枝切法解包裹

枝切法解包裹【计算过程概述】

打算写几篇博客介绍一下做过的东西——枝切法解包裹（二维图像），主要包含以下四个内容：
1、计算过程概述
2、matlab代码实践
3、C++代码实践
4、matlab调用C++代码
（2、3、4步）还没有时间写

所写内容皆是博主结合网上的资料、论文、现有代码总结改进。如有不足敬请批评指正。

注意：文章不是对理论和原理的理解，而是根据现有方法解释为计算步骤，并能够根据步骤写成代码，是工程性的文章。

参考博客：
二维相位解包裹
残点(Residues)
Itoh解包裹）
经典枝切法
枝切法解包裹

---

---

前言

解包裹，做过结构光、光学干涉等领域的人都不陌生。新入门的可以看一下这位博主的文章，当初做本科毕设的时候就是小白一点都不懂，多亏了这位博主的分享。
二维相位解包裹
其中最简单的行列式（逐行逐列）解包裹就不说了，时间解包裹在编写代码上也并不复杂。
本次系列文章只针对二维解包裹的经典枝切法（Goldstein’s branch cut algorithm）

---

一、算法过程

整个算法流程比较清晰，其中输入是包裹的相位图，输出是解包裹之后的相位图，步骤主要包含：
1、计算残（差）点图。
2、根据残（差）点图生成枝切图。（这个名字自己随便起的，本博客暂且这么叫）
3、根据枝切图（①）和包裹相位图（②）解包裹。
值得一提的是，在Unwrap解包裹步骤中，与质量图解包裹法有些类似。（具体在代码实践里面讲）

二、计算残差点图

在路径跟踪法中，理想情况下解包裹是与路径无关的，逐行逐列解包裹即可。但在实际情况中，由于噪声等各种因素影响，不同的路径解包裹的结果是不一样的。
1987年，Ghiglia, Mastin, Romero等，用不连续点（inconsistencies） 来定义导致与路径有关的点；1988年，Goldstein, Zebker, Werner用residues这个术语，也就是残差点，来描述inconsistencies。
残差点的内容见这篇博客:残点(Residues)

具体计算如下图所示，过程描述如下：定义一个2×2的残差点判断窗口，左上角的像素坐标设为(m,n)，按照一定的顺序（一般取顺时针）计算两个像素之间的相位差，然后再对该差值进行解包裹处理，使其回到单周期内，即公式中的W(·)（这里如果不懂，可以看看代码如何实现）。4个像素可以计算4个差值，随后将4个差值求和得到S。

1. 窗口中的相位差均在（-π，π）内，则S应该在(-2π，2π)，说明四个点连续，此时该窗口并无残点，令窗口内所有值的charge/polarity = 0。
2. S 大于2π，定义为正残点，charge/polarity = 1；
3. S 小于-2π，定义为负残点，charge/polarity = -1。

需要注意的点是：
1、残点实际上是一个在图像上2×2个pixel组成的小窗口，算法为了便利，通常将该窗口的左上角作为残点的标记，其他3个像素仍为0。
（这里只是用左上角标记，后续还要用到这个窗口的其他像素）

2、由于使用了2×2的窗口，因此残差点图中最后一列和最后一行不能计算，设定为0。
（假设包裹相位图的维度是500×500，残差点图维度为500×500，但只算了499×499，剩下的边界元素统一为0）

三、计算枝切图

残点的计算实际上是一个2×2窗口内的环路积分，我们可以将这个窗口放大到整张图（假设是500×500），如果没有图中没有残点，那么解包裹区域内的环路积分charge = 0，相位解包裹就是路径无关。
当残点存在时，为了使500×500图中的解包裹路径无关，那么就需要将里面的残点的charge = 0。在上文中我们定义了残点的charge有+1和-1，也就是说当图中残点都被平衡(balanced)之后，相位解包裹自然就和路径无关了。

枝切法的核心就是如何连接残点，形成枝切图！！！！！

连接残点的目标是使得解包裹区域内残点的charge为0，有2种实现方式：
1、将区域内的正残点与负残点连接，这样两个残点即被平衡。
2、将残点连接到图像边缘（边界），与边界连接的残点自动平衡。

如何连接残点是又一定的规则，如下图b和图c，但要求连接残点的branch应该尽可能短，因此图c为其中的一种连接方式。

不同的枝切法连接方式有所不同，Goldstein方法的主要步骤如下：
（ 在上一步中将所有残点的像素位置保存下来，循环每一个残点。）
1、以其中一个残点为中心，放置一个3×3的窗口，然后寻找窗口内的其他残点或边界。如果找到了就连接起来，如果没有找到，就继续增大窗口到5×5、7×7……，一直到该残点的charge = 0。
在放置的过程中一般会出现 4 种情况：

1）窗口内存在边界，则寻找离边界到该残点的最近距离，然后在枝切图（二值化）用1表示连接路径（下图中的第三个蓝色图所示），连接到边界的残点 charge 令其为0。

2）窗口内存在正负残点，两者正好平衡（Balanced)，随后两个残点的charge都等于0。

3）窗口内存在残点，但两者的charge相同，此时下一步将相同的窗口移到另一个残点上，继续搜索残点或者边界。如下面的第二个图，本来是5×5的窗口搜索到了同样charge的残点，下一步将5×窗口

4）其中这一步跟3）的原理是一样的。当搜索到另一个残点，但它的charge在之前的步骤中已经被置为0，连接branch cut，然后将搜索框放置在另一个残点上。

有一点需要提醒的是：由于之前说残点实际上是一个2×2的区域，因此两个残点之间的连接不是一定要经过被标记的两个像素，而是选择两个小区域最短的路径即可，如下图中，正负残点的连接只需要3个像素。（如果经过标记点则需要4个像素，3个像素最少，因此选择该条路径连接）

残点的连接基本上就是以上几种情况。这一步的输入是所有的残点坐标，输出则是一个二值化的branch cut图。

四、解包裹

计算残点是为了连接branch cut，连接branch cut是为了在解包裹的时候避开这些路径。
解包裹的过程其实不难，过程类似Flood Fill算法。
1、先根据之前的branch cut将图片分为不同的区域，每个区域单独解包裹。（在我处理自己的数据时候，连接完branch cut后，整张图还是连通的。因此有时候这步可以跳过）。
2、选择一个非branch上的点作为起始点，标记为已解包裹，将它的4邻域的非branch的点，用Itoh方法解包裹，将解包裹后的点进入队列。
（这个Itoh方法，实际上就是逐行逐列解包裹时，根据两个像素的相位差进行解包裹，具体看这篇博客文章Itoh解包裹）
3、只要队列不为空，首个元素出队，将该点的4邻域未曾入队、非branch且未解包裹的点用Itoh方法解包裹，解完包裹后入队。
4、重复步骤2，直到所有非branch上的点全部解完包裹。
5、对于在branch上的点，根据他们的邻域中已解包裹的点的相位，对它们逐一解包裹。

解包裹示意图如下：

这里入栈的，其实是已经解完包裹的像素坐标。另外，质量图解包裹的示意图如下。
两者的区别是：质量图解包裹入栈的是未解包裹的像素，并且根据像素质量的大小在栈中进行排序，质量大的像素优先解包裹。

参考文献

[1] Ghiglia D C, Pritt M D. Two-dimensional phase unwrapping: theory, algorithms, and software[M]. New York: Wiley, 1998.
[2] Zhao M, Huang L, Zhang Q, et al. Quality-guided phase unwrapping technique: comparison of quality maps and guiding strategies [J]. Applied Optics, 2011, 50(33): 6214-24.