一种面向骨骼CT序列图像的网格模型自适应重建方法与流程

本发明涉及图像处理，尤其涉及一种面向骨骼CT序列图像的网格模型自适应重建方法。

背景技术：

CT扫描作为一种X射线计算机体层成像技术，可以无损且准确的获取活体内部组织器官的影像。在骨骼临床应用中，能否准确重建骨骼的三维网格模型直接关系到后续的骨骼计算机辅助制造、骨骼的有限元分析及3D打印效果。目前将骨骼CT序列图像重建为网格模型的方法主要有：移动立方体算法(Marching Cubes)和基于轮廓集直接三角化的方法。移动立方体算法是逐个处理CT体数据场中的体素，以三角形面片来逼近体素内部的等值面，由于能够产生体素量级的网格，使得骨骼网格模型包含大量的冗余网格，严重影响骨骼生物力学中有限元体网格划分。基于轮廓集直接三角化的方法是通过直接连接相邻轮廓线，从而形成分段的三角网格，然后依次连接所有相邻的轮廓之间的三角网格，最终生成整体的三角网格模型，但是当相邻轮廓有多条轮廓线需要重建时，会出现分支和多元匹配问题，尤其在重建复杂骨骼模型时就会出现拓扑连接错误，同时重建的网格模型含有大量的狭长三角网格，并且不同尺寸的三角网格之间不能合理过渡，表面光顺效果差。

另外，上述两种方法都不能根据骨骼表面曲率特征对网格进行自适应划分，在低曲率表面区域生成了大量小尺寸三角网格，在高曲率表面区域不能进行局部加密网格，不能有效保证网格的离散精度和网格质量，严重影响后续的分析和使用。因此如何利用骨骼表面的曲率自适应重建高质量的三角网格模型，是当前急需解决的问题。

技术实现要素：

发明目的：本发明针对现有技术存在的问题，提供一种面向骨骼CT序列图像的网格模型自适应重建方法。

技术方案：本发明所述的面向骨骼CT序列图像的网格模型自适应重建方法包括：

(1)依次提取每层骨骼CT序列图像中骨骼的二维轮廓数据；

(2)将形成的二维轮廓数据集转化为带有法向量的三维点云；

(3)根据三维点云构建骨骼模型的隐式曲面；

(4)对骨骼隐式曲面进行自适应图元采样；

(5)对采集的自适应图元集进行三角化，形成骨骼网格模型。

其中，步骤(1)具体包括：

依次提取每层骨骼CT序列图像中骨骼的轮廓边缘、细化、连接、轮廓追踪，获取单像素封闭的骨骼二维轮廓，形成二维轮廓数据集。

其中，步骤(2)具体包括：

(2-1)对二维轮廓数据集中每层二维轮廓的内部图像进行边缘填充，以轮廓边缘为约束，对轮廓边缘内部连通区域进行八向连通填充为白色，轮廓外部边缘连通区域填充为黑色，形成轮廓二值图像；

(2-2)对多层轮廓二值图像的序列形成的三维体素集V进行三维高斯滤波，得到滤波后的体素集V'；

(2-3)使用空间Prewitt算子计算二维轮廓的像素点所在的体素位置处的梯度；

(2-4)通过二维轮廓像素点所在的体素位置梯度计算出二维轮廓像素点的单位法向量；

(2-5)依次提取每层图像的像素空间位置pi(x,y,z)，得到三维点云数据P。

其中，步骤(3)具体包括：

(3-1)使用kd-tree把三维点云数据P分解为m个局部子区域，在每个局部子区域上建立局部隐式曲面；

(3-2)使用加权函数将局部隐式曲面进行拼接，形成完整的隐式曲面。

其中，步骤(4)具体包括：

(4-1)对骨骼的隐式曲面进行图元采样，并计算所有采样图元的高斯曲率；

(4-2)利用高斯曲率，使采样图元在骨骼表面自适应分布，得到自适应图元。

其中，步骤(5)具体包括：

(5-1)计算每个自适应图元对应的法向量；

(5-2)根据法向量，使用Delaunay三角剖分对自适应图元集进行三角化，形成自适应分布的骨骼三角网格模型。

有益效果：本发明与现有技术相比，其显著优点是：本发明提出的面向骨骼CT序列图像的网格模型自适应重建方法，根据骨骼表面曲率大小生成不同尺寸的三角网格，在低曲率表面区域生成稀疏的大尺寸三角网格，在高曲率表面区域生成密集的小尺寸三角网格，这样能够更好的表达骨骼局部细节特征，并且不同尺寸的三角网格之间能合理过渡，表面光顺效果好，同时避免了骨骼模型中生成大量冗余的三角网格和狭长三角网格，能有效保证网格的离散精度和网格质量。

附图说明

图1是本发明的方法流程示意图；

图2是多层骨骼CT序列图像；

图3是单层骨骼CT图像；

图4是单层骨骼CT图像的二维轮廓；

图5是单层骨骼轮廓的二值图像；

图6是骨骼的三维点云模型；

图7是骨骼隐式曲面的自适应图元采样；

图8是重建的骨骼的三角网格模型。

具体实施方式

图1是骨骼CT序列图像自适应重建网格模型的方法流程图，如图所示：本实施例包括如下步骤：

S1：骨骼的多层CT序列图像，如图2所示，依次提取每个单层CT序列图像(见图3)的骨骼轮廓边缘、细化、连接、轮廓追踪，获取单像素封闭的骨骼二维轮廓，如图4所示。

S2：将形成的多层二维轮廓数据集转化为带有法向量的三维点云。结合图5和图6说明，具体步骤如下：

1)对封闭的二维轮廓内部图像进行边缘填充，以轮廓边缘为约束，对轮廓边缘内部连通区域进行八向连通填充为白色，轮廓外部边缘连通区域填充为黑色，形成轮廓二值图像，如图5所示；

2)对多层轮廓二值图像的序列形成的三维体素集V进行三维高斯滤波，得到滤波后的体素集V'，其中三维高斯滤波函数表达式如下：

式中[x,y]表示轮廓像素点所在的CT图像平面内坐标，z为切片的层数所在的坐标，σ是标准偏差，在本发明的实施例中，σ取3；

3)使用空间Prewitt算子计算轮廓像素点所在的体素位置处pi(x,y,z)的梯度▽V'(pi)；

4)通过轮廓像素点所在的体素梯度计算出轮廓像素点的单位法向量ni，其计算表达式如下：

5)依次提取每层图像的像素空间位置pi(x,y,z)，得到三维点云数据如图6所示；

S3：根据三维点云构建骨骼模型的隐式曲面，具体步骤如下：

1)对于局部隐式小曲面的建立，首先定义每个局部子区域坐标系(u,v,w)，其中，(u,v)作为局部平面坐标，w作为平面正交方向，使用正交多项式基，计算局部子区域坐标系下的正交多项式hi(u,v)：

式中，β1＝1，j＝1,2,…,k＝1,2…j-1，wj表示局部正交方向，xj表示未知量；

再根据正交多项式hi(u,v)构造出局部拟合函数Qi(x)＝Qi(u,v,w)＝w-hi(u,v)，其中(u,v,w)是x在局部子区域坐标系中的坐标。

对于局部隐式大曲面的建立，使用三维二次曲面函数进行曲面重建，此时局部拟合函数为Qi(x)＝x^TAx+b^Tx+c，其中A是3×3对称矩阵,b是空间向量,c为常系数。

2)使用加权函数将局部隐式曲面进行拼接，形成完整的隐式曲面，具体步骤如下：

每个局部子区域包含在半径为Ri的圆球内，每个圆球内，局部拟合函数Qi(x)通过最小二乘法拟合，通过加权函数将局部拟合函数Qi(x)在欧式空间Ω中进行叠加形成完整隐式曲面f(x)，表达式如下：

x∈Ω

式中，权重函数wi(x)为非均匀有理B样条函数,ci为圆球的球心。

S4：对骨骼隐式曲面进行自适应图元采样，具体步骤如下：

1)计算骨骼模型的隐式曲面f(x)上所有采样图元{qⁱ}的高斯曲率Ki；

2)利用高斯曲率，使采样图元自适应分布，如图7所示。计算采样图元之间的排斥能量Eⁱ和采样图元的移动速度vⁱ，并更新图元的位置，使采样图元根据骨骼表面高斯曲率自适应分布；其中排斥能量Eⁱ的表达式如下：

其中Ki是采样图元qⁱ的高斯曲率，E^ij是采样图元qⁱ和采样图元q^j之间的排斥能量，因为σⁱ≠σ^j，所以E^ij≠E^ji，μ是排斥参数,|r^ij|是采样图元qi和采样图元q^j之间的距离,σⁱ是采样图元间排斥半径，采样图元的移动速度vⁱ为排斥能量Eⁱ的导数，表达为：

S5：对采集的自适应图元集进行三角化，形成骨骼网格模型，具体步骤如下：

1)计算每个自适应图元qⁱ对应的法向量Ni(qⁱ)，其中Ni(qⁱ)＝▽f(qⁱ)，▽f(qⁱ)是图元qⁱ处的梯度值。

2)使用Delaunay三角剖分对图元集{qⁱ}进行三角化，形成自适应分布的三角网格模型，如图8所示。

本发明可以避免骨骼模型中生成大量冗余的三角网格和狭长三角网格，可以根据骨骼表面的曲率自适应重建高质量的三角网格模型，同时不同尺寸的三角网格之间合理过渡，网格光顺效果好。

以上所揭露的仅为本发明一种较佳实施例而已，不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。