CN117078840A - 基于ct图像的髋关节三维建模的自动定量计算方法 - Google Patents
基于ct图像的髋关节三维建模的自动定量计算方法 Download PDFInfo
- Publication number
- CN117078840A CN117078840A CN202310832350.3A CN202310832350A CN117078840A CN 117078840 A CN117078840 A CN 117078840A CN 202310832350 A CN202310832350 A CN 202310832350A CN 117078840 A CN117078840 A CN 117078840A
- Authority
- CN
- China
- Prior art keywords
- point cloud
- point
- hip joint
- head
- neck
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 210000004394 hip joint Anatomy 0.000 title claims abstract description 142
- 238000004364 calculation method Methods 0.000 title claims abstract description 37
- 238000000034 method Methods 0.000 claims abstract description 55
- 210000000689 upper leg Anatomy 0.000 claims abstract description 52
- 210000000588 acetabulum Anatomy 0.000 claims abstract description 37
- 210000003049 pelvic bone Anatomy 0.000 claims abstract description 35
- 238000012952 Resampling Methods 0.000 claims abstract description 13
- 238000001914 filtration Methods 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims abstract description 10
- 238000007781 pre-processing Methods 0.000 claims abstract description 5
- 238000012847 principal component analysis method Methods 0.000 claims abstract description 4
- 239000013598 vector Substances 0.000 claims description 72
- 210000002436 femur neck Anatomy 0.000 claims description 48
- 239000011159 matrix material Substances 0.000 claims description 30
- 210000002391 femur head Anatomy 0.000 claims description 20
- 238000011156 evaluation Methods 0.000 claims description 14
- 230000009467 reduction Effects 0.000 claims description 12
- 206010039227 Rotator cuff syndrome Diseases 0.000 claims description 9
- 238000005520 cutting process Methods 0.000 claims description 9
- 230000005484 gravity Effects 0.000 claims description 9
- 210000004197 pelvis Anatomy 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 230000001629 suppression Effects 0.000 claims description 7
- 238000009826 distribution Methods 0.000 claims description 6
- 239000002356 single layer Substances 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 claims description 4
- 238000003709 image segmentation Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 241000287196 Asthenes Species 0.000 claims description 3
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 230000000877 morphologic effect Effects 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 2
- 210000000988 bone and bone Anatomy 0.000 abstract description 22
- 210000001624 hip Anatomy 0.000 abstract description 5
- 210000004872 soft tissue Anatomy 0.000 abstract description 5
- 230000006870 function Effects 0.000 description 11
- 238000003384 imaging method Methods 0.000 description 8
- 230000002093 peripheral effect Effects 0.000 description 7
- 208000011580 syndromic disease Diseases 0.000 description 6
- 238000002595 magnetic resonance imaging Methods 0.000 description 5
- 210000002414 leg Anatomy 0.000 description 4
- 238000005259 measurement Methods 0.000 description 4
- 238000004804 winding Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000011218 segmentation Effects 0.000 description 3
- 230000000007 visual effect Effects 0.000 description 3
- 238000010146 3D printing Methods 0.000 description 2
- 210000000845 cartilage Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000002591 computed tomography Methods 0.000 description 2
- 238000013135 deep learning Methods 0.000 description 2
- 238000003745 diagnosis Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 230000002980 postoperative effect Effects 0.000 description 2
- 238000000513 principal component analysis Methods 0.000 description 2
- 208000000094 Chronic Pain Diseases 0.000 description 1
- 206010062337 Congenital joint malformation Diseases 0.000 description 1
- 208000004199 Femoracetabular Impingement Diseases 0.000 description 1
- 208000007353 Hip Osteoarthritis Diseases 0.000 description 1
- 206010057178 Osteoarthropathies Diseases 0.000 description 1
- 208000002193 Pain Diseases 0.000 description 1
- 210000003423 ankle Anatomy 0.000 description 1
- 230000001684 chronic effect Effects 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000249 desinfective effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000003814 drug Substances 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000000227 grinding Methods 0.000 description 1
- 208000012285 hip pain Diseases 0.000 description 1
- 230000036244 malformation Effects 0.000 description 1
- 238000003801 milling Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 231100000915 pathological change Toxicity 0.000 description 1
- 230000036285 pathological change Effects 0.000 description 1
- 238000007639 printing Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
- A61B6/02—Devices for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computerised tomographs
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
- A61B6/50—Clinical applications
- A61B6/505—Clinical applications involving diagnosis of bone
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5258—Devices using data or image processing specially adapted for radiation diagnosis involving detection or reduction of artifacts or noise
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/20—Image enhancement or restoration by the use of local operators
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10081—Computed x-ray tomography [CT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30008—Bone
Abstract
本发明涉及一种基于CT图像的髋关节三维建模的自动定量计算方法，所述方法包括步骤如下：步骤1，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作；步骤2，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk，和股骨点云Pm；步骤3，对髋关节的点云数据采用主成分分析方法重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型；步骤4，髋关节三维建模处理等步骤；本发明所述方法的优越效果在于，对髋关节骨性结构和周围软组织结构均能进行精确、直观的三维重现，重建出的髋关节三维模型能够直观地显示髋臼和股骨头的空间关系，全面、精确髋臼双柱、臼顶和关节腔内情况，对髋关节结构的显示和测量更有优势。
Description
技术领域
本发明属于计算机图形学与临床医学学科交叉领域，具体涉及基于CT图像的髋关节三维模型建立，以根据髋关节临床指标自动定量计算的方法。
背景技术
髋关节撞击综合征(Femoroacet abular impingement Syndrome，FAI)，又称股骨髋臼撞击综合征，是一种常见的、可引起慢性髋关节疼痛及髋关节屈曲和内收受限的疾病，是髋关节骨性关节炎常见的病因。髋关节撞击综合征是指股骨头颈交界处与髋臼的频繁撞击所导致的髋关节骨关节病，会使患者的髋关节出现慢性疼痛的症状，髋关节活动功能会随之降低，病情严重的情况下可以致残。
现阶段用于诊断髋关节撞击综合征的影像学检查方式主要有X光片(X-ray Film，也称X线片)检查、计算机断层扫描(Computed Tomography，CT)检查和磁共振(MagneticResonance Imaging，MRI)检查。
X光片检查的益处是操作简单方便、直观，多用于早期筛查，但仅是平面成像，测量的参数有限，无法观察髋关节的细微骨质结构；CT检查的特点是分辨率比较高，能够让微细软骨的显像非常清晰，弥补X射线诊断的一些缺漏，但仍无法观察髋关节软骨、髋盂唇等病理改变；MRI具有多平面成像能力和在不同软组织间良好的对比度，有很高的特异性和敏感性，可显示X线片和CT均无法显示的软组织改变等。
传统的三维重建算法按传感器是否主动向物体照射光源分为主动式和被动式两种方法，主要包括以下几类：第一类是结构光算法，依靠投影仪将编码的结构光投射到被拍摄物体上，由摄像头进行拍摄。由于被拍摄物体上的不同部分相对于相机的距离精度和方向不同，结构光编码的图案的大小和形状也会发生改变。这种变化可以被摄像头捕获，然后通过运算单元将其换算成深度信息，进而获取物体的三维轮廓信息。缺点是容易受环境光干扰，随检测距离增加，其精度也会变差；第二类是激光飞行时间法(Time of flight，TOF)，通过向目标连续发送光脉冲，然后依据传感器接收到返回光的时间或相位差来计算距离目标的距离，这种方式足够的精度需要极为精确的时间测量模块，因此成本相对较高；第三类是三角测距法，依据三角测距原理，不需要较为精密的传感器，三角测距法整体成本较低，在近距离的时候精度较高。但三角测距的测量误差与距离有关，随着测量距离越来越大，测量误差也越来越大，这是由三角测量的原理导致的，不可避免；第四类是单目视觉算法，利用左右相机得到的两幅校正图像找到左右图片的匹配点，根据几何原理恢复出环境的三维信息。该方法难点在于左右相机图片的匹配，匹配不精确会影响最后算法成像的效果。多目视觉采用三个或三个以上摄像机来提高匹配的精度，需要消耗更多的时间，实时性也更差。
在现有技术中，中国发明专利申请号CN201610166101.5公开了一种基于3D打印技术的术前髋关节畸形骨骼模型制作方法，其特征在于，包括如下步骤：(1)对患者髋关节进行X线全长片力线定位，并采集患者髋关节的MRI检查所生成的成像数据；(2)将步骤(1)所述的成像数据输入三维建模软件，通过图像分割、编辑、三维计算处理，提取得到患者原生态的髋关节畸形骨骼数字化三维模型；(3)将所述原生态的髋关节畸形骨骼数字化三维模型输入计算机辅助设计软件进行骨骼分层和定位处理，分别获得需要摘除的畸形骨骼模型图、位于所述畸形骨骼上方的上截骨模型图以及位于所述畸形骨骼下方的下截骨模型图，进一步提取所述的畸形骨骼模型图、上截骨模型图、下截骨模型图的成像数据；(4)将步骤(3)所述的畸形骨骼模型图、上截骨模型图、下截骨模型图的成像数据分别输入三维建模软件，之后利用3D打印设备进行打印，分别得到畸形骨骼模型、上截骨模型和下截骨模型；(5)将步骤(4)所述的上截骨模型、畸形骨骼模型、下截骨模型按照由上至下的顺序进行匹配结合，即得所述的术前髋关节畸形骨骼模型。
又例如，中国发明专利申请号CN201610819251.1公开了一种基于三维重建技术的髋臼侧模型和导板的制备方法，步骤一、将患者术前的髋关节三维CT数据进行加工后建立数字化三维数据模型，重建髋臼侧三维骨性结构，还原疾病原始状态；步骤二、通过分析髋臼位置、评估髋臼状态、明确髋臼周围骨量和计算髋臼周围骨厚度，制定髋臼磨锉和髋臼螺钉置入方案，并对这些导板和模型进行三维建模；步骤三、使用3D打印机制作术前髋臼侧骨性模型、术中磨臼导板、术中髋臼螺钉导板和术后髋臼侧骨性模型；步骤四、上述导板和模型整合并检验，若不合格则重新设计，若合格则将导板和模型消毒后投入手术。
再例如，中国发明专利申请号CN202210362422.8涉及医疗器械技术领域，具体涉及一种全髋关节三维重建系统及方法，包括手术台、扫描结构和支撑结构，所述支撑结构包括第一电缸、支撑环、驱动丝杆、滑动架、第一卷线单元、第二卷线单元、第一套环和第二套环，通过第一电缸调节支撑环的位置，将患病的大腿部分与支撑环卡合，两处脚踝分别与第一套环和第二套环卡合，对病人的双腿进行固定，便于扫描，手术过程中，需要更换腿部姿势，通过驱动丝杆调节滑动架的位置，第二卷线单元和第一卷线单元启动，将另一脚弯折，使得大腿部分向外侧翻动，暴露全髋关节的连接部分，便于后续的锉骨，实现腿部的固定和姿势变换，减轻了医护人员的工作量。
以上公开的发明专利申请与本申请属于同一技术领域，但均存在不能够或不能完全解决对髋关节骨性结构和周围软组织结构进行精确、直观的三维重现，以直观地显示髋臼和股骨头的空间关系，全面、精确髋臼双柱、臼顶和关节腔内情况。
发明内容
为了解决上述问题，本发明提供一种基于CT图像的髋关节三维建模的自动定量计算方法，具体步骤如下：
步骤1，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作；
步骤2，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk，和股骨点云Pm；
步骤3，对髋关节的点云数据采用主成分分析方法(Principal ComponentAnalysis，PCA)重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型；
步骤4，髋关节三维建模处理；
步骤5，提取股骨头部点云Pn通过球面拟合算法，得到股骨颈中心Chead与半径Rhead；
步骤6，基于单层球面截取股骨颈部轮廓线方法提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心Cneck；
步骤7，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线；
步骤8，提取盆骨髋臼外缘线上的点云，拟合成平面S1，用S1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线；
步骤9，计算并输出α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标；
进一步的，步骤1中，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作：
在进行图像分割之前，需要先对原始CT图像进行去噪、二值化、去枕床的预处理操作，先使用中值滤波消除图像的部分噪声，再使用大津算法对灰度图像实现二值化，鉴于盆骨和股骨的点云数据在图像的上半部分而枕床的点云靠近下方，故根据点云的这一分布特征去除枕床部分的点云。
进一步的，步骤2中，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk，和股骨点云Pm：
采用基于边界的方法提取图像轮廓线，分为像素梯度计算、非极大值抑制、滞后阈值处理以及孤立弱边缘抑制，在进行图像轮廓线的提取过程中，使用Sobel算子计算像素梯度：Sobel算子为两个3×3的矩阵Sx、Sy，前者用于计算图像x方向像素梯度矩阵Gx，后者用于计算图像y方向像素梯度矩阵Gy，如下式(1)、(2)：
其中，I为CT灰度图像矩阵，按照下式(3)：
计算得到梯度强度矩阵Gxy，通过梯度强度矩阵Gxy，对当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大即为极值，则保留该像素为边缘点，若非最大，则对其进行抑制，不将其作为边缘点，之后定义一个高阈值和一个低阈值，梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理，在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，本算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合，至此，得到一个最为接近图像真实边缘的边缘图像，最终根据盆骨和股骨提取出的轮廓线长度、位置和形态差异，分割得到盆骨与股骨单独的点云Pk和Pm，其中，k为盆骨点云个数，m为股骨点云个数。
进一步的，步骤3中，对髋关节的点云数据采用PCA重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型：
步骤3.1，基于PCA对髋关节点云进行降维，对于从CT图像中提取的轮廓线，我们通过间隔均匀采样获取初始髋关节点云，为了获得更加光滑的表面，使用常用的数据降维技术PCA对初始髋关节点云进行重采样，通过找到数据中主要方向来降低数据维度，从而实现对点云的重采样，使用PCA能够在髋关节点云保持总体形状不变的同时去除点云中的噪声，且保持髋关节点云上均匀的密度分布，这为下游的重建与定量计算任务奠定了良好的基础，提高了后续髋关节点云处理的效率；
使用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投影到选取的特征向量上，重新计算髋关节点云的坐标来实现点云的降维，获得需要的已经降维的髋关节点云数据集
步骤3.2：无定向髋关节点云的迭代泊松重建：
由于髋关节各处的点法向量变化大，不便于准确计算，而传统的泊松重建方法又需要点的法向量作为输入，在髋关节建模部分选用的是基于无定向点的迭代泊松重建方法，这种方法不需要点的法向量作为输入，而是通过不断迭代来获得各点更加精确的法向量信息，当分配随机法线输入点，泊松重建生成一个表面，在每次迭代中，算法以直接从前一次迭代得到的曲面计算法线样本作为输入点，然后建立质量较好的新髋关节曲面。
进一步的，步骤4中，髋关节三维建模处理步骤如下：
给定一组髋关节无方向点Pu＝{P1,P2,...,Pu}，其中u为无方向点的个数，最大八叉树深度H，收敛阈值δ和筛选的PSR权重γ作为输入，首先构造一个深度最大为H的八叉树，并使用八叉树节点作为样本集O，构造一个kd树以方便样本搜索，对每个样本oi∈O分配一个随机的法向量，之后在每次迭代中，将筛选出的PSR应用于具有当前法线的样本集O，得到一个指标函数x，应用面向八叉树的marching cube算法提取具有等值的等面F，其中，n＝|O|为样本数，对于等面的每一个三角形面，找到最接近它的k个样本，建立三角形面和样本之间的关系，在迭代计算中通过计算与oi相关联的面的法线的平均值来更新每个样本的法线/>最后计算当前迭代和之前迭代的样本之间的法线方差，当法线方差的平均值小于指定的阈值δ时，算法终止，通过迭代泊松重建方法能够得到更加精确的髋关节重建模型。
进一步的，步骤5中，提取股骨头部点云Pn，通过球面拟合算法得到股骨颈中心Chead与半径Rhead：
根据股骨头在整个股骨的占比选取三个点确定一个平面，用平面截取属于股骨头上的点云数据集合Pn，用球面拟合算法对点云集合Pn进行拟合运算，假设球心为(x0,y0,z0),半径为R，(x,y,z)为点云集合Pn中的点：
首先建立一个球面方程：
(x-x0)2+(y-y0)2+(z-z0)2＝R2 ⑷，
建立误差函数：
误差函数E表示所有目标点到拟合球心的误差平方和，如令：
s＝x2+y2+z2-R2 ⑹，
则函数E又能够记作：
当误差函数E为最小值时，球心(x0,y0,z0)是E的最优估计值，E分别对参数x0,y0,z0,s取偏导数，如下式(8)：
通过矩阵进行求解，从而求得x0,y0,z0,s的值，就能够确定出半径Rhead的值，并根据一组在球面上的点拟合出了股骨头球体，得到了股骨头中心Chead和股骨头半径Rhead。
进一步的，步骤6中：基于单层球面截取股骨颈部轮廓线方法，提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心Cneck：
建立球心在股骨头中心Chead(x0,y0,z0)，半径为Rhead+ε的球面方程：(x-x0)2+(y-y0)2+(z-z0)2-(Rhead+ε)2＝0，其中，ε是一个微小偏移量，ε∈[3mm,5mm]，将股骨点云数据代入球面方程，得到球面与股骨颈的相交轮廓线附近的点云数据，将得到的球面与股骨颈的相交轮廓线附近的点云数据拟合成三维空间中的椭圆，取椭圆的中心Cneck为密度中心，连接Chead和Cneck得到股骨颈轴线，记方向向量为n_neck。
进一步的，步骤7中，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线：
在通过高斯曲率寻找股骨头颈相交线时，首先寻找每个点的邻域点，固定距离球状邻域搜索首先设定距离阈值r，通过计算当前点与待定点之间的三维空间距离，将距离小于预设半径的所有待定点视为当前点的球状邻域点，设点p及其周围邻域内m个最近邻点为qi(i＝1,2,3,.,m)，设N为点p的法向量，Mi为点qi的法向量,则p点相对于qi的法曲率kn i估计如下式(9)：
式中:为向量N与向量pqi之间的夹角；β为向量N与向量Mi之间的夹角；|pqi|为p与qi之间的欧氏距离；
根据欧拉公式，法曲率与主曲率有以下关系，如下式(10):
其中：i＝1,2,...,m；k1和k2是主曲率；θi+θ为点p过qi的法截线的切线与主方向的夹角，
令最终求主曲率的问题转化成求解下列方程最小二乘，如下式(11)：
求解最小二乘以获得主曲率k1与k2，而k1×k2即为高斯曲率，得到每个点的高斯曲率值后，通过限定值的大小来找出股骨头头颈相交线附近的点云数据，再通过选出的点云数据拟合成一个空间椭圆，得到股骨头颈相交线。
进一步的，步骤8中，提取盆骨髋臼外缘线上的点云，拟合成平面S1，用S1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线：
在求髋臼外缘线时，首先求出股骨头和股骨颈的半径R1、R2和球心C1、C2以及盆骨的质心O，令股骨头上的每一点P与股骨头球心C1构成一条射线L1，分别计算盆骨上的点Q与L1的距离，若距离小于设定的阈值D，则将该点纳入盆骨球冠点云中，令盆骨球冠上的点与股骨颈中心C2的连线记作L2，盆骨质心与股骨颈中心的连线记作L3，计算L2与L3形成的角度，从中选取最大角度范围内的点构成盆骨外缘线，拟合求得盆骨外缘线所在的平面方程Z＝AX+BY+C，利用平面截取股骨头，该平面与股骨头相交的点即构成髋臼外缘线，至此，得到髋臼外缘线的点云，通过选出的点云数据拟合成一个空间椭圆，得到髋臼外缘线。
进一步的，步骤9中，计算α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标：
步骤9.1，计算α角：α角是股骨头球心Chead和股骨头颈相交线上t点钟方向点T1的连线与股骨头颈轴线n_neck的夹角，连接Chead、T1为向量v1，则α角可记为：
θα＝arccos((v1·n_neck)/(|v1||n_neck|))；
步骤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；
步骤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的夹角；
步骤9.4计算CE角：股骨头的中心为Chead，E点为髋臼最外侧处，经过Chead、E两点的直线与身体中线的平行线之间的夹角即为CE角，E为髋臼外缘线12点钟方向的点，股骨头中心Chead在Z轴上的向量为v6，连接E与Chead记为向量v7，则CE角可记为：θCE＝arccos((v6·v7)/(|v6||v7|))。
与现有技术相比较，本发明的优越效果在于：
1.本发明所述基于深度学习的三维重建方法对髋关节骨性结构和周围软组织结构均能进行精确、直观的三维重现，重建出的髋关节三维模型可以直观地显示髋臼和股骨头的空间关系，全面、精确髋臼双柱、臼顶和关节腔内情况，对髋关节结构的显示和测量更有优势。在得到的髋关节三维模型上进行髋关节三维定量计算以及刚性图像配准，以便于对于治疗髋关节撞击综合征进行客观的手术效果评价
2.本发明所述基于CT图像的髋关节三维建模的自动定量计算方法，通过对髋关节CT图像进行中值滤波、二值化等预处理，通过基于边界的方法提取轮廓线将股骨与盆骨分开，得到股骨与盆骨单独的点云。
3.本发明所述基于CT图像的髋关节三维建模的自动定量计算方法，采用基于无定向点的迭代泊松重建方法，对法向量进行迭代计算，得到更加准确的法向量，在重建之前对点云进行重采样，去除噪声点，建立了精确的髋关节三维模型。
4.本发明所述基于CT图像的髋关节三维建模的自动定量计算方法，基于单层球面截取颈部轮廓线获得股骨颈中心，以及使用点云高斯曲率获得股骨头颈相交线，对判断髋关节撞击综合征的评价指标从三维层面上进行了自动化的计算。
附图说明
图1为本发明所述基于CT图像的髋关节三维建模的自动定量计算方法示意图；
图2表示髋臼部位方向，其中A为股骨头正视图，B为盆骨正视图，C为股骨时钟方向图，D为髋臼时钟方向图；
图3示出股骨髋臼撞击综合征评价指标，其中A为α角，B为股骨头突出指数，C为股骨头颈偏心距，D为CE角。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。
在本实施例中，所述基于CT图像的髋关节三维建模的自动定量计算方法，基于髋关节的CT图像进行分割，将股骨和盆骨分开，提取点云数据，根据从髋关节CT图像提取的点云数据进行髋关节三维模型的建立以及对判断髋关节撞击综合征的四个评价指标从三维层面上进行自动定量计算。基于几何形状描述的评价指标，如α角、CE角、股骨头突出指数、股骨头颈偏心距对股骨髋臼撞击综合征的诊断、术前术后评价有着极其关键的作用。
在本实施例中，以下对本发明有关技术内容进行详细描述与定义：
髋臼部位采用时钟进行方向描述：
临床上将髋臼及股骨头某一侧的解剖面近似为一圆，并依此建立髋臼部位的统一方向描述，如图1所示，A)、B)所呈现的视角规定为人体的正向，则从左(右)侧观察股骨头和髋臼可获得C)、D)中视角，以该侧视角的解剖面可近似视为一圆，其中，D)中髋臼围(也称为髋臼外围轮廓线，该部位与之对应的股骨头上一条轮廓线为髋臼外围股骨头轮廓线，本实施例中不区分两者时统称为外围轮廓线)最高点为12点钟，并为了避免因为从左、右侧观测时3点钟方和9点钟方向不一致，特别规定始终以解剖面上靠近正侧的方向为3点钟方向，远离正侧的方向为9点钟方向，于此，基于髋臼围建立的时钟方向描述建立完成；类似的，C)给出了股骨头剖面上时钟方向描述示意图，对于左、右股骨，与基于髋臼围建立的时钟方向描述使用相同的规定，始终以靠近正侧的方向作为3点钟方向，远离正侧的方向作为9点钟方向。
髋臼部位在世界坐标系下的方向描述：
以CT图像中的世界坐标为参考笛卡尔坐标系，如图2所示，A)中人体直立方向为Z轴正方向，垂直于A)中视角所在平面并指向观察者方向为Y轴正方向，则Z＝X×Y，在给出世界坐标系下方向描述之前，给定以下两个假设前提：
1)外围轮廓线位于三维空间某一平面S上；
2)平面上的外围轮廓线为三维空间椭圆曲线l(θ)；
其中，0≤θ≤2π，令外围轮廓线上正12点位置0＝l(0)＝l(2π)，则给定点钟时外围轮廓线上位置为
需要说明的是，虽然无论从左、右任意侧观察时，始终保持靠近正侧端为3点钟方向，远离正侧端为9点钟方向，但在实际计算中仅为顺时针或逆时针旋转相同角度的关系，即取θ'＝-θ，更为方便的，此处描述时并不作区分。
以下，对FAI评价指标及其相关内容进行定义与说明：
对各FAI评价指标基于形态学的定义：
α角股骨头球心和股骨头非球面起始处的连线与股骨头颈轴线n_neck的夹角，如图3中A)所示，a为股骨头球心，b为股骨颈中心，c为股骨头颈相交线t点钟位置；
股骨头凸出指数：
系股骨头超出髋臼覆盖的百分比，如图3中B)所示，b点为股骨头与股骨颈相交轮廓线t点钟位置，c点为髋臼外缘轮廓线上t点钟位置，a点为b点所在直径与股骨头另一交点；
股骨头颈偏心距：
即测量平行的股骨头前缘最大半径的切线与邻近股骨颈切线的距离，如图3中C)所示，a点、b点分别为股骨头中心、股骨颈中心，过a点、b点所在直线为A，c点为股骨头与股骨颈相交轮廓线t点钟位置，过c点作直线B平行于直线A，在直线A与直线B所在平面上，过a点作直线A垂线与股骨头交于两点，其中最靠近Z正方向的点为d点，过d点与直线A的平行线记为直线C，则股骨头偏心距为直线B与直线C的距离；
CE角C点定义为股骨头的中心：
E点为髋臼最外侧处，经过C、E两点的直线与身体中线的平行线之间的夹角即为CE角，如图3中D)所示，a点为股骨头中心，a点与b点所在直线为Z轴正方向，c点为髋臼外缘轮廓线上t点钟位置，图3中示出的骨髋臼撞击综合征评价指标：A为α角；B为股骨头突出指数；C为股骨头颈偏心距；D为CE角。
下面结合说明书附图1至3，详细描述本发明所述基于CT图像的髋关节三维建模的自动定量计算方法的具体步骤：
步骤1，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作：
在进行图像分割之前，需要先对原始CT图像进行去噪、二值化、去枕床的预处理操作，先使用中值滤波消除图像的部分噪声，再使用大津算法对灰度图像实现二值化，鉴于盆骨和股骨的点云数据在图像的上半部分而枕床的点云靠近下方，故根据点云的这一分布特征去除枕床部分的点云；
步骤2，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk，和股骨点云Pm，其中，k为盆骨点云个数，m为股骨点云个数：
采用基于边界的方法提取图像轮廓线，分为像素梯度计算、非极大值抑制、滞后阈值处理以及孤立弱边缘抑制，在进行图像轮廓线的提取过程中，使用Sobel算子计算像素梯度：Sobel算子为两个3×3的矩阵Sx、Sy，前者用于计算图像x方向像素梯度矩阵Gx，后者用于计算图像y方向像素梯度矩阵Gy，具体形式如下式(1)、(2)：
其中，I为CT灰度图像矩阵，由下式(3)
计算得到梯度强度矩阵Gxy，通过梯度强度矩阵Gxy，对当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大(即为极值)，则保留该像素为边缘点，若不是最大，则对其进行抑制，不将其作为边缘点，之后定义一个高阈值和一个低阈值，梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理，在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，本算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合，至此，得到一个最为接近图像真实边缘的边缘图像，最终根据盆骨和股骨提取出的轮廓线长度、位置和形态差异，分割得到盆骨与股骨单独的点云Pk和Pm；基于计算所得的像素梯度，通过设定阈值划分，利用边缘抑制确定图像轮廓线，分割得到盆骨与股骨单独的点云Pk和Pm；
步骤3，对髋关节的点云数据采用主成分分析方法(Principal ComponentAnalysis，PCA)重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型：
步骤3.1，基于PCA对髋关节点云进行降维，对于从CT图像中提取的轮廓线，我们通过间隔均匀采样获取初始髋关节点云，为了获得更加光滑的表面，使用常用的数据降维技术PCA对初始髋关节点云进行重采样，通过找到数据中主要方向来降低数据维度，从而实现对点云的重采样，使用PCA能够在髋关节点云保持总体形状不变的同时去除点云中的噪声，且保持髋关节点云上均匀的密度分布，这为下游的重建与定量计算任务奠定了良好的基础，提高了后续髋关节点云处理的效率；
使用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投影到选取的特征向量上，重新计算髋关节点云的坐标来实现点云的降维，这样就得到了需要的已经降维的髋关节点云数据集/>
步骤3.2，无定向髋关节点云的迭代泊松重建：
由于髋关节各处的点法向量变化较大，难以准确计算，而传统的泊松重建方法又需要点的法向量作为输入，在髋关节建模部分选用的是基于无定向点的迭代泊松重建方法，这种方法不需要点的法向量作为输入，而是通过不断迭代来获得各点更加精确的法向量信息，当分配随机法线输入点，泊松重建生成一个表面，在每次迭代中，算法以直接从前一次迭代得到的曲面计算法线样本作为输入点，然后建立质量较好的新髋关节曲面；
步骤4，髋关节三维建模处理：
给定一组髋关节无方向点Pu＝{P1,P2,...,Pu}，其中u为无方向点的个数，最大八叉树深度H，收敛阈值δ和筛选的PSR权重γ作为输入，首先构造一个深度最大为H的八叉树，并使用八叉树节点作为样本集O，构造一个kd树以方便样本搜索，对每个样本oi∈O分配一个随机的法向量，之后在每次迭代中，将筛选出的PSR应用于具有当前法线的样本集O，得到一个指标函数x，应用面向八叉树的marching cube算法提取具有等值的等面F，其中，n＝|O|为样本数，对于等面的每一个三角形面，找到最接近它的k个样本，建立三角形面和样本之间的关系，在迭代计算中通过计算与oi相关联的面的法线的平均值来更新每个样本的法线/>最后计算当前迭代和之前迭代的样本之间的法线方差，当法线方差的平均值小于指定的阈值δ时，算法终止，通过迭代泊松重建方法能够得到更加精确的髋关节重建模型；
步骤5，提取股骨头部点云Pn，通过球面拟合算法，得到股骨颈中心Chead与半径Rhead：
根据股骨头在整个股骨的占比选取三个点确定一个平面，用平面截取属于股骨头上的点云数据集合Pn，用球面拟合算法对点云集合Pn进行拟合运算，假设球心为(x0,y0,z0),半径为R，(x,y,z)为点云集合Pn中的点：
建立一个球面方程：
(x-x0)2+(y-y0)2+(z-z0)2＝R2 ⑷，
建立误差函数：
误差函数E表示所有目标点到拟合球心的误差平方和，如令：
s＝x2+y2+z2-R2 ⑹，
则误差函数E又能够记作：
当误差函数E为最小值时，球心(x0,y0,z0)是E的最优估计值，E分别对参数x0,y0,z0,s取偏导数，如下式(8)：
通过矩阵进行求解，求得x0,y0,z0,s的值，确定出半径Rhead的值，并根据一组在球面上的点拟合出了股骨头球体，得到了股骨头中心Chead和股骨头半径Rhead；
步骤6，基于单层球面截取股骨颈部轮廓线方法提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心Cneck：
建立球心在股骨头中心Chead(x0,y0,z0)，半径为Rhead+ε的球面方程：(x-x0)2+(y-y0)2+(z-z0)2-(Rhead+ε)2＝0，其中ε是一个微小偏移量，ε∈[3mm,5mm]，将股骨点云数据代入方程，得到球面与股骨颈的相交轮廓线附近的点云数据。将得到的球面与股骨颈的相交轮廓线附近的点云数据拟合成三维空间中的椭圆，取椭圆的中心Cneck为密度中心，连接Chead和Cneck得到股骨颈轴线，记方向向量为n_neck；
步骤7，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线：
在通过高斯曲率寻找股骨头颈相交线时，首先寻找每个点的邻域点。固定距离球状邻域搜索首先设定距离阈值r，通过计算当前点与待定点之间的三维空间距离，将距离小于预设半径的所有待定点视为当前点的球状邻域点，设点p及其周围邻域内m个最近邻点为qi(i＝1,2,3,.,m)。设N为点p的法向量，Mi为点qi的法向量,则p点相对于qi的法曲率kn i估计如下式(9)：
式中:为向量N与向量pqi之间的夹角；β为向量N与向量Mi之间的夹角；|pqi|为p与qi之间的欧氏距离。
根据欧拉公式,法曲率与主曲率有以下关系，如下式(10):
其中：i＝1,2,...,m；k1和k2是主曲率；θi+θ为点p过qi的法截线的切线与主方向的夹角。
令最终求主曲率的问题可以转化成求解下列方程最小二乘：
求解上述最小二乘以获得主曲率k1与k2,而k1×k2即为高斯曲率，得到每个点的高斯曲率值后，通过限定值的大小来找出股骨头头颈相交线附近的点云数据，然后通过选出的点云数据拟合成一个空间椭圆，得到股骨头颈相交线；
步骤8，提取盆骨髋臼外缘线上的点云，拟合成平面S1，用S1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线：
在求髋臼外缘线时，首先求出股骨头和股骨颈的半径R1、R2和球心C1、C2以及盆骨的质心O，令股骨头上的每一点P与股骨头球心C1构成一条射线L1，分别计算盆骨上的点Q与L1的距离，若距离小于设定的阈值D，则将该点纳入盆骨球冠点云中，令盆骨球冠上的点与股骨颈中心C2的连线记作L2，盆骨质心与股骨颈中心的连线记作L3，计算L2与L3形成的角度，从中选取最大角度范围内的点构成盆骨外缘线，拟合求得盆骨外缘线所在的平面方程Z＝AX+BY+C，利用平面截取股骨头，该平面与股骨头相交的点即构成髋臼外缘线，至此，得到髋臼外缘线的点云，通过选出的点云数据拟合成一个空间椭圆，得到髋臼外缘线；
步骤9，计算α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标；
步骤9.1，计算α角：α角是股骨头球心Chead和股骨头颈相交线上t点钟方向点T1的连线与股骨头颈轴线n_neck的夹角，连接Chead、T1为向量v1，则α角可记为：θα＝arccos((v1·n_neck)/(|v1||n_neck|))；
步骤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；
步骤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的夹角；
步骤9.4，计算CE角：股骨头的中心为Chead，E点为髋臼最外侧处，经过Chead、E两点的直线与身体中线的平行线之间的夹角即为CE角，E为髋臼外缘线12点钟方向的点，股骨头中心Chead在Z轴上的向量为v6，连接E与Chead记为向量v7，则CE角可记为：θCE＝arccos((v6·v7)/(|v6||v7|))；
步骤9.5，输出α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个评价指标的具体数据。
本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明构思和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内，本发明要求保护范围由所附的权利要求书界定。
Claims (10)
1.一种基于CT图像的髋关节三维建模的自动定量计算方法，具体步骤如下：
步骤1，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作；
步骤2，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk，和股骨点云Pm；
步骤3，对髋关节的点云数据采用主成分分析方法重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型；
步骤4，髋关节三维建模处理；
步骤5，提取股骨头部点云Pn通过球面拟合算法，得到股骨颈中心Chead与半径Rhead；
步骤6，基于单层球面截取股骨颈部轮廓线方法提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心Cneck；
步骤7，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线；
步骤8，提取盆骨髋臼外缘线上的点云，拟合成平面S1，用S1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线；
步骤9，计算并输出α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标。
2.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤1中，对髋关节CT图像进行预处理，包括滤波、二值化、去枕床操作：在进行图像分割之前，需要先对原始CT图像进行去噪、二值化、去枕床的预处理操作，先使用中值滤波消除图像的部分噪声，再使用大津算法对灰度图像实现二值化，鉴于盆骨和股骨的点云数据在图像的上半部分而枕床的点云靠近下方，故根据点云的这一分布特征去除枕床部分的点云。
3.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤2中，对预处理后的髋关节CT图像基于轮廓线进行分割，得到盆骨点云Pk和股骨点云Pm：
采用基于边界的方法提取图像轮廓线，分为像素梯度计算、非极大值抑制、滞后阈值处理以及孤立弱边缘抑制，在进行图像轮廓线的提取过程中，使用Sobel算子计算像素梯度：Sobel算子为两个3×3的矩阵Sx、Sy，前者用于计算图像x方向像素梯度矩阵Gx，后者用于计算图像y方向像素梯度矩阵Gy，如下式(1)、(2)：
其中，I为CT灰度图像矩阵，按照下式(3)：
计算得到梯度强度矩阵Gxy，通过梯度强度矩阵Gxy，对当前像素梯度强度与沿正负梯度方向上的相邻像素的梯度强度进行比较，若其最大即为极值，则保留该像素为边缘点，若非最大，则对其进行抑制，不将其作为边缘点，之后定义一个高阈值和一个低阈值，梯度强度低于低阈值的像素点被抑制，不作为边缘点；高于高阈值的像素点被定义为强边缘，保留为边缘点；处于高低阈值之间的定义为弱边缘，留待进一步处理，在高阈值图像中把边缘链接成轮廓，当到达轮廓的端点时，本算法会在断点的8邻域点中寻找满足低阈值的点，再根据此点收集新的边缘，直到整个图像边缘闭合，至此，得到一个最为接近图像真实边缘的边缘图像，最终根据盆骨和股骨提取出的轮廓线长度、位置和形态差异，分割得到盆骨与股骨单独的点云Pk和Pm，其中，k为盆骨点云个数，m为股骨点云个数。
4.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤3中，对髋关节的点云数据采用PCA重采样后，通过基于无定向点的迭代泊松重建进行髋关节三维模型：
步骤3.1，基于PCA对髋关节点云进行降维，对于从CT图像中提取的轮廓线，我们通过间隔均匀采样获取初始髋关节点云，为了获得更加光滑的表面，使用常用的数据降维技术PCA对初始髋关节点云进行重采样，通过找到数据中主要方向来降低数据维度，从而实现对点云的重采样，使用PCA能够在髋关节点云保持总体形状不变的同时去除点云中的噪声，且保持髋关节点云上均匀的密度分布，这为下游的重建与定量计算任务奠定了良好的基础，提高了后续髋关节点云处理的效率；
使用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投影到选取的特征向量上，重新计算髋关节点云的坐标来实现点云的降维，获得需要的已经降维的髋关节点云数据集
步骤3.2：无定向髋关节点云的迭代泊松重建：
无定向髋关节点云的迭代泊松重建不需要点的法向量作为输入，而是通过不断迭代来获得各点更加精确的法向量信息，当分配随机法线输入点，泊松重建生成一个表面，在每次迭代中，算法以直接从前一次迭代得到的曲面计算法线样本作为输入点，然后建立质量较好的新髋关节曲面。
5.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤4中，髋关节三维建模处理步骤如下：
给定一组髋关节无方向点Pu＝{P1,P2,...,Pu}，其中u为无方向点的个数，最大八叉树深度H，收敛阈值δ和筛选的PSR权重γ作为输入，首先构造一个深度最大为H的八叉树，并使用八叉树节点作为样本集O，构造一个kd树以方便样本搜索，对每个样本oi∈O分配一个随机的法向量，之后在每次迭代中，将筛选出的PSR应用于具有当前法线的样本集O，得到一个指标函数x，应用面向八叉树的marching cube算法提取具有等值的等面F，其中，n＝|O|为样本数，对于等面的每一个三角形面，找到最接近它的k个样本，建立三角形面和样本之间的关系，在迭代计算中通过计算与oi相关联的面的法线的平均值来更新每个样本的法线/>最后计算当前迭代和之前迭代的样本之间的法线方差，当法线方差的平均值小于指定的阈值δ时，算法终止，通过迭代泊松重建方法能够得到更加精确的髋关节重建模型。
6.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤5中，提取股骨头部点云Pn，通过球面拟合算法得到股骨颈中心Chead与半径Rhead：
根据股骨头在整个股骨的占比选取三个点确定一个平面，用平面截取属于股骨头上的点云数据集合Pn，用球面拟合算法对点云集合Pn进行拟合运算，假设球心为(x0,y0,z0),半径为R，(x,y,z)为点云集合Pn中的点：
首先建立一个球面方程：
(x-x0)2+(y-y0)2+(z-z0)2＝R2 ⑷，
建立误差函数：
误差函数E表示所有目标点到拟合球心的误差平方和，如令：
s＝x2+y2+z2-R2 ⑹，
则函数E又能够记作：
当误差函数E为最小值时，球心(x0,y0,z0)是E的最优估计值，E分别对参数x0,y0,z0,s取偏导数，如下式(8)：
通过矩阵进行求解，从而求得x0,y0,z0,s的值，就能够确定出半径Rhead的值，并根据一组在球面上的点拟合出了股骨头球体，得到了股骨头中心Chead和股骨头半径Rhead。
7.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤6中，基于单层球面截取股骨颈部轮廓线方法，提取颈部轮廓线，取轮廓线的密度中心作为股骨颈中心Cneck：
建立球心在股骨头中心Chead(x0,y0,z0)，半径为Rhead+ε的球面方程：(x-x0)2+(y-y0)2+(z-z0)2-(Rhead+ε)2＝0，其中，ε是一个微小偏移量，ε∈[3mm,5mm]，将股骨点云数据代入球面方程，得到球面与股骨颈的相交轮廓线附近的点云数据，将得到的球面与股骨颈的相交轮廓线附近的点云数据拟合成三维空间中的椭圆，取椭圆的中心Cneck为密度中心，连接Chead和Cneck得到股骨颈轴线，记方向向量为n_neck。
8.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤7中，计算股骨点云的高斯曲率，获得股骨头与股骨颈的相交处的点云，将点云投影到由点云本身拟合的平面后基于最小二乘法拟合成椭圆得到股骨头颈相交线：
在通过高斯曲率寻找股骨头颈相交线时，首先寻找每个点的邻域点，固定距离球状邻域搜索首先设定距离阈值r，通过计算当前点与待定点之间的三维空间距离，将距离小于预设半径的所有待定点视为当前点的球状邻域点，设点p及其周围邻域内m个最近邻点为qi(i＝1,2,3,.,m)，设N为点p的法向量，Mi为点qi的法向量,则p点相对于qi的法曲率估计如下式(9)：
式中:为向量N与向量pqi之间的夹角；β为向量N与向量Mi之间的夹角；|pqi|为p与qi之间的欧氏距离；
根据欧拉公式，法曲率与主曲率有以下关系，如下式(10):
其中：i＝1,2,...,m；k1和k2是主曲率；θi+θ为点p过qi的法截线的切线与主方向的夹角，
令最终求主曲率的问题转化成求解下列方程最小二乘，如下式(11)：
求解最小二乘以获得主曲率k1与k2，而k1×k2即为高斯曲率，得到每个点的高斯曲率值后，通过限定值的大小来找出股骨头头颈相交线附近的点云数据，再通过选出的点云数据拟合成一个空间椭圆，得到股骨头颈相交线。
9.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤8中，提取盆骨髋臼外缘线上的点云，拟合成平面S1，用S1截取股骨头上的点云数据，拟合成椭圆得到髋臼外缘线：
在求髋臼外缘线时，首先求出股骨头和股骨颈的半径R1、R2和球心C1、C2以及盆骨的质心O，令股骨头上的每一点P与股骨头球心C1构成一条射线L1，分别计算盆骨上的点Q与L1的距离，若距离小于设定的阈值D，则将该点纳入盆骨球冠点云中，令盆骨球冠上的点与股骨颈中心C2的连线记作L2，盆骨质心与股骨颈中心的连线记作L3，计算L2与L3形成的角度，从中选取最大角度范围内的点构成盆骨外缘线，拟合求得盆骨外缘线所在的平面方程Z＝AX+BY+C，利用平面截取股骨头，该平面与股骨头相交的点即构成髋臼外缘线，至此，得到髋臼外缘线的点云，通过选出的点云数据拟合成一个空间椭圆，得到髋臼外缘线。
10.根据权利要求1所述基于CT图像的髋关节三维建模的自动定量计算方法，其特征在于，步骤9中，计算α角、CE角、股骨头突出指数e及股骨头颈偏心距这四个股骨髋臼撞击综合征评价指标：
步骤9.1，计算α角：α角是股骨头球心Chead和股骨头颈相交线上t点钟方向点T1的连线与股骨头颈轴线n_neck的夹角，连接Chead、T1为向量v1，则α角可记为：θα＝arccos((v1·n_neck)/(|v1||n_neck|))；
步骤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；
步骤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的夹角；
步骤9.4计算CE角：股骨头的中心为Chead，E点为髋臼最外侧处，经过Chead、E两点的直线与身体中线的平行线之间的夹角即为CE角，E为髋臼外缘线12点钟方向的点，股骨头中心Chea在Z轴上的向量为v6，连接E与Chead记为向量v7，则CE角可记为：θCE＝arccos((v6·v7)/(|v6||v7|))。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310832350.3A CN117078840A (zh) | 2023-07-07 | 2023-07-07 | 基于ct图像的髋关节三维建模的自动定量计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310832350.3A CN117078840A (zh) | 2023-07-07 | 2023-07-07 | 基于ct图像的髋关节三维建模的自动定量计算方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117078840A true CN117078840A (zh) | 2023-11-17 |
Family
ID=88718250
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310832350.3A Pending CN117078840A (zh) | 2023-07-07 | 2023-07-07 | 基于ct图像的髋关节三维建模的自动定量计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117078840A (zh) |
-
2023
- 2023-07-07 CN CN202310832350.3A patent/CN117078840A/zh active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Douglas | Image processing for craniofacial landmark identification and measurement: a review of photogrammetry and cephalometry | |
| CN107481228B (zh) | 基于计算机视觉的人体背部脊柱侧弯角度测量方法 | |
| Kauffmann et al. | Computer-aided method for quantification of cartilage thickness and volume changes using MRI: validation study using a synthetic model | |
| US8965108B2 (en) | Method and system of automatic determination of geometric elements from a 3D medical image of a bone | |
| US8923584B2 (en) | Method and system of automatic determination of geometric elements characterizing a bone deformation from 3D image | |
| CN108618749B (zh) | 基于便携式数字化眼底照相机的视网膜血管三维重建方法 | |
| CN106997594B (zh) | 一种眼部组织的定位方法及装置 | |
| JP2010279438A (ja) | 画像処理装置、画像処理方法及びコンピュータプログラム | |
| CN111563901B (zh) | 基于磁共振的髋关节图像处理方法及系统、存储介质、设备 | |
| Schumann et al. | An integrated system for 3D hip joint reconstruction from 2D X-rays: a preliminary validation study | |
| US11551371B2 (en) | Analyzing symmetry in image data | |
| Trinh et al. | Accurate measurement of cartilage morphology using a 3D laser scanner | |
| Linney et al. | Three-dimensional visualization of data on human anatomy: diagnosis and surgical planning | |
| Liu et al. | Improve accuracy for automatic acetabulum segmentation in CT images | |
| CN117078840A (zh) | 基于ct图像的髋关节三维建模的自动定量计算方法 | |
| Fischer et al. | Automated morphometric analysis of the hip joint on MRI from the German National Cohort Study | |
| Bhatia et al. | Practical surface patch registration technique | |
| Shui et al. | A Landmark-free Approach for Surface Asymmetry Detection and Profile Drawings from Bilaterally Symmetrical Geometry | |
| Bhatia et al. | A Practical Surface Patch Registration | |
| Nakagawa et al. | Comparison of the depth of an optic nerve head obtained using stereo retinal images and HRT | |
| CN115170452A (zh) | 一种基于空间中心轴的髁突cbct影像数据增强方法 | |
| Liu et al. | Model-based 3D segmentation of the bones of joints in medical images | |
| CN112330723A (zh) | 一种物理到图像/图像到物理的自动配准方法 | |
| EP3948776A1 (en) | Closed surface fitting for segmentation of orthopedic medical image data | |
| CN117958970A (zh) | 一种基于ct与激光口扫的口腔手术实时导航方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |