形状匹配^[1]是计算机视觉中的重要研究内容，它在图像检索、目标识别、图形图像拼接等方面都有重要的应用。在基于形状匹配的目标识别中，由于受传感器视点变化的影响，形状会发生相似变换、形变等。通常，可以通过仿射变换模型来描述这种变化，即对于在不同传感器或不同视角下获得的相同目标的形状，它们之间存在仿射变换的关系。根据所利用的形状信息，仿射变换下的形状匹配分为基于轮廓和区域的两种匹配类型。

基于轮廓的匹配方法通过对轮廓上采样点之间的关系形成描述子进行匹配，例如Mokhtarian等人^[2]提出的曲率尺度空间(CSS)方法。CSS基于多尺度和曲率来表示形状，具有噪声鲁棒性、尺度不变性和旋转不变性，可靠并易于计算。文献[3]对其进行改进，使用曲率积分算子，使其提出的描述子具有除水平剪切外的仿射不变性，但是这种方法不能保证采样点在轮廓上起始位置的一致性，从而使等曲率积分的采样点位置不一致，最终导致匹配误差。因为曲率不具备完全仿射不变性，文献[4]定义了仿射曲率算子，但这种算子计算复杂。Fu等人^[5]在仿射曲率算子的基础上提出一种新的仿射曲率描述子，将四阶的算子降为三阶，降低了计算难度；但这种描述子在拐点处没有定义，而在实际情况中，接近拐点的位置会导致该描述子的值非常大，使得该描述子不稳定。Fan^[6]与Prasad^[7]提出根据形状上的几何元素(直线、曲线、多边形)构造不变量的方法，适合于相对简单的形状。张桂梅等人^[8]提出等面积分割点解决轮廓曲线平滑特征点少而不能被精确描述的问题。文献[9]将弦高点与遗传算法相结合用于图像的仿射配准，但弦高点序列获取和设置的起始点对位置密切相关。蔡慧英等人^[10]使用凸包作为仿射变换的不变点，不仅具有一定的稳健性，而且提取相对简单，但是当凸包的顶点较少时，得到的特征点也较少，造成描述形状能力不足。以上方法都需要确定的、一致性的轮廓起始点，而在图像处理中，采用的方法或尺度不同会导致确定的轮廓曲线起始点不一致。因此，在对曲线匹配前，需要确定一致的起始点。

基于区域的匹配方法利用目标形状的内部信息形成描述子进行匹配，例如，Flusser等人^[11]提出的仿射不变矩，通过对形状计算不同阶次的矩来构造仿射不变量。李迎春等人^[12]将仿射不变矩与神经网络进行结合用于飞机目标的识别，对仿射目标取得了较好的效果。文献[13]利用ICA Zernike矩描述子用于形状检索，该方法首先通过独立成分分析(ICA)将原始形状规范化，规范化的形状具有了尺度不变性；然后进行Zernike变换，使图像具有旋转不变性；最后提取Zernike矩描述子，使得该描述子具有仿射不变性。Wang等人^[14]提出了一种结构简单而且稳定的仿射不变描述子，该描述子源于矩阵变换理论，可用于形状的匹配与检索，但是计算复杂度较高，而且在仿射变换上存在局限性。以上形状表示方法的缺点是他们无法反映物体内部的形状内容，同时这些方法也很难处理包含有两个或更多分离部分的复杂形状。

文献[15]为了获取更丰富的形状信息，结合形状区域与轮廓关系提出一种基于形状投影分布的描述子(SPD)，该描述子将目标点集视为纵向长度片段，投影在重心与轮廓点的连线垂直线上，从而得到形状描述向量。该描述子同时利用了目标点集的空间位置关系与轮廓特征，性能良好，但其匹配采用等间隔采样，对复杂图形效果不好；另外，其采用统计目标所有图像点的方法计算面积比，使得计算用时高。这种方法的优势是它们将形状轮廓的整体信息与局部信息有效地融合在一起，描述形状时更加准确。

针对仿射变换下形状匹配中存在的描述子对形状的描述能力不足，以及描述子计算耗时大的问题，本文改进文献[15]的匹配方法，利用轮廓曲线计算投影面积分布描述子对特征点进行匹配，并对匹配度前5%的目标进行对应的轮廓子曲线匹配。

文献[16]论证了投影不变性：在垂直于两个指定点连线方向上，任一形状点的投影在整个形状投影中的相对位置，仿射变换前后保持不变。利用这个性质可以构造基于投影面积的形状描述子SPD。

如图 1所示，目标所占的区域为$\mathit{\boldsymbol{D}}$，目标的轮廓($\mathit{\boldsymbol{D}}$的边界)记为$E$。对于特征点$P$，以形状质心$o$到$P$的向量${\mathit{\boldsymbol{V}}_{op}}$作${x_p}$轴建立右手直角坐标系$(o, {\mathit{\boldsymbol{x}}_p}, {\mathit{\boldsymbol{y}}_p})$。那么，对形状区域$\mathit{\boldsymbol{D}}$内的目标点向${\mathit{\boldsymbol{y}}_p}$的投影统计为

$ \begin{array}{*{20}{c}} {{h_p}\left( k \right) = }\\ {\# \left\{ {g:{{\boldsymbol{v}}_{og \cdot }}{{\mathit{\boldsymbol{\hat y}}}_p} \in \left[ {{v^p}\left( {k - 1} \right),{v^p}\left( k \right)} \right],g \in \mathit{\boldsymbol{D}}} \right\}}\\ {{v^p}\left( k \right) = v_{\min }^p + k \cdot \left( {v_{\max }^p - v_{\min }^p} \right)/N;k = 1, \cdots ,N} \end{array} $

(1)

式中，$\mathit{\boldsymbol{D}}$代表形状，$g$为$\mathit{\boldsymbol{D}}$中的目标点，${\mathit{\boldsymbol{v}}_{og}}$为点$o$到点$g$的向量，${\mathit{\boldsymbol{\hat y}}_p}$为${\mathit{\boldsymbol{y}}_p}$轴的单位向量，${v^p}\left(k \right)$为形状点在${\mathit{\boldsymbol{y}}_p}$的投影值，${h_p}\left(k \right)$为$\mathit{\boldsymbol{D}}$中在${\mathit{\boldsymbol{y}}_p}$的投影值位于$({v^p}\left({k-1} \right), {v^p}\left(k \right)]$范围内的目标点的个数，$v_{\max }^p$、$v_{\min }^p$分别为投影的最大值与最小值，$N$为投影等间隔划分数目。显然，${h_p}\left(k \right)$的总和为$\mathit{\boldsymbol{D}}$中点的总数$S$，也就是形状面积，为了消除缩放对面积的影响，需要利用$S$对${h_p}\left(k \right)$进行归一化，即

$ {d_p}\left( k \right) = {h_p}\left( k \right)/\sum\limits_{n = 1}^N {{h_p}\left( n \right)} $

(2)

该描述子除了将形状看成柱状直条，忽视了形状的内部结构外，还须对形状区域$\mathit{\boldsymbol{D}}$内的所有目标点进行统计，使得计算耗时非常大。本文对其进行改进，提出快速SPD描述子(简称QSPD)。

如图 2所示，在闭合有向的轮廓序列$\mathit{\boldsymbol{E}}({g_1}, {g_2}, \ldots, {g_k}, {g_{k + 1}}, \ldots, {g_n})$内，以轮廓质心$c$为原点，$c$到轮廓特征点$P$的向量${\mathit{\boldsymbol{V}}_{cp}}$作${\mathit{\boldsymbol{x}}_p}$轴建立右手直角坐标系$(c, {\mathit{\boldsymbol{x}}_p}, {\mathit{\boldsymbol{y}}_p})$。将整个形状区域按照在${\mathit{\boldsymbol{y}}_p}$的投影范围划分为$N$等份，并按照${\mathit{\boldsymbol{y}}_p}$的两侧位置划分为2$N$个区间$\left({d\left(1 \right), d\left(2 \right), \ldots, d\left({2N} \right)} \right)$。由于轮廓是闭合的，前$N$个区间用于表示${\mathit{\boldsymbol{y}}_p}$左侧轮廓的投影，后$N$个区间用于表示${\mathit{\boldsymbol{y}}_p}$右侧轮廓的投影。等分线与轮廓曲线的交点分别记为$G_1, {G_2}, \ldots, {G_t}$，由于轮廓的弯曲不确定，$t$的值并不确定。在图 2中，$N$取10，若将有向轮廓序列看成连续的曲线，则轮廓曲线在$d$(4)区域的投影面积为

$ {{S}_{d\left( 4 \right)}}=\left| \int{\overset\frown{{{G}_{4}}{{G}_{5}}}\text{d}s}+\int{\overset\frown{{{G}_{5}}{{G}_{6}}}\text{d}s}+\int{\overset\frown{{{G}_{6}}{{G}_{7}}}\text{d}s} \right| $

(3)

由于轮廓曲线是有向的，因此$\int{\overset\frown{{{G}_{5}}{{G}_{6}}}\text{d}s}$的符号和$\int{\overset\frown{{{G}_{4}}{{G}_{5}}}\text{d}s}$、$\int{\overset\frown{{{G}_{6}}{{G}_{7}}}\text{d}s}$相反。不管轮廓方向如何，式(3)可以得到轮廓曲线的投影面积。同理，对于$d$(7)，有

$ {{S}_{d\left( 7 \right)}}=\left| \int{\overset\frown{{{G}_{9}}{{G}_{10}}}\text{d}s}+\int{\overset\frown{{{G}_{12}}{{G}_{13}}}\text{d}s} \right| $

(4)

对于离散轮廓序列$\mathit{\boldsymbol{E}}({{g}_{1}}, {{g}_{2}}, \ldots, {{g}_{k}}, {{g}_{k+1}}, \ldots, {{g}_{n}})$，式(3)和式(4)中的积分表示为区域内轮廓相邻两点${{g}_{k}}, {{g}_{k+1}}$与其投影点构成的梯形面积${{s}_{k}}$的累加和。投影面积${{s}_{k}}$计算为

$ {s_k} = 0.5\left( {{y_k} - {y_{k + 1}}} \right)\left( {{\delta _k} + {\delta _{k + 1}}} \right) $

(5)

式中，${{y}_{k}}$为点${{g}_{k}}$在${\mathit{\boldsymbol{y}}_p}$轴上的投影，${{\delta }_{k}}$为点${{g}_{k}}$到${\mathit{\boldsymbol{y}}_p}$轴的距离。

通过式(5)可以得到点集投影面积序列${s_1}, {s_2}, \ldots, {s_{n-1}}$，联合式(3)、(4)对轮廓序列$\mathit{\boldsymbol{E}}({g_1}, \ldots, {g_n})$的投影面积统计为

$ \begin{array}{l} {h_p}\left( k \right) = \sum {\left( {{s_i}} \right)\left\{ {{g_i}:{\mathit{\boldsymbol{v}}_{cg \cdot }}{{\mathit{\boldsymbol{\hat y}}}_p} \in \left[ {{v^p}\left( {k - 1} \right),{v^p}\left( k \right)} \right],} \right.} \\ \left. {{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{v}}_{cg}} \times {{\mathit{\boldsymbol{\hat y}}}_p}} \right) \ge 0} \right\}\\ {h_p}\left( {2k} \right) = \sum {\left( {{s_i}} \right)\left\{ {{g_i}:{\mathit{\boldsymbol{v}}_{cg \cdot }}{{\mathit{\boldsymbol{\hat y}}}_p} \in \left[ {{v^p}\left( {k - 1} \right),{v^p}\left( k \right)} \right],} \right.} \\ \left. {{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{v}}_{cg}} \times {{\mathit{\boldsymbol{\hat y}}}_p}} \right) < 0} \right\}\\ {v^p}\left( k \right) = v_{\min }^p + k\left( {v_{\max }^p - v_{\min }^p} \right)/N,k = 1, \cdots ,N\\ {v^p}\left( 0 \right) = v_{\min }^p \end{array} $

(6)

${h_p}\left(k \right)$为$E$中在${\mathit{\boldsymbol{y}}_p}$的投影位于$({v^p}\left({k-1} \right), {v^p}\left(k \right)]$范围内的点构成的${s_k}$的累加和，${\rm{sgn}}({\mathit{\boldsymbol{v}}_{cg}} \times {\mathit{\boldsymbol{y}}_p})$用来判断投影的方向, $v_{\max }^p$、$v_{\min }^p$分别为投影的最大与最小值。通过式(6)最终可以得到2$N$维度的${h_p}$。对${h_p}\left(k \right)$进行归一化

$ \begin{array}{*{20}{c}} {{d_p}\left( k \right) = {h_p}\left( k \right)/\sum\limits_{n = 1}^{2N} {{h_p}\left( n \right)} ;}&{k = 1, \cdots ,2N} \end{array} $

(7)

如图 3所示，图 3(a)为原目标，通过仿射变换矩阵$\mathit{\boldsymbol{T}} = \{ {\mathit{\boldsymbol{A}}_T}, \mathit{\boldsymbol{b}}\} $将图 3(a)中的骆驼变换到图 3(b)，$\mathit{\boldsymbol{b}}$为平移向量，${\mathit{\boldsymbol{A}}_T}$为仿射矩阵，$\mathit{\boldsymbol{E}}$、$\mathit{\boldsymbol{E'}}$分别为两个目标的轮廓曲线。

仿射变换包含缩放、旋转、平移和水平剪切几种变换，QSPD构造的坐标系是基于轮廓的，对于平移具有不变性。缩放、旋转和水平剪切可以分别表示成变换

$ \begin{array}{l} {\mathit{\boldsymbol{A}}_{{\rm{Tscaling}}}} = \left[ {\begin{array}{*{20}{c}} {{s_x}}&0\\ 0&{{s_y}} \end{array}} \right]\\ {\mathit{\boldsymbol{A}}_{{\rm{Trotation}}}} = \left[ {\begin{array}{*{20}{c}} {\cos \left( \theta \right)}&{ - \sin \left( \theta \right)}\\ {\sin \left( \theta \right)}&{\cos \left( \theta \right)} \end{array}} \right]\\ {\mathit{\boldsymbol{A}}_{{\rm{Tshcaring}}}} = \left[ {\begin{array}{*{20}{c}} 1&k\\ 0&1 \end{array}} \right] \end{array} $

(8)

对于图 3(a)的轮廓点$P$，随机点$g$在${y_p}$上的投影值$v_g^p$计算为

$ v_g^p = {\mathit{\boldsymbol{v}}_{cg}} \cdot {{\mathit{\boldsymbol{\hat y}}}_p} $

(9)

式中，${{\mathit{\boldsymbol{\hat y}}}_p}$为${\mathit{\boldsymbol{y}}_p}$轴的单位向量，${\mathit{\boldsymbol{v}}_{cg}}$为点$c$到点$g$的向量，显然，式(9)也能表达为

$ v_g^p = {\mathit{\boldsymbol{v}}_{c{o_g}}} \cdot {{\mathit{\boldsymbol{\hat y}}}_p} $

(10)

式中，点${o_g}$是点$g$在${y_p}$上的投影，${\mathit{\boldsymbol{v}}_{c{o_g}}}$是点$c$到点${o_g}$的向量。假如点$P$和点$g$分别映射到$\mathit{\boldsymbol{E}}{\rm{'}}$中的点$P′$和点$g′$，即$P' = {\mathit{\boldsymbol{A}}_T}P + \mathit{\boldsymbol{b}}、g' = {\mathit{\boldsymbol{A}}_T}g + \mathit{\boldsymbol{b}}$。根据质心的仿射不变性，可以得出推

$ \begin{array}{*{20}{c}} {v_{g'}^{p'} = {\mathit{\boldsymbol{v}}_{c'g'}} \cdot {{\mathit{\boldsymbol{\hat y}}}_{p'}} = \left( {{{\mathit{\boldsymbol{\hat v}}}_{c'p'}} \times {\mathit{\boldsymbol{v}}_{c'g'}} \times {{\mathit{\boldsymbol{\hat v}}}_{c'p'}}} \right) \cdot {{\mathit{\boldsymbol{\hat y}}}_{p'}} = }\\ {\det \left( {{\mathit{\boldsymbol{A}}_T}} \right) \cdot \frac{{\left| {{\mathit{\boldsymbol{v}}_{c'p'}}} \right|}}{{\left| {{\mathit{\boldsymbol{v}}_{cp}}} \right|}} \cdot \left( {\left( {{{\mathit{\boldsymbol{\hat v}}}_{cp}} \times {\mathit{\boldsymbol{v}}_{cg}} \times {{\mathit{\boldsymbol{\hat v}}}_{cp}}} \right) \cdot {{\mathit{\boldsymbol{\hat y}}}_{p'}}} \right) = }\\ {\det \left( {{\mathit{\boldsymbol{A}}_T}} \right) \cdot \frac{{\left| {{\mathit{\boldsymbol{v}}_{c'p'}}} \right|}}{{\left| {{\mathit{\boldsymbol{v}}_{cp}}} \right|}} \cdot \mathit{\boldsymbol{v}}_g^p} \end{array} $

(11)

式中，$v_{g'}^{p'}$是$g′$在${\mathit{\boldsymbol{y}}_{p'}}$上的投影值，${{\mathit{\boldsymbol{\hat v}}}_{c'p'}}$、${{\mathit{\boldsymbol{\hat v}}}_{cp}}$分别是${\mathit{\boldsymbol{v}}_{c'p'}}$、${\mathit{\boldsymbol{v}}_{cp}}$上的单位向量，${\mathit{\boldsymbol{A}}_T}$为仿射矩阵，$\det \left(\mathit{\boldsymbol{\xi }} \right)$是矩阵$\mathit{\boldsymbol{\xi }}$的行列式，${\left| \mathit{\boldsymbol{\xi }} \right|}$是向量$\mathit{\boldsymbol{\xi }}$的模。根据式(11)，有

$ \begin{array}{l} \left| {{\mathit{\boldsymbol{v}}_{cp}}} \right| \cdot \mathit{\boldsymbol{v}}_{\max }^{p'} = \det \left( {{\mathit{\boldsymbol{A}}_T}} \right) \cdot \left| {{\mathit{\boldsymbol{v}}_{c'p'}}} \right| \cdot \mathit{\boldsymbol{v}}_{\max }^p\\ \left| {{\mathit{\boldsymbol{v}}_{cp}}} \right| \cdot \mathit{\boldsymbol{v}}_{\min }^{p'} = \det \left( {{\mathit{\boldsymbol{A}}_T}} \right) \cdot \left| {{\mathit{\boldsymbol{v}}_{c'p'}}} \right| \cdot \mathit{\boldsymbol{v}}_{\min }^p \end{array} $

(12)

将式(12)代入式(6)得到

$ {v^{p'}}\left( k \right) = \det \left( {{\mathit{\boldsymbol{A}}_T}} \right) \cdot \frac{{\left| {{\mathit{\boldsymbol{v}}_{c'p'}}} \right|}}{{\left| {{\mathit{\boldsymbol{v}}_{cp}}} \right|}} \cdot {\mathit{\boldsymbol{v}}^p}\left( k \right) $

(13)

联合式(11)和式(13)可得

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{v}}_{cp}} \cdot {{\mathit{\boldsymbol{\hat y}}}_p} \in \left( {{v^p}\left( {k - 1} \right),{v^p}\left( k \right)} \right] \Leftrightarrow }\\ {{\mathit{\boldsymbol{v}}_{c'g'}} \cdot {{\mathit{\boldsymbol{\hat y}}}_{p'}} \in \left( {{v^{p'}}\left( {k - 1} \right),{v^{p'}}\left( k \right)} \right]}\\ {{\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{v}}_{cg}} \times {{\mathit{\boldsymbol{\hat y}}}_p}} \right) = {\mathop{\rm sgn}} \left( {{\mathit{\boldsymbol{v}}_{c'g'}} \times {{\mathit{\boldsymbol{\hat y}}}_{p'}}} \right)} \end{array} $

(14)

式(14)表明QSPD的仿射不变性。即${d_p}\left(k \right) = {d_{p\prime }}\left(k \right)$。

CSS^[2]角点具备仿射不变性，本文选其作为备选特征点，使用QSPD描述子结合动态规划算法(蚁群算法)^[17]对两幅图像轮廓的特征点进行匹配。

在利用蚁群算法进行匹配定位之前，需要先定义描述子的相似性度量，本文采用卡方距离(chisquare)作为度量距离。蚁群算法中，蚂蚁每走一步需要对路径的信息素进行更新。针对特征点匹配问题，取匹配适应度函数来更新信息素。当相似性度量采用卡方距离时，适应度函数定义为

$ \begin{array}{*{20}{c}} {R\left( {x,y} \right) = }\\ {\sqrt {\sum\limits_{i = 1}^N {{{\left( {\frac{{{x_i} - {\rm{E}}\left( {{x_i}} \right)}}{{{\rm{E}}\left( {{x_i}} \right)}}} \right)}^2}} + \sum\limits_{i = 1}^N {{{\left( {\frac{{{y_i} - {\rm{E}}\left( {{y_i}} \right)}}{{{\rm{E}}\left( {{y_i}} \right)}}} \right)}^2}} } } \end{array} $

(15)

式中, ${x_i}$为描述子$x$的第$i$个变量的取值，${\rm{E}}\left({{x_i}} \right)$是描述子$x$在第$i$个变量的期望频数，${y_i}$为描述子$y$的第$i$个变量的取值，${\rm{E}}\left({{y_i}} \right)$是描述子$y$在第$i$个变量的期望频数。

将两组特征点视为城市，将特征点匹配转化为旅行商问题，以下为特征点匹配的流程：

1) 蚁群初始化。首先对蚁群算法的基本参数进行初始化，设算法中的蚂蚁个数$M$，最大迭代次数$T$、搜索空间信息素$τ$(0)、挥发系数$ρ$以及信息素强度系数$Q$，$M$个蚂蚁随机分布在两列向量中。

2) 爬行操作。禁忌表${U_{ab{u_k}}}\left({k = 1, 2, \ldots } \right)$记录了蚂蚁走过的路径，$M$个蚂蚁按照转移概率在两列城市间跳转，转移概率为

$ P_{ij}^m\left( t \right) = \frac{{{{\left| {{\tau _{ij}}\left( t \right)} \right|}^\alpha }}}{{\sum {{{\left| {{\tau _{ij}}\left( t \right)} \right|}^\alpha }} }} $

(16)

式中，${\tau _{ij}}\left(t \right)$表示匹配定位参数$i$, $j$上的信息素，$α$为信息启发式因子，表示轨迹的相对重要性，其值越大，则该蚂蚁越倾向于选择其他蚂蚁经过的路径。蚂蚁$k$从$i$构成的城市列按照式(16)跳到$j$构成的城市列，一次跳转对应一个解，即一个匹配位置。

3) 刷新信息素。依照蚂蚁跳转后的位置上的适应度函数更新信息素$\Delta \tau _{ij}^k\left(t \right)$

$ \Delta \tau _{ij}^k\left( t \right) = \left\{ \begin{array}{l} \frac{Q}{{R\left( {x,y} \right)}}\;\;\;\;\;第\;k\;只蚂蚁经过\left( {i,j} \right)\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(17)

式中, $Q$为信息素强度系数。所有蚂蚁多次搜索后会聚集到适应度函数最小的路径，即最佳的匹配位置。

4) 返回步骤2)直到最大搜索次数$T$。

5) 输出最优匹配序列$U(C1, C2)$，同时输出匹配代价$χ$。$C1$、$C2$分别为两个目标轮廓匹配好的特征点序列。

$ \chi = \sum\limits_{t = 1}^n {{\zeta _t}} $

(18)

式中，${\zeta _t}$为特征点$C{1_t}$与特征点$C{2_t}$的描述子的卡方距离，$n$为两个轮廓匹配好的特征点数。

第1节中的粗匹配采用低维度的描述子，缺乏对轮廓曲线的描述，因此需要进一步对子曲线段进行精确匹配，以减少匹配误差。基于投影面积的特征点粗匹配较好的确定了形状仿射变换前后特征点的对应，从而使得按照特征点划分的子曲线段保持对应。文献[18]在对仿射曲线的匹配中采用小波仿射不变描述子得到较好的效果。本文在子曲线段匹配过程中，也采用小波仿射不变函数构造不变量，通过对不变量的相似度度量，达到对子曲线段的精确匹配。

设曲线参数化为点集$\mathit{\boldsymbol{C}} = \left\{ {x\left(t \right), y\left(t \right)} \right\}$，其仿射变换后的点集表示为$\mathit{\boldsymbol{C}}{\rm{'}} = \left\{ {\tilde x\left(t \right), \tilde y\left(t \right)} \right\}$，则点集满足关系

$ \left[ {\begin{array}{*{20}{c}} {\tilde x\left( t \right)}\\ {\tilde y\left( t \right)} \end{array}} \right] = {\mathit{\boldsymbol{A}}_{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {x\left( t \right)}\\ {y\left( t \right)} \end{array}} \right] + \mathit{\boldsymbol{b}} $

(19)

式中，${\mathit{\boldsymbol{A}}_{\rm{T}}}$为变换矩阵，$\mathit{\boldsymbol{b}}$为平移矩阵，若$I\left(t \right)$和${\tilde I\left(t \right)}$分别表示对应曲线上的仿射不变函数，则其具有关系

$ \mathit{\boldsymbol{\tilde I}}\left( t \right) = \mathit{\boldsymbol{I}}\left( t \right){\mathit{\boldsymbol{J}}^w} $

(20)

式中，$\mathit{\boldsymbol{J}} = {\rm{det}}\left(\mathit{\boldsymbol{A}} \right)$，$w$为函数的权值，若$w$=0，则$I$为绝对仿射不变函数。设${W_i}x\left(t \right), {W_j}y\left(t \right)$为曲线$C$在尺度$i$和$j$上的小波系数($i \ne j$)，对式(19)进行小波变换，可得${W_i}{a_0} = {W_j}{b_0} = 0$，则其仿射不变函数为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{i,j}}\left( t \right) = \left| {\begin{array}{*{20}{c}} {{W_i}x\left( t \right)}&{{W_i}y\left( t \right)}\\ {{W_j}x\left( t \right)}&{{W_j}y\left( t \right)} \end{array}} \right| = }\\ {{W_i}x\left( t \right){W_j}y\left( t \right) - {W_i}y\left( t \right){W_j}x\left( t \right)} \end{array} $

(21)

根据$\mathit{\boldsymbol{\tilde I}}\left(t \right)$和$\mathit{\boldsymbol{I}}\left(t \right)$存在的关系，又可得

$ \begin{array}{*{20}{c}} {{{\tilde I}_{i,j}}\left( t \right) = {W_i}\tilde x\left( t \right){W_j}\tilde y\left( t \right) - {W_i}\tilde y\left( t \right){W_j}\tilde x\left( t \right) = }\\ {\det \left( \mathit{\boldsymbol{A}} \right){\mathit{\boldsymbol{I}}_{i,j}}\left( t \right)} \end{array} $

(22)

由式(21)定义的函数，利用正规化消除权值的影响，得到的绝对仿射不变函数

$ {f_{a,b,c,d}}\left( t \right) = \frac{{{{\mathit{\boldsymbol{\tilde I}}}_{a,b}}}}{{{{\mathit{\boldsymbol{\tilde I}}}_{c,d}}}} = \frac{{{\mathit{\boldsymbol{I}}_{a,b}}}}{{{\mathit{\boldsymbol{I}}_{c,d}}}} $

(23)

式中，$a, b, c, d$为不同的尺度。

小波系数容易受到噪声的干扰，而且尺度越高，鲁棒性越低，因此，本文仅提取尺度为1到4的小波系数。对于目标轮廓$O$，设通过特征点划分得到的子曲线段集为$\{ C_O^i|i = 1, 2, \ldots, n\} $，$n$为子曲线段数，本文关于$O$的小波局部不变量的构造过程如下：

1) 在轮廓库$O$的第$i$段子曲线$C_O^i$上等间隔进行采样，采样点集表示为$\mathit{\boldsymbol{p}}_O^i = ({p_1}, {p_2}, \ldots, {p_m})$，$m$为采样点数。

2) 对点集$\mathit{\boldsymbol{p}}_O^i$进行小波分解，得到的小波系数表示为$\{ ({W_j}x\left(t \right), {W_j}y\left(t \right))\left| {j = 1, 2, 3, 4} \right.\} $，$j$为小波分解的尺度。

3) 利用式(23)对点集$p_O^i$中的采样点构造仿射不变量，并按照顺序进行存储，即

$ \mathit{\boldsymbol{I}}_O^i = {\left[ {f_O^i\left( 1 \right) \cdots f_O^i\left( t \right) \cdots f_O^i\left( m \right)} \right]^{\rm{T}}} $

(24)

式中，$f_O^i\left(t \right) = f_{1, 2, 3, 4}^i\left(t \right), m = {\rm{max}}\{ {m_j}\} $。

4) 轮廓$O$的识别向量表示为

$ {\mathit{\boldsymbol{I}}_O} = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{I}}_O^1}&{\mathit{\boldsymbol{I}}_O^2}& \cdots &{\mathit{\boldsymbol{I}}_O^n} \end{array}} \right] $

(25)

设轮廓曲线$O$和$R$，由本文提出的小波局部不变量构造的向量分别为${\mathit{\boldsymbol{I}}_O}$和${\mathit{\boldsymbol{I}}_R}$，对这两条曲线的差异度定义为

$ \xi \left( {{I_O},{I_R}} \right) = \sum\limits_{t = 1}^m {\left\| {{f_O}\left( t \right) - {f_R}\left( t \right)} \right\|} $

(26)

式中，$m$为子曲线上采样点的个数。

在粗匹配阶段，本文采用CSS角点作为备选特征点，采用1.1节方法提取特征点的QSPD描述子，使用400次迭代的蚁群算法对描述子匹配，选取匹配误差较小的50%的特征点对目标轮廓进行划分，设有$n$个匹配特征点，则可将闭合轮廓划分为$n$段轮廓子曲线。图 4为camel的50%特征点匹配图，在图像发生仿射变换的条件下，以这些匹配的特征点划分目标轮廓段，确保目标的子曲线段是仿射对应的。

在目标轮廓被划分后，接着就是对原图$O$和仿射图$O′$子曲线段的匹配，选取粗匹配阶段匹配代价最小的5%的目标进行对应曲线匹配。

图 4目标的子曲线段如图 5所示。图 5(a)(b)分别为$O$和$O′$由特征点所划分出的轮廓段。

对于子曲线段，首先进行等间隔采样，采样间隔为$L/m$，$L$为曲线段的长度，$m$为采样点数，这里取值为40。然后通过小波分解得到的小波系数来构造不变量。小波变换为多尺度，为了降低噪声的干扰，本文选取的尺度为1~4。图 6(a)给出了对$S$中的子曲线段2构造的小波不变量曲线，图 6(b)为对应$S′$中的子曲线段2的小波不变量曲线，其中横轴为维度，纵轴为小波不变量。从曲线中可以看出，$I$和$I′$在某些局部上存在一些区别，这主要是由于曲线数字化过程中所生的影响，但整体相似。

通过对应曲线的匹配，可以得出最后的匹配代价

$ \lambda = \chi \times \left( {{\xi _1} + {\xi _2} + \cdots + {\xi _n}} \right)/n $

(27)

式中，$χ$为算法粗匹配代价，${\xi _k}$为相对应的第$k$段曲线差异度，总共$n$段匹配曲线。

实验所采用的计算机环境为Inetl(R) Core2.66 GHz CPU，4.0 GB内存，编程环境为Matlab2014a。

式(6)中投影等间隔划分数目$N$的取值决定了描述子的维数，这对描述子性能有很大的影响。表 1为描述子的平均计算耗时对比，可以看出QSPD的计算耗时远小于SPD^[15]。

表 2为不同维数情况下的几种方案在MPEG-7中的匹配正确率对比，其中SPD与QSPD为使用DP的描述子匹配方案，SPD-RL为使用松弛标签技术的SPD匹配方案，QSPD-ACO₄₀₀与QSPD-ACO₁₀₀₀分别为使用400与1 000次迭代的蚁群算法方案。结果表明，当划分数目达到一定维度时(表 2中该阈值为10)，描述子性能对划分维数不再敏感。而描述子维度与蚁群算法的迭代次数越高，匹配结果越好。

图 7为部分仿射图像特征点的匹配效果。表 3为几种方案的平均匹配耗时对比。结果表明，描述子的维数和蚁群算法的迭代次数越高，计算耗时也越大。为了减少计算耗时，同时保证匹配正确率，本文使用10维数的QSPD描述子与400次迭代的蚁群算法作为第1阶段的粗匹配算法，并且得到匹配代价cost1。对于匹配代价小于5%的图片，将进行第2阶段的精匹配。

表 4为本文算法两个阶段匹配总的耗时与SPD-RL算法的对比，结果表明本文方案的匹配耗时低于SPD-RL。

本文使用MPEG7数据库原始图像进行测试，该数据库中部分图像如图 8所示。对每幅原始图像，按照仿射模型生成5幅不同的仿射图，其对应的仿射变换矩阵分别为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{T}}_1} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} \cos \left( {0.75{\rm{ \mathsf{ π} }}} \right)\\ \sin \left( {0.75{\rm{ \mathsf{ π} }}} \right) \end{array}&\begin{array}{l} - \sin \left( {0.75{\rm{ \mathsf{ π} }}} \right)\\ \cos \left( {0.75{\rm{ \mathsf{ π} }}} \right) \end{array} \end{array}} \right]}\\ {{\mathit{\boldsymbol{T}}_2} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 1\\ 0 \end{array}&\begin{array}{l} 0.5\\ 2 \end{array} \end{array}} \right]}\\ {{\mathit{\boldsymbol{T}}_3} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 0.6\cos \left( {0.8{\rm{ \mathsf{ π} }}} \right)\\ \sin \left( {0.8{\rm{ \mathsf{ π} }}} \right) \end{array}&\begin{array}{l} - 0.6\sin \left( {0.8{\rm{ \mathsf{ π} }}} \right)\\ 2\cos \left( {0.8{\rm{ \mathsf{ π} }}} \right) \end{array} \end{array}} \right]}\\ {{\mathit{\boldsymbol{T}}_4} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} - 1\\ 0 \end{array}&\begin{array}{l} 1\\ 1 \end{array} \end{array}} \right]}\\ {{\mathit{\boldsymbol{T}}_5} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} 1.5\\ 0 \end{array}&\begin{array}{l} 2\\ - 1 \end{array} \end{array}} \right]} \end{array} $

(28)

式中，${\mathit{\boldsymbol{T}}_{\rm{1}}}$为相似变换(仅发生旋转)，其余的为一般的仿射变换。图 9为camel对应的5个仿射图。由原始图像及其仿射图构成本文测试的图像库。

仿射图像库采用3.1节中生成的图像库。对于算法匹配性能的好坏，通过准确率和查全率两个指标进行评估。以每幅图像为匹配对象进行测试，图 10给出了本文算法与文献[5]、文献[8]、文献[9]、文献[15]、Shape Contex^[19]的查准率—查全率曲线对比图。

综合考虑准确率和查全率，可以用综合评价指标(F-Measure)评价它们，即

$ F = \frac{{\left( {{\alpha ^2} + 1} \right)p \cdot r}}{{{\alpha ^2}\left( {p + r} \right)}} $

(29)

式中，$p$为准确率，$r$为查全率。当$F$较高时则能说明实验方法比较有效。表 5为参数$α$为1时图 10的综合评价指标。

结果表明，在本文创建的仿射图像库中，本文算法的表现优于其他几种算法，说明本文算法能够较好地对仿射形变的目标进行识别。本文算法的描述子具备仿射不变性，第2步匹配方案去除了描述子忽视对轮廓曲线的描述带来的误差，这使得本文算法在仿射目标匹配上表现良好。

本文算法与文献[8]、文献[9]、文献[15]的抗噪性对比实验如图 11所示，1~10分别表示信噪比SNR从55平均递减至15。从图中可以看出本文方法虽然不如文献[15]，但比文献[8]和文献[9]要好。因为本文算法计算描述子时利用轮廓质心，使得描述子具有一定的抗噪性，并且在精匹配中采用的小波变换对噪声具有一定的鲁棒性，所以本文算法在匹配有噪声干扰的形状时相对稳定。本文算法的抗噪性不如文献[15]的原因是，本文算法的计算是基于目标轮廓，更易受噪声干扰，但本文算法又利用了区域信息，因此本文算法仍具有的优良的抗噪性。

本文提出了一种结合轮廓与形状特征的仿射形状匹配算法，该算法的主要贡献在于利用轮廓快速地计算投影面积描述子，并在精匹配阶段采用小波仿射不变量方法对目标轮廓子曲线进行匹配，弥补了基于投影面积比描述子方法缺乏对局部轮廓曲线描述的缺点。本文方法解决了仿射描述子计算速度慢、描述能力不足的问题，并且具有一定的抗噪性。QSPD描述子的严格仿射不变性保证了本文方法对仿射变换形状的适用性。由于QSPD的计算是基于全局的形状信息，因此本文算法并不适用于形状变化大的目标。后续研究将探索如何有效地将本文算法应用到基于轮廓的自然图像和3维图像检索中。