一、算法原理
1、四向最小二乘
2、参考文献
[1]郭媛,吴全,毛琦,等. 基于横向剪切的四向最小二乘相位解包裹算法 [J]. 光电子·激光, 2015, 26 (10): 1953-1959.
二、代码实现
clc;
clear;
close all;
%% --------------------------------初始相位--------------------------------
phi0 = peaks(200)*3;
figure(1);
surf(phi0);
title('山峰3D图')
%% -------------------------------包裹相位---------------------------------
phi = angle(exp(1i*phi0));
figure(2);
%imshow(phi,[]);
imagesc(phi)
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
%title('Wrapped Phase','FontSize',14)
set(figure(2),'name','包裹相位','Numbertitle','off');
axis on
%% ------------------------------相位解包裹--------------------------------
% 对包裹相位求一阶偏微分
[M,N] = size(phi);
deltX = zeros(M,N);
deltY = zeros(M,N);
deltP = zeros(M,N);
deltQ = zeros(M,N);
% deltX(1:M-1,:) = phi(2:M,:)-phi(1:M-1,:);
% deltX = deltX - pi*round(deltX/pi); % 去除差分中的跳跃
deltX(1:M-1,:) = angle(exp(1i*(phi(2:M,:)-phi(1:M-1,:)))); % 公式(2)
deltY(:,1:N-1) = angle(exp(1i*(phi(:,2:N)-phi(:,1:N-1)))); % 公式(3)
deltP(1:M-1,1:N-1) = angle(exp(1i*(phi(2:M,2:N)-phi(1:M-1,1:N-1)))); % 公式(4)
deltQ(1:M-1,1:N-1) = angle(exp(1i*(phi(2:M,1:N-1)-phi(1:M-1,2:N)))); % 公式(5)
% 对包裹相位求二阶偏微分
DdeltX = zeros(M,N);
DdeltY = zeros(M,N);
DdeltU = zeros(M,N);
DdeltV = zeros(M,N);
DdeltX(1:M-1,:) = angle(exp(1i*(deltX(2:M,:)-deltX(1:M-1,:)))); % 公式(10)
DdeltY(:,1:N-1) = angle(exp(1i*(deltY(:,2:N)-deltY(:,1:N-1))));
DdeltU(1:M-1,1:N-1) = angle(exp(1i*(deltP(2:M,2:N)-deltP(1:M-1,1:N-1))));
DdeltV(1:M-1,1:N-1) = angle(exp(1i*(deltQ(2:M,1:N-1)-deltQ(1:M-1,2:N))));
Rou = DdeltX + DdeltY + DdeltU + DdeltV; % 公式(8)
Wij = sqrt(DdeltX.^2 + DdeltY.^2 + DdeltU.^2 + DdeltV.^2); % 公式(9)
k = 0; % 用于调节权重强度,其值由实际情况确定。
Rou_ij = (1+k*Wij).*Rou; % 公式(11)
%% --------------------------DCT求解泊松方程-------------------------------
% 基于离散余弦变换DCT计算ρ_ij
tic
pp = dct2(Rou)+eps;
% 计算Φ_ij在DCT域的精确解
fi = zeros(M,N);
for m = 1:M
for n = 1:N
k1 = 2*cos(pi*(m-1)/M);
k2 = 2*cos(pi*(n-1)/N);
kk = k1 + k2 - 4;
fi(m,n) = pp(m,n)/(kk+eps);
end
end
fi(1,1) = pp(1,1); % 赋值DCT域的Φ11
phs = idct2(fi); % 用iDCT计算解缠相位在空域中的值
toc
phi333 = phs;
figure(3);
surf(Rou,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
zlabel('Phase/Radians','FontSize',14);
%title('lou3','FontSize',14)
set(figure(3),'name','R(x,y) 3D','Numbertitle','off');
figure(4);
imagesc(Rou)
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
set(figure(4),'name','R(x,y) 2D','Numbertitle','off');
figure(5);
surf(phi333,'FaceColor','interp', 'EdgeColor','none','FaceLighting','phong');
camlight left, axis tight
xlabel('X/Pixels','FontSize',14);
ylabel('Y/Pixels','FontSize',14);
zlabel('Phase/Radians','FontSize',14);
%title('Four LSBLS Phase Unwrapping','FontSize',14)
set(figure(5),'name','Four LSBLS 相位解缠结果','Numbertitle','off');