CN101334895A - 一种针对动态增强乳腺磁共振影像序列的影像分割方法 - Google Patents

Abstract

一种针对动态增强乳腺磁共振影像序列的影像分割方法属于乳腺的磁共振影像处理技术领域，其特征在于含有：向计算机输入乳腺的磁共振三维横断面影像序列；把影像分为乳腺与空气界面，以及乳腺与胸腔界面两个部分；用动态阈值控制区域增长的分割处理得到乳房与空气的分界线；用同样方法得到乳腺与胸腔的初始轮廓；用控制水平集的方法得到乳房与胸腔复杂轮廓；再把分割结果拼接合，得到对一个时间点的三维磁共振影像序列，并把它作为下一组三维影像分割的初始位置等步骤。本发明提高了分割速度，解决了水平集算法不易确定初始轮廓和速度函数的问题，实现了复杂、大数据量的动态增强的乳腺磁共振影像的自动分割。

Description

一种针对动态增强乳腺磁共振影像序列的影像分割方法

技术领域

本发明方案属于医学影像的计算机辅助分析的应用领域，具体涉及一种对动态增强的乳腺磁共振影像序列的分割方法。

背景技术

乳腺癌是中老年女性最常见的恶性肿瘤之一，在中国已经高居女性恶性肿瘤的首位，并有逐年上升以及发病年轻化的趋势。目前对乳腺癌的预防尚无良策，因此早期的诊断和对肿瘤性质的判定是提高治愈率和降低死亡率的最有效途径。对乳腺影像的检测分析是成功检测、诊断、治疗乳腺肿瘤的重要依据，近二十年以来，在临床上广泛应用乳腺X线摄片对早期无症状的乳腺肿瘤进行检测，然而，受其成像原理的制约，在乳腺X线摄片中高密度的正常乳腺组织也会呈现与异常组织相接近的高亮度，对乳腺组织较为致密的年轻亚洲女性，其病灶不易被发现，易造成漏诊和误诊。据报道，由于在X光图像上，癌变组织与致密的纤维组织难以分辨，大约40％~50％的绝经前妇女的乳腺X光成像无法提供足够的信息排除乳腺癌的可能，在对经过乳腺保留治疗的妇女进行术后复发的肿瘤检测时，也会遇到同样的问题。

近年发展的影像增强的核磁影像技术，利用磁场和射频场对特定的增强剂进行成像，这种增强剂进入静脉后，将随着血流的增加，毛细血管的渗透性，以及细胞外容量的不同而不同程度地增加乳腺组织在影像上的对比度。由于这些对比度信息与乳腺瘤中的新的血管生成的信息是一致的，所以不同于X光的组织密度成像，影像增强的核磁成像是一种功能成像方式。研究表明，注入组织中的增强剂浓度分布随时间的变化情况与乳腺中的肿瘤病灶性质有特定的关系。对注射影像增强剂前后的乳腺组织影像进行对比，可以清楚地看到影像增强剂在组织内的充盈与洗脱的过程，因此可以对肿瘤组织的血供情况进行判定，从而对肿瘤组织的性质做出判断。目前，影像增强的磁共振影像是检测浸润性乳腺癌最为敏感的方式，敏感度可高达100％。而且它的敏感度不会被纤维组织削弱，在并不致密的乳腺组织中，这种技术也可以检测出隐秘在X光图像中的乳腺癌。因此影像增强的磁共振影像已经成为了乳腺X光成像检查的有力辅助。

对比增强MR成像的典型方法是基于T1加权的回波序列或梯度回波序列，在注射影像增强剂前后多次进行三维成像。通常，首先要获取高空间分辨率的三维磁共振乳腺影像，然后用设定好速率的压力注射器，通过静脉注入已知剂量的增强剂，之后用快速成像磁共振脉冲序列再获取多组高时间分辨率的序列图像。这种成像方式一般会在注射增强剂后产生多组时间间隔小于2秒，空间分辨率从几毫米到1厘米的双侧乳房的序列图像。

与乳腺X线摄片不同，动态影像增强的磁共振影像提供一种四维(4D)的信息，即三维的体积影像体数据，以及该影像三维体数据随时间变化的信息。因此，医生要观察的不再是乳腺X线摄片中所面对的每个病人的4张影像，而是高达每个病人数百张以上的影像。要在数百张影像中进行对比并做出判断，大大增加了医生的工作量。

采用动态增强的乳腺磁共振影像进行诊断的另一困难在于，磁共振影像的对比度和图像亮度值并不代表一个物理量。从制造者到医院中的使用，乳腺磁共振影像的密度和对比度尚没有统一的标准。并且，新的磁共振影像技术的引进(例如更高场强磁共振扫描仪的使用)，将进一步增强影像间的不一致性，给不同的成像数据对比造成困难。由于影像增强的乳腺磁共振影像可以提供的信息量非常巨大，随着成像技术的不断改进，由于时域和空间域分辨率的提高，以及各种辅助信息(如磁共振波谱)的增加，数据量还将不断增大。因此，对这种多维数据的自动预处理和有效地描述信息变得越来越重要，其中，计算机辅助的检测是技术发展的方向。

计算机对影像数据的处理可以分为两类工作：计算机辅助检测(分辨正常和怀疑的组织)和计算机辅助诊断(分辨怀疑组织中的良性和恶性肿瘤)。前者的目的是降低癌变组织未被检出的概率，后者的目的是降低良性肿瘤需要活检的概率。计算机分析的目的是从大量的数据中抽取有用信息，其目标是通过降低相互间对比的难度来增强信息检测的特异性和一致性。本发明的目标主要是解决计算机辅助检测中对乳腺组织的定位问题，通过影像分割处理，去除成像过程中在影像数据中出现的空气、骨骼和胸部脏器等非乳腺组织，以便在进一步的处理中，仅仅对比和处理乳腺组织的影像亮度和对比度信息，降低影像进一步处理的难度。因此本发明的核心技术是提出了一种从大量影像中快速有效地分割出乳腺组织的影像处理方法。

乳腺MRI图像的分割是为了去除两个方面的影响：一方面，成像系统和数据采集过程都会形成噪声点，在进行不同时间上影像的对比时，这些随机噪声被凸显出来，干扰对影像信息的判断；另一方面，胸腔中内脏的活动也会组成不同时间点上影像的差异，影响后续的对比分析。因此要将影像中不属于乳腺组织的像素分割出来，不再进行对比。目前，针对乳腺影像的分割，已经提出了多种方法。其中，Saskia van Engeland在他的研究中提到一种粗略计算乳房位置的方法，他在二维横切图中，用交互方式选取垂直于乳房的平面，在此基础上进行分割。这种算法原理简单，效果稳定，但需人机交互，且其分割误差较大，分割出的乳房部分中缺少了腋下的淋巴部分，而在临床中这部分组织对于判别肿瘤性质很重要。Lina Arbach提出一种基于边界增强的分割方法，将图像分为三部分(空气，乳房，胸腔)，通过线性滤波器，对皮肤与空气，以及胸壁部分进行增强，以实现分割，但在处理复杂影像时，可靠性和速度尚有待提高。另外，对于序列图像的分割，Yoo等人采用区域自适应Snake模型与三维显示相结合的方法对耳蜗影像进行了分割，他们首先对二维切片图像进行分割，并把得到的分割结果作为其邻接切片的初始轮廓，然后将二维分割的结果堆叠成三维数据，并对其进行三维显示，是二维分割与三维显示的结合。但是，这种方法只在相邻图像切片上组织变化不大时有效。在我们处理的在高时间分辨率的乳腺图像序列中，影像数据的空间分辨率较低，相邻图像切片组织的变化较大，因此难于采用这种方法。

综上所述，动态增强的乳腺磁共振影像的分割存在以下几个难点：首先，轮廓较为模糊的单一时间点上的三维数据分割难以自动实现，常常需要通过三维显示和人机交互来进行，增加了工作量；第二，由于影像数据量大，造成分割运算的时间非常长；第三，在影像中胸腔与乳腺的分界不清晰，而且胸腔内脏器的搏动，严重干扰对不同时间点上影像的对比分析，因此将乳腺与胸腔准确分割开十分困难，也十分重要。

为对上述动态增强的乳腺磁共振影像进行分割处理，我们发明了一种多步骤、分层次的影像处理方法，针对乳腺在影像中的形态特点，将每张断层图像进行分区，然后对分区得到的子图像分别进行处理。对分区得到的仅包含乳腺与空气清晰界面的子图像，我们采用简单的阈值控制的区域增长方法进行分割，而对于乳腺与胸腔的模糊边界，我们采用基于水平集(Level set)计算的图像分割方法，迭代地寻找乳腺与胸腔的边界。上述处理方法所形成的处理步骤，构成了一种针对动态增强乳腺磁共振影像序列的影像分割方法与系统。

水平集(Level set)方法主要是从界面传播等研究领域中逐步发展起来的一种数学处理方法，它是处理封闭运动界面随时间演化过程中几何拓扑变化的有效的计算工具。依赖于时间的运动界面的水平集描述首先由Osher和Sethina提出。其主要思想是将移动的界面作为零水平集嵌入高一维的水平集函数中，在图像空间上构造一个静态的、间隔均匀的网格，网格点上的值代表网格点到轮廓，即零水平集的最短距离，这可以理解为地图上的等高线。这样，由闭超曲面的演化方程可以得到水平集函数的演化方程，而嵌入的闭超曲面总是其水平集，最终只要确定零水平集即可确定移动界面演化的结果。Level Set方法自提出以来，已在图像处理和计算机视觉等领域得到了广泛的应用：Sethian运用Level Set进行图像去噪；Mallad应用Level Set进行图像分割；Parogios和Deriche将其应用于纹理分割以及运动目标分割和跟踪等等。其主要优点是可以对任意复杂的形状进行模型化并且隐式地解决拓扑形状的分裂、合并等变化。另外该算法中与曲率相关的速度项可以控制曲线或曲面的光滑特性，因此对于分割具有光滑表面的组织，该算法更加具有优势。

具体讲，Level set方法是一种跟踪曲线或曲面轮廓演化的数学方法。它首先将三维曲面隐含地表达为高维函数的零水平集，即具有相同函数值的点集，通过level set函数的进化来隐含地求解曲面的运动，得到所需要的轮廓。

在处理中，假设随时间演化的零水平集为：

(x，y，t)＝0    (1)

对方程(1)两边求关于时间的偏导数，有

$\frac{\∂\mathrm{\Φ}}{\∂t}+\▿\mathrm{\Φ}(x,y,t)\frac{\∂x}{\∂t}\frac{\∂y}{\∂t}=0---\left(2\right)$

假设F为外法线方向的速度，那么：

$F=\frac{\∂x}{\∂t}\frac{\∂y}{\∂t}*\frac{\▿\mathrm{\Φ}}{|\▿\mathrm{\Φ}|}---\left(3\right)$

因此，我们便得到水平集迭代的基本方程式

${\mathrm{\Φ}}_t+F|\▿\mathrm{\Φ}|=0---\left(4\right)$

其中，Φ为水平集函数，Φ_t为其对时间的导数，Φ的零水平集Γ就表示目标的轮廓，即

Γ(t)＝{x，y|φ(x，y，t)＝0}    (5)

式中表示Level set函数的梯度范数。F为轮廓法线方向上的速度函数，用来控制轮廓的演变运动，当速度为零的时候，轮廓不再演变，就得到了分割的结果。因此，在图像分割中关键的问题是确定合适的速度函数F，以控制水平集迭代中轮廓的演化。一般情况下速度函数F的设计是考虑将图像的特征作为控制曲线边界演变的要素，使得曲线边界向所要寻找的图像特征靠近，因此速度函数中可以包括与图像有关的项(如图像的灰度，灰度梯度等)，以及与曲面的几何形状有关的项(如边界曲面的曲率等)，因此，通过构造水平集演化的速度函数可以将与图像灰度、边界形态等相关的多种因素考虑进去。常用的速度函数一般可以描述为

F＝α·P(X)+β·K(X)    (6)

式中：P(X)为与图像灰度梯度相关的传播速度项，K(X)为与曲率相关的速度项，用来控制轮廓的光滑情况，大小与曲率成正比，α、β为两者的权重系数。

在本发明中采用一种基于阈值区间的三维水平集算法，将图像的灰度值应用一种符号映射函数变换为介于[-1，1]区间的特征图像，作为水平集算法速度函数的一个输入参数，控制轮廓的演化速度。通过选取(6)式中传播速度P(X)的大小与特征图像的函数值比例系数，当特征图像的像素值大于零时，P(X)垂直轮廓向外，反之向内，从而使水平集计算的叠代过程收缩到特定的图像灰度附近的轮廓。

在应用Level Set方法时，有两个突出的问题需要解决：

1.确定初始轮廓，即初始的水平集。

2.确定速度函数。

本发明针对这两个问题，设计了一种自动确定初始轮廓和速度函数的方法。

本发明针对所分割的动态增强乳腺磁共振影像序列的特征，提出了一种从随时间变化的动态增强乳腺影像序列中自动分割出乳腺组织的处理策略，主要包含以下处理步骤：首先，应用分区处理的方法，降低处理对象的复杂性。方法是通过对影像中身体的位置进行估算，并据此将横断面影像分成前后两部分，前面的一半仅包含乳房和空气的影像；然后，通过对这部分影像亮度的统计，得到关于乳腺组织在影像灰度的先验信息；接着，采用基于上述灰度信息，并进行位置信息加权的方法获取水平集的初始轮廓。同时，采用上述灰度信息生成控制水平集演化(膨胀或收缩)的初始化特征图像，从而解决了上述两个应用水平集算法的关键问题；接着，应用基于阈值区间的3D水平集算法分割后一般影像中乳房与胸腔间的复杂轮廓；最后，利用动态增强乳腺磁共振影像数据时间分辨率高的特点，利用对前一时间点上三维影像分割的结果，作为后一时间点上三维影像分割的初始位置，来处理不同时间点上的3D影像序列，从而实现了动态增强乳腺磁共振影像的快速、自动的分割。

发明内容

本发明的目的在于提供一种针对动态增强乳腺磁共振影像序列的影像分割方法，该方法的处理流程可以快速、准确地将动态增强乳腺磁共振影像序列中的乳腺组织分割出来，一方面可以帮助医生在影像对比中排除不相干的像素的干扰，更重要的是可以作为动态增强乳腺磁共振影像序列计算机辅助检测的预处理，有助于提高对乳腺组织肿瘤样病变检测、定位和诊断的效率和精准程度。本发明还提供了实现该方法的系统。

本发明提供的针对动态增强乳腺磁共振影像序列的影像分割方法，其处理步骤包括：

步骤1：输入一个时间点的三维磁共振影像序列，可能包含24～120张切片图像；

步骤2：对横断面的影像进行分区处理，以降低影像的复杂性。方法是有计算机自动对影像中身体的位置进行搜寻，并据此将横断面影像分成前后两部分，前面的一半仅包含乳房和空气背景的影像；

步骤3：采用图像处理中动态阈值控制的区域增长方法，对步骤2得到的前半部分影像进行分割处理；

步骤4：通过对步骤3分割出的前半部分影像中乳腺组织像素灰度的统计，得到关于乳腺组织磁共振影像灰度的先验信息；

步骤5：对步骤2划分出的后半部分影像进行位置相关的灰度加权处理，并基于上述估计的乳腺组织影像灰度信息，进行区域增长，得到乳腺与胸腔间的初始轮廓；

步骤6：利用步骤4估计的乳腺组织影像灰度信息，生成一个控制水平集演化(膨胀或收缩)的特征图像；

步骤7：将步骤5得到的初始轮廓作为水平集初始值，利用步骤6得到的特征图像，用来构造水平集的速度函数，在此基础上用基于阈值区间的三维水平集算法分割步骤2得到的后一半影像中乳房与胸腔间的复杂轮廓；

步骤8：将步骤3和步骤7得到的前半部分与后半部分的图像分割结果进行拼接合成，得到对一个时间点的三维磁共振影像序列的分割结果；

步骤9：输入下一时间点上的三维磁共振影像序列，根据动态增强乳腺磁共振影像数据高时间分辨率的特点，利用对前一时间点上三维影像分割的结果，作为后一时间点上三维影像分割的初始位置(即取代上述步骤5所估算的初始轮廓)，并对所输入的新的影像序列顺次进行步骤2、步骤3、步骤4、步骤6、步骤7的处理、来得到新的时间点上的三维影像序列分割结果，直到将所有时间点上的影像序列分割完毕。

步骤10：将分割结果进行显示，以供医生进行分析，或将最终分割的结果存储在计算机硬盘中，供影像识别软件进行进一步的分析处理。

本发明提供的针对动态增强乳腺磁共振影像序列的影像分割系统，包括图像输入模块、图像分区模块、图像区域增长处理模块、像素灰度统计计算模块、水平集初始化计算模块、水平集迭代计算模块和图像输出模块；

图像输入模块用于接收从MRI影像工作占传来的待分割的乳腺MRI影像，对影像进行排序和分组，并传送给图像分区模块；

图像分区模块用于对输入的图像进行计算，得到图像中人体胸骨前的位置，并将图像分成前、后两个部分(图像)，以便分别处理。分区后的图像分别送入图像区域增长处理模块和水平集初始化计算模块；

图像区域增长处理模块用于对图像分区模块送入的前半部分影像进行分割，通过图像处理的区域增长算法，得到前半部分图像中乳腺组织的像素集合，即前半部分图像的分割结果，并传送给像素灰度统计计算模块；

像素灰度统计计算模块用于统计计算前半部分图像中乳腺组织像素的灰度分布，包括迅速灰度的平均值和分布的标准差，并将计算结果传送给水平集初始化计算模块；

水平集初始化计算模块用于根据前面处理模块送来的乳腺迅速灰度统计结果，以及分区模块送来的后半部分图像，分别计算乳腺与胸腔的初始轮廓和控制水平集演化的特征图像；

水平集迭代计算模块用于对后半部分图像进行优化计算，搜寻乳腺与胸腔的轮廓，其初始值由水平集初始化计算模块得到。

图像输出模块用于对图像区域增长处理模块和水平集迭代计算模块得到的两部分分割的结果进行拼接，并将分割经过显示在计算机屏幕上供医生进行分析，或存储到硬盘，供图像分析软件进一步进行分析。

本发明的主要特点为：

1)算法框架完整，整个算法从一个时间点上的三维影像入手，根据影像特点对影像进行简化处理，并利用处理中得到的关于影像灰度分布的特点，逐步构造搜索复杂边界的算法，并推广到不同时间点上的三维数据体。算法较好地解决了复杂、大数据量的动态增强的乳腺磁共振影像的自动分割问题。

2)所构造的多步骤处理的全套流程，减少了每一步处理的复杂性，使分割速度大大加快，分割算法的稳定性得到提高。

3)将前一步处理得到的图像灰度信息用于估算后一步的初始轮廓，并控制水平集的演化，解决了水平集算法应用中的确定初始轮廓选取和确定速度函数的难题。

附图说明

图1为本发明从MRI输入图像数据和输出分割结果的流程图。

图2为本发明的处理模块的框图。

图3为本发明针对动态增强乳腺磁共振影像序列的影像分割方法的流程图。

图4为实施例中确定影像中身体位置方法的示意图；4(a)为自动搜寻确定图像中分界面的位置示意图，4(b)为乳腺MRI图像沿箭头方向的灰度变化

具体实施方式

下面结合附图和实例对本发明的方法作进一步详细的说明。

如图1所示，本发明所说的输入影像序列是从医用磁共振扫描仪的影像工作站获取的，影像工作站将对被检查者的动态增强的乳腺磁共振影像传送到本发明所说的图像分割系统运行的计算机上，由本发明的针对动态增强乳腺磁共振影像序列的影像分割方法与系统进行处理，得到对序列图像的分割结果。

为实施本发明的处理方法，需要构建如图2所示的处理模块，包括图像输入模块、图像分区模块、图像区域增长处理模块、像素灰度统计计算模块、水平集初始化计算模块、水平集迭代计算模块和图像输出模块。送入的动态增强乳腺磁共振影像序列经过这些模块的处理，可以得到对影像数据集的分割。这些模块的工作需要按照下面的处理步骤逐次进行。

如图3所示，本发明方法包括以下步骤：

(1)从医用磁共振影像工作站输入一个时间点的三维磁共振影像序列。根据检测中的实际情况，一个时间点上获取的动态增强乳腺磁共振三维影像序列中可能包含16～32张横断面影像，或120张左右的人体冠状面的切片影像。

(2)对人体横断面的影像进行分区处理，以降低影像的复杂性。

要对整幅影像进行分割处理，不仅存在大量冗余信息，还容易引起错误的分割，其处理速度也很慢。考虑到动态增强乳腺磁共振成像过程中人体的位置是统一的，病人俯卧在磁共振检查线圈之上。因此形成的横断面影像的形态特征非常明显，影像中下半部分主要呈现了双侧的乳房的横断面，而上半部分则主要呈现人体胸腔部分的横断面。因此我们选取横断面影像入手进行图像的分区处理，利用计算机自动计算和寻找分区的分界位置并将影像分为两部分，使前半部分影像只需要对乳房与空气的分界面进行分割。后半部分则针对乳房与胸腔的不明显分界面进行分割处理。

一幅典型的动态增强乳腺磁共振数据的横切面断层影像如图4a)所示。为了让计算机自动划分图像的区域，沿影像的上下方向，在两侧乳房的中间位置寻找人体组织边界。其方法是沿着图4(a)中箭头所示的位置，统计沿此方向的影像灰度变化，图4(b)给出了沿该直线方向影像像素灰度变化的示意。可以看到空气所对应区域的图像灰度值很低，而人体组织对应像素的灰度快速增长。在图4(a)中顺着箭头方向搜寻灰度的变化率超过给定阈值的位置，就可以确定人体胸骨前组织的边界，并可以将其作为划分两部分影像区域的位置。如图4(b)中横向直线所示。横线以下为前半部分待分割影像区域，即待分割图像(1)，横线以上为后半部分待分割影像(2)。对所有输入的横断面影像均进行这样的划分。

上面所说的给定阈值是这样确定的，首先给定一个大大高出空气所对应区域影像灰度的像素灰度阈值，例如说50。然后沿着图4(b)中曲线从左到右进行搜索，找到该阈值第1次出现的位置，该位置代表从空气中搜索到人体组织的大致位置。然后在从左侧开始到这一位置的区间上，逐次计算两两相邻像素间的灰度变化，即：

I(k)＝I(k+1)-I(k)    (7)

其中k代表沿搜索方向顺次取得的像素位置，而I(k)代表该点的像素灰度。

接下来，对得到的一系列I(k)的值进行排序，并取其中第2大的I(k)值作为给定阈值，顺着图4(a)中箭头方向搜寻相邻像素灰度的变化率超过该给定阈值的位置，作为身体的位置。这样处理的好处是可以非常可靠地获得影像中人体的前后位置。

(3)采用动态阈值控制的区域增长方法对步骤2得到的前半部分影像进行分割处理。

针对乳腺磁共振影像的特点，根据图4(b)中得到的空气中的噪声分布选定阈值。在上一步骤搜索人体组织位置的区间(即计算计算I(k)的区间)内，统计空气区间内所有像素的灰度分布(均值与方差)，然后应用动态阈值控制的区域增长方法，对待分割图像(1)中的乳房与空气进行分割，从而精确地保留乳房部分组织的影像。动态阈值控制的区域增长是指：首先在图像中选定一个种子像素点，然后逐次遍历种子点的图像邻域，检查各个像素的灰度值，如果某个像素处于规定的灰度阈值之内，则将该像素归于与种子点同类，并作为新的种子点，进行一次新的邻域像素寻找过程。否则，将该位置像素设置为背景。通过对所有种子点的遍历，就得到了整个图像的分割结果。本文所说的动态阈值控制的区域增长方法的特点在于：种子点的选择是在步骤2中计算I(k)的区间内，选择区间的中点作为种子像素，并将前面统计得到的空气区间内所有像素的灰度的平均值(mean)和标准差(std)作为区域增长过程中规定灰度阈值的控制条件，对待分割图像(1)的像素进行分类处理。具体的算法为：

从种子点出发，对种子点的图像邻域，逐个检查各个像素的灰度值，如果其灰度与上述灰度平均值不超过三倍的标准差，则认为该像素与种子点属于同一类，否则，认为该像素属于另外一类。这样，最终得到的像素集合可以表示为：

R_air＝{(x，y)|I(x，y)∈mean±3·std}    (8)

其中，R_air表示空气区域像素的集合，I(x，y)表示位于(x，y)位置的像素的灰度，mean和std分别表示统计计算得到的空气区间内像素的灰度均值与方差。这样通过自动确定算法所需要的初始种子点和动态阈值，对不同条件下采集的动态磁共振图像，可以达到稳定的分割结果。

(4)通过对步骤3分割出的乳腺组织亮度的统计，得到关于乳腺组织磁共振影像灰度的先验信息。在待分割图像(1)的分割处理完成后，对所分割出的乳房部分的影像，进一步统计其灰度的均值和方差，从而得到对乳房组织在图像中的灰度分布的估计，作为后续分割处理算法的依据条件。

(5)基于上述估计的乳腺组织影像灰度信息，进行影像中乳腺组织的区域增长计算，得到乳腺与胸腔间的初始轮廓。

在步骤4中，获得了待分割图像(1)中乳房组织的灰度均值和方差，我们以此作为确定胸腔初始轮廓的灰度依据，以处理待分割图像(2)，即图4(a)中横线以上的部分影像。我们采用区域增长的方法获得乳房与胸腔分界面的初始轮廓，其方法是：将种子点放在乳房区域内，即后半部分图像的左下角边界附近，选择一个灰度与乳房灰度均值接近的点，如图4(a)中的点A。选择区域增长的阈值范围为{乳房组织灰度均值±2×乳房组织灰度方差}，获得增长得到的乳房部分的二值图像。对所有输入的横断面影像均进行这样的增长，从而得到由所有横断面切片图像共同构成的三维影像数据集的乳腺与胸腔的初始轮廓。

(6)利用步骤4估计的乳腺组织影像灰度信息，生成控制水平集演化(膨胀或收缩)的特征图像。

在水平集图像分割算法中，边界轮廓的分裂与合并都隐含在了水平集函数的演变中。根据步骤4获得的乳腺部分的像素灰度值，可以自动设置一个阈值区间，从而生成控制水平集演化(膨胀或收缩)的特征图像。基于阈值区间的特征图像的获取是基于所统计的乳腺组织的灰度区间。当像素值位于灰度上、下限[L，U]区间之内时，将其映射到[0，1]的特征幅度，称为前景区域；而位于灰度上、下限[L，U]区间之外的像素值则映射到[-1，0]的特征幅度，称为背景区域。因此基于阈值区间的特征图像的像素值位于[-1，1]区间之内，当某个位置的特征图像的像素值大于0时，选择速度函数使水平集的轮廓扩张，否则使水平集的轮廓收缩。这样，通过特征图像来控制速度函数分别对应动态轮廓不同的传播方向，动态轮廓既可以扩张也可以收缩，降低了水平集搜索中对选择初始轮廓的准确性要求，并且可以适应复杂的轮廓形状，初始轮廓甚至可以部分地位于感兴趣目标之外，而不会造成搜索的失败。

具体讲，利用步骤4获得的乳腺灰度分布的先验信息(乳房组织灰度均值和乳房组织灰度方差)，把乳腺组织所在的灰度区域划为前景，其他灰度的像素则划为背景，生成特征图像，特征图像的映射公式如下：

式中，g(x，y)为原始图像的像素灰度，阈值范围L，U的选择是以步骤4获得的乳房灰度均值及方差为依据的，实际试算中取

U＝255    (10)

L＝乳房灰度均值-2×乳房灰度方差

(7)将第5步骤得到的各图像层片的胸腔初始轮廓作为水平集的初始值，利用步骤6得到的各图像层片的特征图像，构造水平集的速度函数，在此基础上用基于阈值区间的三维水平集算法分割步骤2得到的后一半影像中乳房与胸腔间的复杂轮廓。

具体讲，我们选取(3)式中传播速度P(X)的大小与特征图像值F_prop成正比，当特征图像的像素值大于零时，取P(X)的方向垂直于轮廓向外，反之垂直于轮廓向内，而P(X)幅度与F_prop成正比。从而使水平集计算的迭代过程收缩到特定的图像灰度附近的轮廓上。为进行水平集的演化，还需要设置(6)式中的传播速度和曲率速度的权重系数，用于速度函数的计算，这里，分别取1，50。因此，(6)式演化为：

F＝F_prop+50*K    (11)

其中，K为演化中轮廓上每一点的曲率。根据轮廓的表达式(5)可以计算出曲率 $K=(\frac{\∂\mathrm{\Γ}}{\∂x}\frac{{\∂}^2\mathrm{\Γ}}{\∂y^2}-\frac{\∂\mathrm{\Γ}}{\∂y}\frac{{\∂}^2\mathrm{\Γ}}{\∂x^2})/[{\left(\frac{\∂\mathrm{\Γ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\Γ}}{\∂y}\right)}^2{]}^{3/2},$ x，y为图像的像素点的坐标位置，Γ为各次迭代得到的零水平集，即目标的轮廓。

设定水平集计算的迭代在两种情况下会终止：(a)分割次数到达设置的最大迭代次数；(b)迭代次数增加而轮廓不变或变化非常小。将水平集计算迭代终止得到的轮廓，作为乳腺与胸腔分割的轮廓。从而得到乳腺组织在待分割图像(2)中的分布。

(8)将步骤3得到的前半部分乳腺影像的分割结果，与步骤7得到的后半部分的乳腺影像分割结果进行拼接合成。其方法是进行直接的图像合并，把待分割图像(1)的分割结果和待分割图像(2)的分割结果按照原来划分的位置放回原处，就可以得到对一个时间点上三维磁共振影像序列的分割结果。

(9)输入下一时间点上的又一组三维磁共振影像序列，根据动态增强乳腺磁共振影像数据高时间分辨率的特点，利用对前一时间点上三维影像分割的结果，作为后一时间点上三维影像分割的初始位置(即取代上述步骤5所估算的水平集初始轮廓)，并对所输入的新的影像序列顺次进行步骤2、步骤3、步骤4、步骤6、步骤7的处理、来得到新的时间点上的三维影像序列分割结果。重复这种处理，直到将所有时间点上的影像序列分割完毕。

(10)将分割结果进行显示，以供医生进行分析，或将最终分割的结果存储在计算机硬盘中，供影像识别软件进行进一步的分析处理。

实际影像数据实验的效果：

我们采用1.5T全身磁共振成像仪获取的快速成像人体乳腺增强影像数据进行实验试算，影像的空间分辨率为512*512*24，像素物理间隔为0.625*0.625*5mm³。

在本发明的影像分割计算流程中，最耗时的计算是水平集的迭代计算，我们采用的基于阈值区间的三维水平集迭代算法，其轮廓既可以膨胀也可以收缩，初始轮廓甚至可以部分地位于感兴趣目标之外，因此算法降低了对初始轮廓选择准确性的要求。我们对选取不同的初始轮廓得到的分割结果进行对比。为测试算法，首先人为地在两侧乳房中各选择一个种子点，形成两个球面形初始轮廓，实验中设定了三种不同的半径进行对比，最终水平集的搜索均能够收缩到所要寻找的乳腺与胸腔的分界。可见，采用基于阈值区间的三维水平集算法进行分割时，初始轮廓的选择更具有灵活性，避免了因初始轮廓引起的分割错误。

本发明所说的基于影像灰度和位置信息所获取的初始轮廓，更接近最终的实际边界，这样可大大减小计算中的迭代次数，并避免迭代运算陷入局部极小而失效。我们用10组三维影像数据进行了对比试验。表1为初始轮廓不同时迭代次数的对比。其中，左列为选择一种中等半径的球面形初始轮廓时的迭代次数，右列为按照本发明的方法设置初始轮廓时的迭代次数。选择球面作为初始轮廓时，迭代次数通常为1000步左右，基于灰度和位置信息的方法获取初始轮廓时的迭代次数仅为13步左右。因此本发明的方法极大地提高了分割的速度。

表1不同初始轮廓时迭代次数对比

本发明通过对动态增强乳腺磁共振影像序列的影像的细致分析，发明了一种针对动态增强乳腺磁共振影像序列的影像分割方法。其特点是计算机处理的自动化程度高，方法的可靠性和精确程度都较高，针对随时间变化的动态增强乳腺影像序列的特点，实现了动态增强乳腺磁共振影像的快速、自动的分割。本发明的实现并不局限与上述实例所公开的范围，也可以按照本发明的步骤，采用不同于上述实例的方式来实现上述技术方案。

Claims (1)

1、一种针对动态增强乳腺磁共振影像序列的分割方法，其特征在于，所述方法是在计算机中依次按以下步骤实现的：

步骤(1)在所述计算机中，创建动态增强乳腺磁共振影像序列的分割系统，包括图像输入模块、图像分区模块、图像区域增长处理模块、像素灰度统计计算模块、水平集初始化计算模块、水平集迭代计算模块和图像输出模块；其中，

图像输入模块，用于接收从磁共振MRI影像工作站传来的待分割的乳腺MRI影像，按照扫描序列不同时间段与空间位置对影像进行分组和排序，并传送给图像分区模块；

图像分区模块，用于对输入的图像进行计算，得到图像中人体胸骨前的位置，并将图像分成前、后两个部分图像，以便分别处理，分区后的图像分别送入图像区域增长处理模块和水平集初始化计算模块；

图像区域增长处理模块，用于对图像分区模块送入的前半部分影像进行分割，通过图像处理的区域增长算法，得到前半部分图像中乳腺组织的像素集合，即前半部分图像的分割结果，并传送给像素灰度统计计算模块；

像素灰度统计计算模块，用于统计计算前半部分图像中乳腺组织像素的灰度分布，包括灰度的平均值和方差，并将计算结果传送给水平集初始化计算模块；

水平集初始化计算模块，用于根据前面处理模块送来的乳腺迅速灰度统计结果，以及分区模块送来的后半部分图像，分别计算乳腺与胸腔的初始轮廓和控制水平集演化的特征图像；

水平集迭代计算模块，用于对后半部分图像进行优化计算，搜寻乳腺与胸腔的轮廓，其初始值由水平集初始化计算模块得到；

图像输出模块，用于对图像区域增长处理模块和水平集迭代计算模块得到的两部分分割的结果进行拼接，并将结果显示在计算机屏幕上供医生进行分析，或存储到硬盘，供图像分析软件进一步进行分析。

步骤(2)从医用磁共振影像工作站输入用一次三维扫描获取的动态增强乳腺磁共振三维影像序列中16-32张人体乳腺区域横断面影像；

步骤(3)对横断面的影像进行自动的分区处理，得到分别包含乳腺与空气界面和乳腺与胸腔界面的两部分图像，依次称为待分割影像(1)和待分割影像(2)：

步骤(3.1)以影像上人体两侧乳房的中间位置为寻找路径，顺次统计从身体前方到身体后方垂直方向上图像各像素点的灰度及其相邻像素的灰度值，得到一条以所述寻找方向上的经历的各像素点为横坐标，该各个像素的灰度值为纵坐标的曲线，

步骤(3.2)给定一个高于空气所对应区域影像灰度的像素灰度阈值，然后在所述曲线上从左到右寻找一个像素灰度值超过规定阈值的位置，并将曲线左侧开始位置到该位置设定为搜寻区间，在此区间内寻找灰度变化率序列中取值第2或第3大的位置作为图像中人体组织与空气的分界点，所述的灰度变化率序列为δI(k)＝I(k+1)-I(k)，其中，k为沿所述寻找方向顺次取得的像素点位置，I(k)为k点位置上的像素灰度，

步骤(3.3)通过步骤(3.2)所述分界点作一条沿人体左右方向的直线，把待分割的影像划分为两个子图像，该直线前方为待分割影像(1)，直线后方为待分割影像(2)；

步骤(4)采用图像处理中动态阈值控制的区域增长方法，对步骤(3)得到的待分割影像(1)进行分割处理：

首先，在步骤(3.2)所述的搜寻区间上，逐次统计各个像素点的灰度，得到像素的灰度均值mean和像素灰度的方差std，它们对应的是成像过程中空气区域的成像噪声特性，然后规定区域增长的动态阈值为上述均值mean的上下三倍方差std范围内，即mean±3·std，通过区域增长得到待分割影像(1)中对应空气的区域，而剩余的区域为乳腺组织的图像区域；

步骤(5)对步骤(4)得到的待分割影像(1)中的乳腺区域，计算乳腺组织像素点灰度的均值与方差；

步骤(6)对待分割影像(2)进行处理，方法是基于上述估计的乳腺组织影像灰度信息，进行区域增长，得到乳腺与胸腔间的初始轮廓：

步骤(5)中，获得了待分割影像(1)对应的乳房组织的灰度均值和方差，以此作为确定胸腔初始轮廓的灰度依据，用来通过图像处理中区域增长的方法获得待分割影像(2)中乳腺与胸腔间的初始轮廓：

首先在待分割影像(2)的乳房区域选择种子点，即在待分割影像(2)的左下角边界附近的100×100像素区域内，选择一个灰度值与上述乳房组织灰度均值最接近的像素点，并规定区域增长的阈值范围为乳房组织灰度均值±2×乳房组织灰度方差，然后利用区域增长算法获得待分割影像(2)中乳房部分的区域，对所有输入的横断面影像均进行这样的区域增长运算，得到由所有横断面切片图像共同构成的三维影像数据集的乳腺与胸腔的初始轮廓，作为后续水平集计算的初始轮廓；

步骤(7)利用步骤(5)得到的乳腺组织像素灰度信息，生成一个控制水平集演化速度的特征图像：

步骤(4)获得的待分割影像(1)的乳腺组织像素灰度值，被用来自动设置一个阈值区间，生成控制水平集演化的特征图像，在待分割影像(2)上，当像素值位于规定的灰度上、下限[L，U]区间之内时，将其映射到[0，1]的特征幅度，而位于灰度值上、下限[L，U]区间之外的像素值则映射到[-1，0]的特征幅度。

利用步骤(5)所获得的乳腺灰度信息，即乳房组织灰度的均值和乳房组织灰度的方差，把乳腺组织灰度均值-2×乳房组织灰度的方差作为下限L，把乳腺组织灰度均值+2×乳房组织灰度的方差作为上限U，取U＝255，生成一个特征图像，计算特征图像的映射公式为：

其中：U＝255，L＝乳房灰度均值-2×乳房灰度值方差，g(x，y)为待分割影像(2)的像素灰度值；

步骤(8)将所得到横断面切片二维图像按照扫描位置重构成的三维影像数据集作为待处理三维数据集，步骤(6)得到的初始轮廓按照垂直于二维图像平面方向扫描位置构成三维轮廓，作为水平集初始值，利用步骤(7)得到的对应各个待分割影像(2)的特征图像构造三维的水平集演化的速度函数F，采用图像处理中基于阈值区间的三维水平集图像分割算法，分割通过步骤(3)所得到的子图像2中乳房与胸腔间的复杂轮廓，F＝F_prop+50K，F_prop特征图像，K为演化中轮廓上每一点的曲率；

步骤(9)将步骤(4)和步骤(8)得到的两个子图像的分割结果进行拼接合成，得到对一个时间点的三维磁共振影像序列的分割结果；

步骤(10)输入下一组动态增强的三维磁共振影像序列，并利用前一组三维影像分割的结果，即取代上述步骤(6)所估算的胸腔初始轮廓，构成后组三维影像分割的初始位置。然后对所输入的新影像序列顺次进行步骤(3)、步骤(4)、步骤(5)、步骤(7)、步骤(8)和步骤(9)的处理、来得到新的三维影像序列分割结果，再重复步骤(10)，直到将所有动态增强的磁共振影像序列分割完毕；

步骤(11)将分割结果进行显示，以供医生进行分析，或将最终分割的结果存储在计算机硬盘中，供影像识别软件进行进一步的分析处理。

Images