基于稀疏表示的医学CT图像肾上腺肿瘤分割方法

摘要

本发明属于图像处理技术领域，具体为基于稀疏表示的医学CT图像中肾上腺肿瘤的分割方法。本发明使用训练所得对边界敏感的过完备字典，将二维CT图像的感兴趣区域分解成图像子块并作稀疏表示；对于图像同质区与非同质区两者稀疏分解的第一个系数的绝对值差异明显，选取合适阈值区分该系数，得到对应的图像边界子块，并获得一个二值图像；利用区域生长法在此二值图像上生长出肿瘤的粗糙轮廓，将其作为水平集分割方法的初始轮廓，通过多次迭代，获取肿瘤最终边界。本发明可以大大提高图像中肾上腺肿瘤图像分割的自动化程度，同时降低了水平集分割方法对初始轮廓的依赖程度，使得分割结果更加准确。

说明书

技术领域

本发明属于图像处理技术领域，具体涉及稀疏表示的医学CT图像肾上腺肿瘤分割方法。

背景技术

肾上腺位于腹膜后，其下外侧则与两侧肾的上内侧紧密贴近。它能够分泌肾上腺素、皮质激素等多种调节人体生理功能的重要激素，是体内重要的内分泌器官。一旦肾上腺的某一部位发生肿瘤，相应部位分泌的激素就会失调，引发一些心血管疾病，严重时会威胁生命。

因为肾上腺肿瘤引发的病症与其它肾上腺疾病(如肾上腺皮质增生)相似，但是两者的治疗方案不同。因此，其诊断需要借助于医学影像技术。CT图像由于其高密度分辨率，在成像诊断中有着广泛的应用。但是由于肿瘤病灶组织与周围组织在CT图像中很接近，使得它们之间的边界难于区分。

水平集的方法在医学CT图像分割的问题上有很好的性能[1]-[3]，它具有灵活的约束条件且对(肿瘤)拓扑结构的适应性强。但是水平集的方法对于初始轮廓的要求很高，往往需要人工勾画，较差的初始轮廓容易导致分割结果进入局部极小值。

针对上述问题，本发明利用稀疏表示以及图像块处理的方法，通过区域生长法自动获取肿瘤初始轮廓，该初始轮廓已十分接近于肿瘤的轮廓，用它作为基于局部信息水平集分割方法的初始轮廓，提高了分割的自动化水平，降低了后续水平集分割的难度，获取的最终分割结果更加准确。

发明内容

本发明的目的是提出一种可以有效提高分割准确率与自动化水平的的医学CT图像肾上腺肿瘤分割方法。

本发明提出的医学CT图像中肾上腺肿瘤的分割方法，具体步骤为：

（1）获取一系列包含肿瘤边界图像子块作为训练样本，用K均值聚类奇异值分解(KSVD)的方法，训练出一个对边界敏感的过完备字典；

（2）利用步骤（1）得到的字典，把CT图像上的感兴趣区域(ROI)的每一个图像子块做稀疏表示，获得第一个稀疏系数的绝对值；选取适当的阈值，区分第一个系数值，把图像子块分成同质区和边界区两类，从而把原图像映射成一个二值图像；

（3）用区域生长的方法，从步骤（2）中的二值图像中心点开始，往外生长得到感兴趣区域内肿瘤的粗糙轮廓；

（4）以步骤（3）得到的粗糙轮廓作为水平集分割方法的初始轮廓，用基于局部信息的水平集分割方法获得最终的精确轮廓。

下面就本发明方法的各个步骤涉及的相关技术细节作进一步的具体描述。

关于步骤1，获取一系列包含肿瘤边界图像子块的做法为：首先，对已知肿瘤边界的训练图像I进行分块，利用                                                的滑动窗口，从左到右，从上到下，遍历整幅图像，依次提取图像子块；把N个包含肿瘤边界的图像块保留下来，作为训练样本，并将其展开为的列向量，M=25。其余图像块舍弃。

感兴趣区域是在图1(a)的红框内，图1(b)为放大后的感兴趣区域。

图2是所选取训练样本的示例，红色轮廓为已知的肿瘤轮廓，蓝色框为所选取的训练样本。

所述用K均值聚类奇异值分解(KSVD)的方法，训练出一个对边界敏感的过完备字典，具体做法如下：

设包含肿瘤边界的图像块集合为，N=900。冗余字典D为L个25维原子的集合(L=100)，经过稀疏分解后的系数为X，则包含肿瘤边界的图像块集合Y的稀疏表示为：

                    (1)

其中与分别为l²-范数与l⁰-范数，前者为向量各个元素平方和的次方，后者为向量中非0元素的个数，为稀疏度。式(1)表示，在约束条件下，从稀疏分解后的系数X中找到一系列合适的列向量，使得最终的X能让取到最小值。

式(1)中的图像块集合Y由人工选定的包含边界的图像块组成，因此如果可以找到合适的字典D，使之能够对集合Y中所有元素均能做到稀疏的表示，那么便能够获得一个如上文所述对边界敏感的过完备字典。显而易见，以此字典对包含边界的图像块做稀疏表示，其系数的稀疏性要强于不包含边界(相对匀质)的图像块做稀疏表示的系数。

首先，每一个包含肿瘤边界的图像块样本的稀疏表示为：

             (2)

其中为以离散余弦变换基元取的初始化字典；系数矩阵 ，其第列表示第个训练样本的稀疏分解系数。

接着固定字典D，采用正交匹配跟踪算法找到最优系数X的近似解；然后，利用奇异值分解(SVD)不断更新原子，求得最终的目标字典；设为的第列，为经过第T次的迭代后X的第行，假定更新第个原子，将矩阵相乘的方式从行乘列求和改写为列乘行求和，则式(2)变为：

       (3)

令记录中使用原子的训练样本，即：

                     (4)

定义矩阵为，在处为1，其余全为0，令，则式(3)改为：

                 (5)

对进行奇异值分解： ，利用矩阵的第一列更新原子；经过多次迭代后，对所有原子进行更新，则得到最终的冗余字典D。

重复上述稀疏分解与字典D的更新过程，直至(2)式收敛，获得对边界敏感的字典D。

关于步骤（2），从图3可以看出，包含边界的图像子块与不包含边界的图像子块在纹理和结构上有明显差别，因此两者的稀疏分解系数有较大的差异。所述把原图像映射成一个二值图像的具体过程为：

首先，对CT图像进行感兴趣区域S的人工提取，设，利用的滑动窗口，从左到右，从上到下遍历整幅图像，依次提取图像子块，并将其展开为的列向量；这样，CT图像分成了Z个图像子块，，Z个图像子块的稀疏表示为：

                (6)

其中，D为训练后的字典；系数矩阵，其第列表示第个图像子块的稀疏分解系数。

由于D是基于肿瘤边界训练得到的字典，对边界敏感。将其对CT图像进行稀疏分解，得到的同质区和边缘区的分解系数有很大的不同，特别是第一个系数值的绝对值，差距很大。故比较系数矩阵中每个分解系数的第一个值的绝对值，取适当阈值，判断该的图像子块是否是包含边界的图像块：若的第一个值的绝对值大于阈值，则说明其为边界子块，用逻辑值1代替原图像子块的中心像素值；若的第一个值的绝对值小于阈值，则说明其为同质区子块，用逻辑值0代替原图像子块的中心像素值。在比较了所有图像子块之后，获得了原图像的二值图像的映射。

图4(a)~(f)显示了图3(a)中，1~6号蓝色框的图像子块(包含边界)的稀疏表示结果。图5(a)~(f)显示了图3(b)中，1~6号绿色框的图像子块(不包含边界)的稀疏表示结果。可以明显看出：包含边界的6块，稀疏表示的第一个系数均大于0.30；不包含边界的6块，稀疏表示的第一个系数均小于0.15。

图6例举了原始感兴趣区域与二值映射后的对应图像。

关于步骤（3），从步骤（2）中获得的二值图像后，可以用区域生长的方法获得肿瘤的粗略轮廓。

在保证肿瘤大致位于感兴趣区域中心的前提下，选取感兴趣区域的中心点为初始生长点，以“十”字形的模板向四周生长。若生长点上下左右四邻域内的点值为0，则生长至此点，并将此点作为下一轮生长的种子点之一；若生长点上下左右四邻域内的点值为1，则不生长至此点，中心处的原种子点则不在下一轮继续生长。过程如图7所示。

由于CT图像中肿瘤与周围组织的CT值十分接近，在步骤2中获得的二值映射有时不完全是闭合的。参见图6(a)与(b)。因此，通过生长的方法获得一个肿瘤的粗略轮廓需要一个终止条件。对于一个类圆形的肿瘤，其感兴趣区域的二值映射也是一个类圆形的二值图像。生长过程从1个点开始，第二轮变为4个点，第三轮变为8个点，随后逐渐增加，到达边界后生长点个数开始逐渐下降。若二值图像是闭合的，生长点个数最终下降为0；若二值图像不是闭合的，则生长点个数会在下降到极小值后有再上升的趋势。如图8所示，(c)为从二值图像(a)的中心点开始生长的每轮生长点个数曲线。本方法利用每轮生长点个数的变化作为生长是否终止的判断，当生长点个数达到第一个极小值时，便停止生长。这里，在寻找极小值之前要将每轮的生长点个数曲线作移动平均的处理，目的是为了平滑掉每轮的生长点个数微小波动产生的极小值。图8中，(d)为(c)平滑后结果。图9显示了有无生长终止条件的生长结果差异。(b)为无生长终止条件，(c)为有生长终止条件。

然后，用梯度算子可以十分简便地获得肿瘤的粗略轮廓。该轮廓十分接近肿瘤的真实边缘。

关于步骤（4），获得了肿瘤的初始轮廓之后，用基于局部信息的水平集方法[5]，利用步骤（3）得到的粗略轮廓作为水平集分割的初始轮廓。由于初始轮廓十分接近肿瘤的真实边缘，所以在水平集迭代中，只需设置适当参数，便得到肿瘤的准确边缘。

假设从步骤（3）中已获得初始闭合轮廓C，再定义图像上的点到初始轮廓的有符号距离为，轮廓内的点与轮廓的距离为负，轮廓外的点与轮廓的距离为正。那么轮廓内的点可以用平滑后的Heaviside函数表示为：

         (7)

其中ε是个很小的常数。

显然，轮廓外的点则可以表示成：。再定义局部区域为：

                                         (8)

其中为轮廓上的点，为除外的其它点。表示两者间欧氏距离。

那么对于轮廓上某一个点，定义局部区域能量为：

  (9)

其中表示图像中，点的亮度，与分别为轮廓内与轮廓外局部区域的亮度均值：

                                           (10)

                                          (11)

这里，局部区域的能量可以表示成：

                     (12)

其中有：     (13)

引入轮廓曲率的惩罚项后，水平集能量函数为：

       (14)

其中，λ为权重系数轮廓曲率的惩罚项权重系数，为梯度算子。

由能量函数得到最终的轮廓迭代方程：

                        (15)

其中表示散度。

由式(15)的迭代方程与初始轮廓，最终可获得分割后的最终轮廓。

本发明方法可以大大提高图像中肾上腺肿瘤图像分割的自动化程度，同时降低了基于局部信息水平集分割方法对初始轮廓的依赖程度，使得分割结果更加准确。

附图说明

图1：(a)肾上腺原醛瘤CT图像，(b)感兴趣区域。

图2：包含边界的训练样本(蓝色框)示例图像。

图3：(a)包含边界的图像子块(蓝色框)，编号为1~6，(b)不包含边界的图像子块(绿色框)，编号为1~6。

图4：(a)~(f)对应于图3(a)中1~6号包含边界的图像子块稀疏表示结果。

图5： (a)~(f)对应于图3(b)中1~6号不包含边界的图像子块稀疏表示结果。

图6：(a)与(c)为原CT图像，红色轮廓为肿瘤区轮廓，(b)与(d)分别为(a)与(c)二值映射后的图像。

图7：“十”字形区域生长方法的示意图。

图8：(a)二值映射图像，(b)不闭合边界的区域生长结果，(c)每轮生长点的个数变化图，(d)平滑后的每轮生长点的个数变化图。

图9：(a)肿瘤图像，红色轮廓为肿瘤区轮廓，(b)无终止条件的区域生长结果，(c) 有终止条件的区域生长结果。

图10：分割结果：(a)原图像，(b)手工勾勒的椭圆初始轮廓，(c) 手工勾勒的椭圆初始轮廓分割结果，(d)医生勾勒的肿瘤轮廓，(e)自动提取的初始轮廓，(f) 自动提取的初始轮廓分割结果。

具体实施方式

对本发明提出的分割方法进行实际CT图像测试。分别采用本方法获取的初始轮廓与手工勾画的椭圆轮廓作为水平集分割方法的初始轮廓，对30幅肿瘤图像进行分割。以医生手工勾勒的肿瘤轮廓为金标准，将两种方法的分割结果分别与手工勾勒的肿瘤轮廓作比较。本方法中，二值映射阈值设为0.30。

分割结果比较参数采用重叠面积(area of overlap，AO)和准确率(accuracy，AC)。

重叠面积AO:

准确率AC：

其中TP为真阳性率，FP为假阳性率:

GS、SE分别表示金标准分割与水平集分割结果的点集合，表示集合内点的个数。

分割结果越准确，AO、AC越接近于100%。

如图10所示，本方法获得的初始轮廓更加接近肿瘤形状，最终的分割结果也更加准确。

30幅肾上腺CT肿瘤图像分割结果比较见表1。无论是重叠面积AO还是准确率AC，本方法都优于手工勾画的椭圆轮廓作为水平集分割方法的初始轮廓的分割结果，同时AO、AC这两项参数的标准差也是本方法低，说明本方法更加稳定可靠。整个分割过程中，除了必要的感兴趣区域ROI的选取，其余步骤均能通过计算机完成，不需要人工参与。

综上所述，本发明可以有效地提高CT图像肿瘤分割的自动化程度，同时降低基于局部信息水平集分割方法对初始轮廓的依赖程度，使得分割结果更加准确。

表1：两种方法分别应用于30幅CT图像后的平均TP、FP、AC以及对应的标准差。

。

参考文献

[1] S. Osher, R. Tsai. Level set methods and their applications in image science [J].  Communications in Contemporary Mathematics, 2003, 1(4): 1-20.

[2] T. Chan, L. Vese. Active contours without edges [J]. IEEE Transaction on Image Processing, 2001, 10(2): 266-277.

[3] B.N. Li, C.K. Chui, S. Chang, S.H. Ong. A new unified level set method for semi-automatic liver tumor segmentation on contrast enhanced CT images [J]. Expert Systems with Applications, 2012, 39(10): 9661-9668.

[4] M. Aharon, M. Elad, A. Bruckstein. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation [J]. IEEE Transaction on Signal Processing, 2006, 54(11): 4311-4322.

[5] S. Lankton, A. Tannenbaum. Localizing region-based active contours [J]. IEEE Transaction on Image Processing, 2008, 17(11): 2029-2039.。