多频法(倍频法)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}))+pi; % f=1,加两个负号是为了移动初相
m_wrapPhase = atan2(-(I{8,1}-I{6,1}),-(I{5,1}-I{7,1}))+pi; % f= 8 
h_wrapPhase = atan2(-(I{12,1}-I{10,1}),-(I{9,1}-I{11,1}))+pi; % f = 64

% 可以先求出三个调制度阈值后单独使用
% l_wrapPhase(B1<th)=nan;
% m_wrapPhase(B2<th)=nan;
% h_wrapPhase(B3<th)=nan;

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)
%% 多频法解码
% 条纹频率，需要产生时的条纹频率一一对应
f1 = 1;
f2 = 8;
f3 = 64;
% f = 8 的解包相位
m_unwrapPhase = 2*round((f2./f1 *l_wrapPhase-m_wrapPhase)./(2*pi))*pi+m_wrapPhase;
% f = 64的解包相位
k = round((f3./f2 *m_unwrapPhase-h_wrapPhase)./(2*pi)); % 条纹级次
unwrapPhase = 2*pi*k+h_wrapPhase;
end

generate.m

function generate()
%% 产生条纹，可以直接给DLP4500
W = 912; % 投影仪的分辨率
H = 1140;

f1 = 1; % 条纹频率
f2 =8;
f3 = 64;

p1 = W./f1;
p2 = W./f2;
p3 = W./f3;

[X,~] = meshgrid(1:W,1:H);

a = 0.5;
b = 0.5;

I1 = a+b*cos(2*pi*1./p1*X);
I2 = a+b*cos(2*pi*1./p1*X+pi/2);
I3 = a+b*cos(2*pi*1./p1*X+pi);
I4 = a+b*cos(2*pi*1./p1*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./p2*X);
I6 = a+b*cos(2*pi*1./p2*X+pi/2);
I7 = a+b*cos(2*pi*1./p2*X+pi);
I8 = a+b*cos(2*pi*1./p2*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./p3*X);
I10 = a+b*cos(2*pi*1./p3*X+pi/2);
I11 = a+b*cos(2*pi*1./p3*X+pi);
I12 = a+b*cos(2*pi*1./p3*X+3*pi/2);

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

end

附上百度网盘实际投影仪采集数据与模拟数据:

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