高级迭代法（AIA）和通用相位提取（GPSA）算法对比-----MATLAB

%% *******************苏德志论文仿真********************
%                AIA算法和GPSA算法仿真对比(随机步长，移相未知)

%                作者：James_Ray_Murphy

%                参考文献：高精度干涉测量随机移相技术研究_苏志德

%****************************%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;clc;close all;

M = 4;

for k = 1:M

    delta(k) = (k-1)*pi/2 + (2*rand-1)*pi/4;

end

%delta = [pi/10 pi/2-pi/7 pi+pi/7 pi*3/2-pi/8];

N = 200;

xmax = 1;

ymax = 1;

x = linspace(-xmax,xmax,N);

y = linspace(-ymax,ymax,N);

[X,Y] = meshgrid(x,y);

A = 145;%*exp(-1.8*(X.^2+Y.^2));

B = 100;%*exp(-0.2*(X.^2+Y.^2));

phi = 2*pi*(X.^2+Y.^2);

phi_wrapped = angle(exp(j*phi));

figure(1)

surf(phi,'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('初始相位图','FontSize',14)

set(figure(1),'name','Initial Phase 3D','Numbertitle','off');

I0 = A + B.*cos(phi+delta(1));

I1 = A + B.*cos(phi+delta(2));

I2 = A + B.*cos(phi+delta(3));

I3 = A + B.*cos(phi+delta(4));

figure(2);subplot(221);imshow(I0,[]);title('四步非等步长相位图delta=0+noise');

subplot(222);imshow(I1,[]);title('四步非等步长相位图delta=pi/2+noise');

subplot(223);imshow(I2,[]);title('四步非等步长相位图delta=pi+noise');

subplot(224);imshow(I3,[]);title('四步非等步长相位图delta=3pi/2+noise');

%% *********************GPSA算法相位提取***********************

%                   P = K-1.*L

M = 4;

cos_delta = 0;

sin_delta = 0;

sin_cos_delta = 0;

cos_sin_delta = 0;

cos_2_delta = 0;

sin_2_delta = 0;

L1 = 0; L2 = 0; L3 = 0;

delta_G = [0 pi/2 pi 3*pi/2];

tic

for k = 1:M

    cos_delta = cos_delta + cos(delta_G(k));

    sin_delta = sin_delta + sin(delta_G(k));

    cos_sin_delta = cos_sin_delta + cos(delta_G(k))*sin(delta_G(k));

    sin_cos_delta = sin_cos_delta + sin(delta_G(k))*cos(delta_G(k));

    sin_2_delta = sin_2_delta + sin(delta_G(k))^2;

    cos_2_delta = cos_2_delta + cos(delta_G(k))^2;

end

K = [M cos_delta sin_delta;cos_delta cos_2_delta cos_sin_delta;

    sin_delta sin_cos_delta sin_2_delta];

for i = 1:N

    for j = 1:N

        L1 = I0(i,j) + I1(i,j) + I2(i,j) + I3(i,j); 

        L2 = I0(i,j)*cos(delta_G(1)) + I1(i,j)*cos(delta_G(2)) + I2(i,j)*cos(delta_G(3)) + I3(i,j)*cos(delta_G(4));

        L3 = I0(i,j)*sin(delta_G(1)) + I1(i,j)*sin(delta_G(2)) + I2(i,j)*sin(delta_G(3)) + I3(i,j)*sin(delta_G(4));

        L = [L1 L2 L3];

        P = inv(K)*L';

        %***************分子=0的情况***************

        if (-P(3) == 0) 

            if (P(2) == 0)  

                phi_huifu1(i,j) = NaN;

            elseif(P(2) > 0)

                phi_huifu(i,j) = 0;

            elseif(P(2)<0)

                phi_huifu1(i,j) = pi;

            end  

         end 

         %***************分子>0的情况***************

         if (-P(3) > 0) 

            if (P(2) == 0)  

                phi_huifu1(i,j) = pi/2;

            elseif(P(2) > 0)

                phi_huifu1(i,j) = atan(-P(3)/P(2));

            elseif(P(2) < 0)

                phi_huifu1(i,j) = pi + atan(-P(3)/P(2));

            end  

         end 

         %***************分子<0的情况***************

         if (-P(3) < 0) 

            if (P(2) == 0)  

                phi_huifu1(i,j) = -pi/2;

            elseif(P(2) > 0)

                phi_huifu1(i,j) = atan(-P(3)/P(2));

            elseif(P(2) < 0)

                phi_huifu1(i,j) = -pi + atan(-P(3)/P(2));

            end  

         end     

    end

end

toc

figure(3);imshow(phi_huifu1,[]);title('GPSA非等步长相位提取');

%% ***************************AIA算法相位提取相位提取***********************

%                       P = K-1.*L

%*****************************************************************

% Step1：随机设定初值delta，求相位：

for k = 1:M

    %delta_h(k) = delta(k);

    delta_h(k) = (k-1)*pi/2 + (2*rand-1)*pi/8;

end

    delta_1 = delta_h;

    num = 0;

    tic

while (1)

    cos_delta = 0;

    sin_delta = 0;

    sin_cos_delta = 0;

    cos_sin_delta = 0;

    cos_2_delta = 0;

    sin_2_delta = 0;

    L1 = 0; L2 = 0; L3 = 0;

    for k = 1:M

        cos_delta = cos_delta + cos(delta_h(k));

        sin_delta = sin_delta + sin(delta_h(k));

        cos_sin_delta = cos_sin_delta + cos(delta_h(k))*sin(delta_h(k));

        sin_cos_delta = sin_cos_delta + sin(delta_h(k))*cos(delta_h(k));

        sin_2_delta = sin_2_delta + sin(delta_h(k))^2;

        cos_2_delta = cos_2_delta + cos(delta_h(k))^2;

    end

    K = [M cos_delta sin_delta;cos_delta cos_2_delta cos_sin_delta;sin_delta sin_cos_delta sin_2_delta];

    for i = 1:N

        for d = 1:N

            L1 = I0(i,d) + I1(i,d) + I2(i,d) + I3(i,d); 

            L2 = I0(i,d)*cos(delta(1)) + I1(i,d)*cos(delta(2)) + I2(i,d)*cos(delta(3)) + I3(i,d)*cos(delta(4));

            L3 = I0(i,d)*sin(delta(1)) + I1(i,d)*sin(delta(2)) + I2(i,d)*sin(delta(3)) + I3(i,d)*sin(delta(4));

            L = [L1 L2 L3];

            P = inv(K)*L';

        %***************分子=0的情况***************

            if (-P(3) == 0) 

                if (P(2) == 0)  

                    phi_huifu(i,d) = NaN;

                elseif(P(2) > 0)

                    phi_huifu(i,d) = 0;

                elseif(P(2)<0)

                    phi_huifu(i,d) = pi;

                end  

            end 

         %***************分
4000
子>0的情况***************

            if (-P(3) > 0) 

                if (P(2) == 0)  

                    phi_huifu(i,d) = pi/2;

                elseif(P(2) > 0)

                    phi_huifu(i,d) = atan(-P(3)/P(2));

                elseif(P(2) < 0)

                    phi_huifu(i,d) = pi + atan(-P(3)/P(2));

                end  

            end 

         %***************分子<0的情况***************

            if (-P(3) < 0) 

                if (P(2) == 0)  

                    phi_huifu(i,d) = -pi/2;

                elseif(P(2) > 0)

                    phi_huifu(i,d) = atan(-P(3)/P(2));

                elseif(P(2) < 0)

                    phi_huifu(i,d) = -pi + atan(-P(3)/P(2));

                end  

            end 

        end

    end

 %*****************************************************************

% Step2：根据第一步求得的相位，求出delta值：

    NN = N*N;

    cos_phi = 0;

    sin_phi = 0;

    sin_cos_phi = 0;

    cos_sin_phi = 0;

    cos_2_phi = 0;

    sin_2_phi = 0;

    L11 = 0; L12 = 0; L13 = 0;

    L21 = 0; L22 = 0; L23 = 0;

    L31 = 0; L32 = 0; L33 = 0;

    L41 = 0; L42 = 0; L43 = 0;

    for i = 1:N

        for d = 1:N

            %计算K矩阵

            cos_phi = cos_phi + cos(phi_huifu(i,d));

            sin_phi = sin_phi + sin(phi_huifu(i,d));

            cos_sin_phi = cos_sin_phi + cos(phi_huifu(i,d)).*sin(phi_huifu(i,d));

            sin_cos_phi = sin_cos_phi + sin(phi_huifu(i,d)).*cos(phi_huifu(i,d));

            sin_2_phi = sin_2_phi + sin(phi_huifu(i,d)).^2;

            cos_2_phi = cos_2_phi + cos(phi_huifu(i,d)).^2;

            %计算L矩阵

            L11 = L11 + I0(i,d); 

            L12 = L12 + I0(i,d)*cos(phi_huifu(i,d));

            L13 = L13 + I0(i,d)*sin(phi_huifu(i,d));

            LLL(1,:) = [L11 L12 L13];

            %计算L矩阵

            L21 = L21 + I1(i,d); 

            L22 = L22 + I1(i,d)*cos(phi_huifu(i,d));

            L23 = L23 + I1(i,d)*sin(phi_huifu(i,d));

            LLL(2,:) = [L21 L22 L23];

            %计算L矩阵

            L31 = L31 + I2(i,d); 

            L32 = L32 + I2(i,d)*cos(phi_huifu(i,d));

            L33 = L33 + I2(i,d)*sin(phi_huifu(i,d));

            LLL(3,:) = [L31 L32 L33];

            %计算L矩阵

            L41 = L41 + I3(i,d); 

            L42 = L42 + I3(i,d)*cos(phi_huifu(i,d));

            L43 = L43 + I3(i,d)*sin(phi_huifu(i,d));

            LLL(4,:) = [L41 L42 L43];

        end

    end

    K1 = [NN cos_phi sin_phi;cos_phi cos_2_phi cos_sin_phi;sin_phi sin_cos_phi sin_2_phi];

    delta_1(M) = 0;

    PP(M,3) = 0;

    %*************计算delta值******************

    for k = 1:M

        PP(k,:) = inv(K1)*LLL(k,:)';

     %***************分子=0的情况***************

        if (-PP(k,3) == 0) 

            if (PP(k,2) == 0)  

                delta_1(k) = NaN;

            elseif(PP(k,2) > 0)

                delta_1(k) = 0;

            elseif(PP(k,2)<0)

                delta_1(k) = 2*pi;

            end  

        end 

    %***************分子>0的情况***************

        if (-PP(k,3) > 0) 

            if (PP(k,2) == 0)  

                delta_1(k) = pi/2;

            elseif(PP(k,2) > 0)

                delta_1(k) = atan(-PP(k,3)/PP(k,2));

            elseif(PP(k,2) < 0)

                delta_1(k) = pi + atan(-PP(k,3)/PP(k,2));

            end  

        end 

    %***************分子<0的情况***************

        if (-PP(k,3) < 0) 

            if (PP(k,2) == 0)  

                delta_1(k) = pi;

            elseif(PP(k,2) > 0)

                delta_1(k) = pi*2 +atan(-PP(k,3)/PP(k,2));

            elseif(PP(k,2) < 0)

                delta_1(k) = pi + atan(-PP(k,3)/PP(k,2));

            end  

        end  

    end

    num = num +1;

    if abs(delta_1(2)-delta_1(1)-(delta_h(2)-delta_h(1)))<=10^(-4)&abs(delta_1(3)-delta_1(1)-(delta_h(3)-delta_h(1)))<=10^(-4)&abs(delta_1(4)-delta_1(1)-(delta_h(4)-delta_h(1)))<=10^(-4)

        break;

    end

    delta_h = delta_1;

end

toc

figure(4);imshow(phi_huifu,[]);title('AIA算法相位提取');

% figure(3)

% surf(phi_huifu,'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('AIA相位提取','FontSize',14)

% set(figure(3),'name','Initial Phase 3D','Numbertitle','off');

delta_h

delta_1

delta

%mex Miguel_2D_unwrapper.cpp

Phi = single(phi_huifu);

UnwrappedPhase = Miguel_2D_unwrapper(Phi);

UnwrappedPhase = double(UnwrappedPhase);

figure(5);imshow(UnwrappedPhase,[]);colorbar;

%figure;mesh(UnwrappedPhase);

figure(6)

surf(UnwrappedPhase,'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('解包裹后相位图','FontSize',14)

set(figure(6),'name','Initial Phase 3D','Numbertitle','off');

% Un_phase = Unwrap_Alg(phi_huifu);

% figure;mesh(Un_phase);

%% ******************二维对比图******************

figure(7);

plot(phi_wrapped(N/2,:),'b+');

hold on;plot(phi_huifu1(N/2,:),'r*');

hold on;plot(phi_huifu(N/2,:),'g-.');

xlabel('X/Pixels','FontSize',14);ylabel('Phase/Radians','FontSize',14);%title('Two Methods x = N/2','FontSize',14)

legend('原始包裹相位','GPSA相位提取出的包裹相位','AIA相位提取出的包裹相位')

set(figure(7),'name','Comprasion Between Two Methods x = N/2','Numbertitle','off');