matlab 二维相位解缠

一、算法原理

1、相位解包裹算法的原理

   在相位$\phi$的提取过程中，进行了反正切函数的计算，由于三角函数具有周期性，即对于任意的 $\phi ={\phi }_0+2\pi \cdot n$（${\phi }_0$是相位主值， $n$是整数），$\phi$和${\phi }_0$的正切值相同，因此相位$\phi$被截断在反正切函数的主值范围内，仅获得${\phi }_0$的值。通过进一步判断，可将测量的相位结果扩展到 [$-\pi ,\pi$]。当真实相位值超出 [$-\pi ,\pi$]范围时，相位值受三角函数调制而形成 [$-\pi ,\pi$]的包裹相位分布。在测量中，要获得物体表面的真实信息，就必须获得真实的相位分布。
  解包裹相位$\varphi$和包裹相位$\phi$的关系可以用下式表示：
$$\varphi =\phi +2\pi \cdot n,-\pi \le \phi \le \pi \text{\hspace{1em}(1)}$$
   相位解包裹的过程就是确定每一点的整数值，通过在相应的包裹相位点上加上(或减去)$2\pi$的整数倍，得到真实相位值，解包裹关键在于整数值$n$的获得。

2、 相位解包裹的过程

   以一维的相位包裹图为例来描述相位解包裹的过程如图 1 所示。图 $1(a)$是包裹相位分布，图 $1(b)$是解包裹相位分布。从图 $1(a)$到图 $1(b)$是一个解包裹过程，可以看作是一个映射过程，而从图 $1(b)$到图 $1(a)$则是一个逆映射过程。将图 $1(a)$映射到图 $1(b)$可视为一个逐段平移拼接的过程：从左向右，在第一段保持不动时，第二段向上平移，直至第二段的左端点与第一段的右端点连接起来，这时第二段上所有的点都增加了 $2\pi$；第三段向上平移，当与第二段相连时，第三段上的所有点都增加了$4\pi$；第四段上各点只需增加$2\pi$即可与第三段相连；剩余的不需平移已经与第四段连接了起来。完成整个平移过程后，由图 $1(a)$就得到了图 $1(b)$的解包裹相位分布。如果起始点的真实相位与包裹相位间有差值$\pi$，则将解包裹相位值的每一点加上差值$\pi$，得到起始点以零为基点的真实相位分布，如图 $1(c)$所示。通过对上面的分析和总结，可以对相位解包裹的过程做如下描述：

1. 沿某一方向进行解包裹处理时，设起始点的包裹相位值等于该点的解包裹相位 值，在出现周期性的相位跳变之前，各点的解包裹相位值等于对应的包裹相位值。
2. 如果包裹相位值出现跳变，导致相邻两点的相位差值的绝对值大于$\pi$时，将该点和以后各点加上或减去$2\pi$使相邻点的相位差介于$-\pi$到$\pi$之间。
3. 对上述处理后的各点继续进行扫描，对跳变点重复这样的处理，直至所有的相邻点的差值都在$-\pi$到$\pi$之间，最后得到一个连续的相位分布。

  如果起始点的真实相位和包裹相位有一个差值$d$，那么在进行完上述过程后，将所有的解包裹相位值都加上差值$d$，即得到相位真值。在分析解包裹相位值的相对关系时，则不需要知道这个差值。

3、 二维相位解包裹的数学模型

  二维相位解包裹的行列逐点算法，这种算法以一维相位去包裹的数学模型为基础，将二维相位解包裹运算分解成多次一维解包裹运算。该算法首先对二维相位包裹图的第一列进行一维解包裹处理，然后以第一列的解包裹相位为基准，对每一行进行一维解包裹，从而实现整幅图的相位解包裹处理。对一幅 $M\times N$大小的包裹图，需要经过$M+1$次一维解包裹处理。

4、参考文献

[1]马蛟. 基于数字光栅投影的形貌测量技术研究[D]. 合肥工业大学

二、代码实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% unwrap_2D.m 二维相位解缠
% 参考文献：[1]李斌. 基于径向剪切原理的相位解包裹算法研究[D]. 昆明理工大学, 2013.p11
% 由 CSDN 点云侠 创建于 2024/7/30
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc;
clear;
close all;

%% ----------------------------生成模拟曲面--------------------------------
x = linspace(-2,2,40); % 设置范围与步长
y = linspace(-2,2,40);
[X, Y] = meshgrid(x,y);
Z = X.*X+Y.*Y;
% 绘图
figure
surf(X,Y,Z),xlabel('x'),ylabel('y'),zlabel('z');
title('三维曲面');
colormap spring;
hold on;
%% ----------------------转复数图像并获取缠绕相位--------------------------
exp_data = exp(1i*Z);
phase = angle(exp_data); 
% 绘图
figure
surf(X,Y,phase),xlabel('x'),ylabel('y'),zlabel('z');
title('缠绕相位三维显示');
colormap(autumn(5))
hold on;
% 显示复数图像相位
figure
imagesc(phase);% 相位图
colorbar;
title('缠绕相位二维显示');
hold on;
%% ---------------------------二维相位解缠---------------------------------
P2 = unwrap_2D(phase); % 核心代码已封装为函数，这里直接调用即可!

% 绘图
figure
surf(X,Y,P2),xlabel('x'),ylabel('y'),zlabel('z');
title('相位解缠结果');
colormap jet;
hold on;

三、结果展示