CN103426169B - 一种医学图像的分割方法 - Google Patents

Abstract

本发明公开了一种医学图像的分割算法，涉及数字图像处理技术。本发明的方法具体步骤如下：选择感兴趣区域，布置初始聚类中心点，聚类分割，连通性调整。本发明通过对图像感兴趣区域及附近区域进行分割，将其划分为一系列同质小区域。使得后续过程能够直接在同质区域级别上操作而不是在单个像素点上操作，大大降低了计算量，执行速度快且结果准确。本发明特别适用于计算机控制X射线断层扫描装置或医学放射性影像处理系统中对组织进行定位、测量、识别或分类的操作。

Description

一种医学图像的分割方法

技术领域

本发明涉及数字图像处理技术领域，具体涉及一种医学图像的分割方法。

背景技术

X线计算机断层摄影(Computed tomography，CT)，是用X射线对人体指定部位一定厚度的层面进行扫描，由探测器接收透过该层面的X射线，转变为可见光后，由光电转换变为电信号，再经模拟/数字转换器转为数字，以数值点阵图像的方式存储在计算机中，这就是通常所说的CT图像。CT图像是医学诊断的重要工具，目前的技术水平还无法实现疾病的全自动诊断，医师需要通过观察CT图像对病灶进行诊断。在这种前提下，有效的图像处理算法对于医生能够更加方便的观察与分析CT图像，不仅能够加快医师的诊断速度，提高工作效率，也能够增加医师对疾病的诊断精度。

一套标准CT数据中，一层图像大小为512×512像素，一套数据有几百到千层不等的图像。在医学图像处理软件上医师可以自由观看横断面、冠状面和矢状面图像，进行窗宽窗位调整、图像锐化等操作以方便观察。考虑到在CT图像中，相邻像素具有很强的组织隶属联系，比如如果一颗像素是空气/骨质/血管，那么与它相邻的像素也很可能是空气/骨质/血管，我们认为像素可以进一步合并为更大的同质单元，使医师在诊断时可以直接对组织层面上的体块进行操作，而不仅仅面对一堆单个像素。

在医学图像中，相同组织在空间位置和灰度值上都具有连续性，因此可以使用图像处理中的分割算法对医学图像进行分割，得到同质区域。目前，将图像划分为同质区域的主流算法有以下几种：

分裂与合并(Split and Merge)：对整个图像递归进行分块，对每一块进行判断：当前块不属于同质区域时，对其继续进行分块并对每一块进行判断，否则停止分块，保证得到的每一块都是同质区域。然后考虑相邻块，如果其属于同一组织，则合并成一个较大的块。分裂与合并算法的缺点在于，需要大量的递归分裂与同质区域判断。

均值漂移(Mean Shift)：对图像上每一点执行某种寻找密度中心的算法，算法保证局部同质区域中的所有点最终会被吸引到同一点，收集最终落到同一点的像素标记为同质区域。均值漂移算法的缺点在于，需要对图像中每一点做卷积以判断其密度中心的位置。

分水岭(Watershed)：将图像视为盆地和山脊，从最下方开始往上漫水，当两块盆地中的水连起来时，则在其中增加一座水坝阻止其连通。最终所有的水坝构成图像的一个分割，不同水坝围起的区域成为同质区域。分水岭算法的缺点在于：要求一个不断“漫水”的过程，对于灰度值跨度有4000级、层数有几百层的CT图像序列处理而言计算量过大。

超级像素(SuperPixel)：通过同质区域在图像空间域和灰度值上的邻近性，通过迭代算法将图像划分为许多同质小区域块，称为超级像素。超级像素算法的缺点在于容易过分分割图像，且不能选择感兴趣区域有针对性地进行分割。

综上，目前的算法存在计算量过大，计算结果不准确或者不易修改的缺陷，因此，研究一种快速、准确和易于修改的医学图像分割算法。

发明内容

针对上述现有技术中存在的缺陷或不足，本发明的目的在于，提供一种医学图像的分割方法，该方法把医学图像分割成一系列小块，保证每一小块中的像素属于同一组织区域，使得后续过程能够直接在同质区域级别上操作而不是在单个像素点上操作，大大降低了计算量，执行速度快且结果准确。

为了达到上述目的，本发明采用如下的技术方案予以实现：

一种医学图像的分割算法，该方法包括如下步骤：

步骤一、选择感兴趣区域：

1)设定初始感兴趣区域的灰度值范围为U，初始感兴趣区域是指医学图像中灰度值在范围U中的所有像素点的集合；生成一个空栈S；设定曼哈顿空间距离阈值T₁＝1.5*STEP，T₂＝2.5*STEP，其中，STEP＝round(V^1/3)是一个类块的期望边长的取整，round表示四舍五入取整，V是聚类分割的一个类块中包含的像素点的期望个数；

2)遍历医学图像的所有像素点，如果当前像素点的灰度值属于U，则在栈S中寻找与当前像素点的曼哈顿空间距离d_M小于T₁的点；如果找不到，则将当前像素点放入栈S。如果能找得到或者如果当前像素点的灰度值范围不属于U，则转到下一像素点重复执行步骤2)，直到所有像素点被遍历。

3)遍历栈S中所有像素点，把医学图像中所有的距离栈S中任一像素点的曼哈顿空间距离d_M小于T₂的像素点标记为感兴趣像素点，所有感兴趣像素点组成感兴趣区域ROI。

步骤二、布置初始聚类中心点：

在感兴趣区域中选取初始聚类中心点，保证初始聚类中心点的空间坐标(x，y，z)的每个分量均为STEP的任意整数倍；将初始聚类中心点的灰度值取为其空间坐标(x，y，z)对应的像素点的灰度值，将初始聚类中心点的空间坐标取为对应的x，y，z值。

步骤三、聚类分割：

对感兴趣区域ROI进行niter轮k-means聚类，每一轮聚类过程的具体步骤如下：对ROI中每个像素点的到最近的聚类中心的距离赋值为无穷大；然后遍历所有的聚类中心，对每一个聚类中心，遍历其周围曼哈顿距离d_M小于STEP的所有像素点，如果当前像素点到当前聚类中心的加权距离d小于此像素点存储的到最近的聚类中心的距离，则用当前聚类中心的编号替换该像素点隶属的聚类中心编号，该像素点存储的到最近的聚类中心的距离更新为其到当前聚类中心的加权距离d；

所有的聚类中心遍历结束后，更新每个聚类中心的灰度值为所有隶属于它的像素点的灰度值的平均值，更新每个聚类中心的空间坐标为所有隶属于它的像素的空间坐标的平均值；至此，一轮聚类过程结束；

niter次聚类完成后，隶属于同一聚类中心的像素点为同一类点，总类数即为初始聚类中心数目；

本发明还包括如下其他技术特征：

所述算法还包括步骤四：连通性调整；具体按照如下步骤执行：

记允许类最小体积为包含minV个像素点；为感兴趣区域中每一像素点初始化类标号数组，并统一赋值为未标号；然后对感兴趣区域进行遍历：如果当前像素点的类标号数组的值为未标号，则顺次赋予其一个新类号，并且把与该像素点连通且隶属于相同聚类中心的所有像素点赋予所述的新类号，如果被赋予新类号的类的体积小于minV，则将其中所有像素点的类号置换为与其相邻的任意一类的类号，即将其合并到它周围较大的类中。

所述步骤一中的曼哈顿空间距离是指：记两个像素点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，则X₁和X₂之间的曼哈顿空间距离d_M(X₁,X₂)＝|x₁-x₂|+|y₁-y₂|+|z₁-z₂|。

所述步骤三种的加权距离d是指：记两点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，该两点对应像素点的灰度值分别为I₁和I₂，则X₁和X₂之间的加权距离d为：

$d(X_1,X_2)=\sqrt{\frac{{(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2}{{\mathrm{STEP}}^2}+\frac{{(I_1-I_2)}^2}{{\mathrm{compactness}}^2}}$

式中，compactness是加权系数，取100至400。

所述步骤三中所述的加权系数compactness取400。

所述步骤三中所述的聚类次数niter取5。

与现有的图像分割算法相比较，本发明的方法的有益效果如下：

通过有选择地对医学图像的感兴趣区域进行快速聚类分割，将医学图像划分为一系列同质区域，保证每一个区域中的像素点属于相同的组织，方便后续在组织层面上的操作。同时，本发明的方法根据试验选取了合适的聚类次数，计算量小，执行速度快，分割准确。本发明特别适用于计算机控制X射线断层扫描装置或医学放射性影像处理系统中对组织进行定位、测量、识别或分类的功能中。

附图说明

图1是本发明的医学图像的分割方法的流程图。

图2是选择感兴趣区域的流程图。

图3是布置初始聚类中心的流程图。

图4是一轮聚类过程的流程图。

图5是连通性调整的流程图。

图6是本发明的一个实施例的示意图。其中，(a)为待分割的三维CT图像，(b)中白色区域为初始感兴趣区域，(c)为在初始感兴趣区域中均匀布置的稀疏种子点，(d)为聚类分割的结果，(e)为进行连通性调整后的结果。

具体实施方式

如图1所示，本发明的医学图像的分割方法，处理对象是由多层二维图像构成的三维医学图像，具体步骤包括选取感兴趣区域、布置初始聚类中心点、聚类分割和连通性调整。

步骤一、选择感兴趣区域

在实际应用中，感兴趣的图像区域往往不是整幅医学图像，因此，为了减少算法计算时间，首先要选择医学图像的感兴趣区域(ROI)，后续的聚类分割只在感兴趣区域上进行。当ROI不是全部医学图像时，这一步骤去除了对不感兴趣区域的分割操作，将显著减少算法的运行时间。如图2所示，步骤一具体包括如下步骤：

1)设定初始感兴趣区域的灰度值范围为U，初始感兴趣区域是指医学图像中灰度值在范围U中的所有像素点的集合；生成一个空栈S；设定曼哈顿空间距离阈值T₁＝1.5*STEP，T₂＝2.5*STEP，T₁和T₂的选择根据经验确定。其中，STEP＝round(V^1/3)是聚类分割的一个类块的期望边长的取整，round表示四舍五入取整。V是聚类分割的一个类块中包含像素点的期望个数(也即一个类快的期望体积)，其取值随医学图像的采样率和分割用途而定，一般取1000-4000。

步骤一的目的是将医学图像中初始感兴趣区域中的所有点，以及与初始感兴趣区域在空间上邻近的点一起标记为感兴趣区域参与后续分割。注意，灰度值不在初始感兴趣区域的灰度值范围U中，但是在空间上与初始感兴趣区域邻近的像素点也将成为感兴趣区域中的点。U可以是一个区间，也可以是多个区间的并集，最简单的情形是一个区间：U＝[a,b]，即感兴趣的待分割组织的灰度值范围在区间[a,b]中。如果没有人为指定初始感兴趣灰度值范围，则认为U＝[Imin,Imax]，其中Imin和Imax分别代表医学图像灰度值的最小值和最大值，此时初始感兴趣区域是医学图像的全部。

2)遍历医学图像的所有像素点，如果当前像素点的灰度值属于U，说明当前像素点属于初始感兴趣区域，则在栈S中寻找与当前像素点的曼哈顿空间距离d_M小于T₁的点。记两个像素点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，则X₁和X₂之间的曼哈顿空间距离(即Manhattan距离)d_M(X₁,X₂)＝|x₁-x₂|+|y₁-y₂|+|z₁-z₂|；如果找不到，则将当前像素点放入栈S。如果能找得到或者如果当前像素点的灰度值范围不属于U，则转到下一像素点重复执行步骤2)，直到所有像素点被遍历。

遍历结束时，栈S中存储的像素点的灰度值都在初始感兴趣区域中；栈S中任意两个像素点的曼哈顿空间距离均不小于T₁；对于除栈S之外的任意初始感兴趣区域中的像素点，栈S中至少存在一个像素点与它的曼哈顿空间距离小于T₁。说明栈S中存储的像素点是均匀分布并布满初始感兴趣区域的稀疏像素点。

3)遍历栈S中所有像素点，把医学图像中所有的距离栈S中任一像素点的曼哈顿空间距离d_M小于T₂的像素点标记为感兴趣像素点，所有感兴趣像素点组成感兴趣区域(ROI)。感兴趣区域在空间上可以是一个连通区域，也可以是多个不连通区域的并集，取决于初始感兴趣区域在空间上的分布情况。

在步骤一中，遍历图像、在栈S中存储初始感兴趣区域的“骨架”的目的是为了减少计算量；将栈S中所有像素点的T₂邻域中的所有像素点标记为感兴趣区域而不是直接将初始感兴趣区域标记为感兴趣区域参与分割，是为了获得比初始感兴趣区域更宽的区域作为后续分割区域，吸收初始感兴趣区域附近灰度值不在U中的遗漏的感兴趣像素。

步骤二、布置初始聚类中心点

如图3所示，在感兴趣区域(ROI)中选取初始聚类中心点，保证初始聚类中心点的空间坐标(x，y，z)的每个分量均为STEP的任意整数倍。这样就保证初始聚类中心在医学图像中的分布成正方体晶格状，任意相邻的8个初始聚类中心恰好在体积为V的正方体的8个顶点处，使初始聚类中心“均匀分布”并且“占满”ROI。将初始聚类中心点的灰度值取为其空间坐标(x，y，z)对应的像素点的灰度值，将初始聚类中心点的空间坐标取为对应的x，y，z值。

步骤三、聚类分割

对感兴趣区域(ROI)进行k-means聚类，共进行niter轮聚类；该聚类过程的初始聚类中心即为步骤二中得到的在空间上“均匀分布”的初始聚类中心，随着聚类的进行，聚类中心不断变化，ROI中每个像素点隶属的聚类中心编号也不断变化。

如图4所示，每一轮聚类过程的具体步骤如下：对ROI中每个像素点的到最近的聚类中心的距离赋值为无穷大；然后遍历所有的聚类中心，对每一个聚类中心，遍历其周围曼哈顿距离d_M小于STEP的所有像素点，如果当前像素点到当前聚类中心的加权距离d小于此像素点存储的到最近的聚类中心的距离，则用当前聚类中心的编号替换该像素点隶属的聚类中心编号，该像素点存储的到最近的聚类中心的距离更新为到当前聚类中心的加权距离d。

加权距离d的定义：记两点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，该两点对应像素点的灰度值分别为I₁和I₂，则X₁和X₂之间的加权距离d为：

$d(X_1,X_2)=\sqrt{\frac{{(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2}{{\mathrm{STEP}}^2}+\frac{{(I_1-I_2)}^2}{{\mathrm{compactness}}^2}}$

式中，compactness为加权系数，其意义为控制被STEP归一化的空间欧式距离||X₁-X₂||和灰度值距离|I₁-I₂|之间的权衡。compactness越大，最终聚类结果越规则整齐，具有较大灰度差异的区域越可能被分到一类，compactness越小，最终聚类结果越崎岖，每一类可能具有较大空间距离跨度，但是灰度值较统一。理想的聚类结果应在空间上尽量规则的同时又能把灰度值差异大的像素点分到不同类中，经试验，本发明中取compactness在区间[100,400]中较好。

所有的聚类中心遍历结束后，更新每个聚类中心的灰度值为所有隶属于它的像素点的灰度值的平均值，更新每个聚类中心的空间坐标为所有隶属于它的像素的空间坐标的平均值。至此，一轮k-means聚类结束。

niter轮聚类完成后，隶属于同一聚类中心的像素点为同一类点，总类数即为初始聚类中心数目。过多的聚类次数会增加算法运行时间，过少的聚类次数会使不同特征的像素不能充分分开，经试验，本发明中取niter＝5。

步骤四、连通性调整：

由于k-meas聚类可能产生体积过小的类以及不连通的类，因此需要将体积过小的类进行合并，将不连通的类的每个连通分支分别赋予不同的类号，最终得到一系列连通的，无细小零碎杂质的分类结果。具体执行步骤如下：

如图5所示，记允许类最小体积为包含minV个像素点，minV为合并体积阈值；为感兴趣区域(ROI)中每一像素点初始化类标号数组，并统一赋值为“未标号”。

然后对感兴趣区域进行遍历：如果当前像素点的类标号数组的值为未标号，则顺次赋予其一个新类号，并且把与该像素点连通且隶属于相同聚类中心的所有像素点赋予所述的新类号，如果被赋予新类号的类的体积小于minV，则将其中所有像素点的类号置换为与其相邻的任意一类的类号，即将其合并到它周围较大的类中。

其中，连通定义为：一组像素点B中任意两个像素点可由一条由B中像素点组成的路径相连，使路径上任意相邻像素点的x、y、z方向空间坐标的差值都不大于1。在二维图像中，一个非边缘像素点有8个相邻点，在三维图像中，一个非边缘像素点有26个相邻点。

为了说明本发明的有益效果，发明人给出了如下实施例，该实施例遵循本发明的算法的技术方案，本实验处理对象为296层512×512大小的三维CT图像。选择聚类分割后一个类块中包含的像素点的期望个数V＝2000，空间距离与颜色距离权值compactness＝400，聚类次数niter＝5，合并体积阈值minV＝100，初始感兴趣区间U＝[-24,4000]。试验过程中某一层图像处理步骤如图6所示。其中，(a)为待分割的三维CT图像，(b)中白色区域为初始感兴趣区域，(c)为在初始感兴趣区域中均匀布置的稀疏种子点，(d)为聚类分割的结果，(e)为进行连通性调整，重新分配类号，合并小区域后的最终分割结果。算法的运行结果将不同组织分割到不同的小类块中，并且由于算法只对感兴趣区域进行，运行速度比对全部图像使用算法相比有很大提升，5次聚类在第二次就基本收敛。可以看出，本发明能够对医学图像进行快速、准确地分割。

Claims (5)

1.一种医学图像的分割方法，其特征在于，该方法包括如下步骤：

步骤一、选择感兴趣区域：

1)设定初始感兴趣区域的灰度值范围为U，初始感兴趣区域是指医学图像中灰度值在范围U中的所有像素点的集合；生成一个空栈S；设定曼哈顿空间距离阈值T₁＝1.5*STEP，T₂＝2.5*STEP，其中，STEP＝round(V^1/3)是一个类块的期望边长的取整，round表示四舍五入取整，V是聚类分割的一个类块中包含的像素点的期望个数；

2)遍历医学图像的所有像素点，如果当前像素点的灰度值属于U，则在栈S中寻找与当前像素点的曼哈顿空间距离d_M小于T₁的点；如果找不到，则将当前像素点放入栈S；如果能找得到或者如果当前像素点的灰度值范围不属于U，则转到下一像素点重复执行步骤2)，直到所有像素点被遍历；

3)遍历栈S中所有像素点，把医学图像中所有的距离栈S中任一像素点的曼哈顿空间距离d_M小于T₂的像素点标记为感兴趣像素点，所有感兴趣像素点组成感兴趣区域ROI；

步骤二、布置初始聚类中心点：

在感兴趣区域中选取初始聚类中心点，保证初始聚类中心点的空间坐标(x，y，z)的每个分量均为STEP的任意整数倍；将初始聚类中心点的灰度值取为其空间坐标(x，y，z)对应的像素点的灰度值，将初始聚类中心点的空间坐标取为对应的x，y，z值；

步骤三、聚类分割：

对感兴趣区域ROI进行niter轮k-means聚类，每一轮聚类过程的具体步骤如下：对ROI中每个像素点的到最近的聚类中心的距离赋值为无穷大；然后遍历所有的聚类中心，对每一个聚类中心，遍历其周围曼哈顿距离d_M小于STEP的所有像素点，如果当前像素点到当前聚类中心的加权距离d小于此像素点存储的到最近的聚类中心的距离，则用当前聚类中心的编号替换该像素点隶属的聚类中心编号，该像素点存储的到最近的聚类中心的距离更新为其到当前聚类中心的加权距离d；所述加权距离d是指：记两点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，该两点对应像素点的灰度值分别为I₁和I₂，则X₁和X₂之间的加权距离d为：

$d(X_1,X_2)=\sqrt{\frac{{(x_1-x_2)}^2+{(y_1-y_2)}^2+{(z_1-z_2)}^2}{{\mathrm{STEP}}^2}+\frac{{(I_1-I_2)}^2}{{\mathrm{compactness}}^2}}$

式中，compactness是加权系数，取100至400；

所有的聚类中心遍历结束后，更新每个聚类中心的灰度值为所有隶属于它的像素点的灰度值的平均值，更新每个聚类中心的空间坐标为所有隶属于它的像素的空间坐标的平均值；至此，一轮聚类过程结束；

niter次聚类完成后，隶属于同一聚类中心的像素点为同一类点，总类数即为初始聚类中心数目。

2.如权利要求1所述的医学图像的分割方法，其特征在于，还包括步骤四：连通性调整；具体按照如下步骤执行：

记允许类最小体积为包含minV个像素点；为感兴趣区域中每一像素点初始化类标号数组，并统一赋值为未标号；然后对感兴趣区域进行遍历：如果当前像素点的类标号数组的值为未标号，则顺次赋予其一个新类号，并且把与该像素点连通且隶属于相同聚类中心的所有像素点赋予所述的新类号，如果被赋予新类号的类的体积小于minV，则将其合并到它周围较大的类中。

3.如权利要求1所述的医学图像的分割方法，其特征在于，所述步骤一中的曼哈顿空间距离是指：记两个像素点X₁＝(x₁,y₁,z₁)，X₂＝(x₂,y₂,z₂)，则X₁和X₂之间的曼哈顿空间距离d_M(X₁,X₂)＝|x₁-x₂|+|y₁-y₂|+|z₁-z₂|。

4.如权利要求1所述的医学图像的分割方法，其特征在于，所述步骤三中所述的加权系数compactness取400。

5.如权利要求1所述的医学图像的分割方法，其特征在于，所述步骤三中所述的聚类次数niter取5。