CN110349168B - 一种股骨头ct影像的分割方法 - Google Patents

Abstract

本发明提供一种股骨头CT影像的分割方法，涉及医学图像处理技术领域。该方法使用三维最大类间方差法进行预分割，然后基于图割和形状约束相结合进行自动股骨头精确分割，构建图后，基于Graph cuts模型优化分割结果，基于分层Hough变换检测影像中圆形区，使用SVM对Graph cuts的分割结果进行重新预测、分类，提取邻域梯度特征，分离股骨头和髋臼，以检测出的圆心为种子节点，使用区域生长算法生成股骨头，得到最终的股骨头分割图像。本发明不仅能够有效剔除影像噪声，为Graph cuts模型提供硬约束条件，鲁棒性较好，实现全自动分割股骨头CT影像，还能大大缩短Graph cuts模型的收敛时间，分割出的股骨头边缘完整、细节清晰，分割准确率达到92％。

Description

一种股骨头CT影像的分割方法

技术领域

本发明涉及医学图像处理技术领域，尤其涉及一种股骨头CT影像的分割方法。

背景技术

影像分割在医学影像的定量、定性分析中均扮演着十分重要的角色，它直接影响到计算机辅助诊断系统的后续分析和处理工作。正确分割股骨头的影像，不仅可以通过股骨头的形状来确定患者坏死的程度，还可以通过分割结果近似地求出患者股骨头内部的缺血体积，为股骨头坏死的辅助诊断和分期判断做好准备。目前，股骨头CT影像的分割方法主要包括专家手工分割、计算机交互式分割和全自动分割。手工分割和计算机交互式分割对专家知识和经验要求很高，而且存在不可避免人为误差，同时对海量CT数据进行手工处理是一件耗时的事情，因此，股骨头CT影像的全自动分割具有极大的研究意义和价值。

目前常用的全自动股骨头分割方法主要有阈值法、分水岭法、水平集法和图谱法。其中阈值法单纯依赖图像像素信息进行分割，存在忽略图像噪声和边界处低对比度的缺点，在分割股骨头CT影像时会出现股骨头轮廓不完整、内部存在较大空洞和噪声亮斑过多等问题；分水岭法会受到CT图像中细节纹理和噪声的影响，出现过分割现象，影响分割效果，误分割率高；水平集法对于初始点的选择要求较为苛刻，并且分割速度较慢，分割准确率低；图谱法需要大量的训练样本，当训练样本与测试样本差异过大时，不能完成准确的分割。

发明内容

本发明要解决的技术问题是针对上述现有技术的不足，提供一种股骨头CT影像的分割方法，使用三维最大类间方差法进行预分割，然后基于图割和形状约束相结合进行自动股骨头精确分割，不仅能够有效剔除影像噪声、为Graph cuts模型提供硬约束条件，还能大大缩短Graph cuts模型的收敛时间。

为解决上述技术问题，本发明所采取的技术方案是：

一种股骨头CT影像的分割方法，包括基于三维类间方差的股骨头预分割方法和基于图割和形状约束相结合的自动股骨头精确分割方法；

基于三维类间方差的股骨头预分割方法中，使用三维最大类间方差法进行预分割，根据分割结果，取出骨像素集合中灰度值最高的10％像素和非骨像素集合中灰度值最低的10％非零像素，作为Graph cuts模型的硬约束条件；

基于图割和形状约束相结合的自动股骨头精确分割方法包括以下步骤：

步骤2.1：构建图；

首先把待处理影像转化为符合Graph cuts模型的图网络，具体方法为：设置两个目标节点α和

α值为所有类别为α的像素灰度值的均值，

的值为所有类别不是α的像素灰度值的均值；再将所有像素和两个目标节点α和

分别连接起来，目标节点与像素之间的连线称为t-link；对于所有四相邻的像素对，若经过预分割后像素标签相同，则连接像素，连线称为n-link，若经过预分割后像素标签不同，则在两个像素之间添加辅助节点a，并连接辅助节点和两个像素，再连接辅助节点与

如此完成构图，把一个图像转化为图网络；

然后给图网络中所有的边赋予权值，在图中有四种边，像素之间的连线n-link、目标节点与像素的连线t-link、像素与辅助节点的连线e{p,a}和辅助节点与

的连线

其中p表示像素点；n-link的权值的大小描述像素属于同一标签的概率，使用训练好的PixelsPair_Category树来预测；t-link的权值描述像素属于α和

两类的概率，使用Pixel_Category树预测；e{p,a}的权值等于像素p和a属于同一类的概率，用PixelsPair_Category树预测；

的权值等于辅助节点两端像素的相似性，用PixelsPair_Category树预测；

步骤2.2：基于Graph cuts模型优化分割结果，得到大部分像素分类正确、股骨头的圆形轮廓容易辨认的分割结果；

步骤2.2.1：进行初始化；

利用属于标号集合T₀、T₁的像素形成训练集：

data set1＝{l₁,l₂,…,l_n}，label set＝{0,1}；

其中，l_i表示像素灰度值，n表示像素总个数，label set表示标签集，方法如下：

训练随机森林模型，得到Pixel_Category树；然后用Pixel_Category树对整个图像进行预测，得到整幅图像的标号矩阵initf；然后通过标号集T₀、T₁及其包含的像素，组成像素对(p,q)，如果像素对(p,q)属于同一标号集合，则设定标签为0，否则，如果像素对(p,q)属于不同标号集合，则设定标签为1，最终得到训练集：

data＝{(l₁,l₂),(l₂,l₃),…,(l_n-1,l_n)}，label＝{0,1}；

其中，l_i表示像素灰度值，n表示像素总个数，label表示标签集；使用该训练集来训练随机森林模型得到PixelsPair_Category树，该树用于衡量像素对属于同一标签和不同标签的概率，即度量像素对中两个像素的相似性；

步骤2.2.2：使用α-expansion算法进行迭代；

步骤2.2.2.1：设置迭代的标志“continue”为False；

步骤2.2.2.2：对于每个目标节点标号α，在初始标签initf的基础上实施一次α-expansion操作，得到新的标号矩阵f′；在所有的标号矩阵f′中找到总能量最小的矩阵，即：

当新的标号矩阵f′的总能量E(f′)小于初始标号矩阵initf的总能量E(initf)时，将f′赋值给initf，并且把迭代的标志“continue”修改为True；然后根据新的结果，更新训练集，修正Pixel_Category树；

步骤2.2.2.3：如果continue＝True，则回到步骤2.2.2.1；否则，返回分割结果Seg，得到一阶段分割结果；

步骤2.3：基于分层Hough变换检测影像中圆形区；

确定基于Graph cuts模型的分割方法生成的图像mat2中股骨头近似圆的圆心和半径，使用分层Hough变换检测圆形：第一步，采用传统的梯度法在每层切片上检测圆，从中选出半径最大的圆，该圆在切片上的x、y坐标就是三维股骨头球心的x、y坐标；第二步，利用球心的x、y坐标，在每个切片上估计出圆半径，再利用圆半径计算出球心的z坐标和球半径r，每计算出一组(z,r)，zr平面上对应的(z,r)点的值就加1，默认所有点的初始值为0，最后找出值最大的一组(z,r)，对应的z坐标值就是球心的z坐标，对应的r值就是球半径r；确定球心和球半径之后，利用切片之间的距离计算出其他切片上圆的半径，得到完整的圆；

梯度法检测圆的具体方法为：对每个边界点求梯度，再以该边界点为起点，沿梯度方向作射线，对射线经过的像素点都进行累加，累加值符合阈值要求的点即为圆心；

步骤2.4：股骨头和髋臼分离；

使用支持向量机(SVM)对Graph cuts的分割结果进行重新预测、分类，把像素的灰度值f(x,y)和圆心距离d组成新的二维特征向量[f(x,y),d]，加上骨和非骨的标签0或1组成训练集，训练支持向量机，再用训练后的支持向量机预测整幅图像，得到结果图像mat3；

选择以像素p为中心，7×7的核来提取邻域梯度特征，在像素p的八个方向，若存在灰度值比p大20或以上的像素，则令该方向的特征值为-1，否则为1；只对两个圆心之间的骨类像素进行提取特征、分类，分类的依据是若存在同一直线的两个特征值均为-1，那么该像素为非骨类，否则，仍然为骨类；如此得到股骨头和髋臼分离的图像矩阵mat4；

步骤2.5：得到mat4后，以之前检测出的圆心为种子节点，使用区域生长算法生成股骨头，得到最终的股骨头分割图像mat5。

进一步地，基于三维类间方差的股骨头预分割方法的具体步骤如下：

步骤1.1：提取股骨头CT图像中所有非零像素的坐标存入列表location，计算其邻域均值g(x,y)和邻域中值h(x,y)，与灰度值f(x,y)组成向量

步骤1.2：利用三维最大类间方差法求出向量

组成的邻域均值、邻域中值和灰度值的数据阈值[t^*,s^*,q^*]，并且根据阈值将所有像素分成两类，即骨头和非骨头，骨头标号为1，非骨头标号为0，得到标号矩阵mat1和属于每个标号的像素集合T₀、T₁；

步骤1.3：将标号为1的像素中灰度值最大的10％设置为硬约束条件，即在分割过程中标号一直为1，不发生变化；同样的，标号为0的非零像素中灰度值最小的10％也设置为硬约束条件。

进一步地，所述步骤2.2.2.2中找到在所有的标号矩阵f′中找到总能量最小的矩阵的方法为：

①对图像中的每个像素使用Pixel_Category树来预测像素属于每种标号的概率，得到一个储存所有像素分属每种标号的概率的矩阵，称为Pixel_Prob矩阵；

②收集所有标号为α的像素，计算像素灰度值的均值作为目标节点α的值；

③将每个像素和标号α形成像素对，然后使用“PixelsPair_Category”树来预测每个像素对，也就是每个像素和标号α属于同一类的概率和不同类的概率，得到一个矩阵，称为PixelPair_Prob矩阵；

④依据最大流/最小割算法求出此时图像的最小割，并得到割集All_cuts，此时的标号矩阵为f′；

⑤根据割集All_cuts，使用广度优先遍历方法得到分割后的图像；

⑥计算

和E(initf)。

采用上述技术方案所产生的有益效果在于：本发明提供的股骨头CT影像的分割方法，该方法鲁棒性较好、分割时间短，能实现全自动分割股骨头CT影像，分割出的股骨头边缘完整、细节清晰，分割准确率达到92％(DICE指数)。

附图说明

图1为本发明实施例提供的股骨头CT影像的分割方法流程图；

图2为本发明实施例提供的α-expansion的图网络模型示意图；

图3为本发明实施例提供的邻域梯度特征提取过程示意图；

图4为本发明实施例提供的各步骤分割效果示意图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

如图1所示，本实施例的方法如下所述。

步骤1：基于三维类间方差的股骨头预分割方法；

由于基于Graph Cuts模型的分割方法是一种交互式的分割方法，用户需要提前将某些像素标记为“对象”或“背景”，为分割提供硬约束条件，因此使用三维最大类间方差法进行预分割，根据分割结果，取出骨像素集合中灰度值最高的10％像素和非骨像素集合中灰度值最低的10％非零像素，作为Graph cuts模型的硬约束条件。使用三维最大类间方差法进行预分割，不仅能够有效的剔除影像噪声、为Graph cuts模型提供硬约束条件，还能大大缩短Graph cuts模型的收敛时间。

具体步骤如下：

步骤1.1：提取股骨头CT图像中所有非零像素的坐标存入列表location，计算其邻域均值g(x,y)和邻域中值h(x,y)，与灰度值f(x,y)组成向量

步骤1.2：利用三维最大类间方差法求出向量

组成的邻域均值、邻域中值和灰度值的数据阈值[t^*,s^*,q^*]，并且根据阈值将所有像素分成两类，即骨头和非骨头，骨头标号为1，非骨头标号为0，得到标号矩阵mat1和属于每个标号的像素集合T₀、T₁；

步骤1.3：将标号为1的像素中灰度值最大的10％设置为硬约束条件，即在分割过程中标号一直为1，不发生变化；同样的，标号为0的非零像素中灰度值最小的10％也设置为硬约束条件。

本实施例中，原始股骨头CT图像如图4中的原图所示，预分割后的图像为mat1。

步骤2：基于图割和形状约束相结合的自动股骨头精确分割方法，在得到标号矩阵mat1之后，使用基于Graph Cuts模型的分割方法，在预分割的基础上，将影像中的股骨头和髋臼的整体更准确地分割出来。具体包括以下步骤：

步骤2.1：构建图；

首先把待处理影像转化为符合Graph cuts模型的图网络，具体方法为：设置两个目标节点α和

α的值为所有类别为α的像素灰度值的均值，

的值为所有类别不是α的像素灰度值的均值；再将所有像素和两个目标节点α和

分别连接起来，目标节点与像素之间的连线称为t-link；对于所有四相邻的像素对，若经过预分割后像素标签相同，则连接像素，连线称为n-link，若经过预分割后像素标签不同，则在两个像素之间添加辅助节点a，并连接辅助节点和两个像素，再连接辅助节点与

如此完成构图，把一个图像转化为图网络；如图2所示；

然后给图网络中所有的边赋予权值，在图中有四种边，像素之间的连线n-link、目标节点与像素的连线t-link、像素与辅助节点的连线e{p,a}和辅助节点与

的连线，其中p表示像素点；n-link的权值的大小描述像素属于同一标签的概率，

使用训练好的PixelsPair_Category树来预测；t-link的权值描述像素属于α和

两类的概率，使用Pixel_Category树预测；e{p,a}的权值等于像素p和a属于同一类的概率，用PixelsPair_Category树预测；

的权值等于辅助节点两端像素的相似性，用PixelsPair_Category树预测；边的权值如表1所示，其中，函数D(·)指代Pixel_Category树，函数V(·)指代PixelsPair_Category树，N是相邻像素集，f_p和f_q是像素p和q的标签，P_α是标号为α的像素集合。

表1α-expansion的图网络模型的权值

步骤2.2：基于Graph cuts模型优化分割结果，得到大部分像素分类正确、股骨头的圆形轮廓容易辨认的分割结果；

步骤2.2.1：进行初始化；

利用属于标号集合T₀、T₁的像素形成训练集：

data set1＝{l₁,l₂,…,l_n},label set＝{0,1}；

其中，l_i表示像素灰度值，n表示像素总个数，label set表示标签集，方法如下：

训练随机森林模型，得到Pixel_Category树；然后用Pixel_Category树对整个图像进行预测，得到整幅图像的标号矩阵initf；然后通过标号集T₀、T₁及其包含的像素，组成像素对(p,q)，如果像素对(p,q)属于同一标号集合，则设定标签为0，否则，如果像素对(p,q)属于不同标号集合，则设定标签为1，最终得到训练集：

data＝{(l₁,l₂),(l₂,l₃),…,(l_n-1,l_n)}，label＝{0,1}；

其中，l_i表示像素灰度值，n表示像素总个数，label表示标签集。使用训练集来训练随机森林得到PixelsPair_Category树，该树用于衡量像素对属于同一标签和不同标签的概率，即度量像素对中两个像素的相似性；

本实施例中的图像分辨率较高，训练数据量较大，因此把所有的相邻像素对作为训练集合；

步骤2.2.2：使用α-expansion算法进行迭代；

步骤2.2.2.1：设置迭代的标志“continue”为False；

步骤2.2.2.2：对于每个目标节点标号α，在初始标签initf的基础上实施一次α-expansion操作，得到新的标号矩阵f′；在所有的标号矩阵f′中找到总能量最小的矩阵，即：

当新的标号矩阵f′的总能量E(f′)小于初始标号矩阵initf的总能量E(initf)时，将f′赋值给initf，并且把迭代的标志“continue”修改为True；然后根据新的结果，更新训练集，修正Pixel_Category树；

找到在所有的标号矩阵f′中找到总能量最小的矩阵的方法为：

①对图像中的每个像素使用Pixel_Category树来预测像素属于每种标号的概率，得到一个储存所有像素分属每种标号的概率的矩阵，称为Pixel_Prob矩阵；本实施例中的CT图像分辨率是512×512，两种标号，就可以得到一个(512,512,2)的矩阵，储存了所有像素分属两类的概率；

②收集所有标号为α的像素，计算像素灰度值的均值作为目标节点α的值；

③将每个像素和标号α形成像素对，然后使用“PixelsPair_Category”树来预测每个像素对，也就是每个像素和标号α属于同一类的概率和不同类的概率，得到一个矩阵，称为PixelPair_Prob矩阵；本实施例中，得到一个(512,512,2)的PixelPair_Prob矩阵；

④依据最大流/最小割算法求出此时图像的最小割，并得到割集All_cuts，此时的标号矩阵为f′；

⑤根据割集All_cuts，使用广度优先遍历方法得到分割后的图像；

⑥计算

和E(initf)；

步骤2.2.2.3：如果continue＝True，则回到步骤2.2.2.1；否则，返回分割结果Seg，得到一阶段分割结果；

通过Graph cuts模型，得到了大部分像素分类正确、股骨头的圆形轮廓容易辨认的分割结果mat2，如图4中的mat2图像所示，但是存在着股骨头边缘有缺口、内部有空洞的问题，需要继续优化结果。

步骤2.3：基于分层Hough变换检测影像中圆形区；

由基于Graph cuts模型的分割方法生成的图像mat2可以看出，单纯依靠纹理或像素特征分割得到的结果，边缘线不完整，不准确，只能看出粗略的股骨头形状，要进一步优化结果。因为人体的许多器官组织，尤其是骨骼等硬组织，往往具有比较明显的形状，如果能在分割过程中加入对形状的约束,将更有利于对器的精确分割。由于人体中股骨头近似球形的特点，因此本实施例选择圆形作为每个切片的形状约束。

首先确定基于Graph cuts模型的分割方法生成的图像mat2中股骨头近似圆的圆心和半径，使用分层Hough变换检测圆形，相较于传统的四维空间Hough变换，其时间效率和空间效率都大大提高。分层Hough变换检测圆形的思路是将整个检测过程分位两个步骤：第一步，采用传统的梯度法在每层切片上检测圆，从中选出半径最大的圆，该圆在切片上的x、y坐标就是三维股骨头球心的x、y坐标；本实施例采用传统的梯度法检测圆，即对每个边界点求梯度，再以该边界点为起点，沿梯度方向作射线，对射线经过的像素点都进行累加，累加值符合阈值要求的点即为圆心；第二步，利用球心的x、y坐标，在每个切片上估计出圆半径，再利用圆半径计算出球心的z坐标和球半径r，每计算出一组(z,r)，zr平面上对应的(z,r)点的值就加1(默认所有点的初始值为0)，最后找出值最大的一组(z,r)，对应的z坐标值就是球的z坐标，对应的r值就是球半径；确定球心和球半径之后，利用切片之间的距离计算出其他切片上圆的半径，得到完整的圆。

步骤2.4：股骨头和髋臼分离；

检测到圆后，要解决股骨头轮廓不完整，骨头内部存在较大空洞的问题，由于本实施例是一个二分类问题，并且每幅图像的实验样本较小，因此选择使用支持向量机(SVM)对基于Graph cuts的分割结果mat2进行重新预测、分类，把像素的灰度值f(x,y)和圆心距离d组成新的二维特征向量[f(x,y),d]，加上骨和非骨的标签0或1组成训练集，训练支持向量机，再用训练后的支持向量机预测整幅图像，得到结果图像mat3，如图4中的mat3图像所示；

得到的图像mat3已经解决了股骨头轮廓不完整的问题，内部空洞也会缩小或补足。由于人体中股骨头和髋臼连接紧密，并且医疗检测时患者CT拍摄角度的不同，在CT图像中股骨头和髋臼出现了相连或重叠的现象，难以分割。为了分离股骨头和髋臼，本实施例提取新的像素特征。观察原图可以发现，股骨头和髋臼存在缝隙，并且缝隙像素两边是灰度值较高的骨像素，并且缝隙较窄，因此选取像素的邻域梯度特征作为分离股骨头和髋臼的凭据。因为目的是明确股骨头和髋臼的间隙，这个分类过程只适用于左右两个股骨头圆心之间的像素。

选择以像素p为中心，7×7的核来提取邻域梯度特征，在像素p的八个方向，若存在灰度值比p大20或以上的像素，则令该方向的特征值为-1，否则为1；例如，若p11、p22、p33中存在像素，其灰度值比p的灰度值大20或者更多，那么该方向的特征值q＝-1。又因为提取特征的目的是检测出股骨头和髋臼缝隙的像素，因此只对两个圆心之间的骨类像素进行提取特征、分类，分类的依据是若存在同一直线的两个特征值均为-1，那么该像素为非骨类，否则，仍然为骨类；例如，在图3中，像素p的特征值中，如果z＝-1,e＝-1，那么像素p归为非骨类别。如此得到股骨头和髋臼分离的图像矩阵mat4，如图4中的mat4图像所示。

步骤2.4：得到mat4后，可以发现股骨头和髋臼已经完全分离开，以之前检测出的圆心为种子节点，使用区域生长算法生成股骨头，得到最终的股骨头分割图像mat5，如图4中的mat5图像所示。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (3)

1.一种股骨头CT影像的分割方法，其特征在于：包括基于三维类间方差的股骨头预分割方法和基于图割和形状约束相结合的自动股骨头精确分割方法；

基于三维类间方差的股骨头预分割方法中，使用三维最大类间方差法进行预分割，根据分割结果，取出骨像素集合中灰度值最高的10％像素和非骨像素集合中灰度值最低的10％非零像素，作为Graph cuts模型的硬约束条件；

基于图割和形状约束相结合的自动股骨头精确分割方法包括以下步骤：

步骤2.1：构建图；

首先把待处理影像转化为符合Graph cuts模型的图网络，具体方法为：设置两个目标节点α和

α值为所有类别为α的像素灰度值的均值，

的值为所有类别不是α的像素灰度值的均值；再将所有像素和两个目标节点α和

分别连接起来，目标节点与像素之间的连线称为t-link；对于所有四相邻的像素对，若经过预分割后像素标签相同，则连接像素，连线称为n-link，若经过预分割后像素标签不同，则在两个像素之间添加辅助节点a，并连接辅助节点和两个像素，再连接辅助节点与

如此完成构图，把一个图像转化为图网络；

然后给图网络中所有的边赋予权值，在图中有四种边，像素之间的连线n-link、目标节点与像素的连线t-link、像素与辅助节点的连线e{p，a}和辅助节点与

的连线

其中p表示像素点；n-link的权值的大小描述像素属于同一标签的概率，使用训练好的PixelsPair_Category树来预测；t-link的权值描述像素属于α和

两类的概率，使用Pixel_Category树预测；e{p，a}的权值等于像素p和a属于同一类的概率，用PixelsPair_Category树预测；

的权值等于辅助节点两端像素的相似性，用PixelsPair_Category树预测；

步骤2.2：基于Graph cuts模型优化分割结果，得到大部分像素分类正确、股骨头的圆形轮廓容易辨认的分割结果；

步骤2.2.1：进行初始化；

利用属于标号集合T₀、T₁的像素形成训练集：

data set1＝{l₁，l₂，...，l_n}，label set＝{0，1}；

其中，l_i表示像素灰度值，n表示像素总个数，label set表示标签集，方法如下：

训练随机森林模型，得到Pixel_Category树；然后用Pixel_Category树对整个图像进行预测，得到整幅图像的标号矩阵initf；然后通过标号集T₀、T₁及其包含的像素，组成像素对(p，q)，如果像素对(p，q)属于同一标号集合，则设定标签为0，否则，如果像素对(p，q)属于不同标号集合，则设定标签为1，最终得到训练集：

data＝{(l₁，l₂)，(l₂，l₃)，...，(l_n-1，l_n)}，label＝{o，1}；

其中，l_i表示像素灰度值，n表示像素总个数，label表示标签集；使用该训练集来训练随机森林模型得到PixelsPair_Category树，该树用于衡量像素对属于同一标签和不同标签的概率，即度量像素对中两个像素的相似性；

步骤2.2.2：使用α-expansion算法进行迭代；

步骤2.2.2.1：设置迭代的标志“continue”为False；

步骤2.2.2.2：对于每个目标节点标号α，在初始标签initf的基础上实施一次α-expansion操作，得到新的标号矩阵f′；在所有的标号矩阵f′中找到总能量最小的矩阵，即：

当新的标号矩阵f′的总能量E(f′)小于初始标号矩阵initf的总能量E(initf)时，将f′赋值给initf，并且把迭代的标志“continue”修改为True；然后根据新的结果，更新训练集，修正Pixel_Category树；

步骤2.2.2.3：如果continue＝True，则回到步骤2.2.2.1；否则，返回分割结果Seg，得到一阶段分割结果；

步骤2.3：基于分层Hough变换检测影像中圆形区；

确定基于Graph cuts模型的分割方法生成的图像mat2中股骨头近似圆的圆心和半径，使用分层Hough变换检测圆形：第一步，采用传统的梯度法在每层切片上检测圆，从中选出半径最大的圆，该圆在切片上的x、y坐标就是三维股骨头球心的x、y坐标；第二步，利用球心的x、y坐标，在每个切片上估计出圆半径，再利用圆半径计算出球心的z坐标和球半径r，每计算出一组(z，r)，zr平面上对应的(z，r)点的值就加1，默认所有点的初始值为0，最后找出值最大的一组(z，r)，对应的z坐标值就是球心的z坐标，对应的r值就是球半径r；确定球心和球半径之后，利用切片之间的距离计算出其他切片上圆的半径，得到完整的圆；

梯度法检测圆的具体方法为：对每个边界点求梯度，再以该边界点为起点，沿梯度方向作射线，对射线经过的像素点都进行累加，累加值符合阈值要求的点即为圆心；

步骤2.4：股骨头和髋臼分离；

使用支持向量机(SVM)对Graph cuts的分割结果进行重新预测、分类，把像素的灰度值f(x，y)和圆心距离d组成新的二维特征向量[f(x，y)，d]，加上骨和非骨的标签0或1组成训练集，训练支持向量机，再用训练后的支持向量机预测整幅图像，得到结果图像mat3；

选择以像素p为中心，7×7的核来提取邻域梯度特征，在像素p的八个方向，若存在灰度值比p大20或以上的像素，则令该方向的特征值为-1，否则为1；只对两个圆心之间的骨类像素进行提取特征、分类，分类的依据是若存在同一直线的两个特征值均为-1，那么该像素为非骨类，否则，仍然为骨类；如此得到股骨头和髋臼分离的图像矩阵mat4；

步骤2.5：得到mat4后，以之前检测出的圆心为种子节点，使用区域生长算法生成股骨头，得到最终的股骨头分割图像mat5。

2.根据权利要求1所述的股骨头CT影像的分割方法，其特征在于：所述基于三维类间方差的股骨头预分割方法的具体步骤如下：

步骤1.1：提取股骨头CT图像中所有非零像素的坐标存入列表location，计算其邻域均值g(x，y)和邻域中值h(x，y)，与灰度值f(x，y)组成向量

步骤1.2：利用三维最大类间方差法求出向量

组成的邻域均值、邻域中值和灰度值的数据阈值[t^*，s^*，q^*]，并且根据阈值将所有像素分成两类，即骨头和非骨头，骨头标号为1，非骨头标号为0，得到标号矩阵mat1和属于每个标号的像素集合T₀、T₁；

步骤1.3：将标号为1的像素中灰度值最大的10％设置为硬约束条件，即在分割过程中标号一直为1，不发生变化；同样的，标号为0的非零像素中灰度值最小的10％也设置为硬约束条件。

3.根据权利要求2所述的股骨头CT影像的分割方法，其特征在于：所述步骤2.2.2.2中找到在所有的标号矩阵f′中找到总能量最小的矩阵的方法为：

①对图像中的每个像素使用Pixel_Category树来预测像素属于每种标号的概率，得到一个储存所有像素分属每种标号的概率的矩阵，称为Pixel_Prob矩阵；

②收集所有标号为α的像素，计算像素灰度值的均值作为目标节点α的值；

③将每个像素和标号α形成像素对，然后使用“PixelsPair_Category”树来预测每个像素对，也就是每个像素和标号α属于同一类的概率和不同类的概率，得到一个矩阵，称为PixelPair_Prob矩阵；

④依据最大流/最小割算法求出此时图像的最小割，并得到割集All_cuts，此时的标号矩阵为f′；

⑤根据割集All_cuts，使用广度优先遍历方法得到分割后的图像；

⑥计算

和E(initf)。

Images