基于移动最小二乘曲面的点云对齐（一） 隐式平面的生成

本文重点

　　这次主要介绍一种点云对齐的方法,多视数据最近迭代(ICP)对齐是最常用的点云对齐方法,为了提高对齐的精度及稳定性我们使用一种基于移动最小二乘(MLS)曲面的ICP多视数据对齐方法.该方法无需对数据进行额外的去噪和数据分割.对于优化噪声点的点云对齐可以采用本方法进行点云对齐。

 

1 多视数据相关概念

　　因为扫面仪(传感器)的自身扫描范围限制,工业扫描往往不能一次对完整的模具进行扫描,需要通过多次转换视角进行模具的扫描.由于各视角的点云数据定义于各自的局部坐标系下,我们就需要对这些多视角测量数据进行坐标系转换使其在统一的坐标系下.

---

　　什么是点云对齐？

---

　　点云对齐就是通过坐标平移或旋转将具有相似特征的一组或几组点云转换到统一的全局坐标系下，多视角点云对齐通常有三类方法:

• 利用度量装置自动记录视角间变换,然后通过一组数据作为标准数据,其余数据逆向转换回全局坐标系;
• 被测模具表面贴足够的标识点，通过计算标志点之间的拓扑关系进行点云间的对齐;
• 测量时使不同视角数据拥有足够的数据重叠，通过重叠部分的坐标变换关系，将局部坐标系下的点云转换到全局坐标系下。

 

2 点云对齐的方法

　　直接通过设备记录点云运动的方法成本太高，需要一台记录装置，本文主要针对目前使用最广的先进行手动粗对齐，在用IPC算法精对齐。

　　移动最小二乘（Moving Least-Squares,MLS）曲面逼近方法，不仅能适应任意拓扑形状的点云估计，而且逼近的MLS曲面具有良好的几何性质，可以通过确定逼近误差阈值排除噪点。

---

　　什么是MLS曲面？

---

　　国外研究者Levin将MLS曲面定义为满足对投影映射具有不变性的空间点集合，即

S(x) = {x∈R³|f(x) = x}

　　Amenta在次基础上对MLS曲面投影进行了更精确的定义。即沿法矢量场n(x)满足能量函数e(y,a)的极小值点集合（y为空间点，a为投影方向），并给出了由空间一点x₀到MLS曲面的迭代投影计算方法。

　　如图1所示，从初始点x_o开始沿加权法矢量n(x_o)优化能量函数e(y,a)得到点x₁，其中能量函数e(y,a)的定义为

${\rm{e}}\left( {{\rm{y}},{\rm{a}}} \right) = \mathop \sum \limits_{i = 1}^k {a^T} \cdot {\left( {y - {p_i}} \right)^2}\theta \left( {\left\| {y - {p_i}} \right\|} \right)$   （公式1）

        公式1中：k为输入点y（如图1中的x_i）的邻域点集{pi}的个数（一般k取15）；a取输入点y处的加权法矢量n(x_i),其中n(x_i)可通过邻近点的加权平均得到，即

${\rm{n}}\left( {{{\rm{x}}_i}} \right) = \;\frac{{\mathop \sum \nolimits_{j = 1}^k {n_{{p_j}}}{v_j}}}{{\left\| {\mathop \sum \nolimits_{j = 1}^k {n_{{p_j}}}{v_j}} \right\|}}$   （公式2）

　　公式2中v_j = ɵ(||x_i-p_i||)以及公式1中的v_j = ɵ(||y-p_i||)均为高斯函数，表示为

$\theta \left( {{\rm{xi}} - {\rm{pi}}} \right) = {e^{ - \frac{{{{\left\| {{\rm{xi}} - {\rm{pi}}} \right\|}^2}}}{{{h^2}}}}}$    （公式3）

　　公式3中h为高斯函数的影响因子，取值较小时，MLS曲面比较靠近离散点云，但曲面不光滑；取值较大时，MLS曲面光滑，但在高曲率区域偏差增大，为获取光滑且偏差合适的MLS曲面，h可取点云数据中所有邻近点对间距离的平均值。

图1 计算MLS曲面图解

　　进一步对公式2求关于变量y的偏导，可以得到MLS曲面的隐式定义：

${\rm{g}}\left( {\rm{x}} \right) = {\rm{n}}\left( {\rm{x}} \right)\left( {\frac{{\partial \left( {e\left( {y,n\left( x \right)} \right)} \right)}}{{\partial y}}{|_{{\rm{\hat x}}}}} \right) = 0$     （公式4）

　　由上述MLS曲面的隐式定义及投影计算方法可知，MLS曲面实际为局部邻域点集的平滑逼近，具有很好的连续性及局部计算特性。由于单一或部分视角点云数据通常只覆盖工件部分区域，该区域MLS曲面应为非闭合曲面，其边界应为离散边界点集在MLS曲面上投影所定义(边界确定方法:判断此点是否为k邻域点集的凸多边形顶点).

       由上述定义可知,测量完点云的MLS曲面,如果需要进行点云到曲面的ICP对齐,还需要计算空间一点到MLS平面的最近点以构成ICP对齐的精确对应点对.

 

3 点到曲面的正交投影计算

什么是正交投影？

       正交投影就是在曲面上找一点，使该点与空间某一给定点的连线与曲面在该点处的切平面是垂直的。

       隐式曲面（例如MLS曲面）与参数曲面的区别在于，隐式曲面在对水，云。人体肌肉，不规则模具的建模上有很大优势，但是隐式曲面难以控制精度。在对隐式曲面的求解过程中最重要的环节是给出初始投影点的位置，而求出初始投影点的位置可以归结为点到隐式曲面的正交投影问题。

       求点到曲面的正交投影即是求点到曲面的最短距离,首先需要构造初始点处的一条特殊的法截线,并给出沿着该法截线追踪投影点的二阶泰勒展开式，然后将给定点向初始点处的法截线的曲率圆作投影，提出基于曲率圆的步长控制策略，在此基础上给出了基于梯度的迭代误差矫正方法。

贴士 二阶泰勒展开式: n(x) = n(x₀)+n(x₀)^'+n(x₀)^"/2!

具体计算：

3.1 构造法截线与投影点追踪

       设P(m,n,l)是空间中一给定点，又设空间隐式平面S以F(x,y,z) = 0的形式给出，该隐式平面S上任意一点处法向量为n = $\nabla $F/||$\nabla $F||，其中$\nabla {\rm{\;}} = {\rm{\;}}\frac{\partial }{{\partial x}}{\rm{\hat i}} + \frac{\partial }{{\partial y}}{\rm{\hat j}} + \frac{\partial }{{\partial z}}{\rm{\hat k}}$为Hamilton算子。下面考虑隐式曲面S在点g₀=[x₀，y₀，z₀]处的法截面S₀，S₀为隐式曲面在g₀处的单位法向量n₀、向量g₀P以及g₀点所在平面。

       贴士 哈密顿（Hamilton）算子：$\nabla {\rm{\;}} = {\rm{\;}}\frac{\partial }{{\partial x}}{\rm{\hat i}} + \frac{\partial }{{\partial y}}{\rm{\hat j}} + \frac{\partial }{{\partial z}}{\rm{\hat k}}$；Hamilton算子通过这个运算符形成了一个矢量场，该矢量场反映了标量场的分布。

       贴士 偏导和Hamilton算子的关系：以三维空间为例，偏导是曲面上存在某条曲线使其投影在坐标轴垂直就是该轴的偏导数，由Hamilton公式可知三维空间中三个轴偏导的矢量和为Hamilton的矢量表示。

       贴士 Hamilton算子用于表示梯度：首先明确梯度 ($\nabla $F)是一个向量场，标量场上某一点的梯度指向该标量场增长最快的方向，梯度的长度代表当前点的最大变化率，$\nabla {\rm{f\;}} = {\rm{\;}}\frac{{\partial f}}{{\partial x}}{\rm{\hat i}} + \frac{{\partial f}}{{\partial y}}{\rm{\hat j}} + \frac{{\partial f}}{{\partial z}}{\rm{\hat k}}$。

        梯度解释图：

梯度解释图

        如图2所示，法截面S₀的方程为：

[n₀,Pg₀,(x-x₀,y-y₀,z-z₀)] = 0   （公式5）

图2

        如公式5所示，表示为三个向量的混合积为0，若记G(x,y,z)= [n₀,Pg₀,(x-x₀,y-y₀,z-z₀)],则法截面S₀方程记作G(x,y,z)=0；

       贴士 向量混合积:记做(a,b,c),(b, c, a),(c, a, b)；(axb)·c, (bxc)·a, ( cxa)·b称为向量的混合积，即两个向量叉乘的结果和另外一个向量点乘，两个向量点乘为零说明两个向量垂直。

        我们把法截面S₀和隐式曲面S的交线称为法截线，记为C₀，则法截线的方程可以表示为

${C_0} = \left\{ {\begin{array}{*{20}{c}}{F\left( {{\bf{x}},{\bf{y}},{\bf{z}}} \right) = 0}\\{G\left( {{\bf{x}},{\bf{y}},{\bf{z}}} \right) = 0}\end{array}} \right.$     （公式6）

        设s为法截线C₀的弧长参数，${\rm{t}} = \left[ {\frac{{dx}}{{ds}},\frac{{dy}}{{ds}},\frac{{dz}}{{ds}}} \right]$为沿着法截线C₀的单位切向量，用公式6对弧长参数s求导

${C_0} = \left\{ {\begin{array}{*{20}{c}}{F\left( {{\bf{x}},{\bf{y}},{\bf{z}}} \right)|\hat s = t \cdot \nabla F = 0}\\{G\left( {{\bf{x}},{\bf{y}},{\bf{z}}} \right)|\hat s = t \cdot \nabla G = 0}\end{array}} \right.$v    （方程1）   

图3

        其中$\nabla F,\nabla G$分别为隐式曲面F(x,y,z) = 0和G(x,y,z) = 0在点(x,y,z)处的梯度，t为单位矢量，||t|| = 1，求解方程1（向量叉乘的方向满足右手定则）： 

 ${\rm{t}} = \frac{{\nabla F \times \nabla G}}{{\nabla F \times \nabla G}}$

        设${\rm{c}} = \left[ {\frac{{{d^2}x}}{{d{s^2}}},\frac{{{d^2}y}}{{d{s^2}}},\frac{{{d^2}y}}{{d{s^2}}}} \right]$为法截线C₀在点(x,y,z)处的曲率矢量，将方程1对弧长参数s求导得

\[\left\{ \begin{array}{l}t{{\bf{F}}_2}{t^T} + \nabla F \cdot {c^T} = 0\\t{{\bf{G}}_2}{t^T} + \nabla G \cdot {c^T} = 0\\t{c^T} = 0\end{array} \right.\]   (方程2)

        其中

\[{F_2} = \left[ {\begin{array}{*{20}{c}}{{F_{{\rm{xx}}}}}&{{F_{xy}}}&{{F_{xz}}}\\{{F_{yx}}}&{{F_{yy}}}&{{F_{yz}}}\\{{F_{zx}}}&{{F_{zy}}}&{{F_{zz}}}\end{array}} \right]\]，

${G_2} = \left[ {\begin{array}{*{20}{c}}{{G_{{\rm{xx}}}}}&{{G_{xy}}}&{{G_{xz}}}\\{{G_{yx}}}&{{G_{yy}}}&{{G_{yz}}}\\{{G_{zx}}}&{{G_{zy}}}&{{G_{zz}}}\end{array}} \right]$

       求解方程2得到曲率矢量为

 ${\bf{c}} = \left[ { - {\bf{t}}{F_2}{t^T}, - {\bf{t}}{G_2}{t^T},0} \right]{\left( {W_2^T} \right)^{ - 1}}$

       其中${W_2} = \left[ {\begin{array}{*{20}{c}}{{F_{\rm{x}}}}&{{F_y}}&{{F_z}}\\{{G_x}}&{{G_y}}&{{G_z}}\\{{x_s}}&{{y_s}}&{{z_s}}\end{array}} \right]$。

        进一步，可以求出法截线C₀上任意一点处的曲率，公式6定义的法截线C₀在(x,y,z)处的曲率为
${\rm{k}}\left( {{\rm{x}},{\rm{y}},{\rm{z}}} \right) = \frac{{\left\| {\left( {\left( {\nabla F \times \nabla G} \right) \times \nabla \left( {\nabla F \times \nabla G} \right) \times \left( {\nabla F \times \nabla G} \right)} \right)} \right\|}}{{{{\left\| {\nabla F \times \nabla G} \right\|}^3}}}$   (公式7)

        其中，$\nabla F,\nabla G$代表梯度，$\nabla \left( {\nabla F \times \nabla G} \right)$表示对行向量$\nabla F,\nabla G$的每一个分量求梯度，将得到的三维向量作为矩阵的一个列。

        我们通过沿着隐式曲面S上的初始点g₀处的法截线C₀来追踪P点在隐式曲面S上的正交投影点H，但由于正交投影点H不一定在初始点g₀处的法截线C₀上，因此需要不断调整正交投影点的追踪方向，通过反复构造下一个追踪点处的法截线C₀，使得追踪方向始终指向所要求的目标点H。

        具体步骤：

Step1. 由曲面S上g₀初始点处出发，沿着法截线C₀追踪到点g=[x,y,z].

Step2. 判断点g是否是所需要的正交投影点H，若是，执行下一步；否则，将g作为新的初始点g₀，并构造曲面S在该点处的法截线C₀，转Step1.

Step3. 输出H的值为P点的正交投影坐标，算法结束。

 

---

        如何从g0处出发，沿着法截线追踪到下一个点g(x,y,z)？

---

        如图4所示，由于点g在法截面S₀上，同时也在法截线C₀上，因此向量g₀g可以表示成法截面S₀上的2个线性无关的向量的线性组合。由于法截线C₀在点g₀的单位切向量t₀以及曲率向量c₀相互垂直，并且向量t₀以及曲率向量c₀均位于法截面S₀内，因此向量g₀g可以用上述两个向量线性表示出。利用二阶泰勒公式，g=[x,y,z]可以表示为

${\bf{g}} = {g_0} + {t_0}\Delta s + \frac{1}{2}{c_0}\Delta {s^2}$   （公式8）

        当给出迭代步长Δs时，可以计算出新的点g。

图4 

3.2 如何根据曲率向量设置步长 

        隐式曲面S在正交投影点H处的法向量n与向量PH共线，同时向量PH与经过H点处的任意一条法线是正交的，同时也与法截线在H点处的曲率圆是正交的。若步长Δs恒定时，上述迭代算法并不能保证经过若干次迭代后能得到所要求的目标值H，需要进一步研究步长的问题，即能采用适当的步长控制策略，使得按照公式8计算得到的点能最终达到目标值。在此提出一种基于曲率的变步长控制方法，保证在追踪点g不断逼近H点的过程中，步长能逐渐趋于0，最终找到目标H。

        如图5所示。初始点g₀处的曲率圆、法截线C₀以及空间点P均在法截面S₀上，设法截线C₀在g₀点处的曲率圆的圆心O的坐标为(O_x0, O_y0, O_z0)，则有

(O_x0, O_y0, O_z0) = (x₀, y₀, z₀) + β₀/k₀

        其中β₀=c₀/||c₀||为法截线C’₀在g₀处的主法向量，c₀为g₀处的曲率矢量；k0为法截线C₀在g₀处的曲率。由于该曲率圆位于法截面S₀上，故曲率圆的方程为

 $\left\{ {\begin{array}{*{20}{c}}{{{\left( {X - {O_{x0}}} \right)}^2} + {{\left( {Y - {O_{y0}}} \right)}^2} + {{\left( {Z - {O_{z0}}} \right)}^2} = \frac{1}{{k_0^2}}}\\{\left[ {{n_0},P{g_0},\left( {X - {x_0},Y - {y_0},Z - {z_0}} \right)} \right] = 0}\end{array}} \right.$   (方程3)

图5 点P向法截线在g₀处的曲率圆做投影

        考虑将P点向初始点g₀处的曲率圆作正交投影，设垂足为Q。可以证明：直线PO与方程3所表示的曲率圆交点即为点P在曲率圆上的正交投影点Q。由于直线PO与曲率圆有两个交点Q1和Q2，我们选择离g₀较近的交点作为P点在曲率圆上的投影点，即

\[{\rm{Q}} = \left\{ {\begin{array}{*{20}{c}}{{Q_1},\;\;if\left\| {{g_0}{Q_1}} \right\| < \left\| {{g_0}{Q_2}} \right\|}\\{{Q_2},\;\;\;\;\;\;\;\;else\;}\end{array}} \right.\]

继续学习请转至我的下一篇博客 ：  基于移动最小二乘曲面的点云对齐（二） ——  优化点云对齐迭代点

 

 

【未完待续~】