CN106228567A - 一种基于平均曲率流的脊椎特征点自动识别方法 - Google Patents

Abstract

一种基于平均曲率流的脊椎特征点自动识别方法，首先，对人体脊椎的CT图像进行三维重建，得到所需标注特征点的椎体三维模型，对步骤一得到的椎体模型进行网格划分；其次，手动选取模型上某点，计算在手动选取点周围极小半径为R的球形空间内每个顶点的平均曲率值；最后，选取手动选取点周围平均曲率最大的n个点，对平均曲率最大的n个点，分别与选取点作内积，夹角最小的即为所求的点，能够提高特征点标注的精确度，从而提高数据测量精确度，具有精确度高、易于实现的特点。

Description

一种基于平均曲率流的脊椎特征点自动识别方法

技术领域

本发明属于医学图像处理技术领域，特别涉及一种基于平均曲率流的脊椎特征点自动识别方法。

背景技术

脊椎是人体重要的支撑骨骼，近年来，随着人们的生活和工作方式的改变，脊椎疾病发病率越来越高。由于脊椎结构的复杂性以及脊椎手术本身的特点，目前脊椎手术的难度很高，并且具有很大的危险性。脊椎模型精确的定位和配准依旧是未解决的诊疗难题。

在脊椎模型的定位和配准中，主要任务是特征点的定位。使用纯手动的方法标注特征点，特征点标注不准确，导致测量的数据精确度不高，模型匹配不准确。采用高斯曲率作为选择特征点的参考值的方法，虽然能够在一定程度上提高标注的准确度，但是对于手工选取特征点造成的误差，这种方法的效果有限。因此，对脊椎特征点标注的准确度的研究尤为必要。

发明内容

为了克服上述现有技术的不足，本发明的目的在于提出一种基于平均曲率流的脊椎特征点自动识别方法，能够提高特征点标注的精确度，从而提高数据测量精确度，具有精确度高、易于实现的特点。

为了实现上述目的，本发明采用的技术方案是：

一种基于平均曲率流的脊椎特征点自动识别方法，步骤如下：

步骤一，对人体脊椎的CT图像进行三维重建，得到所需标注特征点的椎体三维模型；

步骤二，对步骤一得到的椎体模型进行网格划分；

步骤三，手动选取模型上某点，计算在手动选取点周围极小半径为R的球形空间内每个顶点的平均曲率值；

步骤四，选取手动选取点周围平均曲率最大的n个点；

步骤五，对平均曲率最大的n个点，分别与选取点作内积，夹角最小的即为所求的点。

所述步骤四和步骤五中的n≤5。

所述步骤三的具体过程如下：

步骤3-1，按照如下步骤计算出椎体模型上每个顶点的泰森多边形区域面积A_M:

a.遍历脊椎三维图像表面的三角面片，获得其三个顶点的坐标v0、v1和v2，并得出边向量e0＝v1-v0，e1＝v2-v1，e2＝v0-v2；

b.根据向量的点积公式，将边向量(e1,e2)、(e2,e0)和(e0,e1)分别带入，可以得到三角面片各个顶角的度数，再根据cell的拓扑关系，计算脊椎图像三维表面上每个点邻域三角形对应的角度和其中#f为三角形面片的个数；

c.计算由这三个顶点组成的三角形的面积A，同理根据cell的拓扑关系可得到每个点附近三角形的面积A_M；

步骤3-2，按照如下公式计算每个顶点的法向量：

$\left\{\begin{array}{c}\left(x_i,y_i,z_i\right)=\underset{c_j\∈k_i-ring\left(v_i\right)}{\mathrm{\Σ}}area\left(c_j\right)\·\\ \end{array}\right.$

上式中，k_i-ring(v_i)表示顶点vi周围的ki环邻域范围，c_j为处于该邻域范围内的每个三角面片，area(c_j)表示c_j的面积，(x_i,y_i,z_i)为计算所得的向量和，表示三角面片c_j的单位法向量；

步骤3-3，利用以下公式遍历计算范围内每个点的平均曲率值：

$k_H\left(p_i\right)=\frac{1}{4A_M}\underset{j\&Element;N\left(i\right)}{\&Sigma;}({\mathrm{cot\&alpha;}}_{ij}+{\mathrm{cot\&beta;}}_{ij})<p_j-p_i,n>$

其中，

α_ij和β_ij是两相对角，它们所在两三角形以其对的边为公共边。

所述步骤五的具体过程如下：

步骤5-1，利用以下公式计算平均曲率值最大的n个点与手动选取点的两个向量内积：

＜p_i-O,p₀-O＞＝x_ix₀+y_iy₀+z_iz₀

其中，O为坐标系原点，p₀为手动选取点，p_i为平均曲率值最大的n个点之一；

步骤5-2，利用以下公式计算两向量的夹角θ：

$\mathrm{\&theta;}=arccos\frac{<p_i-O,p_0-O>}{\left|\right|p_i-O\left|\right|\&CenterDot;\left|\right|p_0-O\left|\right|};$

与手动选取点夹角最小的点即为经过修正后的特征点。

与现有技术相比，本发明具有以下有益效果：

通过对手动选取点进行修正获得特征点，对椎体模型的精度要求不高，容易实现，获取的特征点精度比较高；算法只对每个特征点周围的顶点的平均曲率值进行一次计算，算法复杂度与特征点个数线性相关，效率较高；使用平均曲率值计算得到的特征点，经过对手动选取点的修正，更加趋近于真实特征点；该改进算法能很好地避免由于人为因素产生的偏差，大大提高了特征点标注的准确性，对脊柱高发疾病的诊断和治疗具有一定的参考价值。

附图说明

图1为本发明方法的流程图。

图2(a)为三维重建得到人体脊椎模型正面图，(b)为三维重建得到人体脊椎模型背面图。

图3(a)为模型区域生长分割原始图像，(b)为模型区域生长分割结果图。

图4(a)为腰段L1-L3节模型建立效果图，(b)为腰段L1节腰段L1节。

图5(a)为腰椎L1节模型俯视图，(b)为腰椎L1节模型主视图。

图6为将椎体模型导入到ScanFE中的显示。

图7(a)为腰椎L1-L3节体网格模型实体单元表示，(b)为腰椎L1-L3节体网格模型节点表示。

图8为脊椎腰骶段的右侧面观特征点标记图。

图9为脊椎腰骶段的上面观特征点标记图。

图10(a)为采用高斯曲率选择特征点结果的右视图，(b)为采用高斯曲率选择特征点结果的左视图，(c)为采用高斯曲率选择特征点结果的下视图。

图11(a)为特征点1采用本发明方法结果的俯视图，(b)为特征点1采用本发明方法结果的仰视图，(c)为特征点1采用本发明方法结果的左视图，(d)为特征点1采用本发明方法结果的右视图。

图12(a)为特征点6采用本发明方法结果的俯视图，(b)为特征点6采用本发明方法结果的仰视图，(c)为特征点6采用本发明方法结果的左视图，(d)为特征点6采用本发明方法结果的右视图。

图13(a)为特征点33采用本发明方法结果的俯视图，(b)为特征点33采用本发明方法结果的仰视图，(c)为特征点33采用本发明方法结果的左视图，(d)为特征点33采用本发明方法结果的右视图。

具体实施方式

下面结合实施例对本发明作进一步详细说明。

参见图1，一种基于平均曲率流的脊椎特征点自动识别方法，包括以下步骤：

步骤一，对人体脊椎的CT图像进行三维重建，得到所需标注特征点的椎体三维模型。

参见图2、图3、图4，对脊椎CT数据进行增强和阈值提取，分割出目标椎体的范围。对椎体的分割可以使用simpleware软件的ScanIP模块进行，如图5、图6所示。在模型构建过程中，可以采用图像增强，图像阈值分割和区域生长算法，去噪点等步骤保证模型的准确。

步骤二，对步骤一得到的椎体模型进行网格划分。

设置对椎体进行网格划分的参数，利用simpleware软件的ScanIP模块生成模型的三角面片格式，如图7所示。利用该方式生成的三角面片能够自动优化局部复杂结构的网格划分，在保证整体网格质量的同时，降低网格数量，简化操作流程。

步骤三，在针对算法开发的特征点标注软件中，使用鼠标选取特征点，脊椎采样特征点的定义如图8、图9所示。计算在手动选取点周围极小半径为R的球形空间内每个顶点的平均曲率值k_H。

步骤3-1，按照如下步骤计算出椎体模型上每个顶点的泰森多边形区域面积A_M:

遍历脊椎三维图像表面的三角面片，获得其三个顶点的坐标v0、v1和v2，并得出边向量e0＝v1-v0，e1＝v2-v1，e2＝v0-v2。

根据向量的点积公式，将边向量(e1,e2)、(e2,e0)和(e0,e1)分别带入，可以得到三角面片各个顶角的度数θ_j，再根据cell的拓扑关系，计算脊椎图像三维表面上每个点邻域三角形对应的角度和

计算由这三个顶点组成的三角形的面积A，同理根据cell的拓扑关系可得到每个点附近三角形的面积A_M，计算公式如下：

$A_M=\underset{j\&Element;N\left(i\right)}{\&Sigma;}\frac{1}{2}{\&lsqb;\left|\right|p_j-p_i\left|{|}^2\right||p_{j+1}-p_i|{|}^2-{(p_j-p_i,p_{j+1}-p_i)}^2\&rsqb;}^{\frac{1}{2}}$

步骤3-2，计算每个顶点的法向量计算公式如下：

上式中，k_i-ring(v_i)表示顶点v_i周围的k_i环邻域范围，c_j为处于该邻域范围内的每个三角面片，area(c_j)表示c_j的面积，(x_i,y_i,z_i)为计算出来的向量和，表示三角面片cj的单位法向量，n_vi为顶点v_i的法向量。

步骤3-3，利用以下公式遍历计算曲面上每个点的平均曲率值k_H：

$k_H\left(p_i\right)=\frac{1}{4A_M}\underset{j\&Element;N\left(i\right)}{\&Sigma;}({\mathrm{cot\&alpha;}}_{ij}+{\mathrm{cot\&beta;}}_{ij})<p_j-p_i,n>$

其中，

α_ij和β_ij是两相对角，它们所在两三角形以其对的边为公共边。

步骤四，在手动选取点周围半径为R的球形空间内，选取平均曲率最大的5个点作为候选特征点。

以手动选取点为中心，首先确定一个半径R，对半径R内的球形空间内所有顶点的平均曲率值，进行降序排序，取得平均曲率值最大的5个点作为候选特征点。

步骤五，对5个候选特征点，分别与手动选取点作内积，夹角最小的为所求的点。

利用以下公式分别计算平均曲率值最大的5个点与手动选取点的内积。

＜p_i-O,p₀-O＞＝x_ix₀+y_iy₀+z_iz₀

其中，O为坐标系原点，p₀为手动选取点。

利用以下公式计算以坐标系原点出发的两个向量的夹角θ_i：

${\mathrm{\&theta;}}_i=arccos\frac{<p_i-O,p_0-O>}{\left|\right|p_i-O\left|\right|\&CenterDot;\left|\right|p_0-O\left|\right|}$

由坐标原点到候选点与手动选取点之间的夹角最小的候选点，即为经过算法最终修正的特征点。

对比实验：

为了突出本发明算法的特点，将采用高斯曲率作为选择特征点参考值的标注算法(如图10所示)与本发明方法(如图11、12、13所示)在同一实验模型上对特征点的对比标注情况，如下表1。

表1使用两种不同的算法在同一实验模型上对特征点的标注情况。

从表1中可以看出，对于同一个特征点，平均曲率值均大于高斯曲率值。根据曲率的几何意义得知，曲率越大，模型的几何特征愈发明显，这说明由平均曲率值算法得到的特征点包含更多几何信息，也就是说由平均曲率值算法得到的特征点更加准确，精度更高。

平均曲率流较大的候选顶点与手动选取点作向量内积，获得两顶点向量之间的夹角角度，以夹角角度最小的候选顶点来替代手动选取的特征顶点。使用平均曲率值计算得到的特征点，经过对手动选取点的修正，更加趋近于真实特征点；该改进算法能很好地避免由于人为因素产生的偏差，大大提高了特征点标注的准确性，对脊柱高发疾病的诊断和治疗具有一定的参考价值。

Claims (4)

1.一种基于平均曲率流的脊椎特征点自动识别方法，其特征在于，步骤如下：

步骤一，对人体脊椎的CT图像进行三维重建，得到所需标注特征点的椎体三维模型；

步骤二，对步骤一得到的椎体模型进行网格划分；

步骤三，手动选取模型上某点，计算在手动选取点周围极小半径为R的球形空间内每个顶点的平均曲率值；

步骤四，选取手动选取点周围平均曲率最大的n个点；

步骤五，对平均曲率最大的n个点，分别与选取点作内积，夹角最小的即为所求的点。

2.根据权利要求1所述的一种基于平均曲率流的脊椎特征点自动识别方法，其特征在于，所述步骤四和步骤五中的n≤5。

3.根据权利要求1所述的一种基于平均曲率流的脊椎特征点自动识别方法，其特征在于，所述步骤三的具体过程如下：

步骤3-1，按照如下步骤计算出椎体模型上每个顶点的泰森多边形区域面积A_M:

a.遍历脊椎三维图像表面的三角面片，获得其三个顶点的坐标v0、v1和v2，并得出边向量e0＝v1-v0，e1＝v2-v1，e2＝v0-v2；

b.根据向量的点积公式，将边向量(e1,e2)、(e2,e0)和(e0,e1)分别带入，可以得到三角面片各个顶角的度数，再根据cell的拓扑关系，计算脊椎图像三维表面上每个点邻域三角形对应的角度和其中#f为三角形面片的个数；

c.计算由这三个顶点组成的三角形的面积A，同理根据cell的拓扑关系可得到每个点附近三角形的面积A_M；

步骤3-2，按照如下公式计算每个顶点的法向量：

上式中，k_i-ring(v_i)表示顶点vi周围的ki环邻域范围，c_j为处于该邻域范围内的每个三角面片，area(c_j)表示c_j的面积，(x_i,y_i,z_i)为计算所得的向量和，表示三角面片c_j的单位法向量；

步骤3-3，利用以下公式遍历计算范围内每个点的平均曲率值：

$k_H\left(p_i\right)=\frac{1}{4A_M}\underset{j\&Element;N\left(i\right)}{\&Sigma;}({\mathrm{cot\&alpha;}}_{ij}+{\mathrm{cot\&beta;}}_{ij})<p_j-p_i,n>$

其中，

α_ij和β_ij是两相对角，它们所在两三角形以其对的边为公共边。

4.根据权利要求1或2所述的一种基于平均曲率流的脊椎特征点自动识别方法，其特征在于，所述步骤五的具体过程如下：

步骤5-1，利用以下公式计算平均曲率值最大的n个点与手动选取点的两个向量内积：

＜p_i-O,p₀-O＞＝x_ix₀+y_iy₀+z_iz₀

其中，O为坐标系原点，p₀为手动选取点，p_i为平均曲率值最大的n个点之一；

步骤5-2，利用以下公式计算两向量的夹角θ：

$\mathrm{\&theta;}=arccos\frac{<p_i-O,p_0-O>}{\left|\right|p_i-O\left|\right|\&CenterDot;\left|\right|p_0-O\left|\right|};$

与手动选取点夹角最小的点即为经过修正后的特征点。