CN109272512B - 一种自动分割左心室内外膜的方法 - Google Patents

一种自动分割左心室内外膜的方法 Download PDF

Info

Publication number
CN109272512B
CN109272512B CN201811116609.XA CN201811116609A CN109272512B CN 109272512 B CN109272512 B CN 109272512B CN 201811116609 A CN201811116609 A CN 201811116609A CN 109272512 B CN109272512 B CN 109272512B
Authority
CN
China
Prior art keywords
model
segmentation
image
level set
left ventricle
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811116609.XA
Other languages
English (en)
Other versions
CN109272512A (zh
Inventor
李军华
李林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanchang Hangkong University
Original Assignee
Nanchang Hangkong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanchang Hangkong University filed Critical Nanchang Hangkong University
Priority to CN201811116609.XA priority Critical patent/CN109272512B/zh
Publication of CN109272512A publication Critical patent/CN109272512A/zh
Application granted granted Critical
Publication of CN109272512B publication Critical patent/CN109272512B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00Pattern recognition
    • G06F18/20Analysing
    • G06F18/23Clustering techniques
    • G06F18/232Non-hierarchical techniques
    • G06F18/2321Non-hierarchical techniques using statistics or function optimisation, e.g. modelling of probability density functions
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/12Edge-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/13Edge detection
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10088Magnetic resonance imaging [MRI]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20048Transform domain processing
    • G06T2207/20061Hough transform
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/30Subject of image; Context of image processing
    • G06T2207/30004Biomedical image processing
    • G06T2207/30048Heart; Cardiac

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Probability & Statistics with Applications (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种自动分割左心室内外膜的方法,其步骤如下:1)选取需要处理图像;2)选用优化Mean Shift聚类算法对左心室图像进行预处理;3)用改进hough变换圆检测算法,对左心脏内腔定位,得到内外膜分割初始轮廓;4)用双水平模型分割图像内外膜;5)双水平集模型能量函数公式收敛性检查;6)得到分割内外膜轮廓分割效果图。本发明解决了人为设置分割初始位置问题,提高了内外膜的分割精度,降低了模型分割内外膜的时间,得到无边缘泄露且正则性良好的符合临床定义的内外膜轮廓。

Description

一种自动分割左心室内外膜的方法
技术领域
本发明涉及医学图像处理技术,特别涉及一种自动分割左心室内外膜的方法。
背景技术
近年来,心脏疾病死亡人数的不断增加,引发人们对心脏的更深探索。对左心室外膜准确的提取是进一步研究心脏左心室的重要基础。在临床应用方面,左心室外膜的准确分割能提供一些重要的左心室功能参数,利于医生定量分析左心室功能是否异常;在科学研究方面,左心室的分割是左心室三维建模的基础。一个心脏振动周期的图像数据量大,想在短时间内得到客观、精准的内外膜分割结果,就需要依赖优秀的算法和模型来对左心室进行自动的分割。
当前在医学领域的计算机视觉技术,需解决:
1.对DICOM格式的医学图像(当下医学影像存储、交换及传输的国际标准:CT、核磁共振图像和超声图像都是以该格式进行存储)的调用、去噪、平滑图像处理技术。
2.使用优化Hough变换圆检测模型,对左心室进行定位,利用matlab工具箱函数提取检测圆边缘。
3.利用新提出双阱势函数和各向异性梯度矢量流构建双水平集模型,由Hough变换定位内外膜初始位置,由双水平集模型完成左心室内外膜的分割。
左心室内外膜的分割技术,其方法是利用图像基础信息进行处理,提取出符合临床定义的做心室内外膜。心室内外膜分割是心脏三维模型构建的基础,是心脏四维化显示的关键;提取内外膜的精度越高,对应的心脏功能参数就越准确,就能进一步辅助医生治疗。吸引越来越多国内外优秀学者投入到左心室内外膜分割的研究中来。近些年出现了多种左心室内外膜分割方法,可以大致将他们分为两类:基于主动轮廓模型分割法和基于水平集模型分割法。
主动轮廓模型分割法主要包括基于边缘的主动轮廓模型和基于区域的主动轮廓模型。王元全等利用预测-校正改进的基于边缘主动轮廓模型,用外膜轮廓先验知识对模型预测结果进行校正分割[1];段先华等提出融合滤波的基于区域的主动轮廓模型分割,提高了分割精度[2];Lee H Y等利用区域生长和迭代阈值法分割做心室内膜,用提取的内膜边界引导主动轮廓模型提取外膜[3];朱锴等用结合K均值的自适应核函数带宽的Mean Shift算法分割图像,区域生长法提取目标区域[4]
上述方法都只用到了图像的局部信息,没有考虑到图像全局信息的影响,而且半自动手动校正增加了人为主观因素影响,降低了模型方法分割效率。目前心血管疾病的预防和治愈率能力较弱,迫切的需要一种方法提高左心室内外膜分割精度,辅助医生判断,提高心血管疾病预防能力和治愈率。常用的水平集模型是无需反复初始化的经典水平集模型。王兴家等利用改进的耦合水平集分割算法0,利用心室图像间的结构相似性,定位模型演化的初始轮廓,完成心室膜分割。
水平集算法是一种常用于医疗图像分割的变形模型算法,该模型有较好的拓扑变化能力、不需要精确的初始位置和有稳定的数值近似方法。但当分割图像边缘模糊,提取轮廓边缘正则性差,且分割精度较低;小梁肌和***肌和左心室心肌相连,分布于左心室腔室内,利用模型直接分割,会使分割结果不符合临床上定义小梁肌和***肌归属于心肌的要求。
上述方法中,对于主动轮廓模型分割左心室内外膜有以下几点不足:
①对初始化位置比较敏感,曲线拓扑结构在演化过程不会发生改变。
②基于边缘和区域的模型分割,由于能量函数的非凸性,使曲线在演化过程中容易陷入局部极小值,导致分割失败。
③预测矫正分割方法,手动干预较多,增加算法复杂度和人为主观因素影像,导致分割误差偏大。
而对于水平集模型分割,目前虽然有所改进,但还存在一些问题:
①改进的耦合水平集模型,利用图像之间结构相似性,定位获取左心室腔初始位置,分割出左心室有较好的完整性,但分割出心肌平均相似度较低仅有90%。
②基于形状约束的水平集模型,初始位置需要手动定位,分割效果为满足腔室类圆形状,但导致分割结果错误划分小梁肌和***肌位置。
③模型分别分割内外膜,效率较低。
发明内容
针对现有技术中存在的问题,本发明提供一种自动分割左心室内外膜的方法,以解决DRLSE水平集模型在模糊边界分割的低精确度,边缘正则性差和分割不符合临床要求的问题。
本发明的目的是这样实现的。一种自动分割左心室内外膜的方法,其步骤如下:
1)选取需要处理的左心室MRI图像,并输入图像;
2)图像聚类处理,选用优化Mean Shift聚类算法对左心室MRI图像进行预处理;由Mean Shift向量的基本形式,在d维空间中n个样本点xi,i=1,2,…,n,则:
Figure GDA0003461913530000021
其中:x表示空间中任选的一点,k表示在这n个样本点xi中,有k个点落入sk区域中;
上式中,Sh是一个半径h的球形区域,定义为:
Sh=(y|(y-x)(y-x)T≤h2) (2)
其中:y是在一个半径为h的高维球区域sh,满足公式(2)的集合,x表示空间中样本点;引入核函数和增加样本权重来解决(1)式中存在问题;在核函数中引入高斯核函数:
Figure GDA0003461913530000031
得到改进后的Mean Shift聚类算法形式:
Figure GDA0003461913530000032
其中,G(x)是单位核函数:
Figure GDA0003461913530000033
其中,H是带宽矩阵,即正定的n×n矩阵;(4)式中ω(xi)≥0是添加的每个样本的权重;带宽矩阵H为:
Figure GDA0003461913530000034
其中:h1,h2,hn表示带宽参数;
于是推出Mean Shift概率密度梯度归一化函数:
Figure GDA0003461913530000035
概率密度函数f(x)在d维空间内n个采样点xi,i=1,2,…,n,则f(x)的核密度为:
Figure GDA0003461913530000036
f(x)的梯度
Figure GDA0003461913530000037
的核密度为:
Figure GDA0003461913530000041
定义g(x)=-k′(x),G(x)=g(||x||2),推出:
Figure GDA0003461913530000042
由上式可得,Mean Shift向量就是右边第一个括号中与概率密度梯度函数成正比;则:
Figure GDA0003461913530000043
化简上式,令
Figure GDA0003461913530000044
当且仅当Mh(x)=0时,
Figure GDA0003461913530000045
得到新的圆心位置:
Figure GDA0003461913530000046
定义在以x为圆心,h为半径,总的核密度为:
Figure GDA0003461913530000047
其中,sh,hτ分别是半径为h、τ的高维球区域,C为常数,x、xi为样本点;
Figure GDA0003461913530000048
表示颜色信息,
Figure GDA0003461913530000051
表示空间位置信息;
3)内外膜定位,用改进hough变换圆检测算法,对左心室MRI图像内腔定位,得到左心室MRI图像内外膜分割初始轮廓;
优化Hough变换圆检测算法,在对图像进行找圆定位时,对输入图像进行阈值处理,像素较大点连成区域设为前景区域,累加器先累加前景区域的像素点;利用图像边缘信息进行约束,沿着梯度方向的像素点首先进行投票;
4)用双水平集模型分割左心室MRI内外膜;
用基于ODRLSE模型的0水平集和k水平集融合的双水平集模型,同时分割左心室MRI图像内外膜;双水平集模型的能量函数公式为:
Figure GDA0003461913530000052
其中:μ、λ、α均是大于零的常数;
Figure GDA0003461913530000053
模型距离正则项,
Figure GDA0003461913530000054
是模型边缘约束项,
Figure GDA0003461913530000055
是模型演化速度控制项;
Figure GDA0003461913530000056
保持双水平集模型演化过程中曲线的正则性;
Figure GDA0003461913530000057
融合AGVF约束0水平集和k水平集朝着内外膜边缘及边缘凹陷处演化;
Figure GDA0003461913530000058
约束模型曲线朝目标边缘速率;
定义
Figure GDA0003461913530000059
如下:
Figure GDA00034619135300000510
Figure GDA00034619135300000511
Figure GDA00034619135300000512
其中:p为势函数,p:[0,∞)→R,
Figure GDA00034619135300000515
表示模型轮廓,λ1、λ2、α1、α2均是大于零的常数;且δ和H分别为狄拉克函数和赫维赛德函数:
Figure GDA00034619135300000513
Figure GDA00034619135300000514
其中:ε设置的常数;
边缘指定符公式为:
Figure GDA0003461913530000061
式中:Gσ为标准偏差σ的高斯核函数,I的定义在整个图像区域;
双水平集模型将水平集演化导出为梯度流,用优化距离正则化项和外部能量来最小化能量泛函,从而驱动水平集向内外膜边缘位置移动;
5)对双水平集模型能量函数公式进行收敛性检查,若达到能量最小,分割结束;否则,重复步骤4);
6)得到左心室MRI图像分割内外膜轮廓分割效果图。
本发明解决了人为设置分割初始位置问题,提高了内外膜的分割精度,降低了模型分割内外膜的时间,得到无边缘泄露且正则性良好的符合临床定义的内外膜轮廓。本发明用Mean Shift聚类算法对图像进行预处理,增大了内外膜之间像素的对比度;Hough变换圆检测自动对初始位置进行定位,解决了之前手动定位的繁琐流程;检测初始位置距离内膜较近,减少了曲线演化时间,提高了分割效率;以分割内膜作为外膜分割的初始轮廓,正确处理了小梁肌和***肌的位置归属,使分割结果符合临床定义。
附图说明
图1是左心室内外膜分割流程图;
图2是本发明中待处理的左心室MRI图像;
图3是本发明中Mean Shift聚类算法处理图像的效果图;
图4a是Hough变换圆检测的边缘检测结果图;
图4b是Hough变换圆检测的自适应定位图;
图4c是Hough变换圆检测的输入圆心后的定位图;
图5a是本发明中左心室MRI图像经OTSU算法处理得到的二值图;
图5b是本法明中改进Hough变换圆检测方法中自适应检测的内膜图;
图6a是本发明中双阱势函数对比图;
图6b是本发明中双阱势函数对应扩散速率图;
图7a是图像原始外力场分布图;
图7b是本发明中AGVF外力场的分布图;
图8a是分割原始图;
图8b是图8a的聚类效果图;
图8c是图8b的阈值处理效果图;
图8d是图8c的初始定位图;
图8e是图8a的分割效果图;
图8f是图8e分割效果的三维展示图;
图9a是手动分割效果图;
图9b是Snake模型半自动分割效果图;
图9c是DRLSE模型分割效果图;
图9d是本发明中左心室内外膜分割效果图;
图10a是手动分割效果图;
图10b是Snake模型分割结果图;
图10c是DRLSE模型分割图;
图10d是RSF模型分割图;
图10e是CV模型分割图;
图10f是本发明中左心室内外膜的分割效果图;
图10g是本发明中内外膜分割的三维能量图。
具体实施方式
以下结合附图对本发明作进一步说明。参见图1至图10g,一种自动分割左心室内外膜的方法,其具体步骤如下:
1)选取需要处理的左心室MRI图像(如图2),并输入图像101;
2)图像聚类处理102,选用优化Mean Shift聚类算法对左心室MRI图像进行预处理:
解决模糊边界:
由于图像左心室边界与右心室距离较近,图像底层视觉特征对比度小,为更好得到检测初始位置,图像先用模糊集合处理。由于Mean Shift算法主要应用于计算机视觉和图像处理中的聚类分析,是一种用于定位密度函数最大值的非参数特征空间分析技术。本发明选用优化Mean Shift算法对图像进行预处理。
一般Mean Shift聚类原理:
Mean Shift向量的基本形式,在d维空间中n个样本点xi,i=1,2,…,n,则:
Figure GDA0003461913530000071
上式中,Sh是一个半径h的球形区域,y是球形区域内点,定义为:
Sh=(y|(y-x)(y-x)T≤h2) (2)
这种基本形式存在不足,每点对x的作用是由x的位置决定的,而(1)式中在球形区域的各点对x的作用毫无差异性。
优化Mean Shift聚类算法:
本发明决定在引入核函数和增加样本权重来解决(1)式中存在问题。核函数是一个分段连续、非负、非增函数,引入高斯核函数:
Figure GDA0003461913530000081
得到改进后的Mean Shift聚类算法形式:
Figure GDA0003461913530000082
其中,G(x)是单位核函数:
Figure GDA0003461913530000083
其中,H是带宽矩阵,即正定的n×n矩阵。(4)式中ω(xi)≥0是添加的每个样本的权重。带宽矩阵H为:
Figure GDA0003461913530000084
于是推出Mean Shift概率密度梯度归一化函数:
Figure GDA0003461913530000085
而Mean Shift聚类算法是利用概率密度求得其的局部最优解。对于一个概率密度函数f(x)在d维空间内n个采样点xi,i=1,2,…,n,则f(x)的核密度估计为:
Figure GDA0003461913530000091
f(x)的梯度
Figure GDA0003461913530000098
的核密度估计为:
Figure GDA0003461913530000092
定义g(x)=-k′(x),G(x)=g(||x||2),推出:
Figure GDA0003461913530000093
由上式可得,Mean Shift向量就是右边第一个括号中与概率密度梯度函数成正比。则:
Figure GDA0003461913530000094
化简上式,令
Figure GDA0003461913530000095
当且仅当Mh(x)=0时,
Figure GDA0003461913530000096
得到新的圆心位置:
Figure GDA0003461913530000097
定义在以x为圆心,h为半径,总的核密度为:
Figure GDA0003461913530000101
其中,
Figure GDA0003461913530000102
表示颜色信息,
Figure GDA0003461913530000103
表示空间位置信息;MeanShift聚类算法处理图像效果图(见图3)。
3)用改进hough变换圆检测算法,对左心脏内腔定位,得到内外膜分割初始轮廓103;避免手动初始定位:
传统水平集模型对初始位置比较敏感,为克服这一问题,用优化Hough变换检测圆自动获取初始轮廓。
传统Hough变换圆检测算法原理是利用点与线的的对偶性,将图像空间的线条变为空间的聚焦点,从而检测给定图像是否存在给定性质的曲线。圆的方程为:
(x-a)2+(y-b)2=r2 (16)
其中:(a,b)为圆心,(x,y)是待检测区域坐标,r为半径。首先用算子对图像图2进行边缘检测得到二值图像图4a;Hough变换会将图像空间对应到(a,b,r)参数空间,取和图像平面一样的参数平面,以二值图像中每一个非零点为圆心,用设定的半径在参数平面上画圆,符合(a,b,r)的点进行累加,即对圆覆盖的坐标点进行投票,随着半径步长一直循环运算到设定半径的最大值;最后找出参数平面上的峰值,得到(a,b,r)值及其在图像对应位置显示形状,而后通过调整阈值得到目标圆形。hough变换圆检测算法自适应定位结果见图4b,输入圆心位置定位结果见图4c;
为提高Hough变换圆检测算法自适应定位精确度,解决手动设置圆心位置定位的复杂流程,加快运算效率。优化Hough变换圆检测算法,在对图像进行找圆定位时,对输入图像进行阈值处理(参见图5a),像素较大点连成区域设为前景区域,累加器先累加前景区域的像素点;为避免检测出多个无用圆,利用图像边缘信息进行约束,在有限时间间隔内只允许沿着梯度方向的像素点进行投票。这样无需对预处理图像在进行遍历预算,大大提高了运算效率和定位精度(参见图5b)。原Hough变换圆检测耗时3.209秒,改进后只需要0.392秒就能准确到心室内膜。
4)用双水平模型104分割图像内外膜:
水平集模型:
水平集的方法思想就是用高维度表示低维度,在高维曲面用隐函数等值点来表示平面上闭合曲线轮廓。在传统的水平集方法模型中,出现反复初始化、曲线结果正则性差等问题,Li C提出了一种经典的的距离正则性水平集模型方法,解决了水平集模型分割时初始化问题。根据曲线的运动规律,定义能量泛函为:
ε(φ)=μRp(φ)+εext(φ) (17)
其中:Rp(φ)是定义的距离正则项,μ是大于零的常数,外部能量项εext(φ)由图像先验信息或图像底层视觉特征决定。水平集的距离正则项被定义为:
Figure GDA0003461913530000111
其中:p为势函数,p:[0,∞)→R,
Figure GDA0003461913530000112
是梯度矢量流。传统势函数定义为:
Figure GDA0003461913530000113
偏微分方程:
Figure GDA0003461913530000114
扩散率:
Figure GDA0003461913530000115
提出新双阱势函数:
经典水平集分割模糊边界得到的边界正则性太差,为得到更加平滑的分割轮廓曲线,提高左右心室间模糊边界的分割精度和分割效率,在经典水平集模型中提出新的双阱势函数(NDP,new doublewell potential)代替传统势函数(DP,doublewell potential)和约束模型曲线演化。
提出的双阱势函数:
Figure GDA0003461913530000116
如此定义,知双阱势函数有两个最小值分别是
Figure GDA0003461913530000121
Figure GDA0003461913530000122
扩散率:
Figure GDA0003461913530000123
偏微分方程:
Figure GDA0003461913530000124
且:
Figure GDA0003461913530000125
证明了势函数的扩散速率是有界的。得到一个双阱势的距离正则化项,水平集的演化是一个梯度流,它可以使能量函数最小化;且对水平集函数的梯度模进行约束,避免传统水平集曲线演化中的重新初始化。
由图6a和图6b得:
(1)当
Figure GDA0003461913530000126
r2>r1>0,曲线正向扩散,越扩散
Figure GDA0003461913530000127
越小,曲线运动到
Figure GDA0003461913530000128
NDP需要的迭代次数明显少于DP模型方法;
(2)当
Figure GDA0003461913530000129
r1<r2<0,曲线反向扩散,越扩散
Figure GDA00034619135300001210
越大,NDP模型方法分割效率高于DP分割方法;
(3)当
Figure GDA00034619135300001211
曲线正向扩散,越扩散
Figure GDA00034619135300001212
越小,曲线演化的的结果是
Figure GDA00034619135300001213
Figure GDA00034619135300001214
区域内距模型方法等值线较近,模型扩散速率对曲线演化影响较小;
总之,曲线扩散向着使
Figure GDA00034619135300001215
Figure GDA00034619135300001216
的方向进行。
构造AGVF模型:
由于DRLSE模型在分割凹陷区域时,容易出现欠分割现象,很难正确分割出目标轮廓。因梯度矢量流能将模型分割曲线演化引向边缘深凹处,为使模型有较好分割效果,引入AGVF用以改进模型的外力能量函数,定义图像梯度矢量流为V(x,y)=(u(x,y),v(x,y))。
先对扩散矩阵进行构造,定义结构四通道张量:
Figure GDA00034619135300001217
式中:J是通道张量,
Figure GDA00034619135300001218
是对u求梯度,ux,uy是u对图像x,y方向的一阶偏导,对各个通道进行非线性处理,
Figure GDA0003461913530000131
Figure GDA0003461913530000132
分解非线性结构张量,得到法线和切线方向对应特征值λ12及相应特征向量μ1,μ2。根据特征值对其进行重新定义,即:
Figure GDA0003461913530000133
式中:C1,C2是常数。定义扩散矩阵:
Figure GDA0003461913530000134
根据图像信息构建AGVF模型,图像梯度矢量流为V(x,y)=(u(x,y),v(x,y)),Ui分割区域内第i类信息矩阵,g(x)和h(x)分别是约束分割区域和边缘的函数,
Figure GDA0003461913530000135
K梯度敏感参数,
Figure GDA0003461913530000136
融合图像的各向异性。
Figure GDA0003461913530000137
χ,μ是控制参数。由
Figure GDA0003461913530000138
根据上述扩散矩阵构造过程推出D。令U=Ui,当AGVF能量函数取极小值时,推导出新的扩散模型E,且欧拉方程组为:
Figure GDA0003461913530000139
图7a图像原始外立场凹陷处呈水平状态,图7b AGVF外立场指向凹陷区域深处。对比发现引入AGVF处理图像,扩大了模型的扑捉范围。
双水平集模型:
由于单一水平集曲线演化模型分割心室内外膜边缘,需要进行两次分割,比较耗时,且心室内外膜边缘像素点对比度不同,同一模型分割会造成边界内外膜分割结果出现泄漏,为提高模型分割效率和分割精度,提出用基于ODRLSE模型的0水平集和k水平集融合的双水平集模型,同时分割心室内外膜。
双水平集方法的能量函数:
Figure GDA0003461913530000141
其中:
Figure GDA0003461913530000142
模型距离正则项,
Figure GDA0003461913530000143
是模型边缘约束项,
Figure GDA0003461913530000144
是模型演化速度控制项。
Figure GDA0003461913530000145
保持双水平集模型演化过程中曲线的正则性;
Figure GDA0003461913530000146
融合AGVF约束0水平集和k水平集朝着内外膜边缘及边缘凹陷处演化;
Figure GDA0003461913530000147
约束模型曲线朝目标边缘速率。
定义
Figure GDA0003461913530000148
如下:
Figure GDA0003461913530000149
Figure GDA00034619135300001410
Figure GDA00034619135300001411
且δ和H分别为狄拉克函数和赫维赛德函数:
Figure GDA00034619135300001412
Figure GDA00034619135300001413
边缘指定符:
Figure GDA00034619135300001414
式中:Gσ为标准偏差σ的高斯核函数,I的定义在整个图像区域。
5)双水平集模型能量函数公式收敛性检查,若达到能量最小,分割结束;否则,重复步骤4)105;
6)得到分割内外膜轮廓分割效果图106。
双水平集模型将水平集演化导出为梯度流,用优化距离正则化项和外部能量来最小化能量泛函,从而驱动零水平集向期望位置的运动。距离正则项约束拓扑水平集演化具有唯一的向后和向前扩散效应,它能够维持水平集函数的期望形状,尤其是在水平集附近有符号距离分布。距离正则化消除了传统水平集中需要的重新初始化,简化了运算。选用有限差分法来量化图像梯度,外部能量函数引入AGVF,驱动模型在目标轮廓凹陷也能进行完全分割,减少迭代次数的同时保证应有的计算精度。
本发明首先采用模糊集合对图像进行处理得到特征图像;其次,根据左心室医学图像特征,采用Hough变换圆检测算法对左心室内外膜进行粗略定位,检测结果作为优化水平集模型的初始轮廓,由模型完成左心室内膜分割;最后,利用提取内膜边缘完成对心室外膜的分割。
双阱势函数和各向异性梯度矢量流引入DRLSE模型,融合改进0水平集和K水平集的双水平集模型应用于左心室内外膜分割。观察分割左心室CT图像过程,分割时间大幅度降低,低至9秒内;分割结果边缘无泄漏且正则性良好;分割MRI图像精确度提高明显,分割重叠率D ice提高至0.9529。总体来说,无论是CT图像,还是MRI图像,利用双水平集模型分割左心室内外膜,都提高了边缘分割精度,减少了分割时间,得到符合临床定义的正则化良好的内外膜轮廓。
水平集模型自动分割左心室内外膜步骤如下:
Step 1:选择图像如图8a;
Step 2:对图像进行Mean Shift聚类处理(见图8b);
Step 3:对图像进行OTSU算法自适应处理(见图8c);
Step 4:利用改进Hough变换圆检测,对图像心室内外膜进行定位(见图8d);
Step 5:双水平集模型完成对内外膜的分割,其中模型各参数值是实验所得;
Step 6:得到左心室内外膜分割结果(见图8e),及分割三维效果图(见图8f);
改进算法的优点:融合了Mean Shift聚类算法和Hough圆检测得优化水平集模型算法。增大了模糊边界像素间的对比度;Hough变换圆检测完成了对内膜的初定位和给定模型的初始演化位置;提出新的双阱势函数提高了水平集模型的分割的边缘正则性;各向异性梯度矢量流使得水平集模型在凹凸边界也能得到精准分割;双水平对内外面同时分割,提高了模型分割内外膜的效率。
总结:①Hough变换圆检测对分割图像进行定位,解决了以往模型分割图像初始位置手动输入的问题,减少了人为的主观性影响;
②模型分割初始位置距离提取边缘较近,减少了模型曲线迭代次数,提高了模型分割效率;
③各向异性梯度矢量流的加入,促使曲线能对凹陷区域深入分割,提高了模型的分割精度。
④双水平集同时分割内外膜,提高了分割效率。
实验结果:本发明实验数据集采用了中国台湾中原大学CT心脏影像库和法国鲁昂大学计算机工程与科学MRI心脏影像库,中原大学数据库包含十份心脏的不同振动周期,每个心脏振动周期有387张图像,图像大小256×256;鲁昂大学数据库包含48人心脏序列,每份包含279张图像,图像大小为256×216。程序处理过程中,本发明采用模糊集合、Hough变换和ODRLSE模型方法,试验参数由多次实验所得,设置ε=1.5,μ=0.2,λ1=5,α1=-3.0,λ2=6,α2=-4.0。实验配置处理器为Intel(R)Xeon(R)@1.60GHz 1.60GHz,RAM8GB开发软件平台为Windows7,matlab2014a。
本发明分别用CT图像和MRI图像展示分割效果,其中图9a-图9d是案例一CT图像,图10a-图10f是案例二MRI图像。图9d和图10f是本发明方法的分割效果图。为了更好的证明本发明方法的精准、高效性,做了多组对比试验,即与经典RSF模型、CV模型、基于区域的主动轮廓模型(Snake)、传统DRLSE模型分割左心室外膜结果进行对比。用图9a和图10a手动分割结果作为标准。图9b和图10b是Snake模型分割结果;图9c和图10c是DRLSE模型分割结果;图10d是RSF模型分割结果;图10e是CV模型分割结果。
为更直观看出分割效果,引入三维能量图,定义初始轮廓处的符号距离函数值为零,轮廓内部的符号距离函数定义为正,外部则小于零;中间蓝色线为模型曲线演化最终边缘,图10g是本发明方法分割心脏左心室的三维能量图,能更清晰的看出分割出心室外膜轮廓的效果。
案例一用传统水平集和半自动分割方法与本发明方法同时分割CT图像做对比试验。水平集方法分割模糊边缘时,出现欠分割和过分割现象;基于区域的半自动Snake模型分割方法,结果相对精确,但分割流程相对复杂,人为主观性影响较大,如图9b较粗线为手动校正二次分割结果;本发明方法解决心室膜模糊边界和凹陷区域的分割缺陷,分割出符合左心室外膜临床定义的类圆形边缘。
结合表1案例一数据,本发明提出方法分割外膜迭代次数明显少于DRLSE模型分割,且时间降低至平均8.4833秒,与人机交互二次分割相比,简化了分割流程,大大提高了分割效率;Snake模型方法分割结果Dice平均为0.8934,DRLSE模型方法分割结果Dice平均为0.9081;本发明方法分割结果Dice平均为0.9292,均高于Snake模型和DRLSE模型的分割结果,本发明方法提高了心室外膜分割的精确度。
表1案例一各方法模型分割结果对比
Figure GDA0003461913530000171
其次,为验证本发明方法具有良好的分割效果,分割MRI图像做对比试验时,加入了基于边缘分割的RSF模型方法和基于区域的CV模型方法作对比,由图10b-图10e直接从视觉角度就能看出RSF、CV、Snake、DRLSE模型方法分割左心室外膜时,分割所得边缘正则性差,在模糊边缘处出现过分割现象。对比图10f和表2数据,本发明方法分割的左心室外膜边缘正则性好,精确度高,所得边缘符合临床定义的类圆形。
表2案例二各模型分割结果对比
方法模型 Dice
Snake模型 0.9244
DRLSE模型 0.9438
RSF模型 0.9298
CV模型 0.9125
本发明方法 0.9529
通过上述分析,本发明方法降低分割左心室外膜的耗时,提高了心室外膜模糊边缘和凹陷区域的分割精度,得到边缘正则性良好的轮廓。证实本发明方法无论分割CT心室图像还是MRI心室图像,都能有效的完成左心室外膜的自动分割。
本发明提出基于双阱势模型的左心室内外膜的分割方法,利用优化双阱势函数提高模型的演化速度,引入AGVF的外部能量函数驱动模型曲线引向深凹处边缘,完成左心室外膜的分割。通过多种方法分割效果对比,本发明方法无论分割CT图像还是MRI图像,都提高了外膜模糊边缘和凹陷边界的分割精度和分割效率,并自动分割出符合临床定义的椭圆形心室外膜,具有一定的应用价值。
参考文献:
[1]王元全,汤敏,夏德深等著《基于预测-校正改进Snake的心脏MR图像左心室外轮廓分割》[J],——《模式识别与人工智能》2004,17(1):000047-51。
[2]段先华,吴小俊,夏德深著《带标记线的心脏MRI图像分割方法》[J],——《计算机应用研究》,2007,24(7):295-297。
[3]Lee H Y,Codella N C,Cham M D,et al.《Automatic left ventriclesegmentation using iterative thresholding and an active contour model withadaptation on short-axis cardiac MRI》[J].IEEE《transactions on bio-medicalengineering》,2010,57(4):905-13。
[4]朱锴,付忠良,朱硕等著《基于自适应均值漂移的超声心动图左心室分割方法》[J].——《生物医学工程学杂志》,2018(2)。
王兴家,董利娜,李传富等著《用改进的耦合水平集方法从MSCT中分割左心室》[J].——《中国生物医学工程学报》,2011,30(4):494-499。

Claims (1)

1.一种自动分割左心室内外膜的方法,其步骤如下:
1)选取需要处理的左心室MRI图像,并输入图像;
2)图像聚类处理,选用优化Mean Shift聚类算法对左心室MRI图像进行预处理;
由Mean Shift向量的基本形式,在d维空间中n个样本点xi,i=1,2,…,n,则:
Figure FDA0003317183700000011
其中:x表示空间中任选的一点,k表示在这n个样本点xi中,有k个点落入sk区域中;
上式中,Sh是一个半径h的球形区域,定义为:
Sh=(y|(y-x)(y-x)T≤h2) (2)
其中:y是在一个半径为h的高维球区域sh,满足公式(2)的集合,x表示空间中样本点;
引入核函数和增加样本权重来解决(1)式中存在问题;在核函数中引入高斯核函数:
Figure FDA0003317183700000012
得到改进后的Mean Shift聚类算法形式:
Figure FDA0003317183700000013
其中,G(x)是单位核函数:
Figure FDA0003317183700000014
其中,H是带宽矩阵,即正定的n×n矩阵;(4)式中ω(xi)≥0是添加的每个样本的权重;带宽矩阵H为:
Figure FDA0003317183700000015
其中:h1,h2,hn表示带宽参数;
于是推出Mean Shift概率密度梯度归一化函数:
Figure FDA0003317183700000021
概率密度函数f(x)在d维空间内n个采样点xi,i=1,2,…;n,则f(x)的核密度为:
Figure FDA0003317183700000022
f(x)的梯度
Figure FDA0003317183700000023
的核密度为:
Figure FDA0003317183700000024
定义g(x)=-k′(x),G(x)=g(||x||2),推出:
Figure FDA0003317183700000025
由上式可得,Mean Shift向量就是右边第一个括号中与概率密度梯度函数成正比;则:
Figure FDA0003317183700000026
化简上式,令
Figure FDA0003317183700000027
Figure FDA0003317183700000028
当且仅当Mh(x)=0时,
Figure FDA0003317183700000029
得到新的圆心位置:
Figure FDA0003317183700000031
定义在以x为圆心,h为半径,总的核密度为:
Figure FDA0003317183700000032
其中,sh,hτ分别是半径为h、τ的高维球区域,C为常数,x、xi为样本点;
Figure FDA0003317183700000033
表示颜色信息,
Figure FDA0003317183700000034
表示空间位置信息;
3)内外膜定位,用改进hough变换圆检测算法,对左心室MRI图像内腔定位,得到左心室MRI图像内外膜分割初始轮廓;
优化Hough变换圆检测算法,在对图像进行找圆定位时,对输入图像进行阈值处理,像素较大点连成区域设为前景区域,累加器先累加前景区域的像素点;利用图像边缘信息进行约束,沿着梯度方向的像素点首先进行投票;
4)用双水平集模型分割左心室MRI内外膜;
用基于ODRLSE模型的0水平集和k水平集融合的双水平集模型,同时分割左心室MRI图像内外膜;双水平集模型的能量函数公式为:
Figure FDA0003317183700000035
其中:μ、λ、α均是大于零的常数;
Figure FDA0003317183700000036
模型距离正则项,
Figure FDA0003317183700000037
是模型边缘约束项,
Figure FDA0003317183700000038
是模型演化速度控制项;
Figure FDA0003317183700000039
保持双水平集模型演化过程中曲线的正则性;
Figure FDA00033171837000000310
融合AGVF约束0水平集和k水平集朝着内外膜边缘及边缘凹陷处演化;
Figure FDA00033171837000000311
约束模型曲线朝目标边缘速率;
定义
Figure FDA00033171837000000312
如下:
Figure FDA00033171837000000313
Figure FDA00033171837000000314
Figure FDA0003317183700000041
其中:p为势函数,p:[0,∞)→R,
Figure FDA0003317183700000042
表示模型轮廓,λ1、λ2、α1、α2均是大于零的常数;且δ和H分别为狄拉克函数和赫维赛德函数:
Figure FDA0003317183700000043
Figure FDA0003317183700000044
其中:ε设置的常数;
边缘指定符公式为:
Figure FDA0003317183700000045
式中:Gσ为标准偏差σ的高斯核函数,I的定义在整个图像区域;
双水平集模型将水平集演化导出为梯度流,用优化距离正则化项和外部能量来最小化能量泛函,从而驱动水平集向内外膜边缘位置移动;
5)对双水平集模型能量函数公式进行收敛性检查,若达到能量最小,分割结束;否则,重复步骤4);
6)得到左心室MRI图像分割内外膜轮廓分割效果图。
CN201811116609.XA 2018-09-25 2018-09-25 一种自动分割左心室内外膜的方法 Active CN109272512B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811116609.XA CN109272512B (zh) 2018-09-25 2018-09-25 一种自动分割左心室内外膜的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811116609.XA CN109272512B (zh) 2018-09-25 2018-09-25 一种自动分割左心室内外膜的方法

Publications (2)

Publication Number Publication Date
CN109272512A CN109272512A (zh) 2019-01-25
CN109272512B true CN109272512B (zh) 2022-02-15

Family

ID=65198043

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811116609.XA Active CN109272512B (zh) 2018-09-25 2018-09-25 一种自动分割左心室内外膜的方法

Country Status (1)

Country Link
CN (1) CN109272512B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110148111A (zh) * 2019-04-01 2019-08-20 江西比格威医疗科技有限公司 一种视网膜oct图像中多种视网膜病灶的自动检测方法
CN110232695A (zh) * 2019-06-18 2019-09-13 山东师范大学 基于混合模态图像的左心室图像分割方法及***
CN110288581B (zh) * 2019-06-26 2022-11-04 电子科技大学 一种基于保持形状凸性水平集模型的分割方法
CN112308845B (zh) * 2020-11-03 2021-07-02 赛诺威盛科技(北京)股份有限公司 左心室分割方法、装置及电子设备
CN113160116B (zh) * 2021-02-03 2022-12-27 中南民族大学 左心室内外膜自动分割方法、***及设备
CN112767420B (zh) * 2021-02-26 2021-11-23 中国人民解放军总医院 基于人工智能的核磁影像分割方法、装置、设备和介质
CN113379682B (zh) * 2021-05-21 2022-10-04 郑州大学 一种心脏mri影像耦合水平集分割方法及***
CN113436203A (zh) * 2021-06-02 2021-09-24 亳州学院 一种心室膜自动分割方法
CN113450256B (zh) * 2021-06-30 2023-03-03 重庆理工大学 基于水平集的ct图像主动脉自动分割方法与***
CN114240964B (zh) * 2021-12-09 2023-04-07 电子科技大学 一种心脏磁共振图像心肌区域分割方法及***
CN114693664B (zh) * 2022-04-13 2023-07-04 深圳北芯生命科技股份有限公司 血管超声图像的标注方法、设备及介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105741310A (zh) * 2016-03-21 2016-07-06 东北大学 一种心脏左心室图像分割***及方法
CN106952279A (zh) * 2017-03-31 2017-07-14 福州大学 一种超声心动图分割算法
CN108230342A (zh) * 2017-12-29 2018-06-29 电子科技大学 一种基于心脏解剖结构的左右心室水平集分割方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7715627B2 (en) * 2005-03-25 2010-05-11 Siemens Medical Solutions Usa, Inc. Automatic determination of the standard cardiac views from volumetric data acquisitions
JP2020510463A (ja) * 2017-01-27 2020-04-09 アーテリーズ インコーポレイテッド 全層畳み込みネットワークを利用する自動化されたセグメンテーション

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105741310A (zh) * 2016-03-21 2016-07-06 东北大学 一种心脏左心室图像分割***及方法
CN106952279A (zh) * 2017-03-31 2017-07-14 福州大学 一种超声心动图分割算法
CN108230342A (zh) * 2017-12-29 2018-06-29 电子科技大学 一种基于心脏解剖结构的左右心室水平集分割方法

Also Published As

Publication number Publication date
CN109272512A (zh) 2019-01-25

Similar Documents

Publication Publication Date Title
CN109272512B (zh) 一种自动分割左心室内外膜的方法
Xian et al. Automatic breast ultrasound image segmentation: A survey
CN108898606B (zh) 医学图像的自动分割方法、***、设备及存储介质
Li et al. Brain tumor detection based on multimodal information fusion and convolutional neural network
CN109522908B (zh) 基于区域标签融合的图像显著性检测方法
WO2021203795A1 (zh) 一种基于显著性密集连接扩张卷积网络的胰腺ct自动分割方法
CN109191471A (zh) 基于改进U-Net网络的胰腺细胞图像分割方法
CN111105424A (zh) ***自动勾画方法及装置
CN109300113B (zh) 一种基于改进凸包方法的肺结节辅助检测***及方法
CN111767952B (zh) 一种可解释的肺结节良恶性分类方法
CN113706487A (zh) 基于自监督特征小样本学习的多器官分割方法
CN112215842B (zh) 基于良性甲状腺模板的恶性结节边缘检测图像处理方法
CN110110727B (zh) 基于条件随机场和贝叶斯后处理的图像分割方法
Dong et al. Simultaneous segmentation of multiple organs using random walks
CN111784653A (zh) 基于形状约束的多尺度网络mri胰腺轮廓定位方法
CN107093184A (zh) 一种基于稀疏特征和形状相关性的超声图像序列分割方法
CN111127532B (zh) 基于深度学习特征光流的医学图像形变配准方法及***
Sivanesan et al. Unsupervised medical image segmentation with adversarial networks: From edge diagrams to segmentation maps
CN117036288A (zh) 一种面向全切片病理图像的肿瘤亚型诊断方法
Jiang et al. Segmentation of prostate ultrasound images: the state of the art and the future directions of segmentation algorithms
CN115147640A (zh) 一种基于改进胶囊网络的脑肿瘤图像分类方法
Afshar et al. Lung tumor area recognition in CT images based on Gustafson-Kessel clustering
CN113764101A (zh) 基于cnn的乳腺癌新辅助化疗多模态超声诊断***
CN117115437A (zh) 一种基于区域的多指标多器官医学图像分割模型评估***
CN107492101B (zh) 基于自适应构造最优图的多模态鼻咽肿瘤分割算法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant