生成干涉图的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(RX,RY); 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(i,j)~=0 r = rad/R0; if r~=0 theta = atan2(RX-i,j-RY); end k=k+1; wfront(k)=A(i,j); z(1,k)=1; z(2,k)=r*cos(theta); z(3,k)=r*sin(theta); z(4,k)=2*r^2-1; z(5,k)=r^2*cos(2*theta); z(6,k)=r^2*sin(2*theta); z(7,k)=(3*r^3-2*r)*cos(theta); z(8,k)=(3*r^3-2*r)*sin(theta); end end end orthop=zeros(zj,k); orthop(1,1: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(n,m)=z(n,:)*orthop(m,:)'/(orthop(m,:)*orthop(m,:)'); orthop(n,:)=orthop(n,:)-aa(n,m)*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:zterm,zterm-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-i,j-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(i,j)=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:180,65:180); %figure(3),imshow(uint8(B)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FFT相位提取%%%%%%%%%%%% % %%%%%fft_lvbo%%%%%%%%%%% C=imresize(B,[128 128]); [m,n]=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(Af3_2); Af4_3=fftshift(Af3_3); Aiff5_0=ifft2(Af4_0); Aiff5_1=ifft2(Af4_1); Aiff5_2=ifft2(Af4_2); Aiff5_3=ifft2(Af4_3); %%%%%%%%%%%%%%%%%%%% J0=atan2(imag(Aiff5_0),real(Aiff5_0)); J1=atan2(imag(Aiff5_1),real(Aiff5_1)); J2=atan2(imag(Aiff5_2),real(Aiff5_2)); J3=atan2(imag(Aiff5_3),real(Aiff5_3)); %%%%%%%%%%%%%%%%%%%找中心 x=0; r=ones([2 2]); tp=2; x0=0; for x=3:(m-3) % r=corrcoef(J0(x-1,:),J0(x,:))+corrcoef(J0(x,:),J0(x+1,:)); r=corrcoef(J0(x,:),J0(x+1,:)); if r(1,2)<tp tp=r(1,2); x0=x; end end y=0; r=ones([2 2]); tp=2; y0=0; for y=3:(n-3) % r=corrcoef(J2(:,y-1),J2(:,y))+corrcoef(J2(:,y),J2(:,y+1)); r=corrcoef(J2(:,y),J2(:,y+1)); if r(1,2)<tp tp=r(1,2); y0=y; end end zx0=x0; zy0=y0; x0=x0+0.5; y0=y0+0.5; %%%%%%%%%%% %J=filter0.*J0+filter1.*J1; a1=(y0-1)/(x0-1); b1=(x0-y0)/(x0-1); a2=(n-y0)/(m-x0); b2=(m*y0-n*x0)/(m-x0); a3=(y0-1)/(x0-m); b3=(m*y0-x0)/(m-x0); a4=(y0-n)/(x0-1); b4=(n*x0-y0)/(x0-1); %%%%%%%%%%% flt0=zeros([m n]); flt1=zeros([m n]); flt2=zeros([m n]); flt3=zeros([m n]); x=0;y=0; % for x=(x0+1):m % flt0(x,uint16(b3+a3*x):uint16(b2+a2*x))=1; %%(-1)*a*x ?? % end % for x=1:x0 % flt1(x,uint16(b1+a1*x):uint16(b4+a4*x))=1; % end % for y=(y0+2):n % %flt2((n+2-b*y):(b*y-1),y)=1; % flt2(uint16((y-b4)/a4+2):uint16((y-b2)/a2-2),y)=1; % end % for y=1:(y0-1) % % flt3((b*y+1):(n+1-b*y-1),y)=1; % flt3(uint16((y-b1)/a1+1):uint16((y-b3)/a3-1),y)=1; % end for x=zx0:m flt0(x,uint16(b3+a3*x):uint16(b2+a2*x))=1; %%(-1)*a*x ?? end for x=1:zx0 flt1(x,uint16(b1+a1*x):uint16(b4+a4*x))=1; end for y=zy0:n flt2(uint16((y-b4)/a4):uint16((y-b2)/a2),y)=1; end for y=1:zy0 flt3(uint16((y-b1)/a1):uint16((y-b3)/a3),y)=1; end %%%处理滤波器的重合点 f02=flt0+flt2; f03=flt0+flt3; f12=flt1+flt2; f13=flt1+flt3; f01=flt0+flt1; for x=1:m for y=1:n if f02(x,y)==2 flt2(x,y)=0; end if f03(x,y)==2; flt3(x,y)=0; end if f12(x,y)==2 flt2(x,y)=0; end if f13(x,y)==2; flt3(x,y)=0; end if f01(x,y)==2; flt1(x,y)=0; end end end J=flt0.*J0+flt1.*J1+flt2.*J2+flt3.*J3; figure(6),mesh(double(J)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%j %---基本方法解包裹--- [m n]=size(J); k=double(zeros(m,n)); for j=1:n for h=2:m if (J(h,j)-J(h-1,j))>=pi k(h,j)=k(h-1,j)-1; elseif abs(J(h,j)-J(h-1,j))<pi k(h,j)=k(h-1,j); elseif (J(h,j)-J(h-1,j))<(-pi) k(h,j)=k(h-1,j)+1; end end end for h=1:m for p=2:n if (J(h,p)-J(h,p-1))>=pi k(h,p)=k(h,p-1)-1; elseif abs(J(h,p)-J(h,p-1))<pi k(h,p)=k(h,p-1); elseif (J(h,p)-J(h,p-1))<(-pi) k(h,p)=k(h,p-1)+1; end end end X=J+2*pi*k; figure(7),mesh(double(X)); % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 最小二乘解包裹%%% % %%%%%%%边界处理 % F=spalloc(m*n,m*n,5*m*n); % for w=1:m*n % x=ceil(w/n)-1;y=rem(w,n);% w=i*n+j -->i j--> % if y==0 % y=n; % end % if (x==0) && (y~=1&&y~=n) % F(w,w)=-3;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % F(w,(x+1)*n+y)=1; % % F(w,(x-1)*n+y)=1; % F(w,x*n+y+1)=1; % F(w,x*n+y-1)=1; % end % if (x==m-1) && (y~=1&&y~=n) % F(w,w)=-3;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % % F(w,(x+1)*n+y)=1; % F(w,(x-1)*n+y)=1; % F(w,x*n+y+1)=1; % F(w,x*n+y-1)=1; % end % if (y==1) && (x~=0&&x~=m-1) % F(w,w)=-3;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % F(w,(x+1)*n+y)=1; % F(w,(x-1)*n+y)=1; % F(w,x*n+y+1)=1; % % F(w,x*n+y-1)=1; % end % if (y==n) && (x~=0&&x~=m-1) % F(w,w)=-3;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % F(w,(x+1)*n+y)=1; % F(w,(x-1)*n+y)=1; % % F(w,x*n+y+1)=1; % F(w,x*n+y-1)=1; % end % if (x==0) && (y==1) % &&y~=n) % F(w,w)=-2;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % F(w,(x+1)*n+y)=1; % % F(w,(x-1)*n+y)=1; % F(w,x*n+y+1)=1; % % F(w,x*n+y-1)=1;; % end % if (x==m-1) && (y==1)% &&y~=n) % F(w,w)=-2;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % % F(w,(x+1)*n+y)=1; % F(w,(x-1)*n+y)=1; % F(w,x*n+y+1)=1; % % F(w,x*n+y-1)=1;; % end % if (x==m-1) && (y==n)% &&y~=n) % F(w,w)=-2;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % % F(w,(x+1)*n+y)=1; % F(w,(x-1)*n+y)=1; % % F(w,x*n+y+1)=1; % F(w,x*n+y-1)=1;; % end % if (x==0) && (y==n)% &&y~=n) % F(w,w)=-2;% F(w,(x-1)*n+y)=F(w,(i+1)*n+j) % F(w,(x+1)*n+y)=1; % % F(w,(x-1)*n+y)=1; % % F(w,x*n+y+1)=