CN102800081B - 基于图切割的高抗噪性散斑包裹相位图的展开方法 - Google Patents

Abstract

本发明属于图像处理领域。为解决相位展开算法抗噪性能差，得不到准确解的问题，同时避开噪声的影响，准确提取相位，本发明采取的技术方案是，基于图切割的高抗噪性散斑包裹相位图的展开算法，包括如下步骤：1.针对具体的相位展开问题构造合适的能量函数，建立整数优化的能量最小化模型；2.简化能量最小化模型，将其转化为可迭代求解的0-1优化问题，每次迭代利用图切割理论求解；3.为步骤2中得到的能量函数构造相应的图；4.利用最大流最小割算法确定已建立图的最小割，得到能量函数的最优0-1解，更新能量函数；5.不断迭代来减小能量函数，直到不再减少时终止迭代，得到最优的整数估计，展开相位。本发明主要应用于图像处理。

Description

基于图切割的高抗噪性散斑包裹相位图的展开方法

技术领域

本发明属于图像处理领域，具体涉及剪切散斑干涉术中包裹相位图的展开，即基于图切割的高抗噪性散斑包裹相位图的展开方法。

背景技术

电子剪切散斑干涉技术是20世纪80年代发展起来的一种高精度的光学无损检测技术，广泛应用于粗糙表面的变形测量，具有非接触、高灵敏度和全场实时检测等优点。为了从散斑条纹中得到材料的应变，必须从条纹中提取相位，但是如今的数字散斑测量术中，一般采用相减模式来产生数字散斑条纹图，不可避免地附带有大量的乘性噪声，且通过四步相移法提取相位时，真实相位被包裹在(-π,π]之间，需要将相位展开，得到真实相位，相位展开是相位提取中极其重要的一步，关系到位移、应变等物理量的测量精度。展开包裹的散斑相位图，考虑其大量的乘性噪声，现如今普遍采取先滤波后展开相位的流程，但是散斑条纹图的滤波是非常精巧的任务，过了的话，很容易破坏真实的相位信息，而且针对不同噪声含量的散斑条纹图，滤波策略也得做出相应的改变，这无疑增加了相位展开算法整体的复杂度，故在剪切散斑干涉术中，设计高抗噪性能的相位展开算法十分有意义。

为了能够快速而且准确地展开相位，国内外不同领域的学者提出了众多算法，可以粗略地划分为两类：路径无关法和路径跟踪法。最小二乘法，是与路径无关的全局算法，虽然可以得到平滑解，保证展开相位在全局连续，但是局部的相位展开误差会在全局内传递，无法保证展开相位的相位一致性，即破坏了展开相位与原始相位只可相差2π整数倍的约束条件，得不到精确解。路径跟踪法通过识别残差点后设置枝切，或者引入质量图来指导相位展开的路径，从而避开噪声区，但在散斑包裹相位图的展开中，大量随机的散斑噪声，使得残差点过于密集，很容易形成相位无法展开的孤立区域，无法保证展开相位在全局内连续，以上算法均得不到相对理想的结果。故为了准确展开相位，必须同时保证展开相位的全局连续性和相位一致性。

发明内容

本发明旨在克服现有技术的不足，解决相位展开算法抗噪性能差，得不到准确解的问题，提供一种新型的基于图切割的相位展开算法，避开噪声的影响，准确提取相位，为达到上述目的，本发明采取的技术方案是，基于图切割的高抗噪性散斑包裹相位图的展开方法，包括如下步骤：

1.针对具体的相位展开问题构造合适的能量函数，建立整数优化的能量最小化模型；

2.简化能量最小化模型，将其转化为可迭代求解的0-1优化问题，每次迭代利用图切割理论求解；

3.为步骤2中得到的能量函数构造相应的图，使该图中的割容量表示能量函数；

4.利用最大流最小割算法确定已建立图的最小割，得到能量函数的最优0-1解，更新能量函数；

5.不断迭代来减小能量函数，直到不再减少时终止迭代，得到最优的整数估计，展开相位。

步骤1具体的内容为：四步相移法提取的相位是通过反正切函数得到的，相位被包裹在区间(-π,π]内，相位展开就是将截断的相位恢复成为真实连续的相位，即为：

Φ＝Ψ+2πK      (1)

上式中，Φ是真实相位图，Ψ是包裹相位图，K为整数矩阵；建立一个能量函数来反映全局相位的不连续程度，通过不断迭代整数K来最小化该能量函数，最小的能量函数即对应最优的整数估计

$\hat{K}=\arg \underset{K}{\min }\left\{E\left(K\right)\right\}---\left(2\right)$

能量函数表达为：

$E\left(K\right)=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left(\right|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h|+|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^v\left|\right)---\left(3\right)$

上式中，G为相位值矩阵，(i,j)表示矩阵G第i行，第j列，和分别表示位于相位值矩阵G(i,j)元素水平和竖直方向上的差分，即为：

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h={\mathrm{\φ}}_{\mathrm{ij}-1}-{\mathrm{\φ}}_{\mathrm{ij}}\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v={\mathrm{\φ}}_{i-1j}-{\mathrm{\φ}}_{\mathrm{ij}}\end{array}\right.---\left(4\right)$

上式中，φ_ij为相位图Φ第i行，第j列的元素；φ_ij-1为相位图Φ第i行，第j-1列的元素；φ_i-1j为相位图Φ第i-1行，第j列的元素。

步骤2具体的内容为：每一次的迭代，整数矩阵K中元素k_ij的变化量δ_ij限定为0或1，即t+1次迭代的k_ij可表示为：

$k_{\mathrm{ij}}^{t+1}=k_{\mathrm{ij}}^t+{\mathrm{\δ}}_{\mathrm{ij}}---\left(5\right)$

其中，k_ij的变化量δ_ij∈{0,1}，t+1次迭代时，结合式(1),式(4)和式(5),式(3)中的相位图水平和竖直方向上的差分和可表示为：

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h=2\mathrm{\π}({\mathrm{\δ}}_{\mathrm{ij}-1}-{\mathrm{\δ}}_{\mathrm{ij}})+a^h\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v=2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\end{array}\right.---\left(6\right)$

其中： $a^h=2\mathrm{\π}(k_{\mathrm{ij}-1}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^h,a^v=2\mathrm{\π}(k_{i-1j}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^v,$ 和分别为包裹相位图中水平和竖直方向上的差分。

将式(6)代入式(3)得到：第t+1次迭代时的能量函数E^t+1：

$E(k^t+\mathrm{\δ})=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left[\right|2\mathrm{\π}\left({\mathrm{\δ}}_{\mathrm{ij}-1}{\mathrm{\δ}}_{\mathrm{ij}}\right)+a^h|+|2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\left|\right]---\left(7\right)$

兼顾上式中的水平方向和竖直方向，可以简化为：

$E(k^t+\mathrm{\δ})=\underset{i,j\∈G}{\mathrm{\Σ}}{E}^{\mathrm{ij}}({\mathrm{\δ}}_i,{\mathrm{\δ}}_j)---\left(8\right)$

其中，δ_i,δ_j∈{0,1}，上式中的i,j表示图像中水平或竖直相邻的两个像素点，对整数矩阵K的迭代简化为0-1场δ的迭代。

步骤3的具体内容为：对于图G＝(V,E),由顶点集V和边集E组成，顶点之间由边连接，边被赋予非负的权值，采用有两个终点的图,这两个终点分别被称为源点s和汇点t；把该种图G＝(V,E)中除终点外的点分成两个不相连的子集S和T的过程就是图切割C＝S/T,其中源点s在S集里,汇点t在T集里,切割C相当于标号f，其中f是从一个从顶点集合V-{s,t}到{0,1}的映射：f(v)＝0意味着v∈S，f(v)＝1意味着v∈S；图切割就是用{0,1}给图中的每个顶点赋值；图切割的割容量定义为从集合S到集合T的所有边的权值和，可表达为：

$c(S,T)=\underset{u\∈S,v\∈T,(u,v)\∈E}{\mathrm{\Σ}}c(u,v)---\left(9\right)$

最小割问题就是找到包含最小割容量的切割C。

为了最小化第t+1次迭代时的能量函数E^t+1，构造特定的图，图中的顶点表示能量函数中的变量，变量取值为0或者1，边的权值表示变量的系数，图切割的割容量表示能量函数，确定了最小割就给能量函数赋值，使能量函数最小化；

为其他的单元能量E^ij(δ_i,δ_j)构造相应的子图，然后图切割理论中的可加性定理，将所有的源和汇整合成为单一的源和汇，顶点直接拼接为整体相位的图，这样就构造出了表示整体能量函数的图。

步骤4的具体内容为：根据Ford-Fulkerson定理，确定最小割就等同于计算从源点到汇点的最大流，最小割和最大流是等价的，而图中的最大流可以快速而准确地得到，采用最大流最小割算法max-flow/min-cut先得到图的最大流，进而确定图的最小割。

步骤5具体的内容为：利用求得的能量函数的最小解δ，将其代入式(7),得到t+1次迭代时的能量函数E^t+1,与上一次迭代的能量函数E^t进行比较，不断迭代直到能量函数不再减少，得到最优的整数估计代入下式:

$\mathrm{\Φ}=\mathrm{\Ψ}+2\mathrm{\π}\hat{K}---\left(10\right)$

即在包裹相位图Ψ的基础上叠加2π的倍得到最终的展开相位Φ。

本发明的技术特点及效果：

1.本发明从最优估计的角度解决相位展开问题，建立了相应的理论框架，可以根据不同的相位展开背景，添加相应的约束条件，构造不同的能量函数，相比传统算法，适用范围更广。

2.本发明将图切割理论引入到散斑包裹相位图的展开中，在包裹相位与展开相位相差2π整数倍的约束条件下，建立反映相位全局连续程度的能量函数，为该能量函数构建相应的图，通过一系列0-1优化来最小化能量函数，得到最优的整数估计，从而展开相位。

3.本发明综合了全局算法和路径跟踪算法的优势，获得全局连续相位的同时，可以保证展开相位与原始相位之间很好的一致性。新算法抗噪性能优良，稳定，无需识别残差点，不引入质量图，就可以避开噪声的影响，从散斑包裹相位图中准确提取出真实相位。

附图说明

图1剪切散斑干涉系统中的包裹相位图。

图2四种图切割。图中，(a)(δ_i,δ_j)＝(1,1)，(b)(δ_i,δ_j)＝(1,0)，(c)(δ_i,δ_j)＝(0,1)；(d)(δ_i,δ_j)＝(0,0)。

图3由子图合成的整体图的示意图。

图4算法的流程图。

图5散斑包裹相位图的相位展开结果。图中，(a)为路径跟踪算法的展开结果，(b)为最小二乘法的相位展开结果，(c)图切割算法的相位展开结果。

图6展开相位的再包裹图。图中，(a)为最小二乘算法的展开相位的再包裹图，(b)为图切割算法的展开相位的再包裹图。

具体实施方式

本发明解决的技术问题主要是兼顾展开相位的全局连续性和相位一致性，将相位展开转化为能量最小化问题，并利用图切割不断迭代来取得能量函数的最小解。图切割理论是计算机视觉领域中的一种先进的能量最小化技术，该算法构造合适的图使该图的割容量对应能量函数，求出最小的图切割来最小化能量函数，从而展开相位，本发明具体包括五个部分：

1.针对具体的相位展开问题构造合适的能量函数，建立整数优化的能量最小化模型；

2.简化能量最小化模型，将其转化为可迭代求解的0-1优化问题，每次迭代利用图切割理论求解；

3.为步骤2中得到的能量函数构造相应的图，使该图中的割容量表示能量函数；

4.利用最大流最小割算法确定已建立图的最小割，得到能量函数的最优0-1解，更新能量函数；

5.不断迭代来减小能量函数，直到不再减少时终止迭代，得到最优的整数估计，展开相位。

步骤1具体的内容为：四步相移法提取的相位是通过反正切函数得到的，相位被包裹在区间(-π,π]内，相位展开就是将截断的相位恢复成为真实连续的相位，即为：

Φ＝Ψ+2πK       (1)

上式中，Φ是真实相位图，Ψ是包裹相位图，K为整数矩阵。由上式可知：真实相位与包裹相位的差异为2π的整数倍,相位展开其实就是估计出最优的整数最小化相位图Φ的不连续程度,使其尽可能连续。建立一个能量函数来反映全局相位的不连续程度，通过不断迭代整数K来最小化该能量函数，最小的能量函数即对应最优的整数估计

$\hat{K}=\arg \underset{K}{\min }\left\{E\left(K\right)\right\}---\left(2\right)$

连续性一般由相位图中相邻像素点水平和竖直方向的差分值来表示，能量函数可表达为：

$E\left(K\right)=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left(\right|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h|+|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^v\left|\right)---\left(3\right)$

上式中，G为相位值矩阵，(i,j)为矩阵G第i行，第j列的元素，和分别表示元素(i,j)水平和竖直方向上的差分，即为：

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h={\mathrm{\φ}}_{\mathrm{ij}-1}-{\mathrm{\φ}}_{\mathrm{ij}}\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v={\mathrm{\φ}}_{i-1j}-{\mathrm{\φ}}_{\mathrm{ij}}\end{array}\right.---\left(4\right)$

上式中，φ_ij为相位图Φ第i行，第j列的元素；φ_ij-1为相位图Φ第i行，第j-1列的元素；φ_i-1j为相位图Φ第i-1行，第j列的元素。

步骤2具体的内容为：为了简化能量最小化问题，每一次的迭代，整数K中元素k_ij的变化量δ_ij限定为0或1，即t+1次迭代的k_ij可表示为：

$k_{\mathrm{ij}}^{t+1}=k_{\mathrm{ij}}^t+{\mathrm{\δ}}_{\mathrm{ij}}---\left(5\right)$

其中，δ_ij∈{0,1}，t+1次迭代时，结合式(1),式(4)和式(5),式(3)中的相位图中水平和竖直方向上的差分和可表示为：

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h=2\mathrm{\π}({\mathrm{\δ}}_{\mathrm{ij}-1}-{\mathrm{\δ}}_{\mathrm{ij}})+a^h\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v=2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\end{array}\right.---\left(6\right)$

其中： $a^h=2\mathrm{\π}(k_{\mathrm{ij}-1}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^h,a^v=2\mathrm{\π}(k_{i-1j}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^v,$ 和分别为包裹相位图中水平和竖直方向上的差分。

将式(6)代入式(3)得到：第t+1次迭代时的能量函数E^t+1：

$E(k^t+\mathrm{\δ})=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left[\right|2\mathrm{\π}\left({\mathrm{\δ}}_{\mathrm{ij}-1}{\mathrm{\δ}}_{\mathrm{ij}}\right)+a^h|+|2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\left|\right]---\left(7\right)$

兼顾上式中的水平方向和竖直方向，可以简化为：

$E(k^t+\mathrm{\δ})=\underset{i,j\∈G}{\mathrm{\Σ}}{E}^{\mathrm{ij}}({\mathrm{\δ}}_i,{\mathrm{\δ}}_j)---\left(8\right)$

其中，δ_i,δ_j∈{0,1}，上式中的i,j表示图像中水平或竖直相邻的两个像素点，由上式可以看出，相位展开，可以是不断用0-1场δ进行迭代，逐渐减少能量函数的过程，这样对整数K的迭代简化为0-1场δ的迭代，而图切割算法正好可以有效地解决这个0-1优化问题。

步骤3的具体内容为：对于图G＝(V,E),由顶点集V和边集E组成，顶点之间由边连接，边被赋予非负的权值。本发明中采用有两个终点的图,这两个终点分别被称为源点s和汇点t。把该种类型图中除终点外的点分成两个不相连的子集S和T的过程就是图切割C＝S/T,其中源点s在S集里,汇点t在T集里,切割C相当于标号f，其中f是从一个从顶点集合V-{s,t}到{0,1}的映射：f(v)＝0意味着v∈S，f(v)＝1意味着v∈S。图切割就是用{0,1}给图中的每个顶点赋值。图切割的割容量定义为从集合S到集合T的所有边的权值和，可表达为：

$c(S,T)=\underset{u\∈S,v\∈T,(u,v)\∈E}{\mathrm{\Σ}}c(u,v)---\left(9\right)$

最小割问题就是找到包含最小割容量的切割C。

为了最小化(7)式中的能量函数，可以构造特定的图，图中的顶点表示能量函数中的变量，变量取值为0或者1，边的权值表示变量的系数，图切割的割容量表示能量函数，确定了最小割就给能量函数赋值，使能量函数最小化。

据以上的方法，为其他的单元能量E^ij(δ_i,δ_j)构造相应的子图，然后图切割理论中的可加性定理，将所有的源和汇整合成为单一的源和汇，顶点直接拼接为整体相位的图，这样就构造出了表示整体能量函数的图。

步骤4的具体内容为：根据Ford-Fulkerson定理，确定最小割就等同于计算从源点到汇点的最大流，最小割和最大流是等价的，而图中的最大流可以快速而准确地得到，可用最大流最小割算法(max-flow/min-cut)先得到图的最大流，进而确定图的最小割。

步骤5具体的内容为：利用求得的能量函数的最小解δ，将其代入式(7),得到t+1次迭代时的能量函数E^t+1,与上一次迭代的能量函数E^t进行比较，不断迭代直到能量函数不再减少。得到最优的整数估计代入下式:

$\mathrm{\Φ}=\mathrm{\Ψ}+2\mathrm{\π}\hat{K}---\left(10\right)$

即在包裹相位图Ψ的基础上叠加2π的倍得到最终的展开相位Φ

下面结合附图和实施例进一步说明本发明。

1.算法的具体实施流程

步骤1：新算法进行初始化工作：整数K＝0，迭代次数t＝0,当前t次迭代的相位图Φ^t初始值为包裹相位图Ψ。计算相位图Φ相邻像素点水平和竖直方向上的差分值和

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h={\mathrm{\φ}}_{\mathrm{ij}-1}-{\mathrm{\φ}}_{\mathrm{ij}}\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v={\mathrm{\φ}}_{i-1j}-{\mathrm{\φ}}_{\mathrm{ij}}\end{array}\right.---\left(1\right)$

上式中，φ_ij为相位图Φ第i行，第j列的元素；φ_ij-1为相位图Φ第i行，第j-1列的元素；φ_i-1j为相位图Φ第i-1行，第j列的元素。

代入下式：

$E^t=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left(\right|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h|+|{\mathrm{\Δ\φ}}_{\mathrm{ij}}^v\left|\right)---\left(2\right)$

计算得到t次迭代的能量E^t，建立了相位展开中的能量最小化模型

步骤2：每一次的迭代，整数K中元素k_ij的变化量δ_ij限定为0和1，即K中的元素k_ij可表示为：

$k_{\mathrm{ij}}^{t+1}=k_{\mathrm{ij}}^t+{\mathrm{\δ}}_{\mathrm{ij}}---\left(3\right)$

其中，δ_ij∈{0,1}，则t+1次时的相位图Φ^t+1，相邻像素点水平和竖直方向的差分值和可表示为：

$\left\{\begin{array}{c}{\mathrm{\Δ\φ}}_{\mathrm{ij}}^h=2\mathrm{\π}({\mathrm{\δ}}_{\mathrm{ij}-1}-{\mathrm{\δ}}_{\mathrm{ij}})+a^h\\ {\mathrm{\Δ\φ}}_{\mathrm{ij}}^v=2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\end{array}\right.---\left(4\right)$

其中： $a^h=2\mathrm{\π}(k_{\mathrm{ij}-1}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^h,a^v=2\mathrm{\π}(k_{i-1j}^t-k_{\mathrm{ij}}^t)+{\mathrm{\Δ\ψ}}_{\mathrm{ij}}^v,$ 和分别为包裹相位图中水平和竖直方向上的差分。

对应t+1次迭代时相位图Φ^t+1的能量函数E^t+1为：

$E(k^t+\mathrm{\δ})=\underset{(i,j)\∈G}{\mathrm{\Σ}}\left[\right|2\mathrm{\π}\left({\mathrm{\δ}}_{\mathrm{ij}-1}{\mathrm{\δ}}_{\mathrm{ij}}\right)+a^h|+|2\mathrm{\π}({\mathrm{\δ}}_{i-1j}-{\mathrm{\δ}}_{\mathrm{ij}})+a^v\left|\right]---\left(5\right)$

兼顾上式中的水平方向和竖直方向，可以简化为：

$E(k^t+\mathrm{\δ})=\underset{i,j\∈G}{\mathrm{\Σ}}{E}^{\mathrm{ij}}({\mathrm{\δ}}_i,{\mathrm{\δ}}_j)---\left(6\right)$

其中，δ_i,δ_j∈{0,1}，上式中的i,j表示图像中水平或竖直相邻的两个像素点，这一步得到了能量函数E^t+1的表达式，其中变量δ通过最小化能量函数E^t+1来确定。由上式可以看出，相位展开，可以是不断用0-1场δ进行迭代，逐渐减少能量函数的过程，这样对整数K的迭代简化为0-1场δ的迭代。下一步就是利用图切割理论来最小化能量函数，确定δ。

步骤3：为了最小化(5)式中的能量函数，可以构造特定的图，图中的顶点表示能量函数中的变量，变量取值为0或者1，边的权值表示变量的系数，图切割的割容量表示能量函数，确定了最小割就给能量函数赋值，使能量函数最小化。下面是为(5)式中的能量函数构造特定的图，令式(5)中的每一项可表示为：

$E^{\mathrm{ij}}({\mathrm{\δ}}_i,{\mathrm{\δ}}_j)=E^{\mathrm{ij}}\left(\mathrm{0,0}\right)\left({\stackrel{\‾}{\mathrm{\δ}}}_i\right)\left({\stackrel{\‾}{\mathrm{\δ}}}_j\right)+E^{\mathrm{ij}}\left(\mathrm{0,1}\right)\left({\stackrel{\‾}{\mathrm{\δ}}}_i\right)\left({\mathrm{\δ}}_j\right)+E^{\mathrm{ij}}\left(\mathrm{1,0}\right)\left({\mathrm{\δ}}_i\right)\left({\stackrel{\‾}{\mathrm{\δ}}}_j\right)+E^{\mathrm{ij}}\left(\mathrm{1,1}\right)\left({\mathrm{\δ}}_i\right)\left({\mathrm{\δ}}_j\right)---\left(7\right)$

为了简化表达式：令E^ij(0,0)＝A,E^ij(0,1)＝B,E^ij(1,0)＝C,E^ij(1,1)＝D

则E^ij(δ_i,δ_j)可以表示为：

$E^{\mathrm{ij}}({\mathrm{\δ}}_i,{\mathrm{\δ}}_j)=A+(C-A){\mathrm{\δ}}_i+(D-C){\mathrm{\δ}}_j+(B+C-A-D)\left({\stackrel{\‾}{\mathrm{\δ}}}_i\right)\left({\mathrm{\δ}}_j\right)---\left(8\right)$

由式(5)可知：E^ij(0,0)＝|a|,E^ij(1,1)＝|a|，E^ij(1,0)＝|2π+a|，E^ij(0,1)＝|2π-a|，这样保证了(δ_j)的系数B+C-A-D大于0，下面为C-A>0,D-C<0时的E^ij(δ_i,δ_j)构造相应的图。

构造图中边的权值必须为非负，即得保证E^ij(δ_i,δ_j)表达式中变量的系数为非负，可对(D-C)δ_j做如下变化：

$(D-C){\mathrm{\&delta;}}_j=(C-D)(-{\mathrm{\&delta;}}_j)=(C-D)({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j-1)=(C-D){\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j-C+D---\left(9\right)$

由于常数项不影响最终的最小图切割，舍掉常数项，将式(9)代入式(8)，得到：

$E^{\mathrm{ij}}({\mathrm{\&delta;}}_i,{\mathrm{\&delta;}}_j)=(C-A){\mathrm{\&delta;}}_i+(C-D){\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j+(B+C-A-D)\left({\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_i\right)\left({\mathrm{\&delta;}}_j\right)---\left(10\right)$

这样保证了E^ij(δ_i,δ_j)表达式中变量前的各项系数为非负，可使构造图中边的权值为非负，上式中E^ij(δ_i,δ_j)的构造图以及四种图切割(顶点v₁和v₂的四种标号也即变量对(δ_i,δ_j)四种取值)，如图2所示。

据以上的方法，为其他的单元能量E^ij(δ_i,δ_j)构造相应的子图，然后图切割理论中的可加性定理，将所有的源和汇整合成为单一的源和汇，顶点直接拼接为整体相位的图，这样就构造出了表示式(5)中整体能量函数的图，具体效果如图4所示。

步骤4：根据Ford-Fulkerson定理，确定最小割就等同于计算从源点到汇点的最大流，最小割和最大流是等价的，而图中的最大流可以快速而准确地得到，可用最大流最小割算法(max-flow/min-cut)先得到图的最大流，进而确定图的最小割，这样就得到了能量函数的最小解δ。

步骤5：上一步已经确定了δ，将其代入式(6),得到t+1次迭代时的能量函数E^t+1,与上一次迭代的能量函数E^t进行比较，不断迭代直到能量函数不再减少。得到最优的整数估计代入下式:

$\mathrm{\&Phi;}=\mathrm{\&Psi;}+2\mathrm{\&pi;}\hat{K}---\left(11\right)$

即在包裹相位图Ψ的基础上叠加2π的倍得到最终的展开相位Φ。

图4是算法的具体流程图。

2.算法的处理结果

为了验证新算法的有效性，对剪切散斑干涉系统中的包裹相位图进行相位展开，图1为利用四步相移法获得的包裹相位图，图像大小为512×512像素。新算法的处理结果如图4所示，从图像的视觉效果上，可以看出路径跟踪算法在大量噪声的影响下，虽然可以保证展开相位与包裹相位之间相差2π的整数倍，却很容易形成相位无法展开的孤立区域，无法保证展开相位在全局内连续，而基于图切割的新算法没有出现相位无法展开的孤立区域，保证了展开相位的全局连续。图6是展开相位的再包裹图，其中(a)是最小二乘法的展开相位的再包裹图，(b)是新算法的展开相位的再包裹图，将两者分别与原包裹相位图对比，可以发现：(a)图中的条纹数和原始包裹相位图相比明显减少，而(b)图中两者基本吻合，条纹数没有减少，基于图切割的新算法没有改变相位的动态范围。新算法的相位展开图的再包裹图与原始包裹图的高度一致性体现了新算法的有效性，可以保证展开相位和原始相位的相位一致性，得到精确解。

Claims (1)

1.一种基于图切割的高抗噪性散斑包裹相位图的展开方法，其特征是，包括如下步骤：

(1)针对具体的相位展开问题构造合适的能量函数，建立整数优化的能量最小化模型；

(2)简化能量最小化模型，将其转化为可迭代求解的0-1优化问题，每次迭代利用图切割理论求解；

(3)为步骤(2)中得到的能量函数构造相应的图，使该图中的割容量表示能量函数；

(4)利用最大流最小割算法确定已建立图的最小割，得到能量函数的最优0-1解，更新能量函数；

(5)不断迭代来减小能量函数，直到不再减少时终止迭代，得到最优的整数估计，展开相位；

步骤(1)具体的内容为：

四步相移法提取的相位是通过反正切函数得到的，相位被包裹在区间(-π,π]内，相位展开就是将截断的相位恢复成为真实连续的相位，即为：

Φ＝Ψ+2πK   (1)

上式中，Φ是真实相位图，Ψ是包裹相位图，K为整数矩阵；建立一个能量函数来反映全局相位的不连续程度，通过不断迭代整数矩阵K来最小化该能量函数，最小的能量函数即对应最优的整数估计

$\hat{K}=\arg \underset{K}{\min }\left\{E\left(K\right)\right\}---\left(2\right)$

能量函数表达为：

$E\left(K\right)=\underset{(i,j)\&Element;G}{\mathrm{\&Sigma;}}\left(\right|\mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^h|+|\mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^v\left|\right)---\left(3\right)$

上式中，G为相位值矩阵，(i,j)表示矩阵G第i行，第j列，和分别表示位于相位值矩阵G(i,j)元素水平和竖直方向上的差分，即为：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^h={\mathrm{\&phi;}}_{\mathrm{ij}-1}-{\mathrm{\&phi;}}_{\mathrm{ij}}\\ \mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^v={\mathrm{\&phi;}}_{i-1j}-{\mathrm{\&phi;}}_{\mathrm{ij}}\end{array}\right.---\left(4\right)$

上式中，φ_ij为相位图Φ第i行，第j列的元素；φ_ij-1为相位图Φ第i行，第j-1列的元素；φ_i-1j为相位图Φ第i-1行，第j列的元素；

步骤(2)具体的内容为：

每一次的迭代，整数矩阵K中元素k_ij的变化量δ_ij限定为0或1，即t+1次迭代的k_ij可表示为：

$k_{\mathrm{ij}}^{t+1}=k_{\mathrm{ij}}^t+{\mathrm{\&delta;}}_{\mathrm{ij}}---\left(5\right)$

其中，δ_ij∈{0,1}，t+1次迭代时，结合式(1),式(4)和式(5),式(3)中的相位图水平和竖直方向差分值和可表示为：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^h=2\mathrm{\&pi;}({\mathrm{\&delta;}}_{\mathrm{ij}-1}-{\mathrm{\&delta;}}_{\mathrm{ij}})+a^h\\ \mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{ij}}^v=2\mathrm{\&pi;}({\mathrm{\&delta;}}_{i-1j}-{\mathrm{\&delta;}}_{\mathrm{ij}})+a^v\end{array}\right.---\left(6\right)$

其中： $a^h=2\mathrm{\&pi;}(k_{\mathrm{ij}-1}^t-k_{\mathrm{ij}}^t)+\mathrm{\&Delta;}{\mathrm{\&psi;}}_{\mathrm{ij}}^h,a^v=2\mathrm{\&pi;}(k_{i-1j}^t-k_{\mathrm{ij}}^t)+\mathrm{\&Delta;}{\mathrm{\&psi;}}_{\mathrm{ij}}^v,$ 和分别为包裹相位图中水平和竖直方向上的差分；

将式(6)代入式(3)得到：第t+1次迭代时的能量函数E^t+1：

$E(k^t+\mathrm{\&delta;})=\underset{(i,j)\&Element;G}{\mathrm{\&Sigma;}}\left[\right|2\mathrm{\&pi;}({\mathrm{\&delta;}}_{\mathrm{ij}-1}-{\mathrm{\&delta;}}_{\mathrm{ij}})+a^h|+|2\mathrm{\&pi;}({\mathrm{\&delta;}}_{i-1j}-{\mathrm{\&delta;}}_{\mathrm{ij}})+a^v\left|\right]---\left(7\right)$

兼顾上式中的水平方向和竖直方向，可以简化为：

$E(k^t+\mathrm{\&delta;})=\underset{i,j\&Element;G}{\mathrm{\&Sigma;}}{E}^{\mathrm{ij}}({\mathrm{\&delta;}}_i,{\mathrm{\&delta;}}_j)---\left(8\right)$

其中，δ_i,δ_j∈{0,1}，上式中的i,j表示图像中水平或竖直相邻的两个像素点，对整数矩阵K的迭代简化为0-1场δ的迭代；

步骤(3)的具体内容为：对于图G＝(V,E),由顶点集V和边集E组成，顶点之间由边连接，边被赋予非负的权值，采用有两个终点的图,这两个终点分别被称为源点s和汇点t；把该种图G＝(V,E)中除终点外的点分成两个不相连的子集S和T的过程就是图切割C＝S/T,其中源点s在S集里,汇点t在T集里,切割C相当于标号f，其中f是从一个从顶点集合V-{s,t}到{0,1}的映射：f(v)＝0意味着v∈S，f(v)＝1意味着v∈S；图切割就是用{0,1}给图中的每个顶点赋值；图切割的割容量定义为从集合S到集合T的所有边的权值和，表达为：

$c(S,T)=\underset{u\&Element;S,v\&Element;T,(u,v)\&Element;E}{\mathrm{\&Sigma;}}c(u,v)---\left(9\right)$

最小割问题就是找到包含最小割容量的切割C；

为了最小化第t+1次迭代时的能量函数E^t+1，构造特定的图，图中的顶点表示能量函数中的变量，变量取值为0或者1，边的权值表示变量的系数，图切割的割容量表示能量函数，确定了最小割就给能量函数赋值，使能量函数最小化；

为其他的单元能量E^ij(δ_i,δ_j)构造相应的子图，然后图切割理论中的可加性定理，将所有的源和汇整合成为单一的源和汇，顶点直接拼接为整体相位的图，这样就构造出了表示整体能量函数的图；

步骤(4)的具体内容为：根据Ford-Fulkerson定理，确定最小割就等同于计算从源点到汇点的最大流，最小割和最大流是等价的，而图中的最大流可以快速而准确地得到，采用最大流最小割算法max-flow/min-cut先得到图的最大流，进而确定图的最小割；

步骤(5)具体的内容为：利用求得的能量函数的最小解δ，将其代入式(7),得到t+1次迭代时的能量函数E^t+1,与上一次迭代的能量函数E^t进行比较，不断迭代直到能量函数不再减少，得到最优的整数估计代入下式：

$\mathrm{\&Phi;}=\mathrm{\&Psi;}+2\mathrm{\&pi;}\hat{K}---\left(10\right)$

即在包裹相位图Ψ的基础上叠加2π的倍得到最终的展开相位Φ。