结构光编解码-多频外差代码

	多频外差原理可查看之前的文章《结构光编解码-多频外差》，以下代码为多频外差解相位的代码，该代码需修改图像路径及图像格式，使用的频率也需做相应的替换，替换完后即可正常运行。（欢迎进群交流：874653199）

clc;
clear;

%% 图像路径及图像格式
s1 = 'D:\2_光学测量\6_数据\三频四相\';
s2 = '.bmp';
%% 获取图像宽高，并创建数组用于存储相移图
filename = [s1 '1' s2];
[row,col] = size(imread(filename));
srcImgNums = 12;
images = uint8(zeros(row,col,srcImgNums));
%% 读取相移数据
for i = 1:srcImgNums 
    filename = [s1 num2str(i) s2];
    images(:,:,i) = imread(filename); 
end

%% 解算每个频率的折叠相位
freqImg = uint8(zeros(row,col,4));%uint8
freqWrappedPhase = zeros(row,col,3);
freq = [100 99 90]; 

for i = 1:3
    freqImg = images(:,:,(i-1)*4+1:(i-1)*4+4);
    fenzi = double(freqImg(:,:,4))-double(freqImg(:,:,2));
    fenmu = double(freqImg(:,:,1))-double(freqImg(:,:,3));
    temp_wrapped_phase = atan(fenzi./fenmu);
    
    flagMat = any(fenmu==0,3).*any(fenzi<0,3);
    temp_wrapped_phase = temp_wrapped_phase + 2*pi*flagMat; 
    flagMat = any(fenmu<0,3);
    temp_wrapped_phase = temp_wrapped_phase + pi*flagMat;
    flagMat = any(fenmu>0,3).*any(fenzi<0,3);
    temp_wrapped_phase = temp_wrapped_phase + 2*pi*flagMat;
    
    freqWrappedPhase(:,:,i) = temp_wrapped_phase;
    figure
    hold on
    imshow(freqWrappedPhase(:,:,i),[]);
    imgName=['频率' num2str(freq(i)) '的折叠相位'];
    title(imgName,'color','r')
    hold off
end


%% 计算绝对相位
absPhase = zeros(row,col);%未校正畸变的绝对相位图

[phase_abs13,freq13] = heterodynePhase(freqWrappedPhase(:,:,1),freq(1),freqWrappedPhase(:,:,3),freq(3));
[phase_abs23,freq23] = heterodynePhase(freqWrappedPhase(:,:,2),freq(2),freqWrappedPhase(:,:,3),freq(3));
[phase_abs_abs,freq1323] =heterodynePhase(phase_abs13,freq13,phase_abs23,freq23);
absPhase1 = phase_abs_abs*freq(1);
absPhase2 = round((phase_abs_abs*freq(2)-freqWrappedPhase(:,:,2))/2/pi)*2*pi+freqWrappedPhase(:,:,2);
absPhase3 = round((phase_abs_abs*freq(3)-freqWrappedPhase(:,:,3))/2/pi)*2*pi+freqWrappedPhase(:,:,3);
absPhase = (absPhase1+absPhase2+absPhase3)/sum(freq);

figure
hold on
imshow(absPhase,[]);
title('绝对相位','color','r')
hold off

%% 双频外差函数
function [phaseij_abs, freqij] = heterodynePhase(phasei, freqi, phasej, freqj)% 要求Nfi大于Nfj
freqij = freqi - freqj;
k = freqi/freqij;
phaseij = phasei - phasej;
index = find(phaseij<0);
phaseij(index) = phaseij(index) + 2*pi;
period = round((k*phaseij-phasei)/2/pi);
phaseij_abs = (period*2*pi + phasei)/k;

end

折叠相位展开结果如下：