泊松表面重建算法原理

2022-10-16 00:18 作者:生醫(yī)小王子 0人讀過 | 我要投稿

一、概述

泊松表面重建算法（Poisson Surface Reconstruction）[1]是M. Kazhdan等人于2006年提出的點云表面重建算法，與基于徑向基函數(shù)（Radial Basis Function，RBF）[2]的移動立方體算法（Marching Cube，MC）相似，都是先定義一個隱函數(shù)來表示表面（Surface），根據(jù)點云數(shù)據(jù)求解隱函數(shù)的具體表達(dá)式，然后使用MC方法將表面三角化。但有2點不同，首先使用的隱函數(shù)不同，其次使用了八叉樹來替代3D網(wǎng)格。

隨后，M. Kazhdan等人在2013年提出了屏蔽泊松表面重建算法（Screened Poisson Surface Reconstruction）[3]，他們發(fā)現(xiàn)原方法的表面過于平滑了，所以改進(jìn)了算法，以保證重建表面更貼近于點云。PCL中實現(xiàn)的泊松重建算法即為該改進(jìn)算法。

泊松算法的輸入為帶法向量和空間坐標(biāo)的點云數(shù)據(jù)，輸出為用三角網(wǎng)格表示的表面，其中法向量應(yīng)該保證是準(zhǔn)確的（全都指向物體內(nèi)部或者全指向外部）。在下面的第2和第3章將分別介紹算法的推導(dǎo)過程和實現(xiàn)方法。第4章將介紹屏蔽泊松表面重建算法的不同。

注意：由于我對算法原理也是一知半解，僅知道算法大致流程，所以大家要理解泊松重建算法還是要看原論文。

（專欄有100圖片限制，然后公式也算圖片，也是服了。。）

二、算法推導(dǎo)過程

算法大致思路：找到一個隱函數(shù)表示表面，然后使用移動立方體等表面提取算法來三角化表面

對于物體M，文章定義了一個3D連續(xù)分布的指示函數(shù) $%5Cchi%20_%7BM%7D$ ，物體內(nèi)部的函數(shù)值為1，物體外部函數(shù)值為0，如圖1所示。若能求得該指示函數(shù)，通過MC方法可以很方便地計算出物體表面。
點云法向量是指示函數(shù)梯度的采樣

對指示函數(shù) $%5Cchi%20_%7BM%7D$ 求導(dǎo)計算 $%5Cnabla%7B%5Cchi%20_M%7D$ ，由于物體內(nèi)部和外部區(qū)域都是常數(shù)，這些區(qū)域的梯度為0，僅在表面處有數(shù)值，所以 $%5Cnabla%7B%5Cchi%20_M%7D$ 為指向物體內(nèi)部的表面法向量 $%5Cvec%20V$ 。當(dāng)規(guī)定外部輸入的點云法向量都是指向物體內(nèi)部時，這些點云法向量是 $%5Cwidetilde%20%5Cchi_M$ 的離散采樣，應(yīng)該是有辦法根據(jù)這些離散采樣值估計出 $%5Cvec%20V$ 的分布的。指示函數(shù) $%5Cchi%20_%7BM%7D$ 應(yīng)該保證下面能量函數(shù)是最小的：

$E(%5Cchi)%3D%5Cint%20%7B%5CVert%20%5Cvec%20V-%5Cnabla%5Cchi(p)%5CVert%7D%5E2dp$ ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? (1)
確定表面法向量 $%5Cvec%20V$ 的表達(dá)式
由于 $%5Cchi%20_%7BM%7D$ 是分段常量函數(shù)， $%5Cvec%20V$ 在表面處的數(shù)值是無限大的，不好用點云的法向量表示，故退而求其次，使用平滑濾波器對 $%5Cchi%20_%7BM%7D$ 做卷積，得到平滑的指示函數(shù) $%5Ctilde%20%5Cchi_M$ 。其中 $%5Ctilde%20%5Cchi_M$ 的梯度為 $%5Cchi%20_%7BM%7D$ 的表面法向量平滑后的結(jié)果，如下式所示：
$%5Cnabla(%5Cchi_M*%5Cwidetilde%20F)(q)%3D%5Cint_%7B%5Cpartial%20M%7D%5Cwidetilde%20F_p(q)%5Cvec%20N_%7B%5Cpartial%20M%7D(p)dp$ ? ? ? ? ? ? ? ? ? ?? ? ? ? ?(2)
其中， $%5Cpartial%20M$ 為物體M的表面； $%5Cvec%20N_%7B%5Cpartial%20M%7D(p)$ 為表面的內(nèi)法向量， $p%5Cin%20%5Cpartial%20M$ ； $%5Cwidetilde%20F(q)$ 為平滑濾波器，且滿足 $%5Cwidetilde%20F(q)%3D%5Cwidetilde%20F(q-p)$ ，即模糊權(quán)值與兩點之間的距離有關(guān)。
假設(shè)點云S將表面 $%5Cpartial%20M$ 分成了不同的小塊 $%5CPsi_s$ ，可以使用數(shù)據(jù)點的法向量近似表示 $%5CPsi_s$ 處的整體法向量，故上式（2）可以近似表示為：
$%5Cnabla(%5Cchi_M*%5Cwidetilde%20F)(q)%3D%5Cint_%7B%5Cpartial%20M%7D%5Cwidetilde%20F_p(q)%5Cvec%20N_%7B%5Cpartial%20M%7D(p)dp%5Capprox%20%5Csum_%7Bs%5Cin%20S%7D%7C%5CPsi_s%7C%5Cwidetilde%20F_%7Bs.p%7D(q)s.%5Cvec%20N%5Cequiv%20%5Cvec%20V(q)$ ?(3)
其中， $%7C%5CPsi_s%7C$ 為小塊的面積，s.p和 $s.%5Cvec%20N$ 分別為數(shù)據(jù)點的位置和法向量。 $%5Cwidetilde%20F_%7Bs.p%7D(q)$ 可以選用高斯函數(shù)。當(dāng)點云均勻分布時， $%7C%5CPsi_s%7C$ 為常數(shù)，該常數(shù)不影響最終結(jié)果（后面會說明），故上式（3）可以簡化為：
$%5Cvec%20V(q)%3D%5Csum_%7Bs%5Cin%20S%7D%5Cwidetilde%20F_%7Bs.p%7D(q)s.%5Cvec%20N$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? (4)
注意： $%5Cvec%20V(q)$ 為平滑后的指示函數(shù) $%5Ctilde%20%5Cchi_M$ 的梯度場，最終公式（1）的解也是 $%5Ctilde%20%5Cchi_M$ 。
將公式（1）轉(zhuǎn)換為泊松方程
由于梯度場 $%5Cvec%20V$ 通常是不可積的，我們無法通過方程 $%5Cnabla%20%5Ctilde%5Cchi%3D%5Cvec%20V%20$ 來獲取指示函數(shù)，因此在公式（1）等式兩邊再加上一個散度，得到一個泊松方程：
$%5CDelta%20%5Ctilde%20%5Cchi%3D%5Cnabla%20%5Ccdot%20%7B%5Cvec%20V%7D$ ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(5)

三、算法的實現(xiàn)

1. 為求解泊松方程，需要定義一個函數(shù)空間，即定義一組基函數(shù) $F_o$ ，使得可以用 $F_o$ 的線性組合來有效表示 $%5Ctilde%20%5Cchi$ 和 $%5Cvec%20V$

直接選取平滑函數(shù) $%5Cwidetilde%20F_%7Bs.p%7D(q)$ 作為基函數(shù)。如公式（4）和圖2（a）所示， $%5Cvec%20V$ 已經(jīng)被表示為一組函數(shù)的線性組合了，但文章沒有采用，可能是因為這些基函數(shù)的中心僅位于表面周圍，無法很好表示其他區(qū)域處的指示函數(shù)。注：文章沒有提及這部分內(nèi)容，僅為我自己的想法。
使用3D網(wǎng)格劃分空間，為每個節(jié)點定義一個基函數(shù)。如圖2（b）所示，在每個節(jié)點的中心處定義一個平滑函數(shù)（高斯函數(shù)），并假設(shè)數(shù)據(jù)點均位于節(jié)點中心，那么 $%5Cvec%20V$ 仍可以用公式（4）表示（沒有數(shù)據(jù)點的節(jié)點的權(quán)重為0）。但隨著網(wǎng)格分辨率提升，計算量以三次方的速度遞增，并不實用，文章沒有采用該方法。
使用八叉樹劃分空間，保證所有數(shù)據(jù)點位于葉子節(jié)點，然后在每個節(jié)點中定義一個基函數(shù)，并假設(shè)數(shù)據(jù)點均位于節(jié)點中心，如圖2（c）所示。文章采用了八叉樹來劃分空間。對八叉樹O中的每個節(jié)點o，定義一個具有單位積分性質(zhì)的節(jié)點函數(shù) $F_o$ ，該節(jié)點函數(shù)以節(jié)點中心o.c為中心，以節(jié)點寬度o.w為標(biāo)準(zhǔn)差：

$F_o(q)%3DF%5Cleft(%20%5Cfrac%20%7Bq-o.c%7D%20%7Bo.w%7D%5Cright)%5Cfrac%20%7B1%7D%7B%7Bo.w%7D%5E3%7D$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? (6)
$F(q)%3D%5Cwidetilde%20F%5Cleft(%5Cfrac%20%7Bq%7D%7B2%5ED%7D%5Cright)$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(7)
D為八叉樹的最大深度， $F%5Cleft(%20%5Cfrac%20%7Bq-o.c%7D%20%7Bo.w%7D%5Cright)$ 表示函數(shù)以o.c為中心，o.w為寬度（高斯函數(shù)的標(biāo)準(zhǔn)差）， $%5Ctilde%20%5Cchi$ 表示將函數(shù)寬度轉(zhuǎn)化為 $o.w*2%5ED$ ，即八叉樹的總寬度。由于所有數(shù)據(jù)點在葉子節(jié)點，故在表示 $%5Cvec%20V$ 時，公式（4）可以修改為：
$%5Cvec%20V(q)%3D%7Bo.w%7D%5E3*%5Csum_%7Bs%5Cin%20S%7DF_%7Bo%7D(q)s.%5Cvec%20N$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(8)
注意：這部分分析不保證正確，看看就好。
使用周圍8節(jié)點的線性插值表示數(shù)據(jù)點。作者嫌上面直接用節(jié)點中心表示數(shù)據(jù)點誤差大，故使用周圍8節(jié)點中心的三線性插值來替代公式（4）中的數(shù)據(jù)點，如圖2（d）所示：

$%5Cvec%20V(q)%3D%5Csum_%7Bs%5Cin%20S%7D%5Csum_%7Bo%5Cin%20Ngbr_D(s)%7D%5Calpha_%7Bo%2Cs%7DF_o(q)s.%5Cvec%20N$ ? ? ? ? ? ? ? ? ? ? ? ? ? ?(9)
$Ngbr_D(s)$ 為第D層最接近數(shù)據(jù)點的8個節(jié)點， $%5Calpha_%7Bo%2Cs%7D$ 為三線性插值權(quán)重。公式（9）中的 $%5Cvec%20V(q)$ 與真實表面法向量相差常數(shù)倍。

圖2? （a）直接使用平滑函數(shù)作為基函數(shù)來表達(dá)表面上任意點x，見公式（4）；（b）使用3D網(wǎng)格劃分空間，在節(jié)點中心處定義基函數(shù)，并使用節(jié)點中心替代數(shù)據(jù)點；（c）使用八叉樹劃分空間，在葉子節(jié)點處定義基函數(shù)，并使用節(jié)點中心替代數(shù)據(jù)點；（d）使用八叉樹劃分空間，在葉子節(jié)點處定義基函數(shù)，并使用周圍4節(jié)點（3D是8節(jié)點）的插值來替代數(shù)據(jù)點

2. 在基函數(shù) $F_o$ 定義的函數(shù)空間 $O$ 中描述泊松方程

求得 $%5Cnabla%20%5Ccdot%20%5Cvec%20V$ 在函數(shù)空間 $O$ 中的坐標(biāo)。已知 $%5Cvec%20V$ 可以用 $F_o$ 的線性組合表示，即已知 $%5Cvec%20V$ 在O中的坐標(biāo)，應(yīng)該是可以進(jìn)一步計算出 $%5Cnabla%20%5Ccdot%20%5Cvec%20V$ 在O中的坐標(biāo)v
在函數(shù)空間O中表示 $%5CDelta%20%5Ctilde%20%5Cchi$ 。令 $%5Ctilde%20%5Cchi%3D%5Csum_ox_oF_o$ ，我們是要求解所有 $x_o$ 組成的向量x。為此將 $%5CDelta%20%5Ctilde%20%5Cchi$ 定義為矩陣相乘的形式， $%5CDelta%20%5Ctilde%20%5Cchi%3DLx$ ，對于L中的(o,o')項，滿足下式：
$%5Cpartial%20%5Ctilde%20M%5Cequiv%20%5C%7Bq%5Cin%20R%5E3%7C%5Ctilde%20%5Cchi%20(q)%3D%5Cgamma%20%5C%20with%5C%20%20%5Cgamma%3D%5Cfrac%20%7B1%7D%7B%7CS%7C%7D%5Csum_%7Bs%5Cin%20S%7D%5Ctilde%20%5Cchi%20(s.p)$ ? ? ? ? ? ? (10)
$%3C%5Cfrac%7B%7B%5Cpartial%7D%5E2F_o%7D%7B%5Cpartial%20x%5E2%7D%2CF_%7Bo'%7D%3E$ 表示經(jīng)過拉普拉斯算子運(yùn)算后在基函數(shù) $F_%7Bo'%7D$ 上的投影。

3. 使用Gauss-Seidel方法求解方程組Lx=v

4. 計算表面的數(shù)據(jù)值以提取表面。計算 $%5Ctilde%20%5Cchi$ 在所有數(shù)據(jù)點處的值，選取平均值作為表面的數(shù)據(jù)值。其中 $%5Ctilde%20%5Cchi$ 乘以一個常數(shù)，并不影響最終的表面結(jié)果。

$%5Cpartial%20%5Ctilde%20M%5Cequiv%20%5C%7Bq%5Cin%20R%5E3%7C%5Ctilde%20%5Cchi%20(q)%3D%5Cgamma%20%5C%20with%5C%20%20%5Cgamma%3D%5Cfrac%20%7B1%7D%7B%7CS%7C%7D%5Csum_%7Bs%5Cin%20S%7D%5Ctilde%20%5Cchi%20(s.p)$ ? ? ? ? ? ? (11)

5. 考慮不均勻采樣點的情況

計算每個數(shù)據(jù)點處的點云密度。選取一個深度 $%5Chat%20D$ ，該深度小于八叉樹最大深度 $D$ ，計算點云密度：
$W_%7B%5Chat%20D%7D(q)%5Cequiv%20%5Csum_%7Bs%5Cin%20S%7D%5Csum_%7Bo%5Cin%20Ngbr_%7B%5Chat%20D%7D(s)%7D%7B%5Calpha%7D_%7Bo%2Cs%7DF_o(q)$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? (12)
將權(quán)重加入到 $%5Cvec%20V$ 計算公式中
$%5Cvec%20V(q)%3D%5Csum_%7Bs%5Cin%20S%7D%20%5Cfrac%7B1%7D%7BW_%7B%5Chat%20D%7D(s.p)%7D%5Csum_%7Bo%5Cin%20Ngbr_D(s)%7D%7B%5Calpha%7D_%7Bo%2Cs%7DF_o(q)s.%5Cvec%20N$ ? ? ? ? ? ? ? ?(13)
將權(quán)重加入表面數(shù)據(jù)值計算公式中
$%5Cpartial%20%5Ctilde%20M%5Cequiv%20%5C%7Bq%5Cin%20R%5E3%7C%5Ctilde%20%5Cchi(q)%3D%5Cgamma%5C%20%20with%20%5C%20%5Cgamma%3D%5Cfrac%7B%5Csum%20%5Cfrac%7B1%7D%7BW_%7B%5Chat%20D%7D(s.p)%7D%5Ctilde%20%5Cchi%20(s.p)%7D%7B%5Csum%20%5Cfrac%7B1%7D%7BW_%7B%5Chat%20D%7D(s.p)%7D%7D$ ? ? ? ?(14)

四、屏蔽泊松表面重建算法

重定義指示函數(shù) $%5Cchi_M$ 。定義物體內(nèi)部函數(shù)為1/2，物體外部函數(shù)值為-1/2，令物體表面函數(shù)值為0.
定義一個新的能量函數(shù)值，以確保零值表面是貼近點云的。

原能量函數(shù)如公式（1）所示，新能量函數(shù)如公式（15）所示：

$E(%5Cchi)%3D%5Cint%20%7B%5CVert%20%5Cvec%20V-%5Cnabla%5Cchi(p)%5CVert%7D%5E2dp$ ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?(1)
$E(%5Cchi)%3D%5Cint%7B%5C%7C%5Cvec%20V%20(p)-%5Cnabla%20%5Cchi(p)%20%5C%7C%7D%5E2dp%2B%5Cfrac%20%7B%5Calpha%20%5Ccdot%20Area(P)%7D%7B%5Csum_%7Bp%5Cin%20P%7D%5Comega(p)%7D%5Csum_%7Bp%5Cin%20P%7D%5Comega(p)%7B%5Cchi%7D%5E2(p)$ ? ?(15)
其中 $%5Calpha$ 是平衡梯度擬合和數(shù)值擬合的權(quán)值； $Area(P)$ 為表面的面積，可由之前算到的局部點云密度估計； $%5Comega(p)$ 是每個數(shù)據(jù)點的權(quán)重，在PCL中用數(shù)據(jù)點法向量的模表示，默認(rèn)是1.
使用B樣條函數(shù)替代 $F_o$ 來描述泊松方程
2013年的文章直接使用了B樣條函數(shù)作為基函數(shù)來描述泊松方程。在PCL的實現(xiàn)中，先使用 $F_o$ 計算出每個網(wǎng)格處的 $%5Cvec%20V$ ，然后使用B樣條函數(shù)計算 $%5Cnabla%20%5Cvec%20V$ 、求解線性方程組和提取表面。
屏蔽泊松表面重建算法存在過擬合的情況
當(dāng)點云中存在噪聲，屏蔽泊松算法會過于擬合噪聲點，形成粗糙表面，如圖3所示，而且 $%5Calpha$ 越大，受噪聲影響越大?？赏ㄟ^在包含多個數(shù)據(jù)點的網(wǎng)格中（降低八叉樹深度）構(gòu)建泊松方程來降低噪聲影響。

五、參考文獻(xiàn)

M. Kazhdan, M. Bolitho, H. Hoppe. Poisson Surface Reconstruction. Eurographics Symposium on Geometry Processing. 2006: 61-70
J. C. Carr, R. K. Beatson, J. B. Cherrie, T. J. Mitchell, W. R. Fright, B. C. McCallum, T. R. Evans. Reconstruction and Representation of 3D Objects with Radial Basis Functions. Proceedings of the 28th Annual Conference on Computer? Graphics and Interactive Techniques. 2001: 67-76
M. Kazhdan,?H. Hoppe.?Screened Poisson Surface Reconstruction. ACM Trans. Graph. 32(3): 29

標(biāo)簽：