噪声对二维相位展开的影响

相位展开器通过计算两个连续样本之间的差来检测图像中相位包裹的存在。如果这个差值大于+π，那么相位展开器认为在这个位置有一个包裹。这可能是真正的相位包络，也可能是由于存在噪声而导致的伪包络。为了研究噪声对2D相位展开的影响，让我们将噪声添加到前面图1（a）所示的模拟连续相位图像中。然后我们将对有噪声的相位图像进行包裹。之后，我们将尝试对模拟对象进行相位展开。该过程在如下所示的Matlab程序中实现。噪声方差在这里被设置为值0.4。

从图4中可以看出，如此低水平的添加噪声不会对Itoh展开算法的操作产生不利影响。

1.添加低水平噪声

%This program shows the problems encountered when unwrapping a noisy 2D phase 
%image by using computer simulation 
clc; close all; clear 
N = 512; 
[x,y]=meshgrid(1:N); 
noise_variance = 0.4; 
image1 = 2*peaks(N) + 0.1*x + 0.01*y + noise_variance*randn(N,N); 
figure, colormap(gray(256)), imagesc(image1) 
title('Noisy continuous phase image displayed as visual intensity array') 
xlabel('Pixels'), ylabel('Pixels') 
  
figure  
surf(image1,'FaceColor','interp', 'EdgeColor','none', 'FaceLighting','phong') 
view(-30,30), camlight left, axis tight  
title('Noisy continuous phase image displayed as a surface plot') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 
  
figure, plot(image1(410,:)) 
title('Row 410 of the original noisy continuous phase image') 
xlabel('Pixels'), ylabel('Phase in radians') 
  
%wrap the 2D image 
image1_wrapped = atan2(sin(image1), cos(image1)); 
figure, colormap(gray(256)), imagesc(image1_wrapped) 
title('Noisy wrapped phase image displayed as 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('Noisy wrapped phase image plotted as a surface plot') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 
  
figure, plot(image1_wrapped(410,:)) 
title('Row 410 of the wrapped noisy image') 
xlabel('Pixels'), ylabel('Phase in radians') 
 
%Unwrap the image using the Itoh algorithm: the first method 
%Unwrap the image first by sequentially unwrapping the rows one at a time.  
image1_unwrapped =  image1_wrapped; 
for i=1:N 
image1_unwrapped(i,:) = unwrap(image1_unwrapped(i,:)); 
end  
%Then unwrap all the columns one-by-one  
for i=1:N 
    image1_unwrapped(:,i) = unwrap(image1_unwrapped(:,i)); 
end  
figure, colormap(gray(256)), imagesc(image1_unwrapped) 
title('Unwrapped noisy 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 noisy phase image using the Itoh unwrapper: the first method') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 
  
%Unwrap the image using the Itoh algorithm: the second method 
%Unwrap the image by first sequentially unwrapping all the columns.  
image2_unwrapped =  image1_wrapped; 
for i=1:N 
    image2_unwrapped(:,i) = unwrap(image2_unwrapped(:,i)); 
end  
%Then unwrap all the a rows one-by-one 
for i=1:N 
    image2_unwrapped(i,:) = unwrap(image2_unwrapped(i,:)); 
end 
figure, colormap(gray(256)), imagesc(image2_unwrapped) 
title('Unwrapped noisy image using the 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 noisy phase image using the Itoh algorithm: the second method') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians')

? ? ? a,? ? ? ? ? ? ? b? ? ? ? ? ? ? ? c? ? ? ? ? ? d? ? ? ? ? ?e? ? ? ? ? ? f? ? ? ? ? ? ? ? ? ? g? ? ? ? ?h? ? ? ? ? ? ?

上面是加完噪声之后的：

图5：（a）和（b）一幅有噪声的计算机生成的连续相位图像。（c） &（d）对噪声相位图像进行包裹。（e） &（f）使用Itoh算法的相位展开：第一种方法。（g） &（h）使用Itoh算法的相位展开：第二种方法。噪声方差

1.1代码解读：

1. 生成带噪声的连续相位图像：使用 peaks 函数生成基本相位图像，然后添加线性项 (0.1*x + 0.01*y) 和随机噪声 (noise_variance*randn(N,N))。

2. 显示原始带噪声的连续相位图像：以视觉强度数组和曲面图形式展示这个图像，使用 imagesc 和 surf 函数。

3. 包裹相位图像：使用 atan2(sin(image1), cos(image1)) 将连续相位图像转换为包裹相位图像，处理相位超过 [-π, π] 范围的情况。

4. 显示包裹后的带噪声相位图像：同样以视觉强度数组和曲面图形式展示。

5. 使用Itoh算法解包裹相位图像：展示了两种方法。

   • 第一种方法：先按行（水平方向）逐一解包裹，然后按列（垂直方向）。
   • 第二种方法：先按列逐一解包裹，然后按行。

6. 显示使用Itoh算法解包裹后的图像：分别为两种方法的结果，以视觉强度数组和曲面图形式展示。

1.2 将噪声方差增加到0.6 错误：

当我们将噪声方差增加到0.6的值时，存在问题。在这种情况下，Itoh相位展开算法未能成功地展开该图像。请注意，在展开的相位图像中仍然存在2π不连续性。还要注意，这一次，实现Itoh算法的第一种和第二种方法，现在产生了不同的结果。

在处理带噪声的二维包裹相位图像时，错误累积的问题。具体来说：

1. 错误累积现象：在相位解包裹过程中，尤其是在噪声较多的图像中，很容易产生错误累积。这种错误累积使得解包裹过程变得更加复杂。

2. Itoh算法的应用：使用Itoh算法的第一种方法处理图像时，首先按行（水平方向）逐一解包裹，完成所有行之后，再按列（垂直方向）逐一解包裹。

3. 错误累积的表现：在图中可以观察到错误累积的例子。例如，在图5(e)的第455行出现了一个假的包裹点。这个假包裹点产生了一个2π的错误，这个错误沿着行的方向传播，从假包裹点开始一直到行的末尾。

4. 结果影响：由于这种错误累积，Itoh算法的第一种方法处理后的解包裹相位图像中出现了看起来像是水平线的2π错误。

使用Itoh算法第二种方法处理图像时的错误累积问题。在这种方法中，先逐列解包裹图像，然后逐行处理。例如，图5(g)中第360列出现了一个假包裹点，产生了2π的错误，这个错误沿着该列传播。这种错误在其他列也会发生，导致在最终的解包裹图像中呈现为垂直线状的2π错误。这展示了在处理噪声影响的图像时，错误累积如何影响结果。

2.2D-SRNCP相位展开器

研究人员开发了许多相位展开算法，试图防止误差传播。这些算法中的许多在[2]中进行了解释。此外，在LJMU的通用工程研究所（GERI），我们开发了一种稳健的2D相位展开算法，称为2D-SRNCP相位展开器[3]。我们的算法基于可靠性排序，遵循非连续路径，在处理破坏真实包裹相位图像的噪声方面表现出优异的性能。不要担心它的工作细节，只需将其视为一种非常先进和稳健的展开算法，并将其用作工具即可。您可以通过以下链接在Matlab中下载2D-SRNCP相位展开器http://www.ljmu.ac.uk/geri/90225.htm.

使用2D-SRNCP相位展开器处理图5（c）中所示的包裹相位图像。生成的图像显示为图6（a）中的视觉强度阵列，也显示为如图6（b）所示的3D表面图。将图6（a）和图6（b）分别与图5（a）、图5（b）进行比较表明，我们的算法成功地正确处理了包裹的相位图像，并防止了误差传播。

请注意，2D-SRNCP相位展开器是用C编程语言编写的。此C程序可使用Mex“Matlab可执行文件”动态链接子程序功能从Matlab调用。在调用C代码之前，必须先在Matlab中对其进行编译。要在Matlab中编译C代码，请在Matlab提示下键入以下内容；mex Miguel_2DUnwrapper.cpp下面给出了可用于打开图像的Matlab代码。

%How to call the 2D-SRNCP phase unwrapper from the C language 
%You should have already compiled the phase unwrapper’s C code first 
%If you haven’t, to compile the C code: in the Matlab Command Window type 
%         mex Miguel_2D_unwrapper.cpp 
%The wrapped phase that you present as an input to the compiled C function 
%should have the single data type (float in C)  
WrappedPhase = single(image1_wrapped); 
UnwrappedPhase = Miguel_2D_unwrapper(WrappedPhase); 
figure, colormap(gray(256)) 
imagesc(UnwrappedPhase); 
xlabel('Pixels'), ylabel('Pixels') 
title('Unwrapped phase image using the 2D-SRNCP algorithm') 
  
figure  
surf(double(UnwrappedPhase),'FaceColor','interp', 'EdgeColor','none', 
'FaceLighting','phong') 
view(-30,30), camlight left, axis tight  
title('Unwrapped phase image using the 2D-SRNCP displayed as a surface') 
xlabel('Pixels'), ylabel('Pixels'), zlabel('Phase in radians') 

?

2.1.代码获取?

需要代码的评论区留言，发给你们邮箱。