0x00 To begin with

这篇文章将会包含以下内容：

  • 泊松重建的一些前置数学知识
  • 泊松重建算法原理及推导过程

For reference：

0x01 隐函数表达 Implicit Expression

表达一个三维曲面，例如球心位于原点、半径为 $r$ 的球面，可以使用球的标准方程：

$$ F(x,y,z):x^{2}+y^{2}+z^{2}-r^{2}=0 $$
该方程（形如 $F(x,y,z)=0$）就定义了球面的隐函数表达，它将所有变量都放在等式一侧，隐含了因变量关于自变量的函数关系，在描述球面的情况下，$z$ 关于 $x,y$ 的关系可以显化为：
$$ z=\sqrt{ r^{2}-x^{2}-y^{2} } $$
这就是隐函数的一种显式表达，它描述球面在 $xOy$ 平面以上的部分。如果需要描述整个球面，还需要 $xOy$ 平面以下的部分：
$$ z=-\sqrt{ r^{2}-x^{2}-y^{2} } $$
也就是说显化的隐函数可能无法覆盖所有满足隐函数方程的解集，因为显函数是一一映射，而隐函数可能是一对多或者多对多的映射。

在计算机图形学中，SDF 就是一种常用的可以描述表面的隐函数表达：

$$ \phi(\mathbf{x})=\pm \text{dist}(\mathbf{x},\partial \Omega) $$
SDF 以物体表面为界，空间中任意一点 $\mathbf{x}$ 的 SDF 函数值为与物体表面的带符号距离（物体内部为负、外部为正）。自然的，物体表面可以用 SDF 值为零的水平集描述：
$$ \{\mathbf{x}\,|\,\phi(\mathbf{x})=0\} $$
![使用 SDF 的零水平集描述表面，Jeong et.al](https://raw.gitmirror.com/congyuxiaoyoudao/Picgo-ImageBed/main/TA-interludes/Interlude%205.PoissonRecon/202508201058142.png#center)

0x02 泊松方程 Poisson Equation

泊松表面重建其实不需要对泊松方程有什么很深刻的数学物理了解，笔者也非相关专业，所以就仅凭自己的想法写一点更为直观的解释

泊松方程由法国数学家泊松在对势的拉普拉斯方程作二阶修正时提出，今天更为人熟知的一般形式如下：

$$ \nabla^{2}\phi=-f,\;\phi:\Omega\to \mathbb{R} $$
其中：
  • $\phi$ 为定义在空间 $\Omega$ 上描述标量场（如温度、电势）的未知函数
  • $\nabla^{2}$ 为拉普拉斯算子（也可以表示为 $\Delta$），表示对函数在空间 $\Omega$ 内的二阶导数求和，该算子相当于梯度的散度 $\nabla^{2}=\nabla \cdot \nabla$
  • $f$ 为一个已知函数，描述场量的分布（如热源强度、电荷密度），称为源项

总体来说泊松方程描述的是有源情况下源项如何影响场量的分布（热源影响温度、电荷密度影响电势），而求解泊松方程，就是在给定源项 $f$ 的条件下找到满足等式的 $\phi$。

为了获得期望的解，需要给定一个边界条件 $\phi|_{\partial\Omega}$，其描述了场量函数 $\phi$ 在空间 $\Omega$ 内由于某种条件约束而产生的边界 $\partial \Omega$，通常有两类边界条件：

  • Dirichlet （狄利克雷）边界：$\phi|_{\partial \Omega}=g$，边界的场量已知，比如边界温度已知，边界满足：$\forall p \in \partial \Omega,\;\phi|_{\partial \Omega}(p)=g(p)$
  • Neumann （纽曼）边界：$\frac{\partial\phi}{\partial n}|_{\partial \Omega}=h$，边界上沿法向的净流出量已知，比如边界热通量已知，边界满足：$\forall p \in \partial \Omega,\;\frac{\partial\phi}{\partial n}|_{\partial \Omega}(p)=h(p)$

从更加数学物理的角度看，泊松方程可以对应一个能量泛函：

$$ E[\phi]=\frac{1}{2} \int_{\Omega}|\nabla\phi|^{2}d\Omega-\int_{\Omega}f\phi d\Omega $$
求解泊松方程相当于让泛函 $E[\phi]$ 对任意小变动 $\delta{\phi}$ 的一阶变分为零，计算一阶变分可以得到：
$$ \delta E=\int_{\Omega}(-\Delta\phi-f)d\Omega $$
所以 $-\Delta \phi-f=0 \Leftrightarrow \nabla^{2}\phi=-f$，这就将泊松方程的解和能量泛函存在极小值联系起来，求解泊松方程就是求解有源情况下能量的局部极小稳定态。

0x03 泊松表面重建 Poisson Surface Reconstruction

泊松表面重建（Poisson Surface Reconstruction）是一种隐式表面重建方法，对于重建任务，给定待重建物体 $M$，目标是获取其边界 $\partial M$。输入有向点云（有法向的点集）$S$，可以将 $\forall s \in S$ 视为在真实表面 $\partial M$ 上的均匀采样，每个采样点必须包含 $p$（坐标）和 $\vec{N}$（指向物体内侧的法向）两个属性。定义指示函数 $\chi_{M}:\mathbb{R}^{3}\to\{0,1\}$，描述物体 $M$ 在整个空间中的占用情况，在物体内部指示函数值为 1，否则为 0。 只要能解出这个指示函数，获取边界 $\partial M$ 也就近在咫尺了。

二维泊松重建直观流程，Kazhdan et.al

那么如何求解这个 $\chi_{M}$ 呢？我们注意到指示函数在物体表面法向上的不连续性，在这个前提下考虑指示函数的梯度 $\nabla \chi_{M}$：

  1. 物体内部 $\chi_{M}\equiv1$，所以 $\nabla \chi_{M}=0$
  2. 物体外部 $\chi_{M}\equiv 0$，所以 $\nabla \chi_{M}=0$
  3. 物体表面，$\chi_{M}$ 从 1 突变到 0，$\nabla \chi_{M}$ 的方向为指向物体内部的法线方向，值为无穷大

理想情况下 $\nabla \chi_{M}$ 只在 $\partial M$ 上非零，并且方向与输入点云 $S$ 的法向一致，即在表面附近的采样位置 $s.p$ 上有 $\nabla \chi_{M}\approx\vec{N}$。为了让指示函数可微以使梯度有意义，我们需要对指示函数梯度场进行平滑，由高斯散度定理，平滑后的 $\nabla\tilde{\chi}_{M}$ 等于平滑表面法向得到的向量场：

$$ \nabla(\chi_{M}\,*\,\tilde{F})(q_{0})=\int_{\partial M}\tilde{F}_{p}(q_{0})\vec{N}_{\partial M}(p)dp $$
其中 $\tilde{F}$ 为平滑函数，$*$ 为卷积算符，$q$ 为梯度算子作用的求导变量，$q_{0}$ 为当前评估梯度的位置，$p$ 为积分变量。

由于采样得到的样本点的离散性，无法获得物体表面每点 $p$ 附近的法向，因此将积分转为离散求和，定义 $\mathscr{P}_{s} \subset \partial M$ 为样本点 $s$ 附近的邻域表面，$|\mathscr{P}_{s}|$ 为其面积，且 $\mathscr{P}_{s_{i}}\cap\mathscr{P}_{s_{j}}=\varnothing,\;\cup_{s\in S}\mathscr{P}_{s}=\partial M$，也即 $\{\mathscr{P}_{s_{1}},\mathscr{P}_{s_{2}},\dots,\mathscr{P}_{s_{n}}\}$ 为在 $\partial M$ 上的一个划分。离散化后的平滑操作为：

$$ \begin{aligned} \nabla(\chi_{M}*\tilde{F})(q)&= \sum_{s \in S}\int_{\mathscr{P}_{s}}\tilde{F}_{p}(q)\vec{N}_{\partial M}(p)dp\\ &\approx \sum_{s\in S}|\mathscr{P}_{s}|\tilde{F}_{s.p}(q)s.\vec{N}\equiv\vec{V}(q) \end{aligned} $$
向量场 $\vec{V}$ 为平滑后的指示函数梯度，可以通过平滑表面法向得到。所以最终的问题转化为求解：
$$ \nabla \tilde{\chi}_{M}=\vec{V} $$
直接求解上式需要积分，但 $\vec{V}$ 往往不可积（比如不一定是无旋场），所以通常不存在精确解。那么就退而求其次，使用最小二乘法求近似解，寻找一个标量场 $\phi$，使其梯度 $\nabla \phi$ 为对向量场 $\vec{V}$ 的最佳逼近：
$$ E(\phi)=\min_{\phi}\int||\nabla \phi-\vec{V}||^{2}d\Omega $$
由变分法可证，取最小值时满足泊松方程：
$$ \Delta\tilde{\chi}_{M}=\nabla \cdot \vec{V} $$
同样的，由于 $\chi_{M}$ 的离散性质（只在物体内部为 1，其余为 0），在连续空间中无法解析表达，所以考虑使用一系列基函数的线性组合近似表示：
$$ \chi_{M}=\sum_{i}^n c_{i}F_{i} $$
为此，必须选取合适的基函数以及其所在的函数空间，一个简单的想法是使用 3D 网格，然而这种均匀的体素结构在需要高分辨率的重建任务中需要的函数维度随立方增长，但重建的表面三角形数量却只是平方增长。由于 $\tilde{\chi}_{M}$ 仅在表面 $\partial M$ 附近有值，所以可以考虑一些更为稀疏的表示，自适应八叉树就是一种不错的选择，在 $\partial M$ 附近相对于其它位置深度更深。

定义一棵自适应八叉树 $\mathscr{O}$，其中每个节点 $o\in \mathscr{O}$ 都含有一个函数 $F_{o}$，用这棵八叉树和这些函数来离散化拟合 $\vec{V},\chi_{M}$。

八叉树与基函数 Octree and Basis

在采样点集 $S$ 上，给定最大的八叉树深度 $D$，可以定义一棵最小八叉树，使每个采样点 $s$ 落在某个节点内。

接下来，定义一个固定的、积分为 1 的基函数 $F$，再为每个节点 $o$ 定义经过该节点平移和缩放（即以节点中心为原点，按节点大小缩放）的节点函数 $F_{o}$：

$$ F_{o}(q)\equiv F\left( \frac{q-o.c}{o.w} \right) \frac{1}{{o.w}^{3}} $$
其中 $o.c,o.w$ 分别为节点 $o$ 的中心位置和宽度。由这些节点函数张成的函数空间 $\mathscr{F}_{\mathscr{O},F}\equiv \text{Span}\{F_{o}\}$ 在表达上具有和传统小波的多尺度分辨率特性，在深层节点上拥有高频的函数表达，这些表达在离表面 $\partial M$ 越近时越精确。

为了更有效率地表示向量场 $\vec{V}$，可以把每个采样点的位置替换为其所在叶节点的中心位置，这样每个叶节点 $o$ 只对应一个固定的节点函数 $F_{o}$，每个采样点只对其所在叶节点的节点函数的系数有贡献。通过缩放定义域，可以把 $\tilde{F}$ 定义在一个固定的参考域内（比如 $[-1,1]$），而不需要考虑实际的节点尺度。

$$ F(q)=\tilde{F}\left( \frac{q}{2^{D}} \right) $$
八叉树最大深度 $D$ 对应采样宽度 $2^{-D}$，平滑函数必须近似一个方差数量级为 $2^{-D}$ 的高斯函数（避免过度平滑或数值不稳定），因此 $F$ 应该近似一个单位方差的的高斯函数。论文中作者使用了一个紧支函数（support 有界）来近似：
$$ F(x,y,z)\equiv(B(x)B(y)B(z))^{*n} \; \text{with}\;B(t)= \begin{cases} 1 & |t|<0.5\\ 0 & \text{otherwise} \end{cases} $$
实际上就是盒式滤波做 $n$ 次自卷积，随着 $n$ 增大，$F$ 越来越接近高斯函数并且其 support 范围增加，下图给出了一维情况下盒式滤波函数做 $n$（1-4）次自卷积的形状：

作者最终选择了 $n=3$，此时函数 $F$ 在域 $[-1.5,1.5]^{3}$ 内有 support，对于八叉树的任意一个节点，其节点函数 $F_{o}$ 的 support 最多跨越 $x,y,z$ 三个主轴正负方向各两个邻居节点（$-1.5\leq\frac{q-o.c}{o.w}\leq 1.5\Leftrightarrow q\in[o.c-1.5\cdot o.w,o.c+1.5\cdot o.w]$），即最多有 $5^{3}-1=124$ 个同深度的其它节点的 support 与之重叠，这使得求解散度和拉普拉斯具备了稀疏性。

重写向量场 Vector Field Definition

以上提到的将采样点位置钳制到其所在叶节点的中心位置的方法，显然这样的钳制会造成误差，作者后来使用三线性插值，对与当前采样点 $s$ 最近的八个节点的节点函数值进行插值作为当前采样点的值，最终向量场 $\vec{V}$ 可以写成所有采样点的和：

$$ \vec{V}(q)\equiv \sum_{s\in S}\sum_{o\in\text{Ngbr}_{D}(s)}\alpha_{o,s}F_{o}(q)s.\vec{N} $$
其中，$\text{Ngbr}_{D}(s)$ 为距离当前采样点位置 $s.p$ 最近的八个节点集合，$\alpha_{o,s}$ 为节点 $o$ 的三线性插值权重。由于先前的假设基于采样点分布均匀，因此每个 $\mathscr{P}_{s}$ 的面积 $|\mathscr{P}_{s}|$ 可以视为常数，总体和平滑后的指示函数梯度只差了个缩放系数，不影响重建质量。

但在采样点分布非均匀的情况下，$|\mathscr{P}_{s}|$ 不能视为常数，此时可以根据采样密度自适应的调节每个点的贡献以及核函数宽度，为此首先需要估计一个密度场，并指定 $\hat{D}\leq D$ 为用于估计密度场的八叉树深度：

$$ W_{\hat{D}}(q)\equiv\sum_{s\in S}\sum_{o\in\text{Ngbr}_{\hat{D}}(s)}\alpha_{o,s}F_{o}(q) $$
这种情况下，使用密度场的倒数估计 $|\mathscr{P}_{s}|$：
$$ \vec{V}(q)\equiv \sum_{s\in S} \frac{1}{W_{\hat{D}}(s.p)}\sum_{o\in\text{Ngbr}_{\text{Depth}(s.p)}(s)}\alpha_{o,s}F_{o}(q)s.\vec{N} $$
其中，$\text{Depth}(s.p)(s)$ 代表样本点 $s$ 期望的八叉树深度，通过衡量当前位置的采样密度和平均采样密度 $W$ 并设置：
$$ \text{Depth}(s.p)=\min(D,D+\log_{4}(W_{\hat{D}}(s.p)/W) $$
当 $W_{\hat{D}}(s.p)>W$，即局部密度大于全局密度时，使用最大的八叉树深度（最高精度），保证细粒度的重建；反之，当 $W_{\hat{D}}(s.p)求解泊松方程 Poisson Solution

定义了向量场 $\vec{V}$ 后，我们想要求解 $\tilde{\chi}\in\mathscr{F}_{\mathscr{O},F}$ 以使指示函数的梯度接近 $\vec{V}$，也即求解泊松方程 $\Delta \tilde{\chi}=\nabla\cdot \vec{V}$，但即便 $\tilde{\chi},\vec{V}$ 属于空间 $\mathscr{F}_{\mathscr{O},F}$，$\Delta \tilde{\chi},\nabla \cdot \vec{V}$ 却不一定如此。所以退而求其次，只要求 $\Delta \tilde{\chi}$ 在 $\mathscr{F}_{\mathscr{O},F}$ 上的投影接近于 $\nabla \cdot \vec{V}$ 在 $\mathscr{F}_{\mathscr{O},F}$ 上的投影，也就是解一个最小二乘问题：

$$ \sum_{o \in \mathscr{O}} \left\| \left\langle \Delta\tilde{\chi} - \nabla\cdot\vec{V}, F_o \right\rangle \right\|^2 = \sum_{o \in \mathscr{O}} \left\| \left\langle \Delta\tilde{\chi}, F_o \right\rangle - \left\langle \nabla\cdot\vec{V}, F_o \right\rangle \right\|^2. $$
给定 $|\mathscr{O}|$ 维向量 $v$，其第 $o$ 维坐标 $v_{o}=\left\langle \nabla \cdot \vec{V},F_{o}\right\rangle$，目标为使 $\Delta \tilde{\chi}$ 在每个 $F_{o}$ 上投影组成的向量与 $v$ 尽可能接近。将指示函数写成向量形式 $\tilde{\chi}=\sum_{o}x_{o}F_{o}$，我们要解的就是这个系数向量 $x\in \mathbb{R}^{|\mathscr{O}|}$，然后定义一个 $|\mathscr{O}|\times|\mathscr{O}|$ 的矩阵 $L$，其中第 $(o,o')$ 位置的元素如下：
$$ L_{o,o'}\equiv\left\langle \frac{\partial^{2}F_{o}}{\partial x^{2}},F_{o'}\right\rangle+\left\langle \frac{\partial^{2}F_{o}}{\partial y^{2}},F_{o'}\right\rangle+\left\langle \frac{\partial^{2}F_{o}}{\partial z^{2}},F_{o'}\right\rangle $$
求解 $\tilde{\chi}$ 就相当于求解：
$$ \min_{x\in \mathbb{R}^{\mathscr{O}}}\Vert Lx-v\Vert^{2} $$
而矩阵可以写作：
$$ \begin{aligned} L_{o,o'}=\left\langle \Delta F_{o},F_{o'}\right\rangle&=\int_{\mathbb{R}^{3}}\Delta F_{o}(q)F_{o'}(q)dq\\ &=\int_{\mathbb{R}^{3}}F_{o}(q)\Delta F_{o'}(q)dq\\ &=-\int_{\mathbb{R}^{3}}\nabla F_{o}(q)\nabla F_{o'}(q)dq \end{aligned} $$
显然是对称的，并且由于 $F_{o}$ 是紧支函数，$L$ 中大部分元素都为零，所以 $L$ 也是稀疏的，可以很好的加速求解。

等值面提取 Isosurface Extraction

为了获取重建表面 $\partial \tilde{M}$，首先需要选择一个合适的 isovalue，然后再从指示函数中提取等值面，作者使用了所有采样点位置的指示函数值的平均：

$$ \partial \tilde{M}\equiv\{q\in \mathbb{R}^{3}|\tilde{\chi}(q)=\gamma\}\;\text{with}\;\gamma=\frac{1}{|S|}\sum_{s\in S}\tilde{\chi}(s.p) $$
在采样点分布非均匀的情况下，修改其为采样点位置的指示函数值的加权和：
$$ \partial \tilde{M}\equiv\{q\in \mathbb{R}^{3}|\tilde{\chi}(q)=\gamma\}\;\text{with}\;\gamma=\frac{\sum \frac{1}{W_{\hat{D}}(s.p)}\tilde{\chi}(s.p)}{\sum \frac{1}{W_{\hat{D}}(s.p)}} $$
然后就可以使用 Marching Cubes 等方法提取等值面了。

0x04 Afterwords

泊松重建的算法算是对从 3D 点云中重建水密表面的经典算法，奈何论文作者定义的数学符号比较花哨，对初学者（没错就是我）来说可能比较 “hard to follow”，后续也有许多工作对其进行了各种改进（比如 iPSR 将原始算法对输入点云的有向性依赖去除，以随机初始化的法线经过多次迭代就可重建出高质量表面），但对非水密点云（非闭合表面）的效果仍然没有很好的突破。