SUSAN边缘响应值灰度化转共生矩阵检索

SUSAN边缘响应值灰度化转共生矩阵检索
熊国萍
【摘要】SUSAN初始边缘响应矩阵元素值对应USAN中像素点总数,该统计值与响应矩阵元素分布均体现图像特征.若将该统计值视为灰度值,则响应矩阵中的元素分布特征可以视为灰度矩阵中的纹理特征.提出SUSAN边缘响应值灰度化转共生矩阵检索算法,即先计算图像的SUSAN边缘响应矩阵,再按映射规则转换为灰度矩阵,然后计算灰度共生矩阵及其各特征描述子,最后进行特征检索.实验显示,该算法的查全率与查准率在检索结果数量达到某临界点之后较为满意,且体现一定的仿射变换、亮度变化等不变性与抗噪鲁棒性.
【期刊名称】《计算机工程与应用》
【年(卷),期】2016(052)005
【总页数】7页(P225-230,235)
【关键词】最小吸收核同值区;吸收核同值区;初始边缘响应;灰度化;灰度共生矩阵;图像检索
【作者】熊国萍
【作者单位】上海医疗器械高等专科学校图文信息中心,上海200093
【正文语种】中文
【中图分类】TP301.6;TP391.3
图像检索技术中需要寻求图像特征表达能力强且时空开销小的算法，而特征描述复杂度与算法时空开销往往成正比。

例如，SIFT算子具备尺度与旋转不变性，且在
保质范围内的仿射失真、3D视角变换、添加噪音与照明变化方面提供鲁棒性匹配能力[1]，但是其128维特征向量描述复杂，相应的检索算法实时性不高；灰度直
方图算法时空复杂度低，但却不具备尺度不变性且对噪音与照明变化敏感，特征表达能力低下。

因此，本文试图探索图像特征描述能力较强同时实时性处理能力较高的图像检索新算法。

在图像特征描述方面，角点描述算子因具备局部不变特性而体现较强的图像特征表征能力，该类算法包括Moravec[2]、Harris[3]、SUSAN（SmallestUnivalue SegmentAssimilating Nucleus）[4]、SIFT（Scale Invariant Feature Transform）[1]、SURF（Speeded Up Robust Features）[5]、FAST[6-7]等。

其中，SUSAN算法较Moravec算法具有更好的边缘检测与抗噪性能，同时相对Harris、SIFT、SURF算法而言因不涉及微分运算，故更为简单高效且其边缘算子的性能不受模版尺寸影响。

FAST算法虽然具有高计算效率、高可重复性特点且适合实时应用场景，但其受噪声与阈值影响较大。

此外，文献[4]表明，SUSAN算子提取精度与抗噪能力较早期的边缘提取算子如Prewitt、Sobel、Canny算子更为优越。

因此，本文选择SUSAN算子并提出SUSAN边缘响应值灰度化转共生矩阵检索算法。

SUSAN算法是一种低层图像处理方法，特别适用于边缘检测（一维特征检测）、角点检测（二维特征检测，包括角点、交叉点等）及保结构去噪[4]。

该算法包括
边缘检测与角点检测算法。

本文仅选用其边缘检测算法。

在SUSAN边缘检测算法中，SUSAN边缘探测器是通过使用圆形模板（有时也称为窗口或内核）来实现的，通常半径为3.4像素、包含37像素。

模板将滑过图像的每个像素点，分别以每个像素点为中心，将模板内部的每个像素点的亮度与核（中心）的亮度作比较，比较公式为[4]：
其中，r0是核在二维图像中的位置，r是模板内其他点的位置，I(r)是像素的亮度，
t是亮度差异阈值，c是比较结果输出。

比较完模板内的每个像素之后，将获得一个累计的总值n作为输出：
该总值n正是USAN中的像素个数，亦即USAN面积。

之后，n用于与固定的阈值 g作比较，而g取值为3nmax/4，其中nmax是n最大取值。

因而初始边缘响应可由如下规则创建：
其中，R(r0)为初始边缘响应。

可见，USAN面积越小，边缘响应越大[4]。

至此该算法已能获得较好的结果，然而对c而言可以使用更具稳定性与判断性的公式：
在描述纹理时，不仅要考虑灰度的分布，还要考虑图像中像素的相对位置[8]。

令Q是定义两个像素彼此相对位置的一个算子，并考虑一幅具有L个可能灰度级的图像 f。

令G为一个矩阵，其元素gij是灰度为zi和zj的像素对出现在 f中由Q所指定的位置处的次数，其中1≤i，j≤L。

按此方法形成的矩阵称为灰度级（或灰度）共生矩阵。

当含义明确时，G简称为一个共生矩阵[8]。

图像中可能的灰度基数决定了矩阵G的大小。

对于一幅8比特图像（256个可能的灰度级），G的大小将为256×256。

为了减小计算负担，经常使用的一种方法是，将灰度量化为几段，以便于管理矩阵G的大小。

例如，在256级灰度的情形下，通过令前32个灰度级等于1，令接下来的32个灰度级等于2，以此类推，可完成这种量化工作。

这将得到一个大小为8×8的矩阵[8]。

满足Q的像素对的总数，等于G的元素之和，因此
是满足Q的一个值为(zi，zj)的点对的概率估计。

这些概率的值域为[0，1]，且它们的和为1：
其中，K是方阵G的行（或列）数。

因为G取决于Q，所以选择一个合适的位置算子并分析G的元素，可以检测灰度纹理模式的存在。

首次提出灰度共生矩阵概念的Haralick等在文献[9]中给出14
种表征G的内容的描述子。

文献[10]详细研究了其中6个纹理特征，认为对比度和熵是最重要的两个特征。

本文在实验中采用了如下一组常用的描述子：
相关：
对比度：
一致性（能量）：
同质性：
熵：
其中：
mr是沿归一化后的G的行计算的均值，mc是沿归一化后的G的列计算的均值。

σr和σc是分别沿行和列计算的标准差（方差的平方根）[8]。

4.1 核心思想
图像的SUSAN初始边缘响应一方面包含图像的边缘信息，可直接用于提取图像边缘或进一步提取角点；另一方面，该响应矩阵中的元素值对应USAN中像素点个数统计值（USAN面积），其元素值与元素分布均体现图像特征，若将该统计值视为灰度值，则响应矩阵中的元素分布特征可视为灰度矩阵中的纹理特征。

本算法核心思想正是将SUSAN初始边缘响应矩阵按一定的映射规则转换为灰度矩阵，再计算该灰度矩阵的灰度共生矩阵及其各特征描述子，以此作为检索的相似性度量条件。

4.2 算法步骤
本算法可分为如下步骤：
步骤1图像预处理。

图像预处理包括图像格式的转换与规格化、图像的增强与去噪等操作。

其中，本文采用公式（12）对256色图像灰度化处理：
步骤2使用SUSAN边缘检测算法获得图像的初始边缘响应矩阵。

因SUSAN边缘检测采用37像素圆形模板，故由公式（1）、（2）可知，USAN 面积n(r0)的值域为[0，37]∩N，其中N表示整数集。

公式（3）中，g取值为
3nmax/4，则当nmax=37时，g=27.75，因而R(r0)的值域下取整后为[0，
27]∩N。

当将7×7的37像素圆形模板在原图像滑动时，因本算法不在外围扩展图像，故
窗口中心未覆盖图像外围3像素宽度大小的边缘，所以对于像素矩阵大小为m×n 的图像而言，最终获得的初始边缘响应矩阵大小为(m-3×2)×(n-3×2)，即(m-
6)×(n-6)。

以像素矩阵大小为16×16的Lena图为例，经SUSAN边缘检测后，将获得如图
1所示10×10的边缘响应矩阵。

步骤3 SUSAN边缘响应矩阵元素值灰度化。

初始边缘响应R(r0)的值域下取整后为[0，27]∩N，而[0，27]∩N⊂[0，255]∩N，即采取37像素模板后的SUSAN初始边缘响应矩阵的元素值属于256阶灰度空间范围，故将该矩阵视为灰度矩阵理论上可行。

实际上，在计算灰度共生矩阵时常将灰度量化为几段[8]，如将256阶灰度空间映
射为8阶灰度空间。

类似地，本步骤将值域[0，27]∩N映射至[0，7]∩N，即将
27阶灰度空间映射为8阶灰度空间（5.1节实验显示，该降阶映射转换将有利于
提高本算法的查全率与查准率）。

具体说明如下：
若令fij表示图像第i行、第 j列像素的灰度，m表示图像的高度，n表示图像的
宽度，则图像的灰度矩阵可记为：
其中，0≤i≤m-1，0≤j≤n-1。

若图像经SUSAN边缘检测后的初始边缘响应R(r0)所对应的矩阵记为：
其中，0≤i≤m-7，0≤j≤n-7，初始边缘响应矩阵A经灰度映射转换所得矩阵记为：其中，0≤i≤m-7，0≤j≤n-7，则初始边缘响应矩阵灰度化映射规则可表示为：
因aij∊[0，27]∩N，所以公式（13）可简化为：
则bij∊[0，7]∩N，即B为8阶灰度空间矩阵。

步骤4计算8阶灰度空间矩阵B的灰度共生矩阵及其各特征描述子。

为减小空间开销，本算法采取矩阵分块计算方法，即先按一定分块大小将矩阵B 分成若干子矩阵：
其中，0≤i≤k-1，0≤j≤l-1，然后分别统计 Bij所对应的灰度共生矩阵Gij及其各特征描述子，最后对Gij累计求和得到父矩阵B的灰度共生矩阵G并对Gij的各特征描述子平均化得到父矩阵B的各特征描述子。

在计算灰度共生矩阵时，需要选取定义两个元素彼此相对位置的算子Q。

若取图像中任意一点(x，y)及偏离它的另一点(x+a，y+b)，Q值即为(a，b)。

本算法仅求取水平方向（0°方向）上的灰度共生矩阵及其各特征描述子，其对应偏离点为(x+5，y)，即Q值为（5，0），其中参数a=5的选取是根据样本图片大小通过实验中调节参数观测最佳效果确定的。

实际应用中，可分别计算0°、45°、90°、135°方向上的灰度共生矩阵的特征描述子再加权组合化处理。

步骤5图像检索。

SUSAN响应值灰度化转共生矩阵G的各特征描述子可作为图像的特征向量。

经步骤1～4分别求出示例图像与待检索图像的各特征描述子作为特征向量，则可以两者的欧式距离作为检索的相似性度量依据。

为简化计算，本算法仅选取G的熵[10]作为图像的特征值，以两图的特征值之差的绝对值作为相似性判定条件，该值越小表明示例图像特征与待检索图像特征越接近。

以VC6.0为开发工具，基于Microsoft Windows 7平台设计实验系统。

实验机器配置如下：处理器：Intel®************************；内存（RAM）：4.00 GB；系统类型：64位操作系统。

实验样本库选用科雷尔（Corel）公司收集整理的Corel图像库，该库涵盖多个主题，每个主题包含100幅图像。

Corel图像库原始图像均以24位真彩色JPEG格
式的图片文件存储，且图片大小均为192 px×128 px或128 px×192 px。

实验
中采用Adobe Photoshop软件将其批处理转换为256色BMP格式文件以满足
实验要求。

实验共分三类，分别介绍如下。

5.1 不同主题相同变换检索实验
该类实验测试本文算法在样本为不同主题、相同变换处理情形下的检索性能。

为便于实验统计与分析，该类实验选取Corel库中sun、flowers、horses等10
类主题的原始图片，其中每类主题各选取10张，共计选取100张图片作为测试样本。

本类实验共分四组，分别测试基于GH（Gray Histogram，灰度直方图）、GLCM（Gray Level Co-occurrence Matrix，灰度共生矩阵）、本文算法（未降
阶映射）与本文算法的图像检索性能。

每组实验选用同一种检索方法，但选取三个不同主题的图片分三次测试。

每次实验时，将首先选取某种检索方法，然后在测试样本中选择某一主题的一张图片作为示例图像作相关检索，分别记录系统返回的图像数量T在不同取值下的相关数与P值（Precision，查准率）、R值（Recall，
查全率）及算法时空开销。

其中，本类实验中的相关图片是指同一主题的图片。

对每组实验的三次观测数据求取均值，记为该检索方法的检索性能。

以选取主题分别为sun、buildings、plane的三组实验为例，其检索效果如图2
所示、观测数据如表1所示。

其中，R值随T值增长曲线如图3所示。

该类实验数据体现如下性能特点：
（1）查准率P与查全率R：①当T取值为[20，100]时，本法的P、R值相对均
为最高，特别是当T=50时，本法的R值达到100%；②当T取值为（0，15）时，
GH法的 R值相对为最高，但其平均增长速度不如本法，如图3所示；③无论T
取何值，GLCM法的P、R值相对均为最低。

（2）时空开销：GH法最小，本文算法最大。

与GH、GLCM法相比，本文算法的特点如下：
（1）R值增长最快，如图3所示，且在T值达到某临界点之后，其P、R值相对
最高，优势明显。

原因分析如下：首先，本法采用SUSAN边缘检测算法获得的初始边缘响应矩阵包含图像的边缘与角点特征，相对GH、GLCM具有更强的特征表征能力；其次，在初始边缘响应矩阵元素统计值转灰度化后，其对应的灰度共生矩阵携带之前的元素相对位置特征。

当T值达到某临界点之后，同类图像的同质特
征才相对集中凸显，利用本文算法表征的同一主题的图像特征总体上明显区别于其他主题的图像特征，因而检索出的相关数明显增加，P、R值相对更高。

（2）时空开销最大。

原因分析如下：本文算法相对需要更多的矩阵相关运算，如计算USAN一维矩阵、USAN一维矩阵转二维矩阵等处理，因而需付出更大的时
空开销。

5.2 不同主题不同变换检索实验
该类实验测试本文算法在样本为不同主题、不同变换处理情形下的检索性能。

该类实验首先从Corel图像库中的plane、flowers等5个不同的主题中各选取2张图片，共计选取10张图片。

然后对原图使用Photoshop分别进行如下处理：
对比度增加50、亮度增加75、高宽缩小一半、顺时针旋转45°、顺时针旋转90°、平移（+50 px，+50 px），自由变换（+15°，+15°）。

再对原图使用Matlab
分别进行如下处理：叠加均值为0方差为0.02的高斯噪声、叠加密度为0.04的
椒盐噪声。

若称原图为源图，称经图像处理之后的图像为源图对应的派生图，则经以上处理后，每张源图产生9张派生图，连带10张源图本身共计获得100张图片，以此作为本类实验的测试样本。

本类实验方法与5.1节实验（不同主题相同变换检索实验）相似，区别在于样本的选择与相关数的界定不同。

本类实验中定义源图与其经图像处理之后的派生图互为相关，相关数指检索结果中以示例图像为源图的派生图（包含源图）的图像数。

同法，可得其观测数据如表2所示。

该类实验数据体现如下性能特点：
（1）查准率P与查全率R：①当T取值为[20，100]时，本法的P、R值相对均为最高；②当T取值为（0，15）时，GH法的P、R值相对均为最高，但是其R 值增长速度不如本文算法；③无论T取何值，GLCM法的P、R值相对均为最低。

（2）时空开销：GH法最小，本文算法最大。

进一步分析检索结果及漏检、误检情况，得出如下特点：
（1）不变特性：本文算法在T=50时，除高斯噪音派生图漏检外，几乎能检索出所有其他各种变换处理派生图，可见本文算法相对GH、GLCM法而言，除高斯噪音外，对各变换均表现其不变特性。

原因分析为：本文算法继承了SUSAN边缘检测算法带来的局部不变性特性。

（2）不相关图像：本类实验定义的相关为源图与其派生图，但同一主题的图片具有一定的相似特性，故本类实验检索出的结果中不相关的图像以同一主题的另一图片及其派生图为主。

（3）抗噪：GH法在T取值为[0，50]时无法检索到椒盐噪音与高斯噪音派生图，可见对噪音敏感；本文算法对椒盐噪音不敏感，对高斯噪音敏感。

原因分析为：本文算法继承了SUSAN边缘检测对椒盐噪音抗噪的特性。

（4）自由变换与平移：GH、GLCM法对自由变换与平移敏感，本文算法在该处理上体现其优势。

原因分析为：自由变换与平移没有改变图像的边缘点与角点间的相对位置，而本文算法的初始边缘矩阵灰度化转共生矩阵保留部分相对位置特性。

（5）旋转：GLCM法与本法均对旋转45°变换表现一定的敏感性。

原因分析为：
灰度共生矩阵本身对方向具有一定的敏感性。

（6）对比度与亮度变化：GH、GLCM法均对对比度与亮度变化不敏感，本文算
法对亮度变化体现一定的敏感性。

原因分析为：以上变化不改变灰度与纹理特征，但是对本文算法的预设阈值可能具有一定影响。

5.3 节实验对以上特点进行了进一步验证与分析。

5.3 同一主题不同变换检索实验
该组实验测试本文算法在样本为同一主题、不同变换处理情形下的检索性能。

该类实验样本选自Corel图像库的源图及按5.2节实验预处理后的派生图。

每次实验均选取某一主题的1张源图及其9张派生图，共计10张图片作为测试样本。

本类实验中相关数指检索结果中以示例图像为源图的派生图（包含源图）的图像数。

同法，可得其检索效果如图4所示、观测数据如表3所示。

该类实验进一步验证5.2节实验的特点总结。

5.4 实验总结
实验总结如下：
（1）查准率P与查全率R：5.1节与5.2节实验验证本文算法的R值随T值增长
最快，且当T值达到某临界点后，本文算法的P、R值相对GH、GLCM法更高。

（2）时空开销：5.1节与5.2节实验验证本文算法的时空开销最大，其原因在于
本文算法需要更多的矩阵相关运算。

（3）不变性特征：5.2节与5.3节实验验证本文算法除高斯噪音外，对其他各类
变换如旋转、倾斜、平移、缩放等仿射变换均表现一定的不变特性，其原因在于本文算法继承SUSAN边缘检测的算法特性。

文献[11]认为目前在基于内容的图像和视频检索中原来的典型研究内容基本成熟，期待有新的突破。

SUSAN算子应用于图像配准领域、灰度共生矩阵应用于图像检索领域均已成熟，本文创新在于将SUSAN初始边缘响应矩阵的统计值视为灰度值，
从而将其中元素与元素间的位置关系通过灰度化映射转换到灰度共生矩阵中去研究。

这样既能利用SUSAN算子快速高效且噪音鲁棒、局部不变等优势，又可利用灰度共生矩阵对纹理特征[12]统计的特性。

将SUSAN初始边缘响应矩阵灰度化是算法的关键点与论文的创新点，希望能够抛砖引玉，为相关研究提供创新思路，例如可探索SUSAN矩阵灰度化后与其他检索方法相结合产生优化算法。

然而，实验表明本文算法时空开销较大，要满足实时性应用需要还需提高算法性能。

同时，检索效果受限于初始边缘响应矩阵灰度化后所表现出的纹理特性。

另外，对比经典的GH方法和GLCM方法，本文算法只有当检索返回的图像数量T到达较
大的临界点后才能体现其性能优势，而用户通常关心T较小时的查准率，因此本
文算法的应用范围进一步受限。

这些不足均需在后继相关研究中改善提高或寻找新思路突破。

【相关文献】
[1]Lowe D G.Distinctive image features from scale-invariant
keypoints[J].InternationalJournalofComputerVision，2004，60（2）：91-110.
[2]Moravec H.Obstacle avoidance and navigation in the real world by a seeing robot
rover[D].[S.l.]：Stanford University，1980.
[3]Harris C，Stephens M.A combined corner and edge detector[C]//Proceedings of the
4th Alvey Vision Conference，Manchester，1988：147-151.
[4]Smith S M，Brady J M.SUSAN—a new approach to low level image
processing[J].International Journal of Computer Vision，1997，23（1）：45-78.
[5]Bay H，Ess A，Tuytelaars T，et al.Speeded-Up Robust Features（SURF）[J].Computer Vision and Image Understanding，2008，110（3）：346-359.
[6]Rosten E，Drummond T.Fusing points and lines for high performance
tracking[C]//IEEE InternationalConference on Computer Vision，2005：1508-1515.
[7]Rosten E，Drummond T.Machine learning for high-speed corner
detection[C]//European Conferenceon Computer Vision，2006：430-443.
[8]Gonzalez R C，Woods R E.数字图像处理[M].阮秋琦，阮宇智，译.3版.北京：电子工业出版社，2011：536-538.
[9]Haralick R M，Shanmugam K，Dinstein I.Textural features for image
classification[J].IEEE Trans on Systems，Man，and Cybernetics，1973：610-621.
[10]Baraldi A，Parmiggiani F.An investigation of the textual characteristics associated with gray level cooccurrence matrix statistical parameters[J].IEEE Transactions on Geoscience and Remote Sensing，1995，33（2）：293-304.
[11]章毓晋.中国图像工程：2012[J].中国图象图形学报，2013，18（5）：483-492.
[12]刘丽，匡纲要.图像纹理特征提取方法综述[J].中国图象图形学报，2009，14（4）：622-635.。