基于阈值区间的水平集算法在耳蜗分割中的应用

基于阈值区间的水平集算法在耳蜗分割中的应用 基于阈值区间的水平集算法在耳蜗分割中 的应用 第4O卷第2期 2006年2月 浙 Journa1 江大学(工学版)ofZhejiangUniversity(EngineeringScience)Vo1.40No.2 Feb.2006 基于阈值区间的水平集算法在耳蜗分割中的应用 刁现芬,陈思平,梁长虹.,汪元美 (1.浙江大学生物医学工程及仪器科学学院,浙江杭州310027; 2.深圳安科企业博士后工作站,广东深圳518067;3.广东省人民医院影像科,广东广州510080) 摘要:为了获取独立的三维耳蜗模型,提出了一种交互式,半自动,由粗到细,结合三维显示反馈的耳蜗分割方 法.浏览二维切片,调节图像强度区间并选取感兴趣区域,从颞骨螺旋CT(computedtomography)图像中对耳蜗进 行粗分割;采用基于阈值区问的三维水平集(1evelset)分割算法对耳蜗进行细分割.应用可视化技术,对分割结果 进行实时显示,显示结果作为用户进行参数调节的参考依据,经过多次人机交互,直至获得满意的分割结果.利用 临床颧骨螺旋CT图像进行耳蜗分割实验,结果表明该方法适合于分割结构复杂且表面光滑的目标. 关键词:耳蜗;图像分割,水平集算法;可视化 中图分类号:TNgl1.73文献标识码:A文章编号:1008—973X(2006)02—0262—05 Thresh0ld—baSedlevelsetalgorithmforcochleasegmentation DIAOXian—fen,CHENSi—ping',LIANGChang—hong.,WANGYuan—mei (1.CollegeofBiomedicalEngineeringandinstrumentScience,ZhejiangUniversity,Hangzhou310027,ChinaI 2.AnkeEenterprisePost,DoctorWorkstation,Shenzhen518067,China; 3.DepartmentofRadiology,GuangdongProvincialPeople'SHospital,Guangzhou510080,China) Abstract:Toobtainthe3Dseparatemodelofthecochlea,aninteractive,semiautomatic,coarse—to—fineim- agesegmentationapproachwith3Dvisualfeedbackwasdevelopedforthecochleasegmentation.Thein— tensityrangeandthevolumeofinterestweredefinedtoroughlysegmentthecochleabyviewingthespiral computedtomography(CT)imagevolumeslicebyslice.Thenthelhreshold— based3Dlevelsetalgorithm wasusedtoseparatethecochleafromotheranatomicstructures.Visualizationtechnologywasusedtoren— derthesegmentedresults.Theusercanmodifytheintermediateparametersbasedonthereal—time3Dvis— ualfeedbackuntilgettingsatisfyingresult.ThetestresultsusingspiralCTimagesofthetemporalbone provethatthemethodisapplicableforsegmentingcomplicatedstructureswithsmoothingsurface. Keywords:cochlea;imagesegmentation;levelsetalgorithm;visualization 位于颞骨中的耳蜗是一个螺旋状的骨蜗管,形 状类似蜗牛,有豌豆粒大小.耳蜗的主要功能是把声 波的振动转化为螺旋神经节能够感知的电脉冲信 号.许多学者利用计算机图像处理技术对耳蜗进行 了可视化,几何形状分析等的研究.随着影像学 技术的进步,可以获得三维空间上高分辨率的颞骨 图像数据,利用三维医学图像数据获取更为精确的 耳蜗结构信息成为新的需求.图像分割是图像处理 和计算机视觉领域低层次视觉中的主要问题,同时 也是一个经典难题?.通常使用商业软件,利用现有 的区域增长,阈值区间,边界探测,形态学滤波等算 法对耳蜗进行分割.然而,由于以下几个主要难 点常常使得耳蜗的分割并不完整:首先,由于CT (computertomography)图像存在部分容积效应 收稿日期:2005—03—14.浙江大学(工学版)网址:www.journals.zju.edu.cn/eng 基金项目:广东省科技资助项目(20042270021). 作者简介:刁现芬(1976一),女,山东临沂人,博士生,从事计算机图像处理研究.E—mail:diaoxianfen@asisz.corn 量藩l藐{;戢豁嚣}#豁羲孽}}},鞋H黔配豁棼黪瓢堂黔静嚣}瓤?,靠《豁嚣} 第2期刁现芬,等:基于阈值区间的水平集算法在耳蜗分割中的应用263 (partialvolumeartifacts),导致邻接物质边缘分界 同耳蜗邻接的组织与耳蜗有着相似的灰 不清;其次, 度分布;再次,几毫米大小的耳蜗有着复杂的结构, 因此从颞骨螺旋CT图像中分割出耳蜗极其困难. 自从Osher等人_7制提出水平集(1evelset)算法 以来,该方法逐渐被应用于图像分割,图像建模,计 算机视觉等领域.水平集算法是一种描述曲线或曲 面演化的有力工具,其主要优点是可以对任意复杂 的形状进行模型化并且隐式地解决拓扑形状的分 裂,合并等变化.另外,该算法中与曲率相关的速度 项可以控制曲线或曲面的光滑特性,因此对于分割 具有光滑表面的组织该算法更加具有优势.另外,交 互式分割对分割结果进行实时显示,不断地把中间 结果反馈给用户,用户再根据反馈结果调节算法参 数,这为高效地获取分割目标提供了保证.Hohne 等人..认为交互式分割有望成为图像分割的主流 方法.尽管交互式分割方法已经得到了广泛的应用, 然而具有实时三维显示反馈的交互式三维图像分割 受计算机速度的限制,目前应用范围还不是很广泛. Yoo等人_1妇采用区域自适应Snake模型与三维显 示相结合的方法对耳蜗进行了分割,首先对二维切 片图像进行分割,并把得到的分割结果作为其邻接 切片的初始轮廓,然后将二维分割的结果堆叠成三 维数据,并对其进行三维显示,也就是二维分割与三 维显示的结合.然而,这种方法只在相邻图像切片上 组织变化不大时有效.采用三维水平集算法对耳蜗 进行分割,不仅不受这一假设条件的限制,而且缩短 了处理时间,可以获得更加精确的三维耳蜗模型. 本文首先阐述了水平集算法的基本原理及特征 图像的计算过程,提出了一种从人类颞骨螺旋CT 图像中分割耳蜗的半自动分割方法.该方法与其他 交互式分割方法基本类似,不同的是应用了基于阈 值区间的三维水平集算法,并且应用可视化算法对 分割结果进行了实时三维显示. 1水平集算法原理 1.1基本原理 水平集算法[7吲是一种跟踪曲线或曲面轮廓演 化的数学方法,当用于二维图像分割时称为曲线演 化,三维时则称为曲面演化,这里统称为轮廓演化. 本质上,水平集算法是通过水平集函数(X,,)在微 分方程控制下的演化来控制轮廓的演化.(X,,)的 构造方法为:在图像空间上构造一个静态的,间隔均 匀的网格,网格点上的值代表网格点到轮廓的最短 距离,轮廓内部距离为负,外部为正,给定初始轮廓, 就可以计算出在,一0时刻的水平集函数.在一系列 内力和外力的作用下,水平集函数随时问变化,即网 格上各点的值在变化.在任何时刻,,满足(X,,)一 0的网格点的集合,即函数(X,,)的零水平集(zero levelset)就构成了当前的传播界面,如图1所示.由 此可见,在轮廓演化过程中拓扑结构的任意变化都 隐含在了水平集函数的演变之中. 图1零水平集的概念 Fig.1Conceptofzerolevelset 根据链式规则,得到轮廓演化方程为] +Fll=o.(1) 式中:为梯度算子;F为沿轮廓法线方向的速度 函数,控制着轮廓的膨胀和收缩,一般F包括与图 像相关的速度分量(如梯度信息等)以及与曲线的几 何形状相关的速度分量(如轮廓的曲率等).常用的 速度函数一般可以描述为 F=o?P(X)+p?K(X).(2) 式中:P(X)为与图像相关的传播速度,K(X)为与 曲率相关的速度,a,口为两者的权重系数.传播速度 P(X)使得轮廓在灰度值连续的区域快速扩张,而在 边界处停下来;曲率速度K(X)则用来控制轮廓的 光滑情况,大小与曲率成正比.在光滑的边界处,F 将趋于0.用?,表示曲面演化的时间步长,在三维 情况下式(I)的离散形式为 +F.1l-O.(3) 即 一?,?F?lI.(4) 式中:表示第n个At时刻,图像网格中位于(i, .『,)处的点的距离函数值;At?F?ll表示经 过一个?,时间间隔符号距离值的变化量.式(4)即 为轮廓演化的离散迭代公式.三维网格上计算曲率 K(X)的公式为 一 264浙江大学(工学gt)第4O卷 [虹(+)+(+)+(+)一 2一2:一2J/ (以++).(5) 在算法的实现过程中,需要把图像灰度值或梯 度值应用符号映射函数变换到Eo,1.0]或 [一1.0,1.0]区间,生成所谓的特征图像,并作为水 平集算法的一个输入.当传播速度P(X)的大小与 特征图像成正比,并且特征图像的像素值大于零时, 垂直轮廓向外,反之,向内. 与其他曲线演化方法相比,水平集算法最大的 优势在于它的稳定性以及拓扑无关性.把图像分割 问题转化为求解水平集函数在微分方程控制下随时 间的演化问题,看上去是将二维的问题扩展到了三 维,将三维的问题扩展到了四维,增加了计算复杂 度,但在问题的求解上带来很多优点其最大的优点 就是隐含地解决了轮廓的拓扑变化问题.另外,水平 集算法还具有易于扩展到三维的优点. 1.2特征图像的计算 通常采用的特征图像有基于边界的特征图像和 基于阈值区间的特征图像.基于边界的特征图像与 像素值的大小无关,由相邻像素间的差异大小和符 号映射函数决定;基于阈值区间的特征图像,由像素 值的大小和符号映射函数决定.特征图像决定了与 图像相关的速度分量,对分割结果的好坏有着非常 重要的影响. 基于边界的特征图像的获取过程如下:1)图像 平滑,比如采用高斯滤波;2)计算图像梯度幅值,梯 度幅值的大小反映了像素间的连接紧密程度;3)构 造符号映射函数或采用线性变换,把梯度幅值变换 到[0,1.0]区间,最高的梯度幅值映射到0,最低的 梯度幅值映射到1.0. 基于阈值区间的特征图像的获取过程可以用图 2简单地描述.像素值位于[1ower,upper-]之间,映 射到[0,1.0]区间,称为前景区域;位于区间之外的 像素值映射到[一1.0,0]区间,称为背景区域.基于 阈值区间的水平集算法可以看作是阈值连接分割算 法与水平集算法的结合.使用该算法时,用户需要给 定初始轮廓,设置阈值区问和选择具体的符号映射 函数等. 由于基于边界的特征图像位于[O,l_O]区间,而 传播速度与特征图像成正比,就要求设置的初始轮廓 完全位于感兴趣目标中,因为对于该种情况轮廓只能 膨胀,不能收缩;而基于阈值区间的特征图像的像素 值位于[一1.0,1.O]区间,轮廓既可以膨胀也可以收 缩,因此初始轮廓可以部分地位于感兴趣目标之外. 前景 图2基于阈值区间的符号映射函数 Fig.2Threshold—basedSigmoidmappingfunction 2耳蜗图像分割过程 目前,全自动分割还不能取得令人满意的结果, 本文的目的就是要最大化地利用分割算法的自动处 理功能,使手动操作的工作量最小化.对于耳蜗分割 最为重要的知识依据就是耳蜗的表面非常光滑,之 所以选择水平集算法对耳蜗进行分割,很重要的一 点就是该算法可以控制轮廓的光滑度.采用由粗到 细的分割策略,在粗分割阶段提取感兴趣区域,减少 了细分割的计算量,从而提高了分割速度;在细分割 阶段应用基于阈值区间的水平集分割算法对粗分割 的结果进行改善,最终得到更为精确的分割结果.在 图像分割的两个阶段,均应用可视化技术对分割结 果进行显示,直观的显示结果在整个分割过程中为 用户进行参数调节提供了重要的参考依据. 2.1粗分割 粗分割过程主要是从复杂的颞骨结构中找到耳 蜗,完成感兴趣区域(volumeofinterest,VOI)的提 取和图像插值以及定位出耳蜗所在的强度区间.浏 览颞骨螺旋CT图像切片并利用平行光线投影 通 (parallelraycasting)可视化算法进行三维重建,过不断调节强度区间,不透明度传递函数,颜色映射 函数等,直至耳蜗可见为止,结合3个正交截面(横 断面,冠状面和矢状面)确定出耳蜗所在的区域.对 于断层图像层间距与像素间距不均匀的情况,为了 使水平集分割算法和可视化处理算法更为有效且便 于数字实现,把提取出的VOI区域,采用线性插值 的方法,构成立方体体元的形式.对插值后的V0l 图像数据进行三维显示,通过调节颜色映射函数,不 透明度传递函数,直至耳蜗区域清晰可见,并旋转显 示结果来观察所选择的区域大小是否恰当.也可以 重复以上过程,重新选定感兴趣区域,目的是使所选 择区域更为有效.运用平行光线投影可视化算法绘 制粗分割的结果,如图3所示.显然这一显示结果不 能满足耳蜗研究的需要,更是无法满足人工电子耳 第2期刁现芬,等:基于阈值区间的水平集算法在耳蜗分割中的应用265 图3租分副结果 Fig.3Resultofcoarsesegmentation 蜗植入的术前检查和术后评估的需要.为了满足临 床应用的要求,需要运用更为有效的图像分割算法, 实现感兴趣目标的更为精确的分割. 2.2细分割 在细分割阶段采用基于阈值区间的水平集分割 算法对粗分割的结果进行改善,完成从复杂的颞骨 组织中分割耳蜗的任务.水平集分割算法需要2个 输入:特征图像和初始轮廓.为了便于用户操作,采 用球面作为初始轮廓,用户只需给定球心的位置及 半径就可以完成轮廓的初始化工作.在水平集分割 算法中,轮廓的分裂与合并都隐含在了水平集函数 的演变中,因此可以在不同的位置添加多个初始轮 廓.根据图像灰度值设置强度区间,把耳蜗所在的区 域划为前景,周围的组织划为背景,生成特征图像. 设置传播速度和曲率速度的权重系数,用于速度函 数的计算.在算法迭代过程中,可以通过浏览3个正 交截面方向上的切片来观察轮廓的演变过程,也可 以随时对分割的中间结果进行三维显示,但这将影 响算法的执行速度.如果中间结果就不满足要求,即 可终止本次迭代,重新设置初始轮廓,调节参数后, 再次启动迭代过程.迭代在3种情况下会终止:分割 结果达到要求,用户强行停止迭代过程;迭代次数到 达设置的最大迭代次数;迭代次数增加而轮廓不变 或者变化非常小.在细分割阶段使用了步进立方体 (marchingcubes)可视化算法进行三维显示,这是 由于分割结果中采用了某个整数来标记轮廓及其内 部所有的网格点. 3实验结果 使用GE螺旋CT设备获取一组人体颞骨影像 数据,图像空间分辨率为512×512×65,像素物理 间隔为0.273mm×0.273mm×0.625mm,并且采 用2字节表示1个像素.在粗分割阶段,经过感兴趣 区域提取和线性插值得到物理间隔为0.1mm×0. 1mm×0.1mm,空间分辨率为170×186×182的 预处理图像,然后再对其进行细分割.图4(a)为原 始图像切片之一;图4(b)为提取的感兴趣区域经插 值后并且当阈值区间设为[o,1385]时所生成的特 征图像中的轴向切片,实验中把此时的特征图像作 为水平集算法的输入.图5(a),(c)为设置了3个 初始轮廓后抽取的交于一点的3个正交截面;图5 (d),(f)为迭代172次后的3个正交截面,由图可 以看出在轮廓演变过程中进行了合并.图6(a),(b) 为迭代172次后从两个方向上观察的三维显示结 果,轮廓在耳蜗部分扩张,并未遍历所有的耳蝇区 域;图6(c),(d)为迭代374次后在两个方向上的三 维显示结果,此时轮廓已经遍历了耳蜗的所有区域, 而且已经扩展到周围的组织中.在整个实验中设置 传播权重—1.0,光滑系数一0.57.图中虚线为 标记,z坐标用,三维时的交点为图像的中心. 图5给定初始轮廓后和迭代172次后的正交截面图像 Fig.5Orthogonalplanes(axial,coronal,saggita1) withinitialcontoursandafter17Ziterations 266浙江大学(工学版)第4O卷 囵目 囵圈图6迭代172次和374次后的分割结果 Fig.6Segmentationresultsofl72and374iterations 4结语 将分割算法同可视化技术相结合的半自动分割 系统,由于可以提供三维显示结果用于指导分割过 程,不仅减少了手动分割的工作量,而且提高了分割 效率.由粗到细的分割过程,减少了计算机的工作 量,也使自动分割算法更为有效,结果更加精确.实 验结果表明,基于阈值区间的三维水平集分割算法 非常适用于耳蜗的分割,但由于耳蜗同周围的组织 有着相似的灰度分布,当迭代达到一定次数时,轮廓 将扩展到周围的组织.由于耳蜗结构复杂,通过限定 迭代次数并不能根本解决这一问题,而将初始轮廓 设置在远离耳蜗周围组织的地方,结果会得到一些 改进.另外,设置初始轮廓,参数等的随意性较大,获 得同样的结果需要的迭代次数也不一样,并且随着 轮廓的增加,迭代一次所需要的时间也有所增加. 在粗分割的显示结果中,在部分区域内可以看 到前庭阶和鼓阶的模糊分界线,而经过细分割之后, 虽然耳蜗的轮廓清晰.可见,但未能看到前庭阶和鼓 阶的分界线,这是由于为了使轮廓在演变过程中保 持光滑而平滑掉了耳蜗组织内部的弱边界所致.尽 管如此,采用本文方法还是得到了耳蜗区域的精确 数字描述,在此基础上进行更为细致的分割,分割出 前庭阶和鼓阶,将是接下来的研究工作之一. 参考文献(References): [1]KAWANOA,SELDONHL,CLARKGM.Comput er-aidedthreedimensiona1reconstructioninhumanco ehlearmaps:MeasurementofthelengthsoforganofCorti, outerwall,innerwall,andRosenthal'scanal[J].An— naisofOtology.Rhinoiogy,andLaryngology,1996, 105:701—709. r2]WANGG,VANNIERMW,SKINNERMW,eta1. UnwrappingcochlearimplantsbyspiralCT[J].IEEE TransactionsonBiomedicalEngineering,1996,43(9): 891—9OO. r3]HIMIT,KATAURAA,SAKATAM,eta1.Three~ dimensionalimagingofthetemporalboneusingahelical CTscananditsapplicationinpatientswithcochlearim— plantation口].JournalofOto-Rhino-Laryngoiogy&Its RelatedSpecialties,1996,58(6):298—300. [4]庞全,苏佳,段会龙.基于四叉树和交叉熵的面向对象图 像分割方法[J].浙江大学:工学版,2004,38(12): 1618. 1615— PANGQuan,SUJia.DUANHui—long.Object—oriented imagesegmentationmethodbasedoncross—entropyand quad—tree[J].JournalofZhejiangUniversity:Engineer— ingScience,2004,38(12):1615—1618. [5]KETTENDR,SKINNERMW,WANGG,etaLIn vivomeasuresofcochlearlengthandnucleuscochlear implantarrayinsertiondepth[J].AnnalsofOtoiogy, Rhinoiogy,andLaryngoiogy,1998,175:1—16. 1-6]YOOSK,WANGG,RUBINSTEINJT,eta1.Three~ dimensionalgeometricmodelingofthecochleausing helico—spiralapproximation[J].IEEETransactionson BiomedicalEngineering,2000,47(10):13921402. [73OSHERS,SETHIANJA.Frontspropagatingwith curvaturedependentspeed:AlgorithmsbasedonHamil— ton—Jacobiformulation[J].JournalofComputerPhys一 ' ics,1988,79(1):12—49. [8]MALLADIR,SETHIANJA,VEMURIBC.Shape modelingwithfrontpropagation:Alevelsetapproach[J]. IEEETransctionsonPA??,1995,17(2):158—175. [9]HOHNEKH,HANSONWA.Interactive3Dsegmen tationofMRIandCTvolumesusingmorphologicalop— erations[J].JournalofComputerAssistedTomography, 1992,l6(2)1285—294. [1O]SALVIROONPORNP,ROBATINOA,ZAHAJSZKYJ, eta1.Real—timeinteractivethree-dimensionalsegmentation [J].AcademicRadiology,1998,5(1):49—56. [11]YO0SK,?rANGG,RUBINSlTEINJ,eta1.Semi-au tomaticsegmentationofthecochleausingreal-timevolume renderingandregionaladaptivesnakemodeling[J].Jour— nalofDigitadImaging,2001,14(4):173—181.