马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
对相位进行解包的程序,那错了
clear
n=256;
h=fspecial('gaussian',9,1);
I1=imread('image1.bmp');
% figure(1);
% imshow(I1);
M1=rgb2gray(I1); %彩色图变灰度图像
%imshow(M1);
N1=imcrop(M1,[175 105 390 265]);
% figure(2);
% imshow(N1);
A1=filter2(h,N1(:,:,1));
% figure(3)
% imshow(A1)
I2=imread('image2.bmp');
figure(4)
imshow(I2);
M2=rgb2gray(I2); %彩色图变灰度图像
%imshow(M2);
N2=imcrop(M2,[175 105 390 265]);
% figure(5);
% imshow(N2);
A2=filter2(h,N2(:,:,1));
% figure(6)
% imshow(A2)
I3=imread('image3.bmp');
% figure(7)
% imshow(I3);
M3=rgb2gray(I3); %彩色图变灰度图像
%imshow(M3);
N3=imcrop(M3,[175 105 390 265]);
% figure(8);
% imshow(N3);
A3=filter2(h,N3(:,:,1));
% figure(9)
% imshow(A3)
I4=imread('image4.bmp');
% figure(10)
% imshow(I4);
M4=rgb2gray(I4); %彩色图变灰度图像
%imshow(M4);
N4=imcrop(M4,[175 105 390 265]);
% figure(11);
% imshow(N4);
A4=filter2(h,N4(:,:,1));
% figure(12)
% imshow(A4)
for i=1:256
for j=1:256
A4_2(i,j)=A4(i,j)-A2(i,j);
A1_3(i,j)=A1(i,j)-A3(i,j);
m=double(A4_2(i,j)/A1_3(i,j));
phase(i,j)=atan(m);
end
end
figure(13)
mesh(phase);
width=256;
height=256;
%解包
for i=1:width
for j=1:height
if (A1(i,j)==A3(i,j))
if(A2(i,j)>A4(i,j))
phase(i,j)=pi/2;
elseif(A2(i,j)>A3(i,j))
if A2(i,j)>=A4(i,j)
phase(i,j)=atan(m);
else
phase(i,j)=pi-atan(m);
end
else
if( A2(i,j)>=A4(i,j))
phase(i,j)=2*pi-atan(m);
else
phase(i,j)=2*pi+atan(m);
end
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:256
for j=1:40:256
x(m)=i;y(m)=j;
w(m)=phase(i,j);
m=m+1;n=n+1;
end |