CN105809730A - 一种长骨骨折断面点云提取方法 - Google Patents

Abstract

一种长骨骨折断面点云提取方法，涉及三维医学图像处理领域，具体涉及一种骨折断面点云提取方法。为了解决现有方法中存在的断面区域出现较大起伏时区域生长将会停止导致断面点云数据丢失比较严重的问题，本发明首先对长骨骨折CT图像进行三维重建和面绘制；依据断骨模型表面各个三角形顶点法向量与骨干轴线的夹角的差异，滤除差异较大的点；然后依据是三角面片各个顶点构成三角形网格的联系，采用点云区域生长的方法进行细分割，最终得到断面点云列表。本发明适用于长骨骨折断面点云提取。

Description

一种长骨骨折断面点云提取方法

技术领域

本发明涉及三维医学图像处理领域，具体涉及一种骨折断面点云提取方法。

背景技术

长骨多位于人体四肢上，形状呈长管状，骨干受外力撞击时容易发生骨折。骨折传统治疗手段是开刀手术并植入钢板，但其创口大,骨膜剥离范围广,出血多,愈合速度慢。此外,术后仍需辅助以石膏固定,容易发生关节僵硬。

随着科学技术的发展进步,骨折微创手术应运而生。骨折微创手术相对于传统骨折手术具有创口小、疼痛轻、恢复快、住院时间短、出血少等优点。骨折微创手术通常需要图像引导设备，即C型臂。但是由于C型臂产生的是二维透视图像，当需要准确确定骨折的位置、大小和几何形状等，医生仅通过观察图像是很难做出准确判断，手术操作严重依赖医生的经验，并且X光对医生和患者都具有潜在的辐射危险。

随着计算机图形图像处理技术的发展，三维图像导航成为了可能。断骨经过三维扫描建模，通过图像配准技术实现断面的精确对齐。其中，骨折断面上的点集构建是实现精确配准的关键，为此，世界上很多学者对这一问题提出了解决方案。相关文献中具有代表性的是论文《股骨头坏死与骨折计算机辅助手术》提出了一种股骨断面点集提取方法。该方法使用了基于法向量方法的区域生长。区域生长的条件是各个三角面片顶点单位法向量间的角度差是否小于所设阈值角。该方法处理效果受设定角度阈值影响较大，适用于比较平滑的断面。相反，如果断面区域出现较大起伏，区域生长将会停止，点云数据丢失比较严重，无法完整的提取断面点集。

发明内容

本发明是为了解决现有方法中存在的断面区域出现较大起伏时区域生长将会停止导致断面点云数据丢失比较严重的问题。

一种长骨骨折断面点云提取方法，包括以下步骤：

步骤1、对长骨骨折部位进行CT扫描，读取CT图像，使用MarchingCubes算法进行等值面提取，使用Delaunay算法进行三维重建和面绘制，采用Delaunay算法所绘制的长骨表面模型是由很多三角形构成；每个三角形的顶点会被邻近其他三角形共用；三维重建和面绘制后得到长骨模型及两个列表：顶点列表v和三角形列表f，列表v中每行存储了三角形顶点的坐标，列表f中每行存储一个三角形各个顶点在列表v中的行号；列表v和列表f的关系如图1所示；

步骤2、如图2所示，在长骨模型中部骨干的表面选取一条直线AB，构成骨干方向向量使得直线AB与长骨的中轴线几乎平行；

步骤3、求取的三角形顶点的法向量：长骨模型表面由许多三角形拼接而成，每个三角形的顶点会被邻近其他三角形共用，因此三角形某一顶点的法向量指共享该顶点的所有三角形表面法向量的平均值；

步骤4、粗分割：

由于断面与侧壁的三角形各顶点法向量与骨干方向向量的夹角差异较大，可以求出三角形每个顶点法向量与向量夹角θ；根据经验设定阈值角，大于阈值角的顶点被认为侧壁上的点，小于阈值角的顶点被认为断面的点；将列表v中小于阈值角的点的行号存入列表m中；以列表m为依据，从列表f中提取包含列表m中的点的所有三角形，重新生成三角形列表newf；

步骤5、细分割：

由于长骨模型中部骨干侧壁的局部区域出现起伏，使得所在区域的点与骨干方向向量的夹角小于所设的阈值，出现了少量剩余的侧壁点；为滤除这些侧壁点，从断面上进行细分割；

在长骨模型断面上选取若干个种子点，在列表newf中搜索与种子点邻近的三角形顶点进行种子点生长，最终得到断面点云列表pointcloud。

本发明具有以下效果：

本发明采用粗分割和细分割，能更精确提取长骨骨折断面区域点云，相比《股骨头坏死与骨折计算机辅助手术》所提出的方法，本发明对断面平滑程度没有限制，即使断面区域出现较大起伏时，骨折断面区域生长也不受影响，能够较为完整地提取断面点云用于配准，提取骨折断面点云的完成程度可达95％以上。本发明得到的点云数据可以辅助医生观察骨折断面区域进行手术操作，实现了股骨骨折手术三维图像导航，从而提高医疗诊断和治疗的准确性和可靠性。

附图说明

图1为列表v和和列表f的索引关系；

图2为骨干方向向量示意图；

图3为三角形法向量计算示意图；

图4为顶点p法向量计算示意图；

图5为列表newf生成示意图；

图6为种子点选取过程示意图；

图7为种子点生长过程示意图；

图8为新种子点生成过程示意图；

图9为输出列表pointcloud示意图；

图10为实施例中人工股骨模型的三维绘制效果图；

图11为实施例中骨干方向向量示意图；

图12为断面附近部分股骨效果图；

图13(a)为顶点法向量绘制图；

图13(b)为顶点法向量绘制图的局部放大图；

图14为粗分割结果图；

图15为选取初始种子点示意图；

图16为股骨下半部断面点云图；

图17为股骨上半部断面点云图。

具体实施方式

具体实施方式一：

一种长骨骨折断面点云提取方法，包括以下步骤：

步骤1、对长骨骨折部位进行CT扫描，读取CT图像，使用MarchingCubes算法进行等值面提取，使用Delaunay算法进行三维重建和面绘制，采用Delaunay算法所绘制的长骨表面模型是由很多三角形构成；每个三角形的顶点会被邻近其他三角形共用；三维重建和面绘制后得到长骨模型及两个列表：顶点列表v和三角形列表f，列表v中每行存储了三角形顶点的坐标，列表f中每行存储一个三角形各个顶点在列表v中的行号；列表v和列表f的关系如图1所示；

步骤2、如图2所示，在长骨模型中部骨干的表面选取一条直线AB，构成骨干方向向量使得直线AB与长骨的中轴线几乎平行；

步骤3、求取的三角形顶点的法向量：长骨模型表面由许多三角形拼接而成，每个三角形的顶点会被邻近其他三角形共用，因此三角形某一顶点的法向量指共享该顶点的所有三角形表面法向量的平均值；

步骤4、粗分割：

由于断面与侧壁的三角形各顶点法向量与骨干方向向量的夹角差异较大，可以求出三角形每个顶点法向量与向量夹角θ；根据经验设定阈值角，大于阈值角的顶点被认为侧壁上的点，小于阈值角的顶点被认为断面的点；将列表v中小于阈值角的点的行号存入列表m中；以列表m为依据，从列表f中提取包含列表m中的点的所有三角形，重新生成三角形列表newf；

步骤5、细分割：

由于长骨模型中部骨干侧壁的局部区域出现起伏，使得所在区域的点与骨干方向向量的夹角小于所设的阈值，出现了少量剩余的侧壁点；为滤除这些侧壁点，从断面上进行细分割；

在长骨模型断面上选取若干个种子点，在列表newf中搜索与种子点邻近的三角形顶点进行种子点生长，最终得到断面点云列表pointcloud。

具体实施方式二：

本实施方式步骤3所述求取的三角形顶点法向的具体步骤如下：

步骤3.1、三角形法向量的计算：

如果已知某个三角形的三个顶点，那么可使用叉积计算该三角形的法向量；设三角形的三个顶点分别为：P₁(x₁，y₁，z₁)，P₂(x₂，y₂，z₂)，P₃(x₃，y₃，z₃)；只需构建两个向量和则：和其叉积就是所求三角形法向量，如图3所示；

$\stackrel{\→}{a}\×\stackrel{\→}{b}=\left|\begin{array}{ccc}i& j& k\\ a_x& a_y& a_z\\ b_x& b_y& b_z\end{array}\right|---\left(1\right)$

式中：i、j和k分别为单位向量，

a_x＝x₂-x₁，a_y＝y₂-y₁，a_z＝z₂-z₁(2)

b_x＝x₃-x₁，b_y＝y₃-y₁，b_z＝z₃-z₁(3)

则由公式(1)、(2)和(3)可得：

$\stackrel{\→}{n}=n_x\stackrel{\→}{i}+n_y\stackrel{\→}{j}+n_z\stackrel{\→}{k}---\left(4\right)$

式中：

n_x＝a_yb_z-a_zb_y(5)

n_y＝a_zb_x-a_xb_z(6)

n_z＝a_xb_y-a_yb_x(7)

步骤3.2、顶点法向量计算：

如图4所示，p为要计算法向量的顶点，其周围有r个三角形，设第ti个三角形的单位法向量为ti＝1，2，…，r；定义顶点p的法向量为其周围所有三角形的平均值，并以表示顶点p的法向量，公式如(8)(9)所示：

${\stackrel{\→}{v}}_p=\frac{1}{r}\underset{ti=1}{\overset{r}{\Σ}}{\stackrel{\→}{n}}_{ti}---\left(8\right)$

${\stackrel{\→}{v}}_p=\frac{v_p}{\left|v_p\right|}=\frac{v_{px}\stackrel{\→}{i}+v_{py}\stackrel{\→}{j}+v_{pz}\stackrel{\→}{k}}{\sqrt{{\left(v_{px}\right)}^2+{\left(v_{py}\right)}^2+{\left(v_{pz}\right)}^2}}---\left(9\right)$

其中v_px、v_py、v_pz分别表示法向量三个坐标轴上的分量。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：

本实施方式步骤4所述求出顶点法向量与方向向量夹角θ的余弦公式如下：

$cos\mathrm{\θ}=\frac{\stackrel{\→}{v_p}\·\stackrel{\→}{AB}}{\left|\stackrel{\→}{AB}\right|\·\sqrt{{\left(v_{px}\right)}^2+{\left(v_{py}\right)}^2+{\left(v_{pz}\right)}^2}}---\left(10\right).$

其它步骤及参数与具体实施方式二相同。

具体实施方式四：

本实施方式步骤4所述重新生成三角形列表newf的具体步骤如下：

将小于阈值角的所有顶点认为断面的点；将其在列表v中的行号存入列表m中，记为m(l)，l＝1，2，…；在列表f中遍历搜索含有m(l)的行，即列表f中每一行存储一个三角形的三个顶点，三个顶点中任意一个为m(l)时，则将该行存储至列表newf；此外，在列表newf的第4列设置一个搜索标志位，初始化为0，用于在步骤5细分割时表示此行是否已经被搜索。具体构建列表newf方法如图5所示。

其它步骤及参数与具体实施方式二或三相同。

具体实施方式五：

本实施方式细分割原理：三角形网格中点与点之间具有连通性，即网格中一个点往往被邻近多个三角形所共用，因此通过列表v中一个点的索引值可以在列表newf找到包含该索引值的多行数据(即包含该点的多个三角形)；由于断面点云之间具有联通性，而侧壁上噪声点和断面的点不具有连通性；所以通过在断面上选取种子点，借助连通性不断进行区域生长，便可将侧壁噪声点滤除；

步骤5.1、初始种子点选取：

在长骨模型断面上选取若干个种子点；将这些种子点在列表v中的行号存放至列表D；如图6所示，假设在断面上选取了点(63.4,184,1)，将其在列表v中的行号11存储至列表D，将该点作为在区域生长的一个种子点；

步骤5.2、种子点生长：

取列表D中某一个种子点，假设取第1个种子点D[1]；在列表newf的前3列中搜索含有D[1]的行，将列表newf中对应的整行数据存储到列表T，并将列表newf中该行第4列标志为1，下次不再搜索；如图7所示；

步骤5.3、重复步骤5.2操作，将列表D中所有种子点生长一次；

步骤5.4、新种子点生成：

先将列表D追加至列表D_s中，然后将列表T进行单一化，即同一个数值只保留一个，得到列表T'；最后将列表T'中含有列表D的元素剔除，存储到列表S中，如图8所示；将列表S中的点集赋值给列表D，用于下次搜索；

步骤5.5、循环执行步骤5.2到步骤5.4，直到没有新的种子点生成；输出列表D_s；

步骤5.6、从列表v中提取以列表D_s中元素为行号的点的坐标，并存储至断面点云列表pointcloud，如图9所示；

其中，各个列表功能定义如下，

列表newf：将列表f中的数据经过粗分割后重新生成的表格，其中第4列标志为0，表示此行未被搜索；当标志位变为1时，表示此行已被搜索，以避免重复搜索；

列表D：每次搜索前存储种子点；

列表T：在列表newf中搜索含有列表D中元素的行，将整行数据存储到列表T；

列表T'：列表T单一化，即同一个数值只保留一个，然后存储至列表T'；

列表S：将列表T'中含有列表D的元素剔除，剩余的元素存储至列表S；

列表D_s：存储所有已搜索到的点；具体为每次区域生长时，将搜索到的行号追加至列表D_s中。

其它步骤及参数与具体实施方式一至四之一相同。

实施例

本发明针对长骨断面点云进行提取，以人工股骨模型为例。

将股骨断面点云提取作为目标，通过Matlab编制软件，具体过程如下：

步骤1、对股骨骨折部位进行CT扫描，并使用函数dicomdir读入CT图像，生成三维矩阵存放至列表volsub。

三维重建和面绘制：

a二值化处理：

设定灰度阈值，本例中使用CT值2000为阈值，灰度值小于2000的像素点，标志为0；灰度值大于2000的像素点标志为1。

b重建绘制：

使用isosurface函数提取列表volsub的等值面，得到数据列表f(等值面的面的索引信息)、列表v(等值面的顶点)。三维重建和面绘制结果如图10所示。

步骤2、以下部断骨为例，使用鼠标，在模型表面获取两点A(71，202.5,404),B(72,200.3,208)。计算向量并单位化，近似为(0,0,1)；如图11所示。

步骤3、为减小计算量，提取断面附近部分股骨，如图12。使用公式(1)-(7)求出各个三角形的单位法向量。使用公式(8)和(9)求出各个顶点的法向量，如图13(a)和(b)所示。

3.3角度统计及阈值角

步骤4、使用公式(10)求出三角形每个顶点法向量与向量夹角θ。根据经验选取角度阈值，本例中阈值为70°。

法向量滤波：

筛选出小于70°的点存储到列表e中(列表e记录小于阈值角的顶点，即记录被认为是断面的点)。绘制列表e中各点，结果如图14所示。

步骤5、

选取种子点：本实施例在断面上手动选取16个种子点，存储至列表k，其分布如图15所示。

种子点区域生长：使列表newf进行区域生长，最终得到列表D_s。

输出断面点云：按照列表D_s中的行号，在v中进行索引提取点的坐标，存储至列表pointcloud。

绘制列表pointcloud中的点，如图16所示。

同样方法可得上半部分股骨断面点云，如图17所示。

Claims (5)

1.一种长骨骨折断面点云提取方法，其特征在于包括以下步骤：

步骤1、对长骨骨折部位进行CT扫描，读取CT图像，使用MarchingCubes算法进行等值面提取，使用Delaunay算法进行三维重建和面绘制；得到长骨模型及两个列表：顶点列表v和三角形列表f，列表v中每行存储了三角形顶点的坐标，列表f中每行存储一个三角形各个顶点在列表v中的行号；

步骤2、在长骨模型中部骨干的表面选取一条直线AB，构成骨干方向向量

步骤3、求取的三角形顶点的法向量：长骨模型表面由许多三角形拼接而成，每个三角形的顶点会被邻近其他三角形共用，因此三角形某一顶点的法向量指共享该顶点的所有三角形表面法向量的平均值；

步骤4、粗分割：

求出三角形每个顶点法向量与向量夹角θ；设定阈值角，大于阈值角的顶点被认为侧壁上的点，小于阈值角的顶点被认为断面的点；将列表v中小于阈值角的点的行号存入列表m中；以列表m为依据，从列表f中提取包含列表m中的点的所有三角形，重新生成三角形列表newf；

步骤5、细分割：

在长骨模型断面上选取若干个种子点，在列表newf中搜索与种子点邻近的三角形顶点进行种子点生长，最终得到断面点云列表pointcloud。

2.根据权利要求1所述的一种长骨骨折断面点云提取方法，其特征在于步骤3所述求取的三角形顶点法向的具体步骤如下：

步骤3.1、三角形法向量的计算：

设三角形的三个顶点分别为：P₁(x₁，y₁，z₁)，P₂(x₂，y₂，z₂)，P₃(x₃，y₃，z₃)；构建两个向量和则：和其叉积就是所求三角形法向量；

$\stackrel{\→}{a}\×\stackrel{\→}{b}=\left|\begin{array}{ccc}i& j& k\\ a_x& a_y& a_z\\ b_x& b_y& b_z\end{array}\right|---\left(1\right)$

式中：i、j和k分别为单位向量，

a_x＝x₂-x₁，a_y＝y₂-y₁，a_z＝z₂-z₁(2)

b_x＝x₃-x₁，b_y＝y₃-y₁，b_z＝z₃-z₁(3)

则由公式(1)、(2)和(3)可得：

$\stackrel{\→}{n}=n_x\stackrel{\→}{i}+n_y\stackrel{\→}{j}+n_z\stackrel{\→}{k}---\left(4\right)$

式中：

n_x＝a_yb_z-a_zb_y(5)

n_y＝a_zb_x-a_xb_z(6)

n_z＝a_xb_y-a_yb_x(7)

步骤3.2、顶点法向量计算：

p为要计算法向量的顶点，其周围有r个三角形，设第ti个三角形的单位法向量为ti＝1，2，…，r；定义顶点p的法向量为其周围所有三角形的平均值，并以表示顶点p的法向量，公式如(8)(9)所示：

${\stackrel{\→}{v}}_p=\frac{1}{r}\underset{ti=1}{\overset{r}{\Σ}}{\stackrel{\→}{n}}_{ti}---\left(8\right)$

${\stackrel{\→}{v}}_p=\frac{v_p}{\left|v_p\right|}=\frac{v_{px}\stackrel{\→}{i}+v_{py}\stackrel{\→}{j}+v_{pz}\stackrel{\→}{k}}{\sqrt{{\left(v_{px}\right)}^2+{\left(v_{py}\right)}^2+{\left(v_{pz}\right)}^2}}---\left(9\right)$

其中v_px、v_py、v_pz分别表示法向量三个坐标轴上的分量。

3.根据权利要求2所述的一种长骨骨折断面点云提取方法，其特征在于步骤4所述求出顶点法向量与方向向量夹角θ的余弦公式如下：

$cos\mathrm{\θ}=\frac{\stackrel{\→}{v_p}\·\stackrel{\→}{AB}}{\left|\stackrel{\→}{AB}\right|\·\sqrt{{\left(v_{px}\right)}^2+{\left(v_{py}\right)}^2+{\left(v_{pz}\right)}^2}}---\left(10\right).$

4.根据权利要求3所述的一种长骨骨折断面点云提取方法，其特征在于步骤4所述重新生成三角形列表newf的具体步骤如下：

将小于阈值角的所有顶点认为断面的点；将其在列表v中的行号存入列表m中，记为m(l)，l＝1，2，…；在列表f中遍历搜索含有m(l)的行，即列表f中每一行存储一个三角形的三个顶点，三个顶点中任意一个为m(l)时，则将该行存储至列表newf；此外，在列表newf的第4列设置一个搜索标志位，初始化为0，用于在步骤5细分割时表示此行是否已经被搜索。

5.根据权利要求4所述的一种长骨骨折断面点云提取方法，其特征在于步骤5的具体实施步骤如下：

步骤5.1、初始种子点选取：

在长骨模型断面上选取若干个种子点；将这些种子点在列表v中的行号存放至列表D；

步骤5.2、种子点生长：

取列表D中某一个种子点，假设取第1个种子点D[1]；在列表newf的前3列中搜索含有D[1]的行，将列表newf中对应的整行数据存储到列表T，并将列表newf中该行第4列标志为1，下次不再搜索；

步骤5.3、重复步骤5.2操作，将列表D中所有种子点生长一次；

步骤5.4、新种子点生成：

先将列表D追加至列表D_s中，然后将列表T进行单一化，即同一个数值只保留一个，得到列表T'；最后将列表T'中含有列表D的元素剔除，存储到列表S中；将列表S中的点集赋值给列表D，用于下次搜索；

步骤5.5、循环执行步骤5.2到步骤5.4，直到没有新的种子点生成；输出列表D_s；

步骤5.6、从列表v中提取以列表D_s中元素为行号的点的坐标，并存储至断面点云列表pointcloud；

其中，各个列表功能定义如下，

列表newf：将列表f中的数据经过粗分割后重新生成的表格，其中第4列标志为0，表示此行未被搜索；当标志位变为1时，表示此行已被搜索，以避免重复搜索；

列表D：每次搜索前存储种子点；

列表T：在列表newf中搜索含有列表D中元素的行，将整行数据存储到列表T；

列表T'：列表T单一化，即同一个数值只保留一个，然后存储至列表T'；

列表S：将列表T'中含有列表D的元素剔除，剩余的元素存储至列表S；

列表D_s：存储所有已搜索到的点；具体为每次区域生长时，将搜索到的行号追加至列表D_s中。