PCL MLS論文Computing and Rendering Point Set Surfaces研讀筆記

前言

做點雲處理最有名的庫當屬PCL了。PCL的MLS(Moving Least Squares)模塊用於表面平滑，而其理論基礎便是Computing and Rendering Point Set Surfaces這篇論文。本篇主要探討論文裡的3 Defining the surface - projecting 及5 Generating the representation point set 這兩部份。

本文著重理解算法各步驟的前後關係，並補全了論文當中省略的公式推導。文章中螢光標記的部份表示還沒完全明白，需要再仔細研究。

定義曲面 - 投影

投影步驟

MLS對點雲做平滑，實際上在做的事情是對點雲中的每個點（下文中的$r$），都用它附近的點（下文的 { p 1 , . . . , p N } \{p_1, ..., p_N\} {p1​,...,pN​}）擬合一個曲面（下文的$S_p$），然後再將該點投影到曲面上。但是在這篇論文中，並不是直接尋找曲面，而是先在局部擬合一個"平"面（下文的$H$），然後才去尋找曲面。

尋找點$r$的擬合平面(或說座標系)$H$

論文裡對平面H的描述如下： H = { x ∣ < n , x > − D = 0 , x ∈ R 3 } , n ∈ R 3 , ∥ n ∥ = 1 H = \{x | \left<n, x\right> - D = 0, x \in R^3\}, n \in \R^3, \|n\| = 1 H={x∣⟨n,x⟩−D=0,x∈R3},n∈R3,∥n∥=1。乍看之下有點不明所以。

我們可以先從我們熟悉的平面方程上手：$a{x}_1+b{x}_2+c{x}_3+d=0$。論文裡寫的$⟨n,x⟩-D=0$與上式是其實同一個意思，只要把$n$這個法向量看成$(a,b,c)$，$D$看成$-d$即可。觀察$⟨n,x⟩-D=0$，我們可以發現，一個平面是可以由其法向量及其到原點的距離來定義的(詳見Distance from point to plane)。

要用點$r$附近的點擬合一個平面，我們有無數種選擇，因此這裡我們還需要一個cost function，用來指引什麼樣的平面是我們所偏好的。MLS希望找到的平面用文字說明如下：希望找到一個平面$H$，使得集合 { p 1 , . . . , p N } \{p_1, ..., p_N\} {p1​,...,pN​}裡的所有點$p_i$到平面$H$的距離平方和是最小的。另外，每個點還會被進行加權，離點$r$越近的點權重越高，反之則越低。到此，已經明確了我們的目標：最小化加權距離平方和。

以下即為論文裡的方程式(1)，代表尋找$H$時用的cost function:

$${\Sigma }_{i=1}^N(⟨n,p_i⟩-D{)}^2\theta (\Vert p_i-q\Vert )$$

其中$q$為點$r$在平面$H$上的投影，即$r$沿平面法向量$n$方向移動到平面上的點，它們的關係可以用$q=r+tn$描述，其中$t$為一個常數。

方程式(1)中的$(⟨n,p_i⟩-D{)}^2$代表點$p_i$到平面$H$的距離的平方；$\theta (\Vert p_i-q\Vert )$則是點$p_i$的權重，與$p_i$到$q$的距離成反相關，即距離越遠權重越低。

cost function還有另一種型式，即論文裡的方程式(2)：

$${\Sigma }_{i=1}^N{⟨n,p_i-r-tn⟩}^2\theta (\Vert p_i-r-tn\Vert )$$

$\theta (\Vert p_i-r-tn\Vert )$這一項很好理解，就是把$q$替換成$r+tn$而已。那麼${⟨n,p_i-r-tn⟩}^2$要怎麼理解呢?

先把${⟨n,p_i-r-tn⟩}^2$中的$r+tn$替換回$q$，得到${⟨n,p_i-q⟩}^2=(⟨n,p_i⟩-⟨n,q⟩{)}^2$。因為點$q$在平面上，故$⟨n,q⟩=D$。代回上式：$(⟨n,p_i⟩-⟨n,q⟩{)}^2=(⟨n,p_i⟩-D{)}^2$，得到了與方程式(1)相同型式。因此方程式(2)與方程式(1)是等價的。

cost function中，點$p_i$的權重$\theta (\Vert p_i-r-tn\Vert )$是將$p_i$與$q$的距離代入一個平滑、單調遞減、恆為正值的函數$\theta$中所得到，我們通常選取高斯函數作為$\theta$。另外，這裡代入$\theta$的是$p_i$與$q$的距離而非$p_i$與$r$的距離，這在3.2節會詳述。

注意到當我們把平面$H$定義好了，我們同時也定義了局部座標系，這個座標系的原點即點$q$。接下來在尋找曲面時，就會用到這個局部座標系。

尋找由$r$定義的曲面$S_p$

在這篇論文中，是用一個雙變量的多項式來表達一個曲面的。為何是雙變量呢?這是因為一個曲面有兩個自由度。曲面的數學式寫為$g(x,y)$，為一純量函數。

來看看投影是怎麼進行的：對點$r$來說，它在平面$H$上的投影為$q=r+tn$，在曲面$S_p$上的投影則為$\mathcal{P}(r)=q+g(0,0)n=r+(t+g(0,0))n$。觀察$\mathcal{P}(r)=q+g(0,0)n$，可以發現$$\mathcal{P}(r)$是局部座標系原點$q$沿$n$方向移動$g(0,0)$距離所得到的。我們同時也可以發現，$r$，$q$，$\mathcal{P}(r)$三點位於同一直線上。

對$r$鄰域內的點$p_i$，其在局部座標系下的座標為$(x_i,y_i,f_i)$，其中$f_i=n\cdot (p_i-q)$為點$p_i$在局部座標系下的z座標。點$p_i$在平面$H$上的投影為$q_i=(x_i,y_i,0)$，注意投影後的z座標為0。在曲面$S_p$上的投影為$\mathcal{P}(p_i)=q_i+g(x_i,y_i)n=(x_i,y_i,g(x_i,y_i))$。注意不論是投影到平面或曲面上，x座標及y座標都不會改變。

到了這裡，我們就看得懂論文裡公式(3):

$${\Sigma }_{i=1}^N(g(x_i,y_i)-f_i{)}^2\theta (\Vert p_i-q\Vert )$$

在寫什麼了：它定義了一個cost function，希望找到一個曲面$S_p$(或$g$)，使得$p_i$投影到曲面上之後，能與原來的$p_i$越接近越好(只關注它們的z座標)。公式(3)的第二項即為權重，其定義跟找$H$時相同。

權重函數通常使用高斯函數，即論文裡的公式(4):

$$\theta (d)=e^{-\frac{d^2}{h^2}}$$

其中$d$代表兩點間的距離，$h$則是一個參數，反映了點雲中相鄰點距離的期望值。我們可以透過調整$h$來決定平滑的強度：小的$h$會導致高斯函數下降得很快，所以在平滑時只會考慮附近一小部份的點，所以平滑效果可能不明顯；大的$h$則會考慮較大範圍的點來做平滑，給出較強的平滑效果。

上述"尋找$H$之後尋找$S_p$的過程可由論文中的圖三總結：

projection property

這一小節提到了projection property：$\mathcal{P}(\mathcal{P}(r))=\mathcal{P}(r)$，用白話來說就是投影函數$\mathcal{P}$對點$r$投影一次跟投影多次都會在同一個位置上。

計算權重時用的是$p_i$與$q$的距離而非$p_i$與$r$的距離，這裡說如果使用的是$p_i$與$r$的距離，那麼$\theta$的值將會在法向量方向上變動，這將導致$\mathcal{P}$不再會是一個projection。 $\to$這句話沒看懂，$\theta$的值不是本來就會在法向量方向上變動? 猜想這句話的意思是：假想$r$點沿法向量方向移動，那麼$p_i$與$r$的距離將會隨之變動，反之$p_i$與$q$的距離則不會。我們希望對於"$r$點沿法向量方向移動所得到的直線"上的每一點，與某個$p_i$計算出來的$\theta$值都是相同的。(還是不懂這跟projection property的聯繫)

在前一小節不論是定義$H$或$S_p$的cost function時，裡面都有一項$\theta$權重函數，這裡建議的權重函數為高斯函數：$\theta (d)=e^{(-\frac{d^2}{h^2})}$。

其參數$d$如上面所看到的，應該代入$p_i$與$q$的距離。另外還有個參數$h$，可以把它想成是標準差的角色。$h$越小權重下降得越快，相當於用來擬合的點變少，所以得到的曲面細節會越豐富；反之當$h$越來越大時，權重下降得越慢，相當於用來擬合的點變多，所以曲面會越平滑。

取小$h$與大$h$所產生的平滑效果如論文圖四：

實際計算投影

尋找平面H

觀察方程式(2)，其中$r$與 { p 1 , p 2 , . . . p N } \{p_1,p_2,...p_N\} {p1​,p2​,...pN​}皆是給定的輸入，只有$t$與$n$是未知數，我們的目標便是找到一組$n$跟$t$，使得方程式(2)能最小化。方程式(2)有多個局部最小值，我們要找$t$最小的一個，也就是希望$q$離$r$越近越好(因為 $q=r+tn$)。

此處的思想是：首先估計點$r$的法向量$n$，然後迭代優化$t$與$q$使損失函數最小化。

注：感謝@veryneo的指正，原文的第一到三步已修正成如下形式。

Initial normal estimate in r

首先注意這一步是可選的，如果輸入的點雲已包含法向量，就不需要這一步。

這一步是固定$t=0$，來求法向量$n$的近似(注意$n$的長度被固定為1，所以這一步估計的是法向量$n$的方向)。

將$t=0$代入方程式(2)，$tn$這項消失，並且$\theta$變為$p_i$的函數，跟$n$無關，在此記為${\theta }_i$，是一個常數。方程式(2)至此化簡為${\Sigma }_{i=1}^N{⟨n,p_i-r⟩}^2{\theta }_i$。

接下來是將它寫成向量與矩陣相乘的形式，更正式的名稱為bilinear form。

為了推導出方程式(2)的bilinear form，我們先只關注$i$這一項: ${⟨n,p_i-r⟩}^2{\theta }_i$。

將它展開：

$\begin{array}{rl}{⟨n,p_i-r⟩}^2{\theta }_i& =(n_1(p_{i1}-r_1)+n_2(p_{i2}-r_2)+n_3(p_{i3}-r_3){)}^2{\theta }_i\\ & =((p_{i1}-r_1{)}^2{n}_1^2+(p_{i1}-r_1)(p_{i2}-r_2)n_1{n}_2+(p_{i1}-r_1)(p_{i3}-r_3)n_1{n}_3\\ & +(p_{i2}-r_2)(p_{i1}-r_1)n_2{n}_1+(p_{i2}-r_2{)}^2{n}_2^2+(p_{i2}-r_2)(p_{i3}-r_3)n_2{n}_3\\ & +(p_{i3}-r_3)(p_{i1}-r_1)n_3{n}_1+(p_{i3}-r_3)(p_{i2}-r_2)n_3{n}_2+(p_{i3}-r_3{)}^2{n}_3^2){\theta }_i\end{array}$

到了這裡已經能看出矩陣的雛形了，接下來是把向量跟係數矩陣找出來，寫成向量與矩陣相乘的形式。

我們要找的向量自然是$n$，所以在找係數矩陣時，可以先暫時忽略$n$。找出的係數矩陣為：

$\left[\begin{array}{ccc}(p_{i1}-r_1{)}^2& (p_{i1}-r_1)(p_{i2}-r_2)& (p_{i1}-r_1)(p_{i3}-r_3)\\ (p_{i2}-r_2)(p_{i1}-r_1)& (p_{i2}-r_2{)}^2& (p_{i2}-r_2)(p_{i3}-r_3)\\ (p_{i3}-r_3)(p_{i1}-r_1)& (p_{i3}-r_3)(p_{i2}-r_2)& (p_{i3}-r_3{)}^2\end{array}\right]{\theta }_i$

得到$n$跟係數矩陣後可以自行驗證一下，$n^T[...]n={⟨n,p_i-r⟩}^2{\theta }_i$應該是成立的。

記得剛剛我們只看了$i$這一項，但我們要的應該是${\Sigma }_i^N$，因此還要對上面的係數矩陣中的各項都加上${\Sigma }_i^N$。到了這裡，應該可以很容易地看出，我們計算得到的矩陣跟論文裡的$B$矩陣是一致的： B = { b j k } , B ∈ R 3 × 3 , b j k = Σ i θ i ( p i j − r j ) ( p i k − r k ) B = \{b_{jk}\}, B \in \R^{3 \times 3}, b_{jk} = \Sigma_i \theta_i ({p_i}_j - r_j)({p_i}_k - r_k) B={bjk​},B∈R3×3,bjk​=Σi​θi​(pi​j​−rj​)(pi​k​−rk​)

我們可以來看一下矩陣$B$是由哪些元素組成的，也就是： { p 1 , . . . p N } \{p_1,...p_N\} {p1​,...pN​},$r$及 { θ 1 , . . . θ N } \{\theta_1,...\theta_N\} {θ1​,...θN​}。我們已經知道${\theta }_i$是由$p_i$及$r$決定的常數，所以矩陣$B$實際上由 { p 1 , . . . p N } \{p_1,...p_N\} {p1​,...pN​}及$r$決定，而這些向量全是已知的，所以矩陣$B$也是已知的。

向量$n$及矩陣$B$都找到了之後，就可以將方程式(2)寫成bilinear form:$n^{T}Bn$，回想方程式(2)是一個cost function，所以我們希望找到一個$n$使得該式算出來的值最小化，於是我們的目標變成了：

$$mi{n}_{\Vert n\Vert =1}{n}^{T}Bn$$

此最小化問題的解為矩陣$B$最小特徵值所對應的特徵向量。原因為何?

先來回顧一下矩陣特徵值及特徵向量的物理意義：如果一個向量經過矩陣的變換後方向不變，則稱該向量為特徵向量，其縮放程度則為該特徵向量對應的特徵值。我們知道矩陣$B$對"其最小特徵值對應的特徵向量"所進行的放大程度最小，所以這裡我們選擇該特徵向量的方向作為第一步的解$n$。

回頭看一下矩陣$B$，它其實是一個共變異數/協方差矩陣，這是為什麼呢?

我們可以把$r$看作$p_i$的期望值，即$E[p_i]$。那麼${p_i}_j$與${p_i}_k$的共變異數$cov[{p_i}_j,{p_i}_k]=E[({p_i}_j-r_j)({p_i}_k-r_k)]$，與矩陣$B$的元素$b_{jk}={\Sigma }_i{\theta }_i({p_i}_j-r_j)({p_i}_k-r_k)$兩相對照，發現兩者十分吻合，只要我們將${\Sigma }_i{\theta }_i$這個加權平均的動作看成是取期望值即可。既然知道$b_{jk}=cov[{p_i}_j,{p_i}_k]$，那麼我們也就知道$B$是由$p_i$任兩分量的變異數所組成的3*3共變異數矩陣!

Iterative non-linear minimization

接下來是交替優化$t$與$q$，為什麼不是優化$t$與$n$呢?論文中說道，在$t$固定的情況下，其實我們本來是應該是要優化$n$，相當於在一個球面上尋找更好的$q$。但是為了使用conjugate gradient method來進行優化，所以這裡用的是球面的近似，也就是只在平面$H$上尋找$q$。

固定n求t和固定t優化q這兩步會持續進行，直到所有的參數在某一輪的變化量都小於一個事先定義的閾值$ϵ$才結束。

固定n求t

在$[-h/2,h/2]$的範圍內找一個使損失函數最小化的$t$，計算方法同第二步。

為了求解方便，在論文中限制$t$的範圍為$[-h/2,h/2]$，因為在此處$t$只會有一個局部極小值。 $\to$為什麼?

因為我們知道$t$只有一個局部極小值，所以這一步可以直接設方程式(2)對$t$的偏微分等於0來求出。

欲求${\Sigma }_{i=1}^N{⟨n,p_i-r-tn⟩}^2\theta (\Vert p_i-r-tn\Vert )$對$t$的微分，同樣先只關注$i$這一項，令$A=C^2={⟨n,p_i-r-tn⟩}^2$，所以目前我們所關注的是$A\theta$，其中$A$及$\theta$都是$t$的函數($n$在這一步固定)。

$i$這一項的微分：$\frac{d(A\theta )}{dt}=\frac{dA}{dt}\theta +A\frac{d\theta }{dt}$。

• 先看第一項：$\frac{dA}{dt}=\frac{dA}{dC}\frac{dC}{dt}$，其中$\frac{dA}{dC}=2C$。

  要計算$\frac{dC}{dt}$，我們可以先把$C$拆開來看看，即：

  $\begin{array}{rl}C& =⟨n,p_i-r-tn⟩\\ & =n_1({p_i}_1-r_1-t{n}_1)+n_2({p_i}_2-r_2-t{n}_2)+n_3({p_i}_3-r_3-t{n}_3)\\ & =-(n_1^2+n_2^2+n_3^2)t+⟨n,p_i-r⟩\\ & =-t+⟨n,p_i-r⟩\end{array}$

  所以$\frac{dC}{dt}=-1$。

  所以微分第一項$\frac{dA}{dt}\theta =\frac{dA}{dC}\frac{dC}{dt}\theta =2C\ast (-1)\ast \theta =-2C\theta$，算出來的跟方程式(6)差一個負號?

• 再看第二項：$A\frac{d\theta }{dt}$

  其中$\theta (t)=e^{(-\frac{d^2}{h^2})}$，$d=\Vert p_i-r-tn\Vert$

  $\theta (t)$對$t$的微分：$\frac{d\theta }{dt}=\frac{d\theta }{d(d^2)}\frac{d(d^2)}{d(d)}\frac{d(d)}{dt}$

  其中$\frac{d\theta }{d(d^2)}=\theta \ast (-\frac{1}{h^2})$，$\frac{d(d^2)}{d(d)}=2d$

  $\frac{d(d)}{dt}$的計算較為複雜。我們先將d化簡：

  $$\begin{gathered} d=\Vert p_i-r-tn\Vert \\ =\sqrt{({p_i}_1-r_1-t{n}_1{)}^2+({p_i}_2-r_2-t{n}_2{)}^2+({p_i}_3-r_3-t{n}_3{)}^2} \\ =\sqrt{(n_1^2+n_2^2+n_3^2)t^2-2[n_1({p_i}_1-r_1)+n_2({p_i}_2-r_2)+n_3({p_i}_3-r_3)]t+({p_i}_1-r_1{)}^2+({p_i}_2-r_2{)}^2+({p_i}_3-r_3{)}^2} \\ =\sqrt{t^2-2⟨n,p_i-r⟩t+({p_i}_1-r_1{)}^2+({p_i}_2-r_2{)}^2+({p_i}_3-r_3{)}^2} \end{gathered}$$

  $\frac{d(d)}{dt}=\frac{1}{2d}(2t-2⟨n,p_i-r⟩)=\frac{t-⟨n,p_i-r⟩}{d}=\frac{-C}{d}$

  所以微分第二項$A\ast \frac{d\theta }{dt}=A\ast \theta \ast (-\frac{1}{h^2})\ast 2d\ast (\frac{-C}{d})=2A\theta \frac{C}{h^2}$

算出cost function對$t$的一階導數為$-2C\theta +2A\theta \frac{C}{h^2}=2C(-1+\frac{A}{h^2})\theta$，對所有$i$的偏導數做總和，得到論文中的公式(6)：

$$2{\Sigma }_{i=1}^N⟨n,p_i-r-tn⟩(1+\frac{{⟨n,p_i-r-tn⟩}^2}{h^2})e^{\frac{\Vert p_i-r-tn{\Vert }^2}{h^2}}$$

(兩者差了一個負號?)

之後，可以使用Numerical
Recipes in C: The Art of Scientific Computing (2nd ed.) - 10.3 One-Dimensional Search with FirstDerivatives裡提到的迭代的最小化方法來求解。

固定t優化q

使用conjugate gradient method在平面$H$上找一個更好的$q$。$\leftarrow$這樣子$q=r+tn$還會成立?

計算曲面Sp

回顧一下計算曲面的cost function：${\Sigma }_{i=1}^N(g(x_i,y_i)-f_i{)}^2\theta (\Vert p_i-q\Vert )$。

$H$找到後，$n$跟$t$都固定下來，所以cost function裡的${\theta }_i$及$f_i$也都變為常數。計算cost function對曲面多項式中各係數的一階導數，可以得到一個線性方程組，其方程式數量正好等於多項式係數的個數。這個問題可以使用linear least squares的方法來求解。

生成具代表性的點雲

下採樣

下採樣的核心思想是計算各點的貢獻度，然後逐一移除貢獻度最低的點，直到停止條件被滿足為止。

貢獻度的計算:

• 理想：用$P$擬合一個曲面$S$，再用將$p_i$移除後的點集 { P − p i } \{P-p_i\} {P−pi​}擬合一個曲面 S { P − p i } S_{\{P-p_i\}} S{P−pi​}​，然後計算兩曲面的Hausdorff distance。但是這種做法的計算量太大，因此論文裡採用以下做法。

• 現實：將$p_i$投影到 S { P − p i } S_{\{P-p_i\}} S{P−pi​}​。(然後呢?論文沒說後續怎麼做)猜想是把加上$p_i$投影點的 S { P − p i } S_{\{P-p_i\}} S{P−pi​}​當作$S$，然後計算它們的Hausdorff distance。

得到每個點的貢獻度之後，將點存到一個priority queue裡面。在迭代過程中，每次會移除貢獻度最低的點。當一個點被移除後，我們需要重新計算該點附近的點的貢獻度。

以上過程會一直重複，直到點的數量小於某一閾值或者是各點的貢獻度都大於某一閾值後才停止。

上採樣

上採樣的核心思想是在MLS曲面上畫vonoroi diagram。畫法如下：取平面上任意三個點畫圓，圓要穿過這三個點(但實際上可能穿過更多的點)，並要注意不能讓圓內包含其它的點。這些圓的圓心即為vonoroi diagram的頂點。得到vonoroi diagram之後，我們會挑選擁有最大半徑的圓的圓心，將它投影到MLS曲面上，得到的結果就是一個上採樣的點。

vonoroi diagram的畫法：

• 理想：在MLS曲面上畫vonoroi diagram。
• 現實：隨機從點雲裡挑選一個點，在點的鄰域內擬合一個"平面"，將鄰域內的點都投影到該平面上。然後在該平面上畫vonoroi diagram。

上採樣的過程會一直重複，直到最小圓的半徑小於一個事先定義好的值才停止。