phaseunwrap11.m

相位图像的相位解包算法对一个相位图进行解包裹

clear all
clc
h=fspecial('gaussian',9,1);
%度如图 并求平均值
I1=imread('a01.bmp','bmp');
I1=filter2(h,double(I1(:,:,1)));
I2=imread('a02.bmp','bmp');
I2=filter2(h,double(I2(:,:,1)));
I3=imread('a03.bmp','bmp');
I3=filter2(h,double(I3(:,:,1)));
I4=imread('a04.bmp','bmp');
I4=filter2(h,double(I4(:,:,1)));

% I1=double(I1(201:400,201:500));
% I2=double(I2(201:400,201:500));
% I3=double(I3(201:400,201:500));
% I4=double(I4(201:400,201:500));
I1=double(I1);
I2=double(I2);
I3=double(I3);
I4=double(I4);
[width,height]=size(I1);
% width=200;
% height=300;
%解包
for i=1:width
    for j=1:height
        if (I1(i,j)==I3(i,j))
            if(I2(i,j)>I4(i,j))
                phase(i,j)=pi/2;
            elseif(I2(i,j)<I4(i,j))
                phase(i,j)=3*pi/2;
            else
                phase(i,j)=0;
            end
        elseif (I1(i,j)>I3(i,j))
            if I2(i,j)>=I4(i,j)
                phase(i,j)=atan((I2(i,j)-I4(i,j))/(I1(i,j)-I3(i,j)));
            else
                phase(i,j)=2*pi-atan((I4(i,j)-I2(i,j))/(I1(i,j)-I3(i,j))); 
            end
        else
            if I2(i,j)>=I4(i,j)
                phase(i,j)=pi-atan((I2(i,j)-I4(i,j))/(I3(i,j)-I1(i,j)));
            else
                phase(i,j)=pi+atan((I4(i,j)-I2(i,j))/(I3(i,j)-I1(i,j))); 
            end
        end
    end
end
phase=unwrap(phase);
%消除2pi突变
% midx=fix(width/2);     
% for j=1:width
%     for m=height-1:-1:1
%         a=phase(j,m)-phase(j,m+1);
%         if a>pi
%            for n=m:-1:1
%                phase(j,n)=phase(j,n)-2*pi;
%            end
%         end
%         if a<-pi
%            for n=m:-1:1
%             phase(j,n)=phase(j,n)+2*pi;
%            end
%         end
%     
%         if a==pi
%            if phase(j,m)-phase(j,m-1)>0
%               for n=m:-1:1
%                   phase(j,n)=phase(j,n)-2*pi;
%               end
%            end
%         end
%         if a==-pi
%            if phase(j,m)-phase(j,m-1)<0
%               for n=m:-1:1
%                   phase(j,n)=phase(j,n)+2*pi;
%               end
%            end
%         end  
%     end
% end
% 
% for k=1:height    
%     for m=midx:-1:1
%         c=phase(m,k)-phase(m+1,k);
%         if c>pi
%            for n=m:-1:1
%                phase(n,k)=phase(n,k)-2*pi;
%            end
%         end
%         if c<-pi
%            for n=m:-1:1
%                phase(n,k)=phase(n,k)+2*pi;
%            end
%         end
%     
%         if c==pi
%            if phase(m,k)-phase(m,k-1)>0
%               for n=m:-1:1
%                phase(n,k)=phase(n,k)-2*pi;;
%               end
%            end
%         end
%         if c==-pi
%            if phase(m,k)-phase(m,k-1)<0
%               for n=m:-1:1
%                  phase(n,k)=phase(n,k)+2*pi;
%               end
%            end
%         end 
%     end
%     for m=midx:width
%         b=phase(m,k)-phase(m-1,k);
%         if b>pi
%            for n=m:width
%                phase(n,k)=phase(n,k)-2*pi;
%            end
%         end
%         if b<-pi
%            for n=m:width
%                phase(n,k)=phase(n,k)+2*pi;
%            end
%         end
%     
%         if c==pi
%            if phase(m,k)-phase(m+1,k)>0
%               for n=m:-1:1
%                   phase(n,k)=phase(n,k)-2*pi;;
%               end
%            end
%          end
%          if c==-pi
%             if phase(m,k)-phase(m+1,k)<0
%                for n=m:-1:1
%                    phase(n,k)=phase(n,k)+2*pi;
%                end
%             end
%          end
%      end
% end
phase=phase/(2*pi);
figure,surf(phase)

m=1;n=1;
for i=1:40:200
    for j=1:40:200
        x(m)=i;y(m)=j;
        w(m)=phase(i,j);
        m=m+1;n=n+1;
    end
end
[a b]=size(x);
a=zn(x,y,w,b);
for i=1:200
    for j=1:200
        n=200*(i-1)+j;
        xn(n)=i;yn(n)=j;
        ss=zernikeCT(i,j,36);
        ww(n)=sum(ss.*a(1:35))+a(36);
    end
end
figure,plot3(xn,yn,ww)