基于股骨三维多线外形计算颈干角的方法

何姗姗 刘 浩 唐月红(南京航空航天大学理学院，南京 106)

2(南京航空航天大学机电学院，南京 211106)

引言

近年来，随着我国人口老龄化逐渐达到高峰，股骨颈骨折成为一种常见的骨科疾病，而股骨特征提取则成为当今医学研究的主题之一。由于股骨外形结构复杂多样，所以几何建模成为面向医疗应用的CAD/CAE/CAM 技术的关键。三维几何模型可用于临床应用和病理研究，不仅在三维观察和三维测量中起到重要作用，而且也用于人工假体的数字化设计、智能化虚拟手术、计算机辅助手术导航和三维有限元分析［1］。

股骨近端形态和大小在不同个体间存在明显差异，不同年龄、性别、种族的形态有所不同。前些年，我国股骨假体大多从外国进口，但由于差异性的存在，在外形匹配上存在缺陷。近年来，人工假体设计成为我国医学研究的热点，其中不同个体设计符合个体的假体成为研究的重点，而如何描述个体的股骨近端形态是假体设计的必要条件。

股骨头中心、股骨头直径、颈干角是股骨近端的主要特征，它们为定制合理的股骨假体及医学手术提供条件［2-3］。在医学和形态学上，关于股骨特征的提取都有研究。目前，关于股骨近端特征的提取大多是通过基于X 射线投影的人工测量方法［4-6］，进而获得股骨头中心、颈轴线、干轴线和颈干角等特征数据。但是，对于股骨复杂的三维结构，二维测量具有较大的误差性。文献［7 -8］利用CT 影像和逆向工程，分析股骨近端的三维几何特征，具有很强的突破性和精确性。医学研究表明，股骨颈干角的大小与股骨骨折亦有关联［9］。笔者试图利用数学方法，从股骨几何形态学的角度分析股骨形态，根据不同位置的几何特征来描述股骨，从而寻找计算股骨颈干角的方法。

本研究基于CT 获得髋关节的断层数据，导入医学处理软件进行图像处理和三维建模，获得了股骨的轮廓线并存成IGES 格式，通过Matlab 读入IGES 生成点云数据，然后通过三维最小二乘球拟合来得到股骨头中心和股骨头直径，同时通过二维圆拟合和空间直线拟合来得到干轴线，最后根据已知颈干角的测量数据，通过优化迭代来拟合最佳颈干角，找到颈轴线的位置。

1 股骨外形分析

股骨是人体中最大的长管状骨，可分为一体两端(见图1)。上端朝向内上方，其末端膨大呈球形，叫“股骨头”，股骨头的外下方较细的部分称“股骨颈”。颈与体的夹角称“颈干角”，男性平均132°，女性平均127°。颈体交界处的外侧有一向上的隆起，叫做大转子，其内下方较小的隆起叫做小转子。股骨颈的下方称为股骨体。

图1 股骨上端特征(1—股骨头;2—股骨颈;3—大转子;4—小转子;5—股骨体)Fig． 1 The upper characteristics of femur (1 －femoral head; 2 － femoral neck; 3 － greater trochanter;4 －lesser trochanter;5 － femoral shaft)

髋关节由股骨和髋骨构成，而股骨头正是连接两处2/3 的近似球体，其中最主要的参数是球心坐标和球体半径。在假体设计中，可以利用球拟合来拟合股骨头，补充股骨头的损坏缺陷，从而实现假体与髋臼的良好匹配。

股骨体小转子下端可以近似地看做一个圆柱，其主要参数是中心轴，即为股骨的干轴线。把干轴逐层圆面拟合，再把得到的圆心直线拟合，从而得到干轴线。

对于股骨颈，其主要参数是股骨的颈轴线。根据已得到的股骨干轴和股骨头中心坐标设计算法，推导出颈轴线方程。

2 股骨颈干角计算方法

2.1 股骨的多线外形提取

本研究利用CT 原始的DICOM 文件，导入医学处理软件DataSet 进行基于灰度的图形图像处理，利用阈值分割和区域增长来逐层提取股骨边界轮廓，这个过程如图2 所示。

图2 提取股骨切片边界轮廓Fig．2 Femoral section boundary contour extraction

将逐层提取出的股骨轮廓的多线外形存成IGES 数据，通过Matlab 读入IGES 数据，这些数据由每层 z 坐标相同的若干数据点组成。通过Matlab，可实现对股骨多线外形的光顺、平移、旋转、拟合等分析处理，其中股骨的多线外形数据如图3所示。

图3 股骨多线外形Fig．3 The multiline shape of femur

2.2 股骨颈干轴拟合

2.2.1 股骨头中心和直径的确定文献［8，10 -12］都把股骨头看成一个球体来拟合，笔者也采用该方法。不同的是，上述文献数据为点云数据，在描述股骨外形时数据点数巨大。例如，文献［10］中股骨点云有559 816 个数据点，庞大的数据点将不利于外形编辑，因此分别提出了不同的方法来精简点云数据和剔除误差点［10-11］。而本研究的数据是多线外形数据，是一种特殊的点云，介于有序与无序之间。例如，图3 中的股骨点云数据约有18 684 个数据点，借用多线外形这一特性，可以直接来拟合股骨头。在图3 所示的股骨截面线数据上，观察从股骨体到股骨头的状况，在同一z 坐标截面线从1 条变为2 条再变为1 条。由于原始数据仅由分别位于同一z 坐标的点集组成，找出全部位于股骨头上的数据尚有难度。为了筛选出拟合股骨头的有效数据，选截面线从2 条变为1条开始以上部分的股骨头数据。拟合股骨头的具体方法如下:

设球心坐标为(Ox，Oy，Oz)，半径为R，则球面方程为

令a = - 2Ox，b = - 2Oy，c = - 2Oz，d = Ox2 +Oy2 + Oz2 - R2，则对于数据集(xi，yi，zi)(i = 1，2，…，n)，任一点(xi，yi，zi)到球心的距离的平方为

δ = xi

2 + yi

2 + zi2 + axi+ byi+ czi+ d + R2

求解式(1)，得出a，b，c，d。经过换算，即求出股骨头中心坐标(Ox，Oy，Oz)和股骨头半径R。

同时，对每个采样点(xi，yi，zi)(i = 1，2，…，n)计算拟合残差，有进而计算拟合误差，有

利用式(2)～式(3)计算误差，通过20 个算例计算股骨头的数据点数平均为2 209个，且拟合平均误差为0.434 0 mm，如表1 所示。

可以看出，这种数据的采用方法和拟合方法对整个球体拟合造成的误差相对较少，且具有一定的稳定性，同时在性能上可以通过股骨头直径与正常值的对比来评价，将在文中第3 节具体分析。因此这种方法可以有效地的拟合股骨头。

表1 股骨头拟合实验结果Tab．1 Experimental results of fit femoral head

2.2.2 股骨干轴线的确定

如图1 所示，股骨体小转子下端形似一个圆柱，因此对于股骨干轴的拟合，可把截面线的每一层看做二维圆拟合，求出一组圆心坐标。由于这些圆心不在一条直线上，所以再通过空间直线拟合［13］的方法进行拟合，得到直线方程(即为干轴方程)L:x =az +b，y =cz +d(其中z 为变量)。同时，利用2.2.1中类似的方法，可以计算出干轴直线拟合时的平均误差。

2.2.3 优化求解最佳颈干角

根据现有的关于颈干角测量的研究数据［2，5-6］，颈干角正常值在110° ～140°之间，男性平均为132°，女性平均为127°。颈干角随年龄的增大而减小，儿童的颈干角大于成年人的颈干角。

从股骨形态学上分析，股骨颈类圆柱状，中间细两头粗，因此与颈轴线垂直的平面截颈轴得到类圆形的截面，且截面积最小，而在同处垂直于不精确颈轴的平面截得的横截面面积必定大于真实颈轴截得的面积。从这个理论分析出发，寻找精确颈轴。把颈干角α 在110° ～140°之间迭代搜索最佳颈干角来找寻精确颈轴。由2.2.1 和2.2.2 求出的球心坐标和干轴线确定一平面，在该平面上过股骨头中心O 点做与干轴L 成α 夹角的直线与L 交于B点，如图4 所示。

设B 点坐标为(at +b，ct +d，t)，可知

当α = 110°时，颈干轴交点为B1，其坐标为(at1+ b，ct1+ d，t1);当α =140°时，颈干轴交点为B2，其点坐标为(at2+ b，ct2+ d，t2)。

当α ∈(110°，140°)时，颈干轴交点B 的z 坐标t∈(t1，t2)。OB 为假定颈轴。过线段OB 中点C 做垂直于OB 的平面，截面线数据上分别在每层逐点连线，与平面相交点集D，如图5 所示。

由图5 看出，截面点位于空间的一个平面中，既可以采用空间曲线拟合截面点来计算截面面积，亦可以通过旋转空间平面到与坐标平面平行来计算截面面积。由于曲线拟合和旋转的复杂性以及可能带来的误差，所以可以采用投影的方法来计算截面面积。通过投影面积反求截面面积，使面积的计算变得简单，从而减少了拟合和旋转的复杂性与误差性，具体方法如下:

已知截面数据点存为点集D，原则上D 为椭圆形状，为了求D 所围成图形面积，把D 投影到xoy平面上，得到数据点D1，再利用定积分求解平面点集D1的面积S1，从而反求点集D 所围成的面积S。这样，就能够避免直接拟合椭圆求面积而造成的误差。

图4 颈轴线拟合Fig．4 The fitting of neck axis

图5 截面点数据Fig．5 Data of the cross-section

已知数据点D 所在平面为:

设a = Ox- Bx，b = Oy- By，c = Oz- Bz，则式(4)表示的平面法向量为(a，b，c)。该法向量与z轴夹角的余弦为cosθ，则原数据集D 点集所表示图形的面积为S = S1/cosθ。

这样，就得到截面面积S 是变量t 的函数。当t∈(t1，t2)时，根据横截面面积最小来建立优化模型，即

求解单变量约束最优化式(5)，得到最优解t 。

2.2.4 颈轴线和颈干角的确定:

根据上节求得的最优解，即找到最佳颈干轴交点B 的坐标，这时OB 所在直线就是所求的颈轴线。同时，OB 与干轴线L 的夹角即为该股骨的颈干角。拟合出颈干轴直线，如图6 所示。

图6 拟合出的颈干轴。(a)例1;(b)例2Fig．6 The fitting of neck axis and shaft axis． (a)Example 1;(b)Example 2

3 结果

基于CT 获得髋关节的断层数据，通过导入医学处理软件来进行图像处理和三维建模，获得了股骨的轮廓线并存成IGES 格式，通过Matlab 读入IGES 生成点云数据，再通过优化步骤来得到颈干轴线的位置，同时确定股骨颈干角的大小。该研究主要基于数学方法，通过分析股骨外形结构的角度来计算股骨颈干角。用20 例CT 数据来验证该算法，得到的主要股骨特征数据，如表1 ～表2所示。

表2 颈干角实验结果Tab．2 Experimental results of neck-shaft angle

由于采用的软件在前期截面线的提取尚未完成自动化，所以仍需交互式提取。因提取的复杂性，笔者仅通过20 例CT 数据来进行测量。表1 的数据表明，该方法在拟合股骨头时的平均误差为0.434 0 mm，误差较小，具有很好的拟合效果。在通过圆拟合和直线拟合干轴线时，拟合的平均误差为0.384 1 mm。通过表2 的方差分析可以看出，股骨直径的测量值稳定性较强，同时股骨颈干角测量的平均值是126.325 9°。方差分析表明，本算法具有一定的稳定性，但仍有待减小方差值，进一步的验证工作将在截面线自动提取实现后进行。

4 讨论和结论

利用数学方法，根据股骨外形的几何特征来求解颈干角。具体实验表明，这种方法能快速有效地计算出股骨颈干角，同时可以表示出颈干轴方程。与现有方法相比，有以下优点:第一，为求解颈干角提供了一个新的思路;第二，与二维测量相比，三维测量具有明显的优点;第三，与现有的三维测量方法相比，本方法比文献［8］中的方法更具简便性。文献［8］在求解最小截面时，通过确定一条假定颈轴做一系列截面来确定最小截面，其中并没有说明假定颈轴的具体取法，并且直接对截面进行椭圆拟合会有一定的误差;本研究通过干轴方程和夹角范围确定假定颈轴，从而确定最小截面，相对简单，并且在对截面积的计算时采用投影及定积分求解，具有更多的简便性，可避免不必要的误差性。但是，本方法还存在一定的缺陷，例如，在注重数学方法探讨的同时，需加强实例验证。同时，对误差精度的估计也需进一步探讨。

总之，本研究提出的求解颈干角的方法对于人体假肢的定制提供了一定的参考，并对股骨的变形研究工作提供了前提条件。

［1］ Sun W，Starly B，Nam J，et al． Bio-CAD modeling and its applications in computer-aided tissue engineering［J］. Computer Aided Design，2005，37 (11):1097 -1114.

［2］ 秦岭.基于医学透视图像的股骨上段特征提取［D］. 长沙:中南大学，2010.

［3］ 胡荣慧，巫北海，张绍祥，等. 髋关节CT 三维重建与可视化研究［J］.第三军医大学学报，2007(14):1371 -1373.

［4］ 张怀瑫，郑靖中，杨玉田.国人股骨颈干角及扭转角的测量统计［J］.解剖学报，1982，3:262 -230.

［5］ 谭进.国人股骨的形态学特征及其临床意义［J］. 四川解剖学杂志，2002，10(3):163 -167.

［6］ 郑超，马信龙，马剑雄，等. 大鼠股骨颈干角及前倾角的测量［J］.中国中西医结合外科杂志，2008，14(5):460 -461.

［7］ 李凤春，石福明，周红平，等. 股骨颈干角计算方法的探讨［J］.中国中医骨伤科杂志，1998，6(5):24 -25.

［8］ Mahaisavariya B，Sitthiseripratip K， Tongdee T， et al.Morphological study of the proximal femur:a new method of geometrical assessment using 3 - dimensional reverse engineering［J］.Medical Engineering ＆ Physics，2002，24(9):617 -622.

［9］ 雷田，王岐本，邝满元，等. 颈干角的大小与股骨颈骨折的关系［J］.解剖与临床，2011，16(2):166 -167.

［10］ 宋卫卫，苏铁明，欧宗瑛，等.基于CT 图像反求技术的人体股骨头修复建模［J］. 大连理工大学学报，2009(2):222 -227.

［11］ 陈树广. 基于病例CT 图像的健康股骨头建模研究［J］. 计算机仿真，2009，26(9):172 -175.

［12］ Cibulka MT. Determination and significance of femoral neck anteversion［J］.Physical Therapy，2004，84(6):550 -558.

［13］ 袭杨. 空间直线拟合的一种方法［J］. 齐齐哈尔大学学报，2009，25(2):64 -68.