三频外差法Matlab代码

main.m 

clc,clear,close all
%% 主要参数
N = 4; % 4步相移
isfilter = 0; % 是否高速滤波，可以降低噪音，条纹质量不好可以使用
sig = 0.5; % sigma
I = cell(12,1);
%% 数据准备
fileDir = "模拟数据"; % 存储图片的路径，低频条纹命名序号要小于高频条纹
fileDS=imageDatastore(fileDir);
nums = size(fileDS.Files);
th = 0.1;
for n = 1:nums(1)
    I{n,1}=im2double(imread(fileDS.Files{n}));
    if isfilter==1
        I{n,1}=imgaussfilt(I{n,1},sig);
    end
end

%% 求调制度，具体可以参考相关公式。可以使用通用公式
B = ((I{4,1}-I{2,1}).^2+(I{1,1}-I{3,1}).^2).^0.5/2;
%% 4步相移法求相位
l_wrapPhase = atan2(I{4,1}-I{2,1},I{1,1}-I{3,1});
m_wrapPhase = atan2(I{8,1}-I{6,1},I{5,1}-I{7,1});
h_wrapPhase = atan2(I{12,1}-I{10,1},I{9,1}-I{11,1});

%% 解相
unwrapPhase=decode(l_wrapPhase,m_wrapPhase,h_wrapPhase);
unwrapPhase(B<th)=nan;
imagesc(unwrapPhase),colormap(gray(256))

decode.m代码 

function unwrapPhase=decode(l_wrapPhase,m_wrapPhase,h_wrapPhase)
%%  三频外差法,总共需要投影12张条纹
% l_wrapPhase ：低频包裹相位
% m_wrapPhase ：中频包裹相位
% h_wrapPhase ：高频包裹相位
% 输出：unwrapPhase 高频解包相位

W =912;  % 投影的条纹宽度，本文垂直投影需要使用宽度整个参数。使用宽还是高需要根据条纹编码时使用的参数决定
H =1140; % 未使用

f=[59 64 70]; % 周期是width/f，频率就是912列（相对投影仪的分辨率）中有多少个条纹

wrapPhase1=l_wrapPhase; % 低频
wrapPhase2=m_wrapPhase; % 中频
wrapPhase3=h_wrapPhase; % 高频

% 两两外差
diff12 = wrapPhase2-wrapPhase1;
diff23 =wrapPhase3-wrapPhase2;

diff12(wrapPhase2<wrapPhase1)=diff12(wrapPhase2<wrapPhase1)+2*pi;
diff23(wrapPhase3<wrapPhase2)=diff23(wrapPhase3<wrapPhase2)+2*pi;


diff123 = diff23-diff12;
diff123(diff23<diff12)=diff123(diff23<diff12)+2*pi;

% imagesc(diff12),colormap(gray);
p = H./f;
P12 =p(2)*p(1)/ abs(p(2)-p(1)); 
P23 =p(2)*p(3)/abs(p(3)-p(2));
P123 = P12*P23/abs(P12-P23);

% 逐层求解可以降低放大倍数对外差相位的影响,分母约大外差相位误会将会越大
unwrapPhase12 = 2*round((P123/P12*diff123-diff12)./(2*pi))*pi+diff12;
unwrapPhase23 = 2*round((P12/P23*unwrapPhase12-diff23)/(2*pi))*pi+diff23;

% 求高频解包相位
unwrapPhase = 2*round((P23/p(3)*unwrapPhase23-wrapPhase3)/(2*pi))*pi+wrapPhase3;
% plot(0:length(unwrapPhase(1000,:))-1,unwrapPhase(1000,:),'r');
% imagesc(unwrapPhase),colormap(gray(256))
% mesh(unwrapPhase)
end

generate.m投影条纹生成代码 

%% 产生条纹
W = 912; % 投影仪的分辨率，使用的是DLP4500 
H = 1140;
f = [59 64 70]; % 条纹的频率，和传统正余弦频率有所区别
P = W./f; 
a = 0.5;
b = 0.5;

%% 垂直投影 使用912 
[X,Y]=meshgrid(1:W,1:H);
X = X-1;
Y = Y-1;
I1 =a+b*cos(2*pi*1/P(1)*X);  
I2 =a+b*cos(2*pi*1/P(1)*X+pi/2);  
I3 =a+b*cos(2*pi*1/P(1)*X+pi);  
I4 =a+b*cos(2*pi*1/P(1)*X+3*pi/2);  

imwrite(I1,"I00.bmp");
imwrite(I2,"I01.bmp");
imwrite(I3,"I02.bmp");
imwrite(I4,"I03.bmp");

I5 =a+b*cos(2*pi*1/P(2)*X);  
I6 =a+b*cos(2*pi*1/P(2)*X+pi/2);  
I7 =a+b*cos(2*pi*1/P(2)*X+pi);  
I8 =a+b*cos(2*pi*1/P(2)*X+3*pi/2);  

imwrite(I5,"I04.bmp");
imwrite(I6,"I05.bmp");
imwrite(I7,"I06.bmp");
imwrite(I8,"I07.bmp");

I9 =a+b*cos(2*pi*1/P(3)*X);  
I10 =a+b*cos(2*pi*1/P(3)*X+pi/2);  
I11 =a+b*cos(2*pi*1/P(3)*X+pi);  
I12 =a+b*cos(2*pi*1/P(3)*X+3*pi/2);  

imwrite(I9,"I08.bmp");
imwrite(I10,"I09.bmp");
imwrite(I11,"I10.bmp");
imwrite(I12,"I11.bmp");

附上百度网盘链接,有实际的投影仪采集数据

        链接：https://pan.baidu.com/s/1CKYSuSSsrvAFu5rek5J14A?pwd=1234 
提取码：1234