泊松重建（Poisson Surface Reconstruction）

简介

泊松重建是Michael Kazhdan等在2006年提出的网格重建方法，其文章题目为“Poisson Surface Reconstruction”。
Poisson-Rec的输入是带有法向量属性的点云数据（也可以有RGB信息），输出是三角网格模型，下面来直观的感受一下泊松重建的输入和输出：

表面重建流程：

1、构建八叉树：采用的是自适应的空间网格划分的方法（根据点云的密度调整网格的深度），根据采样点集的位置定义八叉树，然后细分八叉树使每个采样点都落在深度为D的叶节点；
2、设置函数空间：对八叉树的每个节点设置空间函数 F，所有节点函数 F 的线性和可以表示向量场 V，基函数F采用了盒滤波的n维卷积；
3、创建向量场：均匀采样的情况下，假设划分的块是常量，通过向量场 V 逼近指示函数的梯度。采用三次条样插值（三线插值）；
4、求解泊松方程：方程的解采用拉普拉斯矩阵迭代求出；
5、提取等值面：为得到重构表面，需要选择阈值获得等值面；先估计采样点的位置，然后用其平均值进行等值面提取，然后用 Marching Cubes（移动立方体）算法得到等值面。

算法原理

Possion-Rec是一个非常直观的方法。它的核心思想是点云代表了物体表面的位置，其法向量代表了内外的方向。通过隐式地拟合一个由物体派生的指示函数，可以给出一个平滑的物体表面的估计。
给定一个区域 M 及其边界 ∂M ，指示函数 χM 定义为：

这样，把重构 S = ∂M 的问题转换为重构 χM 的问题。
作者给出了将点云及其法向量和 χM 联系起来的公式，下图非常形象地描述了这二者的联系。

这个指示函数是个分段函数，定义模型内部的值大于0，外部的值小于0，而为0的部分即为等值面，提取出来就是几何目标模型的表面。

基本思路
由梯度关系得到采样点和指示函数的积分关系，根据积分关系利用划分块的方法获得点集的向量场，计算指示函数梯度场的逼近，构成泊松方程。根据泊松方程使用矩阵迭代求出近似解，采用移动立方体算法提取等值面，对所测数据点集重构出被测物体的模型，泊松方程在边界处的误差为零，因此得到的模型不存在假的表面框
直接计算梯度场会引起向量场在表面边缘的无穷大值。因此首先用平滑滤波卷积指示函数，然后求平滑函数的梯度场。
高斯散度理论：平滑指示函数的梯度等于平滑表面法向场得到的向量场
由于曲面未知，无法直接计算表面积分，把采样点集划分为小的区域块，通过对所有块的积分求和近似计算。知道向量场V后，可求指示函数。但向量场V不可积，使用最小平方逼近理论求解，应用散度算子得到泊松方程。
泊松表面重建一次性把所有的点都考虑在内，因此对噪声点有很好的弹性。泊松方程允许的层次结构支持局部的基函数，因此对稀疏线性系统的情况有很好的支持。在此基础上描述了多尺度的空间自适应算法，其时间和空间复杂度同重建模型的大小成正比。

从法向量到梯度空间
对于任意点 $p\in \partial M$ ，定义$N⃗\partial M(p)$为向内的表面法向量，$F(q)$为一个平滑滤波器，$Fp(q)=F(q-p)$为$F$沿p方向的平移。因为 $\chi M$ 一般意义上是不好求导的，这里用 $\chi M\ast F$ 的导数来近似
(1)

梯度空间的近似
由于$N⃗\partial M(p)$的分布是未知的，需要通过观测$P=(pi,ni)$来近似。考虑离散点集 $\Omega$，$\partial M$被分割为互不相交的区域 $\wp s,s\in \Omega \wp s,s\in \Omega$。（1）可以转化为积分求和，并将每个小积分近似为常函数（可对照上图中的第二副图），用代表点$s.p$对应的函数值和$\wp s$的面积的积分代替：
(2)

这里希望平滑滤波器$F$能够尽量地窄，不过分平滑数据，同时尽量地宽，使得积分近似能够更准确。高斯滤波器是一种常见的选择。

求解泊松问题
向量空间 $V⃗$和指示函数 $\chi$ 满足
(3)

然而，$V⃗$向量场通常意义上是没法积分的。为了得到(3)式的最小二乘解，将(3)式两边求导，就得到了拉普拉斯方程：
(4)

作者有一个预设，且有一系列数学的证明，证明了指示函数的梯度等于结合表面法线场计算得到的向量场。这也是高斯散度理论。(理论如此，实际是求近似）
倒三角是梯度(希望到这里你还撑的住。)
但是，指示函数不知道，梯度不知道。只能计算出向量场。所以应用散度算子这一媒介。转化成了一个泊松方程。

（倒三角和正三角什么鬼！）其中，

是我们想要求解出的函数，这个正三角就是，梯度的散度（也是拉普拉斯算子）
也就是 $\chi$ 梯度的散度等于向量场的散度，即：

解这个泊松方程就能找到指示函数。方程的解采用拉普拉斯矩阵迭代求出

梯度与散度补充：

3. 散度

具体实现

空间划分
作者采用了一种自适应的网格结构octree来划分空间，并且octree上定义了一个函数空间$Fo$。下图给出了octree在三维空间对一个物体的划分,物体边缘的网格密度远大于远离物体的网格密度。

空间上的基函数选择
如何选择函数空间$Fo$实际上挺有学问的。因为$Fo$一旦给定，定义在这个octree上的向量空间$V⃗$和指示函数$\chi$都会通过$Fo$的线性组合去近似。这样，求解$\chi$就转化为求解$Fo$上的参数组合，进而转化为求解一个线性方程组。
给定octree的深度$D$，作者根据选择了下面的基函数$F$：

$\ast n$代表n次卷积。当$n$趋向于无穷时，$F$趋向于高斯函数，它的定义域也越来越大。当$n=3$时，定义域为$[-1.5,1.5{]}^3$。
octree上某个节点o对应的函数$Fo$定义为：

其中$o.c$是的$o$中心，$o.w$是$o$的宽度。假设根节点（第0层）的宽度为$W$，那么第$d$层节点的宽度为$W/2^d$。这个函数空间和小波空间很像。

Poisson求解
（坚持一下，就快讲完了！）
Poisson求解的方法是$L2$投影。定义octree的节点集合为$O$。向量空间$V⃗$可以近似为：

其中，$NgbrD(s)$是$s$的八个最近邻的叶节点，$\alpha o,s$是三线性插值的权重
虽然$V⃗$和$\chi$都可以在函数空间上表示出来，$\Delta \chi$和$\nabla \cdot V⃗$ 却未必有定义。因此将(4)近似为最小化其在$Fo$上的投影：

令
那么求解$\chi$即是求解$xo$


对上式右边$x=x0$求偏导得，

其中，设octree的节点数为$N$，$N\times N$矩阵$L$在$(o,o\prime )$位置上的值为：

表面提取
用Marching Cubes类似的方法。注意 iso（等值面） 的值取自S个划分的平均，作者还讨论了非均匀采样下的算法。