0 概述

在日常生活中, 三维形状无处不在, 各行各业都涉及三维形状的应用, 包括形状识别^[1-2]、形状分类^[3-4]以及形状检索^[5-6]等。这些应用的基础研究工作就是形状匹配^[7-8], 其本质是一个形状相似性问题, 通过选择合适的相似性度量来定义不同形状的相似程度^[9-10]。

在进行形状匹配时, 要选定一种在不同变换下数值不变的形状描述符, 并计算其数值。同时, 需要选定一种合适的距离, 将一对形状的描述符差值作为其形状相似度。

形状描述符是一种描述形状语义和几何信息的方法, 其对等距、噪声或拓扑等变化鲁棒。形状描述符大致可分为局部描述符和全局描述符2类。局部描述符通常描述形状的局部语义信息, 而全局描述符则描述了形状的整体信息。常用的形状描述符包括:形状上下文(Shape Context, SC)^[11], 梯度直方图^[12]和CSHOT算子^[13]。这些形状描述符在形状发生平移及旋转等刚性变换时仍能保持不变。但是在现实生活中, 除刚性变换外, 3D形状更易发生等距或者近似等距变换。等距变换是指保持曲面上任意曲线长度不变的变换, 例如人弯曲胳膊就是一种等距离变换。在形状发生等距变换后, 其不变性特征称为等距不变性^[14]。

为衡量形状发生等距变换后的相似度, 定义一种谱形状描述符, 并将形状上任意2点的谱形状描述符的距离作为谱距离, 研究表明谱距离具有等距不变性。谱形状描述符主要包括:全局点签名(Global Point Signature, GPS)^[15], 热核签名(Heat Kernel Signature, HKS)^[16]以及双调和签名(Biharmonic Signature, BS)^[17]等。与之相对应的谱距离有通勤时间距离^[18]、扩散距离^[19]以及双调合距离^[20]。在这些谱距离中, 双调合距离能同时调和形状的局部特性和全局特性, 且对参数自由, 不受参数的影响。文献[21]提出的FMPS是一种设计精巧且算法性能较高的匹配算法, 通过计算一对3D非刚性形状之间的泛函映射, 进行3D非刚性形状匹配, 但该算法的时间耗费较高。因此, 本文在双调合距离的基础上, 提出一种三维形状的相似性度量方法。

1 双调合距离

双调合距离是定义在形状上的一种谱距离。在计算双调合距离之前, 首先要计算形状表面的拉普拉斯-贝尔特拉米(Laplace-Baltrami, LB)算子。

1.1 LB算子

LB算子是Laplace算子在黎曼流形上的推广。Laplace算子被定义为形状表面的实值函数f的梯度散度值, 其表达式如下^[22]:

$ \Delta f = \nabla \cdot \nabla f = {\nabla ^2}f = \frac{{{\partial ^2}f}}{{\partial {x^2}}} + \frac{{{\partial ^2}f}}{{\partial {y^2}}} + \frac{{{\partial ^2}f}}{{\partial {z^2}}} $

(1)

将Laplace算子推广到黎曼流形上得到LB算子, 其在局部坐标系中可表示为:

$ \Delta f = \nabla \cdot \nabla f = \frac{1}{{\sqrt G }}\sum\limits_{i,j = 1}^n {{g^{ij}}} \frac{\partial }{{\partial {x^i}}}\left( {\sqrt G {g^{ij}}\frac{{\partial f}}{{\partial {x^i}}}} \right) $

(2)

其中, n为形状的顶点个数, g为度量张量, G为矩阵[g^ij]的行列式。因为本文研究的形状格式为三角网格, 所以三角网格上离散LB算子的定义如下^[23]:

$ Lf\left( {{v_i}} \right) = \sum\limits_{j = 1}^n {{w_{ij}}} \left( {f\left( {{v_i}} \right) - f\left( {{v_j}} \right)} \right) $

(3)

其中, w_ij为系数矩阵中的元素, f(v_i)和f(v_j)分别为作用于函数f的离散LB算子在顶点v_i和v_j的取值, 则式(3)的具体计算过程如下:

$ \begin{array}{l} Lf\left( {{v_i}} \right) = \frac{\partial }{{\partial f\left( {{v_i}} \right)}}\sum\limits_{{v_i} \in Neigh \left( {{v_i}} \right)} {{E_D}\left( {{f_{ith - tri}}} \right)} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\sum\limits_{{v_i} \in Neigh \left( {{v_i}} \right)} {\left( {\cot {\alpha _j} + \cot {\beta _j}} \right)} \cdot \left| {f\left( {{v_i}} \right) - f\left( {{v_j}} \right)} \right| \end{array} $

(4)

其中, α_j、β_j分别表示连接顶点v_i、v_j与边e_ij的对角, Neigh(v_i)表示与v_i相邻的顶点集合。将LB算子进行谱分解, 计算过程如下:

$ {\Delta _M}{\varphi _i} = {\lambda _i}{\varphi _i} $

(5)

其中, λ_i和φ_i分别是LB算子第i个特征值及其对应的特征向量。

1.2 距离计算

测地距离是描述形状属性常用的一种距离, 其测量了曲面上两点之间的最短路径长度。但是, 测地距离只能反映形状的局部属性, 不能反映全局属性, 同时, 其对拓扑噪声敏感, 引入任意小的拓扑变化可能导致较大的测地距离变化。文献[20]通过正则化Laplace算子的特征值, 提出一种新的距离, 即双调和距离。该距离能够很好地平衡形状的局部特征和全局特征, 同时对拓扑变化具有鲁棒性和等距不变性。双调和距离来源于双调和微分方程, 其本质上与通勤时间距离类似。双调和距离的具体思想是将形状映射到一个无限维谱空间, 每一维代表形状上某一点的特征, 则形状上所有顶点的BS值表示为一个n×i维矩阵, 其中, n为形状上采样点的数量, i为LB算子的特征值和特征向量数量, 则BS矩阵可定义为:

$ \mathit{\boldsymbol{BS}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\varphi _1}(1)}}{{{\lambda _1}}}}&{\frac{{{\varphi _2}(1)}}{{{\lambda _2}}}}& \cdots &{\frac{{{\varphi _i}(1)}}{{{\lambda _i}}}}\\ {\frac{{{\varphi _1}(2)}}{{{\lambda _1}}}}&{\frac{{{\varphi _2}(2)}}{{{\lambda _2}}}}& \cdots &{\frac{{{\varphi _i}(2)}}{{{\lambda _i}}}}\\ \vdots & \vdots &{}& \vdots \\ {\frac{{{\varphi _1}(n)}}{{{\lambda _1}}}}&{\frac{{{\varphi _2}(n)}}{{{\lambda _2}}}}& \cdots &{\frac{{{\varphi _i}(n)}}{{{\lambda _i}}}} \end{array}} \right] $

(6)

通过定义任意两点的内积, 可以得到两点间的双调和距离, 如下:

$ B(x,y) = \mathit{\boldsymbol{BS}}(x) \cdot \mathit{\boldsymbol{BS}}(y) $

(7)

B(x, y)是2个无限维向量的内积, 则双调和距离可表示为:

$ {d^2}(x,y) = \int_0^\infty {{B^2}} (x,y){\rm{d}}t $

(8)

为便于计算, 将LB算子进行谱分解, 则双调和距离谱分解后可写为:

$ {d^2}(x,y) = \sum\limits_{i = 1}^{N \to \infty } {\frac{{{{\left( {\varphi \left( x \right) - \varphi \left( y \right)} \right)}^2}}}{{\lambda _i^2}}} $

(9)

2 基于双调和距离的形状相似度

假定形状M上有n个顶点, 如果计算形状上任意2个顶点的双调和距离, 则该距离矩阵为一个n×n维的对称矩阵, 主对角线元素为0, 具体如下:

$ \mathit{\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} 0&{d\left( {1,2} \right)}& \cdots &{d\left( {1,n} \right)}\\ {d\left( {2,1} \right)}&0& \cdots &{d\left( {2,n} \right)}\\ \vdots & \vdots &{}& \vdots \\ {d\left( {n,1} \right)}&{d\left( {n,2} \right)}& \cdots &0 \end{array}} \right] $

(10)

如果应用D进行形状匹配, 其计算复杂度较高, 例如, 当n为10 000时, D为一个非常大的高维矩阵。因此, 在进行形状匹配时需要对D进行降维处理。

2.1 距离矩阵奇异值分解

奇异值分解(Singular Value Decomposition, SVD)是数值分析、统计中最重要的工具之一, 其主要思想是将矩阵分解为若干个矩阵相乘。SVD是提取特征的一种有效方法。本节对D进行奇异值分解并提取特征。

令${\mathit{\pmb{D}}} \in {\mathbb{R}^{m \times n}}$, 则存在正交矩阵${\mathit{\pmb{U}}} \in {\mathbb{R}^{m \times n}}$和${\mathit{\pmb{V}}} \in {\mathbb{R}^{m \times n}}$, 使得:

$ \mathit{\boldsymbol{D}} = \mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}} $

(11)

其中, U、V分别称为D的左、右奇异值矩阵, Σ为一个m×n维的矩阵。

a=diag(Σ)=(σ₁, σ₂, …, σ_n)为矩阵D的奇异值(σ₁≥σ₂≥…≥σ_r>0), σ_i为单奇异值。本文将取降序排列的非零奇异值记为:

$ \mathit{\boldsymbol{a'}} = \left( {{\sigma _1},{\sigma _2}, \cdots ,{\sigma _i}} \right) $

(12)

在实践过程中发现, σ₁、σ₂的值远大于其他单奇异值, 因此将矩阵D的奇异值进行归一化, 如下:

$ {{\mathit{\boldsymbol{a'}}}_{{\rm{norm}}}} = \frac{{\mathit{\boldsymbol{a'}} - \mu }}{v} $

(13)

其中, μ为a′的均值, υ为方差。为了增大不同形状奇异值的差异, 将a′_norm进行指数映射, 从而定义一种新的形状描述符BS_svd:

$ \mathit{\boldsymbol{BS}}\_\mathit{\boldsymbol{svd}}\left( M \right) = {{\rm{e}}^{{{\mathit{\boldsymbol{a'}}}_{{\rm{norm}}}}}} $

(14)

通过定义该描述符, 把形状从复杂的高维矩阵空间映射到简单的一维向量空间, 降低计算复杂度, 提升形状匹配的效率。

2.2 相似度计算

本节基于2.1节定义的形状描述符来计算形状相似度。由于该形状描述符是一个简单的一维向量, 因此直接选取向量的相似性度量作为形状相似度。在所有向量相似度量中, 余弦夹角可以有效地规避个体之间的差异而更注重维度之间的差异, 且能将相似度的范围收敛至[-1, 1]之间, 因此本文选择余弦夹角来进行相似度计算。余弦相似度, 又称为余弦相似性, 其通过计算2个向量的夹角余弦值来衡量个体之间的差异度。余弦值越接近1, 表明2个向量的方向越接近; 余弦值为-1, 表明2个向量的方向完全相反。本文定义的形状相似度计算如下:

$ d\left( {M,N} \right) = \frac{{\sum\limits_1^r {\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M) \times \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} \right)} }}{{\sum\limits_1^r {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M)} \right)}^2}} \times \sum\limits_1^r {{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} \right)}^2}} }} $

(15)

d(M, N)满足以下特性:

1) 非负性, 即d(M, N)≥0, 距离是一个非负的数值。

2) 同一性, 即d(M, M)=0, 形状到自身的距离为0。

3) 对称性, 即d(M, N)=d(N, M), 距离是一个对称函数。

4) 三角不等性, 即d(M, Q)≤d(M, N)+d(N, Q)。

同时, 本文列出其他5种常用的向量相似度计算方法, 并与余弦相似度进行对比。

1) 曼哈顿距离(L1 Distance):

$ d(M,N) = \sum\limits_1^r {\left| {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{s\nu d}}(M) - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{s\nu d}}(N)} \right|} $

(16)

2) 欧氏距离(L2 Distance):

$ \begin{array}{*{20}{l}} {d(M,N) = }\\ {\sqrt {(\mathit{\boldsymbol{BS\_svd}}(M) - \mathit{\boldsymbol{BS\_svd}}(N)) \times {{(\mathit{\boldsymbol{BS\_svd}}(M) - \mathit{\boldsymbol{BS\_svd}}(N))}^{\rm{T}}}} } \end{array} $

(17)

3) 推土机距离(Earth Mover’s Distance):

$ d(M,N) = \sum\limits_{i = 1}^m {\sum\limits_{i = 1}^n {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M) \cdot \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} } $

(18)

4) 相关系数(Correlation Coefficient):

$ d(M,N) = \frac{{ COV \left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M),\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} \right)}}{{\sqrt {D\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M)} \right)} \cdot \sqrt {D\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} \right)} }} $

(19)

5) 切比雪夫距离(Chebyshev Distance):

$ d(M,N) = \mathop {\max }\limits_i \left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(M) - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{S}}_ - }\mathit{\boldsymbol{svd}}(N)} \right) $

(20)

2.3 形状相似度算法

本文形状相似度算法的伪代码如算法1所示。

算法1  形状相似度算法

输入  M, N

输出  d(M, N)

1.for i=1:2

2.compute cotan Laplacian M, N

3.compute eigenvectors and eigenvalues of M, N

4.for j=1:n×n

5.compute biharmonic distance matrix DM

6.end for

7.for k=1:m×m

8.compute biharmonic distance matrix DN

9.end for

10.end for

11.DM=diag(DM), DN=diag(DN)

12.DM=(DM-mean(DM))./std(DM)

13.DN=(DN-mean(DN))./std(DN)

14.aM=exp(DM)

15.aN=exp(DN)

16.d(M, N)=cosine(aM, aN)

形状相似度算法的具体步骤如下:

步骤1  输入一对形状M和N, 读取形状的顶点和面片信息。

步骤2  对形状M和N的LB算子进行谱分解。

步骤3  计算形状M和N的双调和距离矩阵。

步骤4  对双调和距离矩阵进行奇异值分解, 提取特征值, 并将特征值进行归一化处理和指数映射。

步骤5  计算形状M和N的余弦相似度。

3 实验结果与分析

本文在TOSCA2010数据库上进行实验, 该数据库共包括80个形状, 包括dog、david、centaur和horse等, 顶点数为50 000左右, 该数据库适用于等距形状分析, 部分形状如图 1所示。

图 2给出dog、horse、david和centaur单个的形状描述符图及其类别的形状描述符图。从图 2可以发现, 在形状发生等距变化时, 本文定义的形状描述符具有鲁棒性, 能够有效地区分不同类的形状。同时, 对于双调和距离矩阵的单个奇异值呈降序排列时, 数值变化相对较小。为了提高运算速度, 本文中取r为20。

本文使用热力图对本文形状相似度算法与其他相似度算法进行对比。热力图是一种以颜色代替数值变化的效果图, 其清晰明了, 更易看出规律。图 3为本文算法分别使用6种距离时的热力图。

从图 3可以看出, 在使用余弦距离时, 算法的性能较优, 每一类形状内相似度的差较小, 不同类之间的相似度较大, 能有效地区分不同类的形状。

为了验证本文方法的性能, 将该算法与FMPS方法^[21]、SHED方法^[9]进行比较。表 1为3种方法进行4类形状匹配时的耗时情况对比。由表 1可以看出, 与FMPS方法、SHED方法相比, 本文方法的耗时情况介于两者之间。FMPS方法由于需要大量的迭代计算, 因此耗时较长。表 2为使用3种方法进行不同类形状匹配的匹配度对比。由表 2可知, 本文方法、FMPS方法、SHED方法的总体匹配度分别为95.82%、95.82%和75.00%。其原因在于SHED方法不具有等距不变性, 不能很好地区分发生等距变化的形状, 尤其是对于dog和horse的区分度较差, 导致总体匹配度较低。而本文方法与FMPS方法都具有等距不变性, 因此适用于3D非刚性形状。结合时间耗费和形状匹配度对比结果可以看出, 本文方法能够很好地兼顾时间耗费度和形状匹配度。

图 4和图 5分别为使用FMPS方法和SHED方法得到的形状距离矩阵热力图。对比图 3和图 4可知, 本文算法使用余弦距离的结果与FMPS算法的结果较为接近, 而SHED算法对3D非刚性形状的整体区分度较低, 性能较差。

4 结束语

本文在双调和距离的基础上, 提出一种三维形状相似度计算方法, 定义了一种形状双调和距离矩阵奇异值分解的形状描述符。该描述符将复杂的高维矩阵映射为简单的低维向量, 把形状匹配的问题转换为向量距离求解的问题, 并使用余弦距离作为形状相似度。实验结果表明, 与FMPS方法、SHED方法相比, 该方法具有较低的时间耗费和较高的形状匹配度。下一步将引入形状分割理论对部分形状相似度计算进行改进。