基于曲线二叉树的微弱边缘检测方法

王 会，余 阳

(1.成都东软学院 计算机科学与技术学院，四川 成都 611844； 2.成都东软学院 电子商务与信息管理系，四川 成都 611844)

0 引 言

传统的边界检测[1.2]方法多基于微分技术，如基于一阶导数的Roberts算子、sobel算子、Prewitt算子、基于二阶导数的Laplacian算子[3]。这些方法对噪声敏感，检测精度不高。然而，实际图像的信噪比(signal noise ratio，SNR)并不高，如生物医学图像、遥感卫星图像、水下生物图像等[4]。因此，对低SNR图像进行准确地边缘检测一直是研究热点。

为了提高传统梯度算子的边缘检测性能，Duan等[5]提出了改进的canny算子，先对图像进行高斯平滑，再由一阶导数的极大值确定边缘点，但canny算子会引入伪边缘点。基于高斯滤波的二维拉普拉斯(Laplacian of Gaussian，LOG)[6,7]零交叉边缘检测方法应用广泛，该局部探测方法效率很高。然而，在低SNR情况下检测效果较差。为了消除噪声对边缘检测的影响，首先需要图像去噪，然后边缘检测。去噪方法包括双向滤波[8]、各向异性扩散[9]、非局部均值[10]等。然而，这些方法并非为边界探测专门定制，从而使边缘模糊，微弱的边界检测更加困难。

近些年根据自然图像的特点，提出基于多特征和模糊推理的边缘检测方法[11,12]。这类方法虽在一定程度上提高了自然图像的边缘检测准确性，但比较复杂，且对微弱边界检测性能有限。根据噪声与边缘均为高频信息的特性，基于多尺度方法得到应用，最常见的有小波变换[13]。由于小波变换的非奇异性，其不能够适应真实弯曲边界的形状。也有使用匹配滤波器来提高微弱边界的检测性能的，如Tsantis等[14]利用改进的模糊c均值实现超声图像的斑点噪声图像的检测；Wang等[15]利用光束曲线金字塔(beam curve pyramid，BCP)数据结构和次线性算法检测低信噪比边缘。虽然这些方法实现多尺度匹配，但时间复杂度较高，应用困难。

针对传统方法对微弱边界检测性能较低、时间复杂度较高的问题，本文基于人类视觉系统的特性，利用匹配滤波器探测微弱弯曲边界，但仅保留那些在统计学上很重要的弯曲边界。其主要工作包括：①设计多尺度匹配滤波器检测长边缘的微弱边界，提高检测准确性；②采用多线程和自下而上的搜索策略，提高边缘检测速度。

1 本文算法

本文的核心是通过曲线二叉树的方法构建多尺度匹配滤波器。

1.1 基本结构

为了方便阐述本文算法的理论，假设无噪图像I高度与宽度相等，边界为一条非交叉曲线λ，图像像素灰度值表现为阶跃不连续。对于总长度L，通过一系列像素集合P的每条候选曲线λ，其响应向量为

Φ(λ)=[R,L,C,P]

(1)

式中：R为响应值，由曲线λ相对应的匹配滤波器所决定，表示该曲线两边像素和之间的差异性；C=R/m(L)，表示沿着曲线方向的平均对比度，其中m(L)为长度L的匹配滤波器包含的像素总数。令λ1和λ2为两条具有公共端点的曲线，其相对应的响应向量为Φ(λi)=[Ri,Li,Ci,Pi],i=1,2，则其连接的响应向量为

(2)

图1 曲线二叉树

在每个分块V中，对于在其边界∂V上的每一对像素p1和p2，曲线二叉树数据结构CBT存储其相应的响应向量Φ(λ)，每一个响应向量确定p1和p2之间的一个单曲线λ，其很可能为一条边界。之后根据检测阈值，分配边界得分，确定边界曲线。

本文的分块基本原则是：子分块近似等于原分块区域的一半，分割边界的长度近似等于其面积的平方根。

1.2 边界检测算法

算法1给出CBT的构建伪码，其中注意分块V的最大边长为n。算法1中，各种曲线的匹配滤波响应是自下而上的方式计算得到，即从树叶到树根。如果n

然后，根据j+1层的响应向量，计算下一个粗糙层(第j层)的响应。本文通过连接兄弟分块实现该过程，具体如下所述。令分块V1和V2为第j+1层的兄弟分块，V为第j层的母分块，所有像素对遵循p1∈∂V1∩∂V，p2∈∂V2∩∂V。则对于每对像素(p1,p2)，所有像素均在连接线V12=∂V1∩∂V2中。对于每个像素p3∈V12，通过连接p1到p3、p3到p2的两条曲线，计算其曲线响应向量。在所有这些连接曲线中，存储得分最高的那条曲线。该粗糙处理过程的伪码如算法3所示。

算法1：构建曲线二叉树CurveBinaryTree(V)

Ifn

CBT←BottomProcess(V)

else

V1,V2←SubTiles(V)

CBT1←CurveBinaryTree(V1)

CBT2←CurveBinaryTree(V2)

CBT←CoarserProcess(V,V1,V2,CBT1,CBT2)

end

returnCBT

算法2：底层处理BottomProcess(V)

CBT←EmptySet

for ∀p1,p2∈∂V

λ← straight line fromP1toP2

CBT.add(Φ(λ))

end

returnCBT

算法3：粗糙处理CoarserProcess(V,V1,V2,CBT1,CBT2)

CBT←CBT1∪CBT2

if (Basic-Mode==1) then

InterfaceSet←∂V1∩∂V2

else if(Optimized-Mode==1) then

InterfaceSet←BestPixels(∂V1∩∂V2)

end

for ∀p1,p2:p1∈∂V∩∂V1,p2∈∂V∩∂V2

AllResponses←EmptySet

for ∀p3∈InterfaceSet

λ1← curve fromP1toP3in setCBT1

λ2←curve fromP3toP2in setCBT2

Φ(λ)←concatenate(Φ(λ1),Φ(λ2))

AllResponses.add(Φ(λ))

end

CBT.add(AllResponses.bestResponse())

end

returnCBT

本文算法最终输出为模糊边界图像Ei,j=0表示没有边界通过像素(i,j)。Ei,j值越大，表示边界通过的统计可能性越高。该边界图像E的生成具体如下所述。初始设定E的所有像素值均为0，对于每个响应变量Φ(λ)，分配一个重要性得分，取决于其均值对比度C和长度L，一个正得分表示该曲线可能为边界。然后，将该曲线二叉树中所有正得分的曲线从大到小排序。对于该排序列表中的每条曲线λ和每个像素p∈P，设定E(p)记录λ的得分。为了处理正得分的重合曲线，本文应用非极大值抑制方法：如果某些像素p的当前E(p)为正值，则该得分不再减小；如果，P中的大多数像素已经被之前得分较高的曲线标记，则丢弃当前曲线λ，不将其像素添加到E中。

2 算法分析与优化

2.1 算法复杂度分析

t(S)=2t(S/2)+O(S1.5)

(3)

根据主定理(master theorem)，总复杂度t(S)可转换为标准形式，如式(4)所示

t(n)=at(n/b)+f(n)

(4)

下面针对分块过程，计算其具体复杂度。第j=0层即根节点表示整个尺寸为n×n图像，第j=1层为两个尺寸为n×n/2矩形，第j=2层为4个尺寸为n/2×n/2的正方形。所以，当j为偶数时，该层包含尺寸为n/2j/2×n/2j/2的方块；当j为奇数时，该层包括尺寸为n/2(j-1)/2×n/2(j+1)/2的矩形。

在底层，对于每个叶分块，扫描所有从分块的一端开始到另一端结束的直线。这些响应值所需要的运算次数与图像的总像素个数N有关，所以该层运算次数的总数为O(N)。然后，在每个粗糙层j，通过连接j+1层两条短线来构建每条曲线，每条曲线在连接线上有一个公共点。所以当j为偶数时，该层的运行的总次数为

6N×n/2j/2=6N1.5×2-j/2

(5)

当j为奇数时，该层的运行的总次数为

6N×n/2(j+1)/2=6N1.5×2-(j+1)/2

(6)

将各层的运算次数进行求和，得到本文算法的总运算复杂度AC(N)，如式(7)所示

(7)

其中，jb为底层层号，即层数，jeven表示所有偶数层，jodd表示所有奇数层。对式(7)进行转换，得到式(8)，说明本文算法的复杂度小于18N1.5

(8)

2.2 算法复杂度优化

当处理大尺寸图像时，O(N1.5)的运算时间复杂度仍较大，效率较低，因此，需要开发效率更高的算法。

与之前一样，第j层分块V的两个子分块V1和V2，像素对p1、p2∈∂V。优化算法仍然以p1和p2之间最佳响应计算边界曲线。然而并非在连接线V12=∂V1∩∂V2中遍历所有的像素，而是仅考虑包括k个像素的子集。为得到该子集，对每个像素p3∈V12，本文寻找最高响应值的曲线，与之前一样，从∂V1或∂V2开始，至p3结束。然后，仅保留最高响应值的k个像素。

对于6N对起点和终点中的每一对，在连接线上只有最优的k个像素，所以，在第j层共有6Nk条曲线。另外，新增的运算是选择最优的k个像素。由于分块数目等于2j，连接线的长度约等于n/2j/2，所以对于第j层，该选择步骤的运算量由logk·2j·n/2j/2

3 实验结果与分析

本部分实验以visual studio c++(版本2012)为平台，应用多线程方式，处理器：英特尔Core i7 4800 M 2.9 GHz；CUP：16 G；操作系统64为Windows。利用模拟图像和真实图像检验本文算法的性能，并与传统的BCP算法和LOG边缘检测方法进行对比。

3.1 模拟图像

首先比较本文一般算法(算法复杂度为O(N1.5))与本文优化算法(算法复杂度为O(NlogN))的时间情况。对于尺寸128×128的图像，本文一般算法时间为0.25 s，本文优化算法时间为0.19 s；对于尺寸512×512的图像，本文一般算法时间为7.3 s，优化算法为3.1 s。两种算法的边缘检测时间与图像尺寸大小的变化关系，如图2所示。可以看出，本文优化算法所用时间远低于本文一般算法，同时，随着图像尺寸的增大，优化算法节省的时间越多。

图2 本文一般算法和优化算法时间对比

下面对图3(a)所示的模拟图像进行实验。该模拟图像包括狭长三角形、直线、同心圆和“S”形。为了验证本文算法的有效性，对该图像进行高斯加噪处理，使得SNR范围为0～2，间隔为0.2，同时添加影响程度2%的椒盐噪声。本部分利用F-Measure(FM)[16]评价标准量化算法的边界检测性能，定义如下式

(9)

式中：RC(recall)表示查全率，PR(precision)表示查准率。同时，本文一般算法、优化算法同传统的BCP算法和LOG算法相比较。不同算法的FM值随SNR的变化曲线如图4所示，其中SNR=0(纯噪声)、0.8(低SNR)、1.4(中SNR)、2.0(高SNR)的不同算法的FM值见表1。

图3 模拟图像

图4 不同算法FM值随SNR的变化曲线

算法00.81.42.0本文一般算法0.0010.430.850.93本文优化算法0.0010.400.830.92BCP算法0.0030.300.700.85LOG算法0.0050.170.520.68

由图4可以看出，本文两种算法的FM值比传统边缘检测算法有明显提高，本文优化算法的结果稍微低于本文一般算法，但非常接近。说明本文优化算法在大幅度提高运算时间的同时，可以有效保持检测精度。由表1可以看出，当图像中只含有噪声(SRN=0)时，所有算法的FM均很低，即不能探测到边界，并且本文两种算法具有更低的FM值，说明本文算法具有更好的抗噪性能。当图像的SNR为低级(SNR=0.8)时，本文两种算法的FM值为0.4以上，明显高于传统两种方法，说明本文方法检测微弱边界的性能有明显提高。当图像的SNR为中级(SNR=1.4)和高级(SNR=2.0)时，本文两种算法的FM值分别达到0.8和0.9以上，比传统方法有很大提高，说明本文方法具有良好的边缘检测性能。

图3(b)和图3(c)分别给出SNR=1.4和SNR=2.0时，模拟图像的加噪结果。图5和图6分别给出两种模拟加噪图像，本文优化算法和两种传统方法的边缘检测结果。

图5 SNR=1.4的边缘检测结果

图6 SNR=2.0的边缘检测结果

由图5可以看出，当SNR=1.4时，本文方法可以有效检测出图像中含有的边界信息，尤其是可以较好地检测出短直线；传统的BCP算法并不能检测出短直线和完整的圆；传统的LOG算法则不能有效检测到边界。由图6可以看出，当SNR=2.0时，本文算法可良好地检测到图像中的所有边界；传统的BCP算法仍不能检测到短直线，LOG算法的检测效果虽有所提高，但整体仍较差。说明本文算法对于微弱边界的检测效果明显优于传统方法，具有很好的抗噪性能。但是本文算法显得过于注重微弱边界，所以检测到的圆中含有较多短条纹。

3.2 真实图像

为了进一步说明本文算法良好的边缘检测性能，对真实自然图像进行实验。首先对自然景观图像进行边缘检测，然后对生物医学图像进行边缘检测。同时比较本文优化算法和传统BCP算法和LOG算法的检测结果。

3种方法对自然景观的检测结果如图7所示，分别对树干年轮、树叶叶脉、河提和大桥4种自然景观图像进行检测。可以看出，本文算法可以有效的检测出不同自然景观中的细小纹理，尤其是对年轮的检测，几乎可以检测出所有的年轮。而传统方法只能检测部分纹理边缘，对细节检测能力较弱。

图7 自然景观图像检测结果

生物医学图像多属于微观世界，所含的边缘纷繁复杂、细小微弱，对其进行有效边缘检测具有十分重要的意义。本实验分别对神经元、视网膜血管、微生物、细胞结构4种生物医学图像进行检测，检测结果如图8所示。可以看出，本文算法可以良好地检测到传统方法难以检测到的细节。本文算法对生物医学图像细节的良好检测性能，可以为目前生物医学领域的研究做出重要贡献。

图8 生物医学图像检测结果

4 结束语

本文主要研究在噪声图像中实现微弱边界检测的算法。设计多尺度匹配滤波器，检测任意形状和长度的微弱边界，使本文算法对不同应用的图像均可以实现良好的检测；利用多线程和自下而上的搜索策略提高边缘检测速度；同时只保留有效中间元素，对算法进行优化。通过对模拟加噪图像和自然景观、生物医学的真实图像进行实验，本文算法对微弱边界可实现有效的检测。

由于本文算法的焦点多集中于微弱边界的探测，所以对于噪声图像，有时会出现细小的伪边缘。如何更好地处理该种情况，进一步提高边缘检测性能，是下一步研究的重点。