相位解缠Matlab代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%模拟球面  zernike
A=1000000000000000*ones（[115 115]）;
zj=8;
[M0 N0]=size（A）;
coef=zeros（[1 zj]）;
RX = （M0+1）/2;
RY = （N0+1）/2;
R0 =min（RXRY）;
fitted_ball=zeros（[M0 N0]）;
r=0;k=0;
z=0;
wfront=0;
for i = 1:M0
    for j = 1:N0
        rad = sqrt（（i-RX）^2+（j-RY）^2）;%椭圆的问题？
        if rad <= R0;% &&A（ij）~=0
            r = rad/R0;
            if r~=0
                theta = atan2（RX-ij-RY）;
            end
            k=k+1;
            wfront（k）=A（ij）;
            z（1k）=1;
            z（2k）=r*cos（theta）;
            z（3k）=r*sin（theta）;
            z（4k）=2*r^2-1;
            z（5k）=r^2*cos（2*theta）;
            z（6k）=r^2*sin（2*theta）;
            z（7k）=（3*r^3-2*r）*cos（theta）;
            z（8k）=（3*r^3-2*r）*sin（theta）;
        end
    end
end
orthop=zeros（zjk）;
orthop（11:k）=z（1:）;
bb（1）=wfront*orthop（1:）‘/（orthop（1:）*orthop（1:）‘）;
zterm=zj;
for n=2:zterm
    orthop（n:）=z（n:）;
    for m=1:n-1
        aa（nm）=z（n:）*orthop（m:）‘/（orthop（m:）*orthop（m:）‘）;
        orthop（n:）=orthop（n:）-aa（nm）*orthop（m:）;
    end
    bb（n）=wfront*orthop（n:）‘/（orthop（n:）*orthop（n:）‘）;
end
coef（zterm）=bb（zterm）;
for n=1:zterm-1
    coef（zterm-n）=bb（zterm-n）-coef（zterm-n+1:zterm）*aa（zterm-n+1:ztermzterm-n）;
end
coef_ball=coef;
coef_ball（1:3）=0;coef_ball（5:zj）=0;
% coef_qiumian（1）=0;coef_qiumian（4）=0;
r=0;
zz=zeros（[1 zterm]）;
for i = 1:M0
    for j = 1:N0
        rad = sqrt（（i-RX）^2+（j-RY）^2）;
        if rad <= R0
            r = rad/R0;
            if r~=0
                theta = atan2（RX-ij-RY）;
            end
            zz（1）=1;
            zz（2）=r*cos（theta）;
            zz（3）=r*sin（theta）;
            zz（4）=2*r^2-1;
            zz（5）=r^2*cos（2*theta）;
            zz（6）=r^2*sin（2*theta）;
            zz（7）=（3*r^3-2*r）*cos（theta）;
            zz（8）=（3*r^3-2*r）*sin（theta）;
        fitted_ball（ij）=coef_ball*zz‘;
        end
    end
end
figure（1）mesh（fitted_ball）;
 title（‘三维面形图‘‘FontSize‘10）;
 xlabel（‘x轴（pix）‘‘FontSize‘8）;
 ylabel（‘y轴（pix）‘‘FontSize‘8）;
 zlabel（‘面形‘‘FontSize‘8）; 
gst0=255*cos（fitted_ball*20）;
gst1=imresize（（gst0）[230 230]‘bilinear‘）;
figure（2）imshow（uint8（gst1））;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 边界延拓
B=gst1（65:18065:180）;
figure（3）imshow（uint8（B））;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         FFT相位提取%%%%%%%%%%%%
% %%%%%fft_lvbo%%%%%%%%%%%
C=imresize（B[128 128]）;
[mn]=size（C）;
% Af=fft2（double（C））;
Af1=fft2（C）;      
Af2=fftshift（Af1）;
figure（4）imshow（uint8（Af2））;
%%%%%%%lvbo
filter0=ones（[m n]）;
filter0（1:（m/2）:）=0;
filter1=ones（[m n]）;
filter1（（m/2+1）:m:）=0;%%%+1
filter2=ones（[m n]）;
filter2（:1:（n/2））=0;
filter3=ones（[m n]）;
filter3（:（n/2+1）:n）=0;%%%+1
%Af2=abs（Af2）;
Af3_0=filter0.*Af2;
Af3_1=filter1.*Af2;
Af3_2=filter2.*Af2;
Af3_3=filter3.*Af2;
%figure（5）imshow（uint8（Af3_0））;
Af4_0=fftshift（Af3_0）;
Af4_1=fftshift（Af3_1）;
Af4_2=fftshift