clear
clc
close all
%%
z=read_xyz('shiyan1.xyz');
figure;imshow(z);
I0=z(20:440,200:620); %选取有合适的大小
% figure;imshow(I0);
% figure;mesh(I0);
k=I0;
%% 36项系数aa  直角坐标
[m,n]=size(k);
X=-1:2/(m-1):1;
[xa,ya]=meshgrid(X);xx=xa;yy=ya;
[theta1,r1]=cart2pol(xa(:),ya(:));
is_in_circle=(r1<=0.9);   %控制拟合半径
xa=xa(:);ya=ya(:);za=k(:);
xa(~is_in_circle)=[];
ya(~is_in_circle)=[];
za(~is_in_circle)=[];
b=length(xa);
aa=zernike(xa,ya,za,b);
%%
% aa(1:2)=0;
aa(36)=0;
% aa(3)=0;% 离焦项
temp=zeros(b,1);
ss=zernikeCT(xa,ya,35);
for i=1:35
    temp=temp+ss(:,i)*aa(i);
end
temp=temp+aa(36);
result=nan(m,n);
result(is_in_circle)=temp;%result就是你要的三维波面的数据
figure;mesh(xx,yy,result);title('泽尼克拟合');
xlabel('x/pixel');ylabel('y/pixel');zlabel('phase（wave）');
%% 求波相差的PV值及RMS值　
idx=(r1<=1);
[PV RMS]=PV_RMS(result,idx);