泊松表面重建算法原理

2022-10-16 00:18 作者:生医小王子 0人读过 | 我要投稿

一、概述

泊松表面重建算法（Poisson Surface Reconstruction）[1]是M. Kazhdan等人于2006年提出的点云表面重建算法，与基于径向基函数（Radial Basis Function，RBF）[2]的移动立方体算法（Marching Cube，MC）相似，都是先定义一个隐函数来表示表面（Surface），根据点云数据求解隐函数的具体表达式，然后使用MC方法将表面三角化。但有2点不同，首先使用的隐函数不同，其次使用了八叉树来替代3D网格。

随后，M. Kazhdan等人在2013年提出了屏蔽泊松表面重建算法（Screened Poisson Surface Reconstruction）[3]，他们发现原方法的表面过于平滑了，所以改进了算法，以保证重建表面更贴近于点云。PCL中实现的泊松重建算法即为该改进算法。

泊松算法的输入为带法向量和空间坐标的点云数据，输出为用三角网格表示的表面，其中法向量应该保证是准确的（全都指向物体内部或者全指向外部）。在下面的第2和第3章将分别介绍算法的推导过程和实现方法。第4章将介绍屏蔽泊松表面重建算法的不同。

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

（专栏有100图片限制，然后公式也算图片，也是服了。。）

二、算法推导过程

算法大致思路：找到一个隐函数表示表面，然后使用移动立方体等表面提取算法来三角化表面

对于物体M，文章定义了一个3D连续分布的指示函数 $%5Cchi%20_%7BM%7D$ ，物体内部的函数值为1，物体外部函数值为0，如图1所示。若能求得该指示函数，通过MC方法可以很方便地计算出物体表面。
点云法向量是指示函数梯度的采样

对指示函数 $%5Cchi%20_%7BM%7D$ 求导计算 $%5Cnabla%7B%5Cchi%20_M%7D$ ，由于物体内部和外部区域都是常数，这些区域的梯度为0，仅在表面处有数值，所以 $%5Cnabla%7B%5Cchi%20_M%7D$ 为指向物体内部的表面法向量 $%5Cvec%20V$ 。当规定外部输入的点云法向量都是指向物体内部时，这些点云法向量是 $%5Cwidetilde%20%5Cchi_M$ 的离散采样，应该是有办法根据这些离散采样值估计出 $%5Cvec%20V$ 的分布的。指示函数 $%5Cchi%20_%7BM%7D$ 应该保证下面能量函数是最小的：

$E(%5Cchi)%3D%5Cint%20%7B%5CVert%20%5Cvec%20V-%5Cnabla%5Cchi(p)%5CVert%7D%5E2dp$ (1)
确定表面法向量 $%5Cvec%20V$ 的表达式
由于 $%5Cchi%20_%7BM%7D$ 是分段常量函数， $%5Cvec%20V$ 在表面处的数值是无限大的，不好用点云的法向量表示，故退而求其次，使用平滑滤波器对 $%5Cchi%20_%7BM%7D$ 做卷积，得到平滑的指示函数 $%5Ctilde%20%5Cchi_M$ 。其中 $%5Ctilde%20%5Cchi_M$ 的梯度为 $%5Cchi%20_%7BM%7D$ 的表面法向量平滑后的结果，如下式所示：
$%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)$ 为表面的内法向量， $p%5Cin%20%5Cpartial%20M$ ； $%5Cwidetilde%20F(q)$ 为平滑滤波器，且满足 $%5Cwidetilde%20F(q)%3D%5Cwidetilde%20F(q-p)$ ，即模糊权值与两点之间的距离有关。
假设点云S将表面 $%5Cpartial%20M$ 分成了不同的小块 $%5CPsi_s$ ，可以使用数据点的法向量近似表示 $%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$ 分别为数据点的位置和法向量。 $%5Cwidetilde%20F_%7Bs.p%7D(q)$ 可以选用高斯函数。当点云均匀分布时， $%7C%5CPsi_s%7C$ 为常数，该常数不影响最终结果（后面会说明），故上式（3）可以简化为：
$%5Cvec%20V(q)%3D%5Csum_%7Bs%5Cin%20S%7D%5Cwidetilde%20F_%7Bs.p%7D(q)s.%5Cvec%20N$ (4)
注意： $%5Cvec%20V(q)$ 为平滑后的指示函数 $%5Ctilde%20%5Cchi_M$ 的梯度场，最终公式（1）的解也是 $%5Ctilde%20%5Cchi_M$ 。
将公式（1）转换为泊松方程
由于梯度场 $%5Cvec%20V$ 通常是不可积的，我们无法通过方程 $%5Cnabla%20%5Ctilde%5Cchi%3D%5Cvec%20V%20$ 来获取指示函数，因此在公式（1）等式两边再加上一个散度，得到一个泊松方程：
$%5CDelta%20%5Ctilde%20%5Cchi%3D%5Cnabla%20%5Ccdot%20%7B%5Cvec%20V%7D$ (5)

三、算法的实现

1. 为求解泊松方程，需要定义一个函数空间，即定义一组基函数 $F_o$ ，使得可以用 $F_o$ 的线性组合来有效表示 $%5Ctilde%20%5Cchi$ 和 $%5Cvec%20V$

直接选取平滑函数 $%5Cwidetilde%20F_%7Bs.p%7D(q)$ 作为基函数。如公式（4）和图2（a）所示， $%5Cvec%20V$ 已经被表示为一组函数的线性组合了，但文章没有采用，可能是因为这些基函数的中心仅位于表面周围，无法很好表示其他区域处的指示函数。注：文章没有提及这部分内容，仅为我自己的想法。
使用3D网格划分空间，为每个节点定义一个基函数。如图2（b）所示，在每个节点的中心处定义一个平滑函数（高斯函数），并假设数据点均位于节点中心，那么 $%5Cvec%20V$ 仍可以用公式（4）表示（没有数据点的节点的权重为0）。但随着网格分辨率提升，计算量以三次方的速度递增，并不实用，文章没有采用该方法。
使用八叉树划分空间，保证所有数据点位于叶子节点，然后在每个节点中定义一个基函数，并假设数据点均位于节点中心，如图2（c）所示。文章采用了八叉树来划分空间。对八叉树O中的每个节点o，定义一个具有单位积分性质的节点函数 $F_o$ ，该节点函数以节点中心o.c为中心，以节点宽度o.w为标准差：

$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)$ 表示函数以o.c为中心，o.w为宽度（高斯函数的标准差）， $%5Ctilde%20%5Cchi$ 表示将函数宽度转化为 $o.w*2%5ED$ ，即八叉树的总宽度。由于所有数据点在叶子节点，故在表示 $%5Cvec%20V$ 时，公式（4）可以修改为：
$%5Cvec%20V(q)%3D%7Bo.w%7D%5E3*%5Csum_%7Bs%5Cin%20S%7DF_%7Bo%7D(q)s.%5Cvec%20N$ (8)
注意：这部分分析不保证正确，看看就好。
使用周围8节点的线性插值表示数据点。作者嫌上面直接用节点中心表示数据点误差大，故使用周围8节点中心的三线性插值来替代公式（4）中的数据点，如图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层最接近数据点的8个节点， $%5Calpha_%7Bo%2Cs%7D$ 为三线性插值权重。公式（9）中的 $%5Cvec%20V(q)$ 与真实表面法向量相差常数倍。

图2 （a）直接使用平滑函数作为基函数来表达表面上任意点x，见公式（4）；（b）使用3D网格划分空间，在节点中心处定义基函数，并使用节点中心替代数据点；（c）使用八叉树划分空间，在叶子节点处定义基函数，并使用节点中心替代数据点；（d）使用八叉树划分空间，在叶子节点处定义基函数，并使用周围4节点（3D是8节点）的插值来替代数据点

2. 在基函数 $F_o$ 定义的函数空间 $O$ 中描述泊松方程

求得 $%5Cnabla%20%5Ccdot%20%5Cvec%20V$ 在函数空间 $O$ 中的坐标。已知 $%5Cvec%20V$ 可以用 $F_o$ 的线性组合表示，即已知 $%5Cvec%20V$ 在O中的坐标，应该是可以进一步计算出 $%5Cnabla%20%5Ccdot%20%5Cvec%20V$ 在O中的坐标v
在函数空间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$ 表示经过拉普拉斯算子运算后在基函数 $F_%7Bo'%7D$ 上的投影。

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

4. 计算表面的数据值以提取表面。计算 $%5Ctilde%20%5Cchi$ 在所有数据点处的值，选取平均值作为表面的数据值。其中 $%5Ctilde%20%5Cchi$ 乘以一个常数，并不影响最终的表面结果。

$%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. 考虑不均匀采样点的情况

计算每个数据点处的点云密度。选取一个深度 $%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)
将权重加入到 $%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)
将权重加入表面数据值计算公式中
$%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)

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

重定义指示函数 $%5Cchi_M$ 。定义物体内部函数为1/2，物体外部函数值为-1/2，令物体表面函数值为0.
定义一个新的能量函数值，以确保零值表面是贴近点云的。

原能量函数如公式（1）所示，新能量函数如公式（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$ 是平衡梯度拟合和数值拟合的权值； $Area(P)$ 为表面的面积，可由之前算到的局部点云密度估计； $%5Comega(p)$ 是每个数据点的权重，在PCL中用数据点法向量的模表示，默认是1.
使用B样条函数替代 $F_o$ 来描述泊松方程
2013年的文章直接使用了B样条函数作为基函数来描述泊松方程。在PCL的实现中，先使用 $F_o$ 计算出每个网格处的 $%5Cvec%20V$ ，然后使用B样条函数计算 $%5Cnabla%20%5Cvec%20V$ 、求解线性方程组和提取表面。
屏蔽泊松表面重建算法存在过拟合的情况
当点云中存在噪声，屏蔽泊松算法会过于拟合噪声点，形成粗糙表面，如图3所示，而且 $%5Calpha$ 越大，受噪声影响越大。可通过在包含多个数据点的网格中（降低八叉树深度）构建泊松方程来降低噪声影响。

五、参考文献

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