WO2015085757A1

WO2015085757A1 - 一种用于图像匹配、识别、检索的频域形状描述方法

Info

Publication number: WO2015085757A1
Application number: PCT/CN2014/081336
Authority: WO
Inventors: 杨夙
Original assignee: Fudan University
Current assignee: Fudan University
Priority date: 2013-12-14
Filing date: 2014-07-01
Publication date: 2015-06-18
Anticipated expiration: 2016-06-14
Also published as: CN103679200A

Abstract

一种用于图像匹配、识别、检索的频域形状特征描述方法，属于模式识别、图像处理、计算机视觉技术领域。在图像区域划分为若干子区域，分别以各个特征点作为参考点统计其它特征点落在各个子区域的个数，得到直方图，计算直方图的功率谱作为一种形状描述，计算两幅图像的形状描述子之间的相似度进而得到图像相似度，从而实现图像匹配、识别、检索。实验表明，形状上下文谱特征具有较好的旋转和伸缩不变性、对于噪声和变形的鲁棒性、且计算开销较小。

Description

说明书

一种用于图像匹配、识别、检索的频域形状描述方法技术领域

本发明属于模式识别、图像处理、计算机视觉技术领域，具体涉及一种图像形状特征提取方法，可以用于图像匹配、识别、检索。背景技术

图像包括形状、纹理、颜色信息，其中，形状是图像识别、检索所依赖的主要信息，形状特征描述对于图像识别、检索非常重要，【S. Belogie, J. Malik, J. Puzicha： " Shape matching and object recognition using shape contexts IEEE Transactions on Pattern Analysis and Machine Intelligence" , Volume 24, pp. 509-52, 2002】论文中提出了一种称为形状上下文（Shape Contexts) 的形状描述方法，首先提取图像特征点（边缘点），然后以每个特征点作为参考点构造一个网格，统计网格各个区域中包含的特征点数，得到直方图形式的形状描述子。形状上下文是一种对于噪声干扰较为鲁棒的特征，但是形状上下文特征不具有旋转不变性，需要利用特征点所在位置的图像灰度域的切向量作为参考方向以获得旋转不变性，但是这种依赖外在条件的旋转不变性不够稳健，尤其对于二值图像而言，无法获得可以作为参考方向的切向量。因此，发明人在已获得授权的中国发明专利【杨夙：一种通用的用于符号识别的特征描述方法，发明专利，授权时间： 2008年 2月 6日，授权国别：中国，授权号： 200410016733.0】禾口论文【S. Yang: "Symbol recognition via statistical integration of pixel-level constraint histograms： A new descriptor" , IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol. 27, No. 2, pp. 278-281, 2005】中提出了点约束直方图特征，该特征具有内在的旋转和伸缩不变性，且对于噪声和变形较为鲁棒，但是缺点是计算复杂度较高，大约为 0(N3)。发明内容

本发明的目的在于提出一种能以较低的计算开销，获得较好的旋转与伸缩不变性、对噪声和变形的鲁棒性的用于图像匹配、识别、检索的频域形状描述方法，即图像形状特征提取方法。

本发明提出的用于图像匹配、识别、检索的频域形状描述方法，包括以下计算步骤：

(a) 以图像特征点的统计量作为图像形状描述的基本元素；假设一幅图像有 K个特征点，记作 Ρ={Ρ1，Ρ2，一,ΡΚ} , 任选一个特征点 PkeP作为参考点，对其它特征点的空间分布进行统计，得到一个相应的直方图，记作 h(Pk) _; 分别以各个特征点作为参考点，则对应每个特征点分别得到一个直方图，共得到 K个直方图（Pk) k=l, 2, ·'·，Κ}_;

(b)对每个特征点对应的直方图求傅里叶变换；设 h(Pk)的傅里叶变换为 F(Pk)，对矩阵 F(Pk)中的每个元素进行函数 f(. )定义的数学变换，得到 f (F (Pk) ) ,将 { f(F(Pk)) k=l，2，"，，K}作为 K个特征点分别对应的形状描述子;

(c) 以所有特征点的形状描述子 [f(F(Pl))， f(F(P2)), -, f(F(PK))]作为输入进行函数 g(. )定义的数学变换，得到 g(f (F(P1))， f(F(P2)), ···, f(F(PK)))， g(f(F(Pl))， f(F(P2) ), ···, f(F(PK)))是最终用于图像匹配、识别、检索的形状描述子。

上面所述的形状描述方法的计算步骤（a) 、 (b) 、 (c) 各有若干种实现方案。上面所述的形状描述方法的计算步骤（a) 的两种实现方案在下面分别用（al) 和 2) 表示，计算步骤（b) 的两种实现方案在下面分别用（bl) 和（b2) 表示，计算步骤（c) 的两种实现方案在下面分别用（cl) 、 (c2) 表示，具体计算方法如下：

(al)以参考点为中心，将图像最小外接圆所在的空间划分为 MxN的网格，计算落入网格每个区间的特征点的个数得到直方图， M和 N都是自然数；

(a2) 以参考点为中心，将图像最小外接圆所在的空间划分为 M个同心圆环形成的区域，计算落入每个圆环的特征点的个数得到直方图， M是自然数；

(bl) F(Pk)的数学变换 f (. )定义为矩阵 F(Pk)中每个元素的模值的 W次方，设 Fij(Pk)表示矩阵 F(Pk)的第 i行、第 j列的元素，贝 IJ

f (Fij(Pk)) = |Fij(Pk) |W， W是自然数;

(b2) F(Pk)的数学变换 f (. )定义为矩阵 F(Pk)中每个元素的模值的对数值，设 Fij(Pk)表示矩阵 F(Pk)的第 i行、第 j列的元素，贝 IJ

f (Fij(Pk))=log(|Fij(Pk) )； (cl) 以所有特征点的形状描述子 [f(F(Pl))， f(F(P2)), -, f (F (PK) ) ] 作为输入进行的函数 g(.)定义的数学变换为：输出等于输入， g(f(F(Pl))， f(F(P2) )，···， f (F(PK))) = [f (F(P1))， f(F(P2) )，···， f (F (PK) ) ]；

(c2) 以所有特征点的形状描述子 [f(F(Pl))， f(F(P2)), -, f (F (PK) ) ] 作为输入进行的函数 g(.)定义的数学变换为：将所有特征点的对应形状描述子进行矩阵加法， g(f (F (PI)), f (F(P2))，…， f (F(PK)))=f (F(Pl))@f (F(P2))，…， 0f (F (PK) ) , 这里运算符 "@"表示参与运算的各个矩阵的对应位置的元素相加， f(F(Pl))@ f(F(P2))，一， @f(F(PK))表示所有 K个矩阵的下标相同的元素相加。

本发明提出了形状上下文谱特征，点约束直方图可以近似看作形状上下文的自相关，因此其傅里叶变换相当于形状上下文的功率谱，形状上下文谱特征与点约束直方图特征相比计算开销较小，同时继承了点约束直方图特征的旋转和伸缩不变性、以及对噪声和变形的鲁棒性。因此，本发明方法克服了形状上下文和点约束直方图的缺点，具有广泛的应用前景。附图说明

图 1为图像识别系统的组成框图。

具体实施方式

一个图像识别系统通常由以下几个环节组成，图像采集、预处理、特征提取、相似度计算、分类，图像识别系统的目标是从图像数据库中返回与输入图像最相似的图像，整个图像识别系统的组成见图 1。这里，图像采集可以通过相机、扫描仪等各种能够完成物理成像的传感设备完成。

实施例 1:

步骤 1: 对一幅输入图像提取特征点，并计算各个特征点的形状描述子，

^P={P1,P2,〜，PK n{f(F(Pk)) |k=l， 2，···，Κ}分别表示所得到的特征点及其对应的形状描述子，形状描述子的计算步骤如下：

(a)任选一个特征点 PkeP作为参考点，对其它特征点的空间分布进行统计，得到一个相应的直方图，记作 h(Pk); 这里，直方图的具体计算方法如下：以参考点 Pk为中心，将图像最小外接圆所在的空间划分为 MxN的网格，计算落入网格每个区间的特征点的个数得到直方图， M和 N都是自然数；分别以各个特征点作为参考点，则对应每个特征点分别得到一个直方图，共得到 K个直方图 (Pk) k=l, 2, -,K}_;

(b)对每个特征点对应的直方图求傅里叶变换，设 h(Pk)的傅里叶变换为 F(Pk)，对矩阵 F(Pk)中的每个元素进行函数 f(. )定义的数学变换， F(Pk)的数学变换 f (. )定义为矩阵 F(Pk)中每个元素的模值的 W次方，设 Fi j (Pk)表示矩阵 F(Pk)的第 i行、第 j列的元素，则 f (Fij(Pk)) = |Fij(Pk) |W， W=2_; 将 {f(F(Pk)) k=l，2， ···，! }作为 K个特征点分别对应的形状描述子;

步骤 2: 对图像库中任选的一幅图像提取特征点，基于步骤 1所述的形状描述子计算方法计算各个特征点的形状描述子，令9=¼1，92，〜，(^}和

{ f(F(Ql)) |1 = 1，2，···， L}分别表示所得到的特征点及其对应的形状描述子; 步骤 3: 计算输入图像和图像库图像的各形状描述子之间的相似度，记作 {dkl=d(f (F(Pk)), f (F(Q1))) k=l, 2，···，Κ; 1 = 1, 2，···， L}，这里采用内积作为相似度度量，即分别将矩阵 f(F(Pk))和 f(F(Ql))拉直为向量，然后求两个向量的内积；

步骤 4: 按照最近邻原则对点集合 Ρ={Ρ1，Ρ2， ·'·，ΡΚ}和 Q={Q1, Q2, ···, QL} 进行匹配，计算方法如下：进行 m=min{K，L}次迭代，每次迭代得到一个匹配的点对，每次迭代的具体计算步骤如下：（a ) 找到集合0={(^1|1^=1，

2， ···，!； 1=1, 2，…， L}中的最大元素 dst; (b) 将 Ps Qt作为一个匹配对记录到集合 P Q中；（c) 令 dsl=- 00禾卩 dkt=- 00 : k=l, 2，…， K且 1 = 1， 2，…， L; 步骤 5: 采用发明人提出的点集合匹配的校正方法对步骤 4得到的初始的点集合匹配关系 P Q={Pl Ql，P2 Q2，…， Pm Qm}进行校正，得到经过校正的点集合匹配关系 {Pl Ql, P2 Q2, ···, Pn Qn} , n≤m; 令

{(XI, Yl) (xl, yl), (X2, Y2) (x2,y2), ···， (Xn， Yn) ^ (xn， yn) }表示匹配点对 {Pl Ql, P2 Q2, ···, Pn Qn}的坐标值；

步骤 6: 基于步骤 5得到的经过校正的点集合之间的匹配关系计算投影变换，步骤如下： (a)

；这里， β中的 4个

：数定义为： s是尺度缩放比例， α是旋转角度， δχ和 Sy分别是 X和 Υ轴的平

(b) 求 Αβ=1^ 最小二乘解，得到 =^了）-1八；

步骤 7:两幅图像之间相似度定义为 Ρ={Ρ1， Ρ2,…， ΡΚ} 和 Q={Q1， Q2, QL} 之间的相似度，计算步骤如下：

(a) 利用公式」对¼1，92，"，，91}进行变换，设 Qj经过变换得到的坐标值记为 T(Qj)， j=l, 2, ···， L;

(b) 点集合 Ρ={Ρ1，Ρ2， "·，ΡΚ} 和 Q={Q1，Q2，〜，QL}之间的相似度定义为

S(P,Q) = min{— ^ EimHdiP, , T(Qj ))}),-∑ £(min{rf(^. , T(Q

E(d)

上式中：函数， t'是预先设定的一个阈值，这里 t'=30;

d(Pi，T(Qj))表示 Pi和 T(Qj)两点之间的欧几里得距离， i = l，2， -,Κ, j=l，2，…山

步骤 8: 设图像库中有 S幅预存的图像，重复执行步骤 2至步骤 7，分别计算每幅图像与输入图像的相似度，根据最近邻分类原则，按照相似度从大到小的顺序输出图像库中排在前 T位的图像作为识别或者检索的结果， T≤S。

实施例 1中的步骤 1和步骤 2采用了发明人提出的一种图像特征点提取方法实现图像预处理，这里提取二值图像的骨架点作为特征点，具体计算方法如下：

步骤 1: 采用 7x7的高斯核滤波器对图像进行平滑，这里高斯核参数选择 σ=2_;

步骤 2: 采用发明人提出的二值化方法对将图像各个像素点二值化为 1和 0 ， 1和 0分别表示前景点和背景点；

步骤 3: 对二值图像进行边缘点检测，这里边缘点定义为两个相邻点的像素值不相等的点；

步骤 4: 对于每个边缘点，找到包含这个边缘点的像素值连续为 1纵向和横向直线段，令纵向和横向直线段中较短的一个的中点为骨架点；

步骤 5: 删除每个骨架点的一定半径内包含的其它骨架点，这里半径值取 1，如果一个骨架点已经被删除，在后续步骤中将不再扫描它。

上面所述的图像特征点提取方法的步骤 2中所述的二值化方法的具体计算步骤如下：

步骤 21: 假设图像有 n个像素点，将图像所有像素点的灰度值进行排序得到 cl≤c2≤,…，≤cn_;

步骤 22: 令 1 = 1且 1=(^，这里 i = l，2，"'，n; 令 t=cl;

步骤 23: 对于 i = l， 2， ···，!!：计算点（xi， yi)到某个直线的距离 di，这里所指的直线由点（xl， yl)和点（xl， yi)确定；

步骤 24: 令 j=l; 令 l=min{i|ci>0}; 令 tj=cl; 令 rj=(n- I)/n;

I = argmaxiii, I i = 1,2,...,n}

步骤 25: 对于 j=2，3: 令 '' ， tj=cl, rj=(n- I)/n;

步骤 26:对于 j=3， 2, 1：如果 rj大于某个预先设定的阈值 r，这里取 r=2%，则令 t=tj;

步骤 27: 对于 i = l，2，"'，n: 如果 ci≥t，令 bi = l; 否则，令 bi=0; 这里 bl， b2, -,bn表示图像二值化后各点的像素值。

实施例 1中的步骤 5采用发明人提出的点集合匹配的校正方法，具体计算步骤如下：

(1) 邻接矩阵计算；

(2) 图论中最大团问题近似求解；

上面所述的点集合匹配的校正方法中的步骤 1中的邻接矩阵计算的步骤如下：

(11) 假设两个点集合？=^1，？2，〜，？111} 和9=½1，92，"，，9111} 之间初始的匹配关系为 P Q= {Pl Ql, P2 Q2, ···, Pm Qm} , 将邻接矩阵初始化为 C= {ci j=0 I i, j = l, 2, ··· m} ;

(12)计算1?= ^=(1(?1 _]')/(1(91 9_]') |i j=l,2, -,m; i≠j}，这里 d(Pi Pj)表示点 Pi和点 Pj之间的欧几里得距离、 d(Qi Qj) 表示点 Qi和点 Qj之间的欧几里得距离；

(13M {rij=d(Pi Pj)/d(Qi Qj) i, j = l, 2, -,m; i≠j} 按照从小到大的次序排序得到 R1≤R2≤ ≤Rm(m-l); 将 Rk 在 R中的位置记录为

S[k]e{(i, j) |i, j=l,2, -,m; i≠j} k=l 2 ··· m (m- 1)

(14) 将 Rl≤R2≤ ≤Rm(m-l) 分段，分段方法如下：如果存在 n-1 段，其边界对应的下标为 1[1]=1 1[2], 1[3], ··· I[n-1], I [n] =m (m-1)且满足条 RI[i]/RI[i+l]>tARI[i]/RI[i+l]+l≤t, 这里 i = l 2 n- 1 t 是一个接近 1的阈值，则以 1[1]=1 1[2], 1[3], ··· I[n-1], I[n]为边界进行分段;

(15) 从上述步骤 14得到的对 Rl≤R2≤ ≤Rm(m-l) 的 n-1个分段中找到最长的一段，其在 n-1个分段中的对应次序为

* = am max{/r + 11 - /Γ11 i = 1,2,.." n-l) _ία ^

- ，提取该段在 R中的对应下标 {S[k] |k=I[i*]

I[i*]+1, ···, I[i*+1]};

(16) 令 {cS[k]=l |k=I[i*] I[i*]+1 ··· I[i*+1]} ;

上面所述的点集合匹配的校正方法中的步骤 2中的图论中最大团问题近似求解的计算步骤如下：

(21) 将一个完全图的所有节点的下标的集合初始化为 ©={l 2 m}_;将噪声节点对应的下标的集合初始化为空集合 Ψ=Φ，将噪声节点之外的剩余节点的集合初始化为 Ω=Θ

(22) 设置计数器 ' 用于记录有多少节点和节点 i相连接；

(23) 对于 i = l 2， ···，!!!：如果 vi=0，将节点 i加入噪声节点的集合，即 Ψυ{ί}→ψ; 同时将节点从剩余节点的集合 Ω中删除，即 Ω-{ί}→Ω;

(24) 如果 vi=vj对任意的 i≠j都成立，这里 ieQ 且 jeQ，则将 P Q中下标为 Ω所包含元素的对应匹配点对输出并退出；否则转到步骤（25);

(25) 找到0中具有最小连接边数的节点^' ^{= ¾1¾1}^¹^' ¹^^_; 令 Ψυ{ί'}→Ψ 且 Ω- {ί'}→Ω; 找到与节点 i'连接的节点，即「 = W =i^i'2〜_m}_; 对于 j_er，令 vj- l→vj; 转到步骤（24)。

实施例 2:

步骤 1: 与实施例 1的步骤 1相同；

步骤 2: 将步骤 1得到的所有特征点 P={P1，P2，〜，PK}对应的形状描述子 [f(F(Pl))， f(F(P2) ), ···, f(F(PK))]进行矩阵加法得到

G(P)=f (F(Pl))@f (F(P2))，…， @f(F(PK))，这里运算符 "㊉"表示参与运算的各个矩阵的对应位置的元素相加， G (P)表示所有 K个矩阵 f (F (P1) )，

f(F(P2)), -, f(F(PK))的下标相同的元素相加， G(P)是最终用于图像匹配、识别、检索的形状描述子；

步骤 3: 与实施例 1的步骤 2相同；

步骤 4: 将步骤 3得到的所有特征点 Q={Q1， Q2,…， QL}对应的形状描述子 f(F (QD), f(F (Q2) ), ···, f(F (QL) )进行矩阵加法得到

G(Q)=f (F(Ql))0f (F(Q2)), ···, 0f (F(QK))_;

步骤 5: 分别将矩阵 G(P)和 G(Q)拉直为向量，然后求两个向量的内积作为输入图像与图像库图像的相似度

步骤 6: 设图像库中有 S幅预存的图像，重复执行步骤 3至步骤 5，分别计算每幅图像与输入图像的相似度，根据最近邻分类原则，按照相似度从大到小的顺序输出图像库中排在前 T位的图像作为识别或者检索的结果， T≤S。

基于实施例 1和实施例 2描述的方法设计了符号识别程序，并对 GREC2003 图像库（http：〃 www. iapr- tclO. org) 进行了实验，共测试了 6900幅图像。另外，用其它两种特征提取方法代替实施 1采用的特征提取方法进行了对比实验，这两种特征提取方法分别为【杨夙：一种通用的用于符号识别的特征描述方法，发明专利，授权时间： 2008年 2月 6日，授权国别：中国，授权号： 200410016733.0】禾口【S. Belogie, J. Malik, J. Puzicha: "Shape matching and object recognition using shape contexts, IEEE Transactions on Pattern Analysis and Machine Intelligence" , Volume 24, pp. 509-52， 2002】。实验结果中， PLC和 SC分别代表上述两种特征提取方法， SSC代表实施例 1采用的解决方案， SSC0代表实施例 2采用的解决方案。实验结果如下：表 1: 理想图像的识别率（％) (50种模型；级别 1: 5种符号， 5 幅图像; 级别 2: 20种符号， 20幅图像；级别 3: 50种符号， 50种图像）

表 2: 旋转和伸缩图像的识别率 (%) (50种模型；级别 1: 5种符号， 25 幅图像；级别 2: 20种符号， 100幅图像；级别 3: 50种符号， 250幅图像)

表 3: 变形图像的识别率 (%) (50种模型；级别 1: 5种符号， 25幅图像; 级别 2: 15种符号， 75 幅图像）

级别 2 100 100 100 98. 67 级别 1 100 100 100 100

2

级别 2 100 100 100 98. 67 级别 1 100 100 100 100

3

级别 2 98. 67 100 98. 67 98. 67

表 4 : 噪声干扰下图像识别率 (%) (级别 1 : 5种模型， 5种符号， 25幅图像；级别 2 : 20种模型， 20种符号， 100幅图像；级别 3 : 50种模型， 50种符号， 250幅图像）

级别 2 100 100 100 92 级别 3 100 100 100 88 级别 1 100 100 100 96

7 级别 2 97 100 100 85

级别 3 99. 6 100 100 86 级别 1 100 100 100 100

8 级别 2 100 100 100 90

级别 3 100 100 100 94. 8

级别 1 100 100 100 92

9 级别 2 99 100 100 78

级别 3 99. 2 100 100 86 表 5 : 噪声和变形同时出现时的图像识别率（％) (15种模型， 15种符号， 75 幅图像)

/3/: O 980Hil£-S/-S80iAV

Claims

权利要求

1. 一种用于图像匹配、识别、检索的频域形状描述方法，其特征在于包含以下计算步骤：

(b)对每个特征点对应的直方图求傅里叶变换，设 h(Pk)的傅里叶变换为 F(Pk)，对矩阵 F(Pk)中的每个元素进行函数 f(. )定义的数学变换，得到 f (F (Pk) ) ,将 { f(F(Pk)) k=l，2，"，，K}作为 K个特征点分别对应的形状描述子;

2. 根据权利要求 1所述的方法，其特征在于计算步骤（a) 中所述的直方图定义如下：以参考点为中心，将图像最小外接圆所在的空间划分为 MxN的网格，计算落入网格每个区间的特征点的个数得到直方图， M和 N都是自然数。

3. 根据权利要求 1所述的方法，其特征在于计算步骤（a) 中所述的直方图定义如下：以参考点为中心，将图像最小外接圆所在的空间划分为 M个同心圆环形成的区域，计算落入每个圆环的特征点的个数得到直方图， M是自然数。

4. 根据权利要求 1所述的方法，其特征在于计算步骤（b)中所述的 F(Pk) 的数学变换 f (. )定义为矩阵 F (Pk)中每个元素的模值的 W次方，设 Fi j (Pk)表示矩阵 F(Pk)的第 i行、第 j列的元素，则 f(Fij(Pk)) = |Fij(Pk) |W， W是自然数。

5. 根据权利要求 1所述的方法，其特征在于计算步骤（b)中所述的 F(Pk) 的数学变换 f (. )定义为矩阵 F(Pk)中每个元素的模值的对数值，设 Fij (Pk)表示矩阵 F(Pk)的第 i行、第 j列的元素， !lJ f(Fij(Pk))=log(|Fij(Pk) |)₀

6. 根据权利要求 1所述的方法，其特征在于计算步骤（c) 中，所述的以所有特征点的形状描述子 [f(F(Pl))， f(F(P2)), -, f(F(PK))]作为输入进行的函数 g(.)定义的数学变换为：输出等于输入， g(f(F(Pl))， f(F(P2)), -, f (F(PK))) = [f (F(P1))， f(F(P2) )， ···， f(F(PK))]。

7. 根据权利要求 1所述的方法，其特征在于计算步骤（c) 中，所述的以所有特征点的形状描述子 [f(F(Pl))， f(F(P2)), -, f(F(PK))]作为输入进行的函数 g(.)定义的数学变换为：将所有特征点的对应形状描述子进行矩阵加法， g(f(F(Pl))， f (F(P2)), ···, f (F (PK) ) ) =f (F (PI) ) 0f (F (P2) ), ···,

0f (F (PK) ) , 这里运算符 "@"表示参与运算的各个矩阵的对应位置的元素相加， f (F(P1))0 £ (?2))，一，@1" (?1 )表示所有1 个矩阵的下标相同的元素相加。