CN102208117A - 脊椎的三维几何与有限元混合模型的构建方法 - Google Patents

Abstract

本发明提出一种脊椎的三维几何与有限元混合模型的构建方法，属医学图像处理技术领域。构建过程为：输入脊椎的CT图像；对CT图像进行三维重建和三维切割，得到脊椎三维图像集；建立三维几何形态统计模型，首先对脊椎的特征点定义和手动标定，再对脊椎图像进行对齐和配准，然后对样本集进行训练，得到统计模型；生成有限元模型，将统计模型进行导入，首先生成表面网格模型，然后生成体网格模型；该模型可以直接导入到有限元分析软件中进行生物力学分析。本发明既能够精确描述脊椎的几何外形，又保证有限元分析结果的准确性，提高脊椎模型的精确度，使用方便，便于对脊椎部位的形状和受力进行科学计量。可用于外科医学领域中与脊柱和脊椎相关的研究。

Description

脊椎的三维几何与有限元混合模型的构建方法

技术领域

本发明属于医学图像处理技术领域，主要涉及一种脊椎模型的构建方法，具体是一种脊椎的三维几何与有限元混合模型的构建方法，可用于外科医学领域中与脊柱和脊椎相关的研究。

背景技术

目前，与脊椎相关的疾病在我国已成为常见病和高发病，虽然脊椎外科手术近年得到了很大的发展，但是由于脊椎手术本身特点及脊椎的结构复杂，手术难度和危险性很高。

随着计算机图形图像技术在医学各领域内应用的拓展，医学影像设备的不断升级，空间三维定位系统、计算机医学图像处理及三维可视化技术已经深入到医学的各个领域虽然三维医学图像重建可以弥补二维影像的不足，但精确的立体定位和配准仍是未解决的难题。目前手术导航和虚拟手术技术的总体水平尚不能很好地满足临床需求，问题的关键是高精度建模、准确定位与模拟等关键技术尚未很好解决。因此，需要科学而准确的度量数据。

目前，我国还没有脊柱以及各部分脊椎的三维几何统计形态模型，还没建立脊柱和脊椎三维几何数据库，缺乏不同年龄段脊柱的准确度量数据，给脊椎和脊柱的相关医学研究、教学和计算机辅助骨科的精确定位造成了不便。

在Cootes提出主动形状模型(ASM，Active Sharpe Model)后，对利用ASM方法来建立样本的统计模型有了很大的发展，并且有很多人将其拓展到三维形态上。现阶段，采用ASM方法建立的颅骨、颅面模型的技术很成熟，并在颅面复原等方面获得了较好的效果。

ASM方法是基于点分布的，它需要对每一个样本图像的特征点进行标定。不同的特征点标定方法，直接影响到模型建立的准确性。Cootes提出标记点可分为3类：第1类是与模型直接相关的特殊点；第2类是与模型的形状不直接相关的点；第3类是前两类点的中间插值点。一般情况下，选择第1类和第2类点作为标记点。

其中，由于目前我国还没有建立脊柱和脊椎三维几何数据库，缺乏脊椎的准确度量数据，对脊椎的形态分析缺乏科学度量参数。因此，需要建立脊椎的三维几何形态模型。在采用ASM与ICP方法建立的脊椎三维几何形状模型在实际的应用中只能对脊椎的几何形态进行分析，缺乏脊椎软组织与骨骼直接相互作用力，无法对脊椎部位的受力情况进行分析，也不能对脊椎的各个部位和关键点进行科学度量。因此，需要将该统计模型的方法与有限元模型相结合，从而既能描述脊椎的几何形态，又能对脊椎部位的受力情况进行定量分析。

发明的内容

本发明的目的在于克服上述技术存在的缺点和不足，提出一种脊椎的三维几何与有限元混合模型的构建方法，该方法既能够精确描述脊椎的几何外形，又保证有限元分析结果的准确性，为脊椎和脊柱的相关医学研究、教学和计算机辅助骨科的精确定位提供一种科学依据和参考。

本发明的技术方案是：首先对三维重建后的医学CT图像进行特征点定位和标记，从而获得每个脊椎样本的形状矩阵；采用最近点迭代算法(ICP，Iterative Closest Point)算法对样本集中形状矩阵进行对齐和配准；再采用主成分分析法(PCA，Principal Component Analysis)方法对样本进行训练，从而建立脊椎的几何形态统计模型；最后，将模型导入到相关软件生成体网格模型，该网格模型可以直接用来进行有限元分析。

本发明是一种脊椎的三维几何与有限元混合模型的构建方法，其特征在于：包括如下过程：

(1)对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，并保留单个的椎骨模型，经大量处理后建立脊椎三维图像的集合；

(2)对脊椎三维图像的特征点进行定义和手动标定，将每一个脊椎样本标定得到的形状矩阵保存起来，形成形状矩阵样本集；

(3)从样本集中选取一个模型样本，并将剩下每一个的样本都作为待配准样本；

(4)采用ICP方法对模型样本与待配准样本的形状矩阵进行旋转变换和平移变换，通过最小二乘法进行迭代，使二者的对应点之间的距离最小，使两个样本在同一个坐标系中具有可比性，完成模型样本与待配准样本的对齐和配准；

(5)采用PCA方法对样本集进行训练，去除原始数据的线性相关性，从而将原有数据的高维度转化为少数的几个主要变换变量，得到样本变化模型；

(6)将得到的脊椎三维几何形态统计模型网格化，生成有限元分析模型。

主动形状模型(Active Shape Model)最初是由T.F.Cootes等人于1992年提出的，其是一种以训练和统计分析为基础的模型方法。德国马普研究会的Volker Blanz和Thoms Vetter提出基于三维特征点标记的颅面统计模型。目前，运用统计学的原理方法和计算机技术对脊椎几何形状信息的研究方面，国内外相关领域还没有比较成熟的方法。本发明在参考ASM的基本原则和其扩展模型在对三维模型处理第三维信息(也就是坐标系中Z轴方向的坐标值)的方法和统计学原理对脊椎进行建模，得到脊椎在三维空间上的统计知识。它通过对样本标定和训练，使用统计分析的方法建立关于目标平均形状及形变模式的先验模型。它能根据训练集中的形变模式对模型进行灵活形状调整，保证模型在合理的范围变化。

有限元法的基本思想是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的近似解，然后推导求解这个域的总的满足条件，从而得到问题的近似解。本发明将得到的三维几何形态统计模型进行体网格化，构建了脊椎的三维几何与有限元混合模型，可以分别进行详细有限元计算，对脊椎进行受力分析时就可以从每一个元素的受力得到整个脊椎的受力情况，不仅可以观察骨形态细微结构，且可通过随意旋转从任意角度、任意断面显示病变的部位和程度，并可以对不同的部位的受力情况进行详细科学计量，从而为脊椎进行生物力学的评估、统计分析提供科学依据。

本发明的实现还在于：过程(2)中对特征点进行定义和手动标定按照如下过程进行：

参考颅骨特征点的定义方法，根据Cootes提出的特征点分类并结合脊椎的特点，对样本集中脊椎样本的特征点进行定义和手动标定，并将每一个样本标定得到的形状矩阵保存起来，形成形状矩阵集合，定义68个脊椎特征点，每个样本用这68个三维特征点构成的形状矩阵来表示：

$X=\left[\begin{array}{ccccc}{x}_1& x_2& ...& x_i& ...x_n\\ y_1& y_2& ...& y_i& ...y_n\\ z_1& z_2& ...& z_i& ...z_n\end{array}\right]$ {i＝1，2，，n}。

本发明将脊椎的三维图像采用形状矩阵进行保存，把一个三维脊椎图像转化成形状矩阵的方式进行保存，减少了对脊椎图像处理的复杂性，便于对该脊椎图像进行程序处理和分析。

本发明的实现还在于：过程(3)选取模型样本和确定待配准样本按照如下过程进行：

假设在样本集中有标定后的m个脊椎样本，每个脊椎样本的形状矩阵M为：

$M_j=\left[\begin{array}{ccccc}{x}_1& x_2& ...& x_i& ...x_n\\ y_1& y_2& ...& y_i& ...y_n\\ z_1& z_2& ...& z_i& ...z_n\end{array}\right]$ (j＝1，2，...，m)

(i＝1，2，...，n)

首先从样本集中m个脊椎样本选取一个比较理想的，大小、位置都比较标准的样本作为模型样本M_q(1≤q≤m)，然后从其余的m-1个样本中任选一个样本作为待配准样本M_p(1≤p≤m；p≠q)；

为了使这两个样本能够配准起来，就要找到最优的旋转矩阵R和平移向量U，使对齐后的M_p与M_q具有可比性，其中

M′_p＝R*M_p+U，

而旋转矩阵R可以用围绕X，Y，Z坐标轴的连续旋转角度θ_x，θ_y，θ_z表示：

$R=\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_z& {\sin \mathrm{\θ}}_z& 0\\ {-\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_z& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_y& 0& -{\sin \mathrm{\θ}}_y\\ 0& 1& 0\\ {\sin \mathrm{\θ}}_y& 0& {\cos \mathrm{\θ}}_y\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_x& {\sin \mathrm{\θ}}_x\\ 0& -{\sin \mathrm{\θ}}_x& {\cos \mathrm{\θ}}_x\end{array}\right).$

为了方便描述，令Q＝M_q，P＝M_p，则有：

模型样本Q中三维点集为：

$\left\{Q_i\right|Q_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\},$

待配准样本P中三维点集为：

$\left\{P_i\right|P_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\}.$

本发明从样本集中选取一个模型样本，并从剩下的样本中待配准样本，以便将待配准样本进行旋转变换和平移变换。

本发明的实现还在于：其中采用ICP方法完成对齐和配准按照如下过程进行：

4.1初始化迭代次数k＝0；

4.2从待配准样本P中取点

4.3根据约束条件

计算出模型样本Q中的与

对应点

4.4根据约束条件

计算旋转矩阵R；

4.5平移向量U是两个点集之间的重心差异，通过计算R*P^k与Q^k之间的重心差异即得；

4.6第k次迭代配准得到的点集表示为：

$P^{k+1}=\left\{P_i^{k+1}\right|P_i^{k+1}=R^k*P_i^k+U^k,P_i^k\∈P\};$

4.7计算两个点集之间的距离：

$d^{k+1}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|P_i^{k+1}-Q_i^k\left|\right|}^2;$

如果d^k+1小于初始预先给定的匹配距离精度τ或k大于最大迭代次数，终止，否则跳转到4.2，进行下一次迭代。

通过ICP的每一次迭代，使得点集P^k与Q^k之间相对应点的距离缩小，最终得到与模型样本最为接近的匹配样本：

M_p′＝P^k+1。

ICP方法是由Besl和Mckey于1992年提出的。由于ICP配准方法在配准前不需要知道标定点之间的对应关系并且配准精度高，在颅骨的配准中比较普遍和成熟，因此本发明也采用ICP方法进行图像的对齐和配准。ICP的基本思想是从待配准样本到模型样本之间的一个最优的匹配，就是对两个形状矩阵进行旋转变换和平移变换，通过最小二乘法进行迭代，使二者的对应点之间的距离最小。

本发明采用ICP方法对模型样本和待配准样本进行对齐和配准，通过多次迭代使二者之间的距离最小，从而使待配准样本与模型样本在同一坐标系进行比较，这样得到的脊椎三维几何形状统计模型具备较高的精度，能够进行准确的定位，从而更加具有说服力。

本发明的实现还在于：其中过程(5)建立样本变化矩阵按照如下过程进行：

5.1将三维点集M转化为一维的形状向量X：

X_j＝{x₁，y₁，z₁，x₂，y₂，z₂，...，x_i，y_i，z_i...，x_n，y_n，z_n}

；

(j＝1，2，...，m)(i＝1，2，...，n)

5.2计算m个样本的平均模型：

$\stackrel{\‾}{X}=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{X}_i;$

计算每个样本相对于平均样本

的偏移量

$d_{X_i}=X_i-\stackrel{\‾}{X};$

5.3计算样本集的协方差矩阵S：

$S=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{d}_{X_i}{d}_{X^T};$

5.4计算协方差矩阵S的特征值λ和特征向量Φ：

SΦ_k＝λ_kΦ_k；

5.5从计算得到的协方差矩阵S提取前t个特征值b(大小为t*1)，及其对应的t个特征向量组成的矩阵Ψ(大小为3n*t)。从而得到该脊椎样本的形状模型：

$X\≈\stackrel{\‾}{X}+\mathrm{\Ψ}*b$

Ψ＝(Φ₁，Φ₂，...，Φ_t)。

b＝(λ₁，λ₂，...，λ_t)^T

PCA方法的核心思想是利用数学变换的方法，去除原始数据之间的线性相关性，从而将原有数据的高维数转化为少数的几个主要变换变量，是在人脸识别等多个领域广泛使用的一种经典的统计学方法。

本发明将原有数据的高维数转化为少数的几个主要变换变量，得到样本变化模型。通过该样本变化模型，能够将任一脊椎样本转化成为主要的变化向量，既能减少了描述脊椎样本所需要的变量数目，也能保证该脊椎样本的几何形态不发生大的改变。

本发明的实现还在于：其中过程(6)按照如下过程进行：

6.1在软件Simpleware的ScanIP模块的三维视图中，在3D Preview中选择FE，生成表面网格模型；

6.2从ScanIP模块中导出以.sfh为扩展名的文件，再将其导入到ScanFE模块，在Control Panel中点击Mesh选项，在Mesh Options中选择Smoothed，点击Apply，将会得到一个平滑的体网格模型；

6.3在体网格模型的基础上建立有限元分析模型。

本发明的实现还在于：建立有限元分析模型按照如下过程进行：

在体网格导出之前，选择Export选项中的Finite Element Model，然后选择ANSYS选项，在弹出的对话框中点击Configure part materials设置合理的Mass Density，Young’s Modulus和Passion’s Ratio参数，以保证有限元分析结果的准确性。导出的模型可以被直接导入到ANSYS软件中进行有限元分析。

本发明在建立好三维几何形态统计模型之后，就可以对该模型进行网格化了，首先生成的是表面网格模型。根据后续的处理软件如何导入这些值，在simpleware软件的ScanIP模块中合理设置最大密度、杨氏模量和泊松比值。从ScanIP模块导出表面网格模型，将其导入到ScanFE模块，合理设置相关参数就可以生成较为理想的体网格模型。

本发明与现有技术相比具有以下优点：

1.提出一种三维几何与有限元混合模型的构建方法

本发明重建出脊椎的三维统计模型，再转换为有限元网格，既能够精确描述脊椎的几何外形，也能保证有限元分析结果的准确性，使用方便，便于对脊椎部位的形状和受力进行科学计量。

2.在三维重建之后，特征点标定之前首先进行图像分割

由于脊椎的每一块椎骨各异，而且比较难统一，在处理的过程中，首先将脊椎图像的每一块椎骨进行分割后分别进行建模，从而提高脊椎模型的准确性。

3.脊椎特征点的定义和标定，提高了三维几何形态统计模型的准确性

由于脊椎结构的复杂性，在实际中对其特征点的标定还没有统一的标准，国内外医学界也没有脊椎具体的特征点的相关规范和标准。但是在颅脑的特征点标定方面，国际上有一些得到广泛认同的标准。本发明在根据Cootes提出的标记点分类和在应用得到认可的MPEG人脸标记方法的基础上，再结合医学上脊椎的一些常见特征，如椎管相对于椎突在实际外科手术中重要程度高，给出脊椎采样特征点的定义和标定方法。

4.将颅骨的三维几何建模方法推广到脊椎，并对模型进行了改进

目前，我国还缺乏脊椎的三维几何统计形态模型，相关的统计数据库尚未建立，但是采用基于点分布的ASM方法建立的颅骨、颅面模型的技术已经很成熟，并且在颅面复原等方面获得了较好的效果。本发明将这种方法推广到脊椎，对脊椎的第一块腰椎骨建立三维几何统计模型，并针对脊椎的特点对模型进行了改进。

附图说明

图1是本发明模型构建流程图；

图2是三维切割前后的腰椎的三维图像，其中图2(a)是未切割处理的腰椎部分的三维图像，图2(b)是经过切割处理之后的腰椎三维图像；

图3是本发明依据骨骼解剖学中的腰椎右侧面观和上面观的彩绘图像进行特征点标记所得的图像，其中图3(a)是腰椎彩绘图像的右侧面观的特征点标记图像，图3(b)是腰椎彩绘图像的上侧面观的特征点标记图像；

图4是本发明腰椎三维图像的右侧面观和上面观手动标记图像，其中图4(a)是腰椎三维图像的右侧面观的手动标记点图像，图4(b)是腰椎三维图像的上侧面观的手动标记点图像；

图5是本发明配准前后的模型和待配准数据图像，其中图5(a)是腰椎三维图像配准前的模型数据与待配准数据图像的点云示意图，图5(b)是腰椎三维图像配准后的模型数据与待配准数据图像的点云示意图；

图6是本发明腰椎表面网格模型左侧面观和上面观图像，其中图6(a)是腰椎有限元表面网格模型左侧面观图像，图6(b)是腰椎有限元表面网格模型上面观图像；

图7是本发明腰椎体网格模型左侧面观和上面观图像，其中图7(a)是腰椎有限元体网格模型左侧面观图像，图7(b)是腰椎有限元体网格模型上面观图像。

具体实施方式

实施例1.

本发明是一种脊椎的三维几何与有限元混合模型的构建方法，参见图1，脊椎的三维几何与有限元混合模型的构建方法包括如下过程：

过程1，，对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，对该三维模型进行切割并保留单个的腰椎部位的三维模型；通过对大量人体脊椎图像处理，建立脊椎三维图像的集合；为建立脊椎三维几何形态统计模型和有限元分析奠定基础出。

过程2，对脊椎三维图像的特征点进行定义和手动标定，将每一个脊椎样本标定得到的形状矩阵保存起来，形成形状矩阵样本集；从该步骤开始建立脊椎的三维几何形态统计模型；

在实际中对脊椎图像的标定很难有一个统一的标准，国内外医学界也没有脊椎具体的特征点的相关规范和标准。但是在颅脑的特征点标定方面，国际上有一些得到广泛认同的标准，如MPEG-4(Moving Picture Experts Group)的人脸定义参数(FDP)对于人脸特征点的定义。

本发明在根据Cootes提出的标记点分类和在应用得到认可的MPEG人脸标记方法的基础上，再结合医学上脊椎的一些常见特征，如椎管、椎突和刺突等，给出脊椎采样特征点的定义和标定，本发明重点在脊椎腰骶段选取了能够描绘其主要特征的位置点，总共68个三维特征点。

过程3，从样本集中选取一个模型样本，并将剩下每一个的样本都作为待配准样本；

假设在样本集中有标定后的m个脊椎样本，每个脊椎样本的形状矩阵为：

$M_j=\left[\begin{array}{ccccc}{x}_1& x_2& ...& x_i& ...x_n\\ y_1& y_2& ...& y_i& ...y_n\\ z_1& z_2& ...& z_i& ...z_n\end{array}\right]$ (j＝1，2，...，m)

(i＝1，2，...，n)。

首先从样本集中m个脊椎样本选取一个比较理想的，大小、位置都比较标准的样本M_q(1≤q≤m)作为模型样本，然后从其余的m-1个样本中任选一个样本M_p  (1≤p≤m；p≠q)作为待配准样本。

为了使这两个样本能够配准起来，就要找到最优的旋转矩阵R和平移向量U，使对齐后的M_p与M_q具有可比性，其中

M′_p＝R*M_p+U，

而旋转矩阵R可以用围绕X，Y，Z坐标轴的连续旋转角度θ_x，θ_y，θ_z表示：

$R=\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_z& {\sin \mathrm{\θ}}_z& 0\\ {-\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_z& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_y& 0& -{\sin \mathrm{\θ}}_y\\ 0& 1& 0\\ {\sin \mathrm{\θ}}_y& 0& {\cos \mathrm{\θ}}_y\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_x& {\sin \mathrm{\θ}}_x\\ 0& -{\sin \mathrm{\θ}}_x& {\cos \mathrm{\θ}}_x\end{array}\right).$

为了方便描述，令Q＝M_q，P＝M_p，则有：

模型样本Q中三维点集为：

$\left\{Q_i\right|Q_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\},$

待配准样本P中三维点集为：

$\left\{P_i\right|P_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\}.$

过程4，采用ICP方法对模型样本与待配准样本的形状矩阵进行旋转变换和平移变换，通过最小二乘法进行迭代，使二者的对应点之间的距离最小，使两个样本在同一个坐标系中具有可比性，完成模型样本与待配准样本的对齐和配准；

ICP算法的配准过程如下：

4.1初始化迭代次数k＝0；

4.2从待配准样本P中取点

4.3根据约束条件

计算出模型样本Q中的与对应点

4.4根据约束条件

计算旋转矩阵R；

4.5平移向量U是两个点集之间的重心差异，通过计算R*P^k与Q^k之间的重心差异即得；

4.6第k次迭代配准得到的点集表示为：

$P^{k+1}=\left\{P_i^{k+1}\right|P_i^{k+1}=R^k*P_i^k+U^k,P_i^k\∈P\};$

4.7计算两个点集之间的距离：

$d^{k+1}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|P_i^{k+1}-Q_i^k\left|\right|}^2;$

如果d^k+1小于初始预先给定的匹配距离精度τ或k大于最大迭代次数，终止，否则跳转到4.2，进行下一次迭代。

通过ICP的每一次迭代，使得点集P^k与Q^k之间相对应点的距离缩小，最终得到与模型样本最为接近的匹配样本：

M_p′＝P^k+1。

过程5，采用PCA方法对样本集进行训练，去除原始数据的线性相关性，从而将原有数据的高维度转化为少数的几个主要变换变量，得到样本变化模型；

对于配准后的m个样本：

5.1将三维点集M转化为一维的形状向量X：

X_j＝{x₁，y₁，z₁，x₂，y₂，z₂，...，x_i，y_i，z_i...，x_n，y_n，z_n}

；

(j＝1，2，...，m)(i＝1，2，...，n)

5.2计算m个样本的平均模型：

$\stackrel{\‾}{X}=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{X}_i$

计算每个样本相对于平均样本

的偏移量

$d_{X_i}=X_i-\stackrel{\‾}{X};$

5.3计算样本集的协方差矩阵S：

$S=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{d}_{X_i}{d}_{X^T};$

5.4计算协方差矩阵S的特征值λ和特征向量Φ：

SΦ_k＝λ_kΦ_k；

5.5从计算得到的协方差矩阵S提取前t个特征值b(大小为t*1)，及其对应的t个特征向量组成的矩阵Ψ(大小为3n*t)。从而得到该脊椎样本的形状模型：

$X\≈\stackrel{\‾}{X}+\mathrm{\Ψ}*b$

Ψ＝(Φ₁，Φ₂，...，Φ_t)

b＝(λ₁，λ₂，...，λ_t)^T

本发明将68个三维点组成的脊椎样本形状矩阵(维度为：68*3)转换成为少数的几个主要变换变量(维度为：t*1)，既能简化脊椎图像的描述，又能较好的保持脊椎的几何形态。

到此，本发明已经成功建立了脊椎的三维几何形态统计模型。

采用ASM与ICP方法建立的脊椎三维几何形状模型在实际的应用中只能对脊椎的几何形态进行分析，缺乏脊椎软组织与骨骼直接相互作用力，无法对脊椎部位的受力情况进行分析。因此，需要将该统计模型的方法与有限元模型相结合，从而既能描述脊椎的几何形态，又能对脊椎部位的受力情况进行定量分析。

过程6，将得到的脊椎三维几何形态统计模型网格化，生成有限元的表面网格模型和体网格模型。

在建立好三维几何形态统计模型之后，对该模型进行网格化，首先生成表面网格模型。

在软件Simpleware的ScanIP模块的三维视图中，合理设置最大密度、杨氏模量和泊松比值，在3D Preview中选择FE，生成表面网格模型。

综上所述，本发明在建立精确的脊椎三维几何形态统计模型的基础上，对其进行了网格化，所生成的体网格模型能有直接导入到有限元分析软件中进行有限元分析，本模型的构建方法既能够精确描述脊椎的几何外形，又保证有限元分析结果的准确性，便于对脊椎部位的形状和受力进行科学计量。从而对脊柱和脊椎进行生物力学评估。

实施例2

脊椎的三维几何与有限元混合模型的构建方法同实施例1，其中过程1中是对脊椎的三维模型进行切割并保留单个的腰椎部位的三维模型，图2(a)显示的是重建之后分割之前的脊椎的三块相邻的椎骨，可以看出脊椎的各个椎骨形态各异，而且与肋骨联系紧密，即使是相邻的两节椎骨也有很大差异，会在进行特征点标记和模型匹配方面造成了很大的麻烦。因此，在处理的过程中，需要对脊椎图像的每一块椎骨进行分割后，分别进行建模，从而提高脊椎模型的准确性。图2(b)显示的是三维切割后得到的三维腰椎的第一块骨。通过对大量人体脊椎图像处理，建立脊椎三维图像的集合。

过程2是对脊椎采样特征点进行定义和标定，图3(a)和图3(b)是本发明对腰椎三维图像特征点的定义示意图，图3(a)是腰椎彩绘图像的右侧面观的特征点标记图像，图3(b)是腰椎彩绘图像的上侧面观的特征点标记图像。在考虑到脊椎部位的生物力学科学计量的方便，本发明重点在椎管部分的标记点比较多，能够保证椎管部位的精确性。另外，在对三维点的第三维信息，也就是坐标系中Z轴方向的坐标值进行了平均值的修订，从而能够进行较好的配准和对齐。这是本发明针对脊椎特点进行的两个改进。其中68个特征点的定义见表1和表2：

表1腰椎上面观特征点定义

表2腰椎右侧面观特征点定义

由于这些的特征点分布不均匀而且脊椎图像没有一个很确切的边界，在标记软件系统中对图像进行标定时还可以采用手动标定的方法进行辅助标定，这样能够保证模型的准确性，并且简单易行。其中图4(a)是腰椎三维图像的右侧面观的手动标记点图像，图4(b)是腰椎三维图像的上侧面观的手动标记点图像。

其中过程4是对模型样本与待配准样本的对齐和配准，同实施例1，其中图5(a)是腰椎三维图像配准前的模型数据与待配准数据图像的点云示意图，图5(b)是腰椎三维图像配准后的模型数据与待配准数据图像的点云示意图。如图中可见，虽然图5反映是配准前后的三维点的分布示意图，但是可以从图5(b)中看出，三维点的分布比较合理，两个三维图像之间的距离相对较小，较好的实现了待配准数据与模型数据对齐与配准。

其中过程5同实施例1。

其中过程6为将三维几何形态模型生成为有限元体网格模型。在体网格导出之前，选择Export选项中的Finite Element Model，然后选择ANSYS选项，在弹出的对话框中点击Configure part materials设置合理的Mass Density，Young’s Modulus和Passion’s Ratio参数，以保证有限元分析结果的准确性。将从ScanIP模块导出表面网格模型导入到ScanFE模块，合理设置相关参数就可以生成较为理想的体网格模型。

为了保证有限元分析结果的准确性，需要精确的描述脊椎的材料特性，因此在体网格导出之前还要对皮质骨、松质骨、纤维环、髓核等赋予不同的材料属性。导出的网格模型可以直接导入到ANSYS等有限元分析软件中进行有限元分析，同时可以施加负载来对腰椎的受力情况进行分析。图6(a)是腰椎有限元表面网格模型左侧面观图像，图6(b)是腰椎有限元表面网格模型上面观图像；图7(a)是腰椎有限元体网格模型左侧面观图像，图7(b)是腰椎有限元体网格模型上面观图像。从图6和图7中可以看出，该体网格模型既精确描述脊椎的几何外形，又便于对脊椎部位的形状和受力进行分析和度量。

综上所述，本发明提出的是一种脊椎的三维几何与有限元混合模型的构建方法，精确的三维几何形态统计模型为之后有限元分析的准确性奠定了基础。

实施例3

脊椎的三维几何与有限元混合模型的构建方法同实施例1-2，其中过程5得到三维几何形态统计模型如下：

$X\≈\stackrel{\‾}{X}+\mathrm{\Ψ}*b$

Ψ＝(Φ₁，Φ₂，...，Φ_t)

b＝(λ₁，λ₂，...，λ_t)^T

ψ可以表示该样本集中的主要变化模式。对模型的准确性评价可以由ψ对样本集的整体变化模式的解释度来进行衡量。解释度越高说明模型越准确，能够表现整个样本的形状变化情况。

设λ_T为所有特征值λ_k的总和，

则用λ_k/λ_T代表某个模式对整体变化模式的贡献率，ψ对整体变化模式的贡献率可以由它所对应的t个模式对整体模式贡献率的总和来计算。这里表示为解释度α，

$\mathrm{\α}=\frac{\underset{k=1}{\overset{t}{\mathrm{\Σ}}}{\mathrm{\λ}}_k}{{\mathrm{\λ}}_T},0\≤\mathrm{\α}\≤1.$

一般情况下，解释度要求达到95％以上。在实际应用中，可以改变通过改变选取模式个数t来改变解释度α。

三维几何形态统计模型建立好之后还需要具备良好的可用性，才能够在脊椎图像的科学计量以及外科手术中得到应用。

根据某个脊椎样本B计算其主分量参数b：

$b={\mathrm{\Ψ}}^T*(B-\stackrel{\‾}{X})$

通过修改主分量参数b来获得新脊椎数据B：

$B=\stackrel{\‾}{X}+\mathrm{\Ψ}*b$

本发明构建的脊椎三维几何形态统计模型具备较好的准确性和可用性。从而保证了脊椎的三维几何与有限元混合模型的精确度。

Claims (7)

1.一种脊椎的三维几何与有限元混合模型的构建方法，其特征在于：包括如下步骤：

(1)对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，对该三维模型进行切割并保留单个的腰椎部位的三维模型；通过对大量人体脊椎图像处理，建立脊椎三维图像的集合；

(2)对脊椎三维图像的特征点进行定义和手动标定，将每一个脊椎样本标定得到的形状矩阵保存起来，形成形状矩阵样本集；

(3)从形状矩阵样本集中选取一个模型样本，并将剩下每一个的样本都作为待配准样本；

(4)采用ICP方法对模型样本与待配准样本的形状矩阵进行旋转变换和平移变换，通过最小二乘法进行迭代，使二者的对应点之间的距离最小，使两个样本在同一个坐标系中具有可比性，完成模型样本与待配准样本的对齐和配准；

(5)采用PCA方法对形状矩阵样本集进行训练，去除原始数据的线性相关性，从而将原有数据的高维度转化为少数的几个主要变换变量，得到样本变化模型；

(6)将得到的脊椎三维几何形态统计模型网格化，生成有限元的表面网格模型和体网格模型。

2.根据权利要求1所述的脊椎的三维几何与有限元混合模型的构建方法，其特征在于：过程(2)中所述对特征点进行定义和手动标定，按照如下过程进行：

参考颅骨特征点的定义方法，根据Cootes提出的特征点分类并结合脊椎的特点，定义68个脊椎特征点，每个样本用这68个三维特征点构成的形状矩阵来表示：

$X=\left[\begin{array}{ccccc}{x}_1& x_2& ...& x_i& ...x_n\\ y_1& y_2& ...& y_i& ...y_n\\ z_1& z_2& ...& z_i& ...z_n\end{array}\right]$ {i＝1，2，，n}。

3.根据权利要求1所述的脊椎的三维几何与有限元混合模型的构建方法，其特征在于：其中过程(3)按照如下过程进行：

假设在形状矩阵样本集中有标定后的m个脊椎样本，每个脊椎样本的形状矩阵M为：

$M_j=\left[\begin{array}{ccccc}{x}_1& x_2& ...& x_i& ...x_n\\ y_1& y_2& ...& y_i& ...y_n\\ z_1& z_2& ...& z_i& ...z_n\end{array}\right]$ (j＝1，2，...，m)

      (i＝1，2，...，n)

首先从样本集中m个脊椎样本选取一个比较理想的，大小、位置都比较标准的样本作为模型样本M_q(1≤q≤m)，然后从其余的m-1个样本中任选一个样本作为待配准样本M_p(1≤p≤m；p≠q)；

找到最优的旋转矩阵R和平移向量U，使对齐后的M_p′与M_q具有可比性，其中

M′_p＝R*M_p+U，

而旋转矩阵R可以用围绕X，Y，Z坐标轴的连续旋转角度θ_x，θ_y，θ_z表示：

$R=\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_z& {\sin \mathrm{\θ}}_z& 0\\ {-\sin \mathrm{\θ}}_z& {\cos \mathrm{\θ}}_z& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}{\cos \mathrm{\θ}}_y& 0& -{\sin \mathrm{\θ}}_y\\ 0& 1& 0\\ {\sin \mathrm{\θ}}_y& 0& {\cos \mathrm{\θ}}_y\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& {\cos \mathrm{\θ}}_x& {\sin \mathrm{\θ}}_x\\ 0& -{\sin \mathrm{\θ}}_x& {\cos \mathrm{\θ}}_x\end{array}\right),$

令Q＝M_q，P＝M_p，则有：

模型样本Q中三维点集为：

$\left\{Q_i\right|Q_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\},$

待配准样本P中三维点集为：

$\left\{P_i\right|P_i=\left(\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right)i=\mathrm{1,2},...n\}.$

4.根据权利要求1所述的脊椎的三维几何与有限元混合模型的构建方法，其特征在于：过程(4)中所述采用ICP方法完成对齐和配准，按照如下过程进行：

4.1初始化迭代次数k＝0；

4.2从待配准样本P中取点

4.3根据约束条件

计算出模型样本Q中的与

对应点

4.4根据约束条件

计算旋转矩阵R；

4.5平移向量U是两个点集之间的重心差异，通过计算R*P^k与Q^k之间的重心差异即得；

4.6第k次迭代配准得到的点集表示为：

$P^{k+1}=\left\{P_i^{k+1}\right|P_i^{k+1}=R^k*P_i^k+U^k,P_i^k\∈P\};$

4.7计算两个点集之间的距离：

$d^{k+1}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\left|\right|P_i^{k+1}-Q_i^k\left|\right|}^2;$

如果d^k+1小于初始预先给定的匹配距离精度τ或k大于最大迭代次数，终止，否则跳转到4.2，进行下一次迭代；

通过ICP的每一次迭代，使得点集P^k与Q^k之间相对应点的距离缩小，最终得到与模型样本最为接近的匹配样本：

M_p′＝P^k+1。

5.根据权利要求1所述的脊椎的三维几何与有限元混合模型的构建方法，其特征在于：过程(5)中所述建立样本变化模型，按照如下过程进行：

5.1将三维点集M转化为一维的形状向量X：

X_j＝{x₁，y₁，z₁，x₂，y₂，z₂，...，x_i，y_i，z_i...，x_n，y_n，z_n}；

(j＝1，2，...，m)(i＝1，2，...，n)

5.2计算m个样本的平均模型：

$\stackrel{\‾}{X}=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{X}_i;$

计算每个样本相对于平均样本

的偏移量

$d_{X_i}=X_i-\stackrel{\‾}{X};$

5.3计算样本集的协方差矩阵S：

$S=\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{d}_{X_i}{d}_{X^T};$

5.4计算协方差矩阵S的特征值λ和特征向量Φ：

SΦ_k＝λ_kΦ_k；

5.5从计算得到的协方差矩阵S提取前t个特征值b(大小为t*1)，及其对应的t个特征向量组成的矩阵Ψ(大小为3n*t)，从而得到该脊椎样本的形状模型：

$X\≈\stackrel{\‾}{X}+\mathrm{\Ψ}*b$

Ψ＝(Φ₁，Φ₂，...，Φ_t)。

b＝(λ₁，λ₂，...，λ_t)^T

6.根据权利要求1所述的脊椎的三维几何与有限元混合模型的构建方法，其特征在于：其中过程(6)按照如下过程进行：

6.1在软件Simpleware的ScanIP模块的三维视图中，在3D Preview中选择FE，生成表面网格模型；

6.2从ScanIP模块中导出以.sfh为扩展名的文件，再将其导入到ScanFE模块，在Control Panel中点击Mesh选项，在Mesh Options中选择Smoothed，点击Apply，将会得到一个平滑的体网格模型；

6.3在体网格模型的基础上建立有限元分析模型。

7.根据权利要求6所述的建立网格模型及有限元分析模型的构建方法，其特征在于：过程(6)中所述建立有限元分析模型，按照如下过程进行：

在体网格导出之前，选择Export选项中的Finite Element Model，然后选择ANSYS选项，在弹出的对话框中点击Configure part materials设置合理的Mass Density，Young’s Modulus和Passion’s Ratio参数，导出的模型直接导入到ANSYS软件中进行有限元分析。

Images