2021-08-18

相位解缠的简单算例，贡献给大家，matlab可以直接运行！

clc
clear
%% *************初始相位**************
x=0:0.01:pi;
y=0:0.01:pi;
[x0,y0]=meshgrid(x,y);
z0=cos(x0.*y0);
figure(1);
mesh(z0)
phi0=z0;

%%
phi = angle(exp(1i*phi0));         %包裹相位
figure(2);
imshow(phi,[]);
xlabel('X/Pixels','FontSize',14);ylabel('Y/Pixels','FontSize',14);%title('Wrapped Phase','FontSize',14)
set(figure(2),'name','Wrapped Phase','Numbertitle','off');
axis on
%% *************相位解包裹**************
%********************对包裹相位求一阶偏微分**************
[m,n] = size(phi);
phidx = zeros(m,n);
phidy = zeros(m,n);
phidu = zeros(m,n);
phidv = zeros(m,n);
phidx(1:m-1,:) = angle(exp(j*(phi(2:m,:)-phi(1:m-1,:))));
phidy(:,1:n-1) = angle(exp(j*(phi(:,2:n)-phi(:,1:n-1))));
phidu(1:m-1,1:n-1) = angle(exp(j*(phi(2:m,2:n)-phi(1:m-1,1:n-1))));
phidv(1:m-1,1:n-1) = angle(exp(j*(phi(2:m,1:n-1)-phi(1:m-1,2:n))));
%********************对包裹相位求二阶偏微分**************
Rou3 = zeros(m,n);
Rou3dx = zeros(m,n);
Rou3dy = zeros(m,n);
Rou3du = zeros(m,n);
Rou3dv = zeros(m,n);
Rou3dx(1:m-1,:) = angle(exp(j*(phidx(2:m,:)-phidx(1:m-1,:))));
Rou3dy(:,1:n-1) = angle(exp(j*(phidy(:,2:n)-phidy(:,1:n-1))));
Rou3du(1:m-1,1:n-1) = angle(exp(j*(phidu(2:m,2:n)-phidu(1:m-1,1:n-1))));
Rou3dv(1:m-1,1:n-1) = angle(exp(j*(phidv(2:m,1:n-1)-phidv(1:m-1,2:n))));
Rou3 = Rou3dx + Rou3dy +Rou3du + Rou3dv;

figure(3);
surf(Rou3,'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);
imshow(Rou3,[]);
xlabel('X/Pixels','FontSize',14);ylabel('Y/Pixels','FontSize',14);
set(figure(4),'name','R(x,y) 2D','Numbertitle','off');
%***********************DCT求解泊松方程********************
tic
PP3 = dct2(Rou3);
for ii=1:m
    for jj=1:n
        k1=2*cos((ii-1)*pi/(m));
        k2=2*cos((jj-1)*pi/(n));
        KK = k1+k2 + k1*k2 - 8;
        PH3(ii,jj) = PP3(ii,jj)/KK;
    end
end
PH3(1,1) = -(PH3(1,2) + PH3(2,1) - PP3(1,1))/2;
phi333 = idct2(PH3);
toc
phi333 = phi333(1:m,1:n);

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 Phase Unwrapping','Numbertitle','off');