35、分段测地线曲线组合用于交互式图像分割

分段测地线曲线组合用于交互式图像分割

1. 引言

在允许用户预先交互的情况下进行图像中目标对象的提取，通常用户会沿着目标边界提供一组地标点，然后计算这些点之间的最小成本路径以构建封闭轮廓。大多数情况下，最小成本路径可以作为相应最小化问题的全局解被高效找到，例如在无负环的有向图中使用 Dijkstra 算法寻找最短路径。基于图像图模型的交互式分割方法本质上是离散的，包括智能剪刀（或实时线）及其动态扩展，以及基于图像森林变换的河床方法。

2. 曲线扭转量化

将曲线 C 分解为不相交的简单环（Seifert 圆），其中一些简单环是顺时针定向的，而另一些是逆时针定向的。我们将曲线 C 的扭转程度量化为逆时针定向环所界定的面积与曲线 C 总面积的比例，这在一定程度上与解开曲线所需的能量相关。
- 简单环和双环 ：简单环由单个反转段描述，涉及单个交点；双环由两个段组成，涉及两个交点。通过之前介绍的两种交叉点类型来检测这些环。
- 简单反转环：通过负交叉点 (u, v) 检测，且从 u 到 v 没有其他交叉点。
- 双反转环：通过负交叉点 (u, v) 检测，且从 u 到 v 的路径包含正交叉点。

设 SL(C) ⊂[0, 1]² 是描述单反转环中涉及的交点的有序对 (u, v) 的集合，DL(C) ⊂[0, 1]²×[0, 1]² 是描述双反转环中涉及的交点对 ((u1, v1), (u2, v2)) 的集合。SL 和 DL 可以在与曲线长度成线性时间内提取。

当曲线 C 从 s 到 t 的部分 C|s→t 是封闭且简单的时，其内部的有符号面积可以用格林定理表示：
[
\int_{\Omega_{in}(C| {s \to t})} dx = \frac{1}{2} \int {s}^{t} x(u)y’(u) - x’(u)y(u)du = \int_{s}^{t} \frac{C(u)^{\perp} \cdot C’(u)}{2} du
]

最终，曲线 C 的反转环总面积 I[C] 可以通过考虑 SL 和 DL 中的所有简单和双反转环来表示：
[
I[C] = - \sum_{(s,t) \in SL(C)} \int_{s}^{t} \frac{C^{\perp} \cdot C’}{2} du - \sum_{((s_1,t_1),(s_2,t_2)) \in DL(C)} \left( \int_{s_1}^{s_2} \frac{C^{\perp} \cdot C’}{2} du + \int_{t_2}^{t_1} \frac{C^{\perp} \cdot C’}{2} du \right)
]

注意，I[C] 是正的，因为每个反转环的有符号面积是负的。

3. 测地线和分段测地线曲线

为了在给定图像 I : D → Rᵈ 中提取结构，有人提出根据由势函数 P : D → R⁺⁺ 定义的异质各向同性度量来寻找最小长度的曲线。在轮廓提取的上下文中，曲线应位于边缘上。因此，势函数定义为 P = g + w，其中 g : D → R⁺ 是图像 I 在某个尺度 s 下梯度幅度的递减函数：
[
g(x) = \frac{1}{1 + |\nabla G_s * I(x)|}
]
其中 Gs 是高斯核，w ∈ R⁺⁺ 是正则化常数。

给定两点 a 和 b，它们之间的最小长度曲线（测地线路径）γₐ,ₓ 可以通过考虑测地距离图（也称为最小作用图）Uₐ : D → R⁺ 来获得，该图为每个点 x ∈ D 分配从 a 到 x 的最小路径长度。这个图是 Eikonal 方程的唯一粘性解，可以通过快速行进（FM）方法在 O(N log N) 操作中高效计算，其中 N 是网格点的数量。一旦距离图被数值计算出来，就可以通过在 Uₐ 上从 b 开始进行梯度下降，直到到达 a 来提取 a 和 b 之间的最小路径。

然而，最小路径方法可能无法提取所需的曲线，例如当 P 太嘈杂或对比度不足、目标曲线长度过长或正则化常数 w 过高时。

为了解决这个问题，一些方法旨在找到分段测地线，即连接一系列连续地标点或顶点的测地线集合。其中，测地连接活动轮廓（GLAC）模型通过连接一组顶点 V = {vᵢ}₁≤ᵢ≤ₙ 构建的测地线路径来生成。给定曲线连接运算符 ⋓：
[
(C_1 \⋓ C_2)(u) =
\begin{cases}
C_1(2u) & \text{if } 0 \leq u \leq 1/2 \
C_2(2u - 1) & \text{if } 1/2 < u \leq 1
\end{cases}
]
只有当 C₁(1) = C₂(0) 时该运算符才有效，连续测地线组装得到的封闭曲线为 Γ = γᵥ₁,ᵥ₂ ⋓ γᵥ₂,ᵥ₃ ⋓ … ⋓ γᵥₙ₋₁,ᵥₙ ⋓ γᵥₙ,ᵥ₁。需要注意的是，测地线的连接本身不是测地线，因此将其称为分段测地线是合适的。

GLAC 方法的目标是找到生成最小化能量泛函 E 的分段测地线曲线的 n 个顶点序列，E 是边缘项和区域项的加权和。通过类似于贪心算法的方法对 E 关于每个 vᵢ 进行数值最小化。

然而，GLAC 模型存在一个显著的缺点，即其初始状态不一定是简单的封闭曲线。当初始顶点在目标边界周围分布不均匀时，可能会出现测地线在目标边界的特定一侧聚集的情况，导致路径重叠。这是因为每条测地线是独立生成的，得到的分段测地线曲线不依赖于相邻顶点对的访问顺序。这种不良现象可能在测地连接轮廓初始化时就出现，或者在先前初始化良好的轮廓经过多次演化步骤后出现。

下面是一个简单的 mermaid 流程图，展示了 GLAC 模型的基本流程：

graph LR
    A[用户提供地标点] --> B[计算测地线路径]
    B --> C[连接测地线路径形成分段测地线]
    C --> D[最小化能量泛函 E]
    D --> E[输出分段测地线曲线]

4. 可允许路径的组合

为了保留分段测地线曲线的优点，我们研究一种更相关的轮廓构建方法。假设对于每对连续顶点都有一组可能的相关路径，我们可以从每组中选择一条路径并将这些路径组合起来，以构建最佳边界曲线。生成轮廓的相关性通过一个能量泛函来衡量，该泛函结合了现有的轮廓和区域项以及一个新的惩罚自相切和自相交的项，确保轮廓是简单的。

4.1 可允许路径集

给定一对连续顶点 (vᵢ, vᵢ₊₁)，设 Aᵢ 是连接这两个顶点的 Ki 条可允许路径的集合：
[
A_i = {\gamma_{i,j}}_{1 \leq j \leq K_i}
]
我们将其称为可允许集。为了生成这些路径，我们提出一种基于鞍点提取的方法。

测地线路径 γₐ,ₓ 可以通过从 a 传播最小作用，在到达 b 时停止，然后从 b 进行梯度下降来提取。也可以从 b 开始反向操作。第三种方法是同时从 a 和 b 传播，在两个传播前沿相遇的第一个位置停止，然后从相遇位置两侧进行梯度下降，并将得到的两条路径适当组合。

当从两个源点进行传播时，得到的组合作用图 Uₐ,ₓ = min(Uₐ, Uₓ)，两个传播前沿在 Uₐ,ₓ 的鞍点相遇。设 mₐ,ₓ : [0, 1] → D 是中曲线，它沿着与 a 和 b 测地距离相等的点延伸，即 {x | Uₐ(x) = Uₓ(x)}。它在组合作用图上形成一个脊，并且只包含临界点，即 ∇Uₐ,ₓ 在 mₐ,ₓ 上未定义。然而，Uₐ,ₓ 可以沿着中曲线进行微分，从而可以将鞍点定义为 Uₐ,ₓ 沿着 mₐ,ₓ 的局部最小值。

为了在离散化的影响下提取稳健的局部最小值，需要对作用图的值沿着 mₐ,ₓ 进行平滑处理。并非所有鞍点都被保留作为路径构建的起点，实际上，我们将每组可允许路径的数量限制在一个固定阈值以下。

下面是一个表格，总结了不同计算测地线路径的方法：
| 方法 | 步骤 |
| ---- | ---- |
| 从 a 传播 | 从 a 传播最小作用，到达 b 停止，从 b 进行梯度下降 |
| 从 b 反向传播 | 从 b 传播最小作用，到达 a 停止，从 a 进行梯度下降并反转路径 |
| 同时传播 | 同时从 a 和 b 传播，在相遇位置停止，两侧进行梯度下降并组合路径 |

4.2 使用简单性能量组合可允许路径

计算可允许的封闭轮廓可以表述为确定最小化能量 E 的标签序列 {x₁, x₂, …, xₙ}，其中标签 xᵢ 对应于集合 Aᵢ 中选择的路径。能量 E 扩展了 GLAC 中涉及的能量，旨在惩罚具有强重叠部分、自相交、与边缘拟合不佳或分割颜色统计不明显区域的轮廓：
[
E[\Gamma] = \Phi[\Gamma] + \omega_{edge} E_{edge}[\Gamma] + \omega_{region} E_{region}[\Gamma]
]
其中 ωₑd₉ₑ 和 ωᵣₑgᵢₒₙ 是用户定义的参数，控制边缘项和区域项相对于简单性项的相对重要性。简单性项涉及在前面定义的归一化自重叠和扭转度量：
[
\Phi[\Gamma] = |Z_{\Gamma}| - \frac{\sqrt{2}}{\sqrt{2}} + \frac{1}{|\Omega_{in}(\Gamma)|} I[\Gamma]
]
边缘项的定义与之前的文献相同，对于区域项，我们使用 Γ 内外颜色概率分布之间的 Bhattacharyya 系数：
[
E_{region}[\Gamma] = \int_{C} \sqrt{p_{in}(\Gamma, \alpha) p_{out}(\Gamma, \alpha)} d\alpha
]
概率分布函数（PDF）对于给定颜色 α 可以使用基于高斯核的直方图进行估计：
[
p_{in}(\Gamma, \alpha) = \frac{1}{|\Omega_{in}(\Gamma)|} \int_{\Omega_{in}(\Gamma)} G_{\sigma}(\alpha - I(x)) dx
]
[
p_{out}(\Gamma, \alpha) = \frac{1}{|\Omega_{out}(\Gamma)|} \int_{\Omega_{out}(\Gamma)} G_{\sigma}(\alpha - I(x)) dx
]

为了确定最佳的标签序列 {x₁, …, xₙ}，暴力搜索的复杂度为 O(Kⁿ)。为了避免测试所有可能的配置，我们提出一种基于路径特定排序的贪心搜索方法，复杂度为 O(n²K)。

在每个可允许集 Aᵢ 中，路径根据外部性 X 递增排序，外部性 X 是给定路径 C 和从 C(1) 返回 C(0) 的线段所形成的有符号面积：
[
X[C] = \frac{1}{2} \int_{0}^{1} C^{\perp} \cdot C’ du + \frac{1}{2} C(1)^{\perp} \cdot C(0)
]
如果将从 C(0) 到 C(1) 的直线作为参考水平轴，当 C 主要在轴下方时，外部性为负；当 C 主要在轴上方时，外部性为正。由于顶点是顺时针排列的，可允许路径从最内层到最外层排序。

从对应于最内部配置的初始标签 {1, …, 1} 开始，根据局部搜索更改标签，通过迭代测试候选标签。在每次迭代中，给定当前标签序列 S，测试与 S 仅相差一个标签且标签仅增加的候选序列。例如，如果当前序列是 {2, 3, 1}，候选序列将是 {3, 3, 1}、{2, 4, 1} 和 {2, 3, 2}。选择导致最小能量的候选序列作为下一次迭代的基础序列，同时在迭代过程中保留全局最小序列。仅通过增加标签生成候选序列使轮廓单调增长。

下面是一个 mermaid 流程图，展示了贪心搜索选择最佳标签序列的过程：

graph LR
    A[初始化标签序列为最内部配置] --> B[计算当前能量 E]
    B --> C[生成候选标签序列]
    C --> D[计算候选序列的能量 E]
    D --> E[选择最小能量的候选序列]
    E --> F[更新当前标签序列]
    F --> G{是否达到终止条件}
    G -- 否 --> C
    G -- 是 --> H[输出最佳标签序列]

5. 实验与讨论

我们通过在自然彩色图像中恢复对象的封闭边界来验证模型的能力，给定沿目标边界的不同数量的用户提供点。在 Grabcut 数据集上进行实验，并与原始的 GLAC 方法（仅考虑分段测地线曲线，无变形）进行比较，以展示使用可允许路径和额外能量项带来的好处。

在参数选择方面，两种方法都在最有利的配置下进行评估，例如正则化权重 w（两种方法都涉及）和能量权重 ωₑd₉ₑ 和 ωᵣₑgᵢₒₙ（专门针对我们提出的算法）每次都分别进行调整，以实现最相关的分割。对于每个图像，还选择了合适的颜色空间（RGB 或更具感知均匀性的 Lab），这会影响两种方法的势函数 P 以及区域能量中涉及的颜色概率分布函数。颜色分量 I 在 [0, 1] 范围内归一化，正则化权重 w 在 0.01 到 0.1 之间变化，能量权重 ωₑd₉ₑ 和 ωᵣₑgᵢₒₙ 都在 0.5 到 2 之间调整。

实验结果表明：
- 香蕉图像 ：在这种特定的初始点配置下，我们的方法没有比 GLAC 方法在分割上有改进。尽管背景复杂且对象包含许多内部边缘，但边界定义清晰，顶点沿边界均匀分布，原始 GLAC 方法能够提取对象。
- 花朵和玩偶图像 ：当提供的顶点较少时，GLAC 方法在测地线之间表现出强烈的重叠。由于 GLAC 方法没有非重叠约束，它会系统地倾向于具有最低势的轮廓部分，因此平滑的边界段会被忽略，而我们的方法更有可能将这些边界段视为可允许路径的一部分，这证明了我们提出的方法对弱边缘的敏感性较低。
- 陶瓷和泰迪熊图像 ：GLAC 方法由于内部边缘比实际边界更强而在对象内部形成捷径，而我们的算法中的简单性和区域能量项有助于解决这个问题，前者防止重叠，后者有利于内部和外部颜色统计差异较大的区域。

下面是一个表格总结不同图像的实验结果：
| 图像 | 原始 GLAC 表现 | 我们的方法表现 |
| ---- | ---- | ---- |
| 香蕉 | 能提取对象 | 无明显改进 |
| 花朵 | 测地线重叠，忽略平滑边界段 | 对弱边缘不敏感，表现更好 |
| 玩偶 | 测地线重叠，忽略平滑边界段 | 对弱边缘不敏感，表现更好 |
| 陶瓷 | 内部形成捷径 | 解决重叠和捷径问题 |
| 泰迪熊 | 内部形成捷径 | 解决重叠和捷径问题 |

综上所述，我们提出的方法在大多数情况下能够克服 GLAC 方法的缺点，提供更准确和可靠的图像分割结果。通过引入可允许路径和新的能量项，我们的方法能够更好地处理复杂的图像场景，尤其是在处理弱边缘和避免路径重叠方面表现出色。

分段测地线曲线组合用于交互式图像分割

6. 方法优势总结

我们提出的方法相较于原始的 GLAC 方法具有多方面的优势，以下通过列表形式详细呈现：
- 处理弱边缘能力更强 ：在花朵和玩偶图像的实验中，GLAC 方法因无重叠约束，会忽略包含稀疏弱边缘的平滑边界段。而我们的方法将这些边界段视为可允许路径的一部分，对弱边缘的敏感性更低，能更好地适应复杂的图像边界。
- 避免路径重叠和捷径 ：在陶瓷和泰迪熊图像中，GLAC 方法因内部强边缘导致在对象内部形成捷径。我们的算法通过简单性和区域能量项，前者防止路径重叠，后者促进内外颜色统计差异大的区域分割，有效解决了这些问题。
- 轮廓简单性保障 ：引入新的能量项惩罚自相切和自相交，确保生成的轮廓是简单的，避免了 GLAC 方法初始状态可能不是简单封闭曲线的问题。

下面通过一个表格对比两种方法的关键特性：
| 方法特性 | 原始 GLAC 方法 | 我们提出的方法 |
| ---- | ---- | ---- |
| 处理弱边缘 | 易忽略含弱边缘的边界段 | 对弱边缘不敏感，可利用弱边缘 |
| 路径重叠和捷径 | 可能出现重叠和捷径 | 有效避免重叠和捷径 |
| 轮廓简单性 | 不能保证初始轮廓简单 | 确保轮廓简单 |

7. 技术细节回顾

为了更好地理解我们的方法，下面对关键技术细节进行回顾：
- 曲线扭转量化 ：将曲线 C 分解为简单环，通过计算逆时针定向环面积与曲线总面积的比例量化扭转程度，相关公式如下：
- 当曲线 C 从 s 到 t 的部分 C|s→t 封闭且简单时，内部有符号面积：
[
\int_{\Omega_{in}(C| {s \to t})} dx = \frac{1}{2} \int {s}^{t} x(u)y’(u) - x’(u)y(u)du = \int_{s}^{t} \frac{C(u)^{\perp} \cdot C’(u)}{2} du
]
- 曲线 C 的反转环总面积 I[C]：
[
I[C] = - \sum_{(s,t) \in SL(C)} \int_{s}^{t} \frac{C^{\perp} \cdot C’}{2} du - \sum_{((s_1,t_1),(s_2,t_2)) \in DL(C)} \left( \int_{s_1}^{s_2} \frac{C^{\perp} \cdot C’}{2} du + \int_{t_2}^{t_1} \frac{C^{\perp} \cdot C’}{2} du \right)
]
- 可允许路径生成 ：基于鞍点提取生成可允许路径，通过同时从两点传播得到组合作用图，在鞍点相遇。中曲线 mₐ,ₓ 上鞍点为 Uₐ,ₓ 的局部最小值，对作用图值沿中曲线平滑处理，限制可允许路径数量。
- 能量泛函设计 ：能量泛函 E 结合现有轮廓和区域项以及新的简单性项：
[
E[\Gamma] = \Phi[\Gamma] + \omega_{edge} E_{edge}[\Gamma] + \omega_{region} E_{region}[\Gamma]
]
其中简单性项：
[
\Phi[\Gamma] = |Z_{\Gamma}| - \frac{\sqrt{2}}{\sqrt{2}} + \frac{1}{|\Omega_{in}(\Gamma)|} I[\Gamma]
]
区域项使用 Bhattacharyya 系数：
[
E_{region}[\Gamma] = \int_{C} \sqrt{p_{in}(\Gamma, \alpha) p_{out}(\Gamma, \alpha)} d\alpha
]
概率分布函数估计：
[
p_{in}(\Gamma, \alpha) = \frac{1}{|\Omega_{in}(\Gamma)|} \int_{\Omega_{in}(\Gamma)} G_{\sigma}(\alpha - I(x)) dx
]
[
p_{out}(\Gamma, \alpha) = \frac{1}{|\Omega_{out}(\Gamma)|} \int_{\Omega_{out}(\Gamma)} G_{\sigma}(\alpha - I(x)) dx
]
- 贪心搜索算法 ：通过对可允许路径按外部性排序，从最内部配置开始，仅增加标签生成候选序列，选择最小能量序列，使轮廓单调增长。

下面是一个 mermaid 流程图，展示整个方法的完整流程：

graph LR
    A[用户提供地标点] --> B[计算测地线路径和可允许路径集]
    B --> C[定义能量泛函 E]
    C --> D[贪心搜索选择最佳标签序列]
    D --> E[组合路径形成封闭轮廓]
    E --> F[输出分割结果]

8. 实际应用场景

我们的方法在多个实际应用场景中具有重要价值：
- 医学图像分割 ：在医学影像中，器官边界可能存在弱边缘，且不同组织间的对比度可能较低。我们的方法能够准确提取器官的封闭边界，为疾病诊断和治疗提供更可靠的依据。例如，在肺部 CT 图像中，准确分割肺部区域对于检测肺部疾病至关重要。
- 遥感图像分析 ：遥感图像通常包含复杂的地形和地物信息，边界可能不清晰。我们的方法可以有效处理这些复杂情况，实现对土地利用类型、植被覆盖等的准确分割，为资源管理和环境监测提供支持。
- 工业检测 ：在工业生产中，对产品的缺陷检测和质量控制需要精确的图像分割。我们的方法能够适应产品表面的各种纹理和边缘特征，准确识别缺陷区域，提高生产效率和产品质量。

9. 未来发展方向

尽管我们的方法在图像分割领域取得了显著的成果，但仍有一些方面值得进一步研究和改进：
- 自适应参数调整 ：目前参数（如正则化权重 w、能量权重 ωₑd₉ₑ 和 ωᵣₑgᵢₒₙ）需要手动调整。未来可以研究自适应参数调整机制，根据图像的特征自动选择最优参数，提高方法的通用性和易用性。
- 多尺度分析 ：引入多尺度分析可以更好地处理不同大小和复杂程度的对象。通过在不同尺度上提取特征和进行分割，然后融合结果，可以进一步提高分割的准确性。
- 与深度学习结合 ：深度学习在图像分割领域取得了巨大的成功。将我们的方法与深度学习技术相结合，利用深度学习的强大特征提取能力和我们方法的优势，可能会产生更强大的图像分割模型。

下面是一个表格，总结未来发展方向及其预期效果：
| 未来发展方向 | 预期效果 |
| ---- | ---- |
| 自适应参数调整 | 提高方法通用性和易用性 |
| 多尺度分析 | 增强对不同大小和复杂对象的分割能力 |
| 与深度学习结合 | 产生更强大的图像分割模型 |

10. 结论

我们提出的基于分段测地线曲线组合的交互式图像分割方法，通过引入可允许路径和新的能量项，有效克服了原始 GLAC 方法的缺点。在处理弱边缘、避免路径重叠和确保轮廓简单性方面表现出色，通过实验验证了在多个图像数据集上的有效性。未来，通过进一步研究自适应参数调整、多尺度分析和与深度学习结合等方向，有望进一步提升方法的性能，为图像分割领域带来更多的突破和应用。

综上所述，我们的方法为图像分割提供了一种更准确、可靠且通用的解决方案，具有广阔的应用前景和发展潜力。