作者简介:李 璟(1991—),男(汉族),江西省九江市人,深圳大学硕士研究生. E-mail:lijingszu@gmail.com
中文责编:英 子; 英文责编:雨 辰
1)医学超声关键技术国家地方联合工程实验室, 广东省生物医学信息检测与超声成像重点实验室, 深圳大学医学院,深圳 518060; 2)深圳市妇幼保健院超声科,深圳 518028
Li Jing1, Ni Dong1, Li Shengli2, Han Xiao1, Yin Xiaolang1, Wang Tianfu1, and Chen Siping11)National-Regional Key Technology Engineering Laboratory for Medical Ultrasound, Guangdong Key Laboratory for Biomedical Measurements and Ultrasound Imaging, College of Medicine, Shenzhen University, Shenzhen 518060, P.R.China2)Department of Ultrasound, Shenzhen Maternal and Child Healthcare Hospital, Shenzhen 518028, P.R.China
image processing; ultrasound image processing; phase symmetry(PS); edge detection; ellipse fitting; random forests algorithm; measurement of the fetal head circumference
DOI: 10.3724/SP.J.1249.2014.05455
提出一种超声图像中胎儿头围自动测量的新方法.利用机器学习的随机森林(random forests,RFs)算法自动检测感兴趣区域(region of interest,ROI),通过图像局部相位对称(phase symmetry,PS)检测头围边缘,使用非迭代椭圆拟合算法拟合出头围椭圆.与医生手动拟合测量的结果对比,145个头像的平均相对偏差为-3.86 mm,表明该方法可以鲁棒的自动检测胎儿头围.
A novel learning based method is proposed to automatically measure the fetal head circumference in ultrasound image. Firstly, the random forests classifier is utilized to detect the region of interest(ROI). Then, the fetal head contour in the ROI is segmented using the local phase symmetry-based method. Finally, the non-iterative ellipse fitting method is employed to determine an ellipse on the segmented head contour. Experiments on 145 images show that the mean relative deviation between automatic and manual measurements is -3.86 mm, which demonstrates the feasibility of the proposed method.
医学超声成像无辐射、实时且廉价,已成为胎儿产前诊断的首选影像模式.医生主要通过手动测量胎儿的头围、腹围和股骨长等生物参数,预测或确认孕妇的孕周,估计胎儿的尺寸大小和质量,判断其生长发育状况是否良好和有无畸形等.但常用的人工测量方法存在以下问题:① 测量结果的准确性依赖于医生经验; ② 诊断时间长且效率低; ③ 超声医生由于频繁重复劳动容易引起重复性压力损伤.因此,近年来针对胎儿生物参数的自动测量方法的研究被广泛关注.
胎儿头围自动测量的主要困难在于超声图像中通常存在信号衰减、边缘缺失、伪影和斑点噪声(speckle noise)等.在头骨分割和边缘检测方面,Hanna等[1]用形态学开操作去除头围以外的噪声,Lu等[2-3]用K-均值聚类的方法得到头骨高亮部分,这两种方法都基于图像灰度信息,对图像质量依赖较大,不同图像在相同参数情况下难以都达到最佳分割效果.Chalana等[4-6]提出基于活动轮廓模型(active contour models,ACMs)的方法检测头围边缘,但亦有不足:① 要提供种子点(seed point),导致检测不能全自动进行; ② 训练轮廓各点时对手工操作的正确性要求很高. Shan等[7-8]提出基于可形变模型的方法检测头围边缘,但该模型基于Rayleigh分布,没有考虑纹理结构的空间特征[9],边缘检测正确率不高.在椭圆拟合部分,大部分采用经典的Hough变换拟合方法[3,10].然而,Hough变换拟合的缺点在于:① 对于部分缺失的椭圆拟合准确率低,不够鲁棒; ② 耗时较长.
针对传统超声胎儿头围自动测量方法中边缘检测不准确和头围拟合不够鲁棒等问题,本研究采用3个步骤实现胎儿头围的自动测量.首先,利用随机森林(random forests,RFs)算法从超声图像中自动得到胎儿头围的感兴趣区域(region of interest,ROI),从而避免对全图进行处理而影响准确性; 其次,使用相位对称性(phase symmetry,PS)对ROI进行边缘检测,大大减少噪声干扰; 最后,使用非迭代椭圆拟合(ellipse fitting,EllFit)算法进行头围拟合,解决头围图像部分缺失问题.实验结果表明,该算法能够自动准确测量头围.本研究算法的系统结构图如图1.
以往基于监督学习的相关研究[11-12]中,大多采用Adaboost或支持向量机(support vector machine,SVM)分类器来检测解剖结构或者标准切面.近年来, RFs分类器作为一种新型分类器被广泛用于目标检测.RFs分类器相比传统分类器有如下优点:一是通常不会出现过拟合现象,且对异常值和噪声具有很好的鲁棒性; 二是RFs分类器在训练过程中能提供各种输入特征的相对重要性,有利于提高目标检测与分类的准确性.为此,本研究提取训练样本的Haar特征,采用RFs算法训练出用于检测头围ROI的分类器.1个RFs分类器由多棵二值决策树构成.为构建决策树t,其训练数据是原始训练数据的一个随机子集,通过迭代生成新节点,将训练数据集划分成多个有显著区别的子集.当决策树的深度达到设定的最大深度,或者当节点的训练数据集小于一定规模时,节点停止分裂.设Di是输入的某节点的训练数据,通过寻求弱分类器的最优参数来最大化节点分裂前后的信息增益(information gain,IG),并据此将该节点的Di分裂为Dj和Dk. 分裂过程中的IG定义为
IG(Di,Dj,Dk)=E(Di)-(|Dj|)/(|Di|)E(Dj)-
(|Dk|)/(|Di|)E(Dk)(1)
其中, E(Di)、 E(Dj)和E(Dk)分别是数据集i、 j和k的熵.
为训练分类器,从超声图像中分别截取头围ROI图像作为正样本,并从背景图像中随机截取子图作为负样本.部分正负样本间存在 20%~40%的重叠.据统计头围长宽比约为1.2:1.0,故将所有截取的样本图像均缩放至尺寸为 96×80 像素的图像.训练得到分类器后,在超声图像的全图范围内利用RFs分类器定位得到ROI,缩小后续图像处理的范围,从而提高精确度.
自Morrone 等[11]提出用局部相位(local phase)做特征检测,相位信息作为基础特征被应用到各种图像的处理中.本研究使用局部相位对称性来进行胎儿头围的边缘检测.下面对局部相位和相位对称性做基本描述,更多细节见文献[13-14].
正交滤波器通过信号的幅度和相位可计算得到信号的局部属性.本研究使用Log-Gabor滤波器,其频域表达式为
G(ω)=exp{-([lg(ω/ω0)]2)/(2[lg(k/ω0)]2)}(2)
其中,k为尺度参数; ω为滤波器频率; ω0为滤波器的中心频率; 比值k/ω0和滤波器带宽β相关,如β=-(221/2)/((ln2)1/2)×ln(k/ω0).
此时,通过不同尺度的Log-Gabor滤波器可得到空域和频域信号的局部信息.尺度由人为设定的最小波长λmin决定.这里,最小波长λmin、滤波器尺度m和中心频率ω0的关系为
ω0=1/(λmin δm-1)(3)
设Mem(x)=Re{F-1[G(ω)]}和 Mom(x)=Im{F-1[G(ω)]}分别代表Log-Gabor滤波器在尺度为m时的实部和虚部,其中F-1表示傅里叶反变换,则一维信号I(x)在尺度为m时的振幅Am(x)、 局部相位φm(x)和局部能量E(x)依次为
Am(x)=(em(x)2+om(x)2)1/2(4)
φm(x)=arctan((om(x))/(em(x)))(5)
E(x)=([∑mem(x)]2+[∑mom(x)]2)1/2(6)
其中, em(x)=I(x)*Mem(x); om(x)=I(x)* Mom(x),符号*代表卷积运算.经此计算,可获得在I(x)上的每一点x对应的局部相位信息.
利用式(4)至式(6)获得局部相位信息后,可采用Kovesi[15]提出的方法计算不同尺度下每个点的相位对称性值,即
PS(x)=(∑m[|em(x)|+|om(x)|]-T」)/(∑m(em(x)2+om(x)2)1/2+ε)(7)
其中, A」=max(A,0); ε为避免除数为0的小数值; T为噪声阈值[13].
在实际的超声图像中,特征检测包含多个不同方向r. 因此,本研究使用含方向的2维Log-Gabor滤波器进行处理,该滤波器包含了径向分量和角度分量,如
G(ω, φ)=
exp[-(((lg(ω/ω0))2)/(2(lg(k/ω0))2)+((φ-φ0)2)/(2σ2φ))](8)
其中, φ为角度坐标; φ0为滤波器的方向角度; σφ=Δφ/s, 它决定了角度带宽ΔΩ=2σφ(2lg2)1/2, Δφ是相邻角度的间隔, 且Δφ=180°/Nr, Nr决定了分为多少个方向.根据经验设Nr= 6, s=1.2.
通过不同尺度m和不同方向r的二维Log-Gabor滤波器,可得图像中每个点的二维PS值为
PS(x,y)=
(∑r∑m[|erm(x,y)|- |orm(x,y)|]-Tr」)/(∑r∑m(erm(x,y)2+orm(x,y)2)1/2+ε)(9)
图2是边缘检测各步骤结果的示意图.其中,图2(a)是由RFs分类器得到的头围ROI; 图2(b)是对ROI求PS值的结果,PS值相对较大处对应着头骨的边缘部分; 图2(c)是对ROI模板滤波的结果,滤波模板是根据胎儿头围ROI中噪声的分布而设计的,如式(10),其中, e1(x,y)=(x-x0)2/a2+(y-y0)2/b2, e2(x,y)=(x-x0)2/(0.6a)2+(y-y0)2/(0.6b)2, 且(x0, y0)是全图的中心点, a和b分别是图像的长和宽的一半; 图2(d)是对图2(c)进行非极大值抑制,双阈值门限和形态学开运算后得到的头围边缘.
H(x,y)={0, e1(x,y)≥0.9, e2(x,y)≤1.2
1, e1(x,y)<0.9, e2(x,y)>1.2(10)
由于胎儿头骨未完全闭合,以及受超声图像中信号衰减和噪声的影响,图像中头围轮廓常出现或大或小的缺失,故在测量胎儿头围时需要用有效的椭圆拟合算法来拟合出正确的头围结果.
Prasad等[16]提出了基于最小二乘法的非迭代的几何椭圆拟合(EllFit)算法,有效减小了噪声和边缘缺失对拟合结果的影响.
假设2维椭圆的半长轴为a, 半短轴为b, 角度方向为θc, 中心为C(xc,yc), 则几何方程为式(11). 其中, a,b,θc,xc和yc的约束条件C1至C4如式(12).这里,R指实数集; R+指大于0的实数集.
([(x-xc)cosθc-(y-yc)sinθc]2)/(a2)+
([(x-xc)sinθc+(y-yc)cosθc]2)/(b2)=1(11)
{C1: a,b∈R+;
C2: b≤a;
C3: θc∈[0,π);
C4: xc,yc∈R(12)
若设 α=a2sin2θc+b2cos2θc, β=a2cos2θc+b2sin2θc, γ=(a2-b2)sin(2θc),则式(11)可简化为
α(x-xc)2+β(y-yc)2+
γ(x-xc)(y-yc)=a2b2(13)
椭圆上的点Pi(xi,yi)的切线方程为y=mix+ci. 其中, mi为斜率; ci为截距, i=1,2,…,N, 这里N为椭圆上点的数目.假设有一个离椭圆上的点Pi(xi, yi)最近的点P'i(x'i, y'i), 其中x'i, y'i∈Z, 则这两点的距离为
di=|(y'i-mix'i-ci)/((1+m2i)1/2)|(14)
为得到椭圆的各个参数,需使di最小,即
min(di2)=min(((y'i-mix'i-ci)2)/(1+m2i))
约束于C1~C4(15)
式(15)得到极小值的情况是(d2i)/mi=0, 2(d2i)/2mi>0, 根据表达式(d2i)/mi和2(d2i)/2mi易推出y'i-mix'i-ci=0是让d2i最小的最佳方法.但是,由于 P'i(x'i,y'i)并不都在椭圆上,故构造式(16)让ri尽可能接近0.
ri=|y'i-mix'i-ci|(16)
假设P'i(x'i,y'i)与Pi(xi,yi)十分接近,满足式(17)的约束条件时,用P'i(x'i,y'i)代替Pi(xi,yi),通过切线方程得到m'i和c'i,由此证明当m'i/mi→1, c'i/ci→1时, r'i/ri→1, 如式(18).从而,通过最小化差值r'i来间接最小化式(15)中的di.
{|x'i-xi|≤0.5
|y'i-yi|≤0.5(17)
r'i=|y'i-m'ix'i-c'i|(18)
建立矩阵和映射关系如式(19)~式(22)
A=[a b θc xc yc]T, 约束于C1~C4(19)
Φ=[φ1 φ2 φ3 φ4 φ5]T, Φ∈R5(20)
Y=[y'1 y'2 … y'N]T, Y∈ZN(21)
K: A→Φ; X: Φ→Y(22)
其中, A包含了椭圆拟合所需参数; Y包含了各个点的y轴坐标; Φ是一个新定义的过度矩阵,将非线性映射A→Y分成K和X两个部分, 如式(22).根据文献[16],可证映射X是线性的,故需最小化的差值距离r'i为
∑Ni=1(r'i)2==Y-X·Φ=2(23)
这里, =·=为矩阵范数.
如此定义的目的是,利用广义逆矩阵来唯一确定Φ, 从而最小化r', 即
Φ=(XT·X)-1·XT·Y(24)
为获得与式(23)相同的形式,将φ1至φ5写为式(25)的样式, X写为式(26)的矩阵[16]形式.
{φ1=α/β
φ2=γ/β
φ3=2φ1xc+φ2yc
φ4=2yc+φ2xc
φ5=φ1xc2+yc2+φ2xcyc-a2b2/β(25)
X=[
-(x'i)2/y'i -x'i x'i/y'i 1 -(y'i)-1
](26)
为求得含拟合参数的向量A, 建立φ1至φ5对应a, b, θc, xc, yc的逆映射,如
{a=2((φ2φ3φ4-φ24φ1-φ23-φ5(φ22-4φ1))/((φ22-4φ1)[(1-φ1)-((1-φ1)2+φ22)1/2]))1/2
b=2((φ2φ3φ4-φ24φ1-φ23-φ5(φ22-4φ1))/((φ22-4φ1)[(1+φ1)+((1-φ1)2+φ22)1/2]))1/2
θc=-0.5tan-1(-(φ2)/(1-φ1))(27)
xc=(φ2φ4-2φ3)/(φ22-4φ1)
yc=(φ2φ3-2φ1φ4)/(φ22-4φ1)
为保证映射的稳定性,文献[16]将X、 Y和r'i的函数修改为式(28)至(30),从而保证当y'i=0时,不会出现无穷值.
X=[
-x'2i -x'iy'i x'i y'i 1
](28)
Y=[(y'1)2(y'2)2 …(y'N)2]T, Y∈ZN(29)
∑Ni=1[y'i(r'i)2]==Y-X·Φ=2(30)
此外,为避免式(30)出现等式左边因y'i=0而r'i≠0时整个函数为0的情况,可用令(~overx)i=x'i-(~overx)和(~overy)i=y'i-(~overy)的((~overx)i,(~overy)i)代替(x'i, y'i). 其中,
{(~overx)=1/N∑ix'i
(~overy)=1/N∑iy'i(31)
综上,EllFit算法的过程如表1.
利用非迭代的椭圆拟合算法拟合出的 HC椭圆是头骨的中间层,而实际临床中的HC结果是头围外边缘,计算公式为
CH≈π/2×(l+s)(32)
其中, l和s是HC外缘的长轴和短轴.假设图像中头骨的平均宽度为(-overt),则l=2a+(-overt), s=2b+(-overt). 这里,(-overt)值为[3]
(-overt)=1/M∑Mi=1ti=1/M∑Mi=1(Si)/(Li)(33)
其中,Si为各个头骨区域的面积; Li为对应头骨区域的长度; M为头骨区域的个数.各头骨区域Si是通过对头围ROI进行顶帽变换和形态学开运算后得到的面积较大的区域(>80 像素),Li是对Si进行骨架化后的结果.
为验证算法的可行性,首先对头围ROI的检测结果进行评价; 其次将实验结果与医生的手工结果进行定性和定量比较,观察有无显著性差异; 最后,分别对基于相位对称性边缘检测算法和基于形态学的方法,以及非迭代椭圆拟合算法和Hough变换椭圆拟合算法进行定性和定量比较.本研究使用1 048幅超声胎儿头围图像作为训练样本集来得到分类器,用145幅经医生测量过的20~36孕周的超声头围图像作为检测样本集.所有样本均取自深圳市妇幼保健院,通过Siemens Acuson Sequoia 512 超声仪器获得.
本研究利用RFs分类器自动检测头围ROI.以医生手工截取的头围ROI为标准(ground truth),检测结果如表2.表2中,误检指自动检测的ROI与医生手工截取的ROI面积重叠率小于75%,准确检测指自动检测的ROI与医生手工截取的ROI面积重叠率不小于75%.
首先,将实验结果和医生手工结果进行定性比较.图3显示本研究结果和医生手工结果的定性对比,白色虚线是医生手工拟合结果,绿色虚线是本研究自动测量结果.图3(a)为准确测量的典型结果图,图3(b)为测量错误的结果图,是由于靠近头骨右侧的其他解剖结构干扰导致.
其次,将实验结果和医生手工的结果进行了定量比较.定量评价标准如下:
定义临床手动测量值为vm, 自动测量值为va, 则两者之间的相对偏差Er=vm-va, 绝对偏差Ea=|vm-va|. 平均相对偏差(mean relative deviation,MRD)与平均绝对偏差(mean absolute deviation,MAD)分别为
MRD=1/n∑ni=1Eri(34)
MAD=1/n∑ni=1Eai(35)
其中, n为待检测图像幅数.本研究用t检验(t-test)判断自动测量与手动测量之间是否存在显著性差异(当显著性水平α=0.05时, t检验概率P>0.05称差异不显著).
根据以上评价标准,对自动测量的结果进行统计,得到自动测量结果与医生手工测量结果的平均相对偏差为-3.86 mm,标准差为7.18 mm.平均绝对偏差为5.57 mm,标准差为5.94 mm. t检验概率P为0.11,故无统计学差异.
图4绘制了自动测量相对偏差的结果箱线图.从图4可见,下4分位数为-6.90 mm,上4分位数为-0.45 mm,自动测量的整体结果和医生手工拟合的结果接近.
将基于相位对称性的边缘检测算法与基于形态学的方法进行对比.基于形态学的方法是指使用顶帽变换、开运算和骨架化等算法对图像进行边缘检测.2种方法的拟合部分都采用EllFit算法.
首先,通过图5定性比较基于PS的边缘检测算法和基于形态学方法的效果.其中,图5(a)是对一幅头围ROI图像进行形态学处理的结果; 图5(b)是用本研究提出的PS算法对图像进行处理的结果.可见,PS处理后,更好的保留了头围边缘,且很好的抑制了噪声等干扰.但采用形态学的滤波方法噪声干扰依旧很大.图5(c)将两种方法的最终结果绘制在医生手工拟合结果图上做对比.图中白色虚线是医生手工拟合的结果,绿色实线是本研究算法的结果,黄色实线是基于形态学方法的结果.从中可明显观察到本研究算法的结果和医生手工结果更加吻合.
其次,定量比较PS和基于形态学两种边缘检测方法.表3第1、2栏统计了分别使用PS和形态学边缘检测的结果与医生手工拟合结果的MRD、MAD和t检验概率.根据表3可观察到PS的边缘检测方法的MRD和MAD都明显小于形态学的方法所得统计结果.说明用PS对超声胎儿头围的边缘检测效果比用形态学的方法更好.
将EllFit算法与标准Hough变换(Hough transform,HT)拟合算法进行对比,边缘检测部分都使用PS算法.
首先,通过图6定性对比EllFit算法和HT算法.其中,图6(a)为头围ROI; 图6(b)是用PS对(a)图边缘检测后的结果; 图6(c)和(d)是分别使用HT算法和EllFit算法得到的拟合结果,可观察到HT在边缘缺失较严重的情况下不能正确拟合,而EllFit仍能拟合出正确的结果; 图6(e)显示了两种拟合算法的结果与医生手工结果的对比,图6(e)白色虚线是医生手工拟合结果,绿色实线为使用EllFit算法的结果,红色实线是使用HT算法的结果.由图6可明显看到,使用EllFit算法的结果更准确.
其次,定量比较EllFit和HT两种椭圆拟合方法,结果见表3.表3第1和第3栏统计了分别使用EllFit和HT椭圆拟合算法的结果与医生手工结果的MRD,MAD和t检验概率.由表3可见,采用EllFit拟合方法的MRD和MAD都明显小于HT拟合方法所得统计结果.采用HT算法所得的t检验概率P=0.01, 说明所得结果与医生手工结果有显著差异.时间复杂度分析结果显示,HT的复杂度为O(N4)[17], 而EllFit的复杂度为O(N). 实际应用中,采用Matlab软件,经过矩阵化优化, EllFit算法耗时约0.5 s/幅,而HT算法耗时约20 s/幅,故EllFit算法无论从椭圆拟合效果还是耗时方面都优于HT算法.
本研究发现,EllFit算法优于HT算法的原因在于HT是一种聚类算法,它将图像中的每个像素点映射到5维的参数空间中,再找出各个5维空间交汇面的峰值,最终得到椭圆方程各个参数.由于超声胎儿头围图像中受信号衰减、伪影和斑点噪声等影响,采用HT算法难以准确得到头骨中心椭圆线.此外,当对头围ROI进行边缘提取后,进行HT的数据点大大减少,边缘部分缺失且不规则弯曲较多,导致HT拟合效果差.EllFit算法基于最小二乘(least square)法的几何距离最小化.对头围ROI进行边缘提取后,在数据点少而精的情况下更能拟合出准确的头骨中心椭圆线.
综上可得,采用本研究提出的PS边缘检测和EllFit椭圆拟合算法进行超声胎儿头围自动测量是可行的.在台式计算机(CPU为Core i5,3.2 GHz; 内存为4 Gbit)上的运行评估表明,自动测量的速度可达约1 s/幅.因此,本研究实现的超声儿头围自动测量在产业化方面有较好的应用前景.
本研究提出一种超声图像中胎儿头围自动测量的新方法. 利用相位对称性进行边缘检测,有效减小了斑点噪声和伪影的干扰,能较好的提取头骨边缘,而且克服传统方法需初始化的半自动的缺点.采用非迭代的椭圆拟合算法,对于头围边缘部分缺失的情况,拟合准确率较传统拟合算法有显著提高.将所得到的头围结果与医生手工结果进行对比,结果接近且无显著性差异,验证了本研究提出算法的可行性,其速度满足实时性要求.今后,将对相位对称性算法进行参数自适应优化,并研究超声胎儿其他解剖部位的生物参数自动测量,最终将研究成果形成整体软件系统,应用于实际的超声产前检查中.
深圳大学学报理工版
JOURNAL OF SHENZHEN UNIVERSITY SCIENCE AND ENGINEERING
(1984年创刊 双月刊)
主 管 深圳大学
主 办 深圳大学
编辑出版 深圳大学学报理工版编辑部
主 编 阮双琛
国内发行 深圳市邮电局
国外发行 中国国际图书贸易集团有限公司(北京399信箱)
地 址 北京东黄城根北街16号
邮 编 100717
电 话 0755-26732266
0755-26538306
Email journal@szu.edu.cn
标准刊号 ISSN 1000-2618
CN 44-1401/N