数值解法入门：水平集函数的偏导数与曲率离散化实现（GAC/ACWE专用）

前情回顾：在之前的分享《从 PDE 到形态学：解锁高效稳定的曲线与曲面演化新范式》中，我们重点讲解了 GAC/ACWE 算法的形态学解法——用离散形态学算子（膨胀/腐蚀、SIₙ/ISₙ）实现轮廓演化，逻辑简单、无需复杂数学推导。

今天，我们聚焦算法的另一种核心实现方式——数值解法。我们将核心拆解“水平集函数离散化 → 偏导数计算 → 曲率离散化”的完整流程，帮助你快速上手数值解法的核心代码实现，尤其适合需要深入理解水平集理论的新手。

本文所有内容均围绕学术论文中的标准表述展开，结合公式解析、直观理解和代码示例，全程对应“水平集函数 $u$ → 一阶偏导数 → 二阶偏导数 → 曲率”的递进逻辑。吃透这篇，就能轻松搞定数值解法中最核心的“空间离散化”环节。

一、基础铺垫：水平集函数的离散化表示

数值解法与形态学解法的核心区别之一，就是水平集函数 $u$ 的存储形式。形态学解法用布尔值（0/1）表示像素是否在轮廓内，而数值解法需要存储更精细的符号距离信息，因此离散化方式有所不同：

与形态学实现方式一致，我们采用二维数组对水平集函数 $u$ 进行离散化。其中，每个单元格存储一个浮点数值。我们采用中心差分格式对偏微分方程（PDEs）进行空间离散化。网格间距固定为 1，因此一阶导数由公式(33)和(34)给出。

关键细节解析

离散化载体：二维数组（对应图像的像素网格， $i$ 表示行坐标、 $j$ 表示列坐标，与图像坐标系完全一致）；
存储类型：浮点数值（区别于形态学的布尔值）——因为数值解法中， $u$ 通常是符号距离函数（SDF），需要存储像素到轮廓的“带符号距离”（轮廓内正、轮廓外负、轮廓上为 0）；
空间离散化方法：中心差分格式（二阶精度，比前向/后向差分更准确，是数值解法的标准选择）；
网格间距：固定为 1（工程简化惯例，无需额外处理物理间距的缩放，直接用像素坐标计算即可）。

二、第一步：用中心差分计算一阶偏导数（ $u_x$ 、 $u_y$ ）

水平集函数的梯度（ $\nabla u = (u_x, u_y)$ ）是曲率计算、轮廓演化速度计算的基础。而在数值解法中，梯度的分量（一阶偏导数 $u_x$ 、 $u_y$ ）通过中心差分离散计算，对应论文中的公式(33)和(34)：

公式定义（2D 场景）

$u_x(i, j) = \frac{1}{2}(u(i, j+1) - u(i, j-1)) \tag{33}$
$u_y(i, j) = \frac{1}{2}(u(i+1, j) - u(i-1, j)) \tag{34}$

逐公式解析

1. $u_x(i, j)$ ：x 方向（列方向/水平方向）一阶偏导数

计算逻辑：当前像素 $(i, j)$ 的右侧像素 $(i, j+1)$ 减去左侧像素 $(i, j-1)$ ，结果除以 2；
物理意义：描述 $u$ 在水平方向的变化率——值越大，水平方向的轮廓变化越剧烈；
示例：若 $u(i, j+1) = 1.5$ （右侧离轮廓较近）、 $u(i, j-1) = 0.5$ （左侧离轮廓较远），则 $u_x = 0.5$ ，说明水平方向上，像素从左到右逐渐靠近轮廓。

2. $u_y(i, j)$ ：y 方向（行方向/垂直方向）一阶偏导数

计算逻辑：当前像素 $(i, j)$ 的下侧像素 $(i+1, j)$ 减去上侧像素 $(i-1, j)$ ，结果除以 2；
物理意义：描述 $u$ 在垂直方向的变化率——值越大，垂直方向的轮廓变化越剧烈；
关键提醒： $i$ 是行坐标（对应 y 轴），因此 $i+1$ 是下一行、 $i-1$ 是上一行，不要与常规数学坐标系的 x/y 混淆。

优势

中心差分是“前向差分 + 后向差分”的平均，具备二阶精度（前向/后向差分仅一阶），且计算简单——仅需相邻 2 个像素的差值，兼顾精度和效率，是数值解法中偏导数计算的首选。

三、第二步：二阶偏导数离散化（ $u_{xx}$ 、 $u_{yy}$ 、 $u_{xy}$ ）

曲率的显式计算公式需要用到水平集函数的二阶偏导数（ $u_{xx}$ 、 $u_{yy}$ ）和混合二阶偏导数（ $u_{xy}$ ）。结合网格间距 $= 1$ 的简化条件，具体离散公式如下：

1. 二阶偏导数公式（2D 场景）

（1） $u_{xx}(i, j)$ ：x 方向二阶偏导数

$u_{xx}(i, j) = u(i, j+1) + u(i, j-1) - 2u(i, j)$

（2） $u_{yy}(i, j)$ ：y 方向二阶偏导数

$u_{yy}(i, j) = u(i+1, j) + u(i-1, j) - 2u(i, j)$

（3） $u_{xy}(i, j)$ ：混合二阶偏导数

$u_{xy}(i, j) = \frac{1}{4} \left( u(i+1, j+1) - u(i-1, j+1) - u(i+1, j-1) + u(i-1, j-1) \right)$

2. 逐公式解析

（1） $u_{xx}$ 与 $u_{yy}$ ：离散逻辑

本质是“二阶差分”的简化形式，反映“偏导数的变化率”：

记忆技巧： $u_{xx}$ 对应“列方向”的相邻像素（ $i, j \pm 1$ ）， $u_{yy}$ 对应“行方向”的相邻像素（ $i \pm 1, j$ ）。

[!NOTE]
请注意坐标对应关系：在程序实现中， $i$ 通常指代行（垂直）， $j$ 指代列（水平）。因此 $u_{xx}$ （关于水平方向的二阶导）涉及的是 $j \pm 1$ 。

（2） $u_{xy}$ ：混合二阶偏导数

描述“x 方向偏导数随 y 方向的变化”，核心是“对角线像素的差值组合”：

计算逻辑：取当前像素的四个对角线像素，按“右下 - 右上 - 左下 + 左上”的思想（具体符号取决于网格定义）计算，再乘以 $1/4$ 。

3. 边界处理提醒

图像边界像素由于缺少邻域，无法直接计算。工程中的常用方案：

直接将边界像素的二阶偏导数赋值为 0；
用“前向/后向差分”替代中心差分。

四、第三步：曲率的离散化计算（数值解法核心）

在水平集方法中，轮廓的曲率是“平滑项”的核心。曲率描述轮廓的弯曲程度。

1. 曲率的连续数学定义

二维轮廓（零水平集 $u=0$ ）的曲率 $\kappa$ ，标准定义为：
$\kappa = \text{div} \left( \frac{\nabla u}{\|\nabla u\|} \right)$

通俗解释：曲率是“轮廓单位法向量的散度”。

2. 曲率的显式公式（偏导数展开式）

将其展开为一/二阶偏导数的组合，得到数值解法中直接使用的显式公式：
$\kappa = \frac{u_{xx}u_y^2 - 2u_{xy}u_x u_y + u_{yy}u_x^2}{(u_x^2 + u_y^2)^{3/2}}$

符号解析

符号	含义
$\kappa$	当前像素处的曲率值
$u_x, u_y$	一阶偏导数
$u_{xx}, u_{yy}, u_{xy}$	二阶偏导数
$(u_x^2 + u_y^2)^{3/2}$	梯度模长的三次方

3. 公式简化（工程常用）

当 $u$ 是符号距离函数（SDF）时，满足 $\|\nabla u\| = 1$ ，分母变为 1，公式简化为：
$\kappa = u_{xx}u_y^2 - 2u_{xy}u_x u_y + u_{yy}u_x^2$

4. 曲率值的物理意义

曲率 $\kappa$ 的符号	轮廓形态	含义
$\kappa > 0$	凸轮廓	绝对值越大，轮廓越“尖”
$\kappa < 0$	凹轮廓	绝对值越大，凹陷越深
$\kappa = 0$	直线/平面	无弯曲，轮廓平滑

五、完整实现：曲率离散化的 Python 代码示例

import numpy as np

def compute_derivatives(u):
    """
    计算水平集函数 u 的一阶偏导数和二阶偏导数
    """
    u = u.astype(np.float64)
    h, w = u.shape
    
    u_x = np.zeros_like(u)
    u_y = np.zeros_like(u)
    u_xx = np.zeros_like(u)
    u_yy = np.zeros_like(u)
    u_xy = np.zeros_like(u)
    
    # 计算一阶偏导数 (中心差分)
    u_x[:, 1:-1] = 0.5 * (u[:, 2:] - u[:, :-2])
    u_y[1:-1, :] = 0.5 * (u[2:, :] - u[:-2, :])
    
    # 计算二阶偏导数
    u_xx[:, 1:-1] = u[:, 2:] + u[:, :-2] - 2 * u[:, 1:-1]
    u_yy[1:-1, :] = u[2:, :] + u[:-2, :] - 2 * u[1:-1, :]
    
    # 计算混合偏导数
    u_xy[1:-1, 1:-1] = 0.25 * (u[2:, 2:] - u[0, 2:] - u[2:, 0] + u[0:, 0]) # 示例逻辑
    # 修正索引逻辑以匹配标准网格
    # u_xy[i,j] = 0.25 * (u[i+1,j+1] - u[i-1,j+1] - u[i+1,j-1] + u[i-1,j-1])
    term = u[2:, 2:] - u[:-2, 2:] - u[2:, :-2] + u[:-2, :-2]
    u_xy[1:-1, 1:-1] = 0.25 * term
    
    return u_x, u_y, u_xx, u_yy, u_xy

def compute_curvature(u):
    """
    离散化计算曲率（基于 SDF 简化公式）
    """
    u_x, u_y, u_xx, u_yy, u_xy = compute_derivatives(u)
    
    # 简化公式
    curvature = u_xx * (u_y ** 2) - 2 * u_xy * u_x * u_y + u_yy * (u_x ** 2)
    
    # 边界置 0
    curvature[0, :] = 0.0; curvature[-1, :] = 0.0
    curvature[:, 0] = 0.0; curvature[:, -1] = 0.0
    
    return curvature

六、常见问题与注意事项

网格间距的影响：若网格间距 $h \neq 1$ ，需在公式中加入相应的 $1/h$ 或 $1/h^2$ 缩放。
SDF 的重要性：简化公式仅在 $u$ 为 SDF 时成立。若 $u$ 出现退化，必须使用完整公式。

七、总结

从“水平集离散化基础”到“一/二阶偏导数计算”，再到“曲率离散化实现”，我们完整覆盖了数值解法中最关键的步骤。掌握了这些，你就为后续完整实现 GAC/ACWE 算法的数值迭代打下了坚实基础。

如果你对曲率在轮廓演化中的具体应用感兴趣，欢迎继续关注！

从 PDE 到形态学：高效稳定的曲线与曲面演化（7）

从 PDE 到形态学：高效稳定的曲线与曲面演化（7）

数值解法入门：水平集函数的偏导数与曲率离散化实现（GAC/ACWE专用）

一、基础铺垫：水平集函数的离散化表示

关键细节解析

二、第一步：用中心差分计算一阶偏导数（ $u_x$ 、 $u_y$ ）

公式定义（2D 场景）

逐公式解析

1. $u_x(i, j)$ ：x 方向（列方向/水平方向）一阶偏导数

2. $u_y(i, j)$ ：y 方向（行方向/垂直方向）一阶偏导数

优势

三、第二步：二阶偏导数离散化（ $u_{xx}$ 、 $u_{yy}$ 、 $u_{xy}$ ）

1. 二阶偏导数公式（2D 场景）

（1） $u_{xx}(i, j)$ ：x 方向二阶偏导数

（2） $u_{yy}(i, j)$ ：y 方向二阶偏导数

（3） $u_{xy}(i, j)$ ：混合二阶偏导数

2. 逐公式解析

（1） $u_{xx}$ 与 $u_{yy}$ ：离散逻辑

（2） $u_{xy}$ ：混合二阶偏导数

3. 边界处理提醒

四、第三步：曲率的离散化计算（数值解法核心）

1. 曲率的连续数学定义

2. 曲率的显式公式（偏导数展开式）

符号解析

3. 公式简化（工程常用）

4. 曲率值的物理意义

五、完整实现：曲率离散化的 Python 代码示例

六、常见问题与注意事项

七、总结

相关阅读更多精彩内容

友情链接更多精彩内容

从 PDE 到形态学：高效稳定的曲线与曲面演化（7）

数值解法入门：水平集函数的偏导数与曲率离散化实现（GAC/ACWE专用）

一、基础铺垫：水平集函数的离散化表示

关键细节解析

二、第一步：用中心差分计算一阶偏导数（、）

公式定义（2D 场景）

逐公式解析

1. ：x 方向（列方向/水平方向）一阶偏导数

2. ：y 方向（行方向/垂直方向）一阶偏导数

优势

三、第二步：二阶偏导数离散化（、、）

1. 二阶偏导数公式（2D 场景）

（1）：x 方向二阶偏导数

（2）：y 方向二阶偏导数

（3）：混合二阶偏导数

2. 逐公式解析

（1） 与 ：离散逻辑

（2）：混合二阶偏导数

3. 边界处理提醒

四、第三步：曲率的离散化计算（数值解法核心）

1. 曲率的连续数学定义

2. 曲率的显式公式（偏导数展开式）

符号解析

3. 公式简化（工程常用）

4. 曲率值的物理意义

五、完整实现：曲率离散化的 Python 代码示例

六、常见问题与注意事项

七、总结

相关阅读更多精彩内容

友情链接更多精彩内容

二、第一步：用中心差分计算一阶偏导数（ $u_x$ 、 $u_y$ ）

1. $u_x(i, j)$ ：x 方向（列方向/水平方向）一阶偏导数

2. $u_y(i, j)$ ：y 方向（行方向/垂直方向）一阶偏导数

三、第二步：二阶偏导数离散化（ $u_{xx}$ 、 $u_{yy}$ 、 $u_{xy}$ ）

（1） $u_{xx}(i, j)$ ：x 方向二阶偏导数

（2） $u_{yy}(i, j)$ ：y 方向二阶偏导数

（3） $u_{xy}(i, j)$ ：混合二阶偏导数

（1） $u_{xx}$ 与 $u_{yy}$ ：离散逻辑

（2） $u_{xy}$ ：混合二阶偏导数