医学图像分割中的分段常数水平集方法与Mumford - Shah模型研究-

1、 摘要视网膜视神经疾病(如青光眼、糖尿病视网膜病变、老年性黄斑病等是由多种原因引起的视网膜及视觉神经组织病变,并最终导致失明。该疾病在世界范围内有较高的发病率,且早期大多无临床征兆,因此早期的诊断十分重要。其中,视(神经乳头或视盘(optic nerve head, or optic disk的定位、分割、重建与度量是近三十年来视网膜视神经疾病计算机辅助诊断研究的关键技术难题。由于光照不均、血管遮挡等原因,使得对象轮廓拓扑结构复杂多变,现有的大多数方法的分割结果不能满足临床诊断需要。本文重点论述如何将基于分段常数水平集方法(PCLSM的Mumford-Shah 模型应用于视乳头图像杯盘分割中。

2、该方法的特点是:(1它能很好地处理多相图像分割问题;(2只需要一个水平集函数就可以分割多个目标;(3处理过程不需要人工干预,自动化程度较高。本文首先采用基于PCLSM的Mumford-Shah模型来分割视乳头图像的杯盘,从而得到杯盘形态边缘,但这个边缘存在断裂、凹陷;因此,本文进一步结合青光眼视乳头图像杯、盘的先验知识,剔除断裂、凹陷的边缘,提取视杯和视盘边缘真实的特征点;然后,采用光滑样条曲线拟合技术重建被血管遮挡的视杯和视盘部分边缘;最后本文用拟合好的杯盘边缘来重建视乳头杯盘图像,并根据该重建图像自动地提取杯盘比等具有青光眼辅助诊断重要价值的病理特征参数值。不同时期青光眼病人的视乳头图像杯

3、盘重建、分割与度量实验结果表明,该方法能克服噪声污染、血管遮挡、光照不均匀、对比度小、个体间差异大等视网膜图像分割中固有的困难,并能有效重建、分割与度量早、中、晚期和典型欧洲青光眼彩色视乳头图像中的视杯和视盘。与临床结论以及多层Mumford-Shah模型的分割结果对比研究表明基于分段常数水平集Mumford-Shah模型的视乳头杯盘分割、重建与度量的方法是非常有效的。关键词:分段常数水平集方法(PCLSM;Mumford-Shah模型(M-S模型;视乳头图像分割;青光眼;医学图像处理AbstractRetinal optic nerve diseases (such as glaucoma,

4、 diabetic retinopathy, age-related macular disease, etc. are caused by a variety of retinal and neural tissue lesions, ultimately, it leads to blindness. These diseases in the world have high incidence rate, and most of them have no clinical signs at early stage. Therefore, treatment in early stage

5、is very important. The segmentation, reconstruction, and measurement of optic nerve head or optic disc are most challenging problems as for the computer-aided diagnosis of retinal optic nerve diseases in last three decades. However, due to non-uniform illumination, vascular occlusion and other reaso

6、ns, most of existing methods segmentation results may not be satisfying.We mainly focuses on how to use PCLSM based Mumford-Shah modelto segment the images of the optic nerve head. There are some reasons to choose the model: 1, it deals with the multi-phase image well; 2, it segments the image by on

7、ly using a single level set; 3, and there is no need of manual intervention, so it has high degree of automation.Firstly, we can get the cup and the disks morphological edge from segmenting the image of optic nerve head by using the PCLSM based Mumford-Shah model; Secondly, we removed the depressed

8、and broken brink, and retained the edge-points of the cup and disk by combining of the priori knowledge of glaucoma optic nerve head images on cup and disk; And then, we used the technique of smoothing spline curve fitting to reconstruct the edge of the cup and disk obscured partially by vasculars.

9、Finally, we reconstruct the optic nerve heads cup and disk image by using the fitted edge of the cup and disk, and extract the pathological measurement such as cup-disk ratio based on the reconstructed image.The results of the reconstruction, segmentation and measurement of the optic nerve head imag

10、es in different stages demonstrated that this method is able to effectively segment, reconstruct and measure the optic cup and disk in optic nerve head images of different stages glaucoma patients, which are of poor quality, very low contrast, obscure due to blood vessels, and distinct inter-differe

11、nces of individuals. We found that the PCLSM based Mumford-Shah model is very effective by comparing with the clinical conclusions and the segmentation results of Multi-layer Mumford-Shah model.Key Words: Piecewise constant level set methods (PCLSM; Mumford-Shah model(M-S model; Segmentation of opti

12、c nerve head; Glaucoma; Medical image processing目录摘要ABSTRACT第1章绪论 (11.1医学图像分割研究现状 (11.2视乳头图像杯盘分割概述 (41.2.1 研究背景 (41.2.2 现有研究方法 (61.3本文研究内容 (91.4本文的内容安排 (9第2章水平集方法 (112.1引言 (112.2水平集方法概述 (112.3水平集方法的数值实现 (142.4水平集函数演化的快速算法 (162.4.1 窄带法 (162.4.2 快速行进法 (172.5分段常数水平集方法 (17第3章MUMFORD-SHAH模型简介 (213.1引言 (2

13、13.2 经典的M UMFORD-S HAH模型 (223.3简化的M UMFORD-S HAH模型 (233.3.1 Chan-Vese简介 (233.3.2数值解法 (243.4多层M UMFORD-S HAH模型 (263.4.1 模型简介 (263.4.2 水平集逐层迭代算法 (27第4章基于分段常数水平集方法(PCLSM的MUMFORD-SHAH模型 (294.1引言 (294.2基于PCLSM的M UMFORD-S HAH模型 (294.3 基于PCLSM的M UMFORD-S HAH模型的数值解 (314.4 算法分析 (33第5章视乳头图像杯盘分割、重建和度量 (385.1引言

14、 (385.2视乳头图像杯盘分割、重建和度量 (415.2.1 视乳头图像杯盘形态分割与重建 (415.2.2 杯盘边缘特征点的提取 (435.2.3杯盘边缘曲线拟合重建 (445.2.4 杯盘度量 (455.3实验研究 (475.4对比研究 (485.5结论 (51总结与展望 (54参考文献 (56附录A 硕士阶段本人发表的论文 (61附录B 硕士阶段参与的项目 (62致谢 (63第1章绪论1.1医学图像分割研究现状图像分割在自动图像处理领域中一直处于中心地位,一方面,它不同于输入输出都是图像的大多数图像处理技术,图像分割是输入的是图像,而输出是从这些图像中提取出其属性的处理方法;另一方面,

15、图像分割是一种基本的计算机视觉技术,是从图像处理到图像分析的关键步骤。因此,有效、合理的图像分割能够为基于内容的图像检索、对象分析等抽象出十分有用的信息,从而使得更高层的图像理解成为可能。医学图像分割是图像分割的一个分支,它主要指将人体某部位图像尽量按照该部位组织的实际情况分割开来加以标识。医学图像分割是高层次医学图像分析、理解和处理的前提条件,它直接影响到后续的分析、理解和处理工作,因而得到广泛的关注和研究。由于医学图像的分割结果直接影响到诊断、放疗计划的制定,因此比起普通的图像分割,其关键在于分割的准确性,即与实际人体组织的符合程度。虽然国内外有许多大专院校和医院的从业人员都在进行组织分割

16、的研究,提出了许多有效的方法,也有许多软件问世,但分割结果总存在着些许遗憾。目前较为准确并为大家广泛接受的组织分割方法,仍然是通过有丰富解剖知识和临床诊断经验的医生,经过专门训练,在图像上用鼠标直接勾勒出组织结构的边界。当然这个很费时费力,并且取决于分割者的经验和水平。而且对于组织分割的结果评价,目前仍没有一个统一标准,较常用的评估方法有专家目测、体模验证、蒙特利尔神经研究所推出的人体计算机化解剖图谱以及近年出现的可视人体数据库(visible Human Database,vHD,一些专家倾向于用vHD作为组织分割的金标准(gold standard 1。在显微镜发明之后,1895年伦琴发现

17、X射线促使了医学图像第二次重大发展,他使得人们可以无创的探测物体内部,能够避免和减轻患者进行诊断所要承受的痛苦,改变了临床诊断的手段,开创了医学诊断的新纪元。近年来,随着影像技术的发展,现代医学已经越来越离不开影像技术。在现代化的医院里,包含了各式各样的医学成像系统,这些系统被用于临床的很多方面。这些医学成像系统按其成像模式可分为X-线(X照相技术,CT(X线计算断层扫描,MRI(核磁共振,DSA(数字减影血管造影,PET(正电子发射断层扫描,SPCET(单光子发射断层扫描,Ultrasonic Imaging(超声成像,Endoscopy(内窥镜,DIHS(组织切片的数字照相,Thermal

18、 Imaging(热成像,Microscopy(显微镜等.它们提供了丰富的人体2D/3D/4D医学图像,为医生做出正确的诊断和精确的治疗方案提供了科学的方法和平台2,3。随着计算机技术的发展,医学领域也发生了极大的变化,正朝着医疗信息化、智能化的方向发展。计算机图形图像技术与医学影像技术的结合,从根本上改变了医务人员进行诊断的传统方式,例如,在手术过程中,不断获得术中数据,如复合模式图像数据、定位系统的定位坐标、组织器官的位置等,有时候要进行必要的三维显示。在这些应用中,医学图像分割往往是十分关键的一步,是进行医学图像对象提取、定量分析、三维重建等等处理的一个必不可少的关键步骤。比如,如何将心

19、脏瓣膜分割开来,以便可以研究瓣膜的运动特性;如何将脑部CT切片的灰质和白质分割开来,以重构整个脑部的三维结构;如何将病灶部位与健康组织分割开来,以研究病灶的大小和位置等。因此,为了能够准确进行医疗诊断与分析,对目标的边缘和区域进行准确分割显得尤为重要。由于医学图像表现为低信噪比、低对比度、组织特征的可变性,以及不同组织或病灶与背景之间的边界的模糊性,特别是细微结构如血管神经等分布的复杂性,医学图像分割的智能化受到很大的限制,目前还没有一种普遍的方法来完成医学图像分割,临床上大多还处于手动分割,少部分则是半自动分割,费时费力,而且其结果受人的主观因素的影响也较大,因此,提高医学图像的分割自动化程

20、度以及分割的精度是医学图像分割研究一直追求的目标。图像分割是医学图像处理的重要内容之一,也是一个经典难题,从20世纪70年代起,图像分割问题就吸引了很多研究人员为之付出努力。图像分割技术在医学图像处理上有广泛的用途,从成像模式上讲有X-线成像图像处理、C-T成像图像处理、MRI成像图像处理等等,从检测部位来讲,可以分为骨肌、腹盆、头颈、心胸等等。图像分割算法一般是基于亮度值的两个基本特性之一:不连续性和相似性.第一类性质的应用途径是基于亮度的不连续变化的分割图像,比如图像的边缘。基于边缘的分割方法有:并行微分算子、曲面拟合法、基于边界曲线拟合的方法、串行边界查找,和基于变形模型的方法(这种方法

21、又分为两种:参数变形模型方法和几何变形模型方法等。第二类的主要应用途径是依据事先制定的准则将图像分割为相似的区域。基于区域的分割方法有:门限处理、区域生长(Region Growing、区域分裂合并(Split and Merge、分类器和聚类、基于随机场的方法,和其它基于统计学的方法等都是这类方法的实例4。现在还有结合了边缘和区域的方法,比如结合了参数变形模型方法和聚类的方法。医学图像分割具有一般的图像分割的共性,它要满足:(a图像分割要完整,即要对图像中所有像素进行分类;(b两个不同区域没有交集,也即任意一个像素不同时属于两个或两个以上的区域;(c每个分割区域(或类按照某种属性具有均匀性(

22、即相同或相似,这种属性可以是各种图像特征(如灰度、纹理等或效应的某种语义的描述;(d分割后的每个区域之间绝不重合。传统的分割方法有基于边缘的canny方法5、基于区域的Markov方法6,但由于这些方法在模糊边缘、断续边缘、低信噪比、低对比度的处理上存在很多不足,难以满足医学图像分割的要求,因此有必要寻找新的分割方法来解决这些问题。基于形变模型(Deformable Model 7的图像分割方法是近年发展起来的新方法,该方法综合了各种对图像数据的解释和人们的知识,因此,更接近于人类的视觉机理。形变模型包括活动轮廓模型,刚性联接模型和点分布模型。常见的图像分割模型有:基于边界的模型,基于区域的模

23、型,基于边界和区域的分割模型,结合形状约束的分割模型等。基于边界的分割模型是以图像的边界信息来控制变形曲线的运动速度。对于几何活动轮廓模型,Malladi利用图像的边缘强度来控制曲线演化的速度,使曲线停止在图像边缘。该方法对有好的对比度的图像分割效果不错,但当边缘不明显或存在缝隙,则曲线有可能发生边缘泄露现象。Caselles和Yezzi等则从能量最小化角度提出了另一种几何轮廓线的速度函数,叫做测地线活动轮廓模型(Geodesic Active Contour Model 8,9,用以解决Malladi模型的缺陷。该模型对边缘泄露现象有很好的抑制作用,但对大缝隙仍然无能为力。由于基于边界的分割

24、模型的运动方向固定,一旦形变曲线突破了图像的边界就很难返回,从而导致了分割失败。为解决此类问题,很多学者提出了很多方法,如Zhu、Yuille等提出的区域竞争10的思想。近年来基于活动轮廓的一个优秀的模型是由Mumford和Shah提出的结合图像边界和区域的分割模型,即Mumford-Shah模型1113,该模型不需要待分割模型的先验知识。该模型将图像的恢复、区域分割和边缘检测综合到一起,是一个比较完美的模型。之后Chan和Vese根据Mumford-Shah模型又提出了该模型的一个简化模型,也即Chan-Vese模型1416。这两个模型是近年来基于活动轮廓模型中研究最为广泛的模型。Mumfo

25、rd-Shah模型和Chan-Vese模型固然十分优秀,但是他们的求解却颇为困难。近十年来,水平集方法(level set method 1720获得了广泛的研究。该方法基于曲线演化的思想,把N维对象推演到N+1维进行计算,极大推动了非参数化的几何轮廓线模型的研究。将水平集方法和这两个模型相结合,在模型求解中引入曲线演化的思想,取得了非常好的效果。1.2视乳头图像杯盘分割概述1.2.1 研究背景视网膜视神经疾病(如青光眼、糖尿病视网膜病变、老年性黄斑病等是由多种原因引起的视网膜及神经组织病变,并最终导致视功能完全丧失的眼科疾病。该疾病在世界范围内有较高的发病率,且早期大多无临床征兆,因此早期的

26、诊断十分重要。目前,在临床上最常采用的视网膜图像是通过彩色立体眼底视网膜照相机获得的。大量临床病例证明,致盲性眼病早期会在视网膜上有所反映,因此大多数致盲性眼病可通过视网膜(图像进行分析诊断,因此有效的视网膜图像重建与分割方法一直是该研究领域的热点问题2129。在青光眼的诊断中,具有特别分析诊断价值的是视网膜图像中视盘和视杯的形态及其度量(如杯盘比等,因此,如何有效地对视网膜图像中的视神经乳头部分图像进行视盘和视杯的重建、分割与度量是成功开展计算机辅助青光眼诊断的关键。但由于噪声污染、血管遮挡、光照不均匀、对比度小、个体间差异大等视乳头图像分割中固有的原因,使得视(神经乳头或视盘(optic

27、nerve head, or optic disk的定位、分割、重建与度量是近三十年来视网膜视神经疾病计算机辅助诊断研究的关键技术难题3037。视网膜和视神经相关致盲性眼病主要包括青光眼(Glaucoma、糖尿病视网膜病变(Diabetic Retinopathy-DRP、年龄相关性黄斑变性(Age-related Macular Degeneration-AMD三大类,如图1.1所示。 a青光眼b糖尿病视网膜病变c年龄相关性黄斑变性图1.1 常见的视网膜、视神经病统计调查表明视网膜视神经病患者的数目是惊人的,据报告全世界范围内青光眼、糖尿病视网膜病变、年龄相关性黄斑变性、视网膜肿瘤四大致盲性

28、眼病的疑似患者数占总人口的比例高达20%。其中青光眼患者达1%,但早期可疑高眼压患者高达9%;年龄相关性黄斑变性患者占人口总数的7%左右;目前糖尿病发病率为4%,而当患病时间长于10-15年基本上都会伴有糖尿病视网膜病变的发生。我国的统计情况基本与此接近,以青光眼为例,据北京顺义、广东斗门及新会等地区的调查资料显示,非选择人群的原发性青光眼(不包括先天性青光眼患病率为0.52%,以全国13亿人口计,即有676万原发性青光眼患者;而近年来,国外学者根据蒙古和新加坡的流行病学调查资料推算出我国40岁以上人口中青光眼患者高达940万,其中致盲患者达520万!这些致盲性眼病在中国如此高的患病率,对我国

29、人民视功能损害造成了极大的威胁,不仅对患者本身造成极大的痛苦,而且还造成社会经济的巨大负担和人力资源的重大损失,已成为本世纪我国社会和经济发展的一项沉重负担。由于早期视网膜视神经病体现在相关的组织损伤,大多无临床征兆,待病人发现视功能缺陷(视野缺损时,往往已经产生不可逆转的严重视网膜视神经损害并最终引起致盲,丧失了最佳诊疗机会;而早期诊断可及时地对病人进行有效的治疗,从而可避免严重视野缺损,维持正常的视功能。因此对这些致盲性眼病,特别是青光眼和糖尿病视网膜病变的早期诊断具有非常重要的社会及经济意义。为了有效地解决类似彩色视乳头杯、盘形态重建、分割与度量等众多类别生物医学图像分析与理解(计算机视

30、觉的技术难题,同时为普通医疗人员提供科学、先进的视网膜视神经病诊疗手段,为眼科专家提供辅助诊断决策支持,最终在最大程度上实现青光眼等致盲性眼病的早期诊断以达到预防致盲的目的。计算机图像分析是一项迅速发展的技术,对于视网膜图像获取、保存、处理分析和应用有其独特优势,有利于提高视神经组织定量测量的准确性和可重复性。目前视网膜图像可由眼底(立体摄像机、各种激光(断层扫描设备、血管造影设备和超声波成像设备等获取,获取的信息可以是图像或者其他数字化数据(如视网膜组织断层扫描,如图1.2所示:目前,在临床上最常采用的视网膜图像是通过彩色立体眼底视网膜照相机获得的,但从目前公开发表的相关文献可知,在运用计算

31、机进行视网膜图像重建与分割时仅利用了彩色视网膜图像中的亮度或绿色通道等单通道信息,这使某些本可解决的图像分割问题无解。而在青光眼的诊断中,具有特别分析诊断价值的是视网膜图像中视盘和视杯的形态及其度量(如杯盘比等,因此,如何有效地对视网膜图像中的视神经乳头部分图像进行视盘和视杯的重建、分割与度量是成功开展计算机辅助青光眼诊断的关键。但由于噪声污染、血管遮挡、光照不均匀、对比度小、个体间差异大等视乳头图像分割中固有的难题,使得视盘和视杯的重建、分割与度量分析非常困难。视(神经乳头或视盘和视杯(ONH:optic nerve head, or OD:optic disk,OC:optic cup的定

32、位、分割、重建与度量一直是近三十年来视网膜视神经病计算机辅助诊断的关键技术难题。为了解决如上所述的关键技术难题,由湖南大学王耀南教授参与的重大国际合作欧盟第五框架项目“欧洲青光眼预防工程/计算机辅助青光眼预防诊断项目(EGPS project /Glaucoma Prevention by Computer Aided Diagnostics (GlauCADNo QLG1-2000-00651”已取得大量的相关研究成果。目前,由湖南大学刘国才副教 图1.2 多源视网膜图像授负责的国家自然科学基金项目(复杂图像处理中的自由非连续问题及其水平集方法研究,基金号为60872130和两个湖南省自然科

33、学基金项目(基于水平集方法的视网膜图像分割与度量,基金号为06JJ4044,以及生物医学图像处理中的自由非连续问题及其水平集方法研究,基金号为09JJ3118在以上国际合作项目研究的基础上继续深入地开展这方面的研究工作。视乳头是最重要的视觉解剖组织,它是视觉神经和眼部血管的出、入口,其病变直接导致视野的丧失,甚至致盲。视乳头也是定位其它重要视觉解剖组织(如黄斑、中央凹等的重要标志组织。正常人眼的视乳头外观近似为直立的椭圆,长轴弱长于短轴,面积约为2 mm2。在临床上最常采用的彩色立体眼底视网膜图像中,正常人眼的视乳头表现为最亮的黄白色、近似椭圆的轮廓,深红色的主动脉和静脉血管从中穿过,将其分为

34、多个不规则的小区域。通过对视乳头的的定位、分割、重建与度量之后可以得到丰富的病理数据,从而为进一步诊断提供科学定量的依据。1.2.2现有研究方法视网膜血管脉络和视乳头的视杯视盘是视网膜和视神经疾病的两个非常重要的诊断对象之一,至今为止,对于视网膜图像的图像分割也主要集中在这两个方面进行研究。视网膜血管脉络的提取有很多方法,例如采用2维方向匹配滤波器分割方法38、采用曲线结构模型来作为一个确认步骤的多阈值探测方法39、基于曲率评估的数学形态学滤波40、模糊c-均值聚类(Fuzzy c-means Clustering和直方图匹配(Histogram Matching方法41、基于高斯混合模型的贝

35、叶斯分类方法42、脊探测43法、Morlet小波(Morlet Wavelet变换法44等等。其中的2维方向匹配滤波器分割方法和采用曲线结构模型来作为一个确认步骤的多阈值探测方法的运算简单,但是对于对比度小、噪声点多以及眼部损伤的干扰,这些算法常常导致分割不理想甚至失败;模糊c-均值聚类和直方图匹配方法着重在于解决对比度小和眼部损伤可能带来的干扰,但是该方法不是基于全局信息的,去噪的能力较弱;脊探测方法则需要人工指定血管脉络的起点和终点,自动化程度不高;基于高斯混合模型的贝叶斯分类方法以及Morlet小波(Morlet Wavelet变换法的分割结果较为理想,而且能够适用于多种条件的视网膜图像

36、,尤其是Morlet小波(Morlet Wavelet变换法,由于Morlet小波能够较好的去噪以及增强血管脉络,使得其性能更加优异,但是其中的参数选择很麻烦,因此可以结合神经网络的方法可以通过样本来训练参数,样本选择越好,分割的结果也就越理想,然而这样也表现了该方法的缺点,当训练样本选择不好的时候,分割的结果也不理想。人们提出了多种方法来进行视乳头的视杯视盘分割。现有的视乳头定位、分割、重建与度量方法主要基于参数化的活动轮廓(蛇形模型( active contour model, or snakes24、圆形Hough变换31、活动形状模型或点分布模型(active shape model,

37、 or point distribution model32、分水岭变换(watershed transformation33、小波变换(Wavelets45、完全匹配法46,47、水平集方法48、像素特征分类法等49。数学形态学运算主要用在图像预处理阶段移掉视乳头中的血管。实验结果表明基于血管跟踪的视乳头定位最可靠;基于活动轮廓模型的视乳头分割最精确;基于圆形Hough变换的视乳头定位与分割的方法虽然简单,但不能满足临床诊断(特别象青光眼的诊断要求;小波变换和完全匹配法直接且简洁,但需要包含了先验知识的精确模板50。其中,基于活动形状模型和分水岭变换是高精度视乳头分割的可选方法。此外,由于临

38、床诊断对视乳头三维可视化和度量的迫切要求,更困难的视网膜(视乳头图像的三维重建和度量34研究工作才刚开始。至今,比较成功的视乳头分割方法是用于青光眼辅助诊断的参数化活动轮廓(蛇形模型;用于糖尿病视网膜病变辅助诊断的参数化活动轮廓模型、活动形状模型和分水岭变换;用于老年性黄斑病变辅助诊断的圆形Hough变换。但是,基于参数化活动轮廓模型的视乳头分割对其初始的活动轮廓非常敏感。R. Chra, stek a, M. Wolf et.al.(2005基于数学形态学运算、Hough变换和猫定的活动轮廓(蛇模型(an anchored active contour model, or snakes,提出

39、了一种用于Heidelberg retina tomograph (HRT共焦激光断层扫描(confocal scanning laser tomograph 视(神经乳头图像视盘定位和分割的方法。James Lowell, AndrewHunter et.al(2004先采用一个精心设计的滤波模板对彩色立体眼底视网膜图像进行滤波来定位视(神经乳头(中心,然后,采用一个全局的参数化椭圆形变形模型和一个局部的、非刚性的、参数化变形模型联合高精度地确定视盘边缘边界线。该两种方法在一定程度上解决了基于经典的参数化活动轮廓模型的视乳头分割对其初始的活动轮廓非常敏感这一问题,但都非常复杂,难以临床应用。

40、另一个经过修改的形变模型是由Juan Xu, Opas Chutatape等提出的,在文献51中,作者采用基于知识的聚类和平滑的方法(Knowledge-based clustering and smoothing来提取有效的杯盘轮廓曲线。作者首先采用一个基于全局参数的形变模型以及Hough变换来得到初始的边缘轮廓,然后通过基于知识的聚类来把边缘分为确定的边缘点和不确定的边缘点,最后通过基于知识的平滑来连接那些确定的边缘点。这比文献48中用到的椭圆匹配的方法来提取正确的背叛曲线更加准确有效。Huiqi Li, Opas Chutatape(2004 先采用主成份分析(PCA定位彩色立体眼底视网

41、膜图像的视(神经乳头(中心,然后,采用一个改进的活动形状模型(ASM确定视盘边缘边界线,其成功依赖于人工标注的训练样本图像集选取和质量,训练样本中未包括的病例很难正确分割。分水岭变换的成功依赖于基于数学形态学的图像滤波器设计好坏,且难以完全避免对比度小和血管遮挡的彩色立体眼底视网膜图像影响。像素特征分类法49则依赖于大量的正确分割样本,而且初次分类得到的结果边缘毛刺和模糊的地方较多,虽然可以通过后续的处理得到清晰地轮廓,但是该方法显得复杂,而且初始的样本的准确度对分割的结果有较大的影响。采用二阶段(园心和半径圆形Hough变换只能获得激光扫描检眼镜(scanning laser ophthal

42、moscope SLO视网膜图像中近似的视乳头边界,在此视乳头形态参数不直接用于老年性黄斑病变的辅助诊断。在文献48中,为了减少水平集方法带来的运算量大的缺点,作者首先用直方图信息来定位视盘,然后来确定视盘轮廓,然后用用水平集方法来得到可能的杯盘形态边缘,最后采用椭圆匹配的方法来去除由于血管遮挡造成的畸变。通过该方法得到的边缘十分光滑,基本反映了杯盘的边缘特征。然而,上述视乳头分割方法大多需要进行非常复杂的初始化处理过程或图像预处理过程。此外,这些方法都未论述对青光眼辅助诊断至关重要的视杯的分割。由刘国才等提出的多层Mumford-Shah向量模型能够很好的进行彩色视乳头图像杯盘重建、分割和度

43、量2224。该方法基于经典Mumford-Shah模型,先用Mumford-Shah模型分割出背景和目标两个部分,这一步可以称为第一层,然后根据需要,在分割出的目标或背景中继续用该模型进行分割,从而得到第二层的背景和目标,如此层层深入,垂直向下,直到分割出待分割的对象为止。该模型还能够很好的处理具有多个目标的多相图像分割问题。二值Mumford-Shah模型及其相应的Chan-Vese模型将待分割图像划分为目标和背景两个部分,但当图像有多个目标(或多相时,该模型无法分出。根据视乳头图像杯盘分割要求,实际上本文需要处理的是一个具有多相目标的图像。由Xue-Cheng Tai等人提出的基于分段常数

44、水平集函数的Mumford-Shad模型(PCLSM也能够很好处理多相目标图像。该模型仅仅只需要一个水平集来分割多相目标。本文拟运用该方法来分割视乳头图像的杯盘,并将在第4、5章重点讨论该方法。1.3本文研究内容本研究以湖南大学王耀南教授参与的重大国际合作欧盟第五框架项目“欧洲青光眼预防工程/计算机辅助青光眼预防诊断项目(EGPS project /Glaucoma Prevention by Computer Aided Diagnostics (GlauCAD No QLG1-2000-00651”为背景,在刘国才副教授负责的国家自然科学基金项目(复杂图像处理中的自由非连续问题及其水平集方

45、法研究,基金号为60872130和两个湖南省自然科学基金项目(基于水平集方法的视网膜图像分割与度量,基金号为06JJ4044,以及生物医学图像处理中的自由非连续问题及其水平集方法研究, 基金号为09JJ3118的支持下,主要研究基于水平集方法的视乳头杯盘分割方法。具体研究内容主要包括基于分段常数水平集方法(PCLSM 的Mumford-Shah模型及其在视乳头图像杯盘分割中的应用;通过分割结果提取对青光眼诊断具有实际意义的病理数据;通过与临床结论以及与多层模型进行对比分析验证以上方法的有效性。本论文主要做如下工作:(1着重研究基于分段常数水平集方法(PCLSM 5254的Mumford-Sha

46、h模型;(2把基于PCLSM的Mumford-Shah模型应用于视乳头图像杯盘分割中。1.4本文的内容安排本文章节安排如下:第一章简要介绍医学图像和视乳头图像杯盘分割的内容以及目前的研究现状;第二章介绍水平集方法的基本原理和方法;第三章介绍Mumford-Shah模型及其简化的Chan-Vese模型,并简要介绍多层Mumford-Shah模型。第四章详细论述基于分段常数水平集方法(PCLSM的Mumford-Shah模型及其数值解法;第五章论述如何将基于PCLSM的Mumford-Shah模型创造性地应用于视乳头图像杯盘分割中, 在得到初步的分割结果之后,再结合凸集约束的理论以及曲线拟合技术,

47、去除由于血管遮挡和光照不均造成的凹陷和断裂,从而得到了完整的杯盘轮廓曲线,在此基础上提取了对青光眼诊断具有实际意义的病理数据,最后将该结论与文献3537中的多层Mumford-Shah模型的分割结果进行对比研究;第六章总结了基于PCLSM的Mumford-Shah模型在视乳头杯盘分割应用中的优点,并提出不足以及今后继续努力的方向。第2章 水平集方法2.1 引言水平集算法主要是从界面传播等研究领域中逐渐发展起来的,它是处理封闭运动界面随时间演化过程中几何拓扑变化的有效的计算工具,1988年,Osher 和Sethian 在文献17中提出了依赖时间的运动界面的水平集描述,并构造了水平集方程的高精度

48、稳定数值解法,从而形成了一种解决曲线演化的数学工具。所谓水平集技术,就是利用水平集的概念和将求解n维描述的演化转化为求解由于具有n 个实变量的水平集函数的演化而导致的水平集的演化这么一种思想方法来达到某种实用目的的技术。水平集技术可用于图像去噪,图像分割,它可以用来解决一些用常规方法难以解决的数学计算问题。水平集方法处理平面曲线演化问题时不是试图去跟踪演化后的曲线位置,而是遵循一定的规律,在二维固定坐标系中不断更新水平集函数,从而达到演化隐含在水平集函数中的闭合曲线的目的,这种方法的最大的特点是:即使隐含在水平集函数中的闭合曲线发生了拓扑结构的变化(合并或分裂,水平集函数仍然保持为一个有效的函

49、数。由于这个特点,近年来,水平集方法成为了研究的热点,并被应用于图像分割、图像平滑、运动分割以及运动目标追踪,甚至立体视觉和图像恢复。2.2 水平集方法概述该算法的基本思想是把N 维描述视为高一维的(N+1水平集,或者说是把N 维描述视为N+1维变量的水平集函数F 的水平集。这样就把求解N 维描述的演化过程转化为求解关于有N+1维变量的水平集函数F 的演化所导致的水平集的演化过程,通过这种转化,就能推演出水平集在某种规律下依各种具体条件而演化的具体演化方程。也即是说,通过这一规律,就有了从一般到特殊的演绎过程的出发点和依据。先给定一个初始轮廓(0,该初始轮廓不断延法线方向向外或向内以一定的速度

50、演化,则可以得到一个曲线族(t 。而水平集方法把曲线族看成是高一维空间曲面,(t y x 的零水平集0=,其中,(t y x 就称作水平集函数,当轮廓改变拓扑结构时,(t y x 仍然保持连续,这样就从原理上保证了这种方法很容易处理拓扑结构改变的问题。如下图2.1所示,零水平集延z轴演化时分裂和合并,但是,(t y x 仍然保持为一个整体。 图2.1 水平集的曲线演化在经典的水平集方法中,通常采用符号(欧氏距离函数D (x,t ,即从x 到(t 的欧氏距离作为水平集函数, 且当x 位于内部时,取负值;当x 位于外部时,取正值;因此,通常水平集函数定义如下式:=_,(;,0;,(,(x t x

51、D x x t x D t x (2.1 一般情况下,可以采用距离重构方法(Reinitialization 将非符号距离水平集函数水平集函数转化为符号距离水平集函数,以保证水平集函数在(t 附近具有较好几何特性:既不太陡峻,也不太平坦,这样有利于数值处理。距离重构方法之一是求解下列方程的稳态解:,(0,(01|(|sgn(t x x t=+, (2.2 ,1(2sgn(=x H x 式中, H (x为一维的Heaviside 函数。在经典的水平集方法中,要求v N 在(t 的附近具有好的定义,因此常需要将定义于(t 的速度场v N 进行光滑的扩展,以便研究(t 的演化特性。光滑扩展速度场v

52、N 方法之一是求解下列方程的稳态解:.1(2sgn(,(0,(0|(sgn(=+x H x t x v x tN 式中, (2.3水平集处理的基本方程有:考虑零水平集x(t所对应的水平集函数,则有0,(=t t x (2.4 对上式两边对时间t 求偏导,有0,(=+t x t t x t (2.5 假设F 为外法向方向的速度,那么F n tx = 其中|/=n 。这样得到基本方程: 0|=+F t(2.6 除此基本方程式以外,还有其它一些方程式,它们都能从式2.4推导出来。 基本方程式2.6是水平集函数及相应的水平集在法向力F 的推动下的演化方程。 用水平集方法进行图像处理可形成下面的一般性算

53、法:(1设定水平集函数的初态;(2确定动力F 的形式;(3按基本方程推演水平集函数的各状态;(4对于每一水平集函数的状态求解零水平集。尽管具体技术细节千变万化,水平集图像处理的算法万变不离其宗,都是围绕上述一般性算法而展开的。由上面的叙述可以可知,水平集方法实现曲线演化具有如下优点:(1只要速度场v N 光滑,则水平集函数,(t y x 始终保持为一个连续函数,那么其水平集,即演化曲线可以随着,(t y x 的演化很自然的改变其拓扑结构,可以分裂、合并、形成尖角等,如图2.1所示,曲线可以分裂为多个闭合曲线,但,(t y x 仍然保持为一个函数;(2因为,(t y x 演化时始终保持为一个函数

54、,所以很容易实现数值近似算法,可以利用离散网络结构的有限差分方法实现水平集函数,(t y x 的演化;(3演化曲线的内在特征,例如曲率等,可以由,(t y x 直接计算;(4水平集方法很容易扩展到高维,比如可以扩展到三维闭合曲面的演化,可以用于三维图像的分割。2.3 水平集方法的数值实现由上面可以知道,由于,(t y x 演化时始终保持为一个函数,因此可以用离散网格来表达水平集函数,(t y x ,在文献17中详细的说明其数值化方法,包括水平集函数的网格表达,一阶和二阶微分的计算,内向单位法矢以及曲率的计算等等,先简要介绍如下:假设离散网格的间距为h ,时间步长为t ,则在第n 时刻,网格点i

55、j 处的水平集为n ij ,且,(t n jh ih n ij =,则可将2.6式离散化为:01=+n ij ij n ij nij n ij F t (2.7式中,n ij F 表示外法向方向的速度在ij 处n 时刻的值,另外需要注意的是,(t y x 的有限数值差分n ij ij 必须采取适当的形式,正如文献17中说明那样,用以避免常量演化中导致的水平集函数奇异性问题。对于式2.7,它属于Hamilton-Jacobi 方程,可以写成如下形式:0,(=+y x t H (2.8 其中的(H 为哈米顿算子,它需要采用合适的空间导函数值近似方法才能稳定地获得水平集方程的唯一粘性解。涉及水平集函

56、数梯度或曲面法向量数值计算的部分,如水平集方程中的对流运动部分,需要采用迎风差分近似空间导函数,或者引入人工扩散(artificial diffusion 。如果偏微分方程(PDE 中包含足够大的扩散项,不管该扩散来自二阶空间导数,还是来自人工扩散(如Lax-Friedrichs 扩散,则采用迎风差分或者中心差分近似空间导数都能保证数值算法的稳定性,特别如水平集方程中依赖于平均曲率的运动部分,其中涉及的一阶和二阶空间导数都可以采用中心差分近似。式2.8的主要特征是方程具有连续的解,但是导数却不连续。即使对于光滑的初值0,其解也可能产生奇异点,在处理这种问题的时候必须基于波传导理论中的“熵守恒”

57、理论,将该理论运用于偏微分方程的求解上,就是所谓的偏微分方程的迟滞解方法。另外,根据Sethian 在文献55,56中所说的,速度F 必须包含常量速度0F 和基于曲率的速度(K F ,即:(0K F F F +=。其中0F 为轮廓扩张或收缩的基本速度,为速度常量(0.10±=F ,它用来控制曲线按照特定的标记进行均一的扩张或收缩,为正表示向外扩展,为负表示向内收缩;(K F 则与曲率有关。为了是曲线在靠近目标物体的边界时速度逐渐减慢,并最终停止于物体边界,速度F 还应包含速度函数停止因子,(y x C ,该因子包含图像的边界信息,当曲线靠近物体边界的时候,(y x C 使F 趋近于零,从而使曲线能够停止在物体边界上,所以速度函数F 可表示为下式:(,(0K F F y x C F += (2.9将2.9代入2.6可得:+=(,(0K F F y x C t(2.10 对于式2.7中的n ij ij n ij F 计算得到:2/12,2,2,2,0min(0,max(