Poisson 表面重建的核心思想

🧩 一、问题背景

在点云重建中，我们通常有：

• 一组三维点：$\{p_i\}$，在物体表面上；

• 每个点的法向量：$\{n_i\}$。

目标是：

从这些点和法向量重建一个连续的表面 $S$。

🧭 二、Poisson重建的核心思想

Poisson Surface Reconstruction (Kazhdan et al., 2006, 2007) 的核心想法是：

将点云的法向信息看作是某个隐函数的梯度场，然后通过求解泊松方程来恢复这个隐函数。

换句话说：

$$\nabla \chi = \mathbf{v}$$

其中：

• $\chi(\mathbf{x})$：体素网格上的一个 隐式函数（indicator function）；

• $\mathbf{v}(\mathbf{x})$：表示法向场的 向量场。

🧮 三、什么是 $\chi(\mathbf{x})$？

我们定义一个函数 $\chi(\mathbf{x})$ 表示物体内部外部关系：

$$\chi(\mathbf{x}) = \begin{cases} 1, & \text{if } \mathbf{x} \text{ 在物体 } M \text{ 内部}\\ 0, & \text{if } \mathbf{x} \text{ 在物体 } M \text{ 外部} \end{cases}$$

这个函数叫做 indicator function（指示函数） 或 characteristic function。

⚙️ 四、从法向到泊松方程

我们知道在表面上，法向量方向与 $\nabla \chi$ 一致：

$$\mathbf{n} = \frac{\nabla \chi}{\|\nabla \chi\|}$$

因此有：

$$\nabla \chi \approx \mathbf{v}$$

（$\mathbf{v}$ 是点云法向插值得到的连续向量场。）

对上式取散度：

$$\nabla \cdot \nabla \chi = \nabla \cdot \mathbf{v}$$

这就变成 Poisson 方程：

$$\Delta \chi = \nabla \cdot \mathbf{v}$$

🔍 五、求解与提取表面

解出 $\chi(\mathbf{x})$ 后，我们取等值面（通常是 $\chi = 0.5$）：

$$S = \{ \mathbf{x} \mid \chi(\mathbf{x}) = \tau \}$$

这就是重建出的表面。

🧾 六、总结

Poisson 表面重建与有限元（FEM）方法之间的关系

Poisson Surface Reconstruction 本质上是用 有限元方法（Finite Element Method, FEM） 离散求解一个 Poisson 方程。
这听起来很数学，但其实直观地说就是：

“用一堆局部的小积木（基函数）拼出一个连续函数，并让它在平均意义上满足方程。”

一、为什么要用有限元？

我们要求解的偏微分方程（PDE）为(Poisson方程）：

$$\begin{equation} \Delta \chi = \nabla \cdot \mathbf{v} \end{equation}$$

其中 $\Delta$ 是拉普拉斯算子（包含二阶导数）。

其中：

• $\chi(\mathbf{x})$ 是隐函数；
• $\mathbf{v}(\mathbf{x})$ 是由点云法线方向定义的向量场。

计算机不能处理连续函数，只能处理有限维向量，因此必须离散化这个 PDE。
常见的两种离散化思路如下：

有限元法可以看作差分法的更灵活版本——它允许网格大小与形状根据数据密度自适应调整。

二、有限元的基本思想

设要求解的未知函数为 $\chi(\mathbf{x})$。
我们用一组基函数 $\{\phi_i(\mathbf{x})\}$ 近似它：

$$\begin{equation} \chi(\mathbf{x}) \approx \sum_i c_i \, \phi_i(\mathbf{x}), \end{equation}$$

其中：

• $\phi_i(\mathbf{x})$ 为局部形函数（shape function）；
• $c_i$ 为对应系数（待求未知量）；
• 每个 $\phi_i$ 仅在局部单元上非零，因此称为“分片函数（piecewise function）”。

例如在一维中，每个基函数形如三角形：

       /\
      /  \
-----/----\-----
   i-1   i   i+1

在三维中，这些“积木”则对应于立方体、四面体或八叉树节点上的局部基。

三、从强形式到弱形式

Poisson 方程的强形式为：
$$\Delta \chi = \nabla \cdot \mathbf{v}.$$

直接对每个点要求方程成立（强形式）会遇到问题：

• $\Delta \chi$ 要求函数二阶可导；
• 而有限元中的分片线性函数二阶导数为零。

因此我们转化为“弱形式（weak form）”，即方程在积分意义上成立。

对任意测试函数 $\psi(\mathbf{x})$，要求：

$$\begin{equation} \int_\Omega \psi \, \Delta \chi \, d\mathbf{x} = \int_\Omega \psi \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x} \end{equation}$$

这称为“弱形式”或“变分形式”（Weak / Variational Form）。

应用分部积分（Integration by Parts） 可将左式中的二阶导数转为一阶导数。下面是详细过程。

我们利用下面的分部积分（integration by parts）推导详见附录A
$$\begin{equation} \boxed{\displaystyle \int_\Omega ψ\,\Delta χ \,d\mathbf{x} = \int_{\partial\Omega} ψ\,\frac{\partial χ}{\partial n}\,dS \;-\; \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x}} \end{equation}$$
对左边进行替换，得到：

$$\int_{\partial\Omega} ψ\,\frac{\partial χ}{\partial n}\,dS \;-\; \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x} = \int_\Omega \psi \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x}$$

进一步整理后，得到：
$$\begin{equation} \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x} = - \int_\Omega \psi \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x} + \int_{\partial\Omega} ψ\,\frac{\partial χ}{\partial n}\,dS \end{equation}$$

右边第二项就是边界项，假设边界项消失（通常取 Dirichlet 边界或自然边界条件），就得到：

$$\begin{equation} \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x} = - \int_\Omega \psi \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x} \end{equation}$$

这就是 Poisson 方程的弱形式。
这个形式有好处：只含一阶导数，更容易离散。

四、离散化并形成线性系统

令 $\chi = \sum_j c_j \phi_j$，取测试函数 $\psi = \phi_i$，代入公式（6）：

$$\int_{\Omega} \nabla \phi_i \cdot \sum_j c_j \nabla \phi_j \, d\mathbf{x} = - \int_{\Omega} \phi_i \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x}.$$

整理为矩阵形式：

$$\sum_j \underbrace{\left( \int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j \, d\mathbf{x} \right)}_{A_{ij}} c_j = \underbrace{- \int_{\Omega} \phi_i \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x}}_{b_i},$$

即：

$$A \, \mathbf{c} = \mathbf{b}.$$

其中矩阵 $A$ 对称、正定，因此可通过共轭梯度（CG）或多重网格等高效方法求解。

五 左右两边具体如何数值计算

$$\sum_j \underbrace{\left( \int_{\Omega} \nabla \phi_i \cdot \nabla \phi_j \, d\mathbf{x} \right)}_{A_{ij}} c_j = \underbrace{- \int_{\Omega} \phi_i \, (\nabla \cdot \mathbf{v}) \, d\mathbf{x}}_{b_i}$$

1.前提与符号

• $\{\phi_i\}$：全局基（形函数）。每个 $\phi_i$ 有局部支撑（只在若干单元/节点非零）。

• 网格 / 细分单元集合：$\mathcal{T}=\{K\}$。每个单元上可以进行局部积分。

• 目标：组装稀疏刚度矩阵 $A$（$A_{ij}=\int_\Omega \nabla\phi_i\cdot\nabla\phi_j$) 与右端向量 $b$（如前述 $b_i=-\int_\Omega \phi_i(\nabla\cdot v)$ 或等价 $\int_\Omega \nabla\phi_i\cdot v$ 的形式）。

• 最终解线性系统 $A\mathbf{c}=\mathbf{b}$。

通常做法（FEM 风格）：

• 对每个单元 $K$ 计算局部刚度矩阵 $A^K_{ab}=\int_K \nabla\phi_a\cdot\nabla\phi_b\,dV$（这里 a,b 是单元的局部基索引）；

• 把局部矩阵加到全局 $A$（local→global mapping）；

• 对 $b$ 做同样的局部积分与装配。

2.左边：如何计算 $A_{ij}=\int_\Omega \nabla\phi_i\cdot\nabla\phi_j\,d\mathbf{x}$

通用公式（单元分解）

$$A_{ij}=\sum_{K\in\mathcal T}\int_{K} \nabla\phi_i\cdot\nabla\phi_j\,dV$$

只有当 $\operatorname{supp}\phi_i$ 与 $\operatorname{supp}\phi_j$ 在某个单元 $K$ 重叠时该单元才贡献非零，导致 $A$ 稀疏。这个在有原论文中使用basic 样条函数，n=3,支持域是1.5，所以一个Voxel最多有124个neighbor来求积分。

常见情形 1 — 线性单元（tet/tri / 八叉树 piecewise linear）

• 对于一个线性三角形/四面体单元，形函数在单元内是一阶多项式，梯度是常数。

• 因此局部矩阵很简单：

  $$A^K_{ab} = (\nabla\phi_a|_K)\cdot(\nabla\phi_b|_K) \cdot |K|$$

  其中 $|K|$ 是单元体积（或面积），$\nabla\phi_a|_K$ 是常数矢量（可由单元顶点坐标直接算出）。

如何算梯度（四面体例子）：

• 在四面体上，若局部基是“节点为1其余为0”的 hat 函数，梯度可由单元逆雅可比计算（标准 FEM 教材）。

• 计算完每个单元的 $A^K$ 后，把它加到全局矩阵的对应 (global_i, global_j)。

常见情形 2 — 八叉树 / B-spline / higher-order bases

• 若 $\phi$ 是 B-spline/box spline，$\nabla\phi$ 在每个支撑子域可能不是常数，但通常有封闭表达或可以按单元用高斯积分求出：

  $$A^K_{ab}\approx \sum_{q} w_q \, (\nabla\phi_a(x_q)\cdot\nabla\phi_b(x_q))$$

  $x_q$ 为 Gauss 点，$w_q$ 含雅可比权重。

3、右边，如何计算 $b_i = -\int_\Omega \phi_i(\nabla\!\cdot\!v)\,dV$（两种等价实现）

如前面所述，有两种常用做法（更稳定的是第二种）：

方案 A（直接）：先在网格上构建 $\nabla\cdot v$，再积分

1. 在每个单元（或节点）上评估/近似散度 $\nabla\cdot v$（差分或在单元上插值后解析求散度）。

2. 做单元高斯积分（或把散度视作单元常数）：

   $$b_i \approx -\sum_{K}\sum_{q} w_q \, \phi_i(x_q) \, (\nabla\cdot v)(x_q)$$

3. 装配进全局 $b$。

缺点：若 $v$ 是由稀疏点云插值得到且噪声较大，直接求散度不稳定。

方案 B（推荐，等价但数值更好）：使用散度的分部积分，把散度移到基函数上

利用恒等式（把散度换成 $\nabla\phi_i\cdot v$ 加边界项）：

$$-\int_\Omega \phi_i(\nabla\cdot v)\,dV = \int_\Omega \nabla\phi_i\cdot v\,dV \;-\; \int_{\partial\Omega} \phi_i(v\cdot n)\,dS.$$

通常边界项可忽略或单独处理（域足够大或设边界为零）。于是

$$b_i \approx \int_\Omega \nabla\phi_i\cdot v \,dV.$$

这一步非常实用：你只需能在给定点评估 $v(x)$（通过 splatting/RBF/插值），而不必先计算其散度。

局部积分：

$$b_i = \sum_{K}\int_K \nabla\phi_i(x)\cdot v(x)\,dV \approx \sum_{K}\sum_{q} w_q \,\nabla\phi_i(x_q)\cdot v(x_q)$$

六、原文的理解

1️⃣ 指示函数与梯度

• 指示函数（indicator function）本身梯度在表面是无限大，无法直接使用。

• 解决方法是用一个 平滑核（如 Gaussian）去卷积指示函数，得到一个光滑的标量场 $F(\mathbf{x})$。

• 根据数学原理，卷积后梯度可以近似由 点云的法向量分布来表示：

  $$\nabla F(\mathbf{x}) \approx \text{某种基于法向量的向量场}$$

2️⃣ 八叉树表示

• 八叉树的作用是 自适应空间分辨率：点云密度高的地方用更多的节点。

• 每个节点（顶点）对应一个未知系数 $c_i$，表示该位置的函数值。

• 八叉树内部是用基函数的三线性（trilinear）表示，从而可以方便插值， 基函数用 B-spline 表示

3️⃣ 构造梯度场

• 因为真实表面未知，需要在每个八叉树节点附近 splat 点云 patch：

  • 每个点云点（带法向量）会对邻近的节点产生影响，形成 离散向量场。

• 数学上，用 B-spline 基函数的梯度和法向量的内积积分来得到每个节点的贡献：

  $$\mathbf{v}_o = \sum_{p} (\mathbf{n}_p \cdot \nabla B_o(\mathbf{x}_p))$$

• 这就是构造右手边 $b$ 的过程。

4️⃣ 刚度矩阵 $A$

• 刚度矩阵 $A$ 来自 拉普拉斯算子作用在 B-spline 基函数上：

  $$A_{oo'} = \langle \nabla B_o, \nabla B_{o'} \rangle$$

• 本质上是 Poisson 方程的离散化：

  $$\Delta F = \nabla \cdot \mathbf{V} \quad \Rightarrow \quad A \mathbf{c} = \mathbf{b}$$

5️⃣ 求解系数

• 最终解线性系统 $A c = b$，得到每个节点的系数 $c_i$。

• 解出的 $c_i$ 表示 光滑函数的值，然后可以用 marching cubes 或等值面提取得到表面。

其他说明

一、边界项的物理/数学含义

边界项 $\int_{\partial\Omega} ψ\,\frac{\partial χ}{\partial n}\,dS$ 表示“流”（或通量）穿过边界对弱等式的贡献：

• 如果你固定了 χ 在边界（Dirichlet），测试函数通常取为零在边界（消去边界项）。
• 如果你指定了法向导数（Neumann）值，则该边界项会包含已知的边界贡献并进入右端项 $b$。
• 在 Poisson 重建常见做法：把域设大并把 χ 在边界置 0（或用自然边界），实际计算中通常能安全忽略或处理该项。

二、Poisson Surface Reconstruction 的有限元实现

在 Kazhdan 等人的方法（2006, 2007）中：

• 基函数 $\phi_i(\mathbf{x})$ 选用 分片三线性 B-spline；
• 网格结构采用 八叉树（Octree），可自适应细化；
• 向量场 $\mathbf{v}$ 来源于输入点云的法线；
• 方程求得的 $\chi(\mathbf{x})$ 是隐函数场；
• 最终通过提取等值面 $\chi = \text{iso}$ 得到三维表面。

这一过程正是一个 有限元 Poisson 方程求解：

FEM 离散 → 解稀疏线性系统 → 取等值面。

三、与有限差分法的对比

因此，FEM 更适合稀疏、不规则点云的 Poisson 重建。

四、Poisson 图像融合的关系

Poisson 图像融合（Poisson Image Editing）同样解的是 Poisson 方程：

$$\Delta I = \nabla \cdot \mathbf{v},$$

只是定义域变成了 二维图像域，边界条件来自已知图像像素。
Poisson Surface Reconstruction 则是在 三维体素域 上进行相同的数学操作。
两者共享完全相同的 PDE 基础，只是应用场景与基函数不同。

五、直观理解与总结

你可以把 FEM 想象成：

用乐高积木（基函数）拼出一个连续曲面，
调节每块的高度（系数 $c_i$），
让整体尽量满足 Poisson 方程。

最终，有限元法让 Poisson Surface Reconstruction 具备：

• 数学上严格的 PDE 背景；
• 自适应空间分辨率；
• 高效的稀疏矩阵求解特性；
• 更鲁棒的噪声容忍性。

六、小结

有限元方法（FEM）通过基函数近似、弱形式求解，使得连续 Poisson 方程可在离散网格上高效求解。
Poisson Surface Reconstruction 正是利用 FEM 思想在三维空间上重建连续隐函数，从而生成高质量的表面。

参考文献：

• Kazhdan, M., Bolitho, M., & Hoppe, H. (2006). Poisson Surface Reconstruction.
• Kazhdan, M. & Hoppe, H. (2013). Screened Poisson Surface Reconstruction.
• Zienkiewicz, O. C., The Finite Element Method: Its Basis and Fundamentals.
• Botsch, M. et al., Polygon Mesh Processing.

附录A：分部积分

1 对任意标量场 $ψ(\mathbf{x})$ 和 $χ(\mathbf{x})$，有下面恒等式（乘积求散度的向量恒等式） 附录C：

$$\nabla\cdot(ψ\,\nabla χ) = \nabla ψ \cdot \nabla χ + ψ\,\Delta χ$$

这是从乘积求导（类似一维的 $(fg)'=f'g+fg'$）在向量形式下的推广。
把它重排得到：

$$ψ\,\Delta χ = \nabla\cdot(ψ\nabla χ) - \nabla ψ\cdot\nabla χ$$

2 对区域积分并用散度定理（Divergence theorem）

把上式在区域 $\Omega$ 上积分：

$$\int_\Omega ψ\,\Delta χ \,d\mathbf{x} = \int_\Omega \nabla\cdot(ψ\nabla χ)\,d\mathbf{x} \;-\; \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x}$$

利用散度定理（Gauss theorem）附录B把第一个体积分换成边界积分：

$$\int_\Omega \nabla\cdot(ψ\nabla χ)\,d\mathbf{x} = \int_{\partial\Omega} ψ\,(\nabla χ\cdot \mathbf{n})\,dS,$$

其中 $\partial\Omega$ 是域的边界，$\mathbf{n}$ 是外法向，$\nabla χ\cdot\mathbf{n}=\frac{\partial χ}{\partial n}$（法向导数）。

因此合并得：

$$\boxed{\displaystyle \int_\Omega ψ\,\Delta χ \,d\mathbf{x} = \int_{\partial\Omega} ψ\,\frac{\partial χ}{\partial n}\,dS \;-\; \int_\Omega \nabla ψ\cdot\nabla χ\,d\mathbf{x}}$$

可以看到，分部积分后出现边界项＋内积项。

分部积分的意义在于我们把一个二阶方程转换成一阶方程。

附录B 散度定理

对于一个任意的光滑向量场 $\mathbf{F}$：

$$\int_{\Omega} \nabla \cdot \mathbf{F} \, dV = \int_{\partial \Omega} \mathbf{F} \cdot \mathbf{n} \, dS$$

其中：

• $\Omega$：积分域（比如三维空间中的一个体积）
• $\partial \Omega$：这个体积的边界面
• $\mathbf{n}$：边界面上外法向量
• $\mathbf{F}$：任意向量场

左边是一个$\nabla \cdot \mathbf{F}$体积积分（volume integral）
$\nabla \cdot \mathbf{F}$（读作 divergence of F）表示这一点是不是“流的源头或汇点”。

• 如果 $\nabla \cdot \mathbf{F} > 0$：
  表示流体从这个点“源源不断流出” —— 有点像喷泉口；
• 如果 $\nabla \cdot \mathbf{F} < 0$：
  表示流体“流入”这个点 —— 像吸进去；
• 如果 $\nabla \cdot \mathbf{F} = 0$：
  表示进出相等，局部守恒。

右边是一个$\mathbf{F} \cdot \mathbf{n}$边界面积积分（surface integral）
在边界面上，$\mathbf{F} \cdot \mathbf{n}$表示在该点上，流体“穿出”表面的量

• 如果 F 与 n 同方向 → 表示流出；
• 如果反方向 → 表示流入；
• 如果垂直 → 表示不穿出。

整个散度定理的意思是：

“体积内的散度积分 = 表面积上通量积分”

直觉解释一句话总结

“所有点内部流出的总量 = 实际从表面流出去的量”

或者说：

“体积里的源汇之和 = 边界上的通量之和”

这就是散度定理。

附录C 乘积求散度的向量恒等式

这个公式是：

$$\nabla \cdot (\psi \nabla \chi) = \nabla \psi \cdot \nabla \chi + \psi \, \Delta \chi$$

我们要弄懂每个符号是什么意思，然后一步步推导。

符号说明

推导（从最基本定义出发）

我们从三维坐标展开：

$$\nabla \cdot (\psi \nabla \chi) = \frac{\partial}{\partial x}(\psi \frac{\partial \chi}{\partial x}) + \frac{\partial}{\partial y}(\psi \frac{\partial \chi}{\partial y}) + \frac{\partial}{\partial z}(\psi \frac{\partial \chi}{\partial z})$$

用一维乘积法则展开每一项：

$$\frac{\partial}{\partial x}(\psi \frac{\partial \chi}{\partial x}) = \frac{\partial \psi}{\partial x} \frac{\partial \chi}{\partial x} + \psi \frac{\partial^2 \chi}{\partial x^2},$$

同理对 $y, z$ 方向。

把所有加起来：

$$\nabla \cdot (\psi \nabla \chi) = \left( \frac{\partial \psi}{\partial x} \frac{\partial \chi}{\partial x}+ \frac{\partial \psi}{\partial y} \frac{\partial \chi}{\partial y}+ \frac{\partial \psi}{\partial z} \frac{\partial \chi}{\partial z} \right)+ \\ \psi \left(\frac{\partial^2 \chi}{\partial x^2} + \frac{\partial^2 \chi}{\partial y^2} + \frac{\partial^2 \chi}{\partial z^2} \right)$$

用向量符号重新写：

$$\boxed{ \nabla \cdot (\psi \nabla \chi) = \underbrace{\nabla \psi \cdot \nabla \chi}_{\text{梯度点积项}} + \underbrace{\psi \, \Delta \chi}_{\text{乘上拉普拉斯项}} }$$