matlab 解相位_光测力学栅线投影技术-相位求解方法

目前相位求解方法主要有傅里叶方法和相移方法两种，主要有以下区别：

傅里叶方法仅需要一张栅线图像，相移则需要至少三张（两步相移法是一种特殊情况，其相位求解方式严格来说并不是用的相移算法，这个在今后的文章会单独写）
相移法的精度远高于傅里叶方法（这个很好理解，毕竟相移用的信息多）

下面对两种方法进行介绍：

傅里叶方法

傅里叶变换轮廓术是基于离散傅里叶变化的一种相位测量方法，其主要算法流程如图所示。

傅里叶变换轮廓测量方法流程示意图

下面对其理论进行具体介绍。对于被测量物体表面的任意一点处，假设其灰度分布关系有：

$I(x,y)=A(x,y)+B(x,y)\cos [2\pi fx+\varphi (x,y)]\\$

其中

$A$ 为背景灰度，由物体表面的折射率所决定；

$B$ 为光强差，由投影仪栅线的亮度差所决定；

$f$ 为投影栅线的频率。

对于参考平面上的任意一点处，其灰度关系有：

$I(x,y)=A(x,y)+B(x,y)\cos [2\pi fx+{{\varphi }_{\text{0}}}(x,y)]\\$

对于测量物体和参考平面的灰度分布函数，其傅里叶展开式如下所示：

${{I}_{taylor}}(x,y)=r(x,y)\cdot \sum\limits_{N=-\infty }^{\infty }{{{\alpha }_{N}}}\cdot {{e}^{\{i[2\pi fx+N\varphi (x,y)]\}}}\\$

$I_{talor}^{0}(x,y)=\sum\limits_{N=-\infty }^{\infty }{{{\alpha }_{N}}}\cdot {{e}^{\{i[2\pi fx+N{{\varphi }_{\text{0}}}(x,y)]\}}}\\$ 若采用理想的正弦式栅线，其高阶成分主要由物体表面的调制产生，故对于测量来说主要关心其一阶成分。此时可以选取一个窗口（一般为2-D窗口），对其栅线相位变化的方向进行离散傅里叶变化，从而得到其频谱。由于频谱是对称的，故对其去一般分析即可。通过对频谱进行操作提取处其基频分量，如下所示：

${{\hat{I}}_{\text{taylor}}}(x,y)={{\alpha }_{1}}\cdot r(x,y)\cdot {{e}^{\{i[2\pi fx+\varphi (x,y)]\}}}\\$

$\hat{I}_{talor}^{0}(x,y)={{\alpha }_{1}}\cdot {{e}^{\{i[2\pi fx+{{\varphi }_{0}}(x,y)]\}}}\\$

对参考平面的基频分量取其共轭，并将其与被测物体的基频分量相乘有：

$\overline{\hat{I}_{talor}^{0}(x,y)}\centerdot {{\hat{I}}_{\text{taylor}}}(x,y)\text{=}{{\alpha }_{1}}^{2}\cdot r(x,y)\cdot {{e}^{i\cdot \Delta \varphi (x,y)}}\\$

对其取对数即可得到其相位差距

$\Delta \varphi (x,y)$ ：

$\log \left[ \overline{\hat{I}_{talor}^{0}(x,y)}\centerdot {{{\hat{I}}}_{\text{taylor}}}(x,y) \right]\text{=log}\left[ {{\alpha }_{1}}^{2}\cdot r(x,y) \right]+i\Delta \varphi (x,y)\\$

对于真实的测量过程中，多数情况下并不会人为引入真实的参考平面，故直接通过被测物体的基频分量，直接提取出相位信息，再通过假设的参考平面，也可得到物体表面的相位变化信息。

傅里叶轮廓相移术的测量精度较低，但是仅仅利用一张图像即可恢复出物体表面的三维形貌。

相移方法

重写物体表面的灰度关系，如下：

$I(x,y)=A(x,y)+B(x,y)\cos [\theta +\varphi (x,y)]\\$

其中

$\theta$ 为投影相位。对于N步相移方法来说，每一步相移与上一步相移的投影相位相差

$\frac{2\pi }{N}$ 的相位，即以下方程组：

${{I}_{i}}(x,y)=A(x,y)+B(x,y)\cos \left[ \frac{2\pi }{N}\cdot (i-1)+\varphi (x,y) \right]\\$

其中

$i=0,1,2,\ldots ,N-2,N-1$ 。方程组共有N个方程，三个未知量，原理上当

$N\ge 3$ 时即可求解。在实际操作中，为了保证求解的精度，一般取

$N\ge 4$ 使得求解体系变为超定方程求解。以六步相移法为例，采集的图像如下图所示。

六步相移图像（仔细看这几个栅线是有区别的！）

当前方程组为非线性方程，看似无法用直接方法求解。但是注意到方程可以转化为以下形式：

$\begin{align} & {{I}_{i}}(x,y)=A(x,y)+B(x,y)\cdot \cos [\varphi (x,y)]\cdot \cos [\frac{2\pi }{N}\cdot (i-1)] \\ & -B(x,y)\cdot \sin [\varphi (x,y)]\cdot \sin [\frac{2\pi }{N}\cdot (i-1)] \\ & =A(x,y)+{{B}_{1}}(x,y)\cdot \cos [\frac{2\pi }{N}\cdot (i-1)]+{{B}_{2}}(x,y)\cdot \sin [\frac{2\pi }{N}\cdot (i-1)] \end{align}\\$

其中

$i=0,1,2,\ldots ,N-2,N-1$ 。通过简单的三角函数展开，将待求量为

$A(x,y)$ 、

$B(x,y)$ 和

$\varphi (x,y)$ 的非线性方程组，转化为关于待求量

$A(x,y)$ 、

${{B}_{1}}(x,y)$ 和

${{B}_{2}}(x,y)$ 的线性方程组。其中

${{B}_{1}}(x,y)=B(x,y)\cdot \cos [\varphi (x,y)]$ ，

${{B}_{2}}(x,y)=-B(x,y)\cdot \sin [\varphi (x,y)]$ 。即待求方程组可化为以下矩阵模式：

$\left( \begin{matrix} 1 & \cos [0\cdot \frac{2\pi }{N}] & \sin [0\cdot \frac{2\pi }{N}] \\ 1 & \cos [1\cdot \frac{2\pi }{N}] & \sin [1\cdot \frac{2\pi }{N}] \\ 1 & \cos [2\cdot \frac{2\pi }{N}] & \sin [2\cdot \frac{2\pi }{N}] \\ \vdots & \vdots & \vdots \\ 1 & \cos [(N-2)\cdot \frac{2\pi }{N}] & \sin [(N-2)\cdot \frac{2\pi }{N}] \\ 1 & \cos [(N-1)\cdot \frac{2\pi }{N}] & \sin [(N-1)\cdot \frac{2\pi }{N}] \\ \end{matrix} \right)\left( \begin{matrix} A(x,y) \\ {{B}_{1}}(x,y) \\ {{B}_{2}}(x,y) \\ \end{matrix} \right)=\left( \begin{matrix} {{I}_{\text{0}}} \\ {{I}_{1}} \\ {{I}_{2}} \\ \vdots \\ {{I}_{N-2}} \\ {{I}_{N-1}} \\ \end{matrix} \right)\\$

可以简写为：

${{\mathbf{N}}_{(N\times 3)}}{{\mathbf{x}}_{(3\times 1)}}\text{=}{{\mathbf{I}}_{(N\times 1)}}\\$

由最小二乘法求解超定方程的理论，方程的最优解为：

$\mathbf{x=(}{{\mathbf{N}}^{\mathbf{T}}}\mathbf{N}{{\mathbf{)}}^{\mathbf{-1}}}{{\mathbf{N}}^{\mathbf{T}}}\mathbf{I}\\$

由于根据对

${{B}_{\text{1}}}(x,y)$ 与

${{B}_{2}}(x,y)$ 的构造可知：

$\tan [\varphi (x,y)]=-\frac{{{B}_{2}}(x,y)}{{{B}_{1}}(x,y)}\\$

结合超定方程的最小二乘解法，求得

$\varphi (x,y)$ 的最优解为：

$\varphi (x,y)=\arctan \left[ \frac{\sum\limits_{i=0}^{N}{I{}_{i}\cdot \sin (i\cdot \frac{2\pi }{N})}}{\sum\limits_{i=0}^{N}{I{}_{i}\cdot \cos (i\cdot \frac{2\pi }{N})}} \right]\\$

由上图所示的六步相移图像，其所对应的相位图像（已去除背景噪点）如下图所示：

包裹相位

由于投影仪的栅线图案为周期性图像，故相位图中的相位值也具有一定的周期范围

$[-\pi ,\pi ]$ 。为了获取真实的形貌，我们要将离散的周期性相位图转化为连续相位图，这一步一般被称为解包裹。通过每当遇到一个新的周期，对下一个周期补偿若干个

$2\pi$ 或

$-2\pi$ 周期相位，使得两个周期从相位跳跃转化为连续。即如下公式所示：

${\varphi }'(x,y)=\varphi (x,y)+2k\pi \\$

目前准确度较高的解包裹方式为基于质量的解包裹(以后也将介绍各种解包裹算法)，其通过选取图像梯度中较好的相位点向四周进行解包裹，有效的避免了按行解包裹和按列解包裹中的整体相位错位。解包裹的效果如下图所示。

解包裹相位

注意到我们在人脸模型后架设了一块标准平面，已它作为参考平面，利用原始相位减去参考平面相位，即可得到具有人脸的相位高度模型，如下图所示。

相位高度分布

进一步，对系统的各参数进行标定，即可得到人脸模型的真实尺度。

相移轮廓测量技术相对于傅里叶轮廓测量技术，测量精度较高，但是至少需要三张图像才能够恢复出物体表面的三维形貌。属于典型的牺牲测量时间来换取精度的操作。