基于最小二乘解包裹附Matlab代码

最新推荐文章于 2022-10-25 14:33:17 发布

原创最新推荐文章于 2022-10-25 14:33:17 发布 · 449 阅读

6 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #开发语言

物理应用专栏收录该内容

357 篇文章

订阅专栏

该博客介绍了使用MATLAB进行Zernike多项式拟合和相位解包裹的过程。首先，通过模拟球面进行数据生成，然后进行边界延拓和平滑处理。接着，利用FFT进行相位提取，并采用最小二乘法解包裹相位。最后，通过Zernike多项式进行拟合，得到了精确的解包裹结果。整个过程涉及图像处理、傅里叶变换和数值优化等多个技术环节。

1 简介

2 完整代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%模拟球面  zernikeclcclear allclose allA=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    endendorthop=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,:)');endcoef(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);endcoef_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    endendfigure(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));%%%%%%%lvbofilter0=ones([m n]);filter0(1:(m/2),:)=0;filter1=ones([m n]);filter1((m/2+1):m,:)=0;%%%+1filter2=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;  endendy=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;  endend 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; % endfor x=zx0:m    flt0(x,uint16(b3+a3*x):uint16(b2+a2*x))=1;    %%(-1)*a*x ??endfor x=1:zx0    flt1(x,uint16(b1+a1*x):uint16(b4+a4*x))=1;    endfor y=zy0:n    flt2(uint16((y-b4)/a4):uint16((y-b2)/a2),y)=1;  endfor 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    endendJ=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    endendfor 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    endendX=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)=1;%         F(w,x*n+y-1)=1;;%     end%     if (x>0&&x<m-1)&&(y>1&&y<n)%         F(w,w)=-4;%         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% end% % %%--建立泊松方程% dx=zeros(m,n);% dy=zeros(m,n);% bx=(zeros(m,n));% by=(zeros(m,n));% for lx=1:(m-1)%     for ly=1:n%         %dx(lx,ly)=J(lx+1,ly)-J(lx,ly);%dx(lx,:)=..%         if J(lx+1,ly)-J(lx,ly)>=pi%             bx(lx,ly)=0-1; %         elseif J(lx+1,ly)-J(lx,ly)<-pi%             bx(lx,ly)=0+1;%         else%             bx(lx,ly)=0;%         end%         dx(lx,ly)=J(lx+1,ly)-J(lx,ly)+2*pi*bx(lx,ly);%%NOTE: bb的形式！！！need to change%     end% end% for lx=1:m%     for ly=1:(n-1)%         % dy(lx,ly)=J(lx,ly+1)-J(lx,ly);% !tongshang%         if (J(lx,ly+1)-J(lx,ly))>=pi%             by(lx,ly)=0-1; %         elseif (J(lx,ly+1)-J(lx,ly))<-pi%             by(lx,ly)=0+1;%         else%             by(lx,ly)=0;%         end%         dy(lx,ly)=J(lx,ly+1)-J(lx,ly)+2*pi*by(lx,ly);%%NOTE: bb的形式！！！need to change%     end% end% dx(m,:)=0;% dy(:,n)=0;% R=zeros(m,n);% %             R(1,2:n)=dx(1,2:n)-0+dy(1,2:n)-dy(1,1:(n-1));    % %             R(2:m,1)=dx(2:m,1)-dx(1:(m-1),1)+dy(2:m,1)-0;% %             R(1,1)=dx(lx,ly)+dy(lx,ly);% for lx=1:m%     for ly=1:n%         if (lx==1)&&(ly~=1)  %             R(lx,ly)=dx(lx,ly)+dy(lx,ly)-dy(lx,ly-1);    %         elseif (ly==1)&&(lx~=1)%             R(lx,ly)=dx(lx,ly)-dx(lx-1,ly)+dy(lx,ly);%         elseif (ly==1)&&(lx==1)%             R(lx,ly)=dx(lx,ly)+dy(lx,ly);%         else%             R(lx,ly)=dx(lx,ly)-dx(lx-1,ly)+dy(lx,ly)-dy(lx,ly-1);%         end%     end% end% %%%%%%%%%%%%%解方程%%%%!!!!!!!!!!RR=wiener2(R,[10 10]);% for g=1:m%     for t=1:n%         GG((g-1)*n+t)=R(g,t);%     end% end% GT=GG';% % TT=cgs(F,GT,1e-006,100);%300% X=zeros(m,n);% for g=1:m    %     for t=1:n%         X(g,t)=TT((g-1)*n+t);%/2%%%%%2倍关系%     end% end% figure(8),mesh(double(X));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Zernike拟合%%%%%%A=imresize(X,[65 65],'bilinear');zj=16;[M0 N0]=size(A);coef(1:zj) = 0;RX = (M0+1)/2;RY = (N0+1)/2;R0 =min(RX,RY);% fitted(1:M0,1:N0) = 0;fitted_qulijiao=zeros([M0 N0]);% 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);            z(9,k)=6*r^4-6*r^2+1;            z(10,k)=r^3*cos(3*theta);            z(11,k)=r^3*sin(3*theta);            z(12,k)=(4*r^4-3*r^2)*cos(2*theta);            z(13,k)=(4*r^4-3*r^2)*sin(2*theta);            z(14,k)=(10*r^5-12*r^3+3*r)*cos(theta);            z(15,k)=(10*r^5-12*r^3+3*r)*sin(theta);            z(16,k)=20*r^6-30*r^4+12*r^2-1;        end    endendorthop=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,:)');endcoef(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);endcoef_qulijiao=coef;coef_qulijiao(1: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);            zz(9)=6*r^4-6*r^2+1;               zz(10)=r^3*cos(3*theta);            zz(11)=r^3*sin(3*theta);            zz(12)=(4*r^4-3*r^2)*cos(2*theta);            zz(13)=(4*r^4-3*r^2)*sin(2*theta);            zz(14)=(10*r^5-12*r^3+3*r)*cos(theta);            zz(15)=(10*r^5-12*r^3+3*r)*sin(theta);            zz(16)=20*r^6-30*r^4+12*r^2-1;   %             fitted(i,j) =coef*zz';             fitted_qulijiao(i,j)=coef_qulijiao*zz';         end    endendif coef(4)>0    Result_qulijiao=fitted_qulijiao./(4*pi); else    Result_qulijiao=(-1)*fitted_qulijiao./(4*pi);endfigure(9),mesh(double(Result_qulijiao));