王 南,许道云
贵州大学 计算机科学与技术学院,贵阳 550025
左心室MR 图像分割是运用图像分割的方法对左心室内外膜进行准确分割,其结果可以为医生定量地分析左心室收缩功能,提供参数指标。目前,关于左心室的分割方法主要分为图像驱动模型、基于图谱的方法[1]、基于深度学习的方法[2-5]以及主动轮廓模型等。图像驱动模型包括基于阈值、基于区域信息以及像素相似性质分类等方法,该类方法通常具有分割速度快、实现较为简单的特点,但其不能与左心室的先验信息有效结合,对左心室边界弱边界、对比度不高的图像分割能力有限。基于图谱的方法在分割给定的图像时,需要与大型训练集匹配分割图像,该方法的分割结果受配准图像数据集的大小所影响[6]。基于深度学习的方法,在处理左心室图像时,具有鲁棒提取图像特征以及实现左心室全自动分割的优点,但完成该方法的神经网络训练,需要大量图像数据,并且训练时间较长。主动轮廓模型相较于其他方法,其优点是在高噪声的情况下,也能得到较好的分割效果,并且容易与其他几种方法进行结合,分割速度较快。
主动轮廓模型分为参数轮廓模型和几何轮廓模型。参数轮廓模型是Kass 等提出的传统Snake 模型及其改进模型[7-9],其主要是将轮廓曲线向检测物边缘演化的过程与求解能量泛函最小值的过程相关联。该方法在给定适当的初始轮廓后,可以自动收敛到能量最小状态,具有对噪声和对比度不敏感的优点,但其对初始轮廓的位置具有很强的依赖性,在图像的弱边界处容易陷入局部最优。
几何活动轮廓模型是基于曲线演化理论与水平集的方法,将平面曲线隐含的表达为三维曲面函数,通过曲面的演化来隐含地求解曲线演化[10]。C-V模型[11]是几何轮廓模型中一个具有代表性的基于区域信息进行图像分割的算法,其具有易于处理拓扑结构变化、可分割较为复杂的对象、容易实现低维到高维的扩展、对初始轮廓不敏感等优势,但同时具有一定局限性,该模型对强度不均匀图像会产生不良分割结果,需要反复初始水平集,造成计算成本高等限制,众多学者在该算法的基础上进行改进与创新[12-17]。在文献[13]中,卢小鹏等人将目标区域与背景区域的灰度平均值之间的差值添加到C-V模型面积项中,使得C-V模型对亮度不均匀图像的分割能力有所提升。在文献[14]中,Wang等人提出一种将全局和局部统计信息结合的水平集新方法,克服了C-V模型不能成功分割强度不均匀图像的缺点。Murtaza等人在文献[14]基础上,使用变异系数代替方差,提出新的图像变分模型PSM[15],在对图像分割实验中取得更好的分割结果。Li 等人在文献[16]中提出一种对原图像进行拟合的方法。文献[17]中,卢振泰等人提出了一种基于熵和局部邻域信息的高斯约束C-V 模型对左心室进行分割,该模型具有算法执行速度快、避免重新水平集初始化、提高分割精度。文献[18]中,Li 等人提出距离正则化水平集演化(DRLSE)方法在能量函数中增加了一个正则项,以避免在水平集演化过程中出现形状不规则和不稳定性等情况,并且在该模型不需要重新初始化水平集,但该方法对初始轮廓敏感且在对比度不高的图片中演化过程容易陷入局部最小值。在运用DLRSE算法对左心室MR图像内膜进行分割时,会因为受到乳突肌与小梁组织的干扰而造成曲线陷入局部最优的结果,为了避免这种情况的出现,文献[19]中,根据左心室内膜的先验知识,刘宇提出一种保持凸性的水平集演化模型,该方法能保证曲线演化过程中具有良好的凸性,最终演化到左心室内膜边缘,但该算法存在对初始轮廓敏感,对图像弱边缘分割不精确的缺点。
综上所述,在图像分割中DRLSE 算法具有对初始轮廓敏感,演化过程仅依赖图像梯度信息的缺陷。为降低DRLSE 算法对初始轮廓的敏感程度,提升分割弱边缘、强度不均匀图像的能力,提高算法对左心室MR 图像内外的分割精度,提出一种适用于弱边缘信息左心室MR图像的DRLSE分割算法,实现对左心室内膜与外膜分割。在本文的内膜算法中,首先,使用Li提出的图像拟合方法,拟合得到新图像,计算出原图像与新图像的图像差,替换原PSM模型中的图像差,形成PSM新的局部项;其次,将文献[13]中C-V 面积项与上一步得到的PSM 新的局部项添加进 DRLSE 中;再次,将 CPLSE 模型的保持曲线凸性的曲率因子添加到DRLSE 模型中;最终实现对左心室MR 图像的内膜分割。在对左心室外膜进行分割时,首先,运用内膜分割结果,结合左心室心肌各方向厚度基本一致的形状先验,构造外膜形状约束项,形成左心室外膜形状约束项;其次,将外膜形状约束项与内膜分割算法中的PSM新的局部项、文献[13]中的C-V 模型面积项嵌入到DRLSE 算法中;最终实现对左心室MR图像外膜的分割。最后通过实验,对本文提出的模型进行测试,对比实验结果进行分析。
在C-V 模型中,将Ω⊂R2作为图像范围,u:Ω→R作为给出的灰度图像,定义Ωin是在图像u上被闭合曲线C所划分的内部区域,Ωout是在图像u上被闭合曲线C所划分的外部区域。C-V模型的能量泛函定义如下:
其中,μ≥0 ,v≥0 ,λ1>0 ,λ2>0 是各个项的系数,Length(C)是曲线长度项有平滑曲线的作用,Area(inside(C))是Ωin面积项,一般设其系数v=0,后两项为保真项是用于驱动曲线C向真实边界演化。
水平集定义表示如下:
E(C)表示为水平集形式:
其中,c1、c2分别为Ωin、Ωout的灰度平均值,其计算公式如下:
δ是Dirac函数,H是Heaviside函数。
当曲线演化到目标边缘的时,泛函式(3)将达到最小,其最小化可以通过以下的梯度下降流来求解:
C-V模型可以很好地对图像进行分割,但对灰度不均匀、均匀强度重叠的图像分割效果很差。2010 年Wang等人提出LCV模型,该模型对图像灰度不均匀图像有很好的分割效果,但在图像均匀强度一致部分所产生的分割结果较差。
PSM 模型是 2011 年由 Murtaza 等人在 LCV 模型基础上提出,该模型具体如下:
其中,u*=gk∗u-u,gk是高斯核,u*:Ω*→R作为给出的新灰度图像,是曲线C将u*内外两部分区域,m1、m2分别表示区域灰度的平均值。水平集下的PSM公式如下:
其中:
PSM模型的最小化求解过程可以通过以下方程:
PSM 模型在灰度非均匀图像的分割结果要优于LCV、C-V模型。
DRLSE 模型作为一种基于区域的水平集算法,是由Li等人于2010年提出的模型。
定义水平集函数:
DRLSE模型的能量泛函定义为:
其中,μ≥0 ,λ >0 ,α∈R,Rg(φ)是距离正则项,Lg(φ)是加权长度项可以驱使边缘向零水平集向着目标边缘方向演化,Ag(φ)是加权面积项可以改变水平集的演化速度,它们分别被定义为:
其中,P为势函数,是边缘指示函数,Gσ是高斯核;∇ 是梯度算子。在DRLSE方法中,选取:
能量泛函式(14)的最小化可以通过以下梯度下降流来求解:
DRLSE模型允许有更灵活、更有效的初始化,在曲线演化过程中显著地减少了迭代次数与计算时间,但该方法对初始轮廓敏感,仅依靠梯度信息驱动曲线演化,模型容易陷入局部最优。
对左心室图像分割时,DRLSE 方法存在对初始位置敏感且无法有效分割弱边界图像的问题,为改善以上情况,本文运用Li提出的图像拟合方法得到新图像,将新图像与原图像之间的图像差,替换新的PSM 模型中局部项的图像差,将新的局部统计信息与梯度信息结合改进DRLSE模型。为降低新局部项受到图像灰度不均匀情况的干扰,将对分割灰度不均匀图像有一定贡献的文献[13]中C-V的面积项加入到DRLSE模型中,得到初步改进DRLSE,其水平集下公式为:
其中,Hε1(x)是选取与DRLSE中相同的Hε(x) ,Hε2(x)选取与PSM中相同的Hε(x),u*=f-u,m1、m2是曲线划分u*产生区域内外的灰度平均值,拟合图像f的公式为:
其中:
gε是卷积核函数。
运用|∇φ|来消除传统C-V模型中的Dirac函数对检测远边缘的抑制,并将边停函数g加入到演化函数中。基于此,改进的DRLSE 能量泛函方程的最小化求解对应的梯度下降流方程:
在内膜分割中,内膜形状类似于椭圆,而且DRLSE算法受到乳突肌或小梁的影响,内膜分割结果可能是凹形曲线,出现欠分割情况,将保持曲线凸性的曲率因子加入到上述改进的DRLSE 方法中,这里对文献[19]中kmean公式作修改,定义:
最终内膜分割算法水平集最小化演化方程如下:
在分割外膜时,外膜与肌肉组织的灰度接近,在使用DRLSE 模型对图像外膜进行分割时因受到弱边缘、对比度低等因素的影响,造成过分割或欠分割等不好结果,观察图像中左心室的特性,发现左心室心肌各方向厚度基本一致[10],本文设计的外膜形状约束项是采用已分割好的内膜与已知左心室最小厚度Rmin结合先验信息[20],结合先验信息形成形状约束,以实现对演化过程中曲线的约束,具体公式如下:
其中dis是外膜分割曲线c上的点与内膜曲线的最短距离。
最终外膜分割算法水平集最小化演化方程如下:
其中,ω >0 ,γ >0 是两项牵引力的系数。
本文方法的具体步骤如下:
(1)内膜分割
①设定模型参数,对初始水平集φ0,即φt(t=0),最大迭代次数T,0 ≤t≤T;
②根据当前φt,结合公式(4),计算更新当前区域内外平均值c1、c2;
③ 根据公式(20),拟合出图像ft,计算更新u*、m1、m2;
④ 据公式(22)、(23)、(24),计算更新k、s(k)in、kmean;
⑤根据公式(25),更新φt+1;
⑥ 判断t >T,若是,则停止;否则判断φt+1与φt轮廓是否相同,若是,则停止;否则,φt=φt+1,转回②。
内膜分割流程图如图1。
图1 内膜分割算法流程图
(2)外膜分割
①设定模型参数,对初始水平集φ0,即φt(t=0),最大迭代次数T,0 ≤t≤T;
②根据当前φt,结合公式(4),计算更新当前区域内外平均值c1,c2;
③根据公式(20),拟合出图像ft,计算更新u*、m1、m2;
④ 根据公式(26)、(27)、(28),计算更新dmean、s(k)out;
⑤根据公式(29),更新φt+1;
⑥ 判断t >T,若是,则停止;否则判断φt+1与φt轮廓是否相同,若是,则停止;否则,φt=φt+1,转回②。
外膜分割流程图如图2。
图2 外膜分割算法流程图
本文对提出的方法通过实验进行验证,对其结果进行定量分析,并与DRLSE、CPLSE方法进行对比。近年来,基于深度学习的方法在图像分割领域应用成为主流方法,对左心室也有较好的分割效果,选择U-Net 网络加入到对比方法中。实验进行验证所使用的数据集来源于http://www.cse.yorku.ca/~mridataset/的Cardiac MRI dataset,该数据集是从33 个受试者获得共7 980 张的心脏MR 图像,并对图像进行手工分割获得5 011 个分段MR 图像和10 022 个轮廓。每个实验受测对象由20 帧和8~15个切片组。对U-Net网络进行训练时,从数据集中选取4 200 张带有标签的左心室图像,将数据图像分为两部分,总数量10%作为测试数据集,90%作为训练数据集输入到U-Net网络之中进行训练,采用二分类交叉熵作为损失函数,经过训练100 次循环,得到训练好的U-Net网络,其内膜与外膜训练的损失函数值分别为0.001 与0.002。在数据集中选取不同受试者、不同帧、不同切片的左心室MR图像,分别使用DRLSE、DRLSE-shape、CPLSE、训练好的U-Net 网络和本文模型对MR图像进行分割。选择的实验环境为Windows 7 操作系统,CPU为AMD Athlon II X4 640,5.75 GB内存,实验平台是Python3.5。
从数据集中选取不同受试者选取不同帧、不同层、像素为256×256的左心室MR图像进行实验。分别赋予每幅图各自不同的内膜初始化位置,运用本文左心室内膜分割算法对MR图像进行分割,该算法涉及的参数设置为μ=0.04,λ1=2,α=-1.5,λ2=0.1,v=0.001,η=1。为每幅图赋予各自不同的外膜初始化轮廓,运用本文左心室外膜分割算法对MR 图像进行分割,其涉及的参数设置为μ=0.04,λ1=2,α=5,λ2=0.1,v=0.001,η=0.25,γ=5。
(1)给定一张左心室的MR图像,分别运用DRLSE、CPLSE及本文算法,在给与三种方法相同初始化位置的情况下,对MR 图像进行内膜分割,其分割结果如图3,其中每张图的红色虚线为手工轮廓,绿色曲线为算法分割结果。
图3 同一MR图像不同初始位置的分割结果
(2)运用(1)中的三种方法对不同类型左心室MR图片的内膜进行分割,结果如图4。
图4 正常的MR图像分割结果
图5~图7 中,情形 1、2、3 依次代表图像含有乳突肌、弱边界、乳突肌和弱边界的情况。红色虚线为手工轮廓结果,绿色曲线为算法分割结果。
图5 含情形1的MR图像分割结果
图6 含情形2的MR图像分割结果
图7 含情形3的MR图像分割结果
(3)运用DRLSE、DRLSE-shape、U-Net 网络以及本文算法对六个不同左心室图像进行外膜分割,最终MR图像内外膜分割结果如图8。
图8 不同左心室内外膜的分割结果
图8(a)~(f)中的红色曲线代表的是专家手工轮廓;白色虚线代表的是DRLSE算法对MR图像内外膜分割结果;紫色虚线代表的是DRSLE-shape对MR图像外膜分割结果;蓝色曲线代表的是U-Net 网络对MR 图像内外膜的分割结果;绿色曲线代表的是本文算法对左心室内与外膜的分割结果。
(4)将专家手工轮廓、本文算法分割好的内外膜结果进行二值化后得到结果如图9。
为了对本文提出的模型进行定量评价,采用相似系数(Dice Similarity Coefficient,DSC)、平均垂直距离(Average Perpendicular Distance,APD)。DSC是对两个区域重叠度的测量,越接近1,两幅图的重叠度越高;APD 是自动分割所得轮廓曲线上的全部点与手动轮廓曲线的最短垂直距离的平均值,其值越小代表轮廓曲线重合越高;质心距离(Centroid Distance,CD)是自动分割环形质心与手工分割环形质心的欧氏距离。两者指标定义如下:
图9 不同左心室心肌分割结果的二值图
其中,A代表专家已分割轮廓,Aor代表专家手工分割内外膜轮廓之间区域,B代表经过算法自动分割轮廓,Bre代表自动分割内外膜轮廓之间区域,|Ω|代表计算Ω区域内像素点的总数,例如,|A|+|B|代表两个区域内全部像素点个数,|A⋂B|代表两个区域相交部分像素点个数,xi代表B中的像素点,(xor,yor)代表Aor的质心,(xre,yre)代表Bre的质心。
表1 给出的是运用DRLSE、DRLSE-shape、CPLSE、U-Net 网络与本文算法对左心室内外膜分割所用平均时间。
表1 各分割算法的平均时耗s
表2给出的是运用DRLSE、CPLSE和U-Net网络与本文算法分割不同心脏左心室、不同层MR 图像内膜,所得分割结果的DSC与APD值。
表3给出的是运用DRLSE、DRLSE-shape、U-Net网络与本文算法分割不同心脏左心室、不同层MR图像外膜,所得分割结果的DSC和APD值。
表4给出的是,针对不同图像,分别采用本文所提算法对内膜与外膜的分割结果、U-Net网络对内膜与外膜的分割结果、DRLSE对内膜与外膜的分割结果、DRSLE对内膜分割的结果与DRLSE-shape 对外膜分割的结果进行组合得到左心室心肌结果,并将其结果与专家手工分割的心肌结果进行比较,所得到的DSC和CD值。
表2 不同图片内膜分割结果的APD和DSC值
表3 不同图片外膜分割结果的APD和DSC值
表4 不同左心室心肌分割结果的DSC和CD值
由图3知DRLSE、CPLSE依据梯度信息进行MR图像分割时,存在对初始轮廓敏感的问题,本文算法同时依据梯度信息和图像局部区域统计信息来驱动曲线演化,降低了算法对初始位置的依赖性。
由表1可以直观看出各算法的时间复杂度,运用各种方法对左心室内外膜分割,U-Net的平均消耗时间较短,但该方法对网络进行训练时,不仅需要大量图片标签数据,还会消耗较长时间。结合表2、3 可知,相较CPLSE、DRLSE-shape 及本文方法,DRLSE 虽然平均时间最短但其分割效果最不理想。本文算法、CPLSE 及DRLSE-shape方法的运行时间相近,本文算法在运行速度方面表现并不突出,主要因为本文内膜与外膜算法在DRLSE 基础上,两者除了都增加新局部信息及轮廓内外区域灰度平均值的计算之外,内膜算法添加了曲率因子的计算,外膜算法添加了形状约束项的计算,从而影响了内外膜算法的运行速度。
由图4~图7 与表2 知,对比 DRLSE、CPLSE 模型、U-Net网络,本文算法加入图像局部信息,对边缘处有加强的作用,在处理弱边缘、对比度低且有乳突肌影响的左心室图像时,可以得到更高的DSC平均值,本文算法所得轮廓包围的内部区域与专家手工分割所得轮廓包围的内部区域更加重合。分割结果APD平均值相对于对比中的其他模型分割所得APD 平均值更低,本文算法所得的轮廓更加靠近专家手工分割轮廓,有更高的分割精度。
由图8 与表3 知,深度学习方法对左心室MR 图像分割时,可能会无法对图像做出分割。在对DRLSE 算法加入本文先验形状约束项后,提高了DRLSE 方法的分割能力。本文外膜分割算法基于本文内膜分割效果的基础上进行的,内膜分割结果优于DRLSE、CPLSE模型、U-Net网络分割结果,并且在外膜分割方法中依靠了图像差的局部信息,所以本文外膜分割效果是对比方法中最好。
由图9 与表4 知,在对左心室整体的分割结果对比中,U-Net 网络对左心室心肌整体分割的质心CD 平均值最低,更接近手工轮廓质心,但其DSC均值在对比方法中低于本文算法的DSC均值,且本文的质心CD值较传统方法分割结果相较于传统分割方法更接近手工轮廓的质心,本文算法有更好的分割效果。
针对运用DRLSE 模型分割存在弱边界、对比度低的左心室图像时不能得到很好分割效果,为了提高左心室内外膜的分割精度,本文提出改进的DRSLE 算法内膜分割算法和带有形状约束的外膜分割算法。通过区域拟合图像与原图像比较,提供局部区域信息,将改进的C-V 面积项与局部区域信息与保持凸性的方法添加到DRLSE算法中,以此驱动曲线进行演化。实验表明,本文算法在对弱边界、对比度低的左心室图像进行内膜分割时,相对DRLSE、CPLSE算法、U-Net网络有更高的分割精度;在对弱边界外膜分割时,也取得相对稳定、较高的分割精度。但本文算法,需要手动给出初始轮廓,以及调节部分参数,在分割过程中仍依赖人工经验。未来的工作将是,需求新的方法、减少变参数量、降低模型对人工依赖,以达到低成本人工成本同时达到更高的分割精度。