FDDCT.rar_包裹_相位展开_相位解包裹_解包裹算法_解相位

%% 四种最小二乘解包裹算法比较 clear close all O1=peaks(64); O2=0.4*randn(64); O1=O1+O2; U=exp(j.*O1); figure,surf(O1); grid on; axis on; shading interp; ph_wrap=angle(U); figure,imshow(ph_wrap,[]); P=ph_wrap; [M,N]=size(P); %确定包裹相位矩阵的大小 % % M=128; % % N=128; % % fi=zeros(M,N); % % for m=1:M % % for n=1:N % % fi(m,n)=20*exp(-((m-M/2).^2+(n-N/2).^2)/150); % % end % % end % % % fi(x,y)=x*pi/10; % % fii=exp(j*fi); % % ph_wrap=angle(fii); % % % ph_wrap=mod((pi+U),2*pi)-pi;%包裹位相 % % figure,imshow(ph_wrap,[]);%包裹位相图 % % figure,surf(ph_wrap); % % grid off; % % axis off; % % shading interp; % % IM_mask=ones(M,N); % % [residue_charge,residue_sum0]=PhaseResidues(P,IM_mask); % % P=ph_wrap; % % % [unwrapped]=CLUnwrap(P); % % unwrapped=phaseunwrapping(P); % % figure,imshow(unwrapped,[]); % % figure,surf(unwrapped); % % grid on; % % axis on; % % shading interp; % % cc=cos(unwrapped); % % figure,imshow(cc); % % p=exp(j.*unwrapped); % % pp1=angle(p); % % figure,imshow(pp1,[]); % % Mask=ones(M,N); % % UnwrappedFase=LingXingUnwrap2(P,Mask,2,2); % % figure,imshow(UnwrappedFase,[]); % % figure,surf(UnwrappedFase); % % shading interp; % % axis on; % % grid on; % % CC=cos(UnwrappedFase); % % figure,imshow(CC); % clear % close all % O1=imread('9.bmp');%含有样品的像面全息图 % O1=im2double(O1); % % O1=rgb2gray(O1); % figure,imshow(O1); % % O2=fftshift(fft2(ifftshift(O1)));%频谱 % figure,imshow(1000*abs(O2)/max(max(abs(O2)))); % PIN=O2(152:369,223:510); %图9提取频谱中的正或负一级像 % [M,N]=size(PIN); % ZX=fftshift(ifft2(ifftshift(PIN)));%角谱再现 % % ZX1=imag(ZX); % % ZX2=real(ZX); % % ZX1=medfilt2(ZX1,[3,3]); % % ZX2=medfilt2(ZX2,[3,3]); % % ZNCP1=atan(ZX1./(ZX2+eps)); % % figure,imshow(ZNCP1,[]); % ZNCP1=angle(ZX); % P=ZNCP1; % U=ZX; % cc=cos(P); % figure,imshow(cc); %原始DCT法 % clear % close all % O3=imread('12.jpg');%无样品全息图 % % O3=im2double(O3); % O3=rgb2gray(O3); % % figure,imshow(O3); % % O4=fftshift(fft2(ifftshift(O3)));%频谱 % % % 位相重建 % % PIN2=O4(320:490,430:670); % PIN2=O4(395:455,825:915); % ZX2=fftshift(ifft2(ifftshift(PIN2))); % [M,N]=size(PIN2); % ZNCP1=angle(ZX2); % figure,imshow(ZNCP1,[]); % P=ZNCP1; % U=ZX2; % cc=cos(P); % figure,imshow(cc); [r,rsum1]=PhaseResidues(P); Mask=ones(M,N); for m=1:M for n=1:N if r(m,n)==1 | r(m,n)==-1 Mask(m,n)=0; end end end UnwrappedFase=LingXingUnwrap2(P,Mask,2,2); figure,imshow(UnwrappedFase,[]); figure,surf(UnwrappedFase); shading interp; axis on; grid on; dx=zeros(M,N); %预先给出X方向的梯度 dy=zeros(M,N); %预先给出y方向的梯度 m=1:M-1; %通过循环计算梯度 dx(m,:)=P(m+1,:)-P(m,:); max(max(dx)) dx=dx-pi*floor(dx/pi+0.5); %将截断的梯度连接起来 n=1:N-1; dy(:,n)=P(:,n+1)-P(:,n); max(max(dy)) dy=dy-pi*floor(dy/pi+0.5); p=zeros(M,N); %求ρ(x,y) p1=zeros(M,N); p2=zeros(M,N); m=2:M; p1(m,:)=dx(m,:)-dx(m-1,:); n=2:N; p2(:,n)=dy(:,n)-dy(:,n-1); p=p1+p2; %求ρ(x,y) p(1,1)=dx(1,1)+dy(1,1); %边界条件 n=2:N; p(1,n)=dx(1,n)+dy(1,n)-dy(1,n-1); %边界条件 m=2:M; p(m,1)=dx(m,1)-dx(m-1,1)+dy(m,1); %边界条件 m=2:M; w1(m,:)=p1(m,:)-p1(m-1,:); m=2:M; w2(m,:)=p2(m,:)-p2(m-1,:); wij=sqrt(w1.^2+w2.^2); k=0.01; p=(1+k*wij).*p; pp=dct2(p)+eps; %对ρ(x,y)作余弦变换 P_unwrap=zeros(M,N); %求余弦变换后的物体本位Ψ for m=1:M for n=1:N P_unwrap(m,n)=pp(m,n)/(2*cos(pi*(m-1)/M)+2*cos(pi*(n-1)/N)-4+eps); end end P_unwrap(1,1)=pp(1,1); P_unwrap1=idct2(P_unwrap); %再作反余弦变换得到物体相位 figure,imshow(P_unwrap1,[]); figure,surf(P_unwrap1); grid on; axis on; shading interp; % IM_mask=ones(M,N); % [residue_charge,residue_sum1]=PhaseResidues(P_unwrap1,IM_mask); % PP1=cos(P_unwrap1); % figure,imshow(PP1); % Q=unwrap2xin(P); % figure,imshow(Q,[]); % figure,surf(Q); % grid on; % axis on; % shading interp; % PP2=cos(Q); % figure,imshow(PP2); % Q2=unwrap22xin(P); % figure,imshow(Q2,[]); % figure,surf(Q2); % grid on; % axis on; % shading interp; % PP3=cos(Q1); % figure,imshow(PP3); %基于四向的DCT法 dx=zeros(M,N); %预定义X方向的梯度 dy=zeros(M,N); %预定义y方向的梯度 du=zeros(M,N); %预定义左下对角线方向的梯度 dv=zeros(M,N); %预定义右下对角线方向的梯度 m=1:M-1; %通过循环计算梯度 dx(m,:)=P(m+1,:)-P(m,:); dx=dx-pi*floor(dx/pi+0.5); %将截断的梯度连接起来 n=1:N-1; dy(:,n)=P(:,n+1)-P(:,n); dy=dy-pi*floor(dy/pi+0.5); m=2:M-1; n=2:N-1; dv(m,n)=P(m+1,n-1)-P(m,n); n=2:N; dv(1,n)=P(2,n-1)-P(1,n); dv=dv-pi*floor(dv/pi+0.5); m=1:M-1; n=1:N-1; du(m,n)=P(m+1,n+1)-P(m,n); du=du-pi*floor(du/pi+0.5); p=zeros(M,N); %求ρ(x,y) p1=zeros(M,N); p2=zeros(M,N); p3=zeros(M,N); p4=zeros(M,N); m=2:M-1; p1(m,:)=dx(m,:)-dx(m-1,:); n=2:N-1; p2(:,n)=dy(:,n)-dy(:,n-1); m=2:M-1; n=2:N-1; p3(m,n)=dv(m,n)-dv(m-1,n+1); m=2:M-1; n=2:N-1; p4(m,n)=du(m,n)-du(m-1,n-1); p=p1+p2+p3+p4; %求ρ(x,y) p(1,1)=dx(1,1)+dy(1,1)+du(1,1); %边界条件 n=2:N; p(1,n)=dx(1,n)+dy(1,n)-dy(1,n-1)+du(1,n)+dv(1,n); %边界条件 n=2:N-1; p(M,n)=-dx(M-1,n)+dy(M,n)-dy(M,n-1)-dv(M-1,n+1)-du(M-1,n-1); %边界条件 m=2:M; p(m,1)=dx(m,1)-dx(m-1,1)+dy(m,1)-dv(m-1,2)+du(m,1); %边界条件 m=2:M; p(m,N)=dx(m,N)-dx(m-1,N)-dy(m,N-1)+dv(m,N)-du(m-1,N-1); m=2:M; w1(m,:)=p1(m,:)-p1(m-1,:); m=2:M; w2(m,:)=p2(m,:)-p2(m-1,:); m=2:M; w3(m,:)=p3(m,:)-p3(m-1,:); m=2:M; w4(m,:)=p4(m,:)-p4(m-1,:); wij=sqrt(w1.^2+w2.^2+w3.^2+w4.^2); k=0.01; p=(1+k*wij).*p; pp=dct2(p)+eps; %对ρ(x,y)作余弦变换 P_unwrap=zeros(M,N); %求余弦变换后的物体本位Ψ for m=1:M for n=1:N P_unwrap(m,n)=pp(m,n)/(2*cos(pi*(m-1)/M)+2*cos(pi*(n-1)/N)+4*cos(pi*(m-1)/M)*cos(pi*(n-1)/N)-8+eps); end end P_unwrap(1,1)=pp(1,1); P_unwrap2=idct2(P_unwrap); %再作反余弦变换得到物体相位 figure,imshow(P_unwrap2,[]); figure,surf(P_unwrap2); grid on; axis on; shading interp; cc=cos(P_unwrap2); figure,imshow(cc); IM_mask=ones(M,N); % [residue_charge,residue_sum2]=PhaseResidues(P_unwrap2,IM_mask); %基于横向剪切的最小二乘法 ZX=U; U1=zeros(M,N); %初始化横向剪切的再现光场 U2=zeros(M,N); %初始化纵向剪切的再现光场 UOX=zeros(M,N); UOY=zeros(M,N); m=1:M-1; U1(m,:)=ZX(m+1,:); n=1:N-1; U2(:,n)=ZX(:,n+1); UOX(m,n)=U1(m,n)./ZX(m,n); %横向的剪切后光场与原始光场相除 UOY(m,n)=U2(m,n)./ZX(m,n); %纵向的剪切后光场与原始光场相除 ph_wrap1=angle(UOX); %横向位相差 ph_wrap2=angle(UOY); %纵向位相差 %然后利用离散余弦最小二乘解包裹求得结果 p=zeros(M,N); %求ρ(x,y) p1=zeros(M,N); p2=zeros(M,N); m=2:M; p1(m,:)=ph_wrap1(m,:)-ph_wrap1(m-1,:); n=2:N; p2(:,n)=ph_wrap2(:,n)-ph_wrap2(:,n-1); p=p1+p2; %求ρ(x,y) p(1,1)=ph_wrap1(1,1)+ph_wrap2(1,1); %边界条件 n=2:N; p(1,n)=ph_wrap1(1,n)+ph_wrap2(1,n)-ph_wrap2(1,n-1); %边界条件 m=2:M; p(m,1)=ph_wrap1(m,1)-ph_wrap1(m-1,1)+ph_wrap2(m,1); %边界条件 m=2:M; w1(m,:)=ph_wrap1(m,:)-ph_wrap1(m-1,:); m=2:M; w2(m,:)=ph_wrap2(m,:)-ph_wrap2(m-1,:); wij=sqrt(w1.^2+w2.^2); k=0.01; p=(1+k*wij).*p; pp=dct2(p)+eps; %对ρ(x,y)作余弦变换 P_unwrap=zeros(M,N); %求余弦变换后的物体本位Ψ for m=1:M for n=1:N P_unwrap(m,n)=pp(m,n)/(2*cos(pi*(m-1)/M)+2*cos(pi*(n-1)/N)-4+eps); end end P_unwrap(1,1)=pp(1,1); P_unwrap4=idct2(P_unwrap); %再作反余弦变换得到物体相位 figure,imshow(P_unwrap4,[]) figure,surf(P_unwrap4); grid off; axis off; shading interp; %基于四向横向剪切的最小二乘的算法 % [M,N]=size(O1); %确定包裹相位矩阵的大小 U1=zeros(M,N); %初始化横向剪切的再现光场 U2=zeros(M,N); %初始化纵向剪切的再现光场 U3=zeros(M,N);%初始化左下角剪切的再现光场 U4=zeros(M,N);%初始化右下角剪切的再现光场 UOX=zeros(M,N); UOY=zeros(M,N); UOV=zeros(M,N); UOU=zeros(M,N); m=1:M-1; U1(m,:)=U(m+1,:); n=1: