本文重點
這次主要介紹一種點雲對齊的方法,多視數據最近迭代(ICP)對齊是最常用的點雲對齊方法,為了提高對齊的精度及穩定性我們使用一種基於移動最小二乘(MLS)曲面的ICP多視數據對齊方法.該方法無需對數據進行額外的去噪和數據分割.對於優化噪聲點的點雲對齊可以采用本方法進行點雲對齊。
1 多視數據相關概念
因為掃面儀(傳感器)的自身掃描范圍限制,工業掃描往往不能一次對完整的模具進行掃描,需要通過多次轉換視角進行模具的掃描.由於各視角的點雲數據定義於各自的局部坐標系下,我們就需要對這些多視角測量數據進行坐標系轉換使其在統一的坐標系下.
什么是點雲對齊？
點雲對齊就是通過坐標平移或旋轉將具有相似特征的一組或幾組點雲轉換到統一的全局坐標系下，多視角點雲對齊通常有三類方法:
- 利用度量裝置自動記錄視角間變換,然后通過一組數據作為標准數據,其余數據逆向轉換回全局坐標系;
- 被測模具表面貼足夠的標識點，通過計算標志點之間的拓撲關系進行點雲間的對齊;
- 測量時使不同視角數據擁有足夠的數據重疊，通過重疊部分的坐標變換關系，將局部坐標系下的點雲轉換到全局坐標系下。
2 點雲對齊的方法
直接通過設備記錄點雲運動的方法成本太高，需要一台記錄裝置，本文主要針對目前使用最廣的先進行手動粗對齊，在用IPC算法精對齊。
移動最小二乘（Moving Least-Squares,MLS）曲面逼近方法，不僅能適應任意拓撲形狀的點雲估計，而且逼近的MLS曲面具有良好的幾何性質，可以通過確定逼近誤差閾值排除噪點。
什么是MLS曲面？
國外研究者Levin將MLS曲面定義為滿足對投影映射具有不變性的空間點集合，即
S(x) = {x∈R3|f(x) = x}
Amenta在次基礎上對MLS曲面投影進行了更精確的定義。即沿法矢量場n(x)滿足能量函數e(y,a)的極小值點集合（y為空間點，a為投影方向），並給出了由空間一點x0到MLS曲面的迭代投影計算方法。
如圖1所示，從初始點xo開始沿加權法矢量n(xo)優化能量函數e(y,a)得到點x1，其中能量函數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中的xi）的鄰域點集{pi}的個數（一般k取15）；a取輸入點y處的加權法矢量n(xi),其中n(xi)可通過鄰近點的加權平均得到，即
${\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中vj = ɵ(||xi-pi||)以及公式1中的vj = ɵ(||y-pi||)均為高斯函數，表示為
$\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(x0)+n(x0)'+n(x0)"/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在點g0=[x0，y0，z0]處的法截面S0，S0為隱式曲面在g0處的單位法向量n0、向量g0P以及g0點所在平面。
貼士 哈密頓（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所示，法截面S0的方程為：
[n0,Pg0,(x-x0,y-y0,z-z0)] = 0 （公式5）
圖2
如公式5所示，表示為三個向量的混合積為0，若記G(x,y,z)= [n0,Pg0,(x-x0,y-y0,z-z0)],則法截面S0方程記作G(x,y,z)=0；
貼士 向量混合積:記做(a,b,c),(b, c, a),(c, a, b)；(axb)·c, (bxc)·a, ( cxa)·b稱為向量的混合積，即兩個向量叉乘的結果和另外一個向量點乘，兩個向量點乘為零說明兩個向量垂直。
我們把法截面S0和隱式曲面S的交線稱為法截線，記為C0，則法截線的方程可以表示為
${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為法截線C0的弧長參數，${\rm{t}} = \left[ {\frac{{dx}}{{ds}},\frac{{dy}}{{ds}},\frac{{dz}}{{ds}}} \right]$為沿着法截線C0的單位切向量，用公式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]$為法截線C0在點(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]$。
進一步，可以求出法截線C0上任意一點處的曲率，公式6定義的法截線C0在(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上的初始點g0處的法截線C0來追蹤P點在隱式曲面S上的正交投影點H，但由於正交投影點H不一定在初始點g0處的法截線C0上，因此需要不斷調整正交投影點的追蹤方向，通過反復構造下一個追蹤點處的法截線C0，使得追蹤方向始終指向所要求的目標點H。
具體步驟：
Step1. 由曲面S上g0初始點處出發，沿着法截線C0追蹤到點g=[x,y,z].
Step2. 判斷點g是否是所需要的正交投影點H，若是，執行下一步；否則，將g作為新的初始點g0，並構造曲面S在該點處的法截線C0，轉Step1.
Step3. 輸出H的值為P點的正交投影坐標，算法結束。
如圖4所示，由於點g在法截面S0上，同時也在法截線C0上，因此向量g0g可以表示成法截面S0上的2個線性無關的向量的線性組合。由於法截線C0在點g0的單位切向量t0以及曲率向量c0相互垂直，並且向量t0以及曲率向量c0均位於法截面S0內，因此向量g0g可以用上述兩個向量線性表示出。利用二階泰勒公式，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所示。初始點g0處的曲率圓、法截線C0以及空間點P均在法截面S0上，設法截線C0在g0點處的曲率圓的圓心O的坐標為(Ox0, Oy0, Oz0)，則有
(Ox0, Oy0, Oz0) = (x0, y0, z0) + β0/k0
其中β0=c0/||c0||為法截線C’0在g0處的主法向量，c0為g0處的曲率矢量；k0為法截線C0在g0處的曲率。由於該曲率圓位於法截面S0上，故曲率圓的方程為
$\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向法截線在g0處的曲率圓做投影
考慮將P點向初始點g0處的曲率圓作正交投影，設垂足為Q。可以證明：直線PO與方程3所表示的曲率圓交點即為點P在曲率圓上的正交投影點Q。由於直線PO與曲率圓有兩個交點Q1和Q2，我們選擇離g0較近的交點作為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.\]
繼續學習請轉至我的下一篇博客 ： 基於移動最小二乘曲面的點雲對齊（二） —— 優化點雲對齊迭代點
【未完待續~】