多频外差法之三频四相的理论与实现（matlab）

在提取相位时，如果图像内的周期数仅为一个，则相对相位就是绝对相位，则求得的相位主值就是它的相位值。

实际中，我们投射的光栅并不是一个周期，那么在整个测量的空间中就会求得多个相同的相位主值，这时就需要对包裹的相位值进行展开。

相位包裹又是什么意思呢？

相位包裹的主要原因是相移法解相时使用了反正切函数，用atan2函数，得到四个象限的反正切，所以计算的相位都是在(-pi,pi]之间，也就是被包裹。

因而真实的相位值和被包裹的相位主值之间是相差 2npi的，其中 n 为整数,即为

多频外差，字面上意思就是多个频率，利用外差原理。

论文里看到的解释是这样的：

（1）多频外差编码法利用拍频原理，将两个或两个以上周期相近但不相同的相移编码光组合成一组编码光，把解得的不同周期的多个相位作差，将包裹相位的小周期放大为相位差的大周期，直到相位差信号的周期包含整个测量视场。

 这种方法中是用周期表示，也就是条纹的宽度，而且（T2>T1）。

（2） 李中伟博士的博士论文里是这样的：

如果对应的是频率，这里的值应该是小于1的值，即正弦波的个数除以整个图像方向上的像素值，它的周期才会是条纹的宽度，但在图中所对应的并非这种含义，所以我认为这里解释和图示是有矛盾的。

正确图示可能是这样的：

从双频外差可以看出是将小的周期扩展到了一个大的周期内，如果频率差为1可以达到全场范围内只有一个周期的相位的目的，展开仅与频率相位有关，不受颜色的影响。

% 程序开始
clc;
close all;
clear;% 图片的初始化
width = 1024;  
heigth = 768;% 三频率
% 这个可以参见李中伟的博士论文
freq = [70 64 59];  %像素单位为个数，可以看做频率；正弦函数为周期含义% 利用分块矩阵C存储3组共计12张图
% 三种频率，四组相位
C = cell(3,4);  
for i=1:3for j=1:4C{i,j} = zeros(heigth,width);end
end% 利用余弦函数计算12张图的灰度值
% 图像的生成
% 三种频率，四组相位
for i = 1:3 % 对应三种不同的频率for  j = 0:3 % 对应四种相位for k = 1:widthC{i,j+1}(:,k) = 128+127*sin(2*pi*k*freq(i)/width+j*pi/2);endend
end% 对灰度值进行归一化处理
for i = 1:3for j = 1:4C{i,j} = mat2gray(C{i,j});end
end% 显示12张图
% for i = 1:3
%     for j = 1:4
%         n = 4*(i-1)+j;
%         h = figure(n);
%         imshow(C{i,j});
%     end
% end% 初始化三组处理后的图片灰度矩阵
% phi也是分块矩阵
% 存储相位主值图像
phi = cell(3,1);
for i = 1:3phi{i,1} = zeros(heigth,width);
end% 求取相位差
% 计算每种频率对应的相位主值
% 输出三种频率的相位主值，用于相差计算
for i = 1:3 % 对于3组中的每一组图片，每一组相同频率的有四张图片I1 = C{i,1};I2 = C{i,2};I3 = C{i,3};I4 = C{i,4};for g = 1:heigthfor k = 1:width          if     I4(g,k)==I2(g,k)&&I1(g,k)>I3(g,k) %四个特殊位置phi{i,1}(g,k)=0;elseif I4(g,k)==I2(g,k)&&I1(g,k)<I3(g,k) %四个特殊位置phi{i,1}(g,k)=pi; elseif I1(g,k)==I3(g,k)&&I4(g,k)>I2(g,k) %四个特殊位置phi{i,1}(g,k)=pi/2;elseif I1(g,k)==I3(g,k)&&I4(g,k)<I2(g,k) %四个特殊位置phi{i,1}(g,k)=3*pi/2;elseif I1(g,k)<I3(g,k) %二三象限phi{i,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) %第一象限phi{i,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) %第四象限phi{i,1}(g,k)=atan((I4(g,k)-I2(g,k))./(I1(g,k)-I3(g,k)))+2*pi;  endend             end
end% 计算相差
% 保存矩阵，用于多频相差的计算
PH1 = phi{1,1};   %频率1
PH2 = phi{2,1};   %频率2
PH3 = phi{3,1};   %频率3% 初始化相差变量
% 多频相差
PH12 = zeros(heigth,width);
PH23 = zeros(heigth,width);
PH123 = zeros(heigth,width);% 计算相差
% 相差计算
% 解相
for g = 1:heigthfor k = 1:width% 计算第一组和第二组的相差if PH1(g,k)>PH2(g,k)PH12(g,k) = PH1(g,k)-PH2(g,k);elsePH12(g,k) = PH1(g,k)+2*pi-PH2(g,k);end% 计算第二组和第三组的相差if PH2(g,k)>PH3(g,k)PH23(g,k) = PH2(g,k)-PH3(g,k);elsePH23(g,k) = PH2(g,k)+2*pi-PH3(g,k);end% plot(1,k);end
end% 计算最终相差
% 相差图案
% 相位解包裹  相位展开
for g = 1:heigthfor k = 1:width if PH12(g,k)>PH23(g,k)PH123(g,k) = PH12(g,k)-PH23(g,k);elsePH123(g,k) = PH12(g,k)+2*pi-PH23(g,k);endend
end% 显示
figure,imshow(mat2gray(PH12));title('1,2外差');   imwrite(mat2gray(PH12),'12外差.bmp');
figure,imshow(mat2gray(PH23));title('2,3外差');   imwrite(mat2gray(PH23),'23外差.bmp');
figure,imshow(mat2gray(PH123));title('1,2,3外差');imwrite(mat2gray(PH123),'123外差.bmp');
figure,imshow(mat2gray(PH1));title('1相位主值');  imwrite(mat2gray(PH1),'1相位主值.bmp');  
figure,imshow(mat2gray(PH2));title('2相位主值');  imwrite(mat2gray(PH2),'2相位主值.bmp');  
figure,imshow(mat2gray(PH3));title('3相位主值');  imwrite(mat2gray(PH3),'3相位主值.bmp');

如果有问题请在评论区留言，或者qq:857467352联系我交流修改。