CN105139454A - 一种三维ct图像中肝脏三维感兴趣区域的自动提取方法 - Google Patents

Abstract

一种三维CT图像中肝脏三维感兴趣区域的自动提取方法，利用各向异性扩散过程、梯度幅度计算过程及S形曲线的非线性映射过程对人体腹腔三维CT图像进行滤波，提取出初步的三维ROI。再界定三维ROI区域内坐标值和灰度值离散程度阈值，通过对三维ROI区域内各个体素的坐标值和灰度值进行离散程度判断，去除离散程度高的区域，从而更进一步的获得三维ROI。该方法可以提取出包含肝脏信息的三维图像，大幅度降低图像处理的数据量。

Description

一种三维CT图像中肝脏三维感兴趣区域的自动提取方法

技术领域

本发明属于图像处理技术领域，特别涉及一种核医学图像自动提取感兴趣区域(ROI)方法，是相关梯度滤波和依据图像离散程度去除粘连组织的方法用于核医学图像自动提取ROI的技术。

背景技术

现如今，计算机断层扫描(CT)，磁共振(MR)，正电子放射断层扫描(PET)等医学成像技术已经成为医生诊断和治疗的重要辅助工具，尤其在对肿瘤等恶性疾病的诊断中扮演着重要的角色。其中CT成像由于其本身能精确反映解剖信息，信噪比较高，成像分辨率高等优点，成为在肿瘤等疾病的检查与治疗中最常用的成像方式。越来越多的研究者开始对CT图像的感兴趣区域(ROI)进行提取，用以对疾病分析研究及精确治疗。例如，通过提取的三维ROI自动分割出肝脏、肝肿瘤、和肝血管，以辅助医生诊断。因此，提取三维ROI的计算机辅助诊断具有重大的意义。

三维ROI的提取是基于CT成像的计算机辅助诊断的重要环节，近年来，为了实现对CT图像三维ROI的提取，大量研究者进行了深入的研究，现阶段对于CT图像的ROI提取方面已提出诸多方法，以如下几类为主：阈值法、区域增长法、聚类法、形变模型法、水平集、图谱法、图切割和神经网络等。这些方法主要存在以下缺陷：(1)部分方法需要手动交互过程，如选取初始种子点或初始轮廓。(2)部分方法对噪声或初始化位置的选取十分敏感。(3)不同病患的ROI结构存在显著差异、同一病患不同时间ROI的结构变形和密度改变导致很难使用一种分割方法满足不同需求情况。(4)多数方法算法复杂度高，处理时间长，不适合临床使用。

准确、快速、自动化是人们致力研究图像分割的最终目标。但分割的准确性、有效性以及自动化水平也是相互矛盾的。如果能在保证不损失图像有效信息的情况下明显降低图像处理的数据量，将有助于提升图像分割的效率。例如首先在原始图像中提取出ROI区域，仅对ROI区域的数据进行处理，则可以满足上述要求。此时，如何迅速、准确的定位ROI区域，就显得至关重要。

本发明旨在提出一种对CT图像中三维ROI进行准确、快速、自动化提取方法。本发明首先对图像进行各向异性扩散，使图像的均匀区域变得平滑且图像边缘信息予以保留，之后通过高斯核与图像的卷积做一次平滑处理，在进行微分计算得到图像各像素点的梯度幅值，再将此梯度幅值的大小进行映射，得到ROI所属的梯度范围。至此，通过最大连通域提取即可得到初步的ROI。之后对坐标值离散程度高的区域和灰度值离散程度高的区域周围部分粘连组织进行去除。最后通过二值图像的三维包围盒运算得到最终ROI的提取结果。

本发明提出了基于梯度滤波和依据图像离散度进行CT图像三维ROI的自动提取的方法，大幅度降低图像运算数据量，借此能够指导肝脏的快速、准确且自动化程度高的分割，对于指导肝脏肿瘤分割和肝脏血管分割等方面也有重要的意义。另外，经过相关梯度滤波处理的二值化图像，能够作为图割法和水平集法等图像处理过程的初始数据，以便进行其余后续处理。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提出一种基于梯度滤波和依据图像离散程度进行CT图像三维感兴趣区域(ROI)的自动提取的方法，利用各向异性扩散过程、梯度幅度计算过程及S形曲线的非线性映射过程对人体腹腔三维CT图像进行滤波，提取出初步的三维ROI。再界定三维ROI区域内坐标值和灰度值离散程度阈值，通过对三维ROI区域内各个体素的坐标值和灰度值进行离散程度判断，去除离散程度高的区域，从而更进一步的获得三维ROI。该方法可以提取出包含肝脏信息的三维图像，大幅度降低图像处理的数据量。

本发明提出的基于梯度滤波和依据图像离散度进行CT图像三维兴趣域的自动提取方法，包括以下步骤：

1.1.对于一幅三维人体腹腔CT图像I_o进行向下重采样得到图像I'_o；

1.2.对图像I'_o进行梯度滤波，得到包含肝脏梯度范围的三维ROI图像I_s；

1.3.对图像I_s进行二值处理，提取最大连通区域，得到图像I'_lc；

1.4.界定图像坐标离散程度阈值及灰度值离散程度阈值，去除图像I'_lc超过阈值的区域，得到图像I_gd；

1.5.对图像I_gd再次提取最大连通域，对最大连通域进行向上重采样过程，得到图像I'_r，对图像I'_r进行形态学膨胀变换，得到图像I_d；

1.6.对图像I_d进行三维包围盒运算，从原始图像I_o中剪切出最终包含肝脏信息的ROI区域。

上述步骤1.2具体包括如下步骤：

2.1.对向下重采样后图像I'_o进行各向异性扩散处理，得到图像I_g；

2.2.对图像I_g进行梯度幅度计算，得到包含图像梯度幅度信息的图像I_m；

2.3.对图像I_m进行S形曲线非线性映射，得到包含肝脏梯度范围信息的图像I_s；

上述步骤1.3具体包括如下步骤：

3.1.对图像I_s进行二值处理，得到图像I_b；

3.2.对图像I_b进行CT值范围处理，将CT值低于20且高于350的区域CT值置零，得到图像I'_b；

3.3.对图像I'_b提取最大单连通区域，得到图像I_lc；

3.4.对图像I_lc中，CT值非零区域的CT值设为255，得到图像I'_lc；

上述步骤1.4具体包括如下步骤：

4.1.对图像I'_lc进行z轴坐标运算，得到各体素坐标均值m_z和标准差sd_z，以为z轴坐标值离散程度判断依据，与界定离散程度阈值n进行比较，对离散程度超过阈值的区域进行形态学腐蚀变换。重新提取最大单连通区域，再次对此时z轴坐标进行运算，得到新的离散程度值，再次与阈值n比较，如此循环，直到区域的离散程度小于阈值n为止；

4.2.对图像I_c的灰度值进行运算，得到各体素灰度均值m_g和标准差sd_g，对于图像I_c区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围内的部分不做处理；对于图像I_c区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围外部分的灰度值置零，得到图像I_gd。

附图说明

图1：本发明方法的流程图；

图2：待提取的原始三维CT图像(左图为横断面面，中图为矢状状面，右图为冠状面，以下图例类似)；

图3：梯度各向异性扩散过程效果图像；

图4：带平滑过程的梯度幅度计算过程效果图像；

图5：S形曲线的非线性映射过程效果图像；

图6：首次二值化图像后，提取最大单连通域的效果图；

图7：z轴坐标离散程度处理后与图6过程对比效果图(离散程度未超过阈值)；

图8：z轴坐标离散程度处理后与图6过程对比效果图(离散程度超过阈值)；

图9：对经坐标离散处理后图像进行灰度离散度处理过程的前后对比效果图(离散程度未超过阈值)；

图10：对经坐标离散处理后图像进行灰度离散度处理过程的前后对比效果图(离散程度超过阈值)；

图11：再次提取最大连通域，对其进行向上重采样过程效果图；

图12：膨胀后提取的肝脏效果图；

图13：包含肝脏信息的ROI效果图。

具体实施方式

结合附图和实际例子对提取过程进行具体的描述，下分步进行介绍。(由于测试数据选用增强型CT，处理过程中肝脏区域内部有受血管中造影剂影响的孔洞，若使用普通CT数据的处理结果会更佳)。

1.对原始三维CT图像I_o进行向下重采样，x、y、z轴的采样尺度为s_x，s_y，s_z。本方法的优选值分别为s_x＝2，s_y＝2，s_z＝1，得到图像I'_o。

2.使用梯度滤波器对图像进行梯度滤波处理，得到包含肝脏梯度范围的三维ROI图像I_s。具体实施如下，

1)对图像I'_o进行各向异性扩散处理，得到图像I_g，如图3。本方法采用梯度各向异性扩散处理，执行N维标量值图像的经典Perona-Malik各向异性方程，旨在得到均匀密度组织平滑，但边界突出的图像结果。在实现过程中使用曲率各向异性扩散实现这一目的。

Perona-Malik模型为：

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=div\left(c\right(\left|\▿u\right|)\▿u)\\ u(x,y,0)=u_o(x,y)\end{array}\right.$

其中u(x,y,0)＝u₀(x,y)表示输入图像；而G(σ)表示标准差为σ的高斯运算；c(·)＞0称为扩散系数，是单调递减且与图像梯度成反比的函数；

运算符div执行散度运算；是梯度算子，表示进行空间各向全微分运算；

本方法中此时扩散系数为： $C(\▿u)=e^{-{\left(\frac{\left|\▿u\right|}{k}\right)}^2}$

其中包含两个参数：电导系数k和时间参数t，k控制边缘对比处理的灵敏性。本方法中k值优选为3.0，时间参数t值优选为0.0625。

求解Perona-Malik模型的偏微分方程采用数值运算的偏微分方程的迭代求解算法，本发明中优选的迭代次数为4。

2)对图像I_g进行梯度幅度计算，得到包含图像梯度幅度信息的图像I_m，如图4。这一过程计算了每个像素图像的梯度。这个计算过程等同于首先通过将图像和一个离散化的高斯核卷积来平滑图像，再通过微分计算得到像素点梯度值。

设离散化的高斯卷积核H为(2k+1)*(2k+1)维，则高斯卷积核元素的计算方式为：

$H_{i,j}=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-\frac{{(i-k-1)}^2+{(j-k-1)}^2}{2{\mathrm{\σ}}^2}}$

高斯窗口宽度由σ的值决定，其取值范围为1至1.5，此处优选值为1.3。

3)对图像I_m进行S形曲线的非线性映射处理，得到图像I_s，如图5。其目的在于得到肝脏感兴趣区域的梯度范围。为了实现S形曲线的非线性映射的过程，需对四个参数进行设置。接下来的方程表达了S形曲线的非线性映射的灰度变换：

$I^{\′}=(Max-Min)\·\frac{1}{1+e^{-\frac{I-\mathrm{\β}}{\mathrm{\α}}}}+Min$

在上面的方程中，I是输入像素的灰度，I'是输出像素的灰，Min、Max是输出图像的最小值和最大值，α定义了输入灰度范围的宽度，β定义了围绕在范围中心的灰度。本过程中选择以下优选值α＝-0.5，β＝3.0，Min＝0.0，Max＝1.0。

3.对图像I_s进行二值处理，提取最大连通区域，得到图像I'_lc，如图6。具体实施过程如下：

1)对图像I_s进行二值处理，得到图像I_b。若输入图像是归一化图像，阈值下限优选值为0.1，阈值上限优选值为1.0；若输入图像为非归一化图像，阈值下限优选值为25，阈值上限优选值为255。将小于低阈值的图像灰度设为0，将大于高阈值的图像灰度设为255。

2)对图像I_b进行CT值范围处理，将CT值低于20且高于350的区域CT值置零，得到图像I_b'。因为肝脏的CT值优选范围在20至350之间，此操作可去除骨骼及部分脂肪等组织信息或空气等环境信息的影响。

3)对图像I'_b提取最大单连通区域，得到图像I_lc，由于肝脏是人体腹腔内最大的连通组织，由此获得初始的包含肝脏信息的ROI。

4.界定图像坐标离散程度阈值及灰度值离散程度阈值，去除图像I_lc超过阈值的区域，得到图像I'_lc。具体实施过程如下，

1)对图像I'_lc进行z轴坐标运算，得到各体素坐标均值m_z和标准差sd_z，界定z轴坐标值离散程度判断阈值n，如果则不进行操作，如图7；如果则对图像I'_lc进行一次形态学腐蚀变换，得到图像I_be，重新提取最大单连通区域，再次对此时最大单连通区域z轴坐标进行运算，得到坐标均值m'_z和标准差sd'_z，将与n比较，如此循环，直到满足得到图像I_cd，如图8；本方法中，n的优选值为0.4。此步骤将位置不在肝脏区域范围内的脏器(如心脏)予以滤除；

2)对图像I_c的灰度值进行运算，得到各体素灰度均值m_g和标准差sd_g，如果区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围内，则不进行操作，如图9；如果区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围外，则对其灰度值置零，得到图像I_gd，如图10。此步骤将含肝脏信息ROI附近的粘连组织或部分具有灰度差异的器官组织(如肾脏)予以滤除；

5.对图像I_gd再次提取最大连通域，对最大连通域进行向上重采样过程，得到图像I'_r，对图像I'_r进行形态学膨胀变换，得到图像I_d。具体实施过程如下：

1)对图像I_gd再次提取最大连通域；

2)对最大连通域进行向上重采样过程，x、y、z轴的采样尺度为1/s_x，1/s_y，1/s_z；得到图像I_r。

3)图像I_r中，对灰度非零区域的灰度值设定为255。得到图像I'_r，如图11。

4)对图像I'_r进行形态学膨胀变换，得到图像I_d，如图12。本方法进行了三次膨胀变换，变换半径优选值为1，借此膨胀肝脏区域大小，弥补孔洞。

6.针对二值图像I_d进行三维包围盒运算，设定包围盒由六个方向逼近三维ROI，于原始图像I_o中裁剪下最终包含肝脏信息的ROI区域，如图13。图13中含有肝脏信息的三维ROI的数据量为303×273×64，较原始图像I_o数据量有了明显减少(I_o数据量为512×512×79)，其余测试数据中还有明显数据量减少的例子，如原始图像数据量为512×512×129，提取包含肝脏的三维ROI后的数据量减少至217×291×41。

Claims (2)

1.一种三维CT图像中肝脏三维感兴趣区域的自动提取方法，该方法利用各向异性扩散过程、梯度幅度计算过程及S形曲线的非线性映射过程对人体腹腔三维CT图像进行滤波，提取出初步的三维ROI；再界定三维ROI区域内坐标值和灰度值离散程度阈值，通过对三维ROI区域内各个体素的坐标值和灰度值进行离散程度判断，去除离散程度高的区域，从而更进一步的获得三维ROI；该方法可以提取出包含肝脏信息的三维图像，大幅度降低图像处理的数据量；

其特征在于：该方法包括以下步骤，

1.1.对于一幅三维人体腹腔CT图像I_o进行向下重采样得到图像I'_o；

1.2.对图像I'_o进行梯度滤波，得到包含肝脏梯度范围的三维ROI图像I_s；

1.3.对图像I_s进行二值处理，提取最大连通区域，得到图像I′_lc；

1.4.界定图像坐标离散程度阈值及灰度值离散程度阈值，去除图像I′_lc超过阈值的区域，得到图像I_gd；

1.5.对图像I_gd再次提取最大连通域，对最大连通域进行向上重采样过程，得到图像I'_r，对图像I'_r进行形态学膨胀变换，得到图像I_d；

1.6.对图像I_d进行三维包围盒运算，从原始图像I_o中剪切出最终包含肝脏信息的ROI区域；

上述步骤1.2具体包括如下步骤：

2.1.对向下重采样后图像I'_o进行各向异性扩散处理，得到图像I_g；

2.2.对图像I_g进行梯度幅度计算，得到包含图像梯度幅度信息的图像I_m；

2.3.对图像I_m进行S形曲线非线性映射，得到包含肝脏梯度范围信息的图像I_s；

上述步骤1.3具体包括如下步骤：

3.1.对图像I_s进行二值处理，得到图像I_b；

3.2.对图像I_b进行CT值范围处理，将CT值低于20且高于350的区域CT值置零，得到图像I′_b；

3.3.对图像I′_b提取最大单连通区域，得到图像I_lc；

3.4.对图像I_lc中，CT值非零区域的CT值设为255，得到图像I′_lc；

上述步骤1.4具体包括如下步骤：

4.1.对图像I′_lc进行z轴坐标运算，得到各体素坐标均值m_z和标准差sd_z，以为z轴坐标值离散程度判断依据，与界定离散程度阈值n进行比较，对离散程度超过阈值的区域进行形态学腐蚀变换；重新提取最大单连通区域，再次对此时z轴坐标进行运算，得到新的离散程度值，再次与阈值n比较，如此循环，直到区域的离散程度小于阈值n为止；

4.2.对图像I_c的灰度值进行运算，得到各体素灰度均值m_g和标准差sd_g，对于图像I_c区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围内的部分不做处理；对于图像I_c区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围外部分的灰度值置零，得到图像I_gd。

2.根据权利要求1所述的一种三维CT图像中肝脏三维感兴趣区域的自动提取方法，其特征在于：

其实施过程如下，S1对原始三维CT图像I_o进行向下重采样，x、y、z轴的采样尺度为s_x，s_y，s_z；本方法的优选值分别为s_x＝2，s_y＝2，s_z＝1，得到图像I'_o；

S2使用梯度滤波器对图像进行梯度滤波处理，得到包含肝脏梯度范围的三维ROI图像I_s；具体实施如下，

1)对图像I'_o进行各向异性扩散处理，得到图像I_g；本方法采用梯度各向异性扩散处理，执行N维标量值图像的经典Perona-Malik各向异性方程，旨在得到均匀密度组织平滑，但边界突出的图像结果；在实现过程中使用曲率各向异性扩散实现这一目的；

Perona-Malik模型为：

$\left\{\begin{array}{c}\frac{\∂u}{\∂t}=div\left(c\right(\left|\▿u\right|)\▿u)\\ u(x,y,0)=u_o(x,y)\end{array}\right.$

其中u(x,y,0)＝u₀(x,y)表示输入图像；而G(σ)表示标准差为σ的高斯运算；c(·)＞0称为扩散系数，是单调递减且与图像梯度成反比的函数；

运算符div执行散度运算；是梯度算子，表示进行空间各向全微分运算；

本方法中此时扩散系数为：

$C(\▿u)=e^{-{\left(\frac{\left|\▿u\right|}{k}\right)}^2}$

其中包含两个参数：电导系数k和时间参数t，k控制边缘对比处理的灵敏性；本方法中k值优选为3.0，时间参数t值优选为0.0625；

求解Perona-Malik模型的偏微分方程采用数值运算的偏微分方程的迭代求解算法，优选的迭代次数为4；

2)对图像I_g进行梯度幅度计算，得到包含图像梯度幅度信息的图像I_m；这一过程计算了每个像素图像的梯度；这个计算过程等同于首先通过将图像和一个离散化的高斯核卷积来平滑图像，再通过微分计算得到像素点梯度值；

设离散化的高斯卷积核H为(2k+1)*(2k+1)维，则高斯卷积核元素的计算方式为：

$H_{i,j}=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-\frac{{(i-k-1)}^2+{(j-k-1)}^2}{2{\mathrm{\σ}}^2}}$

高斯窗口宽度由σ的值决定，其取值范围为1至1.5，此处优选值为1.3；

3)对图像I_m进行S形曲线的非线性映射处理，得到图像I_s；其目的在于得到肝脏感兴趣区域的梯度范围；为了实现S形曲线的非线性映射的过程，需对四个参数进行设置；接下来的方程表达了S形曲线的非线性映射的灰度变换：

$I^{\′}=(Max-Min)\·\frac{1}{1+e^{-\frac{I-\mathrm{\β}}{\mathrm{\α}}}}+Min$

在上面的方程中，I是输入像素的灰度，I'是输出像素的灰，Min、Max是输出图像的最小值和最大值，α定义了输入灰度范围的宽度，β定义了围绕在范围中心的灰度；本过程中选择以下优选值α＝-0.5，β＝3.0，Min＝0.0，Max＝1.0；

S3对图像I_s进行二值处理，提取最大连通区域，得到图像I′_lc；具体实施过程如下：

1)对图像I_s进行二值处理，得到图像I_b；若输入图像是归一化图像，阈值下限优选值为0.1，阈值上限优选值为1.0；若输入图像为非归一化图像，阈值下限优选值为25，阈值上限优选值为255；将小于低阈值的图像灰度设为0，将大于高阈值的图像灰度设为255；

2)对图像I_b进行CT值范围处理，将CT值低于20且高于350的区域CT值置零，得到图像I′_b；因为肝脏的CT值优选范围在20至350之间，此操作可去除骨骼及部分脂肪等组织信息或空气环境信息的影响；

3)对图像I′_b提取最大单连通区域，得到图像I_lc，由于肝脏是人体腹腔内最大的连通组织，由此获得初始的包含肝脏信息的ROI；

S4界定图像坐标离散程度阈值及灰度值离散程度阈值，去除图像I_lc超过阈值的区域，得到图像I′_lc；具体实施过程如下，

1)对图像I′_lc进行z轴坐标运算，得到各体素坐标均值m_z和标准差sd_z，界定z轴坐标值离散程度判断阈值n，如果则不进行操作；如果则对图像I′_lc进行一次形态学腐蚀变换，得到图像I_be，重新提取最大单连通区域，再次对此时最大单连通区域z轴坐标进行运算，得到坐标均值m'_z和标准差sd'_z，将与n比较，如此循环，直到满足得到图像I_cd；本方法中，n的优选值为0.4；此步骤将位置不在肝脏区域范围内的脏器予以滤除；

2)对图像I_c的灰度值进行运算，得到各体素灰度均值m_g和标准差sd_g，如果区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围内，则不进行操作；如果区域内体素灰度值在[m_g-sd_g,m_g+sd_g]范围外，则对其灰度值置零，得到图像I_gd；此步骤将含肝脏信息ROI附近的粘连组织或部分具有灰度差异的器官组织予以滤除；

S5对图像I_gd再次提取最大连通域，对最大连通域进行向上重采样过程，得到图像I'_r，对图像I'_r进行形态学膨胀变换，得到图像I_d；具体实施过程如下，

1)对图像I_gd再次提取最大连通域；

2)对最大连通域进行向上重采样过程，x、y、z轴的采样尺度为1/s_x，1/s_y，1/s_z；得到图像I_r；

3)图像I_r中，对灰度非零区域的灰度值设定为255；得到图像I'_r；

4)对图像I'_r进行形态学膨胀变换，得到图像I_d；本方法进行了三次膨胀变换，变换半径优选值为1，借此膨胀肝脏区域大小，弥补孔洞；

S6针对二值图像I_d进行三维包围盒运算，设定包围盒由六个方向逼近三维ROI，于原始图像I_o中裁剪下最终包含肝脏信息的ROI区域。