三维重建之多频外差解包裹学习笔记

Posted 2023-02-28 u012507022

 

附matlab多频外差解相位程序

---

%参考博客:https://blog.csdn.net/qq_15295565/article/details/99704919，并进行修改
clc;
close all;
clear;
 
%%  初始化
width  = 1024;  
heigth = 768;

 %三频四步
FREQ = 3;
STEP  = 4;

%freq为三个正弦条纹对应的周期数，λ=width/freq:代表每个条纹的周期(波长)
% 也可理解为频率，三个频率参考李中伟的博士论文
freq = [70 64 59];  
freq12    = 6 ;
freq23    = 5 ;  
freq123  = 1 ;
%% Step0 生成条纹图
% 利用分块矩阵C存储3组共计12张图，三种频率，四步相位
C = cell(FREQ,STEP);  
for i=1:FREQ
    for j=1:STEP
        Ci,j = zeros(heigth,width);
    end
end
 
% 利用余弦函数计算12张图的灰度值，图像的生成
% 三种频率，四组相位
for i = 1:FREQ       % 对应三种频率
    for  j = 0:(STEP) % 对应四步相位
        for k = 1:width
            Ci,j+1(:,k) = 128+127*sin(2*pi*k*freq(i)/width+j*pi/2);
        end
    end
end
 
% 归一化处理
for i = 1:3
    for j = 1:4
        Ci,j = mat2gray(Ci,j);
    end
end
 
% 显示12张图
% for i = 1:FREQ
%     for j = 1:STEP
%         n = STEP*(i-1)+j;
%         h = figure(n);
%         imshow(Ci,j);
%     end
% end
 
%% Step1 求解相位主值

% phi存储相位主值图像
phi = cell(FREQ,1);
for i = 1:FREQ
    phii,1 = zeros(heigth,width);
end
 
% 计算每种频率对应的相位主值
% 输出三种频率的相位主值，用于相差计算
for i = 1:FREQ % 对于3组中的每一组图片，每一组相同频率的有四张图片
     I1 = Ci,1;
     I2 = Ci,2;
     I3 = Ci,3;
     I4 = Ci,4;
      for g = 1:heigth
          for k = 1:width          
            if     I4(g,k)==I2(g,k)&&I1(g,k)>I3(g,k)   %四个特殊位置
                    phii,1(g,k)=0;
            elseif I4(g,k)==I2(g,k)&&I1(g,k)<I3(g,k) %四个特殊位置
                    phii,1(g,k)=pi; 
            elseif I1(g,k)==I3(g,k)&&I4(g,k)>I2(g,k) %四个特殊位置
                    phii,1(g,k)=pi/2;
            elseif I1(g,k)==I3(g,k)&&I4(g,k)<I2(g,k) %四个特殊位置
                    phii,1(g,k)=3*pi/2;
            elseif I1(g,k)<I3(g,k) %二三象限
                    phii,1(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)))+pi;
            elseif I1(g,k)>I3(g,k)&&I4(g,k)>I2(g,k) %第一象限
                    phii,1(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)));
            elseif I1(g,k)>I3(g,k)&&I4(g,k)<I2(g,k) %第四象限
                    phii,1(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)))+2*pi;  
            end %endif
            
          end             
     end
end
 
% 显示
% figure,imshow(mat2gray(PH1));title('1相位主值'); 
% figure,imshow(mat2gray(PH2));title('2相位主值'); 
% figure,imshow(mat2gray(PH3));title('3相位主值'); 

%% Step2 利用外差原理叠加不同频率的相位主值，得到整幅图像只有一个周期的相位
PH1 = phi1,1;   %频率1
PH2 = phi2,1;   %频率2
PH3 = phi3,1;   %频率3

PH12   = zeros(heigth,width);
PH23   = zeros(heigth,width);
PH123 = zeros(heigth,width);
 
for g = 1:heigth
	for k = 1:width
		% 计算第1组和第2组的相差
		if PH1(g,k)>PH2(g,k)
			PH12(g,k) = PH1(g,k)-PH2(g,k);
		else
			PH12(g,k) = PH1(g,k)+2*pi-PH2(g,k);
        end
        
		% 计算第2组和第3组的相差
		if PH2(g,k)>PH3(g,k)
			PH23(g,k) = PH2(g,k)-PH3(g,k);
		else
			PH23(g,k) = PH2(g,k)+2*pi-PH3(g,k);
        end
        
        % 计算第12组和第23组的相差
        if PH12(g,k)>PH23(g,k)
			PH123(g,k) = PH12(g,k)-PH23(g,k);
		else
			PH123(g,k) = PH12(g,k)+2*pi-PH23(g,k);
        end
        
	end
end

%显示
figure,imshow(mat2gray(PH12));title('1,2外差');  
figure,imshow(mat2gray(PH23));title('2,3外差'); 
figure,imshow(mat2gray(PH123));title('1,2,3外差');

%% Step3 相位解包裹；
%根据公式:
%η = floor(((delta_φ*delta_λ/λ)-φ)/2π)，
%Φ = 2πη+φ ，
%由PH123、PH23、PH12反计算出PH1、PH2、PH3的连续相位

%初始化解包裹相位
phUnWrap = cell(3,1);
for i = 1:3
    phUnWrapi,1 = zeros(heigth,width);
end

% 解包裹23：求PH23的绝对相位，deltaPh = PH123，phWrap = PH23
deltaPh = PH123;
phWrap = PH23;
R = freq23/freq123;

Nwrap = floor((deltaPh*R-phWrap)/(2*pi)+0.5);
ph23UnWrap= (2*pi)*Nwrap +phWrap;

% 解包裹12：求PH12的绝对相位，deltaPh = PH123，phWrap = PH12
deltaPh = PH123;
phWrap = PH12;
R = freq12/freq123;

Nwrap = floor((deltaPh*R-phWrap)/(2*pi)+0.5);
ph12UnWrap= (2*pi)*Nwrap +phWrap;

% 解包裹1：求PH1的绝对相位，deltaPh = PH12，phWrap = PH1
deltaPh = ph12UnWrap;
phWrap = PH1;
R = freq(1)/freq12;

Nwrap = floor((deltaPh*R-phWrap)/(2*pi)+0.5);
phUnWrap1,1= (2*pi)*Nwrap +phWrap;

% 解包裹2：求PH2的绝对相位，deltaPh = PH12，phWrap = PH2
deltaPh = ph12UnWrap;
phWrap = PH2;
R = freq(2)/freq12;

Nwrap = floor((deltaPh*R-phWrap)/(2*pi)+0.5);
phUnWrap2,1= (2*pi)*Nwrap +phWrap;

% 解包裹3：求PH2的绝对相位，deltaPh = PH23，phWrap = PH3
deltaPh = ph23UnWrap;
phWrap = PH3;
R = freq(3)/freq23;

Nwrap = floor((deltaPh*R-phWrap)/(2*pi)+0.5);
phUnWrap3,1= (2*pi)*Nwrap +phWrap;

%% show phUnWrap
figure
plot(PH123(100,:),'r');hold on
plot(PH23(100,:),'g');hold on
plot(ph23UnWrap(100,:));
grid on
legend('deltaPh','wrapPH','phUnWrap');
title('解包裹23'); 

figure
plot(PH123(100,:),'r');hold on
plot(PH12(100,:),'g');hold on
plot(ph12UnWrap(100,:));
grid on
legend('deltaPh','wrapPH','phUnWrap');
title('解包裹12'); 

figure
plot(ph12UnWrap(100,:),'r');hold on
plot(PH1(100,:),'g');hold on
plot(phUnWrap1,1(100,:));
grid on
legend('deltaPh','wrapPH','phUnWrap');
title('解包裹1'); 

figure
plot(ph12UnWrap(100,:),'r');hold on
plot(PH2(100,:),'g');hold on
plot(phUnWrap2,1(100,:));
grid on
legend('deltaPh','wrapPH','phUnWrap');
title('解包裹2'); 

figure
plot(ph23UnWrap(100,:),'r');hold on
plot(PH3(100,:),'g');hold on
plot(phUnWrap3,1(100,:));
grid on
legend('deltaPh','wrapPH','phUnWrap');
title('解包裹3'); 

参考：

https://blog.csdn.net/qq_15295565/article/details/99704919
https://blog.csdn.net/gao_summer_cola/article/details/75308618
https://mp.weixin.qq.com/s/eTlkyNrqDk_qrb7eS91ktQ
https://zhuanlan.zhihu.com/p/82521153