CN110246218A - 股骨头三维模型的重建方法、空间骨盆参数测量方法 - Google Patents

Abstract

本发明一种股骨头三维模型的重建方法，包括以下步骤：将包含股骨头的原始二维CT图像进行二值化处理；将含股骨头的二值图采用非梯度最大值抑制算法获得边缘轮廓；利用高斯卷积核得到每个像素的梯度向量；获得的边缘轮廓中的所有像素均为三角形候选像素；计算梯度模值，得到构成一个圆的像素对；对步骤S4中每对像素构成的圆的圆心点进行聚类，识别出股骨头区域；采用聚类算法对S5中识别出的股骨头区域进行聚类，得到三维股骨头坐标，重建股骨头三维模型。本发明还公开一种空间骨盆参数测量方法。本发明具有将二维参数扩展到三维参数，可以拓宽椎体研究的范围，更符合临床实践；消除用户交互需求，保证框架更高效、可靠、准确，技术要求更低的优点。

Description

股骨头三维模型的重建方法、空间骨盆参数测量方法

技术领域

本发明涉及医学图像处理、深度学习技术领域，尤其涉及股骨头三维模型的重建方法、空间骨盆参数测量方法。

背景技术

骨盆是连结脊柱和下肢之间的盆状骨架，由后方的骶、尾骨(脊柱最低的两块骨)和左右两髋骨连接而成的完整骨环。

矢状位脊柱骨盆参数的确定，基于空间股骨头中心定位。目前，空间股骨头中心自动定位是在CT数据自动空间重建股骨头的基础上实现的。CT数据中股骨头位置的检测是重建中最重要的一步。一般来说，空间股骨头中心的自动定位过程分为两个阶段：从CT数据中检测股骨头区域和空间股骨头中心的检测。基于改进的随机采样圆检测算法检测股骨头区域，利用半径滤波获得股骨头空间中心。对图像的一半执行整个过程，检测一个股骨头，然后对图像的另一半执行相同的步骤。

股骨头通常被认为是球形的。在CT切片中，股骨头检测可以转换为圆形检测过程。圆检测的方法很多，但大多数圆检测算法都是基于霍夫变换的，该算法耗时、耗内存。

另外，骨盆参数的测量主要采用人工通过手动在二维图像上测量骨盆参数，即，通过在计算机辅助软件上对二维图像进行手动标定和测量，在测量过程中往往需要在不同切片中进行来回切换来找到准确的测量点；使得测量工作效率低下，技术要求高，结果不可靠。

矢状位脊柱骨盆参数的准确测量是顺利实施脊柱矫形手术重要条件，临床上所用的骨盆参数测量方法以手工测量平面CT图像为主，且整个过程人工干预较多，造成参数估计准确性下降。

发明内容

本发明旨在提供一种基于股骨头三维模型的重建方法、空间骨盆参数测量方法，通过利用病人的CT图像，能完成对病人的包括骨盆入口平面面积等参数的自动测量。

本发明通过以下技术手段实现解决上述技术问题的：一种股骨头三维模型的重建方法，包括以下步骤：

S1、将包含股骨头的原始二维CT图像进行二值化处理，得到含股骨头的二值图；

S2、将含股骨头的二值图采用非梯度最大值抑制算法获得边缘轮廓；

S3、利用高斯卷积核在对边缘轮廓图像在水平和垂直方向上进行卷积，得到每个像素的梯度向量；

S4、获得的边缘轮廓中的所有像素均为三角形候选像素；计算每对像素的梯度模值，得到构成一个圆的像素对；

S5、对步骤S4中每对像素构成的圆的圆心点进行聚类，识别出股骨头区域；

S6、采用聚类算法对S5中识别出的股骨头区域进行聚类，得到三维股骨头坐标，重建股骨头三维模型。

优选地，所述S3采用包括以下方式得到梯度向量：

Grad_x(x,y)＝G_x(x,y)*I_edge

Grad_y(x,y)＝G_y(x,y)*I_edge

其中，Grad_x(x，y)，Grad_x(x，y)分别表示每个像素在水平和垂直方向的梯度值，G_x(x，y)，G_y(x，y)为水平和垂直方向卷积核，卷积核采用sobel算子，I_edge为边缘轮廓图像，水平和垂直两个方向梯度构成一个梯度向量

Grad(x，y)＝(Grad_x(x，y)，Grad_x(x，y))。

优选地，所述S4采用包括以下方式获得构成一个圆的像素对：

比较每一对像素的梯度模值|Grad(x，y)|的大小和构成三角形时两个底角的角度大小；角度相同的被认为是圆上的一对点；该对像素能构成一个圆；

边缘轮廓中两个像素P₁、P₂之间的连接线用表示，以为边构造等腰三角形，圆上像素的梯度向量都是指向圆心的，将归一化梯度向量用表示；构成三角形的两个底角θ₁、θ₂的计算方法如下：

优选地，所述S4采用相似性度量构成一个圆的一对像素的准确率；

相似性度量采用余弦形式，如下：

|cos(θ₁)-cos(θ₂)|≤ε

其中，ε是容错控制项，它的大小控制了等腰三角形识别的准确率，上式成立时，默认对应的P₁、P₂为同一个圆上的两个点，该圆的圆心为等腰三角形的顶点。

优选地，所述S5包括以下步骤：

S51、初始化聚类簇，簇的数量设定为0；将S4得到的所有圆心点作为圆心点集合C；

S52、在圆心点集合C中选择一个圆心点并遍历集合中所有的圆心点，计算当前圆心点与其他圆心点之间的距离；

S53、将当前圆心点和所有使S52的距离小于设定阈值的圆心聚类到同一个簇中；

S54、去除已聚类圆心点并更新圆心点集合C；

S55、直至圆心点集合为空，则完成圆心点聚类；否则，转S52；

统计各个聚类簇中的边长，拥有最多最长边长的聚类簇即为股骨头区域。

优选地，所述S6包括以下步骤：

S61、初始化聚类簇，簇的数量设定为0；

S62、读取识别的S5中的股骨头区域，并遍历所有的簇的中心点，计算股骨头中心点与簇中心点的距离；

S63、若S62的距离小于设定的阈值，将读出的股骨头区域聚类到相应的簇中，簇内股骨头区域数量加1，更新簇内中心为所有股骨头区域中心的平均值，更新最大半径，更新具有最大半径股骨头区域的位置；

S64、若S63不成立，则新建聚类簇，重新依次执行S62、S63；

S65、直至股骨头区域队列为空，则转S66；否则，转S62；

S66、挑选拥有最多股骨头的聚类簇，则最大聚类簇的中心坐标作为三维股骨头中心点的X-Y平面坐标，最大半径所在的位置作为股骨头中心点Z轴坐标，在三维模型中空间股骨头位置绘制球体拟代股骨头，球体半径为在二维平面上识别出的最大半径。

本发明还公开一种基上述的股骨头三维模型的重建方法的空间骨盆参数测量方法，其特征在于，包括以下步骤：

步骤一、将原始二维CT图像进行分割处理，并重建三维骨盆模型；

步骤二、重建股骨头三维模型；

步骤三、建立VGG16网络模型，并对VGG16网络模型进行训练，直至其收敛；

步骤四、将包含L5的原始二维CT图像、S1的原始二维CT图像和股骨头的原始二维CT图像输入至收敛后的步骤三中的VGG16网络模型中，按序识别输入的图像，根据图像类别预测骶骨面位置；

步骤五、对步骤四找到含有骶骨的二维CT图像运行4连通区域识别算法，根据最大连通区域找到用于映射的锚点，即最大连通区域的部分上边缘点；将所述锚点映射回重建的步骤一的三维骨盆模型，生成S1空间骶骨面模型；

步骤六、计算骨盆参数，取股骨头三维模型的空间中心的参数和S1空间骶骨面模型上的参数进行骨盆参数的计算。

优选地，在所述步骤一中，采用基于加权质量评价函数的K-MEANS算法进行二维CT图像分割处理，利用MC算法重建三维骨盆模型。

优选地，所述步骤一采用基于加权质量评价函数的K-MEANS算法进行二维CT图像分割处理，具体方法如下：

首先输入待分割的二维CT图像，之后进行灰度化，使用基于信息熵的迭代算法初始化K个聚簇中心；然后利用下式对图像中的每个像素点计算其到每一个簇的加权距离；

其中，L(P，O_i)表示像素点P和簇i中心像素点O_i之间的加权距离，N是需要分割的像素点的总数量，σ_i是第i个簇的簇内标准差，d为像素点和簇中心的的欧式距离；

接着将图像中的每个像素点划分到与其加权距离最小的簇中，然后重新计算各个簇的聚簇中心，新的聚簇中心为每个簇中所有对象灰度值的平均值，利用下式计算聚簇质量E：

其中，n_i为图像第i个簇中像素点的数量，N是需要分割的像素点的总数量,σ_i是第i个簇的簇内标准差，K表示聚簇中心的数量；

聚簇质量达到了期望值或者达到了预先设定的最大迭代次数，则停止迭代；否则重新迭代聚簇过程；最终根据最后的聚簇结果，将同一个簇中的对象用相同的颜色标记，不同簇中的对象用不同的颜色标记，输出得到分割后的二维CT图像。

优选地，所述步骤六的参数计算包括：

C_mid＝(C_f1+C_f2)/2

其中，C_f1、C_f2分别表示两个股骨头的空间中心，C_p表示骶骨面的空间模型中心，N_p表示骶骨面的空间法向量；C_mid为两个股骨头中心连接线的空间中心；PI_3D表示空间中的骨盆入射角、PT_3D表示空间中的骨盆倾斜角、SS_3D表示空间骶骨倾斜角，z为常量取值为(0,0,1)。

本发明的优点在于：其一，提供一种股骨头三维模型的重建方法，(1)其检测涉及到圆上边缘点的完全覆盖，提高了寻找真实圆的可靠性；(2)计算几何约束不复杂；(3)几何约束的性质可以防止背景噪声和无关纹理的误报。其二，提供一种空间骨盆参数测量方法，这种基于于改进的随机圆检测算法的空间骨盆参数测量方法，不仅提升了估算维度还去除了人工参与，大大提高了测量准确度，保证了脊柱矫形手术的顺利实施。

本发明能够自动地精确实现股骨头的分割定位及骶骨面的识别，与现有技术相比：

(1)将二维参数扩展到三维参数，可以拓宽椎体研究的范围，更符合临床实践；

(2)消除用户交互需求，保证框架更高效、可靠、准确，技术要求更低。

附图说明

图1为本发明实施例2中VGG16网络模型的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，当元件被称为“固定于”另一个元件，它可以直接在另一个元件上或者也可以存在居中的元件。当一个元件被认为是“连接”另一个元件，它可以是直接连接到另一个元件或者可能同时存在居中元件。

实施例1

本实施例公开一种股骨头三维模型的重建方法，包括以下步骤：

S1、将包含股骨头的原始二维CT图像进行二值化处理，得到含股骨头的二值图；

S2、将含股骨头的二值图采用非梯度最大值抑制算法获得边缘轮廓；

S3、利用高斯卷积核在对边缘轮廓图像在水平和垂直方向上进行卷积，得到每个像素的梯度向量；

Grad_x(x,y)＝G_x(x,y)*I_edge

Grad_y(x,y)＝G_y(x,y)*I_edge

其中，Grad_x(x，y)，Grad_x(x，y)分别表示每个像素在水平和垂直方向的梯度值，G_x(x，y)，G_y(x，y)为水平和垂直方向卷积核，卷积核采用sobel算子，I_edge为边缘轮廓图像，水平和垂直两个方向梯度构成一个梯度向量

Grad(x，y)＝(Grad_x(x，y)，Grad_x(x，y))

S4、获得的边缘轮廓中的所有像素均为三角形候选像素；计算每对像素的梯度模值，得到构成一个圆的像素对：

比较每一对像素的梯度模值大小和构成三角形时两个底角的角度大小；角度相同的被认为是圆上的一对点；由于在同一个圆的上的像素的梯度方向指向圆心，且半径相同；因此，圆上的任意一对像素和它们指向的中心位置都可能形成一个等腰三角形。同样，如果一对像素可以构成一个等腰三角形，那么这对像素就能表示一个圆。

边缘轮廓中两个像素P₁、P₂之间的连接线用表示，以为边构造等腰三角形，由于圆上像素的梯度向量都是指向圆心的，所以将归一化梯度向量用表示；构成三角形的两个底角θ₁、θ₂的计算方法如下：

相似性度量采用余弦形式，如下：

|cos(θ₁)-cos(θ₂)|≤ε

其中，ε是容错控制项，它的大小控制了等腰三角形识别的准确率，上式成立时，默认对应的P₁、P₂为同一个圆上的两个点，该圆的圆心为等腰三角形的顶点。优选地，通过选取相对位置距离较大的像素对，能有效抑制检测失误。

S5、对获得的所有圆心点进行聚类，包括以下步骤：

S51、初始化聚类簇，簇的数量设定为0；将S4得到的所有圆心点作为圆心点集合C；

S52、在圆心点集合C中选择一个圆心点并遍历集合中所有的圆心点，计算当前圆心点与其他圆心点之间的距离；

S53、将当前圆心点和所有使S52的距离小于设定阈值的圆心聚类到同一个簇中，本实施例选择的阈值为1，当然，本领域的普通技术人员根据实际情况选择其他的阈值也应该在本发明的保护范围内；

S54、去除已聚类圆心点并更新圆心点集合C；

S55、直至圆心点集合为空，则完成圆心点聚类；否则，转S52；

由于每个等腰三角形都对应一个圆心点，且聚类簇中的圆心距离都相差不大，因此上述步骤实际上是将处于同一个圆上的等腰三角形进行聚类；统计各个聚类簇中等腰三角形的边长即拥有最多边长的聚类簇即为股骨头区域。

S6、采用聚类算法对S5中识别出的股骨头区域进行聚类，得到三维股骨头坐标，重建股骨头三维模型；

所述重建股骨头三维模型方法包括以下步骤：

S61、初始化聚类簇，簇的数量设定为0；

S62、读取识别的S5中的股骨头区域，并遍历所有的簇的中心点，计算股骨头中心点与簇中心点的距离；

S63、若S62的距离小于设定的阈值，本实施例选择的阈值为3，当然，本领域的普通技术人员根据实际情况选择其他的阈值也应该在本发明的保护范围内。将读出的股骨头区域聚类到相应的簇中，簇内股骨头区域数量加1，更新簇内中心为所有股骨头区域中心的平均值，更新最大半径，更新具有最大半径股骨头区域的位置；

S64、若S63不成立，则新建聚类簇，重新依次执行S62、S63；

S65、直至股骨头区域队列为空，则转S66；否则，转S62；

S66、挑选拥有最多股骨头的聚类簇，则最大聚类簇的中心坐标作为三维股骨头中心点的X-Y平面坐标，最大半径所在的位置作为股骨头中心点Z轴坐标，在三维模型中空间股骨头位置绘制球体拟代股骨头，球体半径为在二维平面上识别出的最大半径。

本发明的优点在于：(1)其检测涉及到圆上边缘点的完全覆盖，提高了寻找真实圆的可靠性；(2)计算几何约束不复杂；(3)几何约束的性质可以防止背景噪声和无关纹理的误报。

实施例2

本实施例公开一种空间骨盆参数测量方法，包括以下步骤：

收集病人的原始二维CT图像，组成二维CT图像数据集，原始二维CT图像数据集包括来自于不同地区、不同年龄、不同性别的病人的原始二维CT图像。

步骤一、将原始二维CT图像进行分割处理，并重建三维骨盆模型；

将原始二维CT图像分割处理，为更好地去除杂质，在分割之后，优选对图像做中值滤波处理，然后利用MC算法重建三维骨盆模型。

本发明使用基于加权质量评价函数的K-MEANS算法进行二维CT图像分割处理，阈值选取为192。具体方法如下：

首先输入待分割的二维CT图像，之后进行灰度化，使用基于信息熵的迭代算法初始化K个聚簇中心。然后利用下式对图像中的每个像素点计算其到每一个簇的加权距离。

其中，L(P，O_i)表示像素点P和簇i中心像素点O_i之间的加权距离，N是需要分割的像素点的总数量，σ_i是第i个簇的簇内标准差，d为像素点和簇中心的的欧式距离。

接着将图像中的每个像素点划分到与其加权距离最小的簇中，然后重新计算各个簇的聚簇中心，新的聚簇中心为每个簇中所有对象灰度值的平均值，利用下式计算聚簇质量E：

其中，n_i为图像第i个簇中像素点的数量，N是需要分割的像素点的总数量,σ_i是第i个簇的簇内标准差，K表示聚簇中心的数量。

如果聚簇质量达到了期望值，本实施例选用的期望值为0.8或者达到了预先设定的最大迭代次数，本实施例的最大迭代次数为100次，则停止迭代；否则重新迭代聚簇过程。最终根据最后的聚簇结果，将同一个簇中的对象用相同的颜色标记，不同簇中的对象用不同的颜色标记，输出得到分割后的二维CT图像。此算法能够清晰锐化骨盆图像骨质部分边缘，传统的二值化只是分为两种颜色，本发明每个簇都有一个颜色标记，达到多标记的效果。

步骤二、重建股骨头三维模型，包括以下步骤：

S21、将包含股骨头的原始二维CT图像进行二值化处理，得到含股骨头的二值图；

S2、将含股骨头的二值图采用非梯度最大值抑制算法获得边缘轮廓；

S23、利用高斯卷积核在对边缘轮廓图像在水平和垂直方向上进行卷积，得到每个像素的梯度向量：

Grad_x(x,y)＝G_x(x,y)*I_edge

Grad_y(x,y)＝G_y(x,y)*I_edge

其中，Grad_x(x，y)，Grad_x(x，y)分别表示每个像素在水平和垂直方向的梯度值，G_x(x，y)，G_y(x，y)为水平和垂直方向卷积核，卷积核采用sobel算子，I_edge为边缘轮廓图像，水平和垂直两个方向梯度构成一个梯度向量

Grad(x，y)＝(Grad_x(x，y)，Grad_x(x，y))

S24、获得的边缘轮廓中的所有像素均为三角形候选像素；计算每对像素的梯度模值，具体如下：

比较每一对像素的梯度模值大小和构成三角形时两个底角的角度大小；角度相同的被认为是圆上的一对点；由于在同一个圆的上的像素的梯度方向指向圆心，且半径相同；因此，圆上的任意一对像素和它们指向的中心位置都可能形成一个等腰三角形。同样，如果一对像素可以构成一个等腰三角形，那么这对像素就能表示一个圆。

边缘轮廓中两个像素P₁、P₂之间的连接线用表示，以为边构造等腰三角形，由于圆上像素的梯度向量都是指向圆心的，所以将归一化梯度向量用表示；构成三角形的两个底角θ₁、θ₂的计算方法如下：

相似性度量采用余弦形式，如下：

|cos(θ₁)-cos(θ₂)|≤ε

其中，ε是容错控制项，它的大小控制了等腰三角形识别的准确率，上式成立时，默认对应的P₁、P₂为同一个圆上的两个点，该圆的圆心为等腰三角形的顶点。优选地，通过选取相对位置距离较大的像素对，能有效抑制检测失误。

S25、对获得的所有圆心点进行聚类，包括以下步骤：

S251、初始化聚类簇，簇的数量设定为0；将S24得到的所有圆心点作为圆心点集合C；

S252、在圆心点集合C中选择一个圆心点并遍历集合中所有的圆心点，计算当前圆心点与其他圆心点之间的距离；

S253、将当前圆心点和所有使S252的距离小于设定阈值的圆心聚类到同一个簇中，本实施例选择的阈值为1，当然，本领域的普通技术人员根据实际情况选择其他的阈值也应该在本发明的保护范围内；

S254、去除已聚类圆心点并更新圆心点集合C；

S255、直至圆心点集合为空，则完成圆心点聚类；否则，转S252；

由于每个圆心点都对应一个等腰三角形，且聚类簇中的圆心距离都相差不大，因此上述步骤实际上是将处于同一个圆上的等腰三角形进行聚类；统计各个聚类簇中等腰三角形的边长即拥有最多最长边长的聚类簇即为股骨头区域。

S26、采用聚类算法对S25中识别出的股骨头区域进行聚类，得到三维股骨头坐标，重建股骨头三维模型；

所述重建股骨头三维模型方法包括以下步骤：

S261、初始化聚类簇，簇的数量设定为0；

S262、读取识别的S25中的股骨头区域，并遍历所有的簇的中心点，计算股骨头中心点与簇中心点的距离；

S263、若S262的距离小于设定的阈值，本实施例选择的阈值为3，当然，本领域的普通技术人员根据实际情况选择其他的阈值也应该在本发明的保护范围内。将读出的股骨头区域聚类到相应的簇中，簇内股骨头区域数量加1，更新簇内中心为所有股骨头区域中心的平均值，更新最大半径，更新具有最大半径股骨头区域的位置；

S264、若S263不成立，则新建聚类簇，重新依次执行S262、S63；

S265、直至股骨头区域队列为空，则转S266；否则，转S262；

S266、挑选拥有最多股骨头的聚类簇，则最大聚类簇的中心坐标作为三维股骨头中心点的X-Y平面坐标，最大半径所在的位置作为股骨头中心点Z轴坐标，在三维模型中空间股骨头位置绘制球体拟代股骨头，球体半径为在二维平面上识别出的最大半径。

步骤三、建立VGG16网络模型，并对VGG16网络模型进行训练，直至其收敛。由于VGG16前面几层是卷积层的堆叠，后面几层是全连接层，最后是Softmax层。所有隐层的激活单元都是线性整流函数，同时VGG16使用多个较小卷积核的卷积层代替一个卷积核较大的卷积层，一方面可以减少参数，另一方面相当于是进行了更多的非线性映射，可以增加网络的拟合/表达能力。VGG16网络模型如图1所述：

VGG16模型参数表1如下，输入图像大小为512*512：

表1

本发明的VGG16网络模型采用以下步骤进行训练：VGG16的训练数据集是800张包含L5的原始二维CT图像、S1的原始二维CT图像和股骨头的原始二维CT图像，含L5的二维CT图像的标签为[1,0,0]、含S1的二维CT图像的标签为[0,1,0]、含股骨头的二维CT图像的标签为[0,0,1]。将带有标签的二维CT图像输入建立的VGG16网络模型中，进行训练，直至网络收敛。

本发明利用Adam算法调整VGG16网络模型的权值，直至DRINet网络模型收敛，判定收敛的条件为，收敛函数阈值为0.95；其中，一阶矩估计的指数衰减率β1设定为0.9，二阶矩估计的指数衰减率β2设定为0.999，防除零参数ε设定为1e-8。学习速率设定为1e-3。使用dice系数作为损失函数。

步骤四、将包含L5的原始二维CT图像、S1的原始二维CT图像和股骨头的原始二维CT图像输入至收敛后的步骤三中的VGG16网络模型中，按序识别输入的图像，根据图像类别预测骶骨面位置。

包含L5的原始二维CT图像、S1的原始二维CT图像和股骨头的原始二维CT图像输入至VGG16网络模型中，通过卷积层提取图像特征，图像特征经过全连接层会输出一个向量。向量中的第一个分量表示原始输入图像是含L5的二维CT图像的概率，第二个分量表示原始输入图像是含S1的二维CT图像的概率，第三个分量表示原始输入图像是含股骨头的二维CT图像的概率。由三个分量中最大值的位置确定图像类别，如第一个分量值最大，则此图像为含L5的二维CT图像；如第二个分量值最大，则此图像为含S1的二维CT图像；如第三个分量值最大，则此图像为含股骨头的二维CT图像。

由于骶骨面是在L5椎体下方的一个斜面。然而，L5在空间上遮挡了骶骨面一部分。因此，适合定位的切片图像包含S1和L5空间部分，确定的是L5和S1序列图像之间的连接点一定包含骶骨，本实施例选取L5图像序列的最后一张作为包含骶骨面的二维CT图像。

步骤五、对步骤四找到含有骶骨的二维CT图像运行4连通区域识别算法，根据最大连通区域找到用于映射的锚点，即最大连通区域的部分上边缘点；将所述锚点映射回重建的步骤一的三维骨盆模型，生成S1空间骶骨面模型；

所述步骤五中的4连通区域识别算法包括以下步骤：

S51从找到的含有骶骨的二维CT图像的第一行第一个像素点开始判断此点四邻域中的左方，上方的点的像素值是否为0，如果像素值都为0或不存在上方点和左方点，则此点表示一个新的区域的开始，并给它一个新的标记。

S52如果此点四邻域中的左方点像素值不为0，上方点像素值为0或上方点不存在，则标记此点为左方点的标记值；如果此点四邻域中的左方点像素值为0或左方点不存在，上方点像素值不为0，则标记此点为最上点的标记值。

S53如果此点四邻域中的左方点像素值不为0，上方点像素值不为0，则标记此点为这两个中的最小的标记值，并建立等价标记值对，即上方点像素值的标记和左方点像素值的标记表示了同一个连通域中的部分点。

S54逐行从左向右扫描含有骶骨的二维CT图像上的点，重复S52～S54。

S55根据等价标记值对和含有骶骨的二维CT图像中每个点的标记完成连通域的查找和连通域中点的数量的计算,找到包含点数量最多的连通域即为最大连通域。

所述步骤五中的所述锚点为最大连通域图像上边缘点，通过MC算法将锚点映射回重建的步骤一的三维骨盆模型；采用KD-tree算法在三维骨盆模型中搜索距离锚点最近的三维点；在这些三维点上运行三维连通区域算法，形成S1的空间骶骨模型，

其中，S1的中心和法向量是被测空间平面上所有点的平均位置和法向量。

所述采用KD-tree算法在三维骨盆模型中搜索距离锚点最近的三维点包括以下步骤：

S501构建KD-tree模型；

S5011将重建的三维骨盆数据集合按照X、Y、Z三个方向划分成三个子集合，对每个子集合计算方差，选择具有最大方差的子集合，然后在该子集合上选择中值m做为中心点，用该中心点对该三维骨盆数据集合进行划分，得到两个子集合；同时创建一个树结点，用于存储；

S5012对两个子集合重复S5011步骤的过程，直至所有子集合都不能再划分为止；如果某个子集合不能再划分时，则将该子集合中的数据保存到叶子结点；

S502将锚点Q从根结点开始，按照Q与各个结点的比较结果向下访问Kd-Tree模型，直至达到叶子结点；

其中Q与结点的比较指的是将Q对应于结点中的k维度上的值与m进行比较，若Q(k)<m，则访问左子树，否则访问右子树；达到叶子结点时，计算Q与叶子结点上保存的数据之间的距离，记录下最小距离对应的数据点，记为当前“最近邻点”Pcur和最小距离Dcur。

S503进行回溯操作，找到离Q更近的“最近邻点”；即判断未被访问过的分支里是否还有离Q更近的点，它们之间的距离小于Dcur；

S504如果Q与其父结点下的未被访问过的分支之间的距离小于Dcur，则说明该分支中存在离P更近的数据，进入该结点，进行S501的查找过程，如果找到更近的数据点，则更新为当前的“最近邻点”Pcur，并更新Dcur；

如果Q与其父结点下的未被访问过的分支之间的距离大于Dcur，则说明该分支内不存在与Q更近的点；

回溯的判断过程是从下往上进行的，直到回溯到根结点时已经不存在与P更近的分支为止。

步骤六、计算骨盆参数。取股骨头三维模型的空间中心和S1空间骶骨面模型上的参数进行骨盆参数计算。

C_mid＝(C_f1+C_f2)/2

其中，C_f1、C_f2分别表示两个股骨头的空间中心，C_p表示骶骨面的空间模型中心，N_p表示骶骨面的空间法向量；C_mid为两个股骨头中心连接线的空间中心；PI_3D表示空间中的骨盆入射角、PT_3D表示空间中的骨盆倾斜角、SS_3D表示空间骶骨倾斜角，z为常量取值为(0,0,1)。

本发明包括股骨头识别重建模块、骶骨面识别定位模块、骨盆参数计算模块。首先对原始图像进行阈值处理，得到二值化图像，利用二值化图像及MC算法实现骨盆的三维重建。其次改进的随机圆检测算法在二维原始图像上检测股骨头区域，并保存这些区域。算法对保留区域列表中的区域中心进行聚类，并保持最大的聚类。然后以最大簇的中心为空间股骨头中心，以最大半径为空间半径。然后利用改进的随机圆检测算法找到包含骶骨面的图像，利用连通区域算法，找到最大连通区域上边缘上的点，将其映射到重建好的骨盆图像三维空间中。利用kdtree算法及寻找最邻近点方式实现骶骨面的识别，最终计算三维空间骨盆参数。

本发明提供了一种基于改进的随机圆检测算法的空间骨盆参数测量方法，不仅提升了估算维度还去除了人工参与，大大提高了测量准确度，保证了脊柱矫形手术的顺利实施。

本发明能够自动地精确实现股骨头的分割定位及骶骨面的识别，与现有技术相比：

(1)将二维参数扩展到三维参数，可以拓宽椎体研究的范围，更符合临床实践；

(2)消除用户交互需求，保证框架更高效、可靠、准确，技术要求更低。

需要说明的是，在本文中，如若存在第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

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

Claims (10)

1.一种股骨头三维模型的重建方法，其特征在于，包括以下步骤：

S1、将包含股骨头的原始二维CT图像进行二值化处理，得到含股骨头的二值图；

S2、将含股骨头的二值图采用非梯度最大值抑制算法获得边缘轮廓；

S3、利用高斯卷积核在对边缘轮廓图像在水平和垂直方向上进行卷积，得到每个像素的梯度向量；

S4、获得的边缘轮廓中的所有像素均为三角形候选像素；计算每对像素的梯度模值，得到构成一个圆的像素对；

S5、对步骤S4中每对像素构成的圆的圆心点进行聚类，识别出股骨头区域；

S6、采用聚类算法对S5中识别出的股骨头区域进行聚类，得到三维股骨头坐标，重建股骨头三维模型。

2.根据权利要求1所述肚饿一种股骨头三维模型的重建方法，其特征在于，所述S3采用包括以下方式得到梯度向量：

Grad_x(x,y)＝G_x(x,y)*I_edge

Grad_y(x,y)＝G_y(x,y)*I_edge

其中，Grad_x(x，y)，Grad_x(x，y)分别表示每个像素在水平和垂直方向的梯度值，G_x(x，y)，G_y(x，y)为水平和垂直方向卷积核，卷积核采用sobel算子，I_edge为边缘轮廓图像，水平和垂直两个方向梯度构成一个梯度向量

Grad(x，y)＝(Grad_x(x，y)，Grad_x(x，y))。

3.根据权利要求1所述肚饿一种股骨头三维模型的重建方法，其特征在于，所述S4采用包括以下方式获得构成一个圆的像素对：

比较每一对像素的梯度模值|Grad(x，y)|的大小和构成三角形时两个底角的角度大小；角度相同的被认为是圆上的一对点；该对像素能构成一个圆；

边缘轮廓中两个像素P₁、P₂之间的连接线用表示，以为边构造等腰三角形，圆上像素的梯度向量都是指向圆心的，将归一化梯度向量用表示；构成三角形的两个底角θ₁、θ₂的计算方法如下：

4.根据权利要求1所述的一种股骨头三维模型的重建方法，其特征在于，所述S4采用相似性度量构成一个圆的一对像素的准确率；

相似性度量采用余弦形式，如下：

|cos(θ₁)-cos(θ₂)|≤ε

其中，ε是容错控制项，它的大小控制了等腰三角形识别的准确率，上式成立时，默认对应的P₁、P₂为同一个圆上的两个点，该圆的圆心为等腰三角形的顶点。

5.根据权利要求1所述的一种股骨头三维模型的重建方法，其特征在于，所述S5包括以下步骤：

S51、初始化聚类簇，簇的数量设定为0；将S4得到的所有圆心点作为圆心点集合C；

S52、在圆心点集合C中选择一个圆心点并遍历集合中所有的圆心点，计算当前圆心点与其他圆心点之间的距离；

S53、将当前圆心点和所有使S52的距离小于设定阈值的圆心聚类到同一个簇中；

S54、去除已聚类圆心点并更新圆心点集合C；

S55、直至圆心点集合为空，则完成圆心点聚类；否则，转S52；

统计各个聚类簇中的边长，拥有最多最长边长的聚类簇即为股骨头区域。

6.根据权利要求1所述的一种股骨头三维模型的重建方法，其特征在于，所述S6包括以下步骤：

S61、初始化聚类簇，簇的数量设定为0；

S62、读取识别的S5中的股骨头区域，并遍历所有的簇的中心点，计算股骨头中心点与簇中心点的距离；

S63、若S62的距离小于设定的阈值，将读出的股骨头区域聚类到相应的簇中，簇内股骨头区域数量加1，更新簇内中心为所有股骨头区域中心的平均值，更新最大半径，更新具有最大半径股骨头区域的位置；

S64、若S63不成立，则新建聚类簇，重新依次执行S62、S63；

S65、直至股骨头区域队列为空，则转S66；否则，转S62；

S66、挑选拥有最多股骨头的聚类簇，则最大聚类簇的中心坐标作为三维股骨头中心点的X-Y平面坐标，最大半径所在的位置作为股骨头中心点Z轴坐标，在三维模型中空间股骨头位置绘制球体拟代股骨头，球体半径为在二维平面上识别出的最大半径。

7.一种基于权利要求1-6任一项所述的股骨头三维模型的重建方法的空间骨盆参数测量方法，其特征在于，包括以下步骤：

步骤一、将原始二维CT图像进行分割处理，并重建三维骨盆模型；

步骤二、重建股骨头三维模型；

步骤三、建立VGG16网络模型，并对VGG16网络模型进行训练，直至其收敛；

步骤四、将包含L5的原始二维CT图像、S1的原始二维CT图像和股骨头的原始二维CT图像输入至收敛后的步骤三中的VGG16网络模型中，按序识别输入的图像，根据图像类别预测骶骨面位置；

步骤五、对步骤四找到含有骶骨的二维CT图像运行4连通区域识别算法，根据最大连通区域找到用于映射的锚点，即最大连通区域的部分上边缘点；将所述锚点映射回重建的步骤一的三维骨盆模型，生成S1空间骶骨面模型；

步骤六、计算骨盆参数，取股骨头三维模型的空间中心的参数和S1空间骶骨面模型上的参数进行骨盆参数的计算。

8.根据权利要求7所述的空间骨盆参数测量方法，其特征在于，在所述步骤一中，采用基于加权质量评价函数的K-MEANS算法进行二维CT图像分割处理，利用MC算法重建三维骨盆模型。

9.根据权利要求8所述的空间骨盆参数测量方法，其特征在于，所述步骤一采用基于加权质量评价函数的K-MEANS算法进行二维CT图像分割处理，具体方法如下：

首先输入待分割的二维CT图像，之后进行灰度化，使用基于信息熵的迭代算法初始化K个聚簇中心；然后利用下式对图像中的每个像素点计算其到每一个簇的加权距离；

其中，L(P，O_i)表示像素点P和簇i中心像素点O_i之间的加权距离，N是需要分割的像素点的总数量，σ_i是第i个簇的簇内标准差，d为像素点和簇中心的的欧式距离；

接着将图像中的每个像素点划分到与其加权距离最小的簇中，然后重新计算各个簇的聚簇中心，新的聚簇中心为每个簇中所有对象灰度值的平均值，利用下式计算聚簇质量E：

其中，n_i为图像第i个簇中像素点的数量，N是需要分割的像素点的总数量,σ_i是第i个簇的簇内标准差，K表示聚簇中心的数量；

聚簇质量达到了期望值或者达到了预先设定的最大迭代次数，则停止迭代；否则重新迭代聚簇过程；最终根据最后的聚簇结果，将同一个簇中的对象用相同的颜色标记，不同簇中的对象用不同的颜色标记，输出得到分割后的二维CT图像。

10.根据权利要求7所述的空间骨盆参数测量方法，其特征在于，所述步骤六的参数计算包括：

C_mid＝(C_f1+C_f2)/2

其中，C_f1、C_f2分别表示两个股骨头的空间中心，C_p表示骶骨面的空间模型中心，N_p表示骶骨面的空间法向量；C_mid为两个股骨头中心连接线的空间中心；PI_3D表示空间中的骨盆入射角、PT_3D表示空间中的骨盆倾斜角、SS_3D表示空间骶骨倾斜角，z为常量取值为(0,0,1)。