在实体动画、模型检索等有关三维形状分析和几何处理的许多应用中，模型匹配的准确性对算法起着至关重要的作用. 描述给定模型显著特征的主流方法之一是通过定义模型上各点的特征描述符，即编码该点邻域的多尺度形状信息和几何特征，实现模型整体几何特性的描述. 在众多建立点描述符的方法中，基于扩散几何的方法因具有许多显著的优点成为主流方法之一^[1]. 其中建立点描述符的主要代表工作有热核签名（heat kernel signature，HKS）^[2]和波核签名（wave kernel signature，WKS）^[3]. 从信号处理角度而言，这2种描述符具有统一的数学表达形式，均可以视为一组滤波器作用于拉普拉斯算子的特征值与特征向量. HKS源于一组低通滤波器的使用，主要分析模型的轮廓信息，因此无法实现模型之间点-点匹配时点的精确定位. WKS使用一组带通滤波器，能够获取模型更多的细节信息，实现了比HKS更准确的定位，但由于缺乏低通滤波器捕捉与模型轮廓相关的低频信息，匹配时容易产生大量的匹配噪声. Boscaini等^[4]利用多个几何向量的流形窗口傅里叶变换（WFT）^[5]系数，建立多尺度的点描述符. Melzi等^[6]通过模拟穿过网格点对之间所有路径的离散时间演化序列，创建多尺度的离散时间演化过程描述符（DEP）. 此外，现阶段有许多工作根据不同的应用，利用机器学习建立三维模型点描述符^[7-8].

Hammond等^[9]类比于经典小波在频率域里的性质，基于谱图理论在图域上建立谱图小波变换（spectral graph wavelet transform，SGWT）. 该小波继承了欧氏空间经典小波的许多性质，利用一组低通和高通滤波器实现图上信号的多频带分析. 该小波变换已广泛应用于多个三维模型处理领域，包括分片和校正、模型的稀疏表示和光滑、模型检索以及医学处理等^[10-15].

本文利用谱图小波，创建适用于三维模型点-点匹配的点描述符. 为了充分利用模型的几何特征，将谱图小波构建中使用的图拉普拉斯算子改为流形上的拉普拉斯-贝尔特米算子，利用模型各点处脉冲函数的多尺度小波系数和尺度函数系数描述各点邻域的局部和全局结构信息，建立各点处的描述符. 因为该描述符所使用的滤波器结合HKS和WKS的优势，能够实现对信号的精确重构，这使得所提出的描述符具有很强的分辨力，且形式紧凑、算法鲁棒. 此外，因为算法只需对拉普拉斯矩阵的少量特征值和特征向量进行代数运算，计算效率非常高.

设给定的三维模型X为嵌入 ${{\mathbb{{\bf{R}}}}^3}$ 空间的紧连续二维流形（可能具有边界），将该流形上的平方可积函数空间表示为 ${L^2}(X){\rm{ = }}\left\{ {f:X \to \mathbb{{\bf{R}}},\int_X {f{{(x)}^2}{\rm d}x < \infty } } \right\}$，其中dx为给定黎曼度量下的面积微元. X上的拉普拉斯-贝尔特米算子定义为

${\Delta _X}f(x) = - {\operatorname{div} _X}({\nabla _X}f(x))，.$

式中： ${\nabla _X}f$ 和 ${\operatorname{div} _X}f$ 分别为函数 $f(x) \in {L^2}(X)$ 内蕴的梯度和散度.

在实际应用中，三维模型X通常表示为三角网格 $X = \{ V,E,F\} $，其中V为具有N个顶点的网格顶点集，（i，j）表示连接顶点i和顶点j的边，F为三角面片集合. 将第i个网格顶点的一环邻域记为 $N(i) = $ $\left\{ {j \in V|(i,j) \in E} \right\}$. 向量 ${f}{\rm{ = }}\left[ {{f_i}} \right] \in {\mathbb{{\bf{R}}}^N}$ 表示定义在网格X上的离散信号（函数），其中第i个分量 ${f_i} = f(i)$ 为信号f在网格X第i个顶点处的值.

采用最普遍的Meyer^[16]方法计算拉普拉斯-贝尔特米算子在三角网格X上的离散形式，具体形式为

${\Delta _X}f(i) = \frac{1}{{{a_i}}}\displaystyle\sum\limits_{j \in N(i)} {{\omega _{ij}}[f(i) - } f(j)].$

式中： ${a_i}$ 为顶点i的Voronoi单元面积；ω_ij为权重， ${\omega _{ij}} = \left({{\cot {\alpha _{ij}} + \cot {\beta _{ij}}}}\right)/{2}$； ${\alpha _{ij}}$ 与 ${\beta _{ij}}$ 分别为边（i，j）所对着的两个角，如图1所示.

令面积矩阵A=diag [a_i]，权重矩阵M=[m_ij]，

${m_{ij}} = \left\{ \begin{array}{l} \displaystyle\sum\limits_{k \in N(i)} {{\omega _{ik}},\quad i{\rm{ = }}j}; \\ - {\omega _{ij}},\quad \quad i \ne j{\text{且}}\,j \in N(i);\\ 0,\quad \quad \quad {\text{其他}}, \end{array} \right.$

则拉普拉斯-贝尔特米算子可以表示为 $N \times N$ 的矩阵L=A⁻¹M，称为流形上的拉普拉斯矩阵.

拉普拉斯矩阵L的特征值和特征向量，满足 ${L}{{\phi }_k} =$ $ {\lambda _k}{{\phi }_k},\,\;k = 0,1, \cdots ,N - 1.$ 其中 ${{\phi }_k} = \left[{\phi _k}(i)\right] \in {\mathbb{{\bf{R}}}^N}$. 事实上，L的特征值和特征向量具有调和性质，可以视为流形上信号空间 ${L^2}(X)$ 的频率和傅里叶基^[17]. 在实际计算中，求解L=A⁻¹M的特征系统可以转换为求解广义特征值问题 ${M}{{\phi }_k} = \lambda {A}{{\phi }_k}$. 由于L具有半正定性，特征值 $\displaystyle\{ {\lambda _k}\} _{k = 0}^{N - 1}$ 均为正实数. 将其按升序排列，可得 $0{\rm{ = }}{\lambda _0} < {\lambda _1} \leqslant \cdots \leqslant {\lambda _{N - 1}}$. 对应的特征向量具有关于矩阵A内积的正交性，即

${\left\langle {{{\phi }_i},{{\phi }_j}} \right\rangle _{{A}}} = {{\phi }_i}^{\rm{T}}{{A}}{{\phi }_j} = \left\{ \begin{gathered} 0,\quad i \ne j; \\ 1,\quad i = j. \\ \end{gathered} \right.$

任给信号 ${f}{\rm{ = [}}f{\rm{(}}i{\rm{)]}} \in {\mathbb{{\bf{R}}}^N}$，傅里叶变换系数 $\hat f(k)$ 为 ${f}$ 与特征向量 ${{\phi }_k}$ 的内积，即

(1) $\hat f(k) = {\left\langle {{f},{{\phi }_k}} \right\rangle _{{A}}}{\rm{ = }}{{f}^{\rm{T}}}{A}{{\phi }_k} = \displaystyle\sum\limits_{i = 1}^N {{a_i}f(i)} {\phi _k}(i).$

事实上，由于信号f的能量主要集中在低频部分，实际应用中不需要求解出全部的特征值与特征向量，选择前K（ $K \ll N$）个特征值与特征向量即可满足要求. 类似于经典信号处理，流形上存在反傅里叶变换，即

$f(i) \approx \displaystyle\sum\limits_{k = 0}^K {{{\left\langle {{f},{{\phi }_k}} \right\rangle }_{{A}}}{\phi _k}(i)} ,\quad $

可以高度逼近原始信号.

将上述结果写成矩阵形式，令N×K矩阵 ${\varPhi }{\rm{ = }}$ $\left[ {{{\phi }_1},{{\phi }_2}, \cdots ,{{\phi }_K}} \right]$，向量 $\hat {{f}} = \left[ {\hat f(k)} \right] \in {\mathbb{{\bf{R}}}^N}$ 表示傅里叶系数，则有 $\hat {{f}} = {{\varPhi }^{\rm T}}{Af}$ 和 ${f} = {\varPhi }\hat {{f}}$.

经典小波由欧氏空间内的母小波函数伸缩平移获得，对于定义在三维曲面上的函数，没有直观自然的方式定义这些函数的伸缩和平移变换. 为了克服该问题，Hammond等^[9]类比经典小波在频率域里的滤波特性，基于谱图理论构建图域上的小波. 该小波具有很多与经典小波类似的性质，是分析图域空间信号的强有力工具. 为了充分利用三维模型的几何特征，利用拉普拉斯-贝尔特米算子的特征系统，基于谱图小波的框架构建更适合分析流形上信号的小波，称为谱图小波.

谱图小波函数和小波系数 　给定小波核函数 $g:{\mathbb{{\bf{R}}}^ + } \to {{ \mathbb{{\bf{R}}}}^ + }$，其作用类似于带通滤波器，满足条件 $g(0) = 0,\mathop {{\kern 1pt} \lim }_{x \to \infty } g(x) = 0$. 与经典小波类似，具有不同尺度的谱图小波可以视为脉冲函数经过不同频带带通滤波器的滤波结果.

取拉普拉斯矩阵L的前K个特征值和特征向量，定义尺度为 ${t_l}$，位于三角网格X第i个顶点处的离散小波函数为

${\psi _{{t_l},i}}(m) = \displaystyle\sum\limits_{k = 0}^{K - 1} {g({t_l}{\lambda _k}){{\hat \delta }_i}(k){\phi _k}(m)}. $

式中： ${\hat \delta _i}(k)$ 为位于第i个顶点的单位脉冲函数 ${{\delta }_i}$ 的傅里叶系数. 根据式（1），有

${\hat \delta _i}(k) = {\left\langle {{{\delta }_i},{{\phi }_k}} \right\rangle _{{A}}} = {{\delta }_i}^{\rm{T}}{A}{{\phi }_k} = {a_i}{\phi _k}(i),$

进一步可得

${\psi _{{t_l},i}}(m) = \displaystyle\sum\limits_{k = 0}^{K - 1} {{a_i}g({t_l}{\lambda _k}){\phi _k}(i){\phi _k}(m)} $.

对于任意信号f，小波系数为f与小波函数 ${{\psi }_{{t_l},i}}$ 的内积，即

${W_{{f}}}({t_l},i){\rm{ = }}{\left\langle {{f},{{\psi }_{{t_l},i}}} \right\rangle _{{A}}}{\rm{ = }}\sum\limits_{k = 0}^{K - 1} {{a_i}g({t_l}{\lambda _k})\hat f(k){\phi _k}(i)}. $

尺度函数和尺度函数系数 　为了分析信号的低频部分，构造相应的尺度函数. 选取尺度函数核函数 $h:{\mathbb{{\bf{R}}}^ + } \to {\mathbb{{\bf{R}}}^ + }$，该函数满足条件 $h(0) > 0$ 和 $h(x) $ $ \to 0$ $\left(x \to \infty \right) $. 定义位于网格第i个顶点的尺度函数 ${{\varphi }_i}$ 为

${\varphi _i}(m) = \displaystyle\sum\limits_{k = 0}^{K - 1} {{a_i}h({\lambda _k})} {\phi _k}(i){\phi _k}(m),$

相应的尺度函数系数为

${S_{{f}}}(i){\rm{ = }}{\left\langle {{f},{{\varphi }_i}} \right\rangle _{{A}}} = \sum\limits_{k = 0}^{K - 1} {{a_i}h({\lambda _k})\hat f(k){\phi _k}(i)}. $

谱图小波的尺度函数是为了分析小波函数未能分析的信号低频部分，不需要与小波函数形成类似经典小波里的双尺度关系.

为了方便计算，将网格X上各顶点处的小波系数和尺度函数系数整理成矩阵形式. 令N×1维向量 ${W}_{{f}}^l{\rm{ = }}{W_{{f}}}({t_l}, \cdot )$， ${{S}_{{f}}}{\rm{ = }}({S_{{f}}}(i))$，对角矩阵 ${\varLambda }{\rm{ = }}{\rm diag}\;\left[{\lambda _1},\right.$ $\left.{\lambda _2}, \cdots ,{\lambda _K}\right]$，则有

${W}_{{f}}^l{\rm{ = = }}({A}{\varPhi })g({t_l}{\varLambda })({{\varPhi }^{\rm{T}}}{A}){f},$

${{S}_{{f}}}{\rm{ = }}({A}{\varPhi })h({\varLambda })({\varPhi }{A}){f}.$

核函数的选取 　在谱图小波工具箱（SGW_TOOLBOX）^[9]中提供了多种小波核函数，可以根据应用灵活的选择. 为了使各频率的影响相同，选取墨西哥帽小波核函数，在图4（c）中展示了具有不同尺度 ${t_l}$ 的墨西哥帽小波核函数 $g(x) = x\exp\; ( - x)$ 及其尺度函数核函数

$h(x) = 1.2\exp \left[ { - {{\left( {\displaystyle\frac{x}{{0.4{\lambda _{\min }}}}} \right)}^4} - 1} \right].$

离散尺度的选取 　离散小波尺度 ${t_j}$ 由拉普拉斯矩阵 ${L}$ 的特征值上界 ${\lambda _{\max }}$ 及尺度总数J确定. 最小和最大尺度分别为 ${t_J}{\rm{ = }}1/{\lambda _{\max }}$ 和 ${t_1}{\rm{ = 2}}/{\lambda _{\min }}$，其中 ${\lambda _{\min }} =$ $ {\lambda _{\max }}/u,$ 参数 $u > 0$ 可以由用户根据应用自由选择，其值确定了小波能分析的频率下界. 余下的尺度 ${t_l}\;(l = 2, \cdots ,J - 1)$ 在区间 $[{t_J},{t_1}]$ 上等对数分布. 图2展示了位于模型同一点处的尺度函数和具有不同尺度的小波函数，其中J=7. 随着尺度 ${t_l}\;$ 的变小（即l增大时），谱图小波在模型上的支撑域逐渐变小，这表示谱图小波具备自适应分析高频信息的能力. 事实上，当选定一组尺度 $\left\{ {{t_l}} \right\}_{l = 1}^J$ 之后，函数集 $F = \left\{ {{{\varphi }_i}} \right\}_{i = 1}^N$ $ \cup \left\{ {{{\psi }_{{t_l},i}}} \right\}_{l = 1,i = 1}^{J,N}$ 形成一组框架，边界为 $p = \mathop {\min }_{\lambda \in [0,{\lambda _{\max }}]} G(\lambda )$ 和 $q = \mathop {\max }_{\lambda \in [0,{\lambda _{\max }}]} G(\lambda )$，其中 $G(\lambda ){\rm{ = }}h{(\lambda )^2} + {\displaystyle\sum \nolimits _ {{l}}} {g{{({t_l}\lambda )}^2}} $. 只有当p和q远离零时，才能保证信号f被该组框架稳定重构.

谱图小波系数和尺度函数系数全面分析了信号的各频带信息. 若给定模型第i个顶点处的单位脉冲函数 ${{\delta }_i}$，则该函数的多尺度小波系数和尺度函数系数能够反映模型在该点周围的多尺度结构信息. 将这些多尺度小波系数和尺度函数系数组合成一个向量，可以生成一种新的点描述符，称为谱图小波描述符（spectral graph wavelet descriptor，SGWD），即

${\rm SGWD}(i) = \left( {{S_{{{\delta }_i}}}(i),{W_{{{\delta }_i}}}({t_1},i),{W_{{{\delta }_i}}}({t_2},i), \cdots ,{W_{{{\delta }_i}}}({t_J},i)} \right).$

事实上，根据脉冲函数的傅里叶系数，可得SGWD各分量：

(2) ${W_{{{\delta }_i}}}({t_j},i){\rm{ = }}{\left\langle {{{\delta }_i},{{\psi }_{{t_j},i}}} \right\rangle _{{A}}}{\rm{ = }}\sum\limits_{k = 0}^K {{a_i}^2g({t_j}{\lambda _k}){\phi }_k^2(i)} ,$

(3) ${S_{{{\delta }_i}}}(i){\rm{ = }}{\left\langle {{{\delta }_i},{{\varphi }_i}} \right\rangle _{{A}}}{\rm{ = }}\sum\limits_{k = 0}^K {{a_i}^2h({\lambda _k}){\phi }_k^2(i)} .$

从式（2）、（3）可以推知，SGWD具有如下优良的性质.

1）SGWD具有等距变形不变性. 拉普拉斯矩阵L的特征系统具有模型内蕴等距变形不变性，因SGWD的计算仅涉及L的特征值和特征向量，易知SGWD具有该特性. 图3展示了SGWD的这一性质. 图中，l为每个点SGWD的第l个分量索引. 图3从左至右分别为具有近似等距变形的2个人体模型，最右端图中的实线和虚线分别表示左、右人体模型上相同颜色点的SGWD. 这2种线的高度重合，证实了SGWD具有等距变形不变性.

2）SGWD结合了HKS和WKS的优点，具有更高的分辨能力. 由式（2）、（3）可知，HKS、WKS及SGWD这3种描述符的本质区别在于滤波器，即核函数选取的不同. HKS与WKS的核函数分别为 ${g_{\rm HKS}} = \exp\; ( - t\lambda ),$ ${g_{{\text{WKS}}}} = {C_t}\exp \;\left[ { - {{\left( {\log t - \log \lambda } \right)}^2}/{\sigma ^2}} \right].$如图4所示为具有不同尺度参数t的这3种核函数的图像以及相应的平方和函数 $G(\lambda )$. 图中，λ为特征值，g为相应核函数的值，黑色线表示 $G(\lambda )$，虚线分别表示 $G(\lambda )$ 的上、下界p和q. 从图4可以看出，谱图小波描述符同时使用了一个低通与一组带通滤波器提取模型的多尺度信息，综合了HKS和WKS的优势. 结果显示，只有谱图小波的p和q均远离零，这意味着使用有限的小波系数即可实现信号的稳定重构. 这些优势赋予了SGWD更好的分辨能力和局部定位能力.

3）形式紧凑，计算效率高. SGWD只需计算一组非常简单的脉冲函数的小波系数，无需事先计算更复杂的函数，且小波系数的计算仅涉及拉普拉斯矩阵小部分特征值与特征向量的代数运算，保证了SGWD的高计算效率. 因小波框架能够精确地重构原始信号，这意味着较低维度的多尺度系数可以全面包含模型的信息，使得SGWD结构非常紧凑.

4）鲁棒性. 谱图小波理论^[9]保证了SGWD的鲁棒性. 图6显示了SGWD在经历了多种变形的模型上的鲁棒性.

所有实验均在具有Intel（R）Core^TMi7-4790，CPU主频为3.60 GHz 和16.0 GB RAM 的PC机上开展，使用Matlab R2018a.

数据集：实验采用FAUST^[18]和CAESAR^[19]这2个公共人体模型数据集，模型间点与点之间的准确对应已知. 这2个数据集中的模型因为存在非等距变形，且不同的对象间存在较大的差异，使得模型间点-点匹配非常具有挑战性. 其中FAUST包含10个不同的对象，且每个对象具有10个不同的姿势，共计100个模型，每个模型具有大约7 000个顶点. CAESAR 是最大的商业化人体模型数据集，包含超过4 500个对象且具有同一站立姿势的模型.选取fitted-mesh子数据集中的20个模型进行测试，每个模型大约具有6 000个顶点.

参数设置：在所有实验中，拉普拉斯特征值和特征向量取前K=200个. 理论上J越大，SGWD分辨能力越强，当J到达一定阶段后，描述符分辨力提升较缓但计算效率下降迅速. 如图5所示，计算模型各点处不同J的SGWD，并可视化其与参考点（最左端胳膊处的点）描述符的正则化欧氏距离. 结果显示，J=8时的描述符分辨力比J=4时更强，但和J=16时无明显差异，因此选取J=10.

为了公平比较，以下实验均采用欧氏距离度量描述符空间中点的距离. 将与参考点的描述符具有最短距离的点作为各算法得到的参考点的匹配点. 将准确对应点的内蕴对称点视为正确的匹配点.

评价标准：对于所有的数据集，使用3种标准评估算子的性能，包括累积匹配特征（cumulative match characteristic，CMC），受试者工作特征（receiver operator characteristic，ROC）和匹配质量特征（correspondence quality characteristic ，CQC）. CMC估计在描述符空间的k个最近邻中出现正确对应点的概率. ROC测量在描述符空间里positives and negatives 点对之间距离小于不同阈值的百分比. CQC计算描述符空间里最近邻与准确匹配点的测地距离小于不同r的点的百分比.

实验结果：图6中，最左端的参考模型来自FAUST，选择该模型膝盖处的点作为参考点（圆球标注）. 从左至右为经过多种变形的参考模型，依次为等距变形、非等距变形、光滑、高斯噪声、尖锐噪声以及部分缺失. 计算每个模型各顶点处与参考点SGWD的距离. 如图6所示为正则化后的距离，结果显示，SGWD在各种变形中保持非常高的准确性，证实该描述符具有很好的鲁棒性. 如图7所示为SGWD与HKS、WKS等几种现阶段具有代表性的描述符在数据库FAUST和CAESAR中的性能比较. 图中，C为匹配率，E为测地距离，T为真阳性率，F为假阳性率. 综合CMC、ROC和CQC这3种评价指标可知，SGWD在这2类数据库里都显示了比其他方法更好的匹配结果. DEP和WFT因为算法的要求，计算效率远不及SGWD. 表1给出3个主要模型不同描述符计算时间的统计. 可以看出，本文算法的计算效率优于WFT和DEP.

如图8、9所示为针对不同类别的模型，不同描述符性能的可视化比较结果. 图8选择FAUST中的8个人体模型进行测试. 最左端为参考模型和参考点（圆球标记），计算所有模型顶点处的SGWD，将其与参考点描述符的欧氏距离正则化和可视化. 结果显示，与其他方法相比，利用SGWD寻找到的与参考点距离最小的区域仅发生在准确匹配点的邻域内，且区域面积相对最小. 这表示SGWD具有最高的分辨力和定位能力. 图9选择2类动物模型，最左端为参考模型，右端均为参考模型经过同一等距变形的结果. 在参考模型上选择参考点（圆球标注），计算等距变形模型中所有点与参考点的距离并可视化. 结果显示，SGWD的分辨力和定位能力优于其他算法.

本文利用一类简单函数的谱图小波系数和尺度函数系数，建立流形上的描述符，用以描述模型各点处邻域的局部和全局特征. 该描述符综合了2类代表性描述符的优势，谱图小波的优良性质保证了提出的描述符具有十分强的分辨力和定位能力，形式紧凑，计算效率高. 将谱图小波描述符应用至模型点-点匹配时，获得了比现阶段同类方法更好的性能.

在下一阶段，考虑将SGWD应用至其他几何处理和形状分析领域，以期进一步挖掘SGWD的优势.