matlab_Zernike多项式拟合横向剪切波前重建

clear all; num=1024; x=[1:num]; y=[1:num]; [X1,Y1]=meshgrid(x,y); %%%%s=0.2%%%%%%% Radius=300; CenterX=400; CenterY=500; shear_value=120; lamda=532; rad=2*pi/lamda; X2=(X1-CenterX)/Radius; Y2=(Y1-CenterY)/Radius; s=shear_value/Radius; r=X2.^2+Y2.^2; flag1=(r<=1); flag11=(r>1); r1=X2.^2/(1-s)+Y2.^2/(1-s)^2; flag1_elliptical=(r1<=1); flag11_elliptical=(r1>1); r2=Y2.^2/(1-s)+X2.^2/(1-s)^2; flag2_elliptical=(r2<=1); flag22_elliptical=(r2>1); X=zeros(num,num); Y=zeros(num,num); X(flag1)=X2(flag1); Y(flag1)=Y2(flag1); XX=zeros(num,num); YY=zeros(num,num); XX(flag1_elliptical)=X2(flag1_elliptical); YY(flag1_elliptical)=Y2(flag1_elliptical); a0=rad*3; a1=rad*8; a2=rad*10; a3=rad*15; a4=rad*0.7; a5=rad*0.5; a6=rad*0.2; a7=rad*0.15; a8=rad*0.1; a9=rad*0.01; Z0=1; Z1=Y; Z2=X; Z3=2.*X.*Y; Z4=2.*(X.^2+Y.^2)-1; Z5=X.^2-Y.^2; Z6=3.*X.^2.*Y-Y.^3; Z7=3.*X.^2.*Y+3.*Y.^3-2.*Y; Z8=3.*X.^3-3.*X.*Y.^2-2.*X; Z9=X.^3-3.*X.*Y.^2; phi=a0.*Z0+a1.*Z1+a2.*Z2+a3.*Z3+a4.*Z4+a5.*Z5+a6.*Z6+a7.*Z7+a8.*Z8+a9.*Z9; phi(flag11)=0; Vx=phi; Vx_s=zeros(size(Vx)); Vx_s(shear_value+1:num,:)=Vx(1:num-shear_value,:); V_x=Vx-Vx_s; V_x(flag11_elliptical)=0; % figure,imagesc(V_x);axis square [M,N]=size(Vx); Z11=reshape(Z1',1,M*N); Z22=reshape(Z2',1,M*N); Z33=reshape(Z3',1,M*N); Z44=reshape(Z4',1,M*N); Z55=reshape(Z5',1,M*N); Z66=reshape(Z6',1,M*N); Z77=reshape(Z7',1,M*N); Z88=reshape(Z8',1,M*N); Z99=reshape(Z9',1,M*N); Z1r=Z11; Z2r=Z22-(dot(Z22,Z1r)/dot(Z1r,Z1r))*Z1r; Z3r=Z33-(dot(Z33,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z33,Z1r)/dot(Z1r,Z1r))*Z1r; Z4r=Z44-(dot(Z44,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z44,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z44,Z1r)/dot(Z1r,Z1r))*Z1r; Z5r=Z55-(dot(Z55,Z4r)/dot(Z4r,Z4r))*Z4r-(dot(Z55,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z55,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z55,Z1r)/dot(Z1r,Z1r))*Z1r; Z6r=Z66-(dot(Z66,Z5r)/dot(Z5r,Z5r))*Z5r-(dot(Z66,Z4r)/dot(Z4r,Z4r))*Z4r-(dot(Z66,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z66,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z66,Z1r)/dot(Z1r,Z1r))*Z1r; Z7r=Z77-(dot(Z77,Z6r)/dot(Z6r,Z6r))*Z6r-(dot(Z77,Z5r)/dot(Z5r,Z5r))*Z5r-(dot(Z77,Z4r)/dot(Z4r,Z4r))*Z4r-(dot(Z77,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z77,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z77,Z1r)/dot(Z1r,Z1r))*Z1r; Z8r=Z88-(dot(Z88,Z7r)/dot(Z7r,Z7r))*Z7r-(dot(Z88,Z6r)/dot(Z6r,Z6r))*Z6r-(dot(Z88,Z5r)/dot(Z5r,Z5r))*Z5r-(dot(Z88,Z4r)/dot(Z4r,Z4r))*Z4r-(dot(Z88,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z88,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z88,Z1r)/dot(Z1r,Z1r))*Z1r; Z9r=Z99-(dot(Z99,Z8r)/dot(Z8r,Z8r))*Z8r-(dot(Z99,Z7r)/dot(Z7r,Z7r))*Z7r-(dot(Z99,Z6r)/dot(Z6r,Z6r))*Z6r-(dot(Z99,Z5r)/dot(Z5r,Z5r))*Z5r-(dot(Z99,Z4r)/dot(Z4r,Z4r))*Z4r-(dot(Z99,Z3r)/dot(Z3r,Z3r))*Z3r-(dot(Z99,Z2r)/dot(Z2r,Z2r))*Z2r-(dot(Z99,Z1r)/dot(Z1r,Z1r))*Z1r; Z11=Z1r; Z22=Z2r; Z33=Z3r; Z44=Z4r; Z55=Z5r; Z66=Z6r; Z77=Z7r; Z88=Z8r; Z99=Z9r; VV=reshape(Vx',1,M*N); A=zeros(10,10); A(1,1)=num*num; A(1,2)=sum(Z11); A(1,3)=sum(Z22); A(1,4)=sum(Z33); A(1,5)=sum(Z44); A(1,6)=sum(Z55); A(1,7)=sum(Z66); A(1,8)=sum(Z77); A(1,9)=sum(Z88); A(1,10)=sum(Z99); A(2,1)=sum(Z11); A(2,2)=Z11*Z11'; A(2,3)=Z11*Z22'; A(2,4)=Z11*Z33'; A(2,5)=Z11*Z44'; A(2,6)=Z11*Z55'; A(2,7)=Z11*Z66'; A(2,8)=Z11*Z77'; A(2,9)=Z11*Z88'; A(2,10)=Z11*Z99'; A(3,1)=sum(Z22); A(3,2)=Z22*Z11'; A(3,3)=Z22*Z22'; A(3,4)=Z22*Z33'; A(3,5)=Z22*Z44'; A(3,6)=Z22*Z55'; A(3,7)=Z22*Z66'; A(3,8)=Z22*Z77'; A(3,9)=Z22*Z88'; A(3,10)=Z22*Z99'; A(4,1)=sum(Z33); A(4,2)=Z33*Z11'; A(4,3)=Z33*Z22'; A(4,4)=Z33*Z33'; A(4,5)=Z33*Z44'; A(4,6)=Z33*Z55'; A(4,7)=Z33*Z66'; A(4,8)=Z33*Z77'; A(4,9)=Z33*Z88'; A(4,10)=Z33*Z99'; A(5,1)=sum(Z44); A(5,2)=Z44*Z11'; A(5,3)=Z44*Z22'; A(5,4)=Z44*Z33'; A(5,5)=Z44*Z44'; A(5,6)=Z44*Z55'; A(5,7)=Z44*Z66'; A(5,8)=Z44*Z77'; A(5,9)=Z44*Z88'; A(5,10)=Z44*Z99'; A(6,1)=sum(Z55); A(6,2)=Z55*Z11'; A(6,3)=Z55*Z22'; A(6,4)=Z55*Z33'; A(6,5)=Z55*Z44'; A(6,6)=Z55*Z55'; A(6,7)=Z55*Z66'; A(6,8)=Z55*Z77'; A(6,9)=Z55*Z88'; A(6,10)=Z55*Z99'; A(7,1)=sum(Z66); A(7,2)=Z66*Z11'; A(7,3)=Z66*Z22'; A(7,4)=Z66*Z33'; A(7,5)=Z66*Z44'; A(7,6)=Z66*Z55'; A(7,7)=Z66*Z66'; A(7,8)=Z66*Z77'; A(7,9)=Z66*Z88'; A(7,10)=Z66*Z99'; A(8,1)=sum(Z77); A(8,2)=Z77*Z11'; A(8,3)=Z77*Z22'; A(8,4)=Z77*Z33'; A(8,5)=Z77*Z44'; A(8,6)=Z77*Z55'; A(8,7)=Z77*Z66'; A(8,8)=Z77*Z77'; A(8,9)=Z77*Z88'; A(8,10)=Z77*Z99'; A(9,1)=sum(Z88); A(9,2)=Z88*Z11'; A(9,3)=Z88*Z22'; A(9,4)=Z88*Z33'; A(9,5)=Z88*Z44'; A(9,6)=Z88*Z55'; A(9,7)=Z88*Z66'; A(9,8)=Z88*Z77'; A(9,9)=Z88*Z88'; A(9,10)=Z88*Z99'; A(10,1)=sum(Z99); A(10,2)=Z99*Z11'; A(10,3)=Z99*Z22'; A(10,4)=Z99*Z33'; A(10,5)=Z99*Z44'; A(10,6)=Z99*Z55'; A(10,7)=Z99*Z66'; A(10,8)=Z99*Z77'; A(10,9)=Z99*Z88'; A(10,10)=Z99*Z99'; B=zeros(10,1); B(1,1)=sum(VV); B(2,1)=Z11*VV'; B(3,1)=Z22*VV'; B(4,1)=Z33*VV'; B(5,1)=Z44*VV'; B(6,1)=Z55*VV'; B(7,1)=Z66*VV'; B(8,1)=Z77*VV'; B(9,1)=Z88*VV'; B(10,1)=Z99*VV'; b=ones(10,1); b=A\B; a=[a0,a1,a2,a3,a4,a5,a6,a7,a8,a9]'; c=b-a; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%以上为测试%%%%%%%%%%%%%%%%%%%%%%%%%% Vy_s=zeros(size(Vx)); Vy_s(:,shear_value+1:num)=Vx(:,1:num-shear_value); V_y=Vx-Vy_s; V_y(flag22_elliptical)=0; [M,N]=size(Vx); Z11=reshape(Z1',1,M*N); Z22=reshape(Z2',1,M*N); Z33=reshape(Z3',1,M*N); Z44=reshape(Z4',1,M*N); Z55=reshape(Z5',1,M*N); Z66=reshape(Z6',1,M*N); Z77=reshape(Z7',1,M*N); Z88=reshape(Z8',1,M*N); Z99=reshape(Z9',1,M*N); VVy=reshape(V_y',1,M*N); By=zeros(10,1); By(1,1)=sum(VVy); By(2,1)=Z11*VVy'; By(3,1)=Z22*VVy'; By(4,1)=Z33*VVy'; By(5,1)=Z44*VVy'; By(6,1)=Z55*VVy'; By(7,1)=Z66*VVy'; By(8,1)=Z77*VVy'; By(9,1)=Z88*VVy'; By(10,1)=Z99*VVy'; by=ones(10,1); by=A\By; Mx=zeros(10,6); Mx(:,1)=[s,0,0,0,0,0,0,0,0,0]'; Mx(:,2)=[0,2*s,0,0,0,0,0,0,0,0]'; Mx(:,3)=[s^2,0,2*s,0,0,0,0,0,0,0]'; Mx(:,4)=[0,3,s^2,0,3*s,0,0,0,0,0]'; Mx(:,5)=[3,s^3+s,0,9,s^2,0,3*s,3*s,0,0]'; Mx(:,6)=[s^3,0,3*s^2,0,0,3*s,0,0,0,0]'; My=zeros(10,6); My(:,1)=[s,0,0,0,0,0,0,0,0,0]'; My(:,2)=[-s^2,-2*s,0,0,0,0,0,0,0,0]'; Mx(:,3)=[0,0,2*s,0,0,0,0,0,0,0]'; Mx(:,4)=[0,0,3*s^2,3*s,0,0,0,0,0,0]'; My(:,5)=[-s^3,-3*s^2,0,0,0,3*s,0,0,0,0]'; My(:,6)=[s+3*s^3,9*s^2,0,0,3*s,-3*s,0,0,0,0]'; kx=Mx\b; ky=My\by;