從點雲到網格（三）Poisson重建

Possion重建是Kazhdan等2006年提出的網格重建方法[1]。Possion重建的輸入是點雲及其法向量，輸出是三維網格。Poisson有公開的源代碼[2]。PCL中也有Poisson的實現。

核心思想

Possion重建是一個非常直觀的方法。它的核心思想是點雲代表了物體表面的位置，其法向量代表了內外的方向。通過隱式地擬合一個由物體派生的指示函數，可以給出一個平滑的物體表面的估計。

給定一個區域\(M\)及其邊界\(\partial M\)，指示函數\(\chi_M\)定義為

這樣，把重構\(S = \partial M\)的問題轉換為重構\(\chi_M\)的問題。作者給出了將點雲及其法向量和\(\chi_M\)聯系起來的公式。作者論文中的圖1非常形象地描述了這二者的聯系。

基本思路

對於任意點\(p\in \partial M\)，定義\(\vec{N}_{\partial M}(p)\)為向內的表面法向量，\(\tilde{F}(q)\)為一個平滑濾波器，\(\tilde{F}_p(q)=\tilde{F}(q-p)\)為\(\tilde{F}\)沿\(p\)方向的平移。因為\(\chi_M\)一般意義上是不好求導的，這里用\(\chi_M\ast\tilde{F}\)的導數來近似。

從法向量到梯度空間

梯度空間的近似

由於\(\vec{N}_{\partial M}\)的分布是未知的，需要通過觀測值\(P=\left\lbrace \left( p_i,n_i\right) \right\rbrace\)來近似。考慮離散點集\(\Omega\)，\(\partial M\)被分割為互不相交的區域\(\wp_s,s\in\Omega\)。\((1)\)可以轉化為積分求和，並將每個小積分近似為常函數，用代表點\(s.p\)對應的函數值和\(\wp_s\)的面積的乘積代替。

這里希望平滑濾波器\(\tilde{F}\)能夠盡量地窄，不過分平滑數據，同時盡量地寬，使得積分近似能夠更准確。高斯濾波器是一種常見的選擇。

求解Possion問題

向量空間\(\vec{V}\)和指示函數\(\tilde{\chi}\)滿足

\[\nabla\tilde{\chi}=\vec{V} \tag{3} \]

然而，\(\vec{V}\)通常意義上是沒法積分的（為什么？）。為了得到\((3)\)式的最小二乘解，將\((3)\)式兩邊求導，就得到了拉普拉斯方程

\[\Delta\tilde{\chi}=\nabla\cdot\vec{V} \tag{4} \]

拉普拉斯方程在數學上有很詳細的研究。

實現細節

空間划分

為了解一個偏微分方程問題，首先要將其離散化。最簡單的方法是將空間划分為均勻網格。這種划分非常占內存空間，而且只有邊界附近的值才是我們關心的，大量的空間被浪費了。作者采用了一種自適應的網格結構octree來划分空間，並且octree上定義了一個函數空間\(F_o\)。下圖給出了octree在三維空間對一個物體的划分。物體邊緣的網格密度遠大於遠離物體的網格密度。

http://http.developer.nvidia.com/GPUGems2/elementLinks/37_octree_03.jpg

空間上的基函數選擇

如何選擇函數空間\(F_o\)實際上挺有學問的。因為\(F_o\)一旦給定，定義在這個octree上的向量空間\(\vec{V}\)和指示函數\(\chi\)都會通過\(F_o\)的線性組合去近似。這樣，求解\(\chi\)就轉化為求解\(F_o\)上的參數組合，進而轉化為求解一個線性方程組。

給定octree的深度\(D\)，作者根據選擇了下面的基函數\(F\)。

\(*n\)代表\(n\)次卷積。當\(n\)趨向於無窮時，\(F\)趨向於高斯函數，它的定義域也越來越大。當\(n=3\)時，定義域為\([-1.5,1.5]^3\)。

octree上某個節點\(o\)對應的函數\(F_o\)定義為

\[F_o(q)\equiv F\left(\frac{q-o.c}{o.w}\right)\frac{1}{(o.w)^3} \]

其中\(o.c\)是的\(o\)中心，\(o.w\)是\(o\)的寬度。假設根節點（第0層）的寬度為\(W\)，那么第\(d\)層節點的寬度為\(\frac{W}{2^{d}}\)。這個函數空間和小波空間很像。

Possion求解

Possion求解方法是L2投影（L2 projection）[3]。定義octree的節點集合為\(O\)。向量空間\(\vec{V}\)可以近似為

\[\vec{V}(q)\equiv \Sigma_{s\in \Omega}\Sigma_{o\in Ng(s)}\alpha_{o,s}F_o(q)s.\vec{n} \]

其中\(Ng(s)\)是\(s\)的八個最近鄰的葉節點，\(\alpha_{o,s}\)是三線性插值的權重。

雖然\(\vec{V}\)和\(\tilde{\chi}\)都可以在函數空間上表示出來，\(\Delta\tilde{\chi}\)和\(\nabla\cdot\vec{V}\)卻未必有定義。因此將\((4)\)近似為最小化其在\(F_o\)上的投影

\[\Sigma_{o} \lVert\langle \Delta \tilde{\chi}-\nabla\cdot\vec{V},F_o\rangle\lVert^2=\Sigma_{o} \lVert\langle \Delta \tilde{\chi},F_o\rangle-\langle\nabla\cdot\vec{V},F_o\rangle\lVert^2 \]

令\(\tilde{\chi}=\Sigma_ox_oF_o\)，那么求解\(\tilde{\chi}\)即是求解\(x_o\)。

\[\langle \Delta \tilde{\chi},F_{o'}\rangle=\Sigma_{o}x_{o}\langle \Delta F_o,F_{o'}\rangle \]

\[\Sigma_{o} \lVert\langle \Delta \tilde{\chi},F_o\rangle-\langle\nabla\cdot\vec{V},F_o\rangle\lVert^2 =\Sigma_{o'}\lVert\Sigma_ox_o\langle \Delta F_o,F_{o'}\rangle-\langle\nabla\cdot\vec{V},F_{o'}\rangle\lVert^2 \]

上式右邊對\(x=\lbrace x_o\rbrace\)求偏導，轉化為

\[\min \limits_{x}\left\lVert Lx-v \right\lVert ^2 \]

其中，設octree的節點數為\(N\)，\(N\times N\)矩陣\(L\)在\((o,o')\)位置上的值為

\[L_{o,o'}\equiv \langle\frac{\partial^2 F_o}{\partial x^2},F_{o'}\rangle+\langle\frac{\partial^2 F_o}{\partial y^2},F_{o'}\rangle+\langle\frac{\partial^2 F_o}{\partial z^2},F_{o'}\rangle \]

表面提取

用Marching Cubes類似的方法。注意iso的值取自\(S\)個划分的平均。

作者還討論了非均勻采樣下的算法，在此就不贅述。

Poisson分析

簡單列幾點

• Poisson在邊緣處的銳度比VRIP要好。這是因為VRIP在大的邊緣處TSDF的累加會有平滑效應，而Poisson依據的是法向量，不會引入額外的平滑。
• VRIP是局部方法，每次只更新當前深度圖對應的TSDF。Poisson是全局方法。
• 從個人使用經驗上看，Poisson對於噪聲更加魯棒一些。點雲法向量估計的精度不能太差。
• 如果重建出奇怪的形狀（分層、分塊），請查看原始點雲是否平滑，是否有噪聲，調整生成網格的分辨率以適應點雲。

小結

Poisson是個好方法。

參考文獻

[1]. Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006.
[2]. http://www.cs.jhu.edu/~misha/Code/PoissonRecon/Version8.0/
[3]. http://www.featflow.de/en/software/featflow2/tutorial/tutorial_l2proj.html