CN102509286B - 一种医学图像目标区域勾画方法 - Google Patents

Abstract

一种医学图像目标区域勾画方法，包括以下步骤：（1）预处理；（2）在参考时相上勾画出目标轮廓，并根据该目标轮廓选取若干轮廓控制区域；（3）分别提取所述若干轮廓控制区域的图像纹理特征；（4）根据图像纹理特征，跟踪参考时相上的每个轮廓控制区域在待处理时相上的对应位置；（5）在待处理时相上跟踪得到的每相邻两个轮廓控制区域的中心之间进行插值处理，从而完成待处理时相上目标区域的自动勾画。本发明勾画方法能够根据参考时相上的目标轮廓和轮廓控制区域，快速、准确地自动勾画出其它时相上的目标区域，能够大大减轻工作人员在图像分割中的工作量。

Description

一种医学图像目标区域勾画方法

技术领域

本发明涉及图像处理，特别是利用计算机对4D图像中的目标区域进行勾画的方法。

背景技术

目标靶区勾画是放疗技术中极为重要的步骤，该步的准确度直接决定了放疗计划的精确程度。目前临床中都是由医师手动对目标进行勾画来制定放疗计划。而4D-CT图像实际是由一组各个不同时相的常规CT图像组成的，图像数目众多，通常约为1000～2000张。巨大的图像数目极大地加重了医师勾画的工作量，这种重复勾画不仅浪费了宝贵的人力和精力，而且容易引起医师疲劳而产生人为失误等。

实现目标靶区边界的自动勾画，最直观的思路应该是直接对每个呼吸时相的图像采用图像分割技术来自动提取出目标靶区。区域生长法是目前图像分割中比较常用的方法，需要首先选取出一个或多个像素点作为种子点，然后逐一对种子像素周围邻域内的像素根据事先确定的生长准则进行判定，若满足，则将该像素合并到种子像素集中，将这些新像素当作新的种子像素继续进行上面操作过程，直到再没有满足条件的像素点可被包括进来。然而种子点的选取一般都是人工手动完成的，而对4D-CT的大量图像来说工作量仍然很大。

由于呼吸运动属于半周期运动，具有一定的规律性，而由此引起的目标运动和变形也应服从一定的统计规律。因此，也有学者提出在一些先验知识的基础上建立一个概率模型来描述目标靶区的运动和变形，从而实现对其轮廓的自动跟踪。常用的概率模型有贝叶斯滤波模型（bayesian filter）、卡尔曼滤波模型（kalman filter）、粒子滤波模型（particle filter）等。然而建立一个概率模型需要大量的训练数据，而呼吸运动较大的个体差异性使得模型建立难度进一步加大。并且，当出现不规则呼吸运动时，很难再用经典的统计模型来描述人体内的器官运动和变形。

近年来，基于图像配准的轮廓映射方法成为目标跟踪和分割问题的比较流行的解决方法，其思路是首先在参考图像中手动勾画目标轮廓，然后将轮廓点映射到其他待处理图像上的对应位置。该类方法的关键在于能够得到准确的点点间的映射关系，该映射关系通常都是采用变形配准来实现的。其缺点在于直接对整幅图像进行变形配准，因此计算量很大，而轮廓映射的准确度也会受到离目标轮廓较远区域的影响。

发明内容

本发明的目的是提供一种医学图像目标区域勾画方法，以快速、准确地从医学图像中勾画出目标区域，大大减轻工作人员在图像分割中的工作量。

本发明的总体思路是：手动勾画出参考时相上的目标轮廓，而在其他时相上以该目标轮廓为参考，通过寻找与轮廓控制区域相似的图像小块来实现目标靶区的自动跟踪，最后通过在跟踪得到的每两个相邻轮廓控制区域中心之间进行非线性插值实现对其他时相靶区边界的自动勾画。进一步，可在轮廓跟踪过程中，通过引入一个约束项来保证跟踪到的边界的完整性和连贯性。

本发明的技术方案如下：一种医学图像目标区域勾画方法，包括以下步骤：（1）对图像进行去噪及增强边界特征的预处理；（2）在参考时相上勾画出目标轮廓，并根据该目标轮廓选取若干轮廓控制区域；（3）分别提取所述若干轮廓控制区域的图像纹理特征；（4）根据图像纹理特征，在待处理时相上搜索，跟踪参考时相上的每个轮廓控制区域在待处理时相上的对应位置；（5）在待处理时相上跟踪得到的每相邻两个轮廓控制区域的中心之间进行插值处理，从而完成待处理时相上目标区域的自动勾画。

一种较佳方案中，步骤（1）采用各向异性滤波算法对图像进行预处理。

一种较佳方案中，步骤（2）通过以下方法选取轮廓控制区域：把目标轮廓作为一个多边形，分别以多边形的每个顶点为中心构建轮廓控制区域。更好的一种方案是，当多边形的两个相邻顶点之间的距离大于设定的阈值A时，在该两个相邻顶点之间的目标轮廓上再选取至少一个像素点作为中心构建轮廓控制区域。

一种较佳方案中，步骤（3）采用sobel算子提取轮廓控制区域的图像纹理特征。

一种较佳方案中，步骤（4）根据医学图像对应的器官运动幅度的先验知识，限定待处理时相上每个轮廓控制区域的搜索范围，该搜索范围的中心位置与参考时相上的相应轮廓控制区域的中心位置对应。

步骤（4）中，可以通过以下方案跟踪参考时相上的每个轮廓控制区域在待处理时相上的对应位置：通过公式                                                

计算区域特征向量之间的欧式距离来衡量图像的相似性，其中，Gr _i,k 表示参考时相上第i个轮廓控制区域特征向量的第k个元素，G _j,i,k 表示待处理时相上第i个轮廓控制区域的对应搜索范围内的第j个被搜索区域的第k个元素，M表示特征向量的长度，S _i,j 越小，参考时相上的对应轮廓控制区域与待处理时相上的被搜索区域之间的相似性越大，将相似性最大的被搜索区域作为参考时相上的对应轮廓控制区域在待处理时相上的对应位置。

为了保证跟踪到的边界的完整性和连贯性，本发明进一步可在轮廓跟踪过程中引入一个约束项，具体方法是：步骤（4）中，通过公式计算区域特征向量之间的欧式距离来衡量图像的相似性，其中，Gr _i,k 表示参考时相上第i个轮廓控制区域特征向量的第k个元素，G _j,i,k 表示待处理时相上第i个轮廓控制区域的对应搜索范围内的第j个被搜索区域的第k个元素，M表示特征向量的长度，S _i,j 越小，参考时相上的对应轮廓控制区域与待处理时相上的被搜索区域之间的相似性越大，

同时通过公式

计算相邻轮廓控制区域相关性的变化程度，其中，Cr _i 表示参考时相上第i个轮廓控制区域与待处理时相上对应的相邻轮廓控制区域之间的相关系数，Ct _i,j 表示待处理时相上搜索范围内第i个轮廓控制区域的第j个被搜索区域与参考时相上对应相邻轮廓控制区域之间的相关系数，

    然后通过公式

确定参考时相上的对应轮廓控制区域在待处理时相上的对应位置，其中μ为相似性项S _i,j 与约束项C _i,j 之间的相对权重。

本发明勾画方法能够根据参考时相上的目标轮廓和轮廓控制区域，快速、准确地自动勾画出其它时相上的目标区域，从而能够大大减轻工作人员在图像分割中的工作量。尤其适用于动态器官4D-CT图像的目标区域勾画。

附图说明

图1为本医学图像目标区域勾画方法的框图。

图2为参考时相上轮廓控制区域选取示意图。

图3为实验中CT30时相勾画结果。

图4为实验中CT10、CT30、CT70时相勾画结果。

图5为实验中CT10、CT30时相勾画结果的量化评价结果。

具体实施方式

下面结合附图对本发明做进一步说明。

参照图1，本医学图像目标区域勾画方法包括图像预处理、轮廓控制区域选取、纹理特征提取、轮廓控制区域跟踪和目标轮廓勾画五个步骤，下面分别详细说明。

1. 图像预处理：为了避免图像噪声对后续轮廓跟踪的干扰，首先对图像进行去噪预处理，在去除图像噪声的同时保存并增强图像特征。本实施例采用各向异性扩散技术，将图像去噪过程看成是一个热扩散方程，含有图像噪声的初始图像作为该扩散过程的初始状态，而去噪后的图像则是扩散过程中某个时间点对应的状态。通常选择图像梯度值的倒数作为这一过程的扩散系数，来保证边界像素点上相对较小的扩散以保存边界信息，而非边界像素点上较大扩散以去除图像噪声。

2. 轮廓控制区域选取：图像预处理之后，手动对参考时相上的目标靶区进行勾画。为了避免由于不同勾画者单次操作引起的随机误差，实验中分别请三位医师对同一图像进行两次勾画，取其均值作为最终的参考轮廓。由于呼吸引起的器官运动主要表现为头尾方向上的运动，前后左右方向上的运动相对较小，因此本实施例仅在与冠状面平行的平面上来执行算法。将轮廓近似看成一个多边形，选取该多边形的顶点作为轮廓控制区域的中心。如果相邻两个顶点之间的距离大于设定的阈值A，则需要在两点之间再多选取n个均匀分布的轮廓像素点作为控制区域中心。n的大小取决于这两个顶点之间的距离与阈值A的比值。图2为轮廓控制区域选取的示意图，图中将肺作为目标靶区，其轮廓已被手动勾画出来，见图2中的曲线21，沿曲线21设置的若干小方框22则表示选取的若干轮廓控制区域。可以看出，在轮廓变化较为复杂的区域，选取的控制区域较多较紧密，以便包含足够的反映轮廓变化的信息；而在轮廓较为平坦的部分，选取的控制区域较少较稀疏，从而减少算法的运算量。

阈值A越小，选取的轮廓控制区域越多，跟踪得到的轮廓越平滑，但运算量越大；相反，阈值A越大，轮廓控制区域越少，得到的轮廓越近似多边形，运算量越小。因此，需要合理设置阈值A的大小来达到算法准确度和效率之间很好的平衡。轮廓控制区域的大小对算法的有效性也同样重要，这些区域应该足够大使得能够包含足够的信息来表征该轮廓点，但如果过大则不仅会加重算法的计算负担，而且会引入轮廓距离较远的信息而影响跟踪的准确程度。

    3. 纹理特征提取：该步中，对选取的参考时相上的轮廓控制区域提取图像纹理特征，作为轮廓周围区域的有效表达。考虑到计算复杂性及计算效率，本实施例选用sobel算子来提取图像纹理特征。sobel算子实际上是一个离散差分算子，对每个像素点计算图像灰度值的近似梯度，因此sobel结果反映了该像素点处图像灰度的变化程度，可认为是对该像素邻域的一种描述。最后，对像素各个方向的sobel分量按式（1）结合得到该点处的梯度值，并转化为一维向量作为对应轮廓区域的特征向量。

  ……………………（1）。

4. 轮廓控制区域跟踪：对于参考时相上的每一个轮廓控制区域，通过在其他时相上搜索与其相似的区域来实现轮廓的自动跟踪。由于呼吸运动引起的器官运动主要表现为头尾方向上的运动，左右、前后方向上的器官运动相对较小，因此本实施例中只在与冠状面平行的平面上做二维跟踪取代三维立体跟踪来减少算法的计算量。另外，根据器官运动幅度的先验知识来对其他时相上的搜索区域大小进行限制，从而进一步减少算法的计算量。该搜索区域的中心位置对应于参考时相轮廓控制区域的中心位置。在搜索区域内执行穷尽搜索，计算其中每个小区域与轮廓控制区域之间的相似性。本实施例中，通过计算区域特征向量之间的欧式距离来衡量图像的相似性，其表达式如下：

 ……………………（2）

其中，Gr _i,k 表示参考时相上第i个轮廓控制区域特征向量的第k个元素，G _j,i,k 表示待处理时相上第i个轮廓控制区域的对应搜索范围内的第j个被搜索区域的第k个元素。M表示特征向量的长度。因此，S _i,j 越小，参考时相上的对应轮廓控制区域与待处理时相上的被搜索区域之间的相似性越大。

为了保持轮廓的完整性和连贯性，本实施例在跟踪过程中引入了一个约束项，该约束项考虑了相邻轮廓控制区域之间的相关性，对轮廓可能出现的形状进行约束，其数学表达式如下：

  ……………………（3）

其中，Cr _i 表示参考时相上第i个轮廓控制区域与待处理时相上对应的相邻轮廓控制区域之间的相关系数，Ct _i,j 表示待处理时相上搜索范围内第i个轮廓控制区域的第j个被搜索区域与参考时相上对应相邻轮廓控制区域之间的相关系数。

我们可以这样来理解轮廓跟踪过程。首先，基于相似性的搜索驱使算法将参考轮廓控制区域推动到轮廓运动范围内与其最相似的位置；而当这一相似小块与其周围区域明显不同时，则需要约束项施加一个恢复外力将其重新拉回到轮廓区域的邻域内。因此，待处理时相上的最终轮廓控制区域位置由这两项共同决定，该过程可以用以下的数学表达式来表示：

  ……………………（4）

其中，μ为相似性项与约束项之间的相对权重，决定了跟踪结果分别受图像相似性和轮廓完整性的影响程度。

5. 目标轮廓勾画：跟踪到轮廓控制区域在待处理时相上的对应位置之后，我们将这些区域的中心点认为是轮廓上的像素点。通过在每相邻两点之间进行B样条非线性插值得到的光滑曲线，则被认为是待处理时相上的目标靶区轮廓，从而实现轮廓的自动勾画。

实验：本实验分别用三个人的临床肺部4D-CT图像进行了测试，通过对肺部进行轮廓勾画来初步测试本方法对动态器官的适用性。这三组数据均是通过飞利浦公司Brilliance CT BigBore扫描仪采集的。断层图像格式为DICOM格式，大小为512×512，图像分辨率为0.98×0.98mm²，层厚为5mm。由于肺部周围的脊椎结构随呼吸运动的位移较小，相对肺部运动而言可忽略，而脊椎之间的相似性反而可能会影响轮廓跟踪的准确度，因此，实验选择合适的CT窗来削弱图像中这些脊椎结构的纹理信息。说明：由于图中涉及参考时相上的手动轮廓线直接对应过来的轮廓线、重新手动勾画的轮廓线、以及用本发明方法自动勾画的轮廓线，这些轮廓线有多段重叠部分、并且多次交错，用单一黑白线条结合标号无法清楚地区别，因此，用不同颜色来表示不同的轮廓线。

将呼吸周期分为八个呼吸时相，分别对应4D-CT的CT0～CT70图像。其中CT0表示最大呼气末时相，CT30为最大吸气末时相，因此CT0～CT30为吸气过程，CT40～CT70表示呼气过程。实验中，选择CT0作为参考时相，对该时相的所有断层图像中的肺部轮廓进行手动勾画。当参考时相所有断层都勾画完之后，我们仅在与冠状面平行的平面上来显示肺部轮廓，将其近似为多边形并选择多边形的顶点作为轮廓控制区域的中心。实验中，我们设置距离阈值A为20mm，当相邻顶点之间的距离大于该阈值A时，则在其之间轮廓上再适当选择多个均匀分布的轮廓点作为轮廓控制区域中心。轮廓控制区域的大小设为20mm。文献中报道肺部最大运动位于肺下叶区域，运动幅度为12±2mm，因此我们设置搜索范围为以轮廓点为中心的边长为60mm的方形窗，以保证所有呼吸时相时的轮廓控制区域都始终在搜索范围内。

由于最大吸气末时相与参考时相CT0之间的器官运动及变形最大，因此图3显示了其中一个人最大吸气末CT30时相的肺部勾画结果，分别从横截面、冠状面、矢状面三个视角进行显示。第一列为参考时相即CT0时相，其中曲线31-33表示手动勾画的肺部轮廓。第二列为CT30时相时的轮廓自动勾画结果，其中曲线31’-33’（绿色）为CT0时相手动勾画轮廓在该时相上的对应位置，曲线34-36（红色）为该时相上重新手动勾画的轮廓，曲线37-39（蓝色）则为用本发明方法自动勾画的轮廓。首先，从曲线31’-33’与实际肺部的偏移程度可以看出，肺部运动的幅度较大，且主要发生在肺下部的头尾方向上。而自动勾画的轮廓曲线37-39与手动勾画的轮廓34-36较为相似，大部分都相互重叠，有效地证明了本发明勾画方法的有效性。

另外将CT30时相的勾画结果与CT10、CT70的勾画结果进行对比，来评价本发明勾画方法是否受移位及形变大小的影响。如图4中所示，从左到右依次为CT10、CT30、CT70的勾画结果。从图中可以看出，这三个时相的自动勾画轮廓41-49（蓝色）都与手动勾画轮廓41’-49’（红色）比较吻合。

为了更好的对勾画结果进行评价，这里提出了一个勾画结果的量化标准，即计算目标呼吸时相上自动勾画的轮廓与手动勾画出的轮廓之间的一致性系数（accordance coefficient, AC），其定义如下：

  ……………………（5）

其中，R _m 表示手动勾画轮廓包含的像素集，R _a 表示自动勾画的轮廓所包含的像素集，N(R _m ∩R _a )表示这两个像素集交集中包含的像素个数，N(R _m ∪R _a ) 则表示并集中包含的像素个数。由定义可知，该系数始终小于等于1。如果自动勾画的轮廓与手动勾画的参考轮廓越吻合，其交集和并集所包含的像素个数就越接近，得到的一致性系数就越接近于1；相反，如果自动勾画的轮廓与参考轮廓的差别越大，得到的一致性系数就越小。用该方法对上述实验的勾画结果进行评价，结果如图5所示，可以看出，本发明勾画方法对三个测试者数据的勾画轮廓一致性系数均大于0.9，三组数据CT10时相勾画的一致性系数均值为0.966，CT30时相为0.946，再次说明了本发明勾画方法对动态目标轮廓自动勾画的有效性。而CT30时相的一致性系数小于CT10时相的系数，这是由于CT0与CT30时相之间的变形大于CT0与CT10之间的变形，加大了CT30轮廓区域的图像特征与参考时相上的控制区域特征之间的差别，从而使得算法的准确度有所下降，不过0.946仍然是令人较为满意的结果。

Claims (6)

1.一种医学图像目标区域勾画方法，其特征在于，包括以下步骤：（1）对图像进行去噪及增强边界特征的预处理；（2）在参考时相上勾画出目标轮廓，并根据该目标轮廓选取若干轮廓控制区域；该若干轮廓控制区域沿着所述目标轮廓间隔排列在目标区域的边界处，轮廓控制区域的选取方法为：把所述目标轮廓作为一个多边形，分别以多边形的每个顶点为中心构建轮廓控制区域，当多边形的两个相邻顶点之间的距离大于设定的阈值A时，在该两个相邻顶点之间的目标轮廓上再选取至少一个像素点作为中心构建轮廓控制区域；（3）分别提取所述若干轮廓控制区域的图像纹理特征；（4）根据图像纹理特征，在待处理时相上搜索，跟踪参考时相上的每个轮廓控制区域在待处理时相上的对应位置；（5）在待处理时相上跟踪得到的每相邻两个轮廓控制区域的中心之间进行插值处理，从而完成待处理时相上目标区域的自动勾画。

2.根据权利要求1所述的目标区域勾画方法，其特征在于：步骤（1）中采用各向异性滤波算法对图像进行预处理。

3.根据权利要求1所述的目标区域勾画方法，其特征在于：步骤（3）中采用sobel算子提取轮廓控制区域的图像纹理特征。

4.根据权利要求1所述的目标区域勾画方法，其特征在于：步骤（4）中，根据医学图像对应的器官运动幅度的先验知识，限定待处理时相上每个轮廓控制区域的搜索范围，该搜索范围的中心位置与参考时相上的相应轮廓控制区域的中心位置对应。

5.根据权利要求4所述的目标区域勾画方法，其特征在于：步骤（4）中，通过公式

$S_{i,j}=\sqrt{{\mathrm{\Σ}}_{k=1}^M{(G_{j,i,k}-{\mathrm{Gr}}_{i,k})}^2}$

计算区域特征向量之间的欧式距离来衡量图像的相似性，其中，Gr_i,k表示参考时相上第i个轮廓控制区域特征向量的第k个元素，G_j,i,k表示待处理时相上第i个轮廓控制区域的对应搜索范围内的第j个被搜索区域的第k个元素，M表示特征向量的长度，S_i,j越小，参考时相上的对应轮廓控制区域与待处理时相上的被搜索区域之间的相似性越大，将相似性最大的被搜索区域作为参考时相上的对应轮廓控制区域在待处理时相上的对应位置。

6.根据权利要求4所述的目标区域勾画方法，其特征在于：步骤（4）中，通过公式

$S_{i,j}=\sqrt{{\mathrm{\Σ}}_{k=1}^M{(G_{j,i,k}-{\mathrm{Gr}}_{i,k})}^2}$

计算区域特征向量之间的欧式距离来衡量图像的相似性，其中，Gr_i,k表示参考时相上第i个轮廓控制区域特征向量的第k个元素，G_j,i,k表示待处理时相上第i个轮廓控制区域的对应搜索范围内的第j个被搜索区域的第k个元素，M表示特征向量的长度，S_i,j越小，参考时相上的对应轮廓控制区域与待处理时相上的被搜索区域之间的相似性越大，

同时通过公式

C_i,j=|Ct_i,j-Cr_i|

计算相邻轮廓控制区域相关性的变化程度，其中，Cr_i表示参考时相上第i个轮廓控制区域与待处理时相上对应的相邻轮廓控制区域之间的相关系数，Ct_i,j表示待处理时相上搜索范围内第i个轮廓控制区域的第j个被搜索区域与参考时相上对应相邻轮廓控制区域之间的相关系数，

然后通过公式

j=argmin_jE_i,j=argmin_j(S_i,j+μC_i,j)

确定参考时相上的对应轮廓控制区域在待处理时相上的对应位置，其中μ为相似性项S_i,j与约束项C_i,j之间的相对权重。

Images