0 概述

特征线是计算机图形学、计算机辅助设计及几何建模等领域广泛使用的高层次信息之一, 主要用于曲面设计和编辑, 网格去噪、参数化和变形等应用中。特征线虽易被人眼识别, 但其自动化提取却是一项具有挑战性的任务。特征是一个宽泛的概念, 抽象来说就是被观测物体的某种可度量的启发式属性。在网格处理中常用的特征线包括:尖锐特征线, 脊线与谷线, 轮廓线等。这些线可通过数学语言来刻画, 目前已出现一些特征线检测方法。尖锐特征线是指曲面法向上发生剧烈变化的曲线。文献[1]通过优化拟合误差及尖锐边的个数、位置来构造含有尖锐特征线的细分曲面控制网, 但该方法耗时较长。文献[2]将类似的技术用于重建点云的尖锐特征, 该方法虽易于实现, 但存在阈值参数选取困难、抗噪性差等问题。脊线与谷线是指曲面上沿主曲率方向取到极值的曲线。文献[3]提出先用光滑的隐式曲面来拟合网格模型, 然后利用网格顶点在隐式曲面的投影点来估计相应的曲率张量和曲率导数, 最后通过追踪算法提取脊线与谷线。基于曲率极值的几何本质, 文献[4]提出快速、有效的脊线与谷线提取算法, 然而, 这类算法的结果不可避免地依赖于高阶微分几何量的估计质量。轮廓线是指法向量与视线向量内积为零的曲线。文献[5]引入启发式轮廓线, 即在曲面上沿视线方向的径向曲率为零, 且在切平面上沿视线投影方向的导数大于零的点集。文献[6]采用主成分分析法和局部二次曲面拟合法对网格模型进行微分几何信息估算, 并结合双阈值检测方法提出网格曲面特征线提取算法。文献[7]将图像处理中的Mumford-Shah模型推广到网格曲面上, 由于其利用跟踪的方法提取特征线, 对含噪声的网格特征线提取存在一定的困难。

稀疏表示与优化理论^[8]在信号处理、图像处理、机器学习和统计学等领域被广泛使用。近年来, 有研究者尝试将该方法应用于数字几何处理中的若干问题^[9-10], 例如保特征去噪^[11-12]、网格分割^[13-14]、曲面重建^[15-17]等, 取得了较好的效果。受此启发, 在特征线具有稀疏分布的情况下, 本文使用L₀优化方法提取网格曲面的特征线, 使用L₀范数来抑制噪声的影响, 以提高提取质量。

1 预备知识 1.1 三角网格表示

设M是给定的网格曲面, 可由三元组(V, E, T)表示, 其中, V={v₁, v₂, …, v_|V|}为顶点集、E={e₁, e₂, …, e_|E|}为边集, T={t₁, t₂, …, t_|T|}为面集, 顶点、边和面的个数分别记为| V |、| E |、| T |。网格的每个面可以是三角形、四边形或其他多边形。为保证简单性与一般性, 本文以三角网格为例。网格曲面的几何信息主要通过顶点的位置来定, 拓扑信息由顶点、边和面之间的连接关系构成。若顶点v为边e的端点, 则记为v≺e。类似地, 若边e为三角形t的一条边, 则记为e≺t。

在网格表示中, 由于每条内边被2个三角形所共享, 为方便离散梯度算子的描述, 需要引入边e相对于其邻接三角形t的方向性概念。不失一般性, 假定M是可定向的, 其每个面t按逆时针方向给出外法向。对M中每条边e随机选定一个方向, e≺t, 若边e的定向和三角形t的定向一致, 则记为sgn(e, t)=1;否则sgn(e, t)=－1。

1.2 离散梯度算子

设f(p)是定义于网格M上的函数, 在M的每个面t上取常数值, 即:

$ f\left( \mathit{\boldsymbol{p}} \right)\left| {_{{t_i}}} \right. \equiv {f_i} \in \mathbb{R},\forall \mathit{\boldsymbol{p}} \in {t_i},i = 1,2, \cdots ,\left| T \right| $

则称f(p)是网格M上的分片常函数。显然, 函数f(p)的梯度在每个面t的内部恒为零, 而在面t的边e上可能取非零值, 记为:

$ \nabla f{\left. {(\mathit{\boldsymbol{p}})} \right|_e} = \left\{ {\begin{array}{*{20}{l}} {\sum\limits_{e \prec t} {{f_t}} \operatorname{sgn} (e,t),e \notin \partial M} \\ {0,e \in \partial M} \end{array}} \right. $

其中, ∂M表示网格M的边界。

不失一般性, M上的分片常函数f(p)可由向量 f =(f₁, f₂, …, f_|T|)^T∈ℝ^|T|表示, 此时离散梯度算子▽ f =(▽f(p)|_e1, ▽f(p)|_e2, …, ▽f(p)|_e|E|)^T可写成矩阵形式, 即▽ f ≜ G f, G ∈ℝ^|E|×|T|, 称 G 为离散梯度算子▽的矩阵表示。

1.3 单位外法向量空间

为叙述方便, 除特别说明外, 下文均以输入度量为单位外法向量进行说明。具体地, 对于给定网格M, 在每个面t上计算或读入一个单位外法向量 n _tⁱⁿ=(n_t1ⁱⁿ, n_t2ⁱⁿ, n_t3ⁱⁿ)^T。若记 N ⁱⁿ=(n ₁ⁱⁿ, n ₂ⁱⁿ, …, n_|T| ⁱⁿ)^T∈ℝ^|T|×3, 则 N ⁱⁿ的每一列可视为M上的分片常函数。对于一般维数的输入度量, 可做类似的处理。

在优化单位外法向量时, 需要保持向量的单位性, 引入M上的单位外法向量空间:

$ \boldsymbol{U}_{N}=\left\{\boldsymbol{N} :\left\|\boldsymbol{n}_{t}\right\|_{2}=1, \quad \forall t \in T\right\} $

其中, N =(n₁, n₂, …, n_|T|)^T∈ ^|T|×3, U _N的特征函数记为:

$ \chi(\boldsymbol{N})=\left\{\begin{array}{l}{0, \boldsymbol{N} \in \boldsymbol{U}_{N}} \\ {+\infty, \boldsymbol{N} \notin \boldsymbol{U}_{N}}\end{array}\right. $

1.4 L₀范数

在稀疏表示与优化理论中, L₀范数常被用于描述向量的稀疏度。对于n维线性空间中的列向量 x =(x₁, x₂, …, x_n)^T∈ℝⁿ, L₀范数定义如下:

$ {\left\| \mathit{\boldsymbol{x}} \right\|_0} = \mathop {\lim }\limits_{p \to 0} \left\| \mathit{\boldsymbol{x}} \right\|_p^p = \mathop {\lim }\limits_{p \to 0} \sum\limits_{k = 1}^n {{{\left| {{x_k}} \right|}^p}} = \# \left\{ {k:{x_k} \ne 0} \right\} $

其中, 符号#表示取集合元素的个数, 即‖ x‖₀等于向量 x 中非零分量的个数。特别地, 当n=1时, 若x₁≠0, 有‖ x ‖₀=1;否则, ‖ x ‖₀=0。

在本文中, 对于m维行向量 y =(y₁, y₂, …, y_m)∈ℝ^1×m, 有:

$ {\left\| \mathit{\boldsymbol{y}} \right\|_0} = {\left\| {\left( {\left| {{y_1}} \right| + \left| {{y_2}} \right| + \cdots + \left| {{y_m}} \right|} \right)} \right\|_0} $

即将行向量 y 视为一个整体, 用于刻画非零向量。对于矩阵 A =(a_ij)∈ℝ^n×m, 将其写成 A =(a₁, a₂, …, a_n)^T的形式, 则有:

$ {\left\| \mathit{\boldsymbol{A}} \right\|_0} = \sum\limits_{i = 1}^n {{{\left\| {\mathit{\boldsymbol{a}}_i^{\text{T}}} \right\|}_0}} $

其中, a _i=(a_i1, a_i2, …, a_im)^T, i=1, 2, …, n。

2 基于L₀优化的特征线提取算法 2.1 L₀优化模型

对于给定的网格M, 输入度量 N ⁱⁿ=(n ₁ⁱⁿ, n ₂ⁱⁿ, …, n _|T|ⁱⁿ)^T, 假设算法输出的度量为 N =(n₁, n₂, …, n_|T|)^T, N 的梯度可表示为:

$ \nabla \mathit{\boldsymbol{N}} = \mathit{\boldsymbol{GN}} \triangleq {\left( {{\mathit{\boldsymbol{u}}_1},{\mathit{\boldsymbol{u}}_2}, \cdots ,{\mathit{\boldsymbol{u}}_{\left| E \right|}}} \right)^{\text{T}}} \in {\mathbb{R}^{\left| E \right| \times 3}} $

其中, u _i=(u_i1, u_i2, u_i3)^T, i=1, 2, …, |E|, 稀疏度表示为${\left\| {\nabla \mathit{\boldsymbol{N}}} \right\|_0} = {\sum\limits_{i = 1}^{|E|} {\left\| {\mathit{\boldsymbol{u}}_i^{\rm{T}}} \right\|} _0}$, 用于刻画输出度量在网格边上发生跃变的个数。稀疏度值越小表示特征的稀疏性越高。本文使用分片常数空间, 输入度量在每个面t上取常值或常向量。假定待提取的特征线位于网格的边上, 在特征线具有稀疏分布的情况下建立如下L₀优化模型:

$ \begin{array}{l} \mathop {\min }\limits_{N \in {U_N}} {\left\| {\nabla N} \right\|_0}\\ {\rm{s}}.\;{\rm{t}}.\;\;\left\| {\mathit{\boldsymbol{N}} - {\mathit{\boldsymbol{N}}^{{\rm{in}}}}} \right\|_{\rm{F}}^2 < \varepsilon \end{array} $

(1)

其中, ‖·‖_F表示矩阵的Frobenius范数, ε为误差阈值。从稀疏表示与优化的角度来看, 该模型在允许的阈值范围内对输入度量进行扰动或去噪, 使得输出度量的梯度具有尽可能低的稀疏度。因此, 该模型不仅能检测出特征线, 还可以有效地抑制噪声的不良影响。容易看出, 通过引入参数α及 U _N的特征函数χ(N), 将式(1)改写成如下形式:

$ \begin{array}{l} \mathop {\min }\limits_\mathit{\boldsymbol{N}} {\left\| {\nabla \mathit{\boldsymbol{N}}} \right\|_0} + \alpha \left\| {\mathit{\boldsymbol{N}} - {\mathit{\boldsymbol{N}}^{{\rm{in}}}}} \right\|_{\rm{F}}^2 + \chi \left( \mathit{\boldsymbol{N}} \right) = \\ \;\;\;\;\;\;\mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{n}}_l}} \right\}} \sum\limits_{e = 1}^{\left| E \right|} {{{\left\| {\mathit{\boldsymbol{u}}_e^{\rm{T}}} \right\|}_0}} + \alpha \sum\limits_{t = 1}^{\left| T \right|} {\left\| {{\mathit{\boldsymbol{n}}_t} - \mathit{\boldsymbol{n}}_t^{{\rm{in}}}} \right\|_2^2} + \chi \left( N \right) \end{array} $

(2)

其中, α>0用于调节各项之间的权重, 值得注意的是, { u _e }仅为中间变量, 可利用{ n _t }表示。

2.2 模型求解

由于式(1)中含有非凸、非连续的L₀范数, 因此求解该模型是NP困难问题, 必须使用近似算法。首先, 通过变量分裂技术在每条边上引入辅助变量 h _e=(h_e1, h_e2, h_e3)^T, 将式(1)写成如下等价形式:

$ \begin{array}{l} \mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{n}}_t},{\mathit{\boldsymbol{h}}_e}} \right\}} \sum\limits_{e = 1}^{\left| E \right|} {{{\left\| {\mathit{\boldsymbol{h}}_e^{\rm{T}}} \right\|}_0}} + \alpha \sum\limits_{t = 1}^{\left| T \right|} {\left\| {{\mathit{\boldsymbol{n}}_t} - \mathit{\boldsymbol{n}}_t^{{\rm{in}}}} \right\|_2^2} + \chi \left( N \right)\\ {\rm{s}}.\;{\rm{t}}.\;\;{\mathit{\boldsymbol{h}}_e} = {\mathit{\boldsymbol{u}}_e},\;\;\;\forall e \in E \end{array} $

(3)

然后, 通过二次罚函数方法将其转化成无约束优化问题, 即引入参数β将式(1)进一步近似为如下形式:

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{n}}_t},{\mathit{\boldsymbol{h}}_e}} \right\}} \sum\limits_{e = 1}^{\left| E \right|} {{{\left\| {\mathit{\boldsymbol{h}}_e^{\rm{T}}} \right\|}_0}} + \alpha \sum\limits_{t = 1}^{\left| T \right|} {\left\| {{\mathit{\boldsymbol{n}}_t} - \mathit{\boldsymbol{n}}_t^{{\rm{in}}}} \right\|_2^2} + }\\ {\beta \sum\limits_{e = 1}^{\left| E \right|} {\left\| {{\mathit{\boldsymbol{u}}_e} - {\mathit{\boldsymbol{h}}_e}} \right\|_2^2} + \chi \left( \mathit{\boldsymbol{N}} \right)} \end{array} $

(4)

其中, β为罚因子。当β不断增大时, 式(4)的解逐步逼近式(3)的解。

最后, 利用交替方向优化方法来求解式(4), 即不断交替固定变量{ n _t }与{ h _e }, 逐步迭代求解直到满足预先设定的停止准则为止, 具体见算法1。

算法1  L₀优化模型求解算法

输入   N ⁱⁿ, α, β₀, β_max, ρ, ε

输出   N

步骤1   k←0, β←β₀, N₀← N ⁱⁿ。

步骤2  重复以下步骤直到$\frac{{\left\| {{\mathit{\boldsymbol{N}}_k} - {\mathit{\boldsymbol{N}}_{k - 1}}} \right\|_{\rm{F}}^2}}{{\left| F \right|}} < M$或者β>β_max:

1) 固定 N _k, 求解H-子问题(5)得到 H _k。

2) 固定 H _k, 求解N-子问题(8)得 N_k+1。

3) β←ρβ, k←k+1。

步骤3  返回 N ← N _k。

在算法1中, N ⁱⁿ是输入网格的初始法向, α是优化模型中的权值, β₀、β_max分别是罚因子的初始值以及最大值, ρ是增长因子, ε是终止误差, N 是优化后的网格面法向。H-子问题与N-子问题的具体求解方法如下:

1) H-子问题。当变量 N ={ n _t}固定时, 式(4)可以简化为:

$ \mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{h}}_e}} \right\}} \sum\limits_{e = 1}^{\left| E \right|} {{{\left\| {\mathit{\boldsymbol{h}}_e^{\rm{T}}} \right\|}_0}} + \beta \sum\limits_{e = 1}^{\left| E \right|} {\left\| {{\mathit{\boldsymbol{u}}_e} - {\mathit{\boldsymbol{h}}_e}} \right\|_2^2} $

(5)

上述优化问题可以进一步分解成若干个独立的子问题, 具体如下:

$ \sum\limits_{e = 1}^{\left| E \right|} {\mathop {\min }\limits_{{\mathit{\boldsymbol{h}}_e}} {{\left\| {\mathit{\boldsymbol{h}}_e^{\rm{T}}} \right\|}_0}} + \beta \left\| {{\mathit{\boldsymbol{u}}_e} - {\mathit{\boldsymbol{h}}_e}} \right\|_2^2 $

(6)

对于这些子问题, 有如下命题:

命题1  式(6)的子问题具有如下显式解:

$ {\mathit{\boldsymbol{h}}_e} = \left\{ \begin{array}{l} {\left( {0,0,0} \right)^{\rm{T}}},{\left\| {{\mathit{\boldsymbol{u}}_e}} \right\|^2} \le 1/\beta \\ {\mathit{\boldsymbol{u}}_e},其他 \end{array} \right. $

(7)

证明  记E(h _e)=‖ h_e ^T‖₀+β‖ u _e－ h_e ‖₂², 下面分2种情形进行讨论:

情形1  当‖ u_e ‖²≤1/β时, 若 h _e≠(0, 0, 0)^T, 则E(h _e)=1+β‖ u _e－ h _e‖²≥β‖ u _e‖²+β‖ u _e－ h _e‖²; 若 h _e=(0, 0, 0)^T, 则E(h _e)=β‖ u _e‖²。容易看出, 当 h _e=(0, 0, 0)^T时, 函数E(h _e)取最小值。

情形2 当‖ u _e‖²>1/β时, 有 u _e≠(0, 0, 0)^T。若 h _e≠(0, 0, 0)^T, 则E(h _e)=1+β‖ u _e－ h _e‖²。在此情形下, 当 h _e= u _e时, E(h _e)取最小值1;若 h _e=(0, 0, 0)^T, 则E(h _e)=β‖ u _e‖²>1。容易看出, 当 h _e= u _e时, 函数E(h _e)取最小值。

综上分析, 证毕。

2) N-子问题。当变量 H ={ h _e}固定时, 式(4)可以简化为:

$ \mathop {\min }\limits_{\left\{ {{\mathit{\boldsymbol{n}}_t}} \right\}} \alpha \sum\limits_{t = 1}^{\left| T \right|} {\left\| {{\mathit{\boldsymbol{n}}_t} - \mathit{\boldsymbol{n}}_t^{{\rm{in}}}} \right\|_2^2} + \beta \sum\limits_{e = 1}^{\left| E \right|} {\left\| {{\mathit{\boldsymbol{u}}_e} - {\mathit{\boldsymbol{h}}_e}} \right\|_2^2} + \chi \left( \mathit{\boldsymbol{N}} \right) $

(8)

其中, χ(N)为非线性项, 本文先将其忽略, 即将式(4)简化为式(9)的二次优化问题。

$ \mathop {\min }\limits_{\left\{ {{n_t}} \right\}} \alpha \sum\limits_{t = 1}^{\left| T \right|} {\left\| {{\mathit{\boldsymbol{n}}_t} - \mathit{\boldsymbol{n}}_t^{{\rm{in}}}} \right\|_2^2} + \beta \sum\limits_{e = 1}^{\left| E \right|} {\left\| {{\mathit{\boldsymbol{u}}_e} - {\mathit{\boldsymbol{h}}_e}} \right\|_2^2} $

(9)

然后将所得的{ n _t }单位化来进行处理。若记 H =(h₁, h₂, …, h_{E})^T∈ℝ^|E|×3, N _i= N (:, i)∈ℝ^|T| , N_i ⁱⁿ= N ⁱⁿ(:, i)∈ℝ^|T|, H _i= H (:, i)∈ℝ^|E|, i=1, 2, 3, 则式(9)可写为:

$ \begin{array}{l} \mathop {\min }\limits_{\left\{ {{n_t}} \right\}} \sum\limits_{i = 1}^3 {\left\{ {\alpha \left[ {\mathit{\boldsymbol{N}}_i^{\rm{T}}{\mathit{\boldsymbol{N}}_i} - 2\mathit{\boldsymbol{N}}_i^{\rm{T}}\mathit{\boldsymbol{N}}_i^{{\rm{in}}} + {{\left( {N_i^{{\rm{in}}}} \right)}^{\rm{T}}}\mathit{\boldsymbol{N}}_i^{{\rm{in}}}} \right] + } \right.} \\ \;\;\;\;\;\;\;\;\;\;\left. {\beta \left( {\mathit{\boldsymbol{N}}_i^{\rm{T}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{N}}_i} - 2\mathit{\boldsymbol{N}}_i^{\rm{T}}{\mathit{\boldsymbol{G}}^{\rm{T}}}{\mathit{\boldsymbol{H}}_i} + \mathit{\boldsymbol{H}}_i^{\rm{T}}{\mathit{\boldsymbol{H}}_i}} \right)} \right\} \end{array} $

(10)

其中, G ∈ℝ^|E|×|T|为离散梯度算子▽的矩阵表示。因此, 式(9)的解满足如下线性方程组:

$ \left( {\alpha \mathit{\boldsymbol{I}} + \beta {\mathit{\boldsymbol{G}}^{\rm{T}}}\mathit{\boldsymbol{G}}} \right){\mathit{\boldsymbol{N}}_i} = \alpha \mathit{\boldsymbol{N}}_i^{{\rm{in}}} + \beta {\mathit{\boldsymbol{G}}^{\rm{T}}}{\mathit{\boldsymbol{H}}_i},i = 1,2,3 $

(11)

其中, I 表示|T|阶单位矩阵。可以发现, 上述线性方程组的系数矩阵是稀疏、对称且正定的, 但每次迭代随β变化而改变, 因此采用双共轭梯度算法^[18]来求解。

2.3 稀疏性提升算法

在式(2)中, 参数α可用于调节特征线的稀疏度与逼近误差之间的权重。当特征线的稀疏度较高时, 许多伪特征边会被提取出来。此时, 可通过减小α来降低特征线的稀疏度, 消除多余的伪特征边。但是, 该方法在实践中并不能得到理想的结果, 其主要原因是减小α相当于增大了式(1)中的ε, 使得输出度量的可扰动范围扩大, 即优化模型可行解的搜索区域变大, 从而使L₀优化模型的求解变得更为困难。因此, 本文提出一种迭代的贪婪策略来逐步降低L₀优化问题的解的稀疏度, 具体地, 在进行第l+1次迭代时, 求解如下改进的L₀优化模型:

$ \mathop {\min }\limits_\mathit{\boldsymbol{N}} {\left\| {\nabla \mathit{\boldsymbol{N}}} \right\|_0} + \alpha \left\| {\mathit{\boldsymbol{N}} - {\mathit{\boldsymbol{N}}^l}} \right\|_{\rm{F}}^2 + \chi \left( \mathit{\boldsymbol{N}} \right),l = 1,2, \cdots $

(12)

其中, N ^l是第l次迭代所得的结果, 详见算法2。

算法2  稀疏性提升算法

输入   N ⁱⁿ, α, β₀, β_max, ρ>1, ε, 迭代次数L

输出  最终的特征线Γ

步骤1   l←0, N ⁰← N ⁱⁿ。

步骤2  重复以下步骤, 直到E^l>‖▽ N ^l－1‖₀或者l≥L。

1) 基于 N ^l、α、β₀、β_max、ρ、ε, 调用算法1求解式(12), 得到解 N ^l+1。

2) E^l+1←‖▽ N ^l+1‖₀+α‖ N ^l+1－ N^l ‖_F², l←l+1。

3) 采用第2.4节的方法提取网格曲面的特征线Γ。

容易看出, 利用算法2计算出的式(12)的解序列{ N ^l }满足‖▽ N ^l+1‖₀ < ‖▽ N^l ‖₀, 即解的稀疏度是递减的。通过选取合适的迭代次数L, 可有效提高特征线的提取质量, 提取结果如图 1所示。

2.4 特征线提取

根据式(1), 本文算法待提取的特征线出现在优化后度量值N发生跃变的位置, 即网格M满足如下条件的边集:

$ \mathit{\Gamma } = \left\{ {e\left| {{{\left\| {\mathit{\boldsymbol{u}}_e^{\rm{T}}} \right\|}_0} \ne 0,e \in E} \right.} \right\} $

其中, u_e ^T=▽ N (e, :), e=1, 2, …, |E|。然而, 由于本文采用近似算法求解式(1), 且受到数值计算误差的影响, 不能直接通过比较‖ u_e ^T‖₀与0来判别特征边。因此, 在算法中, 使用如下标准来提取特征边:

$ \mathit{\Gamma } = \left\{ {e\left| {{{\left\| {\mathit{\boldsymbol{u}}_e^{\rm{T}}} \right\|}_2} > \delta } \right.,e \in E} \right\} $

其中, δ为常数。本文所有实验结果均在δ=0.2情形下取得。

3 实验结果与分析

本文使用C++语言和Matlab实现算法研究, 下文实验均在一台配置Intel Core i5 3.1 GHz CPU及8 GB RAM的计算机上进行。

3.1 参数选取

在算法的执行过程中, 有一些参数需要用户设置, 主要包括α、β₀、β_max、ρ、ε、L等。在本文实验中, 部分参数取默认值, 例如β₀=0.1, ρ=1.4, ε=1.0e－7。其余参数依据具体情况进行设置, 通常都在一定的范围内选择, 例如, α∈[2, 40], β_max∈[5 000 β₀, 20 000 β₀], L∈[1, 12], 具体数值见表 1。

3.2 实验结果

为验证本文算法的有效性, 对不同类型的网格模型进行测试, 包括干净的网格曲面、带噪声的网格曲面、带边界的网格曲面以及带颜色属性的网格曲面。此外, 将本文算法与RTSC-1.5算法^[5]、Crest lines算法^[4]和变分算法^[7]进行对比。为确保公平性, 实验中尽可能按照论文作者提供的建议选取相应参数, 以获得最佳实验结果。不同网格模型的提取结果如图 2~图 7所示。

具体分析如下:

1) 干净的网格曲面

图 2、图 3分别是对于Fusse网格模型和Lucy网格模型采用4种不同算法提取特征线的结果。从局部放大图中容易看出, 本文算法提取的特征线质量明显优于其他算法。与其他3种算法不同, 本文算法无须利用跟踪的方法提取特征线, 因此可避免特征线的后处理过程。

2) 带噪声的网格曲面

为了测试噪声对本文算法结果的影响, 按如下步骤人工合成带噪声的网格模型:

步骤1  随机产生向量集{ d _v=(d_v1, d_v2, d_v3)^T|v∈ V}, 其中, 随机变量d_v1、d_v2、d_v3相互独立且均服从高斯分布N(μ, σ²)= ($\sqrt {2{\rm{ \mathit{ π} }} } $ σ)^－1exp^－(x－μ)2/2σ²。本文实验μ恒取0。

步骤2  更新顶点位置${\mathit{\boldsymbol{\tilde p}}_v} = {\mathit{\boldsymbol{p}}_v} + {\mathit{\boldsymbol{d}}_v}$, ∀v∈V, 其中, p _v为顶点v的位置坐标。

步骤3  按更新后的顶点位置重新计算网格的输入度量 N ⁱⁿ。

图 4为对带σ=0.05e的高斯噪声的Trim-star网格模型进行特征线提取的结果, 其中, e为网格M的平均边长。容易看出, 在处理带噪声的数据时, 本文算法的准确性明显高于其他3种算法。图 5对带σ=0.17e高斯噪声的Dodecahedron网格模型的提取结果, 其进一步显示了本文算法对噪声的鲁棒性。

3) 带边界的网格曲面

图 6为本文算法对带边界的网格模型进行特征线提取的结果, 其中虚线表示网格的边界。从图 6可以看出, 本文算法可有效地处理带边界或带洞的网格模型, 且网格缺失的部分并未对本文算法的结果造成明显的影响。

4) 带颜色属性的网格曲面

为验证本文算法的通用性, 对带RGB颜色属性的网格模型进行测试。在实验中输入的度量改为网格每个面上定义的RGB颜色信息, 如图 7(a)和图 7(c)所示。通过本文的L₀优化算法, 优化后的图像及相应的特征线见图 7(b)和图 7(d)。可以看出, 本文算法对于一般的输入度量也能取得较满意的结果。

3.3 算法运行时间

本文算法主要由2个嵌套的循环构成, 其中外循环为稀疏性提升的迭代过程, 内循环为求解L₀模型的交替方向优化算法。容易看出, 求解一系列大型稀疏线性方程组, 即式(11), 是本文算法主要的性能瓶颈, 占用了大部分运行时间。为了测试算法的性能, 对大型的网格模型进行对比实验, 结果如图 8所示。表 1展示了本文算法对于不同模型所需的计算时间, 可以看出, 本文算法的运行速度还不够理想, 其部分原因是本文算法使用C++语言调用Matlab来求解稀疏线性方程组, 未来可对如何加速进行研究。

4 结束语

由于特征线具有稀疏分布, 本文提出一个基于L₀优化的特征线提取算法。对于给定的输入度量, 通过求解L₀优化问题得到优化后的度量, 使之在网格边上发生跃变的总数极小化, 并使优化后的值与优化前的值之间的差异尽可能小。给出基于变量分裂技术与罚函数方法的交替方向优化算法, 并引入一种迭代策略来提升L₀优化问题解的稀疏性, 以取得更高质量的特征线。实验结果表明, 该算法能有效提取网格曲面的特征线, 与变分算法、Crest lines算法、RTSC-1.5算法相比, 该算法提取质量更高, 对带噪声的数据鲁棒性更好。然而, 本文主要使用参数α来调节特征边的稀疏度与逼近误差之间的权重, 未考虑特征线尺度的问题, 同时本文算法的运行速度不够理想, 因此, 引入层次的优化模型并研究收敛速度更快的L₀优化求解算法将是下一步的研究方向。