CN104952109A - 一种基于Harris的脊椎特征点自动识别方法 - Google Patents

Abstract

本发明公开了一种基于Harris的脊椎特征点自动识别方法，通过CT图片进行脊椎三维模型建模和切割；获得要进行特征点识别的特定椎体；对椎体进行网格划分，使得椎体由细小的三角面片构成；对椎体模型的每个顶点，自适应计算每个顶点的阈值环数；计算模型顶点法向量；将顶点及其阈值环数范围内其他顶点的法向量利用主成分分析法获得第一第二特征向量；对第一第二特征向量利用最小二乘法进行函数拟合；对拟合后的函数利用Harris算法计算顶点变化值；对所有顶点的变化值进行降序排序，选出候选特征点；对候选特征点根据一定抑制区间进行非极大值抑制，最后获得该椎体的特征点。可用于外科医学领域中与脊椎相关的研究。

Description

一种基于Harris的脊椎特征点自动识别方法

技术领域

本发明属于医学图像处理技术领域，涉及脊椎部位的特征点识别和标定，特别是一种基于Harris的脊椎特征点自动识别方法。

背景技术

特征点检测是图像处理中热门的研究领域。特征点是指在一定范围内区别于该范围内其他顶点特征的顶点。在三维物体中，特征点可以定义为边，脊，谷，尖点等。特征点识别是指检测并且标记出给定对象上符合特征点定义的特殊点集。将特征点识别应用到医学图像能够为医学诊断提供参考数据，避免人为因素导致的失误。目前在颅骨复原等医学领域已经用到了特征点识别技术。

脊椎是人体的重要支撑骨骼，随着人们生活和工作行为方式的改变，脊椎疾病发病率越来越高，脊椎病被列为“世界十大被忽视健康问题”之一，目前仍然缺乏对脊椎特征点自动识别的研究。特征点识别可以应用到脊椎的诊疗中的病灶定位，椎体测量，虚拟手术等过程中。起到降低手术过程中人为误操作风险，同时提高手术效率，为医生提供准确患者体征数据的作用。

然而由于脊椎结构复杂，目前脊椎特征点识别的主要通过手动和半自动手段进行。手动方法操作步骤主要是依靠操作者的经验，通过特定仪器，直接在给定的椎体上面标记出特征点的位置。半自动主要是在手动操作的基础上，计算机通过高斯曲率等方法对人为手动标记的位置进行微调，从而降低手动特征点标定的误差。这两种方法对操作人员和仪器的要求高，且效率低。因此，研究脊椎特征点自动标定算法尤为必要。

发明内容

针对上述现有技术中存在的缺陷或不足，本发明的目的在于，提供一种基于Harris的脊椎特征点自动识别方法，该方法能够准确、快速、自动、有效地完成脊椎特征点的识别过程。

为了实现上述任务，本发明采用以下技术方案：

一种基于Harris的脊椎特征点自动识别方法，包括以下步骤：

步骤一，对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，对该三维模型进行切割，获得要提取特征点的椎体；

步骤二，对要提取特征点的锥体进行网格划分；

步骤三，获得每个顶点v_i周围阈值环数范围内包括该顶点在内的所有顶点，记为点集set(v_i)；

步骤四，确定网格划分后模型顶点v_i法向量；

步骤五，对点集set(v_i)进行所有顶点法向量主成份分析,得到分析后点集set”(v_i)；

步骤六，对分析后的点集set”(v_i)的第一第二特征向量，利用最小二乘法进行函数拟合；

步骤七，通过Harris算出拟合后的函数顶点变化值，降序排序，选出排序后前10％的顶点作为候选特征点。

步骤八，对候选特征点利用非极大值抑制，得到最终特征点。

进一步地，所述的步骤一的具体过程如下：

脊椎三维模型构建和分割通过商业simpleware软件的ScanIP模块进行，在构建过程中，主要利用到了前期图像增强，中期图像阈值分割和区域生长算法，后期去噪点三个步骤保证了模型的准确。

进一步地，所述的步骤二的具体过程如下：

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

进一步地，所述的步骤三的具体过程如下：

步骤S30，对于网格中的每一个顶点，按照下面的公式计算其自适应邻域面积：

上式中，n是该模型顶点总数，a-ring(v_i)是指网格化后模型中顶点v_i周围a环范围，c_j为顶点v_i周围a环范围内的一个三角面片，area(c_j)为三角面片c_j的面积，a的取值为[1,4]；

步骤S31，按照下述循环算法确定顶点的阈值环数k_i：

a.令k_i＝1；

b.计算顶点v_i周围k_i环内所有三角片面的面积之和d；

c.如果d＜自适应邻域面积，则k_i＝k_i+1，并返回步骤b，否则得到k_i。

步骤S32，把处于顶点v_i阈值环数范围内的所有顶点记为set(v_i)。

进一步地，所述的步骤四中，顶点v_i法向量计算公式如下：

$\stackrel{\→}{n_{v_i}}=\frac{{\mathrm{\Σ}}_{j=1}^{N_i}(\stackrel{\→}{v_J}-\stackrel{\→}{v_m})}{N_i},v_j\∈\mathrm{set}\left(v_i\right).$

其中v_i是待计算法向量的顶点，是顶点v_i的法向量，N_i是v_i邻域范围内的顶点个数，是模型的中心点。

进一步地，所述的步骤五的具体过程如下：

a)算出点集set(v_i)的所有法向量中心点

b)把点集内所有顶点平移到当与三维坐标系原点重合，获得新点集set’(v_i)；

c)对新点集set’(v_i)进行主成份分析，获得三个特征值和对应的特征方程；

d)把最小特征值对应的特征方程旋转到与三维坐标中z轴(0，0，1)重合，set’(v_i)的点进行同样变换，得到set”(v_i)。

进一步地，所述的步骤六的具体过程如下：

对变换后的set”(v_i)，取其x轴和y轴坐标，利用最小二乘法拟合出函数：

$z=f(x,y)=\frac{p_1}{2}{x}^2+p_2\mathrm{xy}+\frac{p_3}{2}{y}^2+p_{4}x+p_{5}y+p_6.$

进一步地，所述的步骤七的具体过程如下：

a.根据拟合后的二维函数，获得Harris的对称矩阵

$E=\left[\begin{array}{cc}\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{{2\mathrm{\σ}}^2}}.f_x{(x,y)}^2\mathrm{dxdy}& \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{\∫e}^{\frac{-(x^2+y^2)}{{2\mathrm{\σ}}^2}}.f_x(x,y).f_y(x,y)\mathrm{dxdy}\\ \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}{\∫e}^{\frac{-(x^2+y^2)}{{2\mathrm{\σ}}^2}}.f_x(x,y).f_y(x,y)\mathrm{dxdy}& \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{{2\mathrm{\σ}}^2}}.f_y{(x,y)}^2\mathrm{dxdy}\end{array}\right].$

其中，f_x为拟合函数在x方向上的导，f_y为拟合函数在y方向上的偏导。(x_i，y_i)为点集set”(v_i)的坐标。σ为高斯窗口大小，设为自适应邻域面积的开方。

b.计算set”(v_i)中，顶点v_i的变化值：

h(x，y)＝det(E)-k*tr(E)²

其中，设λ₁，λ₂是矩阵E的特征值，则det(E)＝λ₁*λ₂，tr(E)＝λ₁+λ₂，k＝0.04。

c.对所有顶点的变化值h(x,y)进行降序排序，选出前10％所对应的作为候选特征点。

进一步地，步骤七中，非极大值抑制的具体过程为：

对于两个顶点距离小于抑制范围r的情况，保留特征值D(v_i)最大的点，舍弃特征值较小的点；特征值D(v_i)最大的点即为特征点。

本发明具有以下技术特点：

1、特征点提取过程全自动，无需手动提取，算法效率较高。通过输入建模好的椎体模型，能够自动标记出椎体模型中特征点的位置；此外，每个顶点只需计算一次，算法复杂度与模型顶点数量成线性关系，算法效率高。

2、针对脊椎的结构特点进行专门优化。已有的三维模型特征点自动标定方法中，顶点向量是定义为顶点与模型中心原点的向量差。当模型中心原点位于模型外部时，这种方式定义出来的顶点向量无法准确地表达顶点的变化程度，从而无法筛选出模型中的特征点。本发明将这种顶点法向量计算方法应用于Harris算法上，有效解决了中心点不在模型内部的问题。

3、实验表明，本发明的方法对脊椎模型，尤其是腰椎模型上进行特征点识别上能够达到较好的效果。

附图说明

图1为本发明的基于Harris的脊椎特征点自动识别方法流程图；

图2为ScanIP模块计算时的示意图；

图3为构建出来的腰椎L3三维模型；

图4为腰椎L3的网格化后的三维模型；

图5为a环邻域范围的说明图；

图6(a)为定义特征点斜俯视图；

图6(b)为定义特征点斜仰视图；

图6(c)为定义特征点俯视图；

图6(d)为定义特征点仰视图；

图6(e)为定义特征点左视图；

图6(f)为定义特征点右视图；

图6(g)为定义特征点后视图；

图6(h)为定义特征点前视图；

图7(a)为本方法结果俯视图；

图7(b)为本方法结果仰视图；

图7(c)为本方法结果左视图；

图7(d)为本方法结果右视图；

图7(e)为本方法结果后视图；

图7(f)为本方法结果前视图；

图8(a)为方法1结果俯视图；

图8(b)为方法1结果仰视图；

图8(c)为方法1结果左视图；

图8(d)为方法1结果右视图；

图8(e)为方法1结果后视图；

图8(f)为方法1结果前视图；

图9(a)为方法2结果俯视图；

图9(b)为方法2结果仰视图；

图9(c)为方法2结果左视图；

图9(d)为方法2结果右视图；

图9(e)为方法2结果后视图；

图9(f)为方法2结果前视图；

以下结合附图和实施例对本发明进行进一步详细说明。

具体实施方式

参见图1，图1分别展示应用本发明的方法在人体腰椎L3上特征点识别的过程和结果。该模型在论文(Wang Y,Yu B,Dai J.Three-DimensionalFinite Element Force Analysis and Simulation of Human Spinal LumbarSegment[J].Journal of Convergence Information Technology,2013,8(4).)中得到验证。

步骤一，对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，对该三维模型进行切割，获得要提取特征点的椎体；

如图2所示，对CT数据进行增强和阈值提取后，分割出目标椎体的范围，得到如图3所示的三维模型。具体构建和分割过程，可通过simpleware软件的ScanIP模块进行进行，在构建过程中，利用前期图像增强，中期图像阈值分割和区域生长算法，后期去噪点三个步骤保证模型的准确。

步骤二，对要提取特征点的锥体进行网格划分；

将单个锥体设置网格相应的参数，利用simpleware软件的ScanIP模块生成模型的三角面片格式。利用该方式生成的三角面片能够自动优化局部复杂结构的网格划分，在保证整体网格质量的同时，降低网格数量，简化操作流程。本例中，利用ScanIP模块将腰椎L3表面数据转化为网格后输出Ascii码格式的STL文件。划分后的网格模型如图4所示。

步骤三，确定划分后网格中每个顶点阈值环数及其点集；

步骤S30，对于网格中的每一个顶点，设定一个自适应邻域面积大小，自适应邻域面积是为了动态调整每个顶点邻域范围的大小；网格化后的模型是由很多三角形面片组成，每个三角形面片有三个顶点，所有的这些顶点都成为模型的顶点。从这些点中挑选出一些有特别特征的点，如突出点，可作为特征点；自适应邻域面积由以下公式计算：

上式中，n是该模型顶点总数，a-ring(v_i)是指网格化后模型中顶点v_i周围a环范围(a-ring)，图5展示了a＝1,2,3时顶点周围a环三角面片的情况。c_j为顶点v_i周围a环范围内的一个三角面片。area(c_j)为三角面片c_j的面积。在实际操作中，a的取值根据模型网格精细程度进行调整，一般取值：[1,4]；本实例中，a取值2，自适应邻域面积为37.3702。自适应邻域面积是为了确定阈值环数k_i。每个顶点参考自适应邻域面积来获得顶点各自的k_i值；

步骤S31，根据步骤S30得到的自适应邻域面积，计算顶点v_i的阈值环数k_i环邻域中k_i的取值，具体计算过程为：

a、令k_i＝1；

b、计算顶点v_i周围k_i环内所有三角片面的面积之和d；

c、如果d＜自适应邻域面积，则k_i＝k_i+1，并返回步骤b，否则得到k_i；

步骤S32，把处于顶点v_i阈值环数范围内的所有顶点记为set(v_i)。

步骤四，计算每个顶点v_i的法向量，计算公式如下：

$\stackrel{\→}{n_{v_i}}=\frac{{\mathrm{\Σ}}_{j=1}^{N_i}(\stackrel{\→}{v_J}-\stackrel{\→}{v_m})}{N_i},v_j\∈\mathrm{set}\left(v_i\right).$

上式中，v_i是待计算法向量的顶点，是顶点v_i的法向量，N_i是v_i邻域范围内的顶点个数，是模型的中心点。

步骤五，对每个顶点v_i的点集set(v_i)所包含所有顶点的法向量进行主成份分析，具体过程为：

a、算出点集set(v_i))的所有法向量中心点

b、把点集内所有顶点平移到当与三维坐标系原点重合，获得新点集set’(v_i)；

c、对新点集set’(v_i)进行主成份分析，获得三个特征值和对应的特征方程；

d、把最小特征值对应的特征方程旋转到与三维坐标中z轴(0，0，1)重合，set’(v_i)的点进行同样变换，得到set”(v_i)。

步骤六，对分析后的点集set”(v_i)的第一第二特征向量，利用最小二乘法进行函数拟合，具体计算过程为：

对变换后的set”(v_i)，取其x轴和y轴坐标，利用最小二乘法拟合出函数：

$z=f(x,y)=\frac{p_1}{2}{x}^2+p_2\mathrm{xy}+\frac{p_3}{2}{y}^2+p_{4}x+p_{5}y+p_6$

步骤七，通过Harris对拟合后的函数算出顶点变化值，降序排序，选出排序后前10％的顶点作为候选特征点。

将每个顶点的变化值代入快速排序算法进行降序排列，得到变化值从大到小的排列结果；由于模型顶点个数在10～4左右，顶点个数多且每个顶点的变化值混乱无序，因此快速排序相比起其他排序算法能够获得更好的时间和空间上的性能。相比起其他特征点的阈值响应方法，利用排序能够快速调整特征点拾取个数，更加适合于脊椎特征点这一新课题的研究。

步骤八，对候选特征点利用非极大值抑制，得到最终特征点。

对于两个顶点距离小于抑制范围r的情况，则保留特征值D(v_i)最大的点，舍弃特征值较小的点；特征值D(v_i)最大的点即为特征点；r根据模型大小设置，一般设置为k_i环邻域的半径；本实例中r＝10。最终筛选出特征点见图7。

对比试验：

为了突出本方法的有效性，展示了传统特征点自动识别算法(Zhou Z D,Ai Q S,Liu Q,et al.A novel feature points selection algorithm for3D triangular mesh models[C]//Signal Processing,2008.ICSP 2008.9thInternational Conference on.IEEE,2008:1023-1026,简称方法1)和(Walter，N.，Aubreton，O.，Laligant，O.:Salient pointcharacterizationfor low resolution meshes.In:ICIP,pp.1512–1515(2008),简称方法2)在本方法实验模型上的应用效果。

参考论文(Jun D，Bin Y,，Ying W.Statistic Model of the Spine inThree-Dimension Geometry[c]//Information Science andEngineering(ISISE),2010 International Symposium on.IEEE，2010:66-70.)定义本实验模型L3上34个特征点，如图6所示。

图7展示了本方法的结果，在自适应面积＝32.2，抑制距离＝10时，正确提取出26个特征点。

图8展示了方法1的结果，a-ring环数取a＝3时，抑制距离＝10时正确提取出23个特征点。

图9展示了方法2的结果，r＝2.0时，抑制距离＝10时正确提取出26个特征点。

从结果可以看出，本发明的方法在针对脊椎特殊模型时，能够达到较好地自动特征点标定效果。方法1除了特征点标定正确率低之外，特征点标定的误差也大于方法2和本方法。其原因在于方法1的顶点法向量计算方式无法适应模型中心点在模型外部的情况。本实验椎体L3的模型中心点处于模型外部，在椎孔中间。

同时，方法2除了在正确率和误差上比方法1的结果好。但是方法2的计算效率低。在计算每个独立球形里面不同三角面片占用体积时，需要根据球形球心的不同而重新计算。

Claims (9)

1.一种基于Harris脊椎特征点自动识别方法，其特征在于，包括以下步骤：

步骤一，对人体脊椎的CT图像进行三维重建，得到脊椎的三维模型，对该三维模型进行切割，获得要提取特征点的椎体；

步骤二，对要提取特征点的锥体进行网格划分；

步骤三，获得每个顶点v_i周围阈值环数范围内包括该顶点在内的所有顶点，记为点集set(v_i)；

步骤四，确定网格划分后模型顶点v_i法向量；

步骤五，对点集set(v_i)进行所有顶点法向量主成份分析,得到分析后点集set”(v_i)；

步骤六，对分析后的点集set”(v_i)的第一第二特征向量，利用最小二乘法进行函数拟合；

步骤七，通过Harris对拟合后的函数算出顶点变化值，降序排序，选出排序后前10％的顶点作为候选特征点；

步骤八，对候选特征点利用非极大值抑制，得到最终特征点。

2.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤一的具体过程如下：

脊椎三维模型构建和分割通过商业simpleware软件的ScanIP模块进行；在构建过程中，主要利用前期图像增强，中期图像阈值分割和区域生长算法，后期去噪点三个步骤，以保证模型的准确性。

3.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤二的具体过程如下：

将单个椎体设置网格相应的参数，利用simpleware软件的ScanIP模块生成模型的三角面片格式，利用生成的三角面片能够自动优化局部复杂结构的网格划分，在保证整体网格质量的同时，降低网格数量，简化操作流程。

4.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤三的具体过程如下：

步骤S30，对于网格中的每一个顶点，按照下面的公式计算其自适应邻域面积：

式中，n是该模型顶点总数，a-ring(v_i)是指网格化后模型中顶点v_i周围a环范围，c_j为顶点v_i周围a环范围内的一个三角面片，area(c_j)为三角面片c_j的面积，a的取值为[1,4]；

步骤S31，按照下述循环算法确定顶点的阈值环数k_i：

a、令k_i＝1；

b、计算顶点v_i周围k_i环内所有三角片面的面积之和d；

c、如果d＜自适应邻域面积，则k_i＝k_i+1，并返回步骤b，否则得到k_i；

步骤S32，把处于顶点v_i阈值环数范围内的所有顶点记为set(v_i)。

5.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤四的具体过程如下：

$\stackrel{\→}{n_{v_i}}=\frac{{\mathrm{\Σ}}_{j=1}^{N_i}(\stackrel{\→}{v_j}-\stackrel{\→}{v_m})}{N_i},v_j\∈\mathrm{set}\left(v_i\right)$

其中v_i是待计算法向量的顶点，是顶点v_i的法向量，N_i是v_i邻域范围内的顶点个数，是模型的中心点。

6.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤五的具体过程如下：

a、算出点集set(v_i)的所有法向量中心点

b、把点集内所有顶点平移到当与三维坐标系原点重合，获得新点集set’(v_i)；

c、对新点集set’(v_i)进行主成份分析，获得三个特征值和对应的特征方程；

d、把最小特征值对应的特征方程旋转到与三维坐标中z轴(0，0，1)重合，set’(v_i)的点进行同样变换，得到set”(vi)。

7.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤六的具体过程如下：

对变换后的set”(v_i)，取其x轴和y轴坐标，利用最小二乘法拟合出函数： $z=f(x,y)=\frac{p_1}{2}{x}^2+p_2\mathrm{xy}+\frac{p_3}{2}{y}^2+p_{4}x+p_{5}y+p_6.$

8.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，所述的步骤七的具体过程如下：

a、根据拟合后的二维函数，获得Harris的对称矩阵

$E=\left[\begin{array}{cc}\frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{2{\mathrm{\σ}}^2}}\·f_x{(x,y)}^2\mathrm{dxdy}& \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{2{\mathrm{\σ}}^2}}\·f_x(x,y)\·f_y(x,y)\mathrm{dxdy}\\ \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{2{\mathrm{\σ}}^2}}\·f_x(x,y)\·f_y(x,y)\mathrm{dxdy}& \frac{1}{\sqrt{2\mathrm{\π}}\mathrm{\σ}}\∫e^{\frac{-(x^2+y^2)}{2{\mathrm{\σ}}^2}}\·f_x{(x,y)}^2\mathrm{dxdy}\end{array}\right]$

其中，f_x为拟合函数在x方向上的导，f_y为拟合函数在y方向上的偏导。(x_i，y_i)为点集set”(v_i)的坐标。σ为高斯窗口大小，设为自适应邻域面积的开方。

b、计算set”(v_i)中，顶点v_i的变化值：

h(x，y)＝det(E)-k*tr(E)²

其中，设λ₁，λ₂是矩阵E的特征值，则det(E)＝λ₁*λ₂，tr(E)＝λ₁+λ₂，k＝0.04。

c、对所有顶点的变化值h(x,y)进行降序排序，选出前10％所对应的作为候选特征点。

9.如权利要求1所述的基于Harris脊椎特征点自动识别方法，其特征在于，步骤七中，非极大值抑制的具体过程为：

对于两个顶点距离小于抑制范围r的情况，保留特征值D(v_i)最大的点，舍弃特征值较小的点；特征值D(v_i)最大的点即为特征点。