相位解包裹(五)枝切法(Goldstein’s branch cut algorithm)

上一篇文章讲了相位解包裹算法的分类,那从这一篇文章开始,就可以开始介绍具体的相位解包裹算法了。

因为介绍过了残点的运算,所以这里首先介绍利用残点的特性去进行相位解包裹算法的一种经典算法,枝切法(Goldstein’s branch cut algorithm)。枝切法的解包裹是和路径有关的,所以其属于空间相位解包裹中的一种。

 

残点的极性(polarity/charge)

从残点的计算我们得知,在一个2*2的小区域里面求相位差的环路积分,结果可能为0,±1(所有结果都除2pi),+1成为正残点(positive),-1成为负残点(negative),这一属性成为残点的极性(polarity)或者charge(实在不知道怎么翻译)

可以得到,相位差的环路积分,实际就是经过多少个正负残点charge的和的2pi倍

\oint{\triangledown \phi \left( r \right)dr=2\pi \times \left( \text{sum of enclosed phase residue charges} \right)}

之前一直说,相位解包裹因为残点存在,是一个与路径有关的过程,那怎么可以让他变成一个与路径无关的过程呢?

从上面的式子可以看出来,要是相位差的环路积分等于0,就需要让charge为0,那怎么才能为0呢,就是即路径经过的区域里,残点的charge被平衡(balance)了,意思是区域中正负残点的数量相等。只要所有区域的正负残点都是平衡的,那相位解包裹的过程自然就和路径没有关系了。

空间相位解包裹的基础就是残点的理论,通过设置一些障碍,例如连接残点的branch,让解包裹的路径绕开这些障碍,解包裹就能顺利进行。

以上讨论都是正常情况,假如认为假如噪声,导致计算的残点不等于±1,这些特殊情况需要另外讨论。

 

枝切法

枝切法就是利用连接残点的branch,让解包裹路径躲开这些区域的相位解包裹的算法。其最重要的贡献还是如何连接残点。

算法的核心:通过枝节(branch)连接距离最近的正负残点,并使它们达到平衡(balance)

 

残点的连接

对于用branch连接残点是有讲究的,如下图(a)从A到B解包裹,有3条路径,明显中间的路径经过了残点的区域,自然得到的结果和另外两条路径会不一样。假如像图(b)一样把正负残点连接使它们平衡,解包裹的路径需要躲开branch,只能走另外两条路的其中一条,不管哪一条,它们的结果都一样。但如果像图(c)一样,通过吧残点和相位图边界连接起来以平衡它们,那么解包裹路径还是会选择中间错误的路。

所以如何 通过branch连接正负残点,是需要一定的法则,如下图,连接残点的branch应该尽可能短。Goldstein等提出的方法,也只是能解决部分残点连接的问题,不是一个完美的方案,在他以后还有很多研究者提出各种各样的方法,要根据情况来决定算法的选取。这里也只讨论Goldstein枝切法一种。

 

Goldstein枝切法残点连接

1、首先当然是需要扫描检测出所有的残点。

2、将一个3*3的搜索窗放在残点的位置,寻找别的残点。如果找到了,将它们用branch连接起来。

    (1)如果另一个残点的极性与当前残点相反,那么就认为这一个branch是已经平衡了,打上“balanced”标记。

    (2)如果另一个残点的极性和当前残点相同,那搜索窗就放在新检测到的这个残点的位置,继续寻找残点,直到找到极性相反的残点并将它们用branch连接起来且他们charge的总和为0(也就是平衡了balanced);或者在搜索窗中再无残点被检测到,如果是后面这种情况,那么搜索窗会增大2,也就是变成5*5的窗,继续从开始的那个残点开始检测别的残点。

 

肥肠重要,一定要注意:

1、当在搜索窗中检测到了残点,都要将它们用branch连接起来。如果检测到的这个残点之前还未被用branch连接过,那它的极性需要加到charge里面来,如果连接过了的就不需要加。(目的就是要用branch上残点的charge总和为0,使它平衡)

2、如果在搜索窗遇到了相位图的边界,那直接把branch和边界连接,同时认定这一branch的charge总和为0。

3、之前讲过,残点实际是一个2*2的小区域,我们标记这个小区域左上角的那个像素为残点,但在残点的连接中,不需要连接被标记的两个像素,而是选择两个小区域最短的路径将他们连接起来,就像下面的图。

 

 

残点连接示例

假设,图中所有残点都已经计算检测出来,并且如图(a),已经连接了一个平衡了的branch,现在从右边的正残点开始检测残点并连接。

1、首先把3*3搜索窗放在右边的正残点上,没有检测到别的残点,当前charge的总和为1。

2、将搜索窗扩大为5*5的搜索窗,没有检测到别的残点当前charge的总和为1。

3、将搜索窗扩大为7*7的搜索窗,如图(b)检测到了箭头指着的正残点,将它们连接起来,由于检测到的这个残点已经连接过了branch,所以当前charge的总和还是为1。

4、7*7的搜索窗内没有检测到新的残点,且charge的总和不为0,将7*7的搜索窗移到箭头所指的残点上继续检测,检测到了相位图的边界,将branch和边界相连,charge的总和为0,这一branch已经平衡,结束检测。

 

 

枝切法相位解包裹

完成残点连接的工作,所有branch都已经平衡了之后,可以开始相位解包裹了,过程类似以flood fill算法。

1、选择一个非branch上的点作为起始点,标记为已解包裹,将它的4邻域的非branch的点,用Itoh方法解包裹,入队

2、只要队列不为空,首个元素出队,将该点的4邻域未曾入队、非branch且未解包裹的点,用Itoh方法解包裹,入队

3、重复步骤2,直到所有非branch上的点全部被解包裹

4、对于在branch上的点,根据他们的邻域中已解包裹的点的相位,对它们逐一解包裹。

 

 

最后,本文的理论和图片,都出自 Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software,我只是当了个翻译工、搬运工,根据自己的理解写出来而已,如果需要更详细的理论,可以自己查阅或者一起交流。

 

参考文献:

[1] Ghiglia D C, Pritt M D. Two-dimensional phase unwrapping: theory, algorithms, and software[M]. New York: Wiley, 1998.

[2] Goldstein R M, Zebker H A, Werner C L. Satellite radar interferometry: Two-dimensional phase unwrapping[J]. Radio science, 1988, 23(4): 713-720.

 

声明:作者水平有限,如文中有错,请务必留言指正。如有学习交流需要,也可通过邮箱zhenyuchung@m.scnu.edu.cn联系我,大家一起讨论学习。

  • 20
    点赞
  • 66
    收藏
    觉得还不错? 一键收藏
  • 22
    评论
Goldstein枝切法相位包裹算法程序的实现可以分为以下几个步骤: 1. 根据相位差值计算出相位差的梯度值。 2. 对梯度值进行高斯滤波,去除噪声。 3. 将滤波后的梯度值映射到[-π,π]范围内。 4. 采用Goldstein枝切法相位进行包裹。 5. 将包裹后的相位进行平滑处理,得到最终的相位包裹结果。 下面是一个简单的Goldstein枝切法相位包裹算法程序的示例: ```python import numpy as np def goldstein_unwrap(phase_diff): """Goldstein枝切法相位包裹算法""" height, width = phase_diff.shape unwrapped_phase = np.zeros((height, width)) # 计算相位差的梯度值 gradient = np.gradient(phase_diff) # 对梯度值进行高斯滤波 sigma = 1 gradient = np.array([np.apply_along_axis(lambda m: np.convolve(m, np.exp(-(np.arange(-5, 6)**2)/(2.*sigma**2))/np.sum(np.exp(-(np.arange(-5, 6)**2)/(2.*sigma**2))), mode='same'), axis=0, arr=row) for row in gradient]) # 将滤波后的梯度值映射到[-π,π]范围内 gradient = np.mod(gradient+np.pi, 2*np.pi)-np.pi # 采用Goldstein枝切法相位进行包裹 for i in range(height): for j in range(width): if j == 0: unwrapped_phase[i,j] = phase_diff[i,j] elif gradient[0][i,j] > np.pi/2: unwrapped_phase[i,j] = unwrapped_phase[i,j-1] - 2*np.pi elif gradient[0][i,j] < -np.pi/2: unwrapped_phase[i,j] = unwrapped_phase[i,j-1] + 2*np.pi else: unwrapped_phase[i,j] = unwrapped_phase[i,j-1] - gradient[0][i,j] # 将包裹后的相位进行平滑处理 for i in range(height): unwrapped_phase[i] = np.convolve(unwrapped_phase[i], np.ones(5)/5, mode='same') return unwrapped_phase ``` 其中,phase_diff是相位差矩阵,可以根据需要自行输入。在实际应用中,还需要根据具体的情况进行参数调整和优化。
评论 22
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值