%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%模擬球面??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