基于CT图像的髋关节三维建模的自动定量计算方法

本发明属于计算机图形学与临床医学学科交叉领域，具体涉及基于ct图像的髋关节三维模型建立，以根据髋关节临床指标自动定量计算的方法。

背景技术：

1、髋关节撞击综合征(femoroacet abular impingement syndrome，fai)，又称股骨髋臼撞击综合征，是一种常见的、可引起慢性髋关节疼痛及髋关节屈曲和内收受限的疾病，是髋关节骨性关节炎常见的病因。髋关节撞击综合征是指股骨头颈交界处与髋臼的频繁撞击所导致的髋关节骨关节病，会使患者的髋关节出现慢性疼痛的症状，髋关节活动功能会随之降低，病情严重的情况下可以致残。

2、现阶段用于诊断髋关节撞击综合征的影像学检查方式主要有x光片(x-ray film，也称x线片)检查、计算机断层扫描(computed tomography，ct)检查和磁共振(magneticresonance imaging，mri)检查。

3、x光片检查的益处是操作简单方便、直观，多用于早期筛查，但仅是平面成像，测量的参数有限，无法观察髋关节的细微骨质结构；ct检查的特点是分辨率比较高，能够让微细软骨的显像非常清晰，弥补x射线诊断的一些缺漏，但仍无法观察髋关节软骨、髋盂唇等病理改变；mri具有多平面成像能力和在不同软组织间良好的对比度，有很高的特异性和敏感性，可显示x线片和ct均无法显示的软组织改变等。

4、传统的三维重建算法按传感器是否主动向物体照射光源分为主动式和被动式两种方法，主要包括以下几类：第一类是结构光算法，依靠投影仪将编码的结构光投射到被拍摄物体上，由摄像头进行拍摄。由于被拍摄物体上的不同部分相对于相机的距离精度和方向不同，结构光编码的图案的大小和形状也会发生改变。这种变化可以被摄像头捕获，然后通过运算单元将其换算成深度信息，进而获取物体的三维轮廓信息。缺点是容易受环境光干扰，随检测距离增加，其精度也会变差；第二类是激光飞行时间法(time of flight，tof)，通过向目标连续发送光脉冲，然后依据传感器接收到返回光的时间或相位差来计算距离目标的距离，这种方式足够的精度需要极为精确的时间测量模块，因此成本相对较高；第三类是三角测距法，依据三角测距原理，不需要较为精密的传感器，三角测距法整体成本较低，在近距离的时候精度较高。但三角测距的测量误差与距离有关，随着测量距离越来越大，测量误差也越来越大，这是由三角测量的原理导致的，不可避免；第四类是单目视觉算法，利用左右相机得到的两幅校正图像找到左右图片的匹配点，根据几何原理恢复出环境的三维信息。该方法难点在于左右相机图片的匹配，匹配不精确会影响最后算法成像的效果。多目视觉采用三个或三个以上摄像机来提高匹配的精度，需要消耗更多的时间，实时性也更差。

5、在现有技术中，中国发明专利申请号cn201610166101.5公开了一种基于3d打印技术的术前髋关节畸形骨骼模型制作方法，其特征在于，包括如下步骤：(1)对患者髋关节进行x线全长片力线定位，并采集患者髋关节的mri检查所生成的成像数据；(2)将步骤(1)所述的成像数据输入三维建模软件，通过图像分割、编辑、三维计算处理，提取得到患者原生态的髋关节畸形骨骼数字化三维模型；(3)将所述原生态的髋关节畸形骨骼数字化三维模型输入计算机辅助设计软件进行骨骼分层和定位处理，分别获得需要摘除的畸形骨骼模型图、位于所述畸形骨骼上方的上截骨模型图以及位于所述畸形骨骼下方的下截骨模型图，进一步提取所述的畸形骨骼模型图、上截骨模型图、下截骨模型图的成像数据；(4)将步骤(3)所述的畸形骨骼模型图、上截骨模型图、下截骨模型图的成像数据分别输入三维建模软件，之后利用3d打印设备进行打印，分别得到畸形骨骼模型、上截骨模型和下截骨模型；(5)将步骤(4)所述的上截骨模型、畸形骨骼模型、下截骨模型按照由上至下的顺序进行匹配结合，即得所述的术前髋关节畸形骨骼模型。

6、又例如，中国发明专利申请号cn201610819251.1公开了一种基于三维重建技术的髋臼侧模型和导板的制备方法，步骤一、将患者术前的髋关节三维ct数据进行加工后建立数字化三维数据模型，重建髋臼侧三维骨性结构，还原疾病原始状态；步骤二、通过分析髋臼位置、评估髋臼状态、明确髋臼周围骨量和计算髋臼周围骨厚度，制定髋臼磨锉和髋臼螺钉置入方案，并对这些导板和模型进行三维建模；步骤三、使用3d打印机制作术前髋臼侧骨性模型、术中磨臼导板、术中髋臼螺钉导板和术后髋臼侧骨性模型；步骤四、上述导板和模型整合并检验，若不合格则重新设计，若合格则将导板和模型消毒后投入手术。

7、再例如，中国发明专利申请号cn202210362422.8涉及医疗器械技术领域，具体涉及一种全髋关节三维重建系统及方法，包括手术台、扫描结构和支撑结构，所述支撑结构包括第一电缸、支撑环、驱动丝杆、滑动架、第一卷线单元、第二卷线单元、第一套环和第二套环，通过第一电缸调节支撑环的位置，将患病的大腿部分与支撑环卡合，两处脚踝分别与第一套环和第二套环卡合，对病人的双腿进行固定，便于扫描，手术过程中，需要更换腿部姿势，通过驱动丝杆调节滑动架的位置，第二卷线单元和第一卷线单元启动，将另一脚弯折，使得大腿部分向外侧翻动，暴露全髋关节的连接部分，便于后续的锉骨，实现腿部的固定和姿势变换，减轻了医护人员的工作量。

8、以上公开的发明专利申请与本技术属于同一技术领域，但均存在不能够或不能完全解决对髋关节骨性结构和周围软组织结构进行精确、直观的三维重现，以直观地显示髋臼和股骨头的空间关系，全面、精确髋臼双柱、臼顶和关节腔内情况。

技术实现思路

1、为了解决上述问题，本发明提供一种基于ct图像的髋关节三维建模的自动定量计算方法，具体步骤如下：

2、步骤1，对髋关节ct图像进行预处理，包括滤波、二值化、去枕床操作；

3、步骤2，对预处理后的髋关节ct图像基于轮廓线进行分割，得到盆骨点云pk，和股骨点云pm；

4、步骤3，对髋关节的点云数据采用主成分分析方法(principal componentanalysis，pca)重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型；

5、步骤4，髋关节三维建模处理；

6、步骤5，提取股骨头部点云pn通过球面拟合算法，得到股骨颈中心chead与半径rhead；

7、步骤6，基于单层球面截取股骨颈部轮廓线方法提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心cneck；

8、步骤7，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线；

9、步骤8，提取盆骨髋臼外缘线上的点云，拟合成平面s1，用s1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线；

10、步骤9，计算并输出α角、ce角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标；

11、进一步的，步骤1中，对髋关节ct图像进行预处理，包括滤波、二值化、去枕床操作：

12、在进行图像分割之前，需要先对原始ct图像进行去噪、二值化、去枕床的预处理操作，先使用中值滤波消除图像的部分噪声，再使用大津算法对灰度图像实现二值化，鉴于盆骨和股骨的点云数据在图像的上半部分而枕床的点云靠近下方，故根据点云的这一分布特征去除枕床部分的点云。

13、进一步的，步骤2中，对预处理后的髋关节ct图像基于轮廓线进行分割，得到盆骨点云pk，和股骨点云pm：

14、采用基于边界的方法提取图像轮廓线，分为像素梯度计算、非极大值抑制、滞后阈值处理以及孤立弱边缘抑制，在进行图像轮廓线的提取过程中，使用sobel算子计算像素梯度：sobel算子为两个3×3的矩阵sx、sy，前者用于计算图像x方向像素梯度矩阵gx，后者用于计算图像y方向像素梯度矩阵gy，如下式(1)、(2)：

15、

16、

17、其中，i为ct灰度图像矩阵，按照下式(3)：

18、

19、计算得到梯度强度矩阵gxy，通过梯度强度矩阵gxy，对当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大即为极值，则保留该像素为边缘点，若非最大，则对其进行抑制，不将其作为边缘点，之后定义一个高阈值和一个低阈值，梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理，在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，本算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合，至此，得到一个最为接近图像真实边缘的边缘图像，最终根据盆骨和股骨提取出的轮廓线长度、位置和形态差异，分割得到盆骨与股骨单独的点云pk和pm，其中，k为盆骨点云个数，m为股骨点云个数。

20、进一步的，步骤3中，对髋关节的点云数据采用pca重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型：

21、步骤3.1，基于pca对髋关节点云进行降维，对于从ct图像中提取的轮廓线，我们通过间隔均匀采样获取初始髋关节点云，为了获得更加光滑的表面，使用常用的数据降维技术pca对初始髋关节点云进行重采样，通过找到数据中主要方向来降低数据维度，从而实现对点云的重采样，使用pca能够在髋关节点云保持总体形状不变的同时去除点云中的噪声，且保持髋关节点云上均匀的密度分布，这为下游的重建与定量计算任务奠定了良好的基础，提高了后续髋关节点云处理的效率；

22、使用pca对髋关节点云进行降维的过程：在三维点云中，髋关节各部位的点云集合表示为p＝{pk}，pk∈rn×m，k∈z*，记为pk中第i行、第j列的值，pk中第j列的均值记为则pk去除中心化之后的点云记为p_newk，且即对髋关节点云中的每个点进行坐标变换，使得点云的重心在原点，之后求解重心在原点的髋关节点云数据集p_new的协方差矩阵，将其记为cov(p_new)，并对协方差矩阵进行特征分解，得到cov(p_new)的特征值和特征向量，将特征值按照从大到小进行排序，选择其中最大的k个作为髋关节点云的主成分，将其对应的k个特征向量分别作为列向量组成特征向量矩阵wn×k，最后计算即将重心在原点的髋关节点云数据集p_new投影到选取的特征向量上，重新计算髋关节点云的坐标来实现点云的降维，获得需要的已经降维的髋关节点云数据集

23、步骤3.2：无定向髋关节点云的迭代泊松重建：

24、由于髋关节各处的点法向量变化大，不便于准确计算，而传统的泊松重建方法又需要点的法向量作为输入，在髋关节建模部分选用的是基于无定向点的迭代泊松重建方法，这种方法不需要点的法向量作为输入，而是通过不断迭代来获得各点更加精确的法向量信息，当分配随机法线输入点，泊松重建生成一个表面，在每次迭代中，算法以直接从前一次迭代得到的曲面计算法线样本作为输入点，然后建立质量较好的新髋关节曲面。

25、进一步的，步骤4中，髋关节三维建模处理步骤如下：

26、给定一组髋关节无方向点pu＝{p1,p2,...,pu}，其中u为无方向点的个数，最大八叉树深度h，收敛阈值δ和筛选的psr权重γ作为输入，首先构造一个深度最大为h的八叉树，并使用八叉树节点作为样本集o，构造一个kd树以方便样本搜索，对每个样本oi∈o分配一个随机的法向量，之后在每次迭代中，将筛选出的psr应用于具有当前法线的样本集o，得到一个指标函数x，应用面向八叉树的marching cube算法提取具有等值的等面f，其中，n＝|o|为样本数，对于等面的每一个三角形面，找到最接近它的k个样本，建立三角形面和样本之间的关系，在迭代计算中通过计算与oi相关联的面的法线的平均值来更新每个样本的法线最后计算当前迭代和之前迭代的样本之间的法线方差，当法线方差的平均值小于指定的阈值δ时，算法终止，通过迭代泊松重建方法能够得到更加精确的髋关节重建模型。

27、进一步的，步骤5中，提取股骨头部点云pn，通过球面拟合算法得到股骨颈中心chead与半径rhead：

28、根据股骨头在整个股骨的占比选取三个点确定一个平面，用平面截取属于股骨头上的点云数据集合pn，用球面拟合算法对点云集合pn进行拟合运算，假设球心为(x0,y0,z0),半径为r，(x,y,z)为点云集合pn中的点：

29、首先建立一个球面方程：

30、(x-x0)2+(y-y0)2+(z-z0)2＝r2 ⑷，

31、建立误差函数：

32、

33、误差函数e表示所有目标点到拟合球心的误差平方和，如令：

34、s＝x2+y2+z2-r2 ⑹，

35、则函数e又能够记作：

36、

37、当误差函数e为最小值时，球心(x0,y0,z0)是e的最优估计值，e分别对参数x0,y0,z0,s取偏导数，如下式(8)：

38、

39、通过矩阵进行求解，从而求得x0,y0,z0,s的值，就能够确定出半径rhead的值，并根据一组在球面上的点拟合出了股骨头球体，得到了股骨头中心chead和股骨头半径rhead。

40、进一步的，步骤6中：基于单层球面截取股骨颈部轮廓线方法，提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心cneck：

41、建立球心在股骨头中心chead(x0,y0,z0)，半径为rhead+ε的球面方程：(x-x0)2+(y-y0)2+(z-z0)2-(rhead+ε)2＝0，其中，ε是一个微小偏移量，ε∈[3mm,5mm]，将股骨点云数据代入球面方程，得到球面与股骨颈的相交轮廓线附近的点云数据，将得到的球面与股骨颈的相交轮廓线附近的点云数据拟合成三维空间中的椭圆，取椭圆的中心cneck为密度中心，连接chead和cneck得到股骨颈轴线，记方向向量为n_neck。

42、进一步的，步骤7中，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线：

43、在通过高斯曲率寻找股骨头颈相交线时，首先寻找每个点的邻域点，固定距离球状邻域搜索首先设定距离阈值r，通过计算当前点与待定点之间的三维空间距离，将距离小于预设半径的所有待定点视为当前点的球状邻域点，设点p及其周围邻域内m个最近邻点为qi(i＝1,2,3,.,m)，设n为点p的法向量，mi为点qi的法向量,则p点相对于qi的法曲率kni估计如下式(9)：

44、

45、式中:为向量n与向量pqi之间的夹角；β为向量n与向量mi之间的夹角；|pqi|为p与qi之间的欧氏距离；

46、根据欧拉公式，法曲率与主曲率有以下关系，如下式(10):

47、

48、其中：i＝1,2,...,m；k1和k2是主曲率；θi+θ为点p过qi的法截线的切线与主方向的夹角，

49、令最终求主曲率的问题转化成求解下列方程最小二乘，如下式(11)：

50、

51、求解最小二乘以获得主曲率k1与k2，而k1×k2即为高斯曲率，得到每个点的高斯曲率值后，通过限定值的大小来找出股骨头头颈相交线附近的点云数据，再通过选出的点云数据拟合成一个空间椭圆，得到股骨头颈相交线。

52、进一步的，步骤8中，提取盆骨髋臼外缘线上的点云，拟合成平面s1，用s1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线：

53、在求髋臼外缘线时，首先求出股骨头和股骨颈的半径r1、r2和球心c1、c2以及盆骨的质心o，令股骨头上的每一点p与股骨头球心c1构成一条射线l1，分别计算盆骨上的点q与l1的距离，若距离小于设定的阈值d，则将该点纳入盆骨球冠点云中，令盆骨球冠上的点与股骨颈中心c2的连线记作l2，盆骨质心与股骨颈中心的连线记作l3，计算l2与l3形成的角度，从中选取最大角度范围内的点构成盆骨外缘线，拟合求得盆骨外缘线所在的平面方程z＝ax+by+c，利用平面截取股骨头，该平面与股骨头相交的点即构成髋臼外缘线，至此，得到髋臼外缘线的点云，通过选出的点云数据拟合成一个空间椭圆，得到髋臼外缘线。

54、进一步的，步骤9中，计算α角、ce角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标：

55、步骤9.1，计算α角：α角是股骨头球心chead和股骨头颈相交线上t点钟方向点t1的连线与股骨头颈轴线n_neck的夹角，连接chead、t1为向量v1，则α角可记为：

56、θα＝arccos((v1·n_neck)/(|v1||n_neck|))；

57、步骤9.2，计算股骨头凸出指数：股骨头凸出指数e是指股骨头超出髋臼覆盖的百分比，b点为股骨头颈相交轮廓线t点钟位置，c点为髋臼外缘轮廓线上t点钟位置，a点为b点所在直径与股骨头另一交点，计算a、b、c所在平面s2，记a、b、c在z轴方向的向量分别为v2,v3,v4，分别计算v2与v3，v2与v4的距离在平面s2上的投影记为d1、d2，则e＝(d1-d2d/)1；

58、步骤9.3计算股骨头颈偏心距：股骨头颈偏心距即测量平行的股骨头前缘最大半径的切线与邻近股骨颈切线的距离，股骨头中心为chead、股骨颈中心为cneck，过chead点、cneck点所在直线为a，即n_neck所在的直线，t2点为股骨头与股骨颈相交轮廓线t点钟位置，过t2点作直线b平行于直线a，在直线a与直线b所在平面上，过chead点作直线a垂线与股骨头交于两点，其中，最靠近z正方向的点为m点，过m点与直线a的平行线记为直线c，连接m与t2为向量v5，股骨头颈偏心距为直线b与直线c的距离，则股骨头颈偏心距：d＝|v5||n_neck|sinβ/|n_neck|，β为向量v5和n_neck的夹角；

59、步骤9.4计算ce角：股骨头的中心为chead，e点为髋臼最外侧处，经过chead、e两点的直线与身体中线的平行线之间的夹角即为ce角，e为髋臼外缘线12点钟方向的点，股骨头中心chead在z轴上的向量为v6，连接e与chead记为向量v7，则ce角可记为：θce＝arccos((v6·v7)/(|v6||v7|))。

60、与现有技术相比较，本发明的优越效果在于：

61、1.本发明所述基于深度学习的三维重建方法对髋关节骨性结构和周围软组织结构均能进行精确、直观的三维重现，重建出的髋关节三维模型可以直观地显示髋臼和股骨头的空间关系，全面、精确髋臼双柱、臼顶和关节腔内情况，对髋关节结构的显示和测量更有优势。在得到的髋关节三维模型上进行髋关节三维定量计算以及刚性图像配准，以便于对于治疗髋关节撞击综合征进行客观的手术效果评价

62、2.本发明所述基于ct图像的髋关节三维建模的自动定量计算方法，通过对髋关节ct图像进行中值滤波、二值化等预处理，通过基于边界的方法提取轮廓线将股骨与盆骨分开，得到股骨与盆骨单独的点云。

63、3.本发明所述基于ct图像的髋关节三维建模的自动定量计算方法，采用基于无定向点的迭代泊松重建方法，对法向量进行迭代计算，得到更加准确的法向量，在重建之前对点云进行重采样，去除噪声点，建立了精确的髋关节三维模型。

64、4.本发明所述基于ct图像的髋关节三维建模的自动定量计算方法，基于单层球面截取颈部轮廓线获得股骨颈中心，以及使用点云高斯曲率获得股骨头颈相交线，对判断髋关节撞击综合征的评价指标从三维层面上进行了自动化的计算。