图像特征点是用来表示图像中具有特殊含义或代表性的节点，如人体的关节点、物体梯度的变化点、目标的感兴趣点等。图像特征点匹配可以对两个图像中的特征点建立一一对应关系，从而为识别图像中物体的位置、姿态、甚至动作打下基础，也使其在对象跟踪、全息影像创作、图像分类、图像对齐等领域处于关键位置^[1-2]。

图论作为数学的一个分支，在20世纪70年代被提出应用到图像匹配问题中，由于需要巨大的运算量，而当时的客观条件无法满足需求，在经历了最初对这个巧妙理论结构的热情后，这方面的研究被搁置了很长一段时间。随着计算机技术的不断发展，近几年，图论理论在图像匹配中的应用再次被提了出来，并且有了不小的进展。

2005年，Berg等人^[3]提出了二维图像匹配算法，将图像匹配问题设定为一个整体二次规划问题。其中的价值函数是基于对应几何模糊点描述符的相似性，例如相应特征点的几何失真等。该算法使基于图论的图像匹配得以实现，并且对异常值也有一定的掌控能力。2008年，Mateus等人^[4]提出了一种解决三维图像匹配的方法，使用权重图描述由像素组成的清晰的形状，图像匹配就可以被归纳为最大子图同构问题。将稠密点通过无监督聚类和EM算法，在正交变换下转换为图，通过用直方图建立的特征函数为图附上权重。在此基础上，通过选择拉普拉斯矩阵的特征函数的最佳子集来寻找两个相对应的$ K $维点集的最好匹配。这样使用谱图理论就将这些图映射在了低维空间中。实验结果表明，该方法可以处理拓扑结构的变化、形状变化和不同的采样率。2011年，Jiang等人^[5]将图像匹配应用于对象跟踪，提出在一个约束线性空间内，同时匹配特征点和估计全局几何变换的线性方程，解决缩放、旋转和变形的视觉形态匹配问题。该方法可以解决包含非常大量的候选特征点的大尺度问题，并且在应对弱特征和混乱的问题时更加稳健。2014年，韩枫^[6]使用图论实现了三维耳廓形状的特征匹配。该工作的目标是利用耳廓形状的唯一性进行个人身份识别，这需要将激光扫描出来的耳廓的三维点云与数据库中的耳廓数据进行匹配，从而进一步计算相似度。该研究使用IsoRank算法对耳廓进行匹配，通过寻找两个需要匹配的耳廓的关键点图的最大整体匹配来找到两个图形节点之间的对应关系，并且通过路径跟踪算法解决非凸问题来提高匹配精度。2015年，Zhou等人^[7]提出了分解图匹配算法(FGM)，通过将1个大的相似矩阵分解为6个小的矩阵来减小运算的复杂度，大大加快了运算速度。并且由于这个算法是一个可以兼容多种策略的框架，所以可以集众家之长来提高准确率和适用范围。并且将路径跟踪算法与迭代最近点法加入其中，形成可变形的图匹配算法(DGM)，初步解决了非刚体图像的特征点匹配问题，对有形变的图像匹配有较高的准确率。同年，Dou等人^[9]用三维图匹配的方法解决了蛋白质三维结构中螺旋线在非刚性变换后的配准问题。在匹配的过程中，将三维螺旋线分解成多段，并构建了相似性函数。匹配的原则是分解后的线段长度和刚性变换的程度。虽然在实验中取得了较好的结果，但该方法面向特定领域，不具备普适性。2016年，Wang等人^[8]提出了Symmetry-aware图匹配策略，旨在解决对称图像的匹配问题，并且达到了一定的效果。该策略提出了一个对称约束，将这个对称约束加入其他图匹配算法中，可以加速具有对称结构的图像的匹配，并且显著提高精度。

综上所述，目前的图匹配算法大多面向二维图像，面向三维图像时，缺少既准确又快速的匹配算法。本文拟将分解图匹配算法(FGM)从现在仅适用于二维图像，推广到适用于三维图像。在众多图匹配算法中选择FGM，是因为在二维图像特征点匹配过程中，FGM表现出了更好的准确率和更快的速度。

面向三维图像的图匹配算法总体流程图如图 1所示，首先将需要匹配的两个三维图像的特征点作为图的节点集(vertices set)。然后通过Delaunay三角剖分算法，将三维特征点相连，得到的边作为图的边集(edges set)，从而建立有向图。令$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $分别表示两个图，通过$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的有向边的起始点，得到边的起点特征矩阵$ {\mathit{\boldsymbol{G}}_1} $和$ {\mathit{\boldsymbol{G}}_2} $；同理，通过$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的有向边的终点，得到边的终点特征矩阵$ {\mathit{\boldsymbol{H}}_1} $和$ {\mathit{\boldsymbol{H}}_2} $。再利用得到的有向图计算需要的三维节点特征和三维边特征，并将各种特征进行归一化处理和整合，从而组成节点特征相似矩阵$ {\mathit{\boldsymbol{K}}_p} $和边特征相似矩阵$ {\mathit{\boldsymbol{K}}_q} $。随后，根据这6个矩阵$ {\mathit{\boldsymbol{G}}_1} $、$ {\mathit{\boldsymbol{G}}_2} $、$ {\mathit{\boldsymbol{H}}_1} $、$ {\mathit{\boldsymbol{H}}_2} $、$ {\mathit{\boldsymbol{K}}_p} $、$ {\mathit{\boldsymbol{K}}_q} $，计算得到节点的对应关系矩阵$ \mathit{\boldsymbol{X}} $，也就是三维图像特征点的匹配关系矩阵。

令需要进行特征点匹配的两个三维图像分别为$ {\mathit{\boldsymbol{I}}_1} $和$ {\mathit{\boldsymbol{I}}_2} $，其特征点集合分别为$ {\mathit{\boldsymbol{N}}_1} $和$ {\mathit{\boldsymbol{N}}_2} $，且$ {\mathit{\boldsymbol{N}}_1} $和$ {\mathit{\boldsymbol{N}}_2} \in {\mathit{\boldsymbol{R}}^3} $。根据图论可知，一个图(graph)的构造，需要节点集和边集。本文将特征点集合$ {\mathit{\boldsymbol{N}}_1} $和$ {\mathit{\boldsymbol{N}}_2} $直接作为图的节点集$ {\mathit{\boldsymbol{V}}_1} $和$ {\mathit{\boldsymbol{V}}_2} $，通过Delaunay三角剖分算法将相邻的节点相连，从而构成边集$ {\mathit{\boldsymbol{E}}_1} $和$ {\mathit{\boldsymbol{E}}_2} $。需要指出的是，这些有向边都是双向的。根据节点集和边集就可以构造与三维图像$ {\mathit{\boldsymbol{I}}_1} $和$ {\mathit{\boldsymbol{I}}_2} $对应的有向图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $，如图 2所示。

令$ {n_1} $和$ {n_2} $分别为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的节点数量，$ {m_1} $和$ {m_2} $分别为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的边数量，定义矩阵$ {\mathit{\boldsymbol{G}}_1} $和$ {\mathit{\boldsymbol{G}}_2} $分别用来描述两个图的边的起始节点，矩阵的行数等于顶点数量，列数等于边数量，每一列表示图中的一条有向边$ e $，并且每一列中只有一个元素的值为1，其他都为0，值为1的元素对应于图中的一个节点$ v $，表示有向边$ e $的起点是$ v $。同理，定义矩阵$ {\mathit{\boldsymbol{H}}_1} $和$ {\mathit{\boldsymbol{H}}_2} $分别用来描述图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的边的终止节点，矩阵的行数等于顶点数，列数等于边数，每一列表示图中的一条有向边$ e $，并且每一列中只有一个元素的值为1，其他都为0，值为1的元素对应于图中的一个节点$ v $，表示有向边$ e $的终点是$ v $。因此，将$ {\mathit{\boldsymbol{G}}} $矩阵与$ {\mathit{\boldsymbol{H}}} $矩阵联合起来就可以表示图中每条边的起始和终止节点。图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $对应的起点特征矩阵$ {\mathit{\boldsymbol{G}}_1} $、$ {\mathit{\boldsymbol{G}}_2} $和终点特征矩阵$ {\mathit{\boldsymbol{H}}_1} $、$ {\mathit{\boldsymbol{H}}_2} $为

$ {\mathit{\boldsymbol{G}}_1} = \left[ {\begin{array}{*{20}{l}} 1&0&0&0&0&0\\ 0&1&0&1&0&0\\ 0&0&1&0&1&0\\ 0&0&0&0&0&1 \end{array}} \right] $

$ {\mathit{\boldsymbol{H}}_1} = \left[ {\begin{array}{*{20}{c}} 0&0&0&1&0&0\\ 1&0&0&0&1&0\\ 0&1&0&0&0&1\\ 0&0&1&0&0&0 \end{array}} \right] $

$ {\mathit{\boldsymbol{G}}_2} = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&1&0\\ 0&0&0&1 \end{array}} \right], \;\;\, {\mathit{\boldsymbol{H}}_2} = \left[ {\begin{array}{*{20}{c}} 0&0&1&0\\ 1&0&0&1\\ 0&1&0&0 \end{array}} \right] $

(1)

令$ {\mathit{\boldsymbol{P}}_1} $与$ {\mathit{\boldsymbol{P}}_2} $分别为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的节点特征矢量，$ {\mathit{\boldsymbol{Q}}_1} $与$ {\mathit{\boldsymbol{Q}}_2} $为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的边特征，则

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{P}}_1} = \left[ {p_1^1, \cdots , p_{{n_1}}^1} \right] \in {{\bf{R}}^{{d_p} \times {n_1}}}}\\ {{\mathit{\boldsymbol{P}}_2} = \left[ {p_1^2, \cdots , p_{{n_2}}^2} \right] \in {{\bf{R}}^{{d_p} \times {n_2}}}}\\ {{\mathit{\boldsymbol{Q}}_1} = \left[ {q_1^1, \cdots , q_{{m_1}}^1} \right] \in {{\bf{R}}^{{d_q} \times {m_1}}}}\\ {{\mathit{\boldsymbol{Q}}_2} = \left[ {q_1^1, \cdots , q_{{m_2}}^1} \right] \in {{\bf{R}}^{{d_q} \times {m_2}}}} \end{array}} \right. $

(2)

式中，$ {{n_1}} $和$ {{n_2}} $分别为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的节点数量，$ {{m_1}} $和$ {{m_2}} $分别为图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的边数量，$ {d_p} $和$ {d_q} $分别为节点特征与边特征的维度，也可以理解为特征的数量。

本文的节点特征选取了度和离心率来表示^[10]。度是指某一节点连接的边的数量；离心率$ \varepsilon $指某一节点与任意节点的最小跳数的最大值，即

$ {\varepsilon _i} = \mathop {\max }\limits_{{n_j} \in \mathit{\boldsymbol{N}}} \left( {{H_{{n_i} \to {n_j}}}} \right) $

(3)

$ H_{A \rightarrow B}=\min\limits _{k \in\{1, \cdots, N-1\}}\left(P_{A \rightarrow B}(k)\right) $

(4)

本文的边节点特征选取了边的长度$ l $和边与$ XOY $平面的夹角$ \theta $来表示^[10]，计算公式为

$ l=\sqrt{\left(x_{1}-x_{2}\right)^{2}+\left(y_{1}-y_{2}\right)^{2}+\left(z_{1}-z_{2}\right)^{2}} $

(5)

$ \begin{array}{c}{\theta=} \\ {\text { arcos } \frac{\sqrt{\left(x_{1}-x_{2}\right)^{2}+\left(y_{1}-y_{2}\right)^{2}}}{\sqrt{\left(x_{1}-x_{2}\right)^{2}+\left(y_{1}-y_{2}\right)^{2}+\left(z_{1}-z_{2}\right)^{2}}}}\end{array} $

(6)

式中，$\left(x_{1}, y_{1}, z_{1}\right)$表示边的起点三维坐标，$\left(x_{2}, y_{2}, z_{2}\right)$表示边的终点的三维坐标。

通过$ {\mathit{\boldsymbol{P}}_1} $、$ {\mathit{\boldsymbol{P}}_2} $可以计算图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的节点特征相似矩阵$ {\mathit{\boldsymbol{K}}_p} $，其行数等于图$ {\mathit{\boldsymbol{A}}_1} $的顶点数，列数等于图$ {\mathit{\boldsymbol{A}}_2} $的顶点数，计算方法是通过$ {\mathit{\boldsymbol{P}}_1} $和$ {\mathit{\boldsymbol{P}}_2} $的对应节点特征的欧氏距离得到，矩阵$ {\mathit{\boldsymbol{K}}_p} $中的元素为$k_{i_{1} i_{2}}^{p}=\phi\left(p_{i_{1}}^{1}, p_{i_{2}}^{2}\right)$，表示图$ {\mathit{\boldsymbol{A}}_1} $中的第$ {i_1} $个节点和图$ {\mathit{\boldsymbol{A}}_2} $中的第$ {i_2} $个节点的相似度。类似地，通过$ {\mathit{\boldsymbol{Q}}_1} $和$ {\mathit{\boldsymbol{Q}}_2} $可以计算图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的边特征相似矩阵$ {\mathit{\boldsymbol{K}}_q} $，行数等于图$ {\mathit{\boldsymbol{A}}_1} $的边数，列数等于图$ {\mathit{\boldsymbol{A}}_2} $的边数，计算方法通过$ {\mathit{\boldsymbol{Q}}_1} $和$ {\mathit{\boldsymbol{Q}}_2} $的对应节点特征的欧氏距离得到，矩阵$ {\mathit{\boldsymbol{K}}_q} $中的元素为$k_{c_{1} c_{2}}^{q}=\phi\left(q_{c_{1}}^{1}, q_{c_{2}}^{2}\right)$，表示图$ {\mathit{\boldsymbol{A}}_1} $中的第$ {c_1} $个边和图$ {\mathit{\boldsymbol{A}}_2} $中的第$ {c_2} $个边的相似度。

有了起点特征矩阵$ {\mathit{\boldsymbol{G}}_1} $与$ {\mathit{\boldsymbol{G}}_2} $、终点特征矩阵$ {\mathit{\boldsymbol{H}}_1} $与$ {\mathit{\boldsymbol{H}}_2} $、节点特征相似矩阵$ {\mathit{\boldsymbol{K}}_p} $以及边特征相似矩阵$ {\mathit{\boldsymbol{K}}_q} $，就可以将图匹配问题转化为求解一个节点的匹配关系矩阵$ \mathit{\boldsymbol{X}} $，矩阵中的元素$x_{i j} \in\{0, 1\}$，若$x_{i j}=1$，表示图$ {\mathit{\boldsymbol{A}}_1} $的第$ i $个节点和图$ {\mathit{\boldsymbol{A}}_2} $的第$ j $个节点相匹配；若$x_{i j}=0$，则表示无法匹配, 并且矩阵$ \mathit{\boldsymbol{X}} $的每行中有且仅有一个元素为1，每列中也有且仅有一个元素为1。这样图$ {\mathit{\boldsymbol{A}}_1} $与$ {\mathit{\boldsymbol{A}}_2} $的匹配问题就归结为寻求

$ {{J}_{\text{gm}}}(X)=\sum\limits_{{{i}_{1}}{{i}_{2}}}{{{x}_{{{i}_{1}}{{i}_{2}}}}}k_{{{i}_{1}}{{i}_{2}}}^{p}+\sum\limits_{\begin{smallmatrix} {{i}_{1}}\ne {{i}_{2}},{{j}_{1}}\ne {{j}_{2}} \\ \text{g}_{{{i}_{1}}{{c}_{1}}}^{1}h_{{{j}_{1}}{{c}_{1}}}^{1}=1n \\ \text{g}_{{{i}_{2}}{{c}_{2}}}^{2}h_{{{j}_{2}}{{c}_{2}}}^{2}=1 \end{smallmatrix}}{{{x}_{{{i}_{1}}{{i}_{2}}}}}{{x}_{{{j}_{1}}{{j}_{2}}}}k_{{{c}_{1}}{{c}_{2}}}^{q} $

(7)

的最大值。若想要式(7)取最大值，则其中的两个部分都需要得到最大值。第1部分中的$k_{i_{1} i_{2}}^{p}$表示两个图中对应节点的相似程度，如果其值很大，表明图$ {\mathit{\boldsymbol{A}}_1} $的$ {i_1} $节点和图$ {\mathit{\boldsymbol{A}}_2} $的$ {i_2} $节点相似程度很高，那么这两个节点很可能就会互相匹配，也就是说$x_{i_{1} i_{2}}=1$，这种情况下，如果$x_{i_{1} i_{2}}=0$，$x_{i_{1} i_{2}}$和$k_{i_{1} i_{2}}^{p}$的乘积就会为0，那么就浪费了这个非常相似的一对节点，也就无法使式(7)取最大值。第2部分中的$k_{c_{1} c_{2}}^{q}$表示两个图中对应边的相似程度，如果其值很大，表明图$ {\mathit{\boldsymbol{A}}_1} $的$ {c_1} $边和图$ {\mathit{\boldsymbol{A}}_2} $的$ {c_2} $边相似程度很高，那么$ {c_1} $和$ {c_2} $边的起始点很可能会互相匹配，而且$ {c_1} $和$ {c_2} $边的终点也很可能会互相匹配，通过式(7)中的$g_{i_{1} c_{1}}^{1} h_{j_{1} c_{1}}^{1}=1$可知，图$ {\mathit{\boldsymbol{A}}_1} $的$ {c_1} $边的起点为节点$ {i_1} $，终点为节点$ {j_1} $。同理，通过式(7)中的$g_{i_{2} c_{2}}^{2} h_{j_{2} c_{2}}^{2}=1$可知，图$ {\mathit{\boldsymbol{A}}_2} $的$ {c_2} $边的起点为节点$ {i_2} $，终点为节点$ {j_2} $，因此可以得到，这种情况下，如果$x_{i_{1} i_{2}}$和$x_{j_{1} j_{2}}$有一个为0，那么$x_{i_{1} i_{2}}、x_{j_{1} j_{2}}、k_{c_{1} c_{2}}^{q}$三者的乘积就会为0，那么就浪费了这个非常相似的一对边，也就无法使式(7)取最大值。通过以上讨论可知，当式(7)取最大值的时候，矩阵$ \mathit{\boldsymbol{X}} $表示图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的节点匹配关系。

为了求解式(7)，定义一个相似关系矩阵$ \mathit{\boldsymbol{K}} $，其元素定义如下

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{k_{{i_1}{i_2}, {j_1}{j_2}}} = \\ \left\{ {\begin{array}{*{20}{l}} {k_{{i_1}{i_2}}^p\;\;\;若{i_1} = {j_1}且{i_2} = {j_2}}\\ {k_{{c_1}{c_2}}^q\;\;\;若{i_1} \ne {j_1}且{i_2} \ne {j_2}且g_{{i_1}{c_1}}^1h_{{j_1}{c_1}}^1g_{{i_2}{c_2}}^2h_{{j_2}{c_2}}^2 = 1}\\ {0\;\;\;\;\;其他} \end{array}} \right.\\ \end{array} $

(8)

矩阵$\mathit{\boldsymbol{K}}$是通过矩阵$ {\mathit{\boldsymbol{K}}_p} $和$ {\mathit{\boldsymbol{K}}_q} $得到的，行数和列数都是$ {{n_1}} \times {{n_2}}$，$ {{n_1}} $和$ {{n_2}} $分别表示两个图$ {\mathit{\boldsymbol{A}}_1} $和$ {\mathit{\boldsymbol{A}}_2} $的顶点数。矩阵$ {\mathit{\boldsymbol{K}}_p} $按列矢量化后的元素，是矩阵$\mathit{\boldsymbol{K}}$的对角线元素；矩阵$\mathit{\boldsymbol{K}}$的第 $\left(a_{1} \times b_{1}, a_{2} \times b_{2}\right)$个元素的值，等于矩阵$ {\mathit{\boldsymbol{K}}_q} $的第$(m, n)$个元素，$ m $、$ n $与${a_1}$、${a_2}$、${b_1}$、${b_2}$的关系推导如下：矩阵$ {\mathit{\boldsymbol{K}}_q} $的第$ (m, n) $个元素表示图$ {\mathit{\boldsymbol{A}}_1} $的第$ m $个边和图$ {\mathit{\boldsymbol{A}}_2} $的第$ n $个边的相似程度，图$ {\mathit{\boldsymbol{A}}_1} $的第$ m $个边，可以由矩阵$ {\mathit{\boldsymbol{G}}_1} $和$ {\mathit{\boldsymbol{H}}_1} $中第$ m $列取值为1的元素来表示，图$ {\mathit{\boldsymbol{A}}_2} $的第$ n $个边，可以由矩阵$ {\mathit{\boldsymbol{G}}_2} $和$ {\mathit{\boldsymbol{H}}_2} $中第$ n $列取值为1的元素来表示，因为$ {\mathit{\boldsymbol{G}}_1} $、$ {\mathit{\boldsymbol{H}}_1} $、$ {\mathit{\boldsymbol{G}}_2} $、$ {\mathit{\boldsymbol{H}}_2} $矩阵中，每列只有一个元素为1，因此这个为1的元素可以标记为矩阵$ {\mathit{\boldsymbol{G}}_1} $中的$({a_1}, m)$，$ {\mathit{\boldsymbol{H}}_1} $中的$({b_1}, m)$, $ {\mathit{\boldsymbol{G}}_2} $中的$({2}, n)$，$ {\mathit{\boldsymbol{H}}_2} $中的$({b_2}, n)$，即$ m $、$ n $与${a_1}$、${a_2}$、${b_1}$、${a_2}$的关系。

有了矩阵$\mathit{\boldsymbol{K}}$，式(7)对节点匹配关系矩阵$ \mathit{\boldsymbol{X}} $的求解可以转换为

$ \mathop {\max }\limits_X {J_{{\rm{gm}}}}(X) = {\mathop{\rm vec}\nolimits} {(\mathit{\boldsymbol{X}})^{\rm{T}}}\mathit{\boldsymbol{K}}{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{X}}) $

(9)

式中，$ {\rm{vec}} $表示将矩阵向量化，对式(9)的求解，借鉴文献[7, 11]的方法，首先将矩阵$\mathit{\boldsymbol{K}}$分解为$ {\mathit{\boldsymbol{G}}_1} $、$ {\mathit{\boldsymbol{G}}_2} $、$ {\mathit{\boldsymbol{H}}_1} $、$ {\mathit{\boldsymbol{H}}_2} $、$ {\mathit{\boldsymbol{K}}_p} $、$ {\mathit{\boldsymbol{K}}_q} $等6个小矩阵，表示为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}} = {\mathop{\rm diag}\nolimits} \left( {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{K}}_p}} \right)} \right) + }\\ {\left( {{\mathit{\boldsymbol{G}}_1} \otimes {\mathit{\boldsymbol{G}}_2}} \right){\mathop{\rm diag}\nolimits} \left( {{\mathop{\rm vec}\nolimits} \left( {{\mathit{\boldsymbol{K}}_q}} \right)} \right){{\left( {{\mathit{\boldsymbol{H}}_2} \otimes {\mathit{\boldsymbol{H}}_1}} \right)}^{\rm{T}}}} \end{array} $

(10)

式中，$ {\rm{diag}} $表示一个对角元素为括号中的向量的对角矩阵，$ \otimes $为克罗内克积。将式(10)代入式(9)，则有

$ \begin{array}{*{20}{c}} {\mathop {\max }\limits_X {J_{{\rm{gm}}}}(\mathit{\boldsymbol{X}}) = }\\ {{\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{K}}_p^{\rm{T}}\mathit{\boldsymbol{X}}} \right) + {\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{K}}_q^{\rm{T}}\left( {\mathit{\boldsymbol{G}}_1^{\rm{T}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{G}}_2^ \circ} \mathit{\boldsymbol{H}}_1^{\rm{T}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{H}}_2}} \right)} \right)} \end{array} $

(11)

式中，${\mathop{\rm tr}\nolimits} \left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{B}}} \right) = {\mathop{\rm vec}\nolimits} {(\mathit{\boldsymbol{A}})^{\rm{T}}}{\mathop{\rm vec}\nolimits} (\mathit{\boldsymbol{B}})$。对矩阵$ {\mathit{\boldsymbol{K}}_q} $进行奇异值分解(SVD)，则式(11)转化为

$ \mathop {\max }\limits_X {J_{{\rm{gm}}}}(\mathit{\boldsymbol{X}}) = {\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{K}}_p^{\rm{T}}\mathit{\boldsymbol{X}}} \right) + \sum\limits_{i = 1}^c {{\mathop{\rm tr}\nolimits} } \left( {\mathit{\boldsymbol{A}}_i^1\mathit{\boldsymbol{XA}}_i^2{\mathit{\boldsymbol{X}}^{\rm{T}}}} \right) $

(12)

式中，$\mathit{\boldsymbol{A}}_i^1 = {\mathit{\boldsymbol{G}}_1}{\mathop{\rm diag}\nolimits} \left( {{u_i}} \right)\mathit{\boldsymbol{H}}_1^{\rm{T}}, \mathit{\boldsymbol{A}}_i^2 = {\mathit{\boldsymbol{G}}_2}{\mathop{\rm diag}\nolimits} \left( {{v_i}} \right)\mathit{\boldsymbol{H}}_2^{\rm{T}}$。

式(12)的求解是一个二次分配问题(QAP)，这类问题的求解一直是数学上的难题，本文采用路径跟踪算法(PFA)^[11]求其近似解。路径跟踪算法是通过将待求目标转化为一系列凹、凸表达的累加，并通过迭代逼近的方式求得最优解。因此，本文需将式(12)转化为凹、凸表达，为此引入常数$ \gamma $，其值为

$ \gamma = \sum\limits_{i = 1}^c {{\mathop{\rm tr}\nolimits} } \left( {\mathit{\boldsymbol{A}}_i^{1{\rm{T}}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{A}}_i^1} \right) + {\mathop{\rm tr}\nolimits} \left( {\mathit{\boldsymbol{A}}_i^2\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{A}}_\mathit{i}^{2{\rm{T}}}} \right) $

(13)

则有

$ {J_{{\rm{vex}}}}(\mathit{\boldsymbol{X}}) = {J_{{\rm{gm}}}}(\mathit{\boldsymbol{X}}) - \frac{\gamma }{2} $

$ {J_{{\rm{cav}}}}(\mathit{\boldsymbol{X}}) = {J_{{\rm{gm}}}}(\mathit{\boldsymbol{X}}) - \frac{\gamma }{2} $

式中，${J_{{\rm{vex}}}}$为凸函数，${J_{{\rm{cav}}}}$为凹函数, 因此，式(12)可以转化为

$ \begin{array}{l} \mathop {\max }\limits_X {J_\alpha }(\mathit{\boldsymbol{X}}) = (1 - \alpha ){J_{{\rm{vex}}}}(\mathit{\boldsymbol{X}}) + \alpha {J_{{\mathop{\rm cav}\nolimits} }}(\mathit{\boldsymbol{X}}) = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{J_{{\rm{gm}}}}(\mathit{\boldsymbol{X}}) + \left( {\alpha - \frac{1}{2}} \right)\gamma \end{array} $

(14)

式中，$\alpha \in[0, 1]$。那么，可以通过迭代方式求解式(14)，即：

1) 将$ \mathit{\boldsymbol{X}} $初始化为$ \mathit{\boldsymbol{X}}_0 $，将$ \alpha $初始化为0;

2) 将$ \mathit{\boldsymbol{X}} $、$ \alpha $代入式(14);

3) 应用Frank-Wolfe’s算法^[12]，计算得到最大值${J_\alpha }\left( {{\mathit{\boldsymbol{X}}^*}} \right)$;

4) 若${J_{{\rm{gm}}}}\left( {{\mathit{\boldsymbol{X}}^*}} \right) > {J_{{\rm{gm}}}}\left( {{\mathit{\boldsymbol{X}}_0}} \right)$，则用${\mathit{\boldsymbol{X}}^*}$取代$ \mathit{\boldsymbol{X}}_0 $，$\alpha=\alpha+\nabla \alpha$，$\nabla \alpha$表示一个小的增量;

5) 重复步骤2), 直到函数收敛或达到指定迭代次数为止。

至此，求解得到的矩阵$ \mathit{\boldsymbol{X}} $，即为两个图的节点匹配矩阵。

实验使用的三维图像数据是从网上(https://free3d.com)免费下载，共下载了9个汽车的三维图像模型，并导入三维空间坐标系中，如图 3所示。随后在这些三维图像上标记特征点，本实验手动标记了102个特征点，特征点集中不包含错误特征点，在每个轮子上标记了4个特征点；在车的前侧与后侧分别标记了4行特征点，每行3个点；在左右增加了一些点来更细致地描述轮廓，如图 4所示。

本实验共分为5个部分，前4个是对手工选取的特征点进行匹配，第5个是对算法自动生成的特征点进行匹配，具体为：1)对相同的图像进行特征点匹配；2)对不同图像进行特征点匹配；3)对相同的图像旋转不同角度后，再进行特征点匹配；4)对相同的图像随机去掉一定数量的节点后，再进行特征点匹配；5)对3D SIFT算法得到的特征点进行匹配。

对9个三维图像模型，各自进行特征点匹配，每个三维图像的匹配准确率如表 1所示。可以看到，在面向相同图像时，基于本文算法的三维图像特征点匹配准确率平均为97.56%，其中同一个三维图像模型的匹配结果如图 5所示。在匹配结果中，绿色线连接的为匹配正确的节点对，红色线连接的为匹配错误的节点对。

对三维图像数据集中的每两个图像之间都进行特征点匹配，因为数据集都是各种车型的三维图像，各种车型之间存在一定的相似性，因此可以进行匹配，此时的匹配依据仍然是三维图起点、终点特征矩阵以及三维图特征相似矩阵，同时又考虑到不同车型的三维图像终归会存在一定程度上的差异性，这会使得一部分本该匹配在一起的特征点因为局部结构的差异性，从而引起一定程度的误匹配。两两之间的匹配结果如表 2所示，平均准确率为76.39%。值得注意的是，表 2中的数据呈现对称的模式，因为从第$ M $个模型向第$ N $个模型进行匹配的结果，与从第$ N $个模型向第$ M $个模型进行匹配的结果是相等的($ 0≤M≤9 $, $ 0≤N≤9 $)。考虑到这是不同车型的匹配效果，在多数情况下，匹配的准确率比较高，特别是轿车与SUV这样差异较大的三维图像数据之间的匹配效果也是很准确的(如图 6所示，图中绿色线连接的为匹配正确的节点对，红色线连接的为匹配错误的节点对)。在102个特征点中，准确匹配85个，误匹配17个(其中车头3个、前轮2个、车身1个、后轮7个、车尾4个)。通过分析图 6中各个角度以及局部放大的匹配结果可知，车头、前轮、后轮特征点的匹配错误是因为存在与待匹配特征点对称的点，例如：左前轮和右前轮的特征点因为图特征非常相似，导致匹配错误。车身和车尾的特征点的匹配错误是因为轿车三维模型与SUV车三维模型存在较大差异，使得待匹配特征点与错误的点在图特征上更相似，导致匹配错误。

为了测试本文算法对三维图像刚性变化情况下的匹配准确率，将三维图像绕$ Y $轴旋转若干角度，然后将旋转前的图像与旋转后的图像进行特征点匹配。三维图像从0°旋转至180°，每隔10°进行一次匹配，共进行18次匹配。旋转角度与准确率的直接关系曲线如图 7所示，可以看到随着旋转角度的增加，特征点匹配的准确率并没有降低，始终维持在95%以上，说明本文算法不受三维图像的刚性旋转的影响。其中一个三维图像数据旋转后的匹配结果如图 8所示。图 8中，绿色线连接的是匹配正确的节点对，红色线连接的是匹配错误的节点对。旋转40°的时候，错误匹配2个特征点；旋转90°的时候，错误匹配3个特征点; 旋转130°的时候，错误匹配2个特征点；旋转180°的时候，错误匹配2个特征点。因为本文选取了图的边与$ XOY $平面的夹角作为图特征矩阵的一部分，该夹角会随着三维图像的旋转而变化，从而造成个别特征点的匹配出现错误。

为了测试本文算法对三维图像非刚性变化情况下的匹配准确率，随机去掉1~30个特征点后再进行匹配。由于每次匹配时，被匹配的三维图像中去掉的特征点是随机的，所以不是所有的特征点都是一一对应的，因此匹配准确率的计算与前3个实验略有不同，由匹配正确的特征点对的数量与可以匹配的特征点对的数量的比值作为准确率。随机去掉若干特征点后，再进行特征点匹配的准确率曲线如图 9所示，可以看到，随着随机去掉特征点的数目的增加，匹配的准确率也在逐渐下降，这是因为待匹配的特征点受到随机删除的影响，使得本该与之最佳匹配的点变成了非最佳，从而导致误匹配。但是即使是去掉了30个特征点的情况下，匹配准确率也在70%以上。说明本文算法受三维图像非刚性变化的影响较小。其中一个三维图像数据在随机去掉10个、20个、30个特征点后的匹配结果如图 10所示，其中绿色线连接的是匹配正确的节点对，红色线连接的是匹配错误的节点对。

之前的验证过程中，特征点都是手工标注的，为了进一步验证本文算法，通过三维SIFT算法^[13]计算得到36个特征点，并应用本文算法对这些三维特征点进行匹配。选取的三维图像仍然是前文所述的9个汽车模型，图 11展示了其中一个的特征点匹配结果，准确率为100%。所有9个三维汽车模型的特征点匹配准确率结果如表 3所示，平均匹配准确率为98.78%。可以验证，本文算法对手工标注的三维特征点和特征提取算法得到的三维特征点的匹配效果都很好。

本文面向三维图像数据，提出了三维图匹配算法，解决了现有图匹配算法仅适用于二维图像特征点匹配的问题。实验表明，本文算法对相同三维图像的特征点匹配有97.56%的平均准确率；对不同三维图像特征点匹配有76.11%的平均准确率；在三维图像有旋转的情况下，有90%以上的平均准确率；在特征点部分缺失的情况下，平均匹配准确率也能达到80%。综上，本文提出的三维图匹配算法可以很好地对三维图像的特征点进行匹配，但是在面向不同三维图像特征点匹配的时候，其准确率还有待提高。今后的研究方向是对刚性及非刚性变换条件下的三维图像做特征点匹配。