多频外差（三频四步相移）结构光三维重建——Matlab仿真

最近想通过三频四步做一下效果对比，发现网上找到的代码不够明了，特此写的仿真代码，虽然也有读图版本，但是个人觉得只需清楚原理即可。

具体数学公式参考李中伟博士论文：

外差原理是指将两种不同频率的相位函数叠加得到一种频率更低相位函数。具体计算我就不再赘述，直接上代码具体问题具体分析。

close all;
clc;
clear;

%% 参数
W = 1000; 
H = 1000;
A = 127;
B = 127;

[x, y] = meshgrid(1:W, 1:H);
def=peaks(W)*10;
figure;
mesh(x,y,def);
title('def');
%%
T1 = 44; 
T2 = 48;
T3 = 52;
T12=T1*T2/(T2-T1);
T23=T2*T3/(T3-T2);
T123=T23*T12/(T23-T12);

%%
I1 = cell(1,4);
I1_r=cell(1,4);
for N=1:4
  for i=1:W
    for j=1:H
        I1{N}(i, j) = A+B*cos((2*pi*j/T1)+((N-1)*2*pi/4));
        I1_r{N}(i, j) = A+B*cos((2*pi*(j+def(i, j))/T1)+((N-1)*2*pi/4));
    end
  end
end

I2 = cell(1,4);
I2_r=cell(1,4);
for N=1:4
  for i=1:W
    for j=1:H
        I2{N}(i, j) = A+B*cos((2*pi*j/T2)+((N-1)*2*pi/4));
        I2_r{N}(i, j) = A+B*cos((2*pi*(j+def(i, j))/T2)+((N-1)*2*pi/4));
    end
  end
end

I3 = cell(1,4);
I3_r=cell(1,4);
for N=1:4
  for i=1:W
    for j=1:H
        I3{N}(i, j) = A+B*cos((2*pi*j/T3)+((N-1)*2*pi/4));
        I3_r{N}(i, j) = A+B*cos((2*pi*(j+def(i, j))/T3)+((N-1)*2*pi/4));
    end
  end
end

% figure;
% imshow(I3_r{4},[]);
%%
%% 重命名
for i = 1:4
    name = sprintf('I1_%d', i);
    assignin('base', name, I1{i});
    name = sprintf('Ir1_%d', i);
    assignin('base', name, I1_r{i});

    name = sprintf('I2_%d', i);
    assignin('base', name, I2{i});
    name = sprintf('Ir2_%d', i);
    assignin('base', name, I2_r{i});

    name = sprintf('I3_%d', i);
    assignin('base', name, I3{i});
    name = sprintf('Ir3_%d', i);
    assignin('base', name, I3_r{i});
end
%% 求包裹相位
phi = cell(1, 3); 
for j = 1:3

    phi{j} = -atan2(eval(['I', num2str(j), '_1']) - eval(['I', num2str(j), '_3']), ...
                    eval(['I', num2str(j), '_4']) - eval(['I', num2str(j), '_2']));

    % phi{j}(phi{j} < 0) = phi{j}(phi{j} < 0) + 2 * pi;
   
    figure('name',sprintf('包裹相位-%d', j)), imshow(phi{j});
    % figure('name',sprintf('相位截面-%d', j)),plot((phi{j}(600,:)));
end
pha = cell(1, 3); 
for j = 1:3

    pha{j} = -atan2(eval(['Ir', num2str(j), '_1']) - eval(['Ir', num2str(j), '_3']), ...
                    eval(['Ir', num2str(j), '_4']) - eval(['Ir', num2str(j), '_2']));

    % pha{j}(pha{j} < 0) = pha{j}(pha{j} < 0) + 2 * pi;
   
    figure('name',sprintf('调制包裹相位-%d', j)), imshow(pha{j});
    % figure('name',sprintf('调制相位截面-%d', j)),plot((pha{j}(600,:)));
end
%% 相位差
% ******原******
P_12=zeros(H,W);
P_23=zeros(H,W);
P_123=zeros(H,W);

for i=1:W    
    for j=1:H
        if(phi{1}(i,j)>=phi{2}(i,j))
            P_12(i,j)=phi{1}(i,j)-phi{2}(i,j);
        else
            P_12(i,j)=phi{1}(i,j)-phi{2}(i,j)+2*pi;
        end

        if(phi{2}(i,j)>=phi{3}(i,j))
            P_23(i,j)=phi{2}(i,j)-phi{3}(i,j);
        else
            P_23(i,j)=phi{2}(i,j)-phi{3}(i,j)+2*pi;
        end


        if(T23<T12)
            if(P_12(i,j)>=P_23(i,j))
               P_123(i,j)=P_12(i,j)-P_23(i,j)+2*pi;
            else
               P_123(i,j)=P_12(i,j)-P_23(i,j);
            end
        else
            if(P_12(i,j)>=P_23(i,j))
               P_123(i,j)=P_12(i,j)-P_23(i,j);
            else
               P_123(i,j)=P_12(i,j)-P_23(i,j)+2*pi;
            end
        end

        
    end
end
 figure('name','相位差图P12'), imshow(P_12,[]);
 % figure('name','相位差图-P12'),plot((P_12(H/2,:)));

 figure('name','相位差图P23'), imshow(P_23,[]);
 % figure('name','相位差图-P23'),plot((P_23(H/2,:)));

 figure('name','相位差图P123'), imshow(P_123,[]);
 % figure('name','相位差图-P123'),plot((P_123(H/2,:)));

% ******调制******
Pa_12=zeros(H,W);
Pa_23=zeros(H,W);
Pa_123=zeros(H,W);

for i=1:W    
    for j=1:H

        if(pha{1}(i,j)>=pha{2}(i,j))
            Pa_12(i,j)=pha{1}(i,j)-pha{2}(i,j);
        else
            Pa_12(i,j)=pha{1}(i,j)-pha{2}(i,j)+2*pi;
        end

        if(pha{2}(i,j)>=pha{3}(i,j))
            Pa_23(i,j)=pha{2}(i,j)-pha{3}(i,j);
        else
            Pa_23(i,j)=pha{2}(i,j)-pha{3}(i,j)+2*pi;
        end
        
        if(T23<T12)
            if(Pa_12(i,j)>=Pa_23(i,j))
               Pa_123(i,j)=Pa_12(i,j)-Pa_23(i,j)+2*pi;
            else
               Pa_123(i,j)=Pa_12(i,j)-Pa_23(i,j);
            end
        else
            if(Pa_12(i,j)>=Pa_23(i,j))
               Pa_123(i,j)=Pa_12(i,j)-Pa_23(i,j);
            else
               Pa_123(i,j)=Pa_12(i,j)-Pa_23(i,j)+2*pi;
            end
        end

    end
end
 figure('name','调制相位差图Pa12'), imshow(Pa_12,[]);
 % figure('name','调制相位差图-Pa12'),plot((Pa_12(H/2,:)));

 figure('name','调制相位差图Pa23'), imshow(Pa_23,[]);
 % figure('name','调制相位差图-Pa23'),plot((Pa_23(H/2,:)));

 figure('name','调制相位差图Pa123'), imshow(Pa_123,[]);
 % figure('name','调制相位差图-Pa123'),plot((Pa_123(H/2,:)));

%% 相位展开
k_12=round(((P_123*T123/T12)-P_12)/(2*pi))*(2*pi);
k=round(((P_12+k_12)*T12/T1-phi{1})/(2*pi))+1;

UP_1=phi{1}+2*pi*k;

figure;
yyaxis left
plot(k(H/2,:));
ylabel('k'); % 设置左侧纵坐标轴标签
yyaxis right
plot(phi{1}(H/2,:));
ylabel('phi1'); % 设置右侧纵坐标轴标签
title('原k_ and 包裹phi{1}');
figure('name','展开相位'), imshow(UP_1,[]);

ka_12=round(((Pa_123*T123/T12)-Pa_12)/(2*pi))*(2*pi);
ka=round(((Pa_12+ka_12)*T12/T1-pha{1})/(2*pi))+1;

UPa_1=pha{1}+2*pi*ka;

figure;
yyaxis left
plot(ka(H/2,:));
ylabel('ka'); % 设置左侧纵坐标轴标签
yyaxis right
plot(pha{1}(H/2,:));
ylabel('pha1'); % 设置右侧纵坐标轴标签
title('调制k_ and 包裹pha{1}');
figure('name','展开相位'), imshow(UPa_1,[]);

%%
phax=UPa_1-UP_1;
figure;
imshow(phax, []);
title('绝对相位差');

phaxx=phax*T1/(2*pi);
figure;
mesh(x, y, phaxx);
title('三维重建');

值得注意的是，我的条纹图的设计是：

for N=1:4
  for i=1:W
    for j=1:H
        I1{N}(i, j) = A+B*cos((2*pi*j/T1)+((N-1)*2*pi/4));
        I1_r{N}(i, j) = A+B*cos((2*pi*(j+def(i, j))/T1)+((N-1)*2*pi/4));
    end
  end
end

 其中相位计算是2*pi*j/T

而变形条纹中的调制高度信息必须在（j+def(i,j)）

其中T可以理解位一个2 周期长度。因为外差最后的一个周期长度T123要大于图像的整体长度。

求出外差相位后，解包裹阶段只需要有由123和12向下解两次，求得1的阶次k即可。

另外，在进行外差相位计算时，要加入判断前后相位的大小。

        if(phi{1}(i,j)>=phi{2}(i,j))
            P_12(i,j)=phi{1}(i,j)-phi{2}(i,j);
        else
            P_12(i,j)=phi{1}(i,j)-phi{2}(i,j)+2*pi;
        end

直接上结果：

待重建物体

调制后的包裹相位：

外差相位：

求解的k：

相位展开：

三维重建的结果：