二维相位展开问题(讨论针对不连续相位展开算法鲁棒性)

作者：Munther Gdeisat博士和Francis Lilley博士

先决条件：为了理解本教程，在阅读本文档之前，您必须已经学习并完成“一维相位展开问题”教程。

有许多应用程序可以生成包裹的相位图像。例如合成孔径雷达（SAR）、磁共振成像（MRI）和条纹图分析。由这些应用程序生成的包裹的相位图像是不可用的，除非它们首先被展开以形成连续的相位图。这意味着，开发一种鲁棒的相位展开算法是所有这些应用的一个重要主题。在这篇文章中，我们不会只在这些应用的特定上下文中讨论相位展开，而是用一般的术语解释二维相位展开问题的概念。

1.二维相位展开简介

我们将如下解释2D相位展开过程。假设我们有一个计算机生成的连续相位图像，它不包含任何相位包裹（2π跳跃）。该图像可以显示为视觉强度阵列，如图1（a）所示。同样的图像也可以绘制为3D表面，如图1（b）所示。该图像（第410行）的单行强度如图1（c）所示。用于生成该相位图像的Matlab代码如下。峰值Matlab函数用于生成连续相位图像。请注意，我们在这里使用术语“连续”不是指模拟信号，而是指不包含任何相位包络的离散1D相位信号或离散2D相位图像。

? ? ? ? ? ? ? ? （a）? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(b)? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(c)

?图1：（a）显示为视觉强度阵列的计算机生成的连续相位图像，（b）绘制为表面的相同图像，（c）相位图像第410行的强度。

%This program is to simulate a continuous phase distribution to act as a dataset 
%for use in the 2D phase unwrapping problem 
clc; close all; clear 
N = 512; 
[x,y]=meshgrid(1:N); 
image1 = 2*peaks(N) + 0.1*x + 0.01*y; 
figure, colormap(gray(256)), imagesc(image1) 
title('Continuous phase image displayed as a visual intensity array') 
xlabel('Pixels'), ylabel('Pixels') 
  
figure  
surf(image1,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong') 
view(-30,30), camlight left, axis tight  
title(' Continuous phase map image displayed as a surface plot') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 
  
figure, plot(image1(410,:)) 
title('Row 410 of the continuous phase image') 
xlabel('Pixels'), ylabel('Phase in radians')

-------------------下面是代码说明---------------------------------?

1. clc; close all; clear：清除MATLAB命令窗口中的内容（clc），关闭所有图形窗口（close all），清除工作空间中的所有变量（clear）。

2. N = 512;：定义变量N，并赋值为512。

3. [x,y]=meshgrid(1:N);：创建一个大小为N×N的网格坐标矩阵x和y。

4. image1 = 2*peaks(N) + 0.1*x + 0.01*y;：生成一个大小为N×N的矩阵image1，其值由MATLAB内置的peaks函数生成的数据，x的线性函数以及y的线性函数组合而成。peaks函数通常用于生成具有峰值和谷值的数据，用于测试图形功能。（）

5. 

6. figure, colormap(gray(256)), imagesc(image1)：创建一个新图形窗口，设置颜色映射为256级灰度，使用imagesc函数显示矩阵image1的图像。这里的图像是以视觉强度数组的形式展示的。

7. title('Continuous phase image displayed as a visual intensity array')：为图像设置标题，“连续相位图像以视觉强度数组形式显示”。

8. xlabel('Pixels'), ylabel('Pixels')：设置x轴和y轴的标签为“Pixels”。

9. figure：创建另一个新的图形窗口。

10. surf(image1,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong')：使用surf函数以表面图的形式显示image1。这里设置了表面颜色插值、边缘颜色为无以及使用Phong光照模型。

11. view(-30,30), camlight left, axis tight：设置视图角度为(-30,30)，加入左侧光源并且让坐标轴紧密适应数据。

12. title(' Continuous phase map image displayed as a surface plot')：为表面图设置标题，“连续相位图映像以表面图形式显示”。

13. xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians')：设置x轴、y轴和z轴的标签，分别为“Pixels”、“Pixels”和“Phase in radians”。

14. figure, plot(image1(410,:))：创建一个新图形窗口，并绘制image1中第410行的一维剖面。

15. title('Row 410 of the continuous phase image')：为这个一维剖面图设置标题，“连续相位图像的第410行”。

16. xlabel('Pixels'), ylabel('Phase in radians')：设置一维剖面图的x轴标签为“Pixels”，y轴标签为“Phase in radians”。

这段代码的目的是为了生成和展示一个连续相位图像的不同视图，包括以灰度图像形式显示的视觉强度数组，以及表面图和一维剖面图。

------------------------------------------------------------------------

现在让我们来包装计算机生成的连续相位图像。执行此任务的Matlab代码如下所示；

%wrap the 2D image 
image1_wrapped = atan2(sin(image1), cos(image1)); 
figure, colormap(gray(256)), imagesc(image1_wrapped) 
title('Wrapped phase image displayed as a visual intensity array') 
xlabel('Pixels'), ylabel('Pixels') 
  
figure  
surf(image1_wrapped,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong') 
view(-30,70), camlight left, axis tight  
title('Wrapped phase image plotted as a surface') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 
  
figure, plot(image1_wrapped(410,:)) 
title('Row 410 of the wrapped phase image') 
xlabel('Pixels'), ylabel('Phase in radians') 

• image1_wrapped = atan2(sin(image1), cos(image1));：这行代码首先计算 image1 中每个元素的正弦和余弦值，然后使用 atan2 函数得到相位的包裹值。包裹相位意味着其值会被限制在 [?π,π] 范围内。这是因为正弦和余弦函数是周期性的，而 atan2 返回的是这两个值的反正切值，这个值在数学上描述了单位圆上的角度，因此也是周期性的。

• figure, colormap(gray(256)), imagesc(image1_wrapped)：创建一个新图形窗口，并使用灰度色彩映射显示包裹相位图像。imagesc 函数会根据数据的最小值和最大值自动调整颜色映射的范围，使得整个数据范围内的细微差别都能在图像中表示出来。

• title('Wrapped phase image displayed as a visual intensity array')、xlabel('Pixels')、ylabel('Pixels')：为图像设置标题和坐标轴标签。

• figure：创建一个新的图形窗口。

• surf(image1_wrapped,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong')：使用 surf 函数创建一个三维表面图，显示包裹相位图像的数据。设置了表面颜色插值、没有边缘颜色和 Phong 光照模型，以便更平滑和真实地呈现图像。

• view(-30,70)、camlight left、axis tight：调整视角、添加光照效果，并使坐标轴紧密包围数据。

• title('Wrapped phase image plotted as a surface')、xlabel('Pixels')、ylabel('Pixels')、zlabel('Phase in radians')：为三维图设置标题和坐标轴标签。

• figure, plot(image1_wrapped(410,:))：创建一个新的图形窗口，并绘制第 410 行的包裹相位数据。这一行的数据以图表的形式表示，能够显示特定行在图像中的相位变化。

• title('Row 410 of the wrapped phase image')、xlabel('Pixels')、ylabel('Phase in radians')：设置图表的标题和坐标轴标签，说明正在显示包裹相位图像的第 410 行的数据。

?包裹后的图像如下所示。

图2：（a）显示为视觉强度阵列的包裹相位图像，（b）绘制为表面的包裹图像，（c）包裹相位图像的行410。

在一维相位解包裹的教程中，相位包裹出现的时候，相位值会在沿线出现多个2π的跳变，形成锯齿波形状。这在图2(c)中有所展示。在二维情况下，我们处理的是以二维数组形式出现的相位图像，相位包裹会表现为轮廓线，如图2(a)所示，我们称这些轮廓线为包裹曲线。这些曲线可能是封闭的，也可能是开放的曲线，在图2(a)中可以看到两种类型的曲线。特别需要注意的是，如果一个开放的曲线进入一个包裹相位图像，那么它也必须离开这个图像。这意味着开放的曲线必须从图像的一边进入并从另一边出去，不能在图像中间断开。

Itoh二维相位解包裹算法主要有两种实现方法。第一种方法是顺序地（一次处理一行）解包裹包裹图像中的行，这会产生一个只进行了部分相位解包裹的中间图像。然后，我们进行类似的过程，但这次是解包裹部分解包裹图像中的所有列。通过这种方式得到的结果是一个解包裹的相位图像，这就是Itoh解包裹算法的第一种实现方式所产生的结果，如图3(a)和(b)所示。执行这项任务的Matlab代码如下。

%Unwrap the image using the Itoh algorithm: the first method is performed 
%by first sequentially unwrapping the all rows, one at a time. 
image1_unwrapped =  image1_wrapped; 
for i=1:N 
    image1_unwrapped(i,:) = unwrap(image1_unwrapped(i,:)); 
end  
%Then sequentially unwrap all the columns one at a time 
for i=1:N 
    image1_unwrapped(:,i) = unwrap(image1_unwrapped(:,i)); 
end  
figure, colormap(gray(256)), imagesc(image1_unwrapped) 
title('Unwrapped phase image using the Itoh algorithm: the first method') 
xlabel('Pixels'), ylabel('Pixels') 
 
figure  
surf(image1_unwrapped,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong') 
view(-30,30), camlight left, axis tight  
title('Unwrapped phase image using the Itoh algorithm: the first method') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 

1. 首先，代码逐行（for i=1:N 循环中的每次迭代代表一行）对图像进行解包裹处理。这是通过MATLAB内置函数 unwrap() 实现的，该函数可以处理一维数组（在这种情况下是图像的每一行）。unwrap() 函数的作用是修改相位数组中的值，消除由于相位超过[?π,π]范围而产生的跳变。

2. 接下来，代码逐列解包裹图像（第二个 for i=1:N 循环），使用的也是 unwrap() 函数。这样就对先前已经逐行解包裹处理过的图像再次进行解包裹处理，但这次是按列进行。

两次解包裹处理完成后，就得到了完全解包裹的相位图像，该图像存储在 image1_unwrapped 变量中。之后，代码使用 imagesc() 函数和 surf() 函数以两种不同的方式展示解包裹后的相位图像：

• imagesc(image1_unwrapped)：以灰度图像的形式显示解包裹的相位图像，其中灰度值代表相位值。

• surf(image1_unwrapped,...)：以三维曲面图的形式显示解包裹的相位图像，其中z轴的高度代表相位值。这种表示方式可以更直观地观察相位的变化。

Itoh算法的第二种实现方法与第一种方法相反。它首先逐列解包裹包裹相位图像中的所有列，这样也会产生一个部分解包裹的图像。然后，它再逐行解包裹这个部分解包裹的图像。使用这种方法得到的解包裹相位图像，就是用第二种Itoh算法实现的结果，如图3(b)所示。这里描述的是一个二维相位图像解包裹的过程，旨在通过先处理列再处理行的方式，来消除包裹效应，以便得到正确的相位信息。

%Unwrap the image using the Itoh algorithm: the second method 
%performed by first sequentially unwrapping all the columns one at a time.  
image2_unwrapped =  image1_wrapped; 
for i=1:N 
    image2_unwrapped(:,i) = unwrap(image2_unwrapped(:,i)); 
end  
%Then sequentially unwrap all the a rows one at a time 
for i=1:N 
    image2_unwrapped(i,:) = unwrap(image2_unwrapped(i,:)); 
end 
figure, colormap(gray(256)), imagesc(image2_unwrapped) 
title('Unwrapped phase image using Itoh algorithm: the second method') 
xlabel('Pixels'), ylabel('Pixels') 
  
figure  
surf(image2_unwrapped,'FaceColor','interp', 'EdgeColor','none', 
'FaceLighting','phong') 
view(-30,30), camlight left, axis tight  
title('Unwrapped phase image using the Itoh algorithm: the second method') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 

图3：使用2D Itoh算法的展开图像；使用（a）和（b）中的第一方法实现并且使用（c）和（d）中的第二方法实现

从上面进行的练习中可以明显看出，Itoh相位展开算法的这两种实现实际上都产生了相同的输出。这是因为这个包裹的相位图像不是真实的，而是一个不包含任何错误的人工数据集。图2（a）所示的包裹相位图像是理想相位图像的一个很好的例子，它不包含任何误差源。我们可以使用任何2D相位展开器轻松处理此图像。如上所述，在这种情况下，我们已经使用2D Itoh算法处理了图像。这是一种非常简单的相位展开算法，仅适用于相位图像几乎无误差的情况。大多数真实世界的应用程序都会生成包含错误的包裹相位图像。在这种情况下，我们需要使用更复杂的二维相位展开器来处理这些图像。在二维相位展开中，有四种误差源使相位展开复杂化。

1.噪音

2.采样不足

3.当连续相位图像包含突然的、突然的相位变化

4.相位提取算法本身产生的误差

在本教程中，我们将仅讨论前三个误差源及其对2D相位展开过程的影响。我们还将解释如何在这三种不同的情况下成功打开图像。第四个误差源取决于用于提取包裹相位的特定算法。读者应该意识到，在执行相位展开时，这是另一个潜在的误差源，然而，关于算法本身对提取的包裹相位的影响的详细讨论不在本教程的范围内，此处不涉及。