论文部分内容阅读
摘要：在3D打印领域中，针对找到受损骨骼待修补部分的问题上，提出一种利用模板进行配准找出需要修补部分的方法。这种模板可以是相似度高的模型，也可以是病人的镜像骨。在完成方法过程中，利用到了最小盒算法和ICP配准，前者是对模型的粗配准，后者是对模型的精配准;在最后进行剔除错误部分和平滑，生成打印模型。实验表明，该方法能完整精确打印出待修补处三维模型，与缺损模型平滑连接。
关键词：3D打印;骨骼修复;三维重建;刚性配准;点云
中图分类号：TP301.6 文献标识码：A
文章编号：1009-3044（2019）31-0231-03
1概述
随着3D打印技术的飞速发展，医学领域和计算机领域的快速结合，也带来了更多优秀人才去加入研究骨骼缺损修复这一领域。多方面领域人才的结合更快的推动了骨骼修复的发展，生物医疗研究人员致力于生产优质的仿生材料去贴合人体构造，而在计算机领域则是构造更完美的模型去进行打印。因此，如何设计更适应骨骼缺陷的三维模型成了如今研究的重点。
在构建契合破损骨骼的三维模型时，关键在于我们要得到精确的形状，并使其能够与正常骨平滑连接。在本文的研究中，我们应需要找到合适的待匹配模板，而这种模板是自于模型库或病人的对侧骨。所有的后续研究都是在这样的基础上进行的。
2总体流程
该修复方法首先重建受损骨骼模型嘲由已知的CT序列得到;接着对模型预处理，选取合适的模板;再通过对应的点云模型求取变换矩阵组，最后进行模型配准和计算缺损部分，具体流程如图1所示。重建过程在win10系统环境下完成，采用
2.1骨骼模型重建
根据原始CT图像的横断面、冠状面及矢状面逐层分割，能更好地提高图像的完整性和准确性。CT数值反应的是组织的密度，我们可以参照正常人体组织CT值表，通过调窗去除大部分非骨骼组织，实现粗分割;然后再利用分割算法[3]进行精分割，最后对分割结果平滑处理。
2.1.1骨骼粗分割
通过查表可知，骨的平扫CT值在150～1000，因此我们对原始数据进行范围规定化，在此过程中，若CT值小于150，则灰度值置为0，;若CT值高于1000，则灰度值置为255;若CT值大于150而小于1000，則灰度值按照（1）得出，在（1）中g（i，j，k）为灰度值，f（i，j，k）为CT值。
g（i，j，k）=（f（i，j，k）-150/（1000-150）*255 （1）
2.1.2骨骼精分割
由于目标具有相同的灰度特征，我们使用Otsu全局阈值分割方法方法是根据图像上灰度值的分布将图像分成背景和前景两部分看待，前景就是我们要按照阈值分割出来的部分。背景和前景的分界值就是我们要求出的阈值。遍历不同的阈值，计算不同阈值下对应的背景和前景之间的类间方差，当类间方差取得极大值时，此时对应的阈值就是OTSU算法所求的阈值。
2.1.3平滑去噪
在进行骨骼精分割以后，图像中也存在少许噪点。此外，骨骼边缘比较粗糙，会直接影响后续点云模型的精度。所以我们还需要进行平滑处理。
步骤1形态学开运算去噪：可选择小半径的圆形结构元，对精分割结果进行开运算。可在保持骨骼不变的情况下去除小噪点，并断开不同骨骼间的伪连接。
步骤2高斯滤波平滑边缘：选择小半径高斯滤波器对去噪结果进行平滑，然后以阈值T=0.5分割平滑结果。
2.2求取缺损部分模型
求取缺损模型首先应该是在适配模型库中找到与缺损模型最相近的完整模型，由于我们需要得到最契合的缺损模型，因此本文提出的方法是在获取模型点云的基础上进行配准的。利用MarchingCubes算法获取体模表面三角片，并计算三角片重心作为点云集。接着进行点云降采样，这样做的目的是减少点云数量，提高配准速率。
2.2.1点云粗配准
在得到匹配模板和缺损模型的点云后，利用最小包围盒算法进行点云的粗配准。最小包围盒也叫外接最小矩阵，是一种求解离散点集最优包围空间的算法，基本思想是用体积稍大且特性简单的几何体（称为包围盒）来近似地代替复杂的几何对象。
首先，获取已知点云的三个主方向，获取质心，计算协方差，获得协方差矩阵，求取协方差矩阵的特征值和特长向量，特征向量即为主方向。接着，利用已获得的主方向和质心，将输入点云转换至原点，且主方向与坐标系方向重回，建立变换到原点的点云的包围盒。最后，给输入点云设置主方向和包围盒，通过输人点云到原点点云变换的逆变换实现。
通俗来说，最小包围盒算法是求解点集XOY平面内投影点的最小面积包围盒，即把空间化问题转化为平面问题来求解。
2.2.2点云精配准
由兰州交通大学杨军等人提出的改进的ICP算法是本文用于点云精配准的主要方法，在该方法中主要在初始值的敏感性以及整体与部分模型的配准精度上做出了改进。选择改进的ICP算法适用于精配准是由于在初始配准后还未达到想要的效果，因此选用KD-Tree搜索大幅度地提高了搜索效率，制定自适应最有阈值则是在经典的ICP算法上提高了配准精度。
假如我们手上的数据有原始点云数据P和经过平移加旋转后的点云Q，那么ICP算法就是求出最优的从P到Q的坐标变换矩阵。从（2）中我们可以看到我们从P到Q的过程，其中旋转为R.平移为T。
2.3计算缺损部分
为解决无法通过点云直接得到缺损部分的问题，我们以点云模型为基准建立初始采样网格，后利用变换矩阵对网格进行逆运算，再对原始体数据进行采样，即可得到配准后的体模型。
把适配模型网格称为模板，缺损模型网格称为实例，那么现在用模板减去实例就可以得到粗结果。当然我们应该在粗结果上进行进行错误部分的提出，其中包括删除扁平、细长、颗粒状部分，最后我们需要利用三维高斯核和球形中值滤波器对剩余部分进行平滑处理。
3实验
在本实验中选用羊的胫骨作为实验对象，数据集一共83张516KB的CT图像，CT数据如图2所示。读取整个序列，并根据dicom参数（ResaleSlope、RescaleIntercept、WindowCenter、WindowWidth、ShceLocationl进行调窗与调序，获得数值为0～255的灰度图像;接着利用阈值分割提取所有骨骼，利用连通分量筛选目标骨骼;再根据dicom参数（PixelSpacing、SliceTh-ickness），重建序列体数据，过程示例如图3所示。
在选取到适匹配模板后，那就可以生成缺损模型和匹配到的完整模型的点云，如图4所示，经过点云的粗配准即最小包围盒算法以及精配准ICP算法后的点云如图5所示，在经过这一步操作后我们可以得到变换矩阵，这一矩阵的得到是为了下一步做铺垫。
在得到缺损模型前，需要用采样网格对点云进行包围，在图6中三维坐标下竖直采样网格代表缺损模型，倾斜采样网格代表匹配到的完整模型，在已知变换矩阵的条件下对采样网格进行逆运算，从而得到配准后模型。
那么在这样的情况下，就可以进行两模型的直接相减得到粗略的缺损模型，如图7所示，在经过剔除错误部分和平滑后，就可以生成打印模型。MarchingCubes面绘制算法获取三角片，并以光顺优化获得更好的三角片，最后以sTL的格式导出，结果的缺损模型如图8所示。
关键词：3D打印;骨骼修复;三维重建;刚性配准;点云
中图分类号：TP301.6 文献标识码：A
文章编号：1009-3044（2019）31-0231-03
1概述
随着3D打印技术的飞速发展，医学领域和计算机领域的快速结合，也带来了更多优秀人才去加入研究骨骼缺损修复这一领域。多方面领域人才的结合更快的推动了骨骼修复的发展，生物医疗研究人员致力于生产优质的仿生材料去贴合人体构造，而在计算机领域则是构造更完美的模型去进行打印。因此，如何设计更适应骨骼缺陷的三维模型成了如今研究的重点。
在构建契合破损骨骼的三维模型时，关键在于我们要得到精确的形状，并使其能够与正常骨平滑连接。在本文的研究中，我们应需要找到合适的待匹配模板，而这种模板是自于模型库或病人的对侧骨。所有的后续研究都是在这样的基础上进行的。
2总体流程
该修复方法首先重建受损骨骼模型嘲由已知的CT序列得到;接着对模型预处理，选取合适的模板;再通过对应的点云模型求取变换矩阵组，最后进行模型配准和计算缺损部分，具体流程如图1所示。重建过程在win10系统环境下完成，采用
2.1骨骼模型重建
根据原始CT图像的横断面、冠状面及矢状面逐层分割，能更好地提高图像的完整性和准确性。CT数值反应的是组织的密度，我们可以参照正常人体组织CT值表，通过调窗去除大部分非骨骼组织，实现粗分割;然后再利用分割算法[3]进行精分割，最后对分割结果平滑处理。
2.1.1骨骼粗分割
通过查表可知，骨的平扫CT值在150～1000，因此我们对原始数据进行范围规定化，在此过程中，若CT值小于150，则灰度值置为0，;若CT值高于1000，则灰度值置为255;若CT值大于150而小于1000，則灰度值按照（1）得出，在（1）中g（i，j，k）为灰度值，f（i，j，k）为CT值。
g（i，j，k）=（f（i，j，k）-150/（1000-150）*255 （1）
2.1.2骨骼精分割
由于目标具有相同的灰度特征，我们使用Otsu全局阈值分割方法方法是根据图像上灰度值的分布将图像分成背景和前景两部分看待，前景就是我们要按照阈值分割出来的部分。背景和前景的分界值就是我们要求出的阈值。遍历不同的阈值，计算不同阈值下对应的背景和前景之间的类间方差，当类间方差取得极大值时，此时对应的阈值就是OTSU算法所求的阈值。
2.1.3平滑去噪
在进行骨骼精分割以后，图像中也存在少许噪点。此外，骨骼边缘比较粗糙，会直接影响后续点云模型的精度。所以我们还需要进行平滑处理。
步骤1形态学开运算去噪：可选择小半径的圆形结构元，对精分割结果进行开运算。可在保持骨骼不变的情况下去除小噪点，并断开不同骨骼间的伪连接。
步骤2高斯滤波平滑边缘：选择小半径高斯滤波器对去噪结果进行平滑，然后以阈值T=0.5分割平滑结果。
2.2求取缺损部分模型
求取缺损模型首先应该是在适配模型库中找到与缺损模型最相近的完整模型，由于我们需要得到最契合的缺损模型，因此本文提出的方法是在获取模型点云的基础上进行配准的。利用MarchingCubes算法获取体模表面三角片，并计算三角片重心作为点云集。接着进行点云降采样，这样做的目的是减少点云数量，提高配准速率。
2.2.1点云粗配准
在得到匹配模板和缺损模型的点云后，利用最小包围盒算法进行点云的粗配准。最小包围盒也叫外接最小矩阵，是一种求解离散点集最优包围空间的算法，基本思想是用体积稍大且特性简单的几何体（称为包围盒）来近似地代替复杂的几何对象。
首先，获取已知点云的三个主方向，获取质心，计算协方差，获得协方差矩阵，求取协方差矩阵的特征值和特长向量，特征向量即为主方向。接着，利用已获得的主方向和质心，将输入点云转换至原点，且主方向与坐标系方向重回，建立变换到原点的点云的包围盒。最后，给输入点云设置主方向和包围盒，通过输人点云到原点点云变换的逆变换实现。
通俗来说，最小包围盒算法是求解点集XOY平面内投影点的最小面积包围盒，即把空间化问题转化为平面问题来求解。
2.2.2点云精配准
由兰州交通大学杨军等人提出的改进的ICP算法是本文用于点云精配准的主要方法，在该方法中主要在初始值的敏感性以及整体与部分模型的配准精度上做出了改进。选择改进的ICP算法适用于精配准是由于在初始配准后还未达到想要的效果，因此选用KD-Tree搜索大幅度地提高了搜索效率，制定自适应最有阈值则是在经典的ICP算法上提高了配准精度。
假如我们手上的数据有原始点云数据P和经过平移加旋转后的点云Q，那么ICP算法就是求出最优的从P到Q的坐标变换矩阵。从（2）中我们可以看到我们从P到Q的过程，其中旋转为R.平移为T。
2.3计算缺损部分
为解决无法通过点云直接得到缺损部分的问题，我们以点云模型为基准建立初始采样网格，后利用变换矩阵对网格进行逆运算，再对原始体数据进行采样，即可得到配准后的体模型。
把适配模型网格称为模板，缺损模型网格称为实例，那么现在用模板减去实例就可以得到粗结果。当然我们应该在粗结果上进行进行错误部分的提出，其中包括删除扁平、细长、颗粒状部分，最后我们需要利用三维高斯核和球形中值滤波器对剩余部分进行平滑处理。
3实验
在本实验中选用羊的胫骨作为实验对象，数据集一共83张516KB的CT图像，CT数据如图2所示。读取整个序列，并根据dicom参数（ResaleSlope、RescaleIntercept、WindowCenter、WindowWidth、ShceLocationl进行调窗与调序，获得数值为0～255的灰度图像;接着利用阈值分割提取所有骨骼，利用连通分量筛选目标骨骼;再根据dicom参数（PixelSpacing、SliceTh-ickness），重建序列体数据，过程示例如图3所示。
在选取到适匹配模板后，那就可以生成缺损模型和匹配到的完整模型的点云，如图4所示，经过点云的粗配准即最小包围盒算法以及精配准ICP算法后的点云如图5所示，在经过这一步操作后我们可以得到变换矩阵，这一矩阵的得到是为了下一步做铺垫。
在得到缺损模型前，需要用采样网格对点云进行包围，在图6中三维坐标下竖直采样网格代表缺损模型，倾斜采样网格代表匹配到的完整模型，在已知变换矩阵的条件下对采样网格进行逆运算，从而得到配准后模型。
那么在这样的情况下，就可以进行两模型的直接相减得到粗略的缺损模型，如图7所示，在经过剔除错误部分和平滑后，就可以生成打印模型。MarchingCubes面绘制算法获取三角片，并以光顺优化获得更好的三角片，最后以sTL的格式导出，结果的缺损模型如图8所示。