CN111028356A - 基于非凸非光滑二阶正则项和稀疏保真项的优化方法 - Google Patents

Abstract

本发明公开了一种基于非凸非光滑二阶正则项和稀疏保真项的优化方法，该方法提出了一个包含非凸非光滑二阶正则项和稀疏保真项的优化模型，用以对几何格网的面法向量进行优化；然后根据优化后的面法向量更新顶点位置进而实现对几何格网的优化。模型中的正则项能够极大化地保留原始格网的几何特征，稀疏保真项能使模型处理多种类型的噪音，如冲量噪音、高斯噪音及混合噪音等。与现有的几何格网模型质量优化技术相比，本发明提出的方法能够显著提升几何格网模型的可视化质量，并且具有计算速度快、计算精度高的优点。

Description

基于非凸非光滑二阶正则项和稀疏保真项的优化方法

技术领域

本发明涉及几何建模领域，更具体地说，涉及一种基于非凸非光滑二阶正则项和稀疏保真项的优化方法，用于对三维网格质量进行优化。

背景技术

目前，几何格网作为一种最常用的三维数据，已经被广泛运用在测绘、地质和地理等领域。近年来，通过激光扫描仪扫描获取及和格网的方式也逐渐普及，但受采集仪器精度以及外界环境等影响，获取到的点云数据不可避免的会含有大量的冗余和噪音，由此生成的几何格网必然也包含噪音。这些噪音不仅会严重影响模型的可视化质量，还会导致在后续的几何处理操作中产生误差，因此，必须对几何格网模型进行优化处理。

长期以来，国内外众多学者对几何格网模型的质量优化问题进行了大量的研究，也提出了很多有效的优化方法，其中典型的有基于滤波器的方法、基于压缩感知的格网优化算法与基于机器学习的方法。虽然这几种方法能达到很好的优化效果，但仍存在一些局限性，主要体现在：

1.基于滤波器的方法，可分为各向同性滤波与各向异性滤波两种类型。各向同性的滤波方法虽然能够有效地优化几何格网模型的质量，但与此同时会使模型的几何特征变得模糊；各向异性的滤波方法在处理高频噪音时也会使模型的几何特征变得模糊。

2.基于压缩感知的格网优化算法，这类方法的基本思想是通过优化全局能量函数来优化几何格网模型的质量，虽然使用该方法能够在质量优化的同时很好地保持模型的几何特征，但会在平滑区域产生阶梯效应的问题。

3.基于机器学习的方法，虽然该方法能够在不考虑潜在曲面的几何特征和噪音模式的前提下有效地优化模型质量，但这种方法的性能取决于训练数据集的完整性，并且十分耗时。

发明内容

针对上述技术问题，本发明提出了一个包含非凸非光滑二阶正则项和稀疏保真项的优化模型，用以对几何格网的面法向量进行优化；然后根据优化后的面法向量更新顶点位置进而实现对几何格网的优化。模型中的正则项能够极大化地保留原始格网的几何特征，稀疏保真项能够使模型处理多种类型的噪音，如冲量噪音、高斯噪音及混合噪音等。由于非凸非光滑模型本身求解的困难性，在求解时，首先需要采用迭代加权L1-范数最小化(IRL1)算法将其转化为一系列凸模型进行求解，然后采用增广拉格朗日乘子法(ALM)和变量分割方法对转化的一系列凸模型进行迭代求解，以得到优化后的面法向量；最后，根据优化后的面法向量更新顶点位置，得到质量优化后的几何格网模型。

具体的，本发明为解决其技术，所采用的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，用于对三维网格质量进行优化，包含如下步骤：

S1、数据准备，具体包括如下子步骤：

S11、获取三维几何格网模型的顶点、边和面，分别用{v_i：i＝1，2，...，V}，{e_j：j＝1，2，...，E}，{τ_k：k＝1，2，...，T}表示，其中V，E和T分别代表三维几何格网模型顶点、边和面的数量；

S12、根据拓扑关系，计算并存储三维几何格网模型中所有顶点、边和面的邻域信息；

S13、根据表达式

计算模型的面法向量，(v_i，v_j，v_k)为三角形τ中按逆时针方向排列的三个顶点；

S2、根据求得的面法向量

设定基于非凸非光滑二阶正则项和稀疏保真项的优化模型：

式中，第一项

为正则项，第二项∑_τarea(τ)||N_τ-Nτ为保真项，

α为优化参数，l是连接三角形τ的重心与某一顶点的线段，len(l)为连接三角形τ的重心与某一顶点的线段l的长度，area(τ)为三角形τ的面积，

为基于线段l的各向异性二阶差分算子，N_τ表示面法向量，

表示的是含有噪音的几何格网中三角形τ的面法向量，p表示p范数，p∈(0，1)；

S3、采用迭代加权L1-范数最小化算法将所述优化模型转化为一系列凸模型；

S4、采用增广拉格朗日乘子法和变量分割方法对转化后的一系列凸模型进行迭代求解，以得到优化后的面法向量；

S5、采用格网重建算法，根据三维几何格网模型的初始顶点坐标与优化后的面法向量更新顶点位置，得到质量优化后的三维几何格网模型。

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S2中，基于线段l的各向异性二阶差分算子定义为：

其中，w_e＝exp(-||N_e，1-N_e，2||²/2σ²)，e＝e⁺、e^-，e⁺和e^-为该顶点的两条临边，τ⁺为基于公共边e⁺和三角形τ相邻的三角形，τ^-为基于公共边e^-和三角形τ相邻的三角形，N_e，1和N_e，2分别表示有公共边e的两个面法向量，σ表示权函数中高斯核的大小。

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S3中，采用迭代加权L1-范数最小化算法将所述优化模型转化为一系列凸模型具体是指基于下述公式迭代求解N^k+1：

式中，

∈是一个防止分母为0的常量。

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S4具体包括：

S41、初始化面法向量N^-1、权重W^-1、第一循环控制变量k以及阈值δ；

S42、通过下述公式迭代求解N^k+1：

具体包括如下子步骤S421-S427：

S421、引入辅助变量X、Y，构造步骤S3中得到的一系列凸模型的增广拉格朗日函数：

式中，

为拉格朗日乘子，r_x、r_y为惩罚因子，

X_l表示l所对应的X的分量，T_τ表示τ所对应的Y的分量，Λ_x，l表示l所对应的Λ_x的分量，Λ_y，τ表示τ所对应的Λ_y的分量；

S422、初始化拉格朗日乘子

N^k，-1＝X^k，-1＝Y^k，-1，第二循环控制变量c；其中，Λ_x、Λ_y、X、Y以及N的右上标的第一位表k相关的变量，第二表示c相关的变量；

S423、固定拉格朗日乘子

和辅助变量X^k，c-1、Y^k，c-1，忽略非线性项ψ(N)，优化法向量N^k，c：

将得到的结果投影到单位球表面，进行归一化处理；

S424、固定拉格朗日乘子

辅助变量Y^k，c-1、法向量N^k，c，求解X^k：

S425、固定拉格朗日乘子

辅助变量X^k，c-1、法向量N^k，c，求解Y^k：

S426、更新拉格朗日乘子

为

为

S427、判断是否满足条件∑_τarea(τ)||N^k，c-N^k，c-1||＜δ或迭代次数c是否大于或等于第二预设次数，若是则令当前法向量N^k等于N^k，c，进入步骤S43，否则返回步骤S423；

S43、根据当前法向量

更新权重

式中∈是一个防止分母为0的常量：

S44、判断是否满足条件∑_τarea(τ)||N^k-N^k-1||＜δ或迭代次数k是否大于或者等于第一预设次数，若是则输出优化的法向量，否则将k更新为k+1并返回步骤S42。

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S424中进行求解时：

由于X中的每个元素X_l相互独立，因此，逐一求解每个X_l：

X_l有封闭解

式中ζ_x，l表示l对应的ζ_x的分量，

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S425中进行求解时：

由于Y中的每个元素T_τ相互独立，因此，逐一求解每个T_τ：

T_τ有封闭解

式中ζ_y，τ表示τ对应的ζ_y的分量，

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S423中所述将得到的结果投影到单位球表面，进行归一化处理是通过Eigen库中的LDLT方法实现。

进一步地，在本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法中，步骤S5中采用的格网重建算法为：

式中(v_i，v_j，v_k)表示三角形τ中按逆时针方向排列的三个顶点，

表示含有噪音的几何格网的顶点位置。

实施本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，具有以下有益效果：本发明公开了一种有效的几何格网模型质量优化方法，该方法提出了一个包含非凸非光滑二阶正则项和稀疏保真项的优化模型，用以对几何格网的面法向量进行优化；然后根据优化后的面法向量更新顶点位置进而实现对几何格网的优化，模型中的正则项能够极大化地保留原始格网的几何特征，稀疏保真项能使模型处理多种类型的噪音，如冲量噪音、高斯噪音及混合噪音等。与现有的几何格网模型质量优化技术相比，本发明提出的方法能够显著提升几何格网模型的可视化质量，并且具有计算速度快、计算精度高的优点。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法的流程图；

图2为三角形τ基于线段l(蓝色表示)的各向异性二阶差分，图中红色的圆点为三角形τ的重心；

图3(a)为含有冲量噪音的模型图，图3(b)为本发明方法的优化结果图；

图4(a)为含有高斯噪音的模型图，图4(b)为本发明方法的优化结果图；

图5(a)为含噪音的CAD模型图，图5(b)为本发明方法的优化结果图；

图6(a)为含噪音的non-CAD模型图，图6(b)为本发明方法的优化结果图；

图7(a)为含噪音的扫描模型图，图7(b)为本发明方法的优化结果图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

参考图1，图1是本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法的流程图。本发明的基于非凸非光滑二阶正则项和稀疏保真项的优化方法用于对三维网格质量进行优化，包含如下步骤：

S1、数据准备，具体包括如下子步骤：

S11、获取三维几何格网模型的顶点、边和面，分别用{v_i：i＝1，2，...，V}，{e_j：j＝1，2，...，E}，{τ_k：k＝1，2，...，T}表示，其中V，E和T分别代表三维几何格网模型顶点、边和面的数量；

S12、根据拓扑关系，计算并存储三维几何格网模型中所有顶点、边和面的邻域信息；

S13、根据表达式

计算模型的面法向量，(v_i，v_j，v_k)为三角形τ中按逆时针方向排列的三个顶点。

S2、根据求得的面法向量

设定基于非凸非光滑二阶正则项和稀疏保真项的优化模型：

式中，第一项

为正则项，第二项∑_τarea(τ)||N_τ-Nτ为保真项，

α为优化参数，l是连接三角形τ的重心与某一顶点的线段，len(l)为连接三角形τ的重心与某一顶点的线段l的长度，area(τ)为三角形τ的面积，

为基于线段l的各向异性二阶差分算子，N_τ表示面法向量，

表示的是含有噪音的几何格网中三角形τ的面法向量，p表示p范数，p∈(0，1)。

参考图2，基于线段l的各向异性二阶差分算子定义为：

其中，w_e＝exp(-||N_e，1-N_e，2||²/2σ²)，e＝e⁺、e^-，e⁺和e^-为该顶点的两条临边，τ⁺为基于公共边e⁺和三角形τ相邻的三角形，τ^-为基于公共边e^-和三角形τ相邻的三角形，N_e，1和N_e，2分别表示有公共边e的两个面法向量，σ表示权函数中高斯核的大小。

S3、采用迭代加权L1-范数最小化算法(IRL1)将所述优化模型转化为一系列凸模型。由于步骤S2的优化模型本身求解的困难性，在求解时采用IRL1将其转化为一系列图模型，基于下述公式迭代求解N^k+1：

式中，

∈是一个防止分母为0的常量。

S4、采用增广拉格朗日乘子法和变量分割方法对转化后的一系列凸模型进行迭代求解，以得到优化后的面法向量。

S5、采用格网重建算法，根据三维几何格网模型的初始顶点坐标与优化后的面法向量更新顶点位置，得到质量优化后的三维几何格网模型。步骤S5中采用的格网重建算法为：

式中(v_i，v_j，v_k)表示三角形τ中按逆时针方向排列的三个顶点，

表示含有噪音的几何格网的顶点位置。

关于上述步骤，步骤S4具体包括：

S41、初始化面法向量N^-1＝0、权重W^-1＝0、第一循环控制变量k＝0以及阈值δ＝1e-6；

S42、通过下述公式迭代求解N^k+1：

具体包括如下子步骤S421-S427：

S421、引入辅助变量X、Y，构造步骤S3中得到的一系列凸模型的增广拉格朗日函数：

式中，

为拉格朗日乘子，r_x、r_y为惩罚因子，

X_l表示l所对应的X的分量，Y_τ表示τ所对应的Y的分量，Λ_x，l表示l所对应的Λ_x的分量，Λ_y，τ表示τ所对应的Λ_y的分量。

S422、初始化拉格朗日乘子

N^k，-1＝X^k，-1＝Y^k，-1，第二循环控制变量c＝0；其中，Λ_x、Λ_y、X、Y以及N的右上标的第一位表k相关的变量，第二表示c相关的变量。

S423、固定拉格朗日乘子

和辅助变量X^k，c-1、Y^k，c-1，忽略非线性项ψ(N)，优化法向量N^k，c：

将得到的结果投影到单位球表面，进行归一化处理，该过程可通过Eigen库中的LDLT方法实现。

S424、固定拉格朗日乘子

辅助变量Y^k，c-1、法向量N^k，c，求解X^k：

由于X中的每个元素X_l相互独立，因此，逐一求解每个X_l：

X_l有封闭解

式中ζ_x，l表示l对应的ζ_x的分量，

S425、固定拉格朗日乘子

辅助变量X^k，c-1、法向量N^k，c，求解Y^k：

由于Y中的每个元素Y_τ相互独立，因此，逐一求解每个Y_τ：

Y_τ有封闭解

式中ζ_y，τ表示τ对应的ζ_y的分量，

S426、更新拉格朗日乘子

为

为

S427、判断是否满足条件∑_τarea(τ)||N^k，c-N^k，c-1||＜δ或迭代次数c是否大于或等于第二预设次数，例如50次，若是则令当前法向量N^k等于N^k，c，进入步骤S43，否则返回步骤S423。

S43、根据当前法向量

更新权重

式中∈是一个防止分母为0的常量：

S44、判断是否满足条件∑_τarea(τ)||N^k-N^k-1||＜δ或迭代次数k是否大于或者等于第一预设次数，例如20次，若是则输出优化的法向量，否则将k更新为k+1并返回步骤S42。

采用C++语言实现上述算法，代码环境配置具体为：Microsoft VisualStudio2010；CGAL-4.9；Eigen-3.2.8。参数格式为：

HOMF<input><output>(<α><β><r><NITRS><VITRS>)其中，<input><output>分别为输入的噪声网格和输出的结果网格，(<α><β><r><NITRS><VITRS>)为算法运行参数。

本发明已经过严格的测试和验证，图3(a)、图3(b)、图4(a)、图4(b)分别展示了本发明方法对含有冲量噪音和高斯噪音的模型的质量优化效果，图5(a)、图5(b)、图6(a)、图6(b)和图7(a)、图7(b)分别展示了本发明方法在CAD模型、non-CAD模型和扫描模型上的质量优化效果。表明了该方法不仅能够在处理多种类型噪音的同时极大化地保留原始格网的几何特征，而且还不会产生阶梯效应，这显著地提升了几何格网模型的可视化质量。

上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (8)

1.一种基于非凸非光滑二阶正则项和稀疏保真项的优化方法，用于对三维网格质量进行优化，其特征在于，包含如下步骤：

S1、数据准备，具体包括如下子步骤：

S11、获取三维几何格网模型的顶点、边和面，分别用{v_i：i＝1，2，...，V}，{e_j：j＝1，2，...，E}，{τ_k：k＝1，2，...，T}表示，其中V，E和T分别代表三维几何格网模型顶点、边和面的数量；

S12、根据拓扑关系，计算并存储三维几何格网模型中所有顶点、边和面的邻域信息；

S13、根据表达式

计算模型的面法向量，(v_i，v_j，v_k)为三角形τ中按逆时针方向排列的三个顶点；

S2、根据求得的面法向量

(设定基于非凸非光滑二阶正则项和稀疏保真项的优化模型：

式中，第一项

为正则项，第二项∑_τ area(τ)||N_τ-Nτ为保真项，CN＝N∈R3×T：

α为优化参数，l是连接三角形τ的重心与某一顶点的线段，len(l)为连接三角形τ的重心与某一顶点的线段l的长度，area(τ)为三角形τ的面积，

为基于线段l的各向异性二阶差分算子，N_τ表示面法向量，

表示的是含有噪音的几何格网中三角形τ的面法向量，p表示p范数，p∈(0，1)；

S3、采用迭代加权L1-范数最小化算法将所述优化模型转化为一系列凸模型；

S4、采用增广拉格朗日乘子法和变量分割方法对转化后的一系列凸模型进行迭代求解，以得到优化后的面法向量；

S5、采用格网重建算法，根据三维几何格网模型的初始顶点坐标与优化后的面法向量更新顶点位置，得到质量优化后的三维几何格网模型。

2.根据权利要求1所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S2中，基于线段l的各向异性二阶差分算子定义为：

其中，w_e＝exp(-||N_e，1-N_e，2||²/2σ²)，e＝e⁺、e^-，e⁺和e^-为该顶点的两条临边，τ⁺为基于公共边e⁺和三角形τ相邻的三角形，τ^-为基于公共边e^-和三角形τ相邻的三角形，N_e，1和N_e，2分别表示有公共边e的两个面法向量，σ表示权函数中高斯核的大小。

3.根据权利要求1所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S3中，采用迭代加权L1-范数最小化算法将所述优化模型转化为一系列凸模型具体是指基于下述公式迭代求解N^k+1：

式中，

∈是一个防止分母为0的常量。

4.根据权利要求3所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S4具体包括：

S41、初始化面法向量N^-1、权重W^-1、第一循环控制变量k以及阈值δ；

S42、通过下述公式迭代求解N^k+1：

具体包括如下子步骤S421-S427：

S421、引入辅助变量X、Y，构造步骤S3中得到的一系列凸模型的增广拉格朗日函数：

式中，Λ_x，

为拉格朗日乘子，r_x、r_y为惩罚因子，

X_l表示l所对应的X的分量，Y_τ表示τ所对应的Y的分量，Λ_x，l表示l所对应的Λ_x的分量，Λ_y，τ表示τ所对应的Λ_y的分量；

S422、初始化拉格朗日乘子

N^k，-1＝X^k，-1＝Y^k，-1，第二循环控制变量c；其中，Λ_x、Λ_y、X、Y以及N的右上标的第一位表k相关的变量，第二表示c相关的变量；

S423、固定拉格朗日乘子

和辅助变量X^k，c-1、Y^k，c-1，忽略非线性项ψ(N)，优化法向量N^k，c：

将得到的结果投影到单位球表面，进行归一化处理；

S424、固定拉格朗日乘子

辅助变量Y^k，c-1、法向量N^k，c，求解X^k：

S425、固定拉格朗日乘子

辅助变量X^k，c-1、法向量N^k，c，求解Y^k：

S426、更新拉格朗日乘子

为

为

S427、判断是否满足条件∑_τarea(τ)||N^k，c-N^k，c-1||＜δ或迭代次数c是否大于或等于第二预设次数，若是则令当前法向量N^k等于N^k，c，进入步骤S43，否则返回步骤S423；

S43、根据当前法向量

更新权重

式中∈是一个防止分母为0的常量：

S44、判断是否满足条件∑_τ area(τ)||N^k-N^k-1||＜δ或迭代次数k是否大于或者等于第一预设次数，若是则输出优化的法向量，否则将k更新为k+1并返回步骤S42。

5.根据权利要求4所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S424中进行求解时：

由于X中的每个元素X_l相互独立，因此，逐一求解每个X_l：

X_l有封闭解

式中ζ_x，l表示l对应的ζ_x的分量，

6.根据权利要求4所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S425中进行求解时：

由于Y中的每个元素Y_τ相互独立，因此，逐一求解每个Y_τ：

Y_τ有封闭解

τ，式中ζ_y，τ表示τ对应的ζ_y的分量，

7.根据权利要求4所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S423中所述将得到的结果投影到单位球表面，进行归一化处理是通过Eigen库中的LDLT方法实现。

8.根据权利要求1所述的基于非凸非光滑二阶正则项和稀疏保真项的优化方法，其特征在于，步骤S5中采用的格网重建算法为：

式中(v_i，v_j，v_k)表示三角形τ中按逆时针方向排列的三个顶点，

表示含有噪音的几何格网的顶点位置。

Images