结构光系统标定（三）基于多项式拟合的结构光系统标定

结构光系统标定的方法实际上还挺多的，我只讨论分析我读过且实现过的两种，两种方法都比较经典、应用较多、速度较快且操作简便的。这篇文章会讨论第一种，基于多项式拟合的结构光系统标定。当然这名字是我自己给它取的，如果不合适的话，也请见谅。

基于多项式拟合的结构光系统标定

按我的理解来说，基于多项式拟合的结构光系统标定，是将本需要直接测量的结构光系统中的长度值，隐式地表达在一个多项式中，同时引入经典的摄像机小孔成像模型，结合张正友摄像机标定法，从而获得待标定的相位——深度映射中的相关参数，以恢复被测物体的三维形貌。

实际上，研究者对基于多项式拟合的结构光系统标定提出了许多种不同的系统模型和标定方法，这里要讲的结构光系统标定方法是由达飞鹏和盖绍彦提出的，参考文献[1]是原文，[2]后来他们出的书，其实是一样的。这个模型和标定方法具有较高的速度和可操作性。我在这尽可能详细地按照我的理解来分析他们的模型，以及实际操作计算中，并不需要做的文章里面一样复杂，可以用些工具简化中间步骤的计算，所以有些计算会和文章不一样，但道理是一样的。

这个系统模型要引入摄像机小孔成像模型，以及张正友标定，在这文章里我就默认各位读者都熟悉成像模型和相机标定了，这里直接用，而不会去介绍为什么，有需要就自己网上搜吧，这方面的解读多得数都数不过来。

系统模型

该方法主要建立摄像机坐标系和被测物体上某点的相位、相位和世界坐标系之间的关系。该系统模型区别于标准结构光系统模型，对摄像机和投影仪的摆放位置没有严格的要求。

${{\Omega }_{\text{W}}}\left( {{O}_{w}}{{X}_{w}}{{Y}_{w}}{{Z}_{w}} \right)$ 、 ${{\Omega }_{\text{C}}}\left( {{O}_{c}}{{X}_{c}}{{Y}_{c}}{{Z}_{c}} \right)$ 分别表示世界坐标系和摄像机坐标系， $o{{u}_{c}}{{v}_{c}}$ 组成摄像机成像面坐标系， $O_{w}^{'}$ 为投影中心， $A$ 为被测物体上一点。其中世界坐标系 ${{\Omega }_{\text{W}}}\left( {{O}_{w}}{{X}_{w}}{{Y}_{w}}{{Z}_{w}} \right)$ 是根据投影仪的位置建立的，参考平面 ${{O}_{w}}{{X}_{w}}{{Y}_{w}}$ 平行于投影平面， ${{Y}_{w}}$ 轴平行于正弦条纹的条纹方向， ${{Z}_{w}}$ 轴过投影中心。在该模型中，通过旋转矩阵 $R=\left[ \begin{matrix} {{r}_{11}} & {{r}_{12}} & {{r}_{13}} \\ {{r}_{21}} & {{r}_{22}} & {{r}_{23}} \\ {{r}_{31}} & {{r}_{32}} & {{r}_{33}} \\ \end{matrix} \right]$ 和平移矩阵 $t=\left[ \begin{matrix} {{t}_{x}} \\ {{t}_{y}} \\ {{t}_{z}} \\ \end{matrix} \right]$ 描述世界坐标系 ${{\Omega }_{\text{W}}}$ 和摄像机坐标系 ${{\Omega }_{\text{C}}}$ 之间的转换关系，通过小孔成像模型描述成像平面和摄像机坐标系之间的关系，如：点 $A$ 在世界坐标系中坐标为 $\left( {{X}_{w}},{{Y}_{w}},{{Z}_{w}} \right)$ ，在摄像机坐标系中坐标为 $\left( {{X}_{c}},{{Y}_{c}},{{Z}_{c}} \right)$ ，在成像平面上的坐标为 $\left( {{u}_{c}},{{v}_{c}} \right)$ ，那么

$\left[ \begin{matrix} {{X}_{c}} \\ {{Y}_{c}} \\ {{Z}_{c}} \\ \end{matrix} \right]=\left[ \begin{matrix} R & t \\ \end{matrix} \right]\left[ \begin{matrix} {{X}_{w}} \\ {{Y}_{w}} \\ {{Z}_{w}} \\ 1 \\ \end{matrix} \right]$

$s\left[ \begin{matrix} {{u}_{c}} \\ {{v}_{c}} \\ 1 \\ \end{matrix} \right]=\left[ \begin{matrix} {{f}_{cx}} & 0 & {{c}_{cx}} \\ 0 & {{f}_{cy}} & {{c}_{cy}} \\ 0 & 0 & 1 \\ \end{matrix} \right]\left[ \begin{matrix} {{X}_{c}} \\ {{Y}_{c}} \\ {{Z}_{c}} \\ \end{matrix} \right]$

以上是几个坐标系之间的转换关系，接下来是通过几何关系，找出相位和世界坐标系之间的关系。

$A'$ 为点 $A$ 在参考平面 ${{O}_{w}}{{X}_{w}}{{Y}_{w}}$ 上的投影， ${{O}_{w}}'A$ 与 ${{O}_{w}}{{X}_{w}}{{Y}_{w}}$ 平面交于点 $B$ ， $A''$ 、 $B'$ 分别为 $A'$ 、点 $B$ 在 ${{X}_{w}}$ 轴上的投影。易知， $AA'//{{O}_{w}}'{{O}_{w}}$ ，则 $\vartriangle BAA'\sim \vartriangle B{{O}_{w}}'{{O}_{w}}$ ，有

$\frac{{{O}_{w}}A'}{{{O}_{w}}B}=\frac{{{O}_{w}}{{O}_{w}}'-A'A}{{{O}_{w}}{{O}_{w}}'}$

同时， $A'A''//BB'$ ，则 $\vartriangle OA'A''\sim \vartriangle OBB'$ ，有

$\frac{{{O}_{w}}A'}{{{O}_{w}}B}=\frac{{{O}_{w}}A''}{{{O}_{w}}B'}$ ，

联立以上两式，可得

$\frac{{{O}_{w}}A''}{{{O}_{w}}B'}=\frac{{{O}_{w}}{{O}_{w}}'-A'A}{{{O}_{w}}{{O}_{w}}'}$ ，

上式中， ${{O}_{w}}{{O}_{w}}'=l$ ， $AA'$ 和 ${{O}_{w}}A''$ 分别为点 $A$ 在世界坐标系中坐标 ${{Z}_{w}}$ 和 ${{X}_{w}}$ ；由于 $B'$ 为点 $B$ 在 ${{X}_{w}}$ 轴上的投影，且 ${{Y}_{w}}$ 轴平行于正弦条纹的条纹方向，所以 $B'$ 和 $B$ 的相位相等，也等于点 $A$ 的相位（摄像机拍到物体A点的相位，实际等于参考平面上B点的相位），因此

$OB'=\frac{\lambda }{2\pi }\left( \Phi -{{\Phi }_{0}} \right)$

上式中， $\lambda$ 为光栅的周期， ${{\Phi }_{0}}$ 为原点 ${O}_{w}$ 处的相位。将式（2-53）代入式（2-52）中，可得到相位和点在世界坐标系的坐标 $\left( {{X}_{w}},{{Y}_{w}},{{Z}_{w}} \right)$ 之间的关系，

$\Phi =\frac{\frac{2\pi l}{\lambda }{{X}_{w}}-{{\Phi }_{0}}{{Z}_{w}}+l{{\Phi }_{0}}}{l-{{Z}_{w}}}$

结合上面世界坐标系和摄像机坐标系的转换关系，可得相位 $\Phi$ 和点 $A$ 在摄像机坐标系的坐标 $\left( {{X}_{c}},{{Y}_{c}},{{Z}_{c}} \right)$ 间的关系

$\Phi =\frac{{{a}_{1}}{{X}_{c}}+{{a}_{2}}{{Y}_{c}}+{{a}_{3}}{{Z}_{c}}+{{a}_{4}}}{{{a}_{5}}{{X}_{c}}+{{a}_{6}}{{Y}_{c}}+{{a}_{7}}{{Z}_{c}}+{{a}_{8}}}$

式中， $\left\{ \begin{array}{*{35}{l}} \begin{array}{*{35}{l}} {{a}_{1}}=\frac{2\pi \left( {{r}_{11}}-{{\Phi }_{0}}{{r}_{31}} \right)l}{\lambda } \\ {{a}_{2}}=\frac{2\pi \left( {{r}_{12}}-{{\Phi }_{0}}{{r}_{32}} \right)l}{\lambda } \\ {{a}_{3}}=\frac{2\pi \left( {{r}_{13}}-{{\Phi }_{0}}{{r}_{33}} \right)l}{\lambda } \\ {{a}_{4}}=\frac{2\pi \left( {{t}_{x}}-{{\Phi }_{0}}{{t}_{z}} \right)l}{\lambda }+{{\Phi }_{0}}l \\ \end{array} \\ \begin{array}{*{35}{l}} {{a}_{5}}=-\frac{2\pi {{r}_{31}}l}{\lambda } \\ {{a}_{6}}=-\frac{2\pi {{r}_{32}}l}{\lambda } \\ {{a}_{7}}=-\frac{2\pi {{r}_{33}}l}{\lambda } \\ {{a}_{8}}=-\frac{2\pi {{t}_{z}}l}{\lambda }+l \\ \end{array} \\ \end{array} \right.$

得到以上关系后，可以看出，系统标定过程就是在求 ${{a}_{1}}-{{a}_{8}}$ 八个参数的过程。

结构光系统标定步骤

从上面推导的式子里可以看出来，要求 ${{a}_{1}}-{{a}_{8}}$ 八个参数，也就是解一组八元一次方程，需要得到至少8组相位和在摄像机坐标系的坐标的数据才能解出这个方程。可见有两类数据要获取，一个是相位，一个是三维坐标，我是通过下面的方式获得的：

1.获取相位，我会用PMP方法获取相位，用Multi-frequency相位解包裹获得连续相位图；

2.一系列在摄像机坐标系的坐标，我会通过检测图像中的棋盘角点，经过张正友相机标定，利用标定的外参，计算这些角点在摄像机坐标系中的三维坐标。

以上用到的方法，在这里都默认大家已经熟悉了，这里不会再去讲，如果不清楚，PMP和相位解包裹可以翻我以往的文章，或者自己去看文献，张氏相机标定可以看网上的讲解。

下面我就列出标定的各个步骤：

步骤一，将棋盘摆放在被测范围内的任意位置，用摄像机拍摄棋盘，检测棋盘角点（OpenCV或者MATLAB都有对应的函数），记录各角点在成像平面中坐标 $\left( {{u}_{ci}},{{v}_{ci}} \right)\left( i=1,2,3\cdots ,n \right)$ 。

步骤二，保持棋盘不动，利用投影仪，将正弦条纹图样投影在棋盘上，获取此时的相位图，记录各角点位置的相位 ${{\Phi }_{i}}\left( i=1,2,3\cdots ,n \right)$ 。文章中用的是标定板，因为标定板还挺贵的，所以我就直接用棋盘了。

步骤三，多次调整棋盘的摆放位置，重复步骤一和步骤二，获得多组角点坐标 $\left( {{u}_{ci}},{{v}_{ci}} \right)$ 和对应的相位 ${{\Phi }_{i}}$ 。为更精确标定摄像机，我建议获取9张或以上棋盘图像，否则后果自负哦。

步骤四，通过张正友摄像机标定方法标定摄像机，获得摄像机内参 ${{A}_{c}}=\left[ \begin{matrix} {{f}_{cx}} & 0 & {{c}_{cx}} \\ 0 & {{f}_{cy}} & {{c}_{cy}} \\ 0 & 0 & 1 \\ \end{matrix} \right]$ 及各幅棋盘图像对应的外参 $\left[ \begin{matrix} {{R}_{i}} & {{t}_{i}} \\ \end{matrix} \right]$ ，根据外参可以计算出各幅棋盘图像中的角点 $\left( {{u}_{ci}},{{v}_{ci}} \right)$ 对应在摄像机坐标系下的坐标 $\left( {{X}_{ci}},{{Y}_{ci}},{{Z}_{ci}} \right)$ 。我的理解和这个标定方法提出的文章主要区别就在 $\left( {{X}_{ci}},{{Y}_{ci}},{{Z}_{ci}} \right)$ 的计算不一样，用张世标定的数据可以很快求出所需的坐标，但文章里面的计算相对复杂，但道理其实类似的，感兴趣可以读一读文章。

步骤五，利用步骤三得到的棋盘角点位置的相位 ${{\Phi }_{i}}$ 和步骤四得到的棋盘角点在摄像机坐标系下的坐标 $\left( {{X}_{ci}},{{Y}_{ci}},{{Z}_{ci}} \right)$ ，就构造多项式

${{a}_{1}}{{X}_{c}}+{{a}_{2}}{{Y}_{c}}+{{a}_{3}}{{Z}_{c}}+{{a}_{4}}-{{a}_{5}}\Phi {{X}_{c}}-{{a}_{6}}\Phi {{Y}_{c}}-{{a}_{7}}\Phi {{Z}_{c}}-{{a}_{8}}\Phi =0$

把所有数据写成矩阵形式

$\left[ \begin{array}{*{35}{l}} {{X}_{c1}} & {{Y}_{c1}} & {{Z}_{c1}} & 1 & -{{\Phi }_{1}}{{X}_{c1}} & -{{\Phi }_{1}}{{Y}_{c1}} & -{{\Phi }_{1}}{{Z}_{c1}} & -{{\Phi }_{1}} \\ {{X}_{c2}} & {{Y}_{c2}} & {{Z}_{c2}} & 1 & -{{\Phi }_{2}}{{X}_{c2}} & -{{\Phi }_{2}}{{Y}_{c2}} & -{{\Phi }_{2}}{{Z}_{c2}} & -{{\Phi }_{2}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ {{X}_{cn}} & {{Y}_{cn}} & {{Z}_{cn}} & 1 & -{{\Phi }_{n}}{{X}_{cn}} & -{{\Phi }_{n}}{{Y}_{cn}} & -{{\Phi }_{n}}{{Z}_{cn}} & -{{\Phi }_{n}} \\ \end{array} \right]\left[ \begin{array}{*{35}{l}} {{a}_{1}} \\ {{a}_{2}} \\ {{a}_{3}} \\ {{a}_{4}} \\ {{a}_{5}} \\ {{a}_{6}} \\ {{a}_{7}} \\ {{a}_{8}} \\ \end{array} \right]=\left[ \begin{array}{*{35}{l}} 0 \\ \vdots \\ 0 \\ \end{array} \right]$

就能通过最小二乘法求出 ${{a}_{1}}-{{a}_{8}}$ 八个参数。

从相位到三维坐标

完成系统标定后即可得到 ${{a}_{1}}-{{a}_{8}}$ 八个参数，当实际测量时，摄像机获取的图像上任意一点 $\left( {{u}_{c}},{{v}_{c}} \right)$ 以及该点对应的相位 $\Phi$ ，怎么才能求出该点在摄像机坐标系下的三维坐标 $\left( {{X}_{c}},{{Y}_{c}},{{Z}_{c}} \right)$ ，说实话，这个推导过程困扰了我时间最长。

前面求参数的过程和程序我花的时间不是特别多，就解决了。然而可能是我的数学能力太差了，推导最后一步的公式一直卡着，放弃了很多次，有了新点子推了很多次，最后受一篇古老文章启发才知道，原来是摄像机小孔成像模型的地方卡住我了，换个形式就好办多了。

这个推导，文章一般都不会详细写，我还问了一些做出来的博主，我得到的回答都是这个很好推的，你再想清楚成像模型就好了。当时是真崩溃······

以下我就直接列出我的推导出来吧，希望能帮到和我一样的数学渣渣、结构光小白，在迷茫中有人解惑真的会感动哭。

已知任意一点 $\left( {{u}_{c}},{{v}_{c}} \right)$ 以及该点对应的相位 $\Phi$ ，还有摄像机的内参，可以慢慢构造下面的矩阵。

所以最后的三维坐标 $\left( {{X}_{c}},{{Y}_{c}},{{Z}_{c}} \right)$ 和任意一点 $\left( {{u}_{c}},{{v}_{c}} \right)$ 以及该点对应的相位 $\Phi$ 的关系就是下面这个式子了。

$\left[ \begin{matrix} {{X}_{c}} \\ {{Y}_{c}} \\ {{Z}_{c}} \\ \end{matrix} \right]={{\left[ \begin{matrix} {{f}_{cx}} & 0 & {{c}_{cx}}-{{u}_{c}} \\ 0 & {{f}_{cy}} & {{c}_{cy}}-{{v}_{c}} \\ {{a}_{1}}-{{a}_{5}}\Phi & {{a}_{2}}-{{a}_{6}}\Phi & {{a}_{3}}-{{a}_{7}}\Phi \\ \end{matrix} \right]}^{\text{-}1}}\left[ \begin{matrix} 0 \\ 0 \\ {{a}_{8}}\Phi -{{a}_{4}} \\ \end{matrix} \right]$

实验中间步骤及结果

最后稍微放一点标定过程中的一些中间步骤图片，和最后的标定结果。系统不同，大家的数据也不同，所以这部分意义不大，读者们也可以略过。

PMP所用的条纹：我用MATLAB生成4个频率、纵向的三步相移条纹图样，尺寸为720×500像素，频率分别是1/400、1/80、1/20、1/6

棋盘的角点检测结果（用的是OpenCV的角点检测函数）

标定的相机内参 ${{A}_{c}}=\left[ \begin{matrix} 3319.37 & 0 & 453.755 \\ 0 & 3303.42 & 198.200 \\ 0 & 0 & 1 \\ \end{matrix} \right]$ ，标定后的八个参数 $\left[ \begin{matrix} \begin{matrix} {{a}_{1}} \\ {{a}_{2}} \\ {{a}_{3}} \\ {{a}_{4}} \\ \end{matrix} \\ \begin{matrix} {{a}_{5}} \\ {{a}_{6}} \\ {{a}_{7}} \\ {{a}_{8}} \\ \end{matrix} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \text{6}\text{.35764e-03} \\ \text{1}\text{.81497e-04} \\ \text{1}\text{.21451e-04} \\ \text{9}\text{.99980e-01} \\ \end{matrix} \\ \begin{matrix} \text{9}\text{.57739e-07} \\ \text{1}\text{.44972e-07} \\ \text{3}\text{.89189e-06} \\ \text{-3}\text{.26908e-05} \\ \end{matrix} \\ \end{matrix} \right]$ 。

重建结果可以用直接用深度图来显示，不过应该要阈值化后才能正常显示，图（a）。

我还用PCL库做了3D点云可视化，图（b）。

参考文献：

[1] 盖绍彦, 达飞鹏. 一种新的相位法三维轮廓测量系统模型及其标定方法研究[J]. 自动化学报, 2007, 33(9): 902-910.

[2] 达飞鹏, 盖绍彦. 光栅投影三维精密测量[M]. 北京: 科学出版社, 2011.

[3] 刘顺涛, 骆华芬, 陈雪梅, 等. 结构光测量系统的标定方法综述[J]. 激光技术, 2015, 39(2): 252-258.

声明：本文中的图片、公式等都和作者本人的毕业论文有关，如果抄袭盗用吃亏是你自己哦~实在有需要，这里建议你自己码一下公式，自己做一下相关的图哦亲。作者水平有限，如文中有错，请务必留言指正。如有学习交流需要，也可通过邮箱zhenyuchung@m.scnu.edu.cn联系我，大家一起讨论学习。