CN104616289A - 一种3d ct图像中骨组织的移除方法及系统 - Google Patents

Abstract

本发明公开了一种3D CT图像中骨组织的移除方法及系统，读取3D CT图像数据序列，利用聚类法对其进行划分；对每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量；构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对超体素进行分类并标记类别；依据超体素标记值，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。本发明使用对3D CT数据进行过分割后的超体素作为基本操作单元，无需对不同部位的骨组织模板进行学习训练，同时也避免了由于模板不匹配而引起的骨组织提取不准确问题。

Description

一种3D CT图像中骨组织的移除方法及系统

技术领域

本发明属于目标检索及提取技术领域，具体涉及一种3D CT图像中骨组织的移除方法及系统。本发明用于对CT影像设备所采集的3D CT数据进行自动骨组织的提取及移除。

背景技术

相对于X射线影像，CT(计算机断层成像)能为医生提供关于目标器官完整的、高分辨的3D信息，其在目前已经成为临床观测的一个重要工具。随着数字化临床技术的发展，有效的计算机辅助诊断及决策系统对医生合理及快速诊断及诊疗方案的制定具有重要意义。在诸多应用中，对不同的器官进行提取及分割是对CT图像进行精确解译的一个重要步骤。大量的基于像素、边界及子区域的分割及分类方法可以被用于3D CT图像数据得到可接受的分割结果，如阈值(thresholding)，分水岭(watershed)、区域生长(region growing)、图切(graph-cut)、水平集(level-set)方法等。但通常情况下，3D CT图像数据规模交大，如腹部CT序列通常大于300(含有不少于512*512*300的像素数)，这给数据在计算机中的存储及处理带来一定负担。同时，各种不同的成像因素，如放射剂量，CT图像的重构方法，扫描的层厚等，会对最终的CT图像引入大量的伪影及噪声，这也会使各种解译方法的使用带来一定的难度。

在通过CT数据对心脏和血管的可视化及病理变化研究中，准确及快速的骨组织分割及提取是计算机辅助医学工程中一个非常重要的研究课题。与其他器官的提取及分类不同的是，由于骨组织具有复杂的结构及位置，因此利用简单的机器学习方法对CT数据逐层进行自动骨提取是非常困难的。由于在CT成像中，骨组织比其他的软组织具有较高的CT值，因此阈值方法通常被广泛引用提取骨组织。虽然阈值类方法比较简单且易于实现，但通过一个CT值对骨头组织及软组织进行分类是不合理的，特别是对于血管壁和骨头组织，他们都具有相近的CT值。相对于基于像素的处理方法，基于过分割子区域的处理方法具有一定的优势，特别是在实际的应用系统中，其拥有比基于像素的处理方法更小的存储的计算代价，但该类方法对图像的噪声极为敏感，同时对细小的及具有低对比度的骨结构的检测能力有限。通过将分割问题建模为一个能量最小化过程，形变模型可以精确提取对应目标的边界，特别是对于部分软组织，如肾脏，心脏，肺部等，但对于部分具有非匀质及扩散性边界的骨区域来说(如盆骨及股骨区域)，该类方法对边界的定位及提取仍不能达到较为满意的结果。同时，形变模型对初始化过程较为敏感，较差的初始化方法会使得对该模型的求解陷入局部最优值。 除了上面所提到的各种方法用于骨组织提取之外，统计形状类方法通过大量标记数据对待提取的目标的平均形状进行学习。统计形状类方法能有效抑制成像过程中所产生的伪影及噪声的影响，但提取性能受所学的平均形状的精确性的影响，同时，对骨组织数据来说，形状的标记及训练过程需要极大的人力及时间代价。

发明内容

针对上述现有技术中存在的缺陷与不足，本发明的一个目的在于，提供一种3D CT图像中骨组织的移除方法及系统。

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

一种3D CT图像中骨组织的移除方法，具体包括如下步骤：

步骤1：读取3D CT图像数据序列，利用聚类法对其进行划分，并对划分后得到的类作为超体素进行标号；

步骤2：对步骤1得到的每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量；

步骤3：构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对步骤1得到的超体素进行分类并标记类别；

步骤4：依据步骤3中所得到的超体素标记值y，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。

进一步的，所述步骤1的具体步骤如下：

步骤1.1：读入3D CT图像数据序列I＝{p₁,p₂,…,p_N}，其中，pi表示3D CT图像中的第i个像素；p_i＝(g(p_i),x_i,y_i,z_i)，其中，g(p_i)表示像素p_i的CT值，idx(p_i)＝[x_i,y_i,z_i]表示像素p_i在3D CT图像数据中的位置坐标；N表示3D CT图像中像素的数目；

步骤1.2：输入骨组织CT值阈值区间T＝[t_l,t_h]、期望划分的超体素的数目K及迭代次数ITER；

步骤1.3：在3D CT图像数据序列中，随机选取K个聚类中心C_k＝{g(p_k),idx(p_k)}；其中，g(p_k)表示3D CT图像中的第k个像素p_k的CT值；idx(p_k)表示像素p_k在3D CT图像数据中的位置坐标；

步骤1.4：通过下式对每个聚类中心判断，得到聚类更新标记集U＝{u₁,u₂,…,u_K}；

$u_k=\left\{\begin{array}{c}1,\mathrm{if}{t}_l\≤g{\mathrm{Min}}_k\≤t_h\\ 0,\mathrm{othrewise}\end{array}\right.$

其中，gMin_k表示在第k个聚类中心C_k的S*S*S邻域中像素的CT值的最小值，S取2～3；u_k表示第k个聚类中心C_k的标记，k＝1～K；

步骤1.5：如果第k个聚类中心C_k的标记u_k＝1，则对该聚类中心执行步骤1.6；

步骤1.6：利用下式计算聚类中心C_k的M*M*M的邻域内每个像素对应的梯度值，M＝3～7；并将梯度值最小的像素作为新的聚类中心C_k；

$\▿\left(I\right)={\left|\right|{\▿}_x^2+{\▿}_y^2+{\▿}_z^2\left|\right|}_2^2,s.t.\left\{\begin{array}{c}{\▿}_x=g(x+1,y,z)-g(x-1,y,z)\\ {\▿}_y=g(x,y+1,z)-g(x,y-1,z)\\ {\▿}_z=g(x,y,z+1)-g(x,y,z-1)\end{array}\right.$

其中，g(x,y,z)表示的是在空间位置(x,y,z)处的CT值，同时，▽_x,▽_y及▽_z表示CT值在3D CT图像的中坐标系三个维度上的梯度值；

步骤1.7：对每个聚类中心C_k进行如下处理：在聚类中心C_k的2S*2S*2S邻域内选择与该聚类中心像素距离最近的像素，其中，S＝20～30；所述距离的计算根据下式进行：

$d=|g\left(p_i\right)-g\left(p_j\right)|+\frac{{\left|\right|\mathrm{idx}\left(p_i\right)-\mathrm{idx}\left(p_j\right)\left|\right|}^2}{M}$

其中，‖idx(p_i)-idx(p_j)‖²表示两个像素p_i与p_j在3D CT图像数据中的位置坐标之间的欧氏距离；M表示聚类中心对应的超体素中包含像素的数目；

然后将该距离最近的像素归入当前聚类中心所在的类，同时将该像素从之前类中移除，最后将聚类中心C_k使用其所属超体素的质心点替换；

步骤1.8：将步骤1.7执行ITER次，超体素划分过程结束，此时得到划分后的K个超体素，对它们进行顺序标号。

进一步的，所述步骤1.2中，所述骨组织CT值阈值T＝[960,1400]，迭代次数Iter＝3～8。

进一步的，所述步骤2的具体步骤如下：

步骤21，计算得到每个超体素中所有像素的CT值直方图H_G、CT值差值直方图H_D及CT值累积直方图H_C；然后根据CT值直方图H_G求取灰度均值、灰度方差、最大值、最小值、上四分位值、中位值、下四分位值、众值、熵、能量、对比度、二阶距、三阶距、四阶距和平滑性共15个特征值；根据CT值差值直方图H_D求取熵、能量、对比度、二阶距和平滑性共5个特征值；根据CT值累积直方图H_C求取熵、能量、对比度、二阶距和平滑性共5个特征值；其中，熵Et、能量Eg、对比度Ct及k阶中心距M^k的计算公式如下：

$\begin{array}{c}\mathrm{Et}=-{\mathrm{\Σ}}_{i=1}^M{H}^{\left(i\right)}*\log \left(H^{\left(i\right)}\right)\\ \mathrm{Eg}={\mathrm{\Σ}}_{i=1}^M{\left(H^{\left(i\right)}\right)}^2\\ \mathrm{Ct}={\mathrm{\Σ}}_{i=1}^M{i}^2*H^{\left(i\right)}\\ M^k={\mathrm{\Σ}}_{i=1}^M(i-\stackrel{\‾}{i})*H^{\left(i\right)}\end{array}---\left(4\right)$

其中，H表示CT值直方图H_G、CT值差直方图H_D或CT值累积直方图H_C，H⁽ⁱ⁾表示某直方图H中第i个位置对应的纵坐标值，M表示直方图的统计尺度；

同时，计算每种直方图的最小值、上四分位、中值、下四分位数及最大值；

步骤22，对步骤1得到的每个超体素中的每个像素，在像素的3*3*3邻域提取对应的3DSobel、梯度边缘滤波器、高斯平滑滤波器、Laplacian、线性锐化滤波器的响应值以及局部均值、局部中值、局部极大值和局部极小值共9个特征值；其中，局部均值是指上述像素的3*3*3邻域中所有26个像素的CT值的算术平均值；局部中值是指26个像素由小到大排序后位于中部的值，局部极大值和局部极小值同理也是指对该邻域的统计结果。

通过上述步骤得到第i个像素p_i的像素级特征f(p_i)后，通过下式得到超体素在像素级的特征描述：

f(V_s)^j＝avg(f(p_i)^j),1≤j≤N

其中，f(p_i)表示第i个像素的像素级特征；f(V_s)^j表示第s个超体素Vs中所有像素的像素级特征f(p_i)的均值；N表示该超体素中所包含的像素的数目；

步骤23，两个相邻的超体素间的邻域关系通过该两个超体素的CT值均值、CT值均值比、中心距离、两个超体素的CT值直方图H_G相关度、CT值差直方图H_D相关度及CT值累积直方图H_C相关度共6个特征表示，两个超体素的某种直方图的相关度采用如下式所示的卡方距离进行描述：

${\mathrm{\χ}}^2(h_1,h_2)=\frac{1}{2}{\mathrm{\Σ}}_{i=1}^M\frac{{(h_1^{\left(i\right)}-h_2^{\left(i\right)})}^2}{h_1^{\left(i\right)}+h_2^{\left(i\right)}}$

其中，及分别表示直方图h₁及h₂的第i个元素；M表示直方图的统计尺度。

经过上述步骤，得到对每个超体素进行描述的特征向量。

进一步的，所述步骤2的具体步骤如下：

步骤31，分类器训练：从步骤1得到的超体素中随机选取超体素构造训练集(x_i,y_i)，1<i<TrNum，其中，x_i表示训练集中第i个超体素对应的特征向量，yi表示超体素的类别标记，yi＝1表示超体素为骨组织，yi＝0表示超体素为非骨组织；TrNum表示训练集中超体素数目；使用构造的训练集训练分类器；

步骤32，对于步骤1得到的所有超体素中任一超体素，将步骤2提取得到的该超体素对应的特征向量x输入训练后的分类器，则得到该特征向量对应的标记值y，如果y＝1，则表示该超体素为骨组织，y＝0为非骨组织。

本发明的另一个目的在于，提供一种3D CT数据中骨组织的提取及移除系统，所述系统 包括如下模块：

图像数据读取模块，用以读取3D CT图像数据序列，利用聚类法对其进行划分，并对划分后得到的类作为超体素进行标号；

超体素特征描述模块，用以对步骤1得到的每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量；

分类器训练及超体素分类模块，用以构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对步骤1得到的超体素进行分类并标记类别；

骨组织移除模块，用以依据步骤3中所得到的超体素标记值y，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。

进一步的，所述图像数据读取模块具体实现如下流程的功能：

步骤1.1：读入3D CT图像数据序列I＝{p₁,p₂,…,p_N}，其中，pi表示3D CT图像中的第i个像素；p_i＝(g(p_i),x_i,y_i,z_i)，其中，g(p_i)表示像素p_i的CT值，idx(p_i)＝[x_i,y_i,z_i]表示像素p_i在3D CT图像数据中的位置坐标；N表示3D CT图像中像素的数目；

步骤1.2：输入骨组织CT值阈值区间T＝[t_l,t_h]、期望划分的超体素的数目K及迭代次数ITER；

步骤1.3：在3D CT图像数据序列中，随机选取K个聚类中心C_k＝{g(p_k),idx(p_k)}；其中，g(p_k)表示3D CT图像中的第k个像素p_k的CT值；idx(p_k)表示像素p_k在3D CT图像数据中的位置坐标；

步骤1.4：通过下式对每个聚类中心判断，得到聚类更新标记集U＝{u₁,u₂,…,u_K}；

$u_k=\left\{\begin{array}{c}1,\mathrm{if}{t}_l\&le;g{\mathrm{Min}}_k\&le;t_h\\ 0,\mathrm{othrewise}\end{array}\right.$

其中，gMin_k表示在第k个聚类中心C_k的S*S*S邻域中像素的CT值的最小值，S取2～3；u_k表示第k个聚类中心C_k的标记，k＝1～K；

步骤1.5：如果第k个聚类中心C_k的标记u_k＝1，则对该聚类中心执行步骤1.6；

步骤1.6：利用下式计算聚类中心C_k的M*M*M的邻域内每个像素对应的梯度值，M＝3～7；并将梯度值最小的像素作为新的聚类中心C_k；

$\&dtri;\left(I\right)={\left|\right|{\&dtri;}_x^2+{\&dtri;}_y^2+{\&dtri;}_z^2\left|\right|}_2^2,s.t.\left\{\begin{array}{c}{\&dtri;}_x=g(x+1,y,z)-g(x-1,y,z)\\ {\&dtri;}_y=g(x,y+1,z)-g(x,y-1,z)\\ {\&dtri;}_z=g(x,y,z+1)-g(x,y,z-1)\end{array}\right.$

其中，g(x,y,z)表示的是在空间位置(x,y,z)处的CT值，同时，▽_x,▽_y及▽_z表示CT值在3D CT图像的中坐标系三个维度上的梯度值；

步骤1.7：对每个聚类中心C_k进行如下处理：在聚类中心C_k的2S*2S*2S邻域内选择与 该聚类中心像素距离最近的像素，其中，S＝20～30；所述距离的计算根据下式进行：

$d=|g\left(p_i\right)-g\left(p_j\right)|+\frac{{\left|\right|\mathrm{idx}\left(p_i\right)-\mathrm{idx}\left(p_j\right)\left|\right|}^2}{M}$

其中，‖idx(p_i)-idx(p_j)‖²表示两个像素p_i与p_j在3D CT图像数据中的位置坐标之间的欧氏距离；M表示聚类中心对应的超体素中包含像素的数目；

然后将该距离最近的像素归入当前聚类中心所在的类，同时将该像素从之前类中移除，最后将聚类中心C_k使用其所属超体素的质心点替换；

步骤1.8：将步骤1.7执行ITER次，超体素划分过程结束，此时得到划分后的K个超体素，对它们进行顺序标号。

进一步的，所述步骤1.2中，所述骨组织CT值阈值T＝[960,1400]，迭代次数Iter＝3～8。

进一步的，所述超体素特征描述模块具体实现如下模块的功能：

步骤21，计算得到每个超体素中所有像素的CT值直方图H_G、CT值差值直方图H_D及CT值累积直方图H_C；然后根据CT值直方图H_G求取灰度均值、灰度方差、最大值、最小值、上四分位值、中位值、下四分位值、众值、熵、能量、对比度、二阶距、三阶距、四阶距和平滑性共15个特征值；根据CT值差值直方图H_D求取熵、能量、对比度、二阶距和平滑性共5个特征值；根据CT值累积直方图H_C求取熵、能量、对比度、二阶距和平滑性共5个特征值；其中，熵Et、能量Eg、对比度Ct及k阶中心距M^k的计算公式如下：

$\begin{array}{c}\mathrm{Et}=-{\mathrm{\&Sigma;}}_{i=1}^M{H}^{\left(i\right)}*\log \left(H^{\left(i\right)}\right)\\ \mathrm{Eg}={\mathrm{\&Sigma;}}_{i=1}^M{\left(H^{\left(i\right)}\right)}^2\\ \mathrm{Ct}={\mathrm{\&Sigma;}}_{i=1}^M{i}^2*H^{\left(i\right)}\\ M^k={\mathrm{\&Sigma;}}_{i=1}^M(i-\stackrel{\&OverBar;}{i})*H^{\left(i\right)}\end{array}---\left(4\right)$

其中，H表示CT值直方图H_G、CT值差直方图H_D或CT值累积直方图H_C，H⁽ⁱ⁾表示某直方图H中第i个位置对应的纵坐标值，M表示直方图的统计尺度；

同时，计算每种直方图的最小值、上四分位、中值、下四分位数及最大值；

步骤22，对步骤1得到的每个超体素中的每个像素，在像素的3*3*3邻域提取对应的3DSobel、梯度边缘滤波器、高斯平滑滤波器、Laplacian、线性锐化滤波器的响应值以及局部均值、局部中值、局部极大值和局部极小值共9个特征值；其中，局部均值是指上述像素的3*3*3邻域中所有26个像素的CT值的算术平均值；局部中值是指26个像素由小到大排序后位于中部的值，局部极大值和局部极小值同理也是指对该邻域的统计结果。

通过上述步骤得到第i个像素p_i的像素级特征f(p_i)后，通过下式得到超体素在像素级的特征描述：

f(V_s)^j＝avg(f(p_i)^j),1≤j≤N

其中，f(p_i)表示第i个像素的像素级特征；f(V_s)^j表示第s个超体素Vs中所有像素的像素级特征f(p_i)的均值；N表示该超体素中所包含的像素的数目；

步骤23，两个相邻的超体素间的邻域关系通过该两个超体素的CT值均值、CT值均值比、中心距离、两个超体素的CT值直方图H_G相关度、CT值差直方图H_D相关度及CT值累积直方图H_C相关度共6个特征表示，两个超体素的某种直方图的相关度采用如下式所示的卡方距离进行描述：

${\mathrm{\&chi;}}^2(h_1,h_2)=\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^M\frac{{(h_1^{\left(i\right)}-h_2^{\left(i\right)})}^2}{h_1^{\left(i\right)}+h_2^{\left(i\right)}}$

其中，及分别表示直方图h₁及h₂的第i个元素；M表示直方图的统计尺度。

经过上述步骤，得到对每个超体素进行描述的特征向量。

进一步的，所述分类器训练及超体素分类模块具体实现如下流程的功能：

步骤31，分类器训练：从步骤1得到的超体素中随机选取超体素构造训练集(x_i,y_i)，1<i<TrNum，其中，x_i表示训练集中第i个超体素对应的特征向量，yi表示超体素的类别标记，yi＝1表示超体素为骨组织，yi＝0表示超体素为非骨组织；TrNum表示训练集中超体素数目；使用构造的训练集训练分类器；

步骤32，对于步骤1得到的所有超体素中任一超体素，将步骤2提取得到的该超体素对应的特征向量x输入训练后的分类器，则得到该特征向量对应的标记值y，如果y＝1，则表示该超体素为骨组织，y＝0为非骨组织。

与传统基于模板匹配方法相比，由于本发明使用对3D CT数据进行过分割后的超体素作为基本操作单元，无需对不同部位的骨组织模板进行学习训练，同时也避免了由于模板不匹配而引起的骨组织提取不准确问题；其次，本发明相对于基于像素的骨组织提取及移除方法，其计算及存储需求将极大降低，同时由于每个超像素的处理过程大部分上是相互独立的，其可以很容易被扩展到GPU并行平台上进行分布式并行处理以进一步提高性能。最后，由于本发明所使用超体素作为基本操作单元，因此所提取的视觉特征描述会比简单的基于像素的处理方法的描述能力更强。

附图说明

图1是本发明的3D CT图像中骨组织的移除方法的流程图。

图2是3D CT图像中骨组织的移除方法的结构示意图。

图3是从部分标记3D CT数据中获取骨区域/非骨区域及所有超体素的CT值分布直方图。

图4是像素级特征提取示例。其中，图4(a)是对像素x0所定义的一个3*3*3邻域；图4(b)是某一296层CT数据中的第51层CT图像数据；图4(c)--图4(e)是第51层CT数据在3D Sobel滤波器在3个不同的方向上的响应；图4(f)--图4(h)是3D线性锐化，Laplacian及Gaussian滤波器在第51层CT数据上的结果。

图5是训练病例Pat_1的超像素标记结果。其中，图5(a)原始CT数据，图5(b)是从图5(a)中移除骨组织后的剩余部分的CT数据。

图6是对CT数据通过步骤1的tSLIC分割得到的超体素结果。其中，图6(a)是病例Pat_1CT数据中的207层数据；图6(b)是超体素是没有进行任何阈值预处理的SLIC分割结果；图6(c)-图6(e)是阈值区间分别设置为960,1024及1200时的超体素分割结果。在图6(c)-图6(e)中，空白区域是没有任何超体素被划分。

图7是使用本发明的方法SLIC及tSLIC对不同的3D CT数据序列进行划分性能比较，其中，图7(a)显示的是划分超体素的数目，图7(b)显示的是划分超体素所使用时间。

图8是本发明实施案例1，对训练集中的病例3D CT数据Pat_2,Pat_3,Pat_4的骨提取及移除结果。其中，图8(a)—图8(c)是原始3D CT数据，图8(d)—图8(f)是通过本发明所提取的骨数据图像，图8(g)—图8(i)是从原始数据中移除骨组织后的数据图像。

图9是本发明实施案例2：对部分测试3D CT数据进行骨提取及移除结果。其中，第一行的四个图像是原始3D CT数据图像，第2行是从第一行中的图像中提取的骨组织，第3行是进行骨组织移除之后的3D CT数据。图9中每一列表示一个病例。

具体实施方式

相关技术术语定义：

1、超体素：超体素是对3D CT图像数据像素集的一个划分，3D CT图像I上所定义的超体素满足下式：

同时，位于同一个超体素中的像素点必须满足空间上的邻域关系及在CT测量值的近邻关系。其中，S表示的是从3D CT图像I中划分得到的超体素的数目。

在将3D CT图像I进行超体素划分后，骨组织提取及移除任务可视为一个二分类问题：

$L:I\&RightArrow;T=\left\{\mathrm{0,1}\right\},\&ForAll;i,L\left(V_i\right)\&Element;T.---\left(7\right)$

其中，1表示骨组织，0表示非骨组织。

因此，骨组织的提取及移除任务可抽象为寻找一个合理的分类器L，分类器L用于对从3D CT图像I中所提取的超体素进行准确快速的标记。

2、随机森林：在机器学习中，随机森林是一个包含多个决策树的分类器，其对样本类别的预测结果是由个别树输出的类别的众数而定的。在随机森林中，对每一棵树的叶子节点都通过骨组织及非骨组织的后验概率进行标记。同时，随机森林中的每一个中间节点都包含一个将视觉特征空间的最佳分割。在将超体素V的视觉特征向量用x表示之后，骨组织提取及移除问题可通过随机森林描述为一个二分类问题，即：L(x)＝0表示超体素被标识为非骨头组织，L(x)＝1则表示超体素被标记为骨组织，其中L表示的是由随机森林所构成的分类器。

本发明的主要思路是：充分利用超体素在数据处理中的优势，首先通过基于阈值，对3DCT数据序列进行快速超体素划分；在超体素划分之后，通过分析骨组织与其他软组织的CT值不同表现，从3个方面，即像素级，超体素级及超体素邻域关系对超体素进行特征描述；最后将骨组织提取及移除任务视为分类问题后，通过随机森林对超体素特征进行分类标记以进行骨提取及移除。

如图1所示，本发明的3D CT图像中骨组织的移除方法，具体包括如下步骤：

步骤1：读取3D CT图像数据序列，利用聚类法对其进行划分，并对划分后的类作为超体素进行标号。具体操作如下：

步骤1.1：读入3D CT图像数据序列I＝{p₁,p₂,…,p_N}，其中，pi表示3D CT图像中的第i个像素；p_i＝(g(p_i),x_i,y_i,z_i)，其中，g(p_i)表示像素p_i的CT值，idx(p_i)＝[x_i,y_i,z_i]表示像素p_i在3D CT图像数据中的位置坐标；N表示3D CT图像中像素的数目；

发明人通过对大量3D CT图像中骨组织及非骨组织的CT值进行统计，统计结果如图3所示。从中可以观测到骨组织的CT值分布在特定的区间内，因此在超体素的划分过程中，仅仅关注位于该特定区间的CT图像数据，则能够有效降低超体素数目及划分时间代价。

步骤1.2：输入骨组织CT值阈值区间T＝[t_l,t_h]、期望划分的超体素的数目K及迭代次数ITER；本发明中，T＝[960,1400]，迭代次数Iter＝3～8，K是根据3D CT图像数据的规模确定的，每个超体素包含10*10*10个像素；

步骤1.3：在3D CT图像数据序列中，随机选取K个聚类中心C_k＝{g(p_k),idx(p_k)}；其中，g(p_k)表示3D CT图像中的第k个像素p_k的CT值；idx(p_k)表示像素p_k在3D CT图像数据中的位置坐标；

步骤1.4：通过下式对每个聚类中心判断，得到聚类更新标记集U＝{u₁,u₂,…,u_K}；

$u_k=\left\{\begin{array}{c}1,\mathrm{if}{t}_l\&le;g{\mathrm{Min}}_k\&le;t_h\\ 0,\mathrm{othrewise}\end{array}\right.$

其中，gMin_k表示在第k个聚类中心C_k的S*S*S邻域中像素的CT值的最小值；(S取2～3)u_k表示第k个聚类中心C_k的标记，k＝1～K；

步骤1.5：如果第k个聚类中心C_k的标记u_k＝1，则对该聚类中心执行步骤1.6；

步骤1.6：利用下式计算聚类中心C_k的M*M*M的邻域内每个像素对应的梯度值，M＝3～7；并将梯度值最小的像素作为新的聚类中心C_k；

$\&dtri;\left(I\right)={\left|\right|{\&dtri;}_x^2+{\&dtri;}_y^2+{\&dtri;}_z^2\left|\right|}_2^2,s.t.\left\{\begin{array}{c}{\&dtri;}_x=g(x+1,y,z)-g(x-1,y,z)\\ {\&dtri;}_y=g(x,y+1,z)-g(x,y-1,z)\\ {\&dtri;}_z=g(x,y,z+1)-g(x,y,z-1)\end{array}\right.$

其中，g(x,y,z)表示的是在空间位置(x,y,z)处的CT值，同时，▽_x,▽_y及▽_z表示CT值在3D CT图像的中坐标系三个维度上的梯度值；

步骤1.7：对每个聚类中心C_k进行如下处理：在聚类中心C_k的2S*2S*2S邻域内选择与该聚类中心像素距离最近的像素，其中，S＝20～30；所述距离的计算根据下式进行：

$d=|g\left(p_i\right)-g\left(p_j\right)|+\frac{{\left|\right|\mathrm{idx}\left(p_i\right)-\mathrm{idx}\left(p_j\right)\left|\right|}^2}{M}$

其中，‖idx(p_i)-idx(p_j)‖²表示两个像素p_i与p_j在3D CT图像数据中的位置坐标之间的欧氏距离；M表示聚类中心对应的超体素中包含像素的数目。

然后将该距离最近的像素归入当前聚类中心所在的类(即所在的超体素)，同时将该像素从之前类中移除，最后将聚类中心C_k使用其所属超体素的质心点替换；

步骤1.8：将步骤1.7执行ITER次，超体素划分过程结束，此时得到划分后的K个超体素，对它们进行顺序标号；

在上述步骤1中，已将3D CT图像中所有的像素点都与其最近的聚类中心进行关联之后，原始聚类中心则被更新为属于当前类的像素的位置均值。在进行若干次迭代之后，对在空间上不相邻，但被标记为同类的像素重新设置聚类标记。最后，记录最终的聚类中心及标记为对应超体素的位置及超体素标记。

步骤2：对步骤1得到的每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量。具体操作如下：

步骤21，计算得到每个超体素中所有像素的CT值直方图H_G、CT值差值直方图H_D及CT值累积直方图H_C；然后根据CT值直方图H_G求取灰度均值、灰度方差、最大值、最小值、上四分位值、中位值、下四分位值、众值、熵、能量、对比度、二阶距、三阶距、四阶距和平滑性共15个特征值；根据CT值差值直方图H_D求取熵、能量、对比度、二阶距和平滑性 共5个特征值；根据CT值累积直方图H_C求取熵、能量、对比度、二阶距和平滑性共5个特征值；其中，熵Et、能量Eg、对比度Ct及k阶中心距M^k的计算公式如下：

$\begin{array}{c}\mathrm{Et}=-{\mathrm{\&Sigma;}}_{i=1}^M{H}^{\left(i\right)}*\log \left(H^{\left(i\right)}\right)\\ \mathrm{Eg}={\mathrm{\&Sigma;}}_{i=1}^M{\left(H^{\left(i\right)}\right)}^2\\ \mathrm{Ct}={\mathrm{\&Sigma;}}_{i=1}^M{i}^2*H^{\left(i\right)}\\ M^k={\mathrm{\&Sigma;}}_{i=1}^M(i-\stackrel{\&OverBar;}{i})*H^{\left(i\right)}\end{array}---\left(4\right)$

其中，H表示CT值直方图H_G、CT值差直方图H_D或CT值累积直方图H_C，H⁽ⁱ⁾表示某直方图H中第i个位置对应的纵坐标值，M表示直方图的统计尺度(本实施例中为0-4095)；

同时，计算每种直方图的最小值、上四分位、中值、下四分位数及最大值；

步骤22，对步骤1得到的每个超体素中的每个像素，在像素的3*3*3邻域提取对应的3DSobel、梯度边缘滤波器、高斯平滑滤波器、Laplacian、线性锐化滤波器的响应值以及局部均值、局部中值、局部极大值和局部极小值(部分特征滤波器响应如图4所示)共9个特征值；其中，局部均值是指上述像素的3*3*3邻域中所有26个像素的CT值的算术平均值；局部中值是指26个像素由小到大排序后位于中部的值，局部极大值和局部极小值同理也是指对该邻域的统计结果。

通过上述步骤得到第i个像素p_i的像素级特征f(p_i)后，通过下式得到超体素在像素级的特征描述：

f(V_s)^j＝avg(f(p_i)^j),1≤j≤N

其中，f(p_i)表示第i个像素的像素级特征(即3D Sobel梯度边缘滤波器等)；f(V_s)^j表示第s个超体素Vs中所有像素的像素级特征f(p_i)的均值；N表示该超体素中所包含的像素的数目；

步骤23，两个相邻的超体素间的邻域关系通过该两个超体素的CT值均值、CT值均值比、中心距离、两个超体素的CT值直方图H_G相关度、CT值差直方图H_D相关度及CT值累积直方图H_C相关度共6个特征表示，两个超体素的某种直方图的相关度采用如下式所示的卡方距离进行描述：

${\mathrm{\&chi;}}^2(h_1,h_2)=\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^M\frac{{(h_1^{\left(i\right)}-h_2^{\left(i\right)})}^2}{h_1^{\left(i\right)}+h_2^{\left(i\right)}}$

其中，及分别表示直方图h₁及h₂的第i个元素；M表示直方图的统计尺度。

最终，本实施例使用一个长度为40的特征向量对每个超体素进行描述，特征向量x中的每个分量所描述的特征值如表1所示。

表1对超体素进行描述的特征向量

步骤3：构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对步骤1得到的K个超体素进行分类并标记类别；本发明中，分类器采用随机森林。具体步骤如下：

步骤31，分类器训练：构造训练集(x_i,y_i)，1<i<TrNum，其中，x_i表示训练集中第i个超体素对应的特征向量，yi表示超体素的类别标记，yi＝1表示超体素为骨组织，yi＝0表示超体素为非骨组织；TrNum表示训练集中超体素数目，其值根据需要选取，值越大，分类器的分类结果越精确，选取时综合考虑结果的精确性以及训练时间。构造训练集的超体素是从步骤1得到的标号后的超体素中随机选取的。

本实施例中，用构造的训练集训练一个含有100棵树，每棵树深度为10的随机森林。该实施例给出了如表2所示的4个病例的CT数据，例如，对于表2中的病例Pat_1，其通过步骤1得到的超体素共有169026个，在步骤3中，从步骤1得到的所有超体素中，随机选取了60414个超体素构成训练集，其中已知5594个超体素为骨组织，将它们标记为1，有54820个超体素为非骨组织，将它们标记为0；对病例Pat_1的标记数据如图3所示。

表2标记的骨区域及非骨区域超体素数目

由表1可看到，四个病例共有445223个超体素，其中被标记为骨组织的超体素共有7103个，约占总超体素的1.60％；被标记为非骨组织的超体素有62043个，约占总体超体素数目的13.94％。可见，非骨组织和骨组织所占比例区别较大，因此，由于该分布不平衡，相差10倍左右，因此采用随机森林进行分类能够得到更加准确的分类结果。

步骤32，对于步骤1得到的所有超体素中任一超体素，将步骤2提取得到的该超体素对应的特征向量x输入训练后的随机森林，则得到该特征向量对应的标记值y，如果y＝1，则表示该超体素为骨组织，y＝0为非骨组织。

步骤4：依据步骤3中所得到的超体素标记值y，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。

为了验证本发明的有效性，发明人对3D CT图像骨提取及移除系统进行测试。测试平台为：winXP SP3，Intel(R)Dual-Core CPU E66003.06GHZ,3.0G内存。测试用3D CT序列信息如表3所示。

表3实例测试用3D CT数据序列信息

实施例1：

本实施例对测试病例Pat_4,Pat_5,Pat_2对超体素分割过程进行测试，对每个病例进行超体素分割结果如图7所示。在超体素分割之后，部分超体素划分结果如图6所示。由图7及图6可见，随着阈值T＝[t_l,t_h]的提高，超体素数量急剧减少，同时，超体素的分割耗时也越来越少。

实施例2：

本实施例对7组病例进行骨提取及移除测试。在该测试实例中，病例Pat_2，Pat_3及Pat_4中的部分超体素被标记用于随机森林分类器训练，具体超体素标记信息如表1所示。骨提取及去除结果如图8所示。对病例Pat_5，Pat_6，Pat_7及Pat_8的骨组织提取及移除结果如图9所示。由图8、图9可见，本发明的系统无需任何对骨组织形状的先验知识即能得到较满意的提取及移除结果。

Claims (10)

1.一种3D CT图像中骨组织的移除方法，其特征在于，具体包括如下步骤：

步骤1：读取3D CT图像数据序列，利用聚类法对其进行划分，并对划分后得到的类作为超体素进行标号；

步骤2：对步骤1得到的每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量；

步骤3：构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对步骤1得到的超体素进行分类并标记类别；

步骤4：依据步骤3中所得到的超体素标记值y，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。

2.如权利要求1所述的3D CT图像中骨组织的移除方法，其特征在于，所述步骤1的具体步骤如下：

步骤1.1：读入3D CT图像数据序列I＝{p₁,p₂,…,p_N}，其中，pi表示3D CT图像中的第i个像素；p_i＝(g(p_i),x_i,y_i,z_i)，其中，g(p_i)表示像素p_i的CT值，idx(p_i)＝[x_i,y_i,z_i]表示像素p_i在3D CT图像数据中的位置坐标；N表示3D CT图像中像素的数目；

步骤1.2：输入骨组织CT值阈值区间T＝[t_l,t_h]、期望划分的超体素的数目K及迭代次数ITER；

步骤1.3：在3D CT图像数据序列中，随机选取K个聚类中心C_k＝{g(p_k),idx(p_k)}；其中，g(p_k)表示3D CT图像中的第k个像素p_k的CT值；idx(p_k)表示像素p_k在3D CT图像数据中的位置坐标；

步骤1.4：通过下式对每个聚类中心判断，得到聚类更新标记集U＝{u₁,u₂,…,u_K}；

$u_k=\left\{\begin{array}{cc}1,& \mathrm{if}{t}_l\&le;{\mathrm{gMin}}_k\&le;t_h\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中，gMin_k表示在第k个聚类中心C_k的S*S*S邻域中像素的CT值的最小值，S取2～3；u_k表示第k个聚类中心C_k的标记，k＝1～K；

步骤1.5：如果第k个聚类中心C_k的标记u_k＝1，则对该聚类中心执行步骤1.6；

步骤1.6：利用下式计算聚类中心C_k的M*M*M的邻域内每个像素对应的梯度值，M＝3～7；并将梯度值最小的像素作为新的聚类中心C_k；

$\&dtri;\left(I\right)={\left|\right|{\&dtri;}_x^2+{\&dtri;}_y^2+{\&dtri;}_z^2\left|\right|}_2^2,s.t.\left\{\begin{array}{c}{\&dtri;}_x=g(x+1,y,z)-g(x-1,y,z)\\ {\&dtri;}_y=g(x,y+1,z)-g(x,y-1,z)\\ {\&dtri;}_z=g(x,y,z+1)-g(x,y,z-1)\end{array}\right.$

其中，g(x,y,z)表示的是在空间位置(x,y,z)处的CT值，同时，及表示CT值在3D CT图像的中坐标系三个维度上的梯度值；

步骤1.7：对每个聚类中心C_k进行如下处理：在聚类中心C_k的2S*2S*2S邻域内选择与该聚类中心像素距离最近的像素，其中，S＝20～30；所述距离的计算根据下式进行：

$d=|g\left(p_i\right)-g\left(p_j\right)|+\frac{{\left|\right|\mathrm{idx}\left(p_i\right)-\mathrm{idx}\left(p_j\right)\left|\right|}^2}{M}$

其中，||idx(p_i)-idx(p_j)||²表示两个像素p_i与p_j在3D CT图像数据中的位置坐标之间的欧氏距离；M表示聚类中心对应的超体素中包含像素的数目；

然后将该距离最近的像素归入当前聚类中心所在的类，同时将该像素从之前类中移除，最后将聚类中心C_k使用其所属超体素的质心点替换；

步骤1.8：将步骤1.7执行ITER次，超体素划分过程结束，此时得到划分后的K个超体素，对它们进行顺序标号。

3.如权利要求2所述的3D CT图像中骨组织的移除方法，其特征在于，所述步骤1.2中，所述骨组织CT值阈值T＝[960,1400]，迭代次数Iter＝3～8。

4.如权利要求1所述的3D CT图像中骨组织的移除方法，其特征在于，所述步骤2的具体步骤如下：

步骤21，计算得到每个超体素中所有像素的CT值直方图H_G、CT值差值直方图H_D及CT值累积直方图H_C；然后根据CT值直方图H_G求取灰度均值、灰度方差、最大值、最小值、上四分位值、中位值、下四分位值、众值、熵、能量、对比度、二阶距、三阶距、四阶距和平滑性共15个特征值；根据CT值差值直方图H_D求取熵、能量、对比度、二阶距和平滑性共5个特征值；根据CT值累积直方图H_C求取熵、能量、对比度、二阶距和平滑性共5个特征值；其中，熵Et、能量Eg、对比度Ct及k阶中心距M^k的计算公式如下：

$\begin{array}{c}\mathrm{Et}=-{\mathrm{\&Sigma;}}_{i=1}^M{H}^{\left(i\right)}*\log \left(H^{\left(i\right)}\right)\\ \mathrm{Eg}={\mathrm{\&Sigma;}}_{i=1}^M{\left(H^{\left(i\right)}\right)}^2\\ \mathrm{Ct}={\mathrm{\&Sigma;}}_{i=1}^M{i}^2*H^{\left(i\right)}\\ M^k={\mathrm{\&Sigma;}}_{i=1}^M(i-\stackrel{\&OverBar;}{i})*H^{\left(i\right)}\end{array}---\left(4\right)$

其中，H表示CT值直方图H_G、CT值差直方图H_D或CT值累积直方图H_C，H⁽ⁱ⁾表示某直方图H中第i个位置对应的纵坐标值，M表示直方图的统计尺度；

同时，计算每种直方图的最小值、上四分位、中值、下四分位数及最大值；

步骤22，对步骤1得到的每个超体素中的每个像素，在像素的3*3*3邻域提取对应的3DSobel、梯度边缘滤波器、高斯平滑滤波器、Laplacian、线性锐化滤波器的响应值以及局部均值、局部中值、局部极大值和局部极小值共9个特征值；其中，局部均值是指上述像素的3*3*3邻域中所有26个像素的CT值的算术平均值；局部中值是指26个像素由小到大排序后位于中部的值，局部极大值和局部极小值同理也是指对该邻域的统计结果。

通过上述步骤得到第i个像素p_i的像素级特征f(p_i)后，通过下式得到超体素在像素级的特征描述：

f(V_s)^j＝avg(f(p_i)^j),1≤j≤N

其中，f(p_i)表示第i个像素的像素级特征；f(V_s)^j表示第s个超体素Vs中所有像素的像素级特征f(p_i)的均值；N表示该超体素中所包含的像素的数目；

步骤23，两个相邻的超体素间的邻域关系通过该两个超体素的CT值均值、CT值均值比、中心距离、两个超体素的CT值直方图H_G相关度、CT值差直方图H_D相关度及CT值累积直方图H_C相关度共6个特征表示，两个超体素的某种直方图的相关度采用如下式所示的卡方距离进行描述：

${\mathrm{\&chi;}}^2(h_1,h_2)=\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^M\frac{{(h_1^{\left(i\right)}-h_2^{\left(i\right)})}^2}{h_1^{\left(i\right)}+h_2^{\left(i\right)}}$

其中，及分别表示直方图h₁及h₂的第i个元素；M表示直方图的统计尺度。

经过上述步骤，得到对每个超体素进行描述的特征向量。

5.如权利要求1所述的3D CT图像中骨组织的移除方法，其特征在于，所述步骤2的具体步骤如下：

步骤31，分类器训练：从步骤1得到的超体素中随机选取超体素构造训练集(x_i,y_i)，1<i<TrNum，其中，x_i表示训练集中第i个超体素对应的特征向量，yi表示超体素的类别标记，yi＝1表示超体素为骨组织，yi＝0表示超体素为非骨组织；TrNum表示训练集中超体素数目；使用构造的训练集训练分类器；

步骤32，对于步骤1得到的所有超体素中任一超体素，将步骤2提取得到的该超体素对应的特征向量x输入训练后的分类器，则得到该特征向量对应的标记值y，如果y＝1，则表示该超体素为骨组织，y＝0为非骨组织。

6.一种3D CT数据中骨组织的提取及移除系统，其特征在于，所述系统包括如下模块：

图像数据读取模块，用以读取3D CT图像数据序列，利用聚类法对其进行划分，并对划分后得到的类作为超体素进行标号；

超体素特征描述模块，用以对步骤1得到的每个标号后的超体素分别从超体素级、像素级和超体素邻域关系级进行描述，得到对每个超体素进行描述的特征向量；

分类器训练及超体素分类模块，用以构造超体素训练集，利用其中的超体素对应的特征向量对分类器进行训练，然后用训练后的分类器对步骤1得到的超体素进行分类并标记类别；

骨组织移除模块，用以依据步骤3中所得到的超体素标记值y，如果某一超体素被标记为骨组织，则从原始CT数据中对其进行移除，否则保留。

7.如权利要求1所述的3D CT数据中骨组织的提取及移除系统，其特征在于，所述图像数据读取模块具体实现如下流程的功能：

步骤1.1：读入3D CT图像数据序列I＝{p₁,p₂,…,p_N}，其中，pi表示3D CT图像中的第i个像素；p_i＝(g(p_i),x_i,y_i,z_i)，其中，g(p_i)表示像素p_i的CT值，idx(p_i)＝[x_i,y_i,z_i]表示像素p_i在3D CT图像数据中的位置坐标；N表示3D CT图像中像素的数目；

步骤1.2：输入骨组织CT值阈值区间T＝[t_l,t_h]、期望划分的超体素的数目K及迭代次数ITER；

步骤1.3：在3D CT图像数据序列中，随机选取K个聚类中心C_k＝{g(p_k),idx(p_k)}；其中，g(p_k)表示3D CT图像中的第k个像素p_k的CT值；idx(p_k)表示像素p_k在3D CT图像数据中的位置坐标；

步骤1.4：通过下式对每个聚类中心判断，得到聚类更新标记集U＝{u₁,u₂,…,u_K}；

$u_k=\left\{\begin{array}{cc}1,& \mathrm{if}{t}_l\&le;{\mathrm{gMin}}_k\&le;t_h\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中，gMin_k表示在第k个聚类中心C_k的S*S*S邻域中像素的CT值的最小值，S取2～3；u_k表示第k个聚类中心C_k的标记，k＝1～K；

步骤1.5：如果第k个聚类中心C_k的标记u_k＝1，则对该聚类中心执行步骤1.6；

步骤1.6：利用下式计算聚类中心C_k的M*M*M的邻域内每个像素对应的梯度值，M＝3～7；并将梯度值最小的像素作为新的聚类中心C_k；

$\&dtri;\left(I\right)={\left|\right|{\&dtri;}_x^2+{\&dtri;}_y^2+{\&dtri;}_z^2\left|\right|}_2^2,s.t.\left\{\begin{array}{c}{\&dtri;}_x=g(x+1,y,z)-g(x-1,y,z)\\ {\&dtri;}_y=g(x,y+1,z)-g(x,y-1,z)\\ {\&dtri;}_z=g(x,y,z+1)-g(x,y,z-1)\end{array}\right.$

其中，g(x,y,z)表示的是在空间位置(x,y,z)处的CT值，同时，及表示CT值在3D CT图像的中坐标系三个维度上的梯度值；

步骤1.7：对每个聚类中心C_k进行如下处理：在聚类中心C_k的2S*2S*2S邻域内选择与该聚类中心像素距离最近的像素，其中，S＝20～30；所述距离的计算根据下式进行：

$d=|g\left(p_i\right)-g\left(p_j\right)|+\frac{{\left|\right|\mathrm{idx}\left(p_i\right)-\mathrm{idx}\left(p_j\right)\left|\right|}^2}{M}$

其中，||idx(p_i)-idx(p_j)||²表示两个像素p_i与p_j在3D CT图像数据中的位置坐标之间的欧氏距离；M表示聚类中心对应的超体素中包含像素的数目；

然后将该距离最近的像素归入当前聚类中心所在的类，同时将该像素从之前类中移除，最后将聚类中心C_k使用其所属超体素的质心点替换；

步骤1.8：将步骤1.7执行ITER次，超体素划分过程结束，此时得到划分后的K个超体素，对它们进行顺序标号。

8.如权利要求7所述的3D CT数据中骨组织的提取及移除系统，其特征在于，所述步骤1.2中，所述骨组织CT值阈值T＝[960,1400]，迭代次数Iter＝3～8。

9.如权利要求1所述的3D CT数据中骨组织的提取及移除系统，其特征在于，所述超体素特征描述模块具体实现如下模块的功能：

步骤21，计算得到每个超体素中所有像素的CT值直方图H_G、CT值差值直方图H_D及CT值累积直方图H_C；然后根据CT值直方图H_G求取灰度均值、灰度方差、最大值、最小值、上四分位值、中位值、下四分位值、众值、熵、能量、对比度、二阶距、三阶距、四阶距和平滑性共15个特征值；根据CT值差值直方图H_D求取熵、能量、对比度、二阶距和平滑性共5个特征值；根据CT值累积直方图H_C求取熵、能量、对比度、二阶距和平滑性共5个特征值；其中，熵Et、能量Eg、对比度Ct及k阶中心距M^k的计算公式如下：

$\begin{array}{c}\mathrm{Et}=-{\mathrm{\&Sigma;}}_{i=1}^M{H}^{\left(i\right)}*\log \left(H^{\left(i\right)}\right)\\ \mathrm{Eg}={\mathrm{\&Sigma;}}_{i=1}^M{\left(H^{\left(i\right)}\right)}^2\\ \mathrm{Ct}={\mathrm{\&Sigma;}}_{i=1}^M{i}^2*H^{\left(i\right)}\\ M^k={\mathrm{\&Sigma;}}_{i=1}^M(i-\stackrel{\&OverBar;}{i})*H^{\left(i\right)}\end{array}---\left(4\right)$

其中，H表示CT值直方图H_G、CT值差直方图H_D或CT值累积直方图H_C，H⁽ⁱ⁾表示某直方图H中第i个位置对应的纵坐标值，M表示直方图的统计尺度；

同时，计算每种直方图的最小值、上四分位、中值、下四分位数及最大值；

步骤22，对步骤1得到的每个超体素中的每个像素，在像素的3*3*3邻域提取对应的3DSobel、梯度边缘滤波器、高斯平滑滤波器、Laplacian、线性锐化滤波器的响应值以及局部均值、局部中值、局部极大值和局部极小值共9个特征值；其中，局部均值是指上述像素的3*3*3邻域中所有26个像素的CT值的算术平均值；局部中值是指26个像素由小到大排序后位于中部的值，局部极大值和局部极小值同理也是指对该邻域的统计结果。

通过上述步骤得到第i个像素p_i的像素级特征f(p_i)后，通过下式得到超体素在像素级的特征描述：

f(V_s)^j＝avg(f(p_i)^j),1≤j≤N

其中，f(p_i)表示第i个像素的像素级特征；f(V_s)^j表示第s个超体素Vs中所有像素的像素级特征f(p_i)的均值；N表示该超体素中所包含的像素的数目；

步骤23，两个相邻的超体素间的邻域关系通过该两个超体素的CT值均值、CT值均值比、中心距离、两个超体素的CT值直方图H_G相关度、CT值差直方图H_D相关度及CT值累积直方图H_C相关度共6个特征表示，两个超体素的某种直方图的相关度采用如下式所示的卡方距离进行描述：

${\mathrm{\&chi;}}^2(h_1,h_2)=\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^M\frac{{(h_1^{\left(i\right)}-h_2^{\left(i\right)})}^2}{h_1^{\left(i\right)}+h_2^{\left(i\right)}}$

其中，及分别表示直方图h₁及h₂的第i个元素；M表示直方图的统计尺度。

经过上述步骤，得到对每个超体素进行描述的特征向量。

10.如权利要求1所述的3D CT图像中骨组织的移除方法，其特征在于，所述分类器训练及超体素分类模块具体实现如下流程的功能：

步骤31，分类器训练：从步骤1得到的超体素中随机选取超体素构造训练集(x_i,y_i)，1<i<TrNum，其中，x_i表示训练集中第i个超体素对应的特征向量，yi表示超体素的类别标记，yi＝1表示超体素为骨组织，yi＝0表示超体素为非骨组织；TrNum表示训练集中超体素数目；使用构造的训练集训练分类器；

步骤32，对于步骤1得到的所有超体素中任一超体素，将步骤2提取得到的该超体素对应的特征向量x输入训练后的分类器，则得到该特征向量对应的标记值y，如果y＝1，则表示该超体素为骨组织，y＝0为非骨组织。