用改进的耦合水平集方法从MSCT中分割左心室

王兴家 董利娜 李传富 范 亚 冯焕清*（中国科学技术大学电子科学与技术系，合肥 3007）

2（安徽中医学院第一附属医院影像中心，合肥 230031）

引言

随着多层螺旋CT（multi-slice spiral CT，MSCT）技术的快速发展，64层乃至更高层数（如128、256）的MSCT装备进了许多医院，320层的机器也已问世。它们通过多排探测器技术来提供高扫描速度、高时间分辨率的体数据，同时采用心电门控技术来同步扫描，做到每个心动周期扫描一周，能显著减少心跳的伪影；通过指定不同的同步位置，还能获得心脏不同搏动时相的体数据。例如，以0.1s间隔，取R峰及前后共10个位置为扫描起始点，各心动周期成像一次，就可获得表征一个心动周期不同搏动时相的10个CT数据集，实现所谓的“四维CT成像”。对该4D-CT数据集做进一步的处理和分析，就可能形成心、肺、冠状动脉等器官的动态图像，进行心、肺功能的更深入研究。

四维CT数据处理的第一步是针对各个时相的数据集进行心脏或肺的分割，工作量很大，必须实现自动分割，而且算法要能保证较高的分割精度，这是后续研究顺利进行的根本保证。由于左心室在心脏功能研究中的重要地位，因此本研究以它为例，介绍从4D-CT心脏数据中分割出左心室腔和心肌的方法。在MRI心脏成像领域已有不少左心室分割的研究工作，但针对MSCT的左心室分割的报道较少，已有的研究基于主动形状模型或主动表现模型，需要大量训练数据集，且不适合病变个体［1-2］。采用嵌入先验知识的改进耦合水平集（improved couple level set，ICLS）模型，从MSCT心脏数据中精确提取左心室腔和心肌。下面在介绍水平集模型的基础上，重点讨论如何针对左心室分割的特点和先验知识，分别对水平集耦合因子、水平集起始轮廓、边缘检测因子进行改进的思路，并给出采用特定心相MSCT数据集进行验证的结果。

1 原理和算法

1.1 水平集模型

水平集（level set）方法首先由Osher和Sethian提出，并用于曲线演化模型［3］。Caselles与Malladi等将该方法引入图像处理领域，提出了几何主动轮廓模型，使初始曲线向能量泛函的极小值演化［4-5］，利用水平集方法，使模型能自适应被检物体的拓扑变化，且与参数无关。近年来，出现了许多基于水平集的图像分割算法［6-9］。

在水平集模型中，用C表示演化曲线，C（t）=｛（x，y）|φ（t，x，y）=0｝是水平集函数φ（t，x，y）的零水平集，其中初始化轮廓为｛（x，y）|φ0（x，y）=0｝。假设I为待分割图像，为使外部能量驱动零水平集函数向目标边缘移动，构造图像梯度的递减函数为边界检测因子g，使得g在图像边缘的值近似为0。通常将g定义为

式中，Gσ是标准差为σ的高斯函数，*表示卷积运算。

边界检测因子g在图像分割中扮演极其重要的角色。从式（1）可看出，在图像梯度较大区域，g值近似为0，在灰度一致区域为正值。采用文献［7］的WRLS模型为原始水平集模型，该模型包含了额外的能量项，用于控制水平集函数和符号距离函数的位置偏差，避免计算量大且复杂的重新初始化。能量函数由惩罚项作为内部能量P（φ）和外部能量εg，λ，ν（φ）构成，有

式中，λ（λ＞0）和ν为常数，δ为狄拉克函数，H为海维塞德函数。

外部能量项εg，λ，ν（φ）驱动零水平集曲线向目标边界移动，同时内部能量项μP（φ）约束水平集函数趋近符号距离函数，从而保持水平集函数的平滑性。最小化能量函数，得到水平集演化方程，即

下面将针对MSCT图像的特点以及左心室的先验知识改进WRLS模型，使原始模型能适用于左心室的定位和精确分割。

1.2 水平集耦合因子

在结构上，左心室可看成由类圆形的腔和类环形的心肌构成，因此采用了不同的水平集函数φ1、φ2来表示左心室的腔和心肌。为避免两水平集函数发生重叠，定义一个耦合水平集（coupled level set，CLS）因子，它的正负和权重控制着水平集曲线的方向和演化速度，即

式中，d表示从左心室重心到内、外心膜的平均距离，w表示平均距离的容许范围。

在实际应用中，可根据左心室的结构和大小对d、w取值。通过耦合水平集，水平集函数在内、外心膜的对应区域演化而不相互影响。

1.3 起始轮廓的初始化

1.3.1 左心室定位

由于分辨率很高，而且Z轴分辨率与层片内的分辨率一致，即数据是各向同性的，所以MSCT数据具有结构连续性的显著特征，即相邻层片间的结构变化是连续的，结构差异不大，灰度具有较大的相似性［10］。根据结构连续特征，提出MSCT左心室的定位算法和粗轮廓提取算法，其中粗轮廓用于水平集起始函数的初始化。

设I是心脏某心相的左心室短轴MSCT的3D数据集，它包含N个层片In：Ω⊂R2→R+，n∈［1，2，…，N］。选取层片Is（s∈［1，2，…，N］）作为序列I的定位模板，其中Is需满足条件：左心室腔面积在层片中较大（腔面积要比主动脉等大血管面积大），左心室腔有较光滑的边界，腔周围有较小的组织噪声。由于左心室腔的灰度值较大，边缘较明显，所以可根据阈值法和形态学操作获得左心室腔的定位初始模板。

1）利用阈值法和形态学模型，去除胸腔的肺、骨头、大血管等组织，初步获取定位模板图像中的左心室腔粗轮廓模板M，即

式中，T1和T2是去除左心室腔以外组织的低阈值和高阈值，*表示二值形态学操作，结构元素SE1用于填充左心室腔中的小孔，Bwselect用于提取阈值分割结果，max用于提取出面积最大的左心室腔。

2）通过层片与模板相乘，提取出左心室腔的粗轮廓ΓtM为

式中，*表示矩阵对应点的值相乘。

沿着左心室的短轴从腔由大变小的方向，使用初始定位模板来获取邻近层片左心室腔的粗轮廓，层片In的腔室掩板可由层片In-1得到，并利用图1所示流程更新定位模板，从而获取整个数据集中层片的左心室腔粗轮廓，完成左心室的定位。

图1 左心室腔粗轮廓的获取流程Fig.1 The flowchart to extract coarse contour of LV cavity

1.3.2 起始轮廓的初始化

水平集曲线演化之前必须确定初始化轮廓的位置、面积，不合适的起始位置、形状和大小可能导致错误的分割结果。在大多数水平集模型中，初始轮廓的形状是圆形、矩形等规则形状，与待分割目标形状无任何关系。本模型则按如下方法、利用左心室腔粗轮廓来定义起始形状，有

式中，Γtn是左心室的粗轮廓，area（Γtn）表示左心室腔的区域，Centroid（ΓLVn-1）表示上一层片左心室的质心，*表示二值形态学运算，结构因子SE3、SE4用于平滑区域轮廓。

利用结构连续性和前一层片的分割结果，可以精确地定位腔室和心肌的起始轮廓，而且由于起始轮廓近似于待分割目标边缘，因此有效降低了水平集演化时间。

1.4 基于先验区域的边缘检测因子

边缘检测因子又称为停止项，文献［8］将其改进为

式中，β是由高斯函数Gσ和引导因子λ卷积构成的矩阵，*表示矩阵对应点的值相乘。

引导因子定义为式中，Θ表示给定的感兴趣区域。

引导因子可以把先验知识融合到曲线的演化当中。由式（11）可知，当曲线演化到感兴趣区域时，引导因子退化为常数，g还原成式（1）的边缘检测因子，曲线保持正常的演化。当引导因子为0时，曲线会越过不感兴趣的区域。由于MSCT图像层片间左心室的结构和位置相对不变，假设Γn-1是上一层片图像分割的左心室腔区域，则下一层片边缘和层片n-1相近。本模型中设Γn-1为感兴趣区域，则分割室腔的水平集引导函数可定义为

同样，可以设置心肌的感兴趣区域。通过重新设计引导函数，水平集演化方程在每个层片上都有对应的检测因子和停止函数，降低了曲线的演化泄露和曲线在目标区域内局部极值上的迭代。

1.5 改进耦合水平集模型

把融合区域先验模型的边缘检测因子、水平集耦合因子加入水平集演化方程，便可得到ICLS模型的演化方程为

步骤1：通过式（7）和式（8），获取层片In的左心室腔粗轮廓。

步骤3：通过式（10）和式（12），更新边缘检测因子gn。

步骤5：计算|φt+1n-φtn|＜θ，θ是自定义的小常数。如果对比水平集演化前后的面积差满足该条件，便停止演化，否则继续步骤4。

2 实验方法

在应用ICLS模型分割前，必须对获取的MSCT图像数据进行有效的预处理。首先对原始CT数据作非线性灰度变换，把一定范围内的灰度值映射到特定范围内，使得在该灰度范围内的图像更加便于后续图像处理［11］。此外，CT数据会受到数据采集系统的电子噪声和重建噪声等的影响。为提高分割质量，必须去除其中的局部噪声和斑点噪声，同时要保留待分割目标的边缘。采用Paris等提出的快速双边滤波器［12］来去除CT图像噪声，与传统的高斯滤波器相比，双边滤波器能较好地保持图像边缘，而高斯滤波后的左心室边缘却比较模糊（见图2）。

图2 滤波效果。（a）高斯滤波器；（b）双边滤波器Fig.2 Filtering results.（a）Gaussian filter；（b）bilateral filter

实验数据采自GE公司的64排螺旋CT，采用心电门控技术获取的3D心脏MSCT数据，分辨率512像素×512像素×297像素；以及来自Philips公司的256排MSCT的8例3D心脏数据集，分辨率为512像素×512像素×331像素。首先选用64排3D心脏数据集，初步验证ICLS模型的分割效果。ICLS模型的参数设置为μ=0.04，λ=5，ν=3，ε=1.5，θ=10，曲线演化步长Δt=5。为了进一步验证算法对三维数据分割的适用性和有效性，保持ICLS模型的参数不变，对8例256排MSCT三维数据集进行左心室心肌和血腔的分割与提取。

为了精确比较并评价自动分割结果与手工分割结果的相似度，用ICLS模型自动提取该数据集的短轴位［13］左心室腔和心肌，然后将自动分割的结果与手工勾勒的结果进行比较。采用如下方法［14］进行相似度（Similarity，SM）对比：设CT层片经计算机分割处理后，得到目标面积S1、手工勾勒轮廓包围的面积S2和共同包围的面积S，相似度的计算公式为

SM的范围为0～1。当S1与S2完全重合时值为1，完全不重合时值为0；当相似度数值越接近于1，表明算法的分割结果越接近于手工分割，算法的准确性就越高。

3 结果

图3是原水平集模型（WRLS）和改进耦合水平集模型（ICLS）起始轮廓的对比，后者起始轮廓已近似于待分割轮廓，并能保证左心室的精确定位。图4是两种模型在相同参数下左心室的分割结果，可见在灰度近似区域ICLS模型减少了边缘泄露。

图3 起始轮廓。（a）WRLS模型；（b）ICLS模型Fig.3 The initial contour.（a）WRLS model；（b）ICLS model

图4 分割结果。（a）WRLS模型；（b）ICLS模型Fig.4 Segmentation results.（a）WRLS model；（b）ICLS model

图5是相似度分析后的结果，其中左心室腔的平均相似度0.950 3，标准偏差0.043 5；左心室心肌的平均相似度0.920 6，标准偏差为0.026 8。在腔和包围腔的心肌之间，灰度对比度较为明显，边缘较为清晰，腔的定位和分割效果较好。心肌和邻近组织的灰度相似度很高，心肌边缘不是特别明显。不过，除极个别层片外，其他层片的分割情况都较好，自动分割与手工分割的结果基本吻合。图6显示的是对分割结果进行3D表面重建后的结果，可以看出分割结果在左心房短轴位具有很好的同一性和完整性。

表1是ICLS模型与分割对心肌的分割结果，用式（14）计算获得平均相似度。可得ICLS模型的血腔分割结果的平均相似度为0.959 2，平均标准偏差为0.019 0；心肌分割结果的平均相似度为0.906 3，平均标准偏差为0.025 6。可以看出，对于不同来源的三维左心室数据，ICLS模型同样具有较高的分割精度，对心肌的平均分割准确率可以达到95％以上，对心肌的分割准确率可以达到90％以上，说明ICLS算法具有普遍的适用性。因此，ICLS模型可以用于提取整个3D+t的MSCT心脏数据，这为后续多心相的心功能分析研究打下良好的基础。

图5 分割相似度。（a）左心室腔；（b）心肌Fig.5 Similaritiy of Segmentation results.（a）LV cavity；（b）myocardium

图6 左心室三维重建结果。（a）不同角度下的心肌；（b）不同角度下的血腔Fig.6 3D surface reconstruction results of LV.（a）myocardium from different views；（b）cavity from different views

4 结论

针对MSCT图像左心室定位与分割的特定应用，本研究提出了融合左心室先验形状和结构连续性的耦合水平集模型——ICLS模型。它首先根据结构连续性和阈值法定位，提取出左心室腔的粗轮廓；然后利用左心室先验知识以及结构连续性，指导后续水平集曲线的演化。ICLS模型对MSCT数据集的实验表明，它能实现左心室腔和心肌的自动分割，与其他基于单幅梯度边缘、无耦合因子、无先验知识的水平集方法比较，结果更准确和完整，还提高了曲线演化速度并减少了边缘泄露，这为进一步分析左心室运动提供良好的数据保证。

表1 ICLS模型分割结果与手工分割结果的相似度与标准偏差Tab.1 The similarity and standard deviation between ICLS model and manual segmentation results

［1］Fritz D，Kroll J，Dillmann R，et al.Automatic 4D-segmentation of the left ventricle in cardiac-CT-data［C］.//Pluim JPW，Reinhardt JM.，eds.Medical Image 2007：Imaging Proceeding.Bellingham：SPIE-INT Soc Optical Engineering，2007，6512：N5123.

［2］Garreau M，Simon A，Boulmie D，et al.Assessment of left ventricular function in cardiacmsct imaging by a 4D hierarchical surface-volume matching process［J］.International Journal of Biomedical Imaging，2006，1-10.

［3］OsherS，Sethian JA.Fronts propagating with curvature dependent speed：algorithms based on the Hamilton-Jacobi formulation［J］.Journal of Computational Physics，1988，79：12-49.

［4］Caselles V，Catte F，Coll T，et al.A geometric method for active contours［J］.Numeriche Meathematik，1993，66：1-31.

［5］Malladi R，Sethian JA，Vemuri BC.Shape modeling with front propagation：a level set approach［J］.IEEE Trans Pattern Analysis and Machine Intelligence，1995，17（2）：158-175.

［6］Chan T，Vese L.Active contours without edges［J］.IEEE Transactions on Image Processing，2001，10（2）：266-277.

［7］Li CM，Xu CY，Gui C，et al.Level set evolution without reinitialization：a new variational formulation［C］.//Schmid C，Soatto S，Tomasi C.，eds.IEEE International Conference on Computer Vision and Pattern Recognition.Los Alamitos：IEEE Computer Soc，2005，1：430-436.

［8］Zhang H，Morrow P，McClean S，et al.Incorporating feature based priors into the geodesic active contour model and its application in biomedical imagery［C］//McDonald J，Markham C，Ghent J，eds.International Machine Vision and Image Processing Conference.Los Alamitos：IEEE Computer Soc，2007：67-74.

［9］Lynch M，Ghita O，Whelam PF.Left-ventricle myocardium segmentation using a coupled level set with a priori knowledge［J］.Computerized Medical Imaging and Graphics，2006，30：255-262.

［10］Liu Junwei.，Feng Huanqing，Zhou Yingyue，et al.A novel automatic extraction method of lung texture tree from HRCT image［J］.Acta Automatica Sinica，2009，35（4）：345-349.

［11］Heinemann EG.Simultaneous brightness induction as a function of inducing and test field luminances［J］.Journal of Experimental Psychology，1955，50（2）：89-96.

［12］Paris S，Durand F.A fast approximation of the bilateral filter using a signal processing approach［C］//Bischof H，Leonardis A.，eds.European conference on Computer Vision.Berlin：Springer-Verlag，2006，81（1）：24-52.

［13］American Heart Association Writing Group on Myocardial Segmentation and Registration for Cardiac Image.Standardized myocardial segmentation and nomenclature fortomographic imaging of the heart［J］.Circulation，2002，105：539-542.

［14］周颖玥，冯焕清，李传富，等.一种从胸HRCT图像序列分割肺的自动化方法［J］.北京生物医学工程，2008，27（1）：6-10.