三频外差相位展开

之前的文章介绍了PMP的基本原理，我们了解到待测物体的高度信息包含在条纹的相位中。下面介绍下如何从条纹中提取连续的相位。

$$\psi (x,y)=-atan(\frac{\sum I_i(x,y)\ast sin({\delta }_i)}{\sum I_i(x,y)\ast cos({\delta }_i)})\text{\hspace{1em}(1)}$$

从相机拍摄回来的条纹图案灰度值直接提取到相位如公式(1)所示，具体推导过程请参考我之前的文章。因为atan函数的原因，条纹相位被截断在$[-\pi ,\pi )$。如下图所示，我们的目标是把下图的三角锯齿波展开成红色虚线的连续相位。

把截断相位展开为连续相位的过程叫相位展开，一般文献称公式(1)的相位为包裹相位(wrapped phase)，我觉得这个名字太难听了，我一般称之为截断相位。相位展开的方法由很多，大体上可分为空间相位展开方法和时间相位展开方法，这里不对他们的做过多的介绍，感兴趣的可以查阅相关文献。本文介绍的三频外差相位展开法就属于时间相位展开方法的一种。

外差原理

假设有两组不同周期（频率）的条纹$T_1$，$T_2$，其中，$T_1$ < $T_2$，则两组不同频率条纹差分的等效周期可表示为：$T_{12}=\frac{T_1\ast T_2}{T_2-T_1}$。使用相移法得到两组截断相位${\psi }_1$，${\psi }_2$，把两组截断相位差值得到等效截断相位${\psi }_{12}$，如果两个条纹的等效周期能完全覆盖整个测量范围，即$T_{12}\ge width$，则等效相位将不是截断相位，而是连续相位，即不需要进行相位展开，然后根据等效相位即可求得截断相位的级次，从而对截断相位进行展开。

三频外差

由外差原理可知，只要选择合适的条纹周期，就可以方便的展开截断相位。但是实际上双频外差的条纹周期比较大，即在整个测量范围内具有较少的条纹。由前面的分析可知，条纹周期大，测量系统的分辨率降低，并且最大测量深度也会降低。为了提高系统测量范围以及精度，人们拓展双频外差到三频外差。具体相位展开过程如下，$\psi$表示截断相位，$\varphi$表示展开后的连续相位：

1. 选择合适的三种条纹周期$T_1$，$T_2$，$T_3$，满足，$T_1<T_2<T_3<T_{12}<T_{23}<T_{123}$，且$T_{123}\ge width$，有：

$$\{\begin{array}{l}{T}_{12}=\frac{T_1{T}_2}{T_2-T_1}\\ T_{23}=\frac{T_2{T}_3}{T_3-T_2}\\ T_{123}=\frac{T_{12}{T}_{23}}{T_{23}-T_{12}}\end{array}\text{\hspace{1em}(2)}$$

2. 按照公式(1)分别获得三种频率条纹的截断相位，并缩放到$[0,2\pi )$区间；
3. 按公式(2)差分截断相位，并缩放到$[0,2\pi )$：

$$\{\begin{array}{l}{\psi }_{12}={\psi }_1-{\psi }_2\\ {\psi }_{23}={\psi }_2-{\psi }_3\\ {\psi }_{123}={\psi }_{12}-{\psi }_{23}\end{array}\text{\hspace{1em}(3)}$$

4. 因为$T_{123}\ge width$，所以公式(2)中的截断相位${\psi }_{123}$实际上是连续相位，记为${\varphi }_{123}$，根据${\varphi }_{123}$展开截断相位${\psi }_{23}$：

$$k_{23}=floor(\frac{{\varphi }_{123}\ast T_{123}/T_{23}-{\psi }_{23}}{2\pi }+0.5)\text{\hspace{1em}(4)}$$

$${\varphi }_{23}={\psi }_{23}+2\pi k\text{\hspace{1em}(5)}$$

5. 使用公式(4)，(5)递归求解${\varphi }_{12},{\varphi }_3,{\varphi }_2,{\varphi }_1$；

$$\{\begin{array}{l}{k}_{12}=floor(\frac{{\varphi }_{23}\ast T_{23}/T_{12}-{\psi }_{12}}{2pi}+0.5)\\ {\varphi }_{12}={\psi }_{12}+2\pi k_{12}\\ k_3=floor(\frac{{\varphi }_{12}\ast T_{12}/T_3-{\psi }_3}{2\pi }+0.5)\\ {\varphi }_3={\psi }_3+2\pi k_3\\ k_2=floor(\frac{{\varphi }_3\ast T_3/T2-{\psi }_2}{2\pi }+0.5)\\ {\varphi }_2={\psi }_2+2\pi k_2\\ k_1=floor(\frac{{\varphi }_2\ast T_2/T_1-{\psi }_1}{2\pi }+0.5)\\ {\varphi }_1={\psi }_1+2\pi k_1\end{array}$$

6. 选择频率较高的条纹绝对相位${\varphi }_1$作为系统的绝对连续相位；

仿真代码

% 频率，注意，不是条纹周期
f1 = 70;
f2 = 64;
f3 = 59;

T1 = 1/f1;
T2 = 1/f2;
T3 = 1/f3;

T12 = T1 * T2 / (T2 - T1);
T23 = T2 * T3 / (T3 - T2);
T123 = T12 * T23 / (T23 - T12);

f12 = 1 / T12;
f23 = 1 / T23;
f123 = 1 / T123;

width = 1024;
height = 1140;
N = 4;
I1 = zeros(height, width, N);
I2 = zeros(height, width, N);
I3 = zeros(height, width, N);

[x, ~] = meshgrid(linspace(0, width, width), linspace(0, height, height));
A = 135;
B = 115;
for i = 1 : N
    I1(:,:,i) = A + B * cos(2 * pi * f1 * x + 2 * pi * (i - 1) / N) - 5 + 5 * rand(height, width);
    I2(:,:,i) = A + B * cos(2 * pi * f2 * x + 2 * pi * (i - 1) / N) - 5 + 5 * rand(height, width);
    I3(:,:,i) = A + B * cos(2 * pi * f3 * x + 2 * pi * (i - 1) / N) - 5 + 5 * rand(height, width);
end

h = fspecial('gaussian', 5, 5);
% 滤波
for i= 1 : N
    I1(:,:,i) = imfilter(I1(:,:,i), h);
    I2(:,:,i) = imfilter(I2(:,:,i), h);
    I3(:,:,i) = imfilter(I3(:,:,i), h);
end 

% 截断相位
phi1 = atan2(I1(:,:,4) - I1(:,:,2), I1(:,:,1) - I1(:,:,3));
phi2 = atan2(I2(:,:,4) - I2(:,:,2), I2(:,:,1) - I2(:,:,3));
phi3 = atan2(I3(:,:,4) - I3(:,:,2), I3(:,:,1) - I3(:,:,3));

% 截断相位映射到0~2pi
phi1(phi1 < 0) = phi1(phi1 < 0) + 2 * pi;
phi2(phi2 < 0) = phi2(phi2 < 0) + 2 * pi;
phi3(phi3 < 0) = phi3(phi3 < 0) + 2 * pi;

%%
% 外差
phi12 = phi1 - phi2;
phi12(phi12 < 0) = phi12(phi12 < 0) + 2 * pi;
phi23 = phi2 - phi3;
phi23(phi23 < 0) = phi23(phi23 < 0) + 2 * pi;
phi123 = phi12 - phi23;
phi123(phi123 < 0) = phi123(phi123 < 0) + 2 * pi;

%%
% 解相位
% 1. 根据phi123求解phi23的绝对相位
k23 = floor((phi123 * T123 / T23 - phi23)/(2 * pi) + 0.5);
uph23 = phi23 + 2 * pi * k23;

% 2. 根据phi23的绝对相位解phi12绝对相位x
k12 = floor((uph23 * T23 / T12 - phi12)/(2 * pi) + 0.5);
uph12 = phi12 + 2 * pi * k12;

% 3. 根据phi12的绝对相位解phi3绝对相位
k3 = floor((uph12 * T12 / T3 - phi3)/(2 * pi) + 0.5);
uph3 = phi3 + 2 * pi * k3;

% 4. 根据phi3绝对相位解phi2绝对相位
k2 = floor((uph3 * T3 / T2 - phi2)/(2 * pi) + 0.5);
uph2 = phi2 + 2 * pi * k2;

% 根据phi2绝对相位解phi1绝对相位
k1 = floor((uph2 * T2 / T1 - phi1)/(2 * pi) + 0.5);
uph1 = phi1 + 2 * pi * k1;