CN103198521A - 一种面向个性化设计的股骨头表面三维模型重构方法 - Google Patents

Abstract

本发明属于医学图像处理领域，特别涉及一种面向个性化设计的股骨头表面三维模型重构方法。该方法包括下述步骤：对输入的髋关节CT序列图像中每一幅二维CT图像内的股骨头和髋臼进行轮廓提取与内部填充；利用像素间的连通性实现每一幅二维CT图像内髋臼的分割操作；利用已分割完毕的髋臼序列图像，进行髋臼三维模型重构；利用空间辐射线法，进行髋臼凹面数据点采样；采用基于一般二次曲面约束的椭球面拟合方法，实现个性化股骨头表面重构。该方法能够为股骨头假体的个性化设计提供精确的数据模型。

Description

一种面向个性化设计的股骨头表面三维模型重构方法

技术领域

本发明属于医学图像处理领域，特别涉及一种面向个性化设计的股骨头表面三维模型重构方法。 

背景技术

股骨头缺血性坏死是一种世界性疾病，应用人工股骨头置换手术以挽救髋关节功能的人数日趋增多，但由于假体松动常常使人工股骨头失效，其主要原因在于假体与髋臼之间的不稳定匹配。有效提高股骨头假体与髋臼之间的长期稳定性是提高人工股骨头置换手术远期效果的关键。 

股骨头假体与髋臼匹配效果不佳，是由人体器官的个体差异性所造成的。不同患者的髋关节形状各不相同，而目前国内外生产的人工股骨头假体的规格、种类有限，因此为患者植入的人工股骨头假体大都不能与髋臼形成解剖匹配，接触面积小，假体与髋臼之间形成的是几个点或者局部区域的接触，而不是大范围的面接触，导致假体与髋臼接触面的应力分布极不合理，从而显著影响假体的稳定性，同时使患者术后股骨活动不自如。单个股骨头假体的使用寿命仅有3～15年，给股骨头坏死患者带来巨大生活影响和经济压力。因此，根据患者自身的个体器官参数，将坏死的股骨头表面重构为与髋臼良好匹配状态下的三维模型，为“个体化”、“定做式”的股骨头假体生产提供几何数据参数，已成为当前股骨头坏死的康复治疗中亟待解决的问题。 

针对股骨头表面模型的个性化三维重构，国内学者提出了一些解决方法，发表的文献主要包括：《西安交通大学学报》的“基于反求工程的个体匹配化骨骼制造方法的研究”；《大连理工大学学报》的“基于CT图像反求技术的人体股骨头修复建模”。此类技术主要利用CT（或MRI）体数据图像对坏死的股骨头 进行三维重构，然后根据生物力学分析结果对所得到的模型进行曲面离散化二次重构。此类方法的精度相比传统假体设计方法有所提高，但其设计周期较长，重构过程较复杂。 

国外学者针对股骨头表面个性化三维重构所提出的方法包括：在《International Journal of Computer Assisted Radiology and Surgery》的“Articular surface remodeling of the hip after periacetabular osteotomy”和“An integrated approach for reconstructing a surface model of the proximal femur from sparse input data and a multi-resolution point distribution model:an in vitro study”。此类技术主要利用校准后X光片中的稀疏点数据，采用多节分辨率的点分布模型（MR-PDM）来重建个性化患者的股骨近端表面模型。此类方法的表面重构精度进一步提高，但其重构结果的稳定性较低，对重构过程的数据环境提出了苛刻的要求。 

发明内容

本发明的目的正是为了克服现有技术的不足，解决如何快速利用患者自身生物几何参数，进行股骨头表面个性化稳定重构的问题，而提出一种面向个性化设计的股骨头表面三维模型重构方法。该方法所得到的股骨头表面模型能够与髋臼实现良好的匹配状态，具有鲁棒性较强、计算过程方便等优点，能够为股骨头假体的个性化设计提供精确的数据模型。 

本发明的目的是通过下述技术方案实现的。 

一种面向个性化设计的股骨头表面三维模型重构方法，所述方法包括以下步骤： 

①输入由常规医用CT机所采集的髋关节CT序列图像； 

②对已输入的序列图像中每一幅二维CT图像内的股骨头和髋臼进行轮廓 提取与内部填充；先利用Sobel算子提取二维CT图像内股骨头和髋臼的轮廓，再利用扫描线法实现股骨头和髋臼轮廓的内部填充； 

③在已填充的股骨头和髋臼中，利用像素间的连通性实现每一幅二维CT图像内髋臼的分割操作，具体步骤如下： 

1）在股骨头和髋臼封闭区域内任取两个像素，如果这两个像素有连通的路径，则将这两个像素归为同一类，如果这两个像素没有连通的路径，则将这两个像素归为两类； 

2）遍历股骨头和髋臼封闭区域内的所有像素，将所有像素归为两类； 

3）舍弃属于股骨头封闭区域内的像素集合，保留属于髋臼封闭区域内的像素集合，实现髋臼的分割操作； 

④利用已分割完毕的髋臼序列图像，进行髋臼三维模型重构； 

⑤在已得到的髋臼三维模型上，利用空间辐射线法，进行髋臼凹面数据点采样，包括如下步骤： 

1）在髋臼三维模型的凹面内选择一个空间坐标点，以此坐标点为中心，向全周空间发射辐射直线； 

2）每条辐射直线与髋臼三维模型相交时即获得一个交点，所有交点的集合即为髋臼凹面数据采样点； 

⑥利用已得到的髋臼凹面数据点，采用基于一般二次曲面约束的椭球面拟合方法，实现个性化股骨头表面重构，具体步骤如下： 

1）将一般二次曲面约束为椭球曲面，设一般二次曲面的方程为： 

ax²+by²+cz²+2fyx+2gxz+2hxy+2px+2qy+2rz+d=0     （1） 

设： 

I=a+b+c     （2） 

J=ab+bc+ac-f²-g²-h     （3） 

$K=\left|\begin{array}{ccc}a& h& g\\ h& b& f\\ g& f& c\end{array}\right|---\left(4\right)$

对于任意一个椭球面都有α≥4，使得αJ-I²>0，选用α=4来将二次曲面约束为椭球曲面； 

2）利用髋臼凹面数据采样点，采用基于代数距离的最小二乘法进行椭球曲面拟合，获得公式（1）中的所有未知系数； 

3）设X=(x,y,z)^T，将椭球面一般二次曲面方程化归为矩阵形式： 

d+(2p,2q,2r)X+X^TKX=0     （5） 

对于矩阵K，使用Jacobi正交变换法求其特征根(λ₁,λ₂,λ₃)，及3×3的特征向量R，令 

(a₁,a₂,a₃)=(2p,2q,2r)R     （6） 

有 

$a_0=d-\frac{1}{4}(\frac{{a_1}^2}{{\mathrm{\λ}}_1}+\frac{{a_2}^2}{{\mathrm{\λ}}_2}+\frac{{a_3}^2}{{\mathrm{\λ}}_3})---\left(7\right)$

得到椭球面的三个主轴的长度 

$a=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_1}}$ $b=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_2}}$ $c=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_3}}---\left(8\right)$

设椭球面三个主轴的旋转角为(α,β,γ)，得到椭球面中心坐标系与主坐标系的变换关系 

X′=X₀+R₁(α)R₁(β)R₁(γ)X     （9） 

$X_0={(\frac{a_1}{{2\mathrm{\λ}}_1},\frac{a_2}{{2\mathrm{\λ}}_2},\frac{a_3}{{2\mathrm{\λ}}_3})}^T---\left(10\right)$

R₁(α)R₁(β)R₁(γ)=R^T     （11） 

附图说明

图1为本发明股骨头表面三维模型重构方法流程图。 

图2为本发明实施例中所涉及的髋关节CT图像。 

图3为本发明实施例中利用Sobel算子提取二维CT图像内股骨头和髋臼轮廓的示意图。 

图4为本发明实施例中利用扫描线法进行股骨头和髋臼轮廓的内部填充的示意图。 

图5为本发明实施例中实现股骨头和髋臼轮廓的内部填充后的示意图。 

图6为本发明实施例中在股骨头和髋臼区域内利用像素间的连通性进行像素分类过程的示意图。 

图7为本发明实施例中遍历股骨头和髋臼封闭区域内的所有像素，将所有像素归为两类后的示意图。 

图8为本发明实施例中舍弃属于股骨头封闭区域内的像素集合，保留属于髋臼封闭区域内的像素集合，实现髋臼的分割操作后的示意图。 

图9为本发明实施例中利用已分割完毕的髋臼序列图像，进行髋臼三维模型重构的示意图。 

图10为本发明实施例中在已得到的髋臼三维模型上，利用空间辐射线法，进行髋臼凹面数据点采样的示意图。 

图11为本发明实施例中利用已得到的髋臼凹面数据点，采用基于一般二次曲面约束的椭球面拟合方法，实现个性化股骨头表面重构后的示意图。 

图12为本发明实施例中个性化股骨头与髋臼良好匹配的示意图。 

具体实施方式

下面结合附图对本发明的较佳实施例进行详细阐述，以使本方明的优点和特征能够更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定。 

请参考图1，本发明提供一种面向个性化设计的股骨头表面三维模型重构方法，该方法包括下述步骤：输入髋关节CT序列图像；股骨头和髋臼边缘检测与内部填充；基于连通性的髋臼分割；髋臼三维模型重构；髋臼凹面数据点采样；个性化股骨头椭球面拟合。下面将详细介绍该一种面向个性化设计的股骨头表面三维模型重构方法。 

请参考图2～12，下面将具体介绍该一种面向个性化设计的股骨头表面三维模型重构方法。 

输入由常规医用多层螺旋CT机所采集到的病态髋关节序列CT图像。 

首先，对已输入的序列图像中每一幅二维CT图像内的股骨头和髋臼进行轮廓提取与内部填充，具体步骤为： 

利用Sobel算子提取二维CT图像内股骨头和髋臼的轮廓。定义Sobel算子为： 

s(i,j)=|Δ_xf|+|Δ_yf|=|(f(i-1,j-1)+2f(i-1,j)+f(i-1,j+1))-(f(i+1,j-1)+2f(i+1,j)+f(i+1,j+1))|+|(f(i-1,j-1)+2f(i,j-1)+f(i+1,j-1))-(f(i-1,j+1)+2f(i,j+1)+f(i+1,j+1))| 

其卷积算子为： 

${\mathrm{\Δ}}_{x}f=\left[\begin{array}{ccc}-1& 0& 1\\ -2& 0& 2\\ -1& 0& 1\end{array}\right],$ ${\mathrm{\Δ}}_{y}f=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]$

适当选取阈值T（本例中T=160），作如下判断：s(i,j)>T时，(i,j)为阶跃状边缘点，{s(i,j)}为边缘图像。 

在股骨头和髋臼轮廓线所在的图像中，从上至下，逐行生成横向扫描线，每条扫面线与股骨头和髋臼轮廓相交时会产生若干个交点，每条扫描线的第N个交点与第N+1（N=1,3,5,7,9,……）个交点之间可判定为股骨头和髋臼轮廓内部，其余为外部。据此，利用扫描线法实现了股骨头和髋臼轮廓的内部填充。 

然后，在已填充的股骨头和髋臼封闭区域内任取两个像素，以其中一个像素作为种子点，向8邻域方向进行迭代式扩展生长，若最终生长后的区域内包含另一像素，则判定这两个像素是连通的，即这两个像素有连通的路径，则将这两个像素归为同一类，如果这两个像素没有连通的路径，则将这两个像素归为两类。 

遍历股骨头和髋臼封闭区域内的所有像素，最终可将所有像素归为两类。 

舍弃属于股骨头封闭区域内的像素集合，保留属于髋臼封闭区域内的像素集合，即可实现髋臼的分割操作。据此，利用像素间的连通性实现了每一幅二维CT图像内髋臼组织的分割操作。 

利用已分割完毕的髋臼序列图像，采用Marching Cubes算法（设阈值为60），进行髋臼三维模型重构，得到髋臼组织的三维模型。 

在已得到的髋臼三维模型上，在髋臼窝凹面内选取一个空间坐标点（初始的采样点中心点），以该中心点为发射中心，向全周空间范围发射辐射直线（两条辐射线间夹角为20°），每条辐射直线与髋臼三维模型相交时即可获得一个交点，所有交点的集合即为髋臼凹面数据采样点集。 

最后，利用已得到的髋臼凹面数据点，采用基于一般二次曲面约束的椭球面拟合方法，实现个性化股骨头表面重构，具体步骤为： 

在直角三维坐标系X-Y-Z下，一般二次曲面方程 

ax²+by²+cz²+2fyx+2gxz+2hxy+2px+2qy+2rz+d=0 

设： 

I=a+b+c 

J=ab+bc+ac-f²-g²-h 

$K=\left|\begin{array}{ccc}a& h& g\\ h& b& f\\ g& f& c\end{array}\right|$

由二次曲面理论可知，I,J,K为二次曲面在坐标变换下的不变量，当J>0,I×K>0时，该二次曲面一定为椭球。当短轴至少为长轴的一半时，必有4J-I²>0，因此4J-I²>0是二次曲面为椭球的必要条件。对于股骨头曲面来说，其大致形状是与球体类似的，也就是说一般来说股骨头椭球面的短轴至少为长轴的一般，因此选用α=4，即4J-I²>0来将二次曲面约束为椭球曲面对于本例来说是可行的。 

设待进行拟合的髋臼凹面数据采样点集为 

$\{p_i(x_i,y_i,z_i){\}}_{i=1}^n$

根据点集中每个点的坐标，构建矩阵 

$x_i={(x_i^2,y_i^2,x_i^2,{2y}_i{z}_i,{2x}_i{y}_i,{2x}_i,{2y}_i,{2z}_i,1)}^T$

D=(X₁,X₂,X₃,...,X_n) 

由前述一般二次曲面方程，设 

v=(a,b,c,f,g,h,p,q,r,d)^T

在椭球条件约束下的基于代数距离最小二乘法可以用公式描述为 

minDv²subject tokJ-I²=1 

其中kJ-I²=1可以用矩阵等价表示为 

v^TCv=1 

$C=\left(\begin{array}{cc}{C}_1& 0_{6\×4}\\ 0_{4\×6}& 0_{4\×4}\end{array}\right)$

$C_1=\left(\begin{array}{cccccc}-1& \frac{k}{2}-1& \frac{k}{2}-1& 0& 0& 0\\ \frac{k}{2}-1& -1& \frac{k}{2}-1& 0& 0& 0\\ \frac{k}{2}-1& \frac{k}{2}-1& -1& 0& 0& 0\\ 0& 0& 0& -k& 0& 0\\ 0& 0& 0& 0& -k& 0\\ 0& 0& 0& 0& 0& -k\end{array}\right)$

使用拉格朗日乘子法，解方程组 

$\left\{\begin{array}{c}{\mathrm{DD}}^{T}v=\mathrm{\λCv}\\ v^T\mathrm{Cv}=1\end{array}\right.$

在k>3时，该方程组中的第一个方程只有一个解，即为C^-1DD^T中唯一一个正项特征值的特征向量。同时注意到C是一个稀疏矩阵，在行列号大于6时矩阵中的值都为0，因此将DD^T分解为 

${\mathrm{DD}}^T=\left(\begin{array}{cc}{S}_{11}& S_{12}\\ {S_{12}}^T& S_{22}\end{array}\right)$

$v=\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)$

S₁₁，S₁₂，S₂₂分别为6×6，6×4，4×4的矩阵，v₁，v₂分别为6行和4行的列向量，那么要求解的方程组能够进一步简化为 

$\left\{\begin{array}{c}(S_{11}-{\mathrm{\λC}}_1)v_1+S_{12}{v}_2=0\\ S_{12}^T{v}_1+S_{22}{v}_2=0\end{array}\right.$

因为待拟合的髋臼凹面数据采样点集是不共面的，所以S₂₂是可逆的，因此解方程组可以得到 

$v_2=-S_{22}^{-1}{S}_{12}^T{v}_1$

代入方程组，从而得到 

$(S_{11}-S_{12}{S}_{22}^{-1}{S}_{12}^T)v_1={\mathrm{\λC}}_1{v}_1$

在实际求解过程中，取k=4，这样能够保证

是可逆的。而当C₁是可逆的时候，v₁即为唯一一个正项特征值对应的特征向量。同时将v₁的解带入式

即可求解前述一般二次曲面方程中的所有未知系数。 

下面将椭球的一般二次曲面方程化归为标准形式，设X=(x,y,z)^T，那么一般二次曲面方程表示为矩阵形式为 

d+(2p,2q,2r)X+X^TKX=0 

对于矩阵K，使用Jacobi正交变换法求其特征根(λ₁,λ₂,λ₃)，及3×3的特征向量R，于是有 

$R^T\mathrm{KR}=\left(\begin{array}{ccc}{\mathrm{\λ}}_1& & \\ & {\mathrm{\λ}}_2& \\ & & {\mathrm{\λ}}_3\end{array}\right)$

同时设Y=R^TX，连同上式，代入d+(2p,2q,2r)X+X^TKX=0，可得 

$d+(a_1,a_2,a_3)Y+Y^T\left(\begin{array}{ccc}{\mathrm{\λ}}_1& & \\ & {\mathrm{\λ}}_2& \\ & & {\mathrm{\λ}}_3\end{array}\right)Y=0$

(a₁,a₂,a₃)=(2p,2q,2r)R 

由于R为特征向量，坐标X，Y之间的关系为 

$X=Y+{(\frac{a_1}{{2\mathrm{\λ}}_1},\frac{a_2}{{2\mathrm{\λ}}_2},\frac{a_3}{{2\mathrm{\λ}}_3})}^T$

代入 $d+(a_1,a_2,a_3)Y+Y^T\left(\begin{array}{ccc}{\mathrm{\λ}}_1& & \\ & {\mathrm{\λ}}_2& \\ & & {\mathrm{\λ}}_3\end{array}\right)Y=0,$ 可得 

$a_0+X^T\left(\begin{array}{ccc}{\mathrm{\λ}}_1& & \\ & {\mathrm{\λ}}_2& \\ & & {\mathrm{\λ}}_3\end{array}\right)X=0$

$a_0=d-\frac{1}{4}(\frac{{a_1}^2}{{\mathrm{\λ}}_1}+\frac{{a_2}^2}{{\mathrm{\λ}}_2}+\frac{{a_3}^2}{{\mathrm{\λ}}_3})$

与椭球标准二次曲面的矩阵方程比较，得到椭球的三个主轴长度为 

$a=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_1}},$ $b=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_2}},$ $c=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_3}}$

设椭球三个轴上的旋转角为(α,β,γ)，得到椭球面中心坐标系与主坐标系的变换关系为 

X′=X₀+R₁(α)R₁(β)R₁(γ)X 

$X_0={(\frac{a_1}{{2\mathrm{\λ}}_1},\frac{a_2}{{2\mathrm{\λ}}_2},\frac{a_3}{{2\mathrm{\λ}}_3})}^T$

R₁(α)R₁(β)R₁(γ)=R^T

据此，计算得到了目标椭球面的三个轴长、中心坐标以及方位角。 

本发明所述的一种面向个性化设计的股骨头表面三维模型重构方法，在普通个人电脑上即可实施，无需特殊设置的支持，这使得本发明很容易推广。以上所述实施例仅表达了本发明的实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。 

Claims (1)

1.一种面向个性化设计的股骨头表面三维模型重构方法，其特征包括以下步骤：

①输入由医用CT机所采集的髋关节CT序列图像；

②对已输入的序列图像中每一幅二维CT图像内的股骨头和髋臼进行轮廓提取与内部填充；先利用Sobel算子提取二维CT图像内股骨头和髋臼的轮廓，再利用扫描线法实现股骨头和髋臼轮廓的内部填充；

③在已填充的股骨头和髋臼中，利用像素间的连通性实现每一幅二维CT图像内髋臼的分割操作，具体步骤如下：

1）在股骨头和髋臼封闭区域内任取两个像素，如果这两个像素有连通的路径，则将这两个像素归为同一类，如果这两个像素没有连通的路径，则将这两个像素归为两类；

2）遍历股骨头和髋臼封闭区域内的所有像素，将所有像素归为两类；

3）舍弃属于股骨头封闭区域内的像素集合，保留属于髋臼封闭区域内的像素集合，实现髋臼的分割操作；

④利用已分割完毕的髋臼序列图像，进行髋臼三维模型重构；

⑤在已得到的髋臼三维模型上，利用空间辐射线法，进行髋臼凹面数据点采样，包括如下步骤：

1）在髋臼三维模型的凹面内选择一个空间坐标点，以此坐标点为中心，向全周空间发射辐射直线；

2）每条辐射直线与髋臼三维模型相交时即获得一个交点，所有交点的集合即为髋臼凹面数据采样点；

⑥利用已得到的髋臼凹面数据点，采用基于一般二次曲面约束的椭球面拟合方法，实现个性化股骨头表面重构，具体步骤如下：

1）将一般二次曲面约束为椭球曲面，设一般二次曲面的方程为：

ax²+by²+cz²+2fyx+2gxz+2hxy+2px+2qy+2rz+d=0    （1）

设：

I=a+b+c    （2）

J=ab+bc+ac-f²-g²-h    （3）

$K=\left|\begin{array}{ccc}a& h& g\\ h& b& f\\ g& f& c\end{array}\right|---\left(4\right)$

对于任意一个椭球面都有α≥4，使得αJ-I²>0，选用α=4来将二次曲面约束为椭球曲面；

2）利用髋臼凹面数据采样点，采用基于代数距离的最小二乘法进行椭球曲面拟合，获得公式（1）中的所有未知系数；

3）设X=(x,y,z)^T，将椭球面一般二次曲面方程化归为矩阵形式：

d+(2p,2q,2r)X+X^TKX=0    （5）

对于矩阵K，使用Jacobi正交变换法求其特征根(λ₁,λ₂,λ₃)，及3×3的特征向量R，令

(a₁,a₂,a₃)=(2p,2q,2r)R    （6）

有

$a_0=d-\frac{1}{4}(\frac{{a_1}^2}{{\mathrm{\λ}}_1}+\frac{{a_2}^2}{{\mathrm{\λ}}_2}+\frac{{a_3}^2}{{\mathrm{\λ}}_3})---\left(7\right)$

得到椭球面的三个主轴的长度

$a=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_1}}$ $b=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_2}}$ $c=\sqrt{-\frac{a_0}{{\mathrm{\λ}}_3}}---\left(8\right)$

设椭球面三个主轴的旋转角为(α,β,γ)，得到椭球面中心坐标系与主坐标系的变换关系

X′=X₀+R₁(α)R₁(β)R₁(γ)X    （9）

$X_0={(\frac{a_1}{{2\mathrm{\λ}}_1},\frac{a_2}{{2\mathrm{\λ}}_2},\frac{a_3}{{2\mathrm{\λ}}_3})}^T---\left(10\right)$

R₁(α)R₁(β)R₁(γ)=R^T    （11）。

Images