近年来，商标检索受到广泛关注。如，刘婷等^[1]定义了一种ZM与颜色空间的商标检索方案，ZM是一种良好的全局特征描述符，但采取的颜色空间对商标区域的方向变化不显著，无法准确地获取局部特征。栾咏红等^[2]定义了一种结合SIFT特征与多元词袋的图像检索方案，但该方法只统计了局部特征，图像数量较大时，需要计算的数据量巨大，实时性不足，并且缺乏对全局特征进行相似度度量，算法的鲁棒性不强。Yan等^[3]定义了一种颜色和空间描述符的商标检索方法，该方法能够根据颜色特征较好地完成检索，但对于颜色相似的商标效果不佳，此外，在颜色量化过程中计算量较大。刘剑英^[4]提出了一种角点描述与区域特征的商标检索算法，该方法通过双级匹配技术，提高检索精确度，但该算法主要依赖商标的局部特征，忽略全局特征，易使特征匹配出现偏差，限制了检索准确性。Nigam等^[5]提出基于SIFT与色度-饱和度-明度(Hue-Saturation-Value，HSV)的商标检索算法，该算法综合了形状与颜色特征，提高了商标检索能力，但SIFT提取过程依赖局部区域的梯度方向，易使特征匹配出现误差，在扭曲与颜色差异性较小的情况下效果不佳。

针对上述问题，本文提出一种结合图像全局特征和局部特征的商标图像检索(Trademark Image Retrieval，TIR)有效解决方案。全局特征捕捉整体形状，而局部特征表示商标的细节信息。本算法将ZM作为全局特征，SIFT为局部特征。首先，从商标中提取ZM向量的相异值，并对其进行排序，数据库中的大多数不同图像将被排除并构建候选图像；然后，利用SIFT特性对查询图像与候选图像进行精细匹配，输出检索结果；最后，通过实验对所提算法进行验证。

1 Zernike矩(ZM)

ZM的核心是在极坐标上表示一组Zernike多项式。设强度函数f(r, θ)，其重复长度为q的p阶二维ZM表示如下^[6]：

${Z_{pq}} = \frac{{p + 1}}{\pi }\int_0^1 {\int_{ - \pi }^\pi {f(r, \theta )} } r{\rm{d}}r{\rm{d}}\theta $

(1)

式中：θ为旋转角度；r为半径，$\left| r \right| \leqslant 1$。ZM多项式${V_{pq}}(r, \theta )$描述如下：

${V_{pq}}(r, \theta ) = {R_{pq}}(r){{\rm{e}}^{ - \hat {\rm{j}}q\theta }}$

(2)

式中j为一个虚数，$\hat {\rm{j}} = \sqrt { - 1} $。真实值径向多项式${R_{pq}}(r)$定义如下：

${R_{pq}}(r) = \sum\limits_{k = 0}^{\frac{{p - \left| q \right|}}{2}} {{{( - 1)}^k}} \frac{{(p - k)!}}{{k!(\frac{{p + \left| q \right|}}{2} - k)!(\frac{{p - \left| q \right|}}{2} - k)!}}{r^{p - 2k}}$

(3)

式中$0 \leqslant \left| q \right| \leqslant p$，$p - \left| q \right|$为偶数，$p \geqslant 0$。

由于ZM是在(r, θ)表示的，$\left| r \right| \leqslant 1$。因此，Zernike多项式的求解需在域$(x, y) \in {R^2}$进行线性变换。图像尺寸为N×N的ZM的离散形式表示如下：

${Z_{pq}} = \frac{{p + 1}}{\pi }\sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {f(x, y){V_{pq}}(x, y)} } $

(4)

式中x²+y²≤1。

使用ZM的目的是将图像展开为一系列正交基，因此形状表示的精确度取决于展开中使用的矩的数量^[7]。对于人眼来说，五阶以上的系数太小，无法可靠地测量，可忽略高阶系数。因此，本文计算了前4阶ZM，从而得到最有效、最可靠的全局形状测量结果。表 1给出了前4阶Zernike的多项式。

根据前4阶ZM的值，用Euclidean距离计算任意2幅商标的相异值，定义如下：

$d = \sqrt {\sum\limits_{i = 1, 2, \cdots , 15} {{{({\mathit{\boldsymbol{Z}}_{\rm{S}}}(i) - {\mathit{\boldsymbol{Z}}_{\rm{O}}}(i))}^2}} } $

(5)

式中${\mathit{\boldsymbol{Z}}_{\rm{S}}}(i), {\mathit{\boldsymbol{Z}}_{\rm{O}}}(i)$分别为查询商标与数据库商标中第i个ZM特征向量。根据式(5)得到的相异值选择最合适的候选商标，并进行排序。

2 SIFT算子

SIFT对旋转、缩放、亮度等保持不变性，SIFT可在大量特征库进行快速、准确的匹配，并且容易与其他的特征向量联合使用^[8]。

2.1 尺度空间极值检测

SIFT计算的第一步是利用Gaussian的差分来检测尺度空间中关键点位置。而尺度空间可由不同的Gaussian卷积得到：

$L(x, y, \sigma ) = G(x, y, \sigma )I(x, y)$

(6)

式中：I(x, y)为输入图像；$\sigma $为关键点的尺度，也是Gaussian的标准差，表示模糊度；$L(x, y, \sigma )$为图像的Gaussian尺度空间。$G(x, y, \sigma )$为Gaussian核函数^[9]：

$G(x, y, \sigma ) = \frac{1}{{2\pi {\sigma ^2}}}{e^{ - ({x^2} + {y^2})/2{\sigma ^2}}}$

(7)

$D(x, y, \sigma )$主要用于检测尺度空间中关键点位置，使用2个图像之间的差异来计算：

$D(x, y, \sigma ) = (G(x, y, k\sigma ) - G(x, y, \sigma )) + I(x, y) = L(x, y, k\sigma ) - L(x, y, \sigma )$

(8)

式中k为Gaussian尺度空间的比例因子。

假设一个Gaussian金字塔包含多个组，每一组含有多个层。不同层的$\sigma $不一样，邻近2层间的尺度差异为比例因子k。如果每组有S层，则$k = {2^{1/S}}$。为了检测$D(x, y, \sigma )$的局部极值，需将每一点与同一尺度下的8个邻域点进行对比。如果这个值是所有这些点的最小或最大值，则该点就是一个极值，极值被用作筛选关键点。

2.2 关键点定位

通过上述过程得到的点可能并不是实际的极值点，图 1描述了2D离散空间计算的极值点与连续空间极值点的区别。

为增强关键点的精准性，对得到的点拟合优化^[10]：

$D(X) = D + \frac{{\partial {D^{\rm{T}}}}}{{\partial X}}X + \frac{1}{2}{X^{\rm{T}}}\frac{{{\partial ^2}D}}{{\partial {X^2}}}X$

(9)

式中：$X = {(x, y, \sigma )^{\rm{T}}}$；D为高斯核函数$D(x, y, \sigma )$。

对式(9)求导并等于零，计算每个关键点的Laplacian值，可得到极值点的位置为：

$\hat X = - \frac{{{\partial ^2}{D^{ - 1}}}}{{\partial {X^2}}} \times \frac{{\partial D}}{{\partial X}}$

(10)

对应极值点方程的值可表示为：

$D(\hat X) = D + \frac{1}{2} \times \frac{{\partial {D^{\rm{T}}}}}{{\partial X}}\hat X$

(11)

式中$\hat X = {(x, y, \sigma )^{\rm{T}}}$表示偏移值。当它在$(x, y, \sigma )$上的移位 > 0.5时，说明中心远离了相邻点，需更改点的位置。

采用的拟合函数具有强烈的边缘响应，需消除不稳定点。通常的方法是获取特征点处的Hessian矩阵H^[11]：

$\mathit{\boldsymbol{H}} = \left[ \begin{gathered} {D_{xx}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {D_{xy}} \\ {D_{xy}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {D_{yy}} \\ \end{gathered} \right]$

(12)

式中${D_{xx}}{\kern 1pt} , {D_{xy}}, {D_{yy}}$表示候选点相邻的差分。$\alpha = {\lambda _{\max }}$表示H最大值；$\beta = {\lambda _{\min }}$表示H最小值。

随后，计算矩阵H的迹Tr(H)与行列式Det(H)：

$\left\{ \begin{gathered} Tr(\mathit{\boldsymbol{H}}) = {D_{xx}} = + {D_{yy}} = \alpha + \beta \\ D{\rm{et}}(\mathit{\boldsymbol{H}}) = {D_{xx}} + {D_{xx}} - D_{xy}^2 = \alpha \beta \\ \end{gathered} \right.$

(13)

设$\gamma = \frac{\alpha }{\beta }$表示$\alpha = {\lambda _{\max }}$和$\beta = {\lambda _{\min }}$的比值，则：

$\frac{{Tr{{(\mathit{\boldsymbol{H}})}^2}}}{{D{\rm{et}}(\mathit{\boldsymbol{H}})}} = \frac{{{{(\alpha + \beta )}^2}}}{{\alpha \beta }} = \frac{{{{(\gamma + 1)}^2}}}{\gamma }$

(14)

式(14)中计算的结果与α, β的比值有关。当α, β值相同时，其值最小，并随$\gamma $的增加而变大。因此，为判断主曲率处于${T_\gamma }$下，只需判断以下公式是否满足：

$\frac{{Tr{{(\mathit{\boldsymbol{H}})}^2}}}{{D{\rm{et}}(\mathit{\boldsymbol{H}})}} > \frac{{{{({T_\gamma } + 1)}^2}}}{{{T_\gamma }}}$

(15)

式中${T_\gamma }$代表曲率阈值。

如果满足式(15)，则去除该点，反之保留。通常情况下${T_\gamma }$=10。

2.3 方向分配

为使特征保持旋转不变性，给每一个特征增加一个方向表示。对于每个图像样本L(x, y)，在这个尺度上，梯度大小G(x, y)和方向θ(x, y)表示如下：

$G(x, y) = \sqrt {{{\left( {L(x + 1, y) - L(x - 1, y)} \right)}^2} + {{\left( {L(x, y + 1) - L(x, y - 1)} \right)}^2}} $

(16)

$\theta (x, y) = {\rm{arc}}\;\tan \left( {L(x + 1, y) - L(x - 1, y)} \right)/\left( {L(x, y + 1) - L(x, y - 1)} \right)$

(17)

梯度直方图将0~360°划分为36组，每组为10°。如图 2所示，峰值方向表示点的主方向。为提高鲁棒性，将高于0.8倍的主方向峰值作为某点的候选方向。对此，在相同位置与尺度可形成多个不同方向，从而提高匹配的稳定性。

2.4 关键点描述

通过一组向量来描述一个关键点，使其具有光照、视角不变性等。这个描述符不仅能表示关键点，同时还含有邻域中的其他信息，可有助于提升特征点正确匹配的几率。SIFT是关键点邻域梯度计算的表示。将关键点邻域区域划分不同块。描述符利用4×4的窗口，测量8个方向的梯度，产生4×4×8=128个向量。其计算过程如下：

1) 确定描述符区域：关键点描述符与其所在的尺度有关，因此，对梯度的计算应在Gaussian上完成。将某点邻域分割为d×d(d=4)个子块，每个子块视为一个种子，每个种子有8个方向。

2) 将坐标变换为关键点的方位，如图 3所示。变换后新坐标如下：

$\left( \begin{gathered} {x'} \\ {y'} \\ \end{gathered} \right) = \left( \begin{gathered} \cos \theta {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - \sin \theta \\ \sin \theta {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cos \theta \\ \end{gathered} \right)\left( \begin{gathered} x \\ y \\ \end{gathered} \right)$

(18)

3) 将数据点分配到相应的子块内，将子块内的梯度分派到8个方向上，测算其权重w：

$w = m(a + x, b + y){{\rm{e}}^{ - \frac{{{{(x')}^2} + {{(y')}^2}}}{{2{{(0.5d)}^2}}}}}$

(19)

式中：a和b均为常量；d为子块的长度；m为子块。

4) 测量每个种子8个方向的梯度。图 4为直方图，对于红色点，对0行3列的作用设为dr，对1行3列的作用为1-dr，邻近两列的作用为dc，邻近两行的作用为do，累加在每个方向上的梯度可表示如下：

$weight = wd{r^k}{(1 - dr)^{1 - k}}d{c^m}{(1 - dc)^{1 - m}}d{o^n}{(1 - do)^{1 - n}}$

(20)

式中k, m, n为常量，其值为0或1。

5) 通过计算4×4×8=128个梯度作为关键点的特征值。特征值生成后，为消除光照变化的干扰，对其执行归一化，生成的向量为$\mathit{\boldsymbol{H}} = ({h_1}, {h_2}, \cdots , {h_{128}})$，归一化表示为$L = ({l_1}, {l_2}, \cdots , {l_{128}})$。

6) 对于不均匀光照导致个别方向的梯度巨大，通过设置限值(归一化后取0.2)排除大的梯度，提高特征的判断力。

7) 特征描述向量进行排序。

2.5 SIFT关键点匹配

假设关键点p和q的特征描述符分别是S_p和S_q，则基于Euclidean距离测量可定义如下^[12]：

$d = \sqrt {\sum\limits_{i = 1, 2, \; \cdots , 128} {{{\left[ {{S_p}(i) - {S_q}(i)} \right]}^2}} } $

(21)

3 本文商标检索算法

本文商标检索系统由离线数据库构建部分和在线图像检索部分组成。离线数据库构建部分旨在通过提取数据库中每个图像的特征集来确保高检索效率。在脱机方式下，将特征集及其相应的图像存储在数据库中。因此，当向系统显示查询图像时，系统不必对每个数据库图像执行在线特征提取。其详细的检索过程如下：

Step1：分别提取查询图像A的ZM, SIFT特征。

Step2：对于数据库B中的任意图像，设B的特征集为F_b。

Step3：系统根据查询图像A的特征集与数据库B中存储的图像的特征集之间的ZM特征来测量相似性。

Step4：根据Euclidean距离的大小，按升序对ZM相似度排序，创建候选商标集。

Step5：设图像A的特征为F_a，根据Euclidean距离，对于每个元素K_i，使用最优节点优先(Best Bin First，BBF)搜索算法^[13]计算出k最近邻。假设f₁, f₂是返回的2个最近邻结果。

Step6：取距第一个最近邻点的距离与距第二个最近邻点的距离之比d₁/d₂来确定关键点f₁是否与K_i匹配。

Step7：对所有K_i重复1~6过程，并计算F_b中匹配的关键点n的数量。

Step8：重复4~7，计算查询图像A与候选商标集中匹配的关键点数目。

Step9：对匹配的关键点进行排序，并返回结果，完成检索，算法结束。

4 实验与分析

为了验证所提算法的检索性能，随机选择100幅常见的标准商标，将所选的商标进行一系列的变换，以生成一个含有1 200幅图像的商标库。为体现该算法的先进性，分别选择文献[1]、文献[2]、文献[3]作为对照组。图 5是以“长安商用”为对象，将标准商标进行一系列的变换，生成了9张变换商标。

4.1 评价指标

为准确评价算法性能，引入当前流行的Precision-Recall与F值，Precision表示返回结果中正确商标比例，Recall是指所关注的结果中召回商标的比例，定义如下^[14]：

${\rm{Precision}} = \frac{{{T_{\rm{p}}}}}{{{T_{\rm{r}}}}}\quad {\rm{Recall}} = \frac{{{T_{\rm{p}}}}}{{{T_{\rm{s}}}}} $

(22)

式中：T_p表示相关商标；T_r表示检索结果中的图像总数；T_s表示数据库中相似商标的数量。

F值是Precision-Recall的综合评价，根据式(22)可得F值^[15]：

$F = \frac{{2PR}}{{P + R}}$

(23)

式中：P代表返回结果中正确商标的数量；R代表召回商标的数量。

4.2 实验结果

利用所提算法，在生成的数据库中进行了6次实验，统计每次实验的检索商标数量(见表 2)。从表 2可看出，在6组实验中，平均Precision为91.66%，平均Recall为76.39%。表 3是以“长安”商标为对象的Zernike多项式计算结果，其中，取ρ=0.5, θ=30°。根据前4阶ZM的值，用Euclidean距离计算任意商标的相异值d=0.091 2。表 4是以“美的”商标为对象，Zernike多项式计算的结果，其中ρ=0.8, θ=60°。根据前4阶ZM的值，用Euclidean距离计算任意商标的相异值d=0.076 6。

图 6、图 7为对照文献[1]~文献[3]算法与本文算法进行商标检索结果，为便于标记，记录为实验一、实验二。其中图 6(a)、图 7(a)为检索对象，图 6(b)~(e)、图 7(b)~(e)依次为文献[1]~文献[3]与本文算法返回的前12个结果。图 6(b)中前9张是相关商标，后3张是非相关商标，图 6(c)、(d)中前10张是相关商标，后2张是非相关的。而所提算法返回的12张商标中含有11张相关商标，1张非相关。图 7(b)、(c)返回9张相关商标，3张不相关商标。图 7(d)返回了10张相关商标，2张非相关商标。图 7(e)中返回了11张相关商标，1张非相关。从这4种检索结果中看出，所提算法检索性能最优，准确性高，对各种变换形式下的商标能够有效识别判断。主要原因是所提算法是分别从查询图像中提取ZM和SIFT特征。ZM能够有效描述物体整体形状，具有旋转不变性，SIFT对旋转、缩放、亮度等保持不变性，能够在大量特征库进行快速、准确的匹配。根据ZM特征进行度量，进行排序，形成候选商标集。再利用SIFT对查询图像与候选图像完成更精准匹配，并将结果返回。而文献[1]采取的颜色空间对方向、大小等变化较弱，无法有效准确地捕捉对象的局部特征。文献[2]的方法只统计了局部特征，并且缺乏全局特征进行相似度量，算法的鲁棒性不强。文献[3]的方法在颜色量化中无法准确有效识别颜色特征，空间描述符对噪声和模糊表达能力不强，导致该算法检索效果不太理想。

图 8、图 9为图 6、图 7实验中的Precision-Recall、F统计分析。在实验一中，当Recall为0.5时，文献[1]~文献[3]与本文算法的Precision值分别为0.61, 0.65, 0.72, 0.81。在实验二中，当Recall为0.5时，文献[1]、文献[2]、文献[3]与本文算法的Precision值依次为0.56, 0.62, 0.70, 0.78。从P-R曲线图中可看出，本文算法具有良好的查准率与查全率。所提算法曲线平稳，在相同返回商标数量时，得到的F值最高，表明算法性能最优。

5 结论

本文提出了一种全局特征与局部特征相结合的商标图像检索方案。利用ZM作为全局特征，能够有效捕捉形状本质；选择SIFT作为局部特征，能够准确描述商标的内部细节。通过计算商标中ZM，并对其进行排序，筛选非全局相似商标，形成候选商标集。对于不均匀光照，饱和度变化导致某些方向的梯度值过大情况，通过设置限值并进行归一化处理，可有效消除影响。再利用SIFT特性对查询图像与候选图像进行精细匹配，返回检索商标。实验结果表明，该方法提高了检索性能，获得了较好的Precision-Recall值与F值。