CN107330281B - 全自动骨折钢板模型个性化重构方法 - Google Patents
全自动骨折钢板模型个性化重构方法 Download PDFInfo
- Publication number
- CN107330281B CN107330281B CN201710542904.0A CN201710542904A CN107330281B CN 107330281 B CN107330281 B CN 107330281B CN 201710542904 A CN201710542904 A CN 201710542904A CN 107330281 B CN107330281 B CN 107330281B
- Authority
- CN
- China
- Prior art keywords
- point
- model
- points
- broken bone
- bone
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 78
- 229910000831 Steel Inorganic materials 0.000 title claims abstract description 32
- 239000010959 steel Substances 0.000 title claims abstract description 28
- 208000010392 Bone Fractures Diseases 0.000 claims abstract description 133
- 210000000988 bone and bone Anatomy 0.000 claims abstract description 60
- 239000013598 vector Substances 0.000 claims abstract description 49
- 239000011159 matrix material Substances 0.000 claims abstract description 33
- 206010017076 Fracture Diseases 0.000 claims abstract description 31
- 238000000926 separation method Methods 0.000 claims abstract description 8
- 208000006670 Multiple fractures Diseases 0.000 claims abstract description 7
- 238000005452 bending Methods 0.000 claims abstract description 5
- 230000008719 thickening Effects 0.000 claims abstract description 4
- 230000008569 process Effects 0.000 claims description 27
- 241000270295 Serpentes Species 0.000 claims description 19
- 230000000694 effects Effects 0.000 claims description 15
- 230000011218 segmentation Effects 0.000 claims description 15
- 238000000605 extraction Methods 0.000 claims description 13
- 230000000644 propagated effect Effects 0.000 claims description 10
- 238000012935 Averaging Methods 0.000 claims description 7
- 230000009466 transformation Effects 0.000 claims description 7
- 238000009616 inductively coupled plasma Methods 0.000 claims description 5
- 230000007480 spreading Effects 0.000 claims description 5
- 230000005489 elastic deformation Effects 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 2
- 238000010586 diagram Methods 0.000 description 14
- 238000002474 experimental method Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000000338 in vitro Methods 0.000 description 3
- 241001442234 Cosa Species 0.000 description 2
- 241001347978 Major minor Species 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 208000014674 injury Diseases 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 208000027418 Wounds and injury Diseases 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000000740 bleeding effect Effects 0.000 description 1
- 230000006378 damage Effects 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
- 238000013519 translation Methods 0.000 description 1
- 230000008733 trauma Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Graphics (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Generation (AREA)
Abstract
本发明公开了一种全自动骨折钢板模型个性化重构方法，包括如下步骤：断骨图像分离步骤：断骨三维重构步骤：断骨轴线提取步骤：生成所述断骨一和断骨二三维模型的初始样本矩阵Z；对样本矩阵进行中心化，得到矩阵X；计算该矩阵X的协方差矩阵C和协方差矩阵C的特征值，选取最大特征值对应的特征向量作为所述断骨一和断骨二的断骨轴线；断骨截面分割步骤：在所述的断骨一和断骨二的三维模型上提取特征点；将特征点按长度缓冲值以及网格模型划分为首尾两部分；分别形成断骨一和断骨二首尾特征点的共4个集合；断骨配准步骤：虚拟钢板预弯步骤：在拼接好断骨断面附件进行点选，确定钢板模型尺寸和形状；计算每个三角面片的法向量值，每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位，最终得到模拟钢板三维数据。
Description
技术领域
本发明涉及全自动骨折钢板模型个性化重构方法。涉及专利分类号G06计算；推算；计数G06T一般的图像数据处理或产生G06T7/00图像分析，例如从位像到非位像。
背景技术
目前骨折手术普遍采用人工复位和伤肢内固定相结合的方法，这种方法存在的问题是创伤大、出血多、容易引发神经血管损伤等并发症。因此，我们可以利用计算机对断骨模型进行虚拟拼接，从而在术前得到钢板的各种几何参数。然而目前的断骨虚拟拼接方法存在很多问题，比如断面分割不准确、预配准不精确，这些问题都很大程度地影响了断骨模型配准的精确度。
现有技术中主要的技术手段主要包括：根据断面法向量与轴线的夹角来进行断面分割以及单纯根据断骨的轴线进行预配准的方法。
根据断面法向量与轴线的夹角进行断面分割的方法，鲁棒性差，遇到结构复杂的断面时不能做到很好地分割，而且需要手动点选，不能做到自动分割。
单纯将断骨的主轴线进行对齐并不能满足预配准的要求，可能会出现断骨的截面方向相背离的情况，而且截面形状也无法做到大致吻合，很容易造成精配准的不准确。
发明内容
本发明要解决的技术问题是：一种全自动骨折钢板模型个性化重构方法，包括如下步骤：
断骨图像分离步骤：
对采集的断骨CT图片序列进行至少包括二值化的预处理；提取二值化后 CT图片序列中的每副图像的断骨轮廓；
根据该断骨轮廓，应用snake模型完成CT图片中的断骨一和断骨二图像区域分离；分别得到断骨一的图片序列集和断骨二的图片序列集；
断骨三维重构步骤：
以二值化后的CT图片序列，三维空间中的相邻像素点为顶点，构建多个正方体体元；
根据图像中的灰度值，将图像中的像素点划分为实点和虚点，将所述正方体体元中既包含实点和虚点的体元定义为边界体元；取边界体元中连接实点和虚点的边的重点作为等值面三角形的顶点，以此方法在所述的边界体元中构建等值面三角形；
提取所述两个断骨的图片序列集中的全部等值面三角形，分别组成断骨一和断骨二的三维模型；
断骨轴线提取步骤：
生成所述断骨一和断骨二三维模型的初始样本矩阵Z；对样本矩阵进行中心化，得到矩阵X；计算该矩阵X的协方差矩阵C和协方差矩阵C的特征值，选取最大特征值对应的特征向量作为所述断骨一和断骨二的断骨轴线；
断骨截面分割步骤：
在所述的断骨一和断骨二的三维模型上提取特征点；将特征点按长度缓冲值以及网格模型划分为首尾两部分；分别形成断骨一和断骨二首尾特征点的共4 个集合；
利用所述的4个特征点集合，建立包围盒；通过分析比较包围盒三轴方差值，确定表示断骨断面所在的2个包围盒，并进行包围盒对准；对准后，进行点集搜索，得到两个断骨模型截面的点集P1和P2；
断骨配准步骤：
该步骤包括以通过坐标轴旋转，拉近所述断骨一和断骨二距离，以及完成截面大致对齐的预配准步骤和基于ICP算法对所述两断骨模型点集进行精确配准的精确配准步骤；
虚拟钢板预弯步骤：
在拼接好断骨断面附件进行点选，确定钢板模型尺寸和形状；记录点选三角平面值，选中范围内所有表面三角面片，计算每个三角面片的法向量值，每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位，最终得到模拟钢板三维数据。
作为优选的实施方式，所述的snake模型完成CT图像中的断骨图像自动分离具体过程如下：
分别选取CT图像集的第一张和最后一张的轮廓线作为正向和反向分割的初始曲线；
在轮廓线的控制点上定义snake能量函数
其中，snake能量函数中前两项：
为内部力，用来控制轮廓线的弹性形变；第三项
为定义的曲线的外部力，
表示v点处的灰度梯度，用来控制曲线的位置与局部特征吻合；
设置每张图像中snake算法的迭代次数，迭代过程中轮廓的控制点在局部范围内搜索，使总能量Etotal趋近于极小值；
从CT图像集的首部至尾部迭代来分离出断骨一，再反向重新迭代分离出断骨二。
更进一步的，分离断骨图像和背景之后，还具有图像去噪步骤：
将重合的像素提取出来，存入像素点集S；
对S中的第i个像素点，分别获取相邻10张图像中该像素邻近位置(5x5 的矩形区域)中像素点的灰度值g，若g>140，则记录为有效灰度值，对全部10
张图片中的有效灰度值求和得到Gi，定义像素点集S与该图像集的相关值
其中，n为S中像素点的数量；
分别计算断骨一和断骨二两个图像集中重合部分的相关值M1和M2，在相关值小的图像集中舍去重合的像素点集。
作为优选的实施方式，对断骨进行三维重构的过程如下：
以图像集的相邻像素点为顶点构建w×h×n个正方体体元w：图像宽度，h：图像高度，n为图像数量；
将图像中灰度值大于阈值的像素点规定为实点，反之为虚点，顶点中既包含实点又包含虚点的体元为边界体元；
取边界体元中连接实点和虚点的边的中点作为等值面三角形的顶点，以此方法在边界体元中构建等值面三角形；
分别提取两个图像集中所有的等值面三角形构成两个断骨三维模型。
更进一步的，构建得到断骨模型后，还具有模型降噪步骤：
随机选择模型中一个未遍历过的面片作为起始面片，用递归的方法对相邻面片进行蔓延；
若蔓延结束时，蔓延到的面片数量超过阈值T；
将蔓延到的面片全部加入已蔓延集合中，否则，加入到孤立面片集合中；
重复上述步骤直到所有面片均被蔓延，然后删除孤立面片集合中所有的面片及其顶点，去噪完成。
作为优选的实施方式，断骨截面分割步骤中特征点提取过程下：
计算每个点的曲率α和每个面片的法向量
并计算
与主轴线l的夹角余弦值β，当满足以下两个条件中任意一个时，即提取该点为断面分割的特征点：该点的曲率α>0.9；该点所在的面片β>0.35。
更进一步的，特征点提取后，还具有特征点去噪步骤：
移动两个断骨模型的空间位置，使断骨轴线与x轴平行，方便后续沿坐标轴划分网格；
利用网格法寻找噪音点，将整个模型点集划分为a×a×b个网格，沿主轴线方向进行b等分，将整个网格模型建立一个三维坐标系，从(0,0,0)到(a,a,b)每个网格都可以用一个三维坐标(xi,yi,zi)表示，并统计每个网格中点的数量ni；
设定网格中点数目的阈值t，t＝5，若0<ni<t，则将此网格i列为待定网格；
遍历所有网格，若网格j满足nj≥t，则计算此网格j与当前待定网格i的直线平方距离d2＝(xi-xj)2+(yi-yj)2+(zi-zj)2，遍历完成后找出数值最小的5个平方距离并求平均值
设定平均距离的阈值r，若
则网格i中的点为好点，否则，为噪音点；
重复上述步骤，在待定网格中找出所有噪音点并去除，去噪效果如图8所示。
更进一步的，特征点划分步骤如下：
设定一个长度缓冲值m＝l/2.2166，l为断骨模型沿主轴线方向的长度，m的作用是为后续的点集划分过程提供一个弹性缓冲，并根据断骨长度的不同提供不同的缓冲值；
利用模型去噪过程中构造的a×a×b的网格模型，沿主轴线进行了b等分，每个单位长度含有a×a个网格，统计这a2个网格中点的总数目n，若n小于阈值r，则m-1，若n大于阈值r，则初始化缓冲值m；
上述过程结束之后，进行下一组a2个网格的检测，循环此过程直到m<0或者遍历完全部b个网格组，经过上述过程遍历到的点组成一个新的点集；
按上述步骤从0到b-1遍历一遍得到点集S1，再反向遍历一遍得到点集S2，此时两个点集可能会有重合部分，即
l1l2分别是两个点集的长度，l是断骨模型的长度；
将两个点集的重合部分进行了均分，均分方法：
l0＝(l1+l2-l)/2
l1′＝l-l2+l0
l2′＝l-l1+l0
l1′和l2′是首尾两个点集最终的长度。
建立包围盒。利用包围盒算法分别计算四个点集(每个断骨模型包含首尾两
个点集)的OBB包围盒，如图10，并记录长轴、中轴和短轴的方向和长度；
截面包围盒自动选取。求两个包围盒三轴方差的公式为：
l1l2l3代表包围盒1的长中短轴的长度；
l1′l2′l3′代表包围盒2的长中短轴的长度；
将模型一中的两个包围盒与模型二中的两个包围盒进行两两比较，共计算得到4个方差值，选取方差最小的一组，即为两个断面所在的包围盒，如图11 即为提取出来的两个包围盒。
作为优选的实施方式，所述的包围盒对准过程如下：判断两个断骨模型的贴近一端是否是断面，若不是则进行调整，使断面相互贴近；
提取两个包围盒的三轴在y-o-z平面上的投影方向；
在y-o-z平面上计算长轴和中轴的角分线方向向量
以x轴增大的方向为正方向，将短轴统一为正方向，提取短轴在y-o-z平面上的投影向量；
固定模型一，将模型二以1°为最小单位绕主轴线进行旋转遍历。设旋转的角度为α，则模型二从初始位置进行旋转的旋转矩阵为：
旋转后的向量为
V1′＝AV1
V2′＝AV2
求旋转变换后，两个包围盒对应轴之间的夹角余弦值；
D＝cosa+cosb
对断骨2进行360°的旋转遍历，D的值最大时对应的角度即为包围盒配准时断骨2所需旋转的角度。
作为优选的实施方式，截面点集搜索步骤，
对包围盒对准后的模型重新计算面片法向量；提取包含该点的所有面片的法向量vi＝(xi,yi,zi)，将所有法向量的长度化为1；
设该数据点处的法向量vp＝(xp,yp,zp)；
n为包含该数据点的面片数量，yp，zp同理；遍历提取出来的所有数据点，计算该点处的法向量与断骨轴线的夹角β，若cosβ>0.35则直接将该点提取至断面点集，否则进行如下步骤；
设置坐标距离阈值r1＝20.0，向量阈值r2＝1.0，对于当前研究的数据点p，在另一个断骨模型中遍历包围盒内数据点p′，同时满足以下三个条件则视该点p′为点p的匹配点：
坐标的平方距离小于阈值，即(x1-x2)2+(y1-y2)2+(z1-z2)2<r1； (xp+xp′)2+(yp+yp′)2+(zp+zp′)2<r2，(xp,yp,zp)(xp′,yp′,zp′)分别是点p和p′处的法向量，点p和p′处的法向量夹角为γ，cosγ>0；
若对点p，能找到至少5个p′满足上述条件，则将p加入至断面点集。(5) 分别对断骨模型1和断骨模型2中的包围盒内数据点进行上述遍历，最终提取得到两个断面/截面点集P1和P2。
附图说明
为了更清楚的说明本发明的实施例或现有技术的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做一简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明的算法流程图
图2为本发明中CT图像自动分离效果示意图
图3为本发明中二值化及轮廓提取效果示意图
图4为本发明中对断骨进行三维重构的结果示意图
图5为本发明断骨轴线提取效果示意图
图6为本发明特征点提取结果示意图
图7为本发明中噪音点示意图
图8为本发明中特征点去噪效果示意图
图9为本发明特征点划分结果示意图
图10为本发明中建立的包围盒示意图
图11为本发明中截面包围盒自动提取示意图
图12为本发明中包围盒对准示意图
图13为本发明中截面点集搜索结果示意图
图14为本发明中断骨预配准步骤示意图
图15为本发明算法步骤中断骨精配准结果示意图
图16为本发明算法中点选钢板控制点示意图
图17为本发明最终拟合的钢板示意图
具体实施方式
为使本发明的实施例的目的、技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：
如图1-17所示：全自动骨折钢板模型个性化重构方法本方法基于snake模型进行断骨CT图像的自动分离。分离效果如图2所示。
图像二值化。导入的CT图像为灰度图，设置一个阈值T，大于T的像素灰度值设为255，小于等于T的像素灰度值设为0，经过大量实验，T的值为130 时效果最佳，二值化的效果如图3。
提取轮廓。在应用snake算法之前进行轮廓提取的好处有两点：一、可以用最初的轮廓线作为snake初始包围曲线，无需人工勾画轮廓线，实现了全自动分割；二、只专注于断骨外部轮廓，避免了内部像素对snake模型分割效果的影响。此处采用基于OpenCV的轮廓提取技术来提取断骨的外轮廓，提取效果如图3。
应用snake模型自动分割。(1)分别选取CT图像集的第一张和最后一张的轮廓线作为正向和反向分割的初始曲线；(2)在轮廓线的控制点上定义snake能量函数
其中，snake能量函数中前两项：
为内部力，用来控制轮廓线的弹性形变；第三项
为定义的曲线的外部力，
表示v点处的灰度梯度，用来控制曲线的位置与局部特征吻合。
三个力共同作用使能量函数趋于极小化，最终可以分割出图像的局部特征；(3)能量函数中α决定曲线的伸展能力，β决定曲线的弯曲能力，经过大量的实验，取α＝0.05，β＝0.25；(4)设置每张图像中snake算法的迭代次数为1000次，迭代过程中轮廓的控制点在局部范围内搜索，使总能量Etotal趋近于极小值。Snake演变过程是从CT图像集的首部至尾部迭代来分离出断骨一，再反向重新迭代分离出断骨二。
图像去噪。分割结果中的噪音主要出现在骨骼断面的部分，可能会因为分割不完全而导致两个断骨的图像数据出现重合的部分，因此在分割完成之后进行去噪处理是十分必要的。去噪的目的是将重合部分进行划分，避免在两个断骨数据中出现耦合，方法如下：(1)将重合的像素提取出来，存入像素点集S；(2) 对S中的第i个像素点，分别获取相邻10张图像中该像素邻近位置(5×5的矩形区域)中像素点的灰度值g，若g>140，则记录为有效灰度值，对全部10张图片中的有效灰度值求和得到Gi，定义像素点集S与该图像集的相关值。
n为S中像素点的数量；(3)分别计算断骨一和断骨二两个图像集中重合部分的相关值M1和M2，在相关值小的图像集中舍去重合的像素点集。
基于MC算法的断骨三维重构
本方法基于MarchingCubes算法对断骨的CT图像进行三维重构，重构的结果如图4所示。
利用MC算法进行三维重构，方法如下：(1)以图像集的相邻像素点为顶点构建w×h×n个正方体体元(w：图像宽度，h：图像高度，n：图像数量)；(2)将图像中灰度值大于150的像素点规定为实点，反之为虚点，顶点中既包含实点又包含虚点的体元为边界体元，边界体元也是三维重构中的重点研究对象；(3) 取边界体元中连接实点和虚点的边的中点作为等值面三角形的顶点，以此方法在边界体元中构建等值面三角形；(4)分别提取两个图像集中所有的等值面三角形构成两个断骨三维模型。
三维模型去噪。重构得到的三维模型难免会有噪音，噪音的来源主要是体外孤立的面片，去噪的方法是用递归来寻找体外孤立面片然后去除，方法如下： (1)随机选择模型中一个未遍历过的面片作为起始面片，用递归的方法对相邻面片进行蔓延；(2)若蔓延结束时，蔓延到的面片数量超过阈值T(经过大量实验，设置阈值T为10)，则将蔓延到的面片全部加入已蔓延集合中，否则，加入到孤立面片集合中；(3)重复上述步骤直到所有面片均被蔓延，然后删除孤立面片集合中所有的面片及其顶点，去噪完成。
基于PCA算法的断骨轴线提取
本方法利用PCA算法对两个断骨模型提取轴线，提取效果如图5所示。
计算点集中心点坐标。分别计算断骨模型中所有点的X、Y、Z坐标的平均值，得到断骨点集中心点的坐标p0(x0,y0,z0)。
特征中心化。将断骨点集的空间坐标构建为3×n矩阵
n为点的总数量，分别计算矩阵中每个点的坐标与中心点坐标p0的差值，得到更新后的矩阵A，即
计算协方差矩阵。将矩阵A与其转置矩阵A′相乘，得到协方差矩阵M，即 M＝AA′
求协方差矩阵M的特征值和特征向量。最大的特征值对应的向量即为该模型点集的主方向，即断骨模型主轴线的方向向量。
基于包围盒的断骨截面自动分割
提取特征点。计算每个点的曲率α和每个面片的法向量
并计算
与主轴线l的夹角余弦值β，当满足以下两个条件中任意一个时，即提取该点为断面分割的特征点：(1)该点的曲率α>0.9；(2)该点所在的面片β>0.35。特征点提取的结果如图6所示。
特征点去噪。噪音的来源主要是模型中部零星的特征点，如图7所示，这些噪音会对后续的包围盒建立造成很大误差，因此需要进行去噪处理，方法如下：
(1)移动两个断骨模型的空间位置，使断骨轴线与x轴平行，方便后续沿坐标轴划分网格；(2)利用网格法寻找噪音点，将整个模型点集划分为a×a×b个网格(沿主轴线方向进行b等分)，将整个网格模型建立一个三维坐标系，从(0,0,0) 到(a,a,b)每个网格都可以用一个三维坐标(xi,yi,zi)表示，并统计每个网格中点的数量ni；(3)设定网格中点数目的阈值t(经过大量实验，设置t为5)，若0<ni<t，则将此网格i列为待定网格；(4)遍历所有网格，若网格j满足nj≥t，则计算此网格j与当前待定网格i的直线平方距离d2＝(xi-xj)2+(yi-yj)2+(zi-zj)2，遍历完成后找出数值最小的5个平方距离并求平均值
(5)设定平均距离的阈值r(本方法中设为15.0)，若
则网格i中的点为好点，否则，为噪音点；(6)重复上述步骤，在待定网格中找出所有噪音点并去除，去噪效果如图8所示。
特征点划分。划分特征点的目的是将去噪后的特征点集划分为首尾两个部分，便于后续自动分离出断面点集。划分方法如下：(1)设定一个长度缓冲值 m＝l/2.2166，l为断骨模型沿主轴线方向的长度，m的作用是为后续的点集划分过程提供一个弹性缓冲，并根据断骨长度的不同提供不同的缓冲值；(2)仍然利用去噪过程中构造的a×a×b的网格模型，沿主轴线进行了b等分，每个单位长度含有a×a个网格，统计这a2个网格中点的总数目n，若n小于阈值r(设为3)，则m-1，若n大于阈值r，则初始化缓冲值m。上述过程结束之后，进行下一组 a2个网格的检测，循环此过程直到m<0或者遍历完全部b个网格组，经过上述过程遍历到的点组成一个新的点集；(3)按上述步骤从0到b-1遍历一遍得到点集S1，再反向遍历一遍得到点集S2，此时两个点集可能会有重合部分(即l1+l2>l， l1l2分别是两个点集的长度，l是断骨模型的长度)，会影响后续包围盒的准确度，因此本方法将两个点集的重合部分进行了均分，均分方法：
l0＝(l1+l2-l)/2
l1′＝l-l2+l0
l2′＝l-l1+l0
l1′和l2′是首尾两个点集最终的长度，特征点划分的结果如图9所示。
建立包围盒。利用包围盒算法分别计算四个点集(每个断骨模型包含首尾两个点集)的OBB包围盒，如图10，并记录长轴、中轴和短轴的方向和长度。
截面包围盒自动选取。求两个包围盒三轴方差的公式为
l1l2l3代表包围盒1的长中短轴的长度
l1′l2′l3′代表包围盒2的长中短轴的长度
将模型一中的两个包围盒与模型二中的两个包围盒进行两两比较，共计算得到4个方差值，选取方差最小的一组，即为两个断面所在的包围盒，如图11 即为提取出来的两个包围盒。
包围盒对准。为了减弱包围盒的准确度对配准精度的影响，提高轴线配准的作用，因此包围盒对准是建立在轴线对准的基础上，即两个断骨模型绕轴线进行旋转来达到最佳的匹配效果。方法如下：(1)先判断两个断骨模型的贴近一端是否是断面，若不是则进行调整，使断面相互贴近(2)提取两个包围盒的三轴在y-o-z平面上的投影方向(此时两个模型轴线与x轴重合，只研究包围盒在y-o-z 平面上的配准)，为了解决某些情况下因为长轴和中轴长度接近而出现的颠倒问题，本方法不直接用三轴进行对准，而是采用长轴和中轴的角分线，再综合短轴来进行对准。在y-o-z平面上计算长轴和中轴的角分线方向向量
以x轴增大的方向为正方向，将短轴统一为正方向，
提取短轴在y-o-z平面上的投影向量+
固定模型一，将模型二以1°为最小单位绕主轴线进行旋转遍历。设旋转的角度为α，则模型二从初始位置进行旋转的旋转矩阵为
旋转后的向量为
V1′＝AV1
V2′＝AV2
(4)求旋转变换后，两个包围盒对应轴之间的夹角余弦值，
D＝cosa+cosb
对断骨2进行360°的旋转遍历，D的值最大时对应的角度即为包围盒配准时断骨2所需旋转的角度，包围盒对准结果如图12所示。
截面点集搜索。将包围盒进行对准的目的是拉近两个截面中对应点之间的距离，便于利用搜索的方法进行截面提取。(1)对包围盒对准后的模型重新计算面片法向量；(2)将两个包围盒中所有的数据点提取出来作为研究对象，计算这些数据点处的法向量，方法是：提取包含该点的所有面片的法向量vi＝(xi,yi,zi)。
将所有法向量的长度化为1，
yizi同理。设该数据点处的法向量vp＝(xp,yp,zp)
(3)遍历提取出来的所有数据点，计算该点处的法向量与断骨轴线的夹角β，若cosβ>0.35则直接将该点提取至断面点集，否则进行步骤(4)；(4)设置坐标距离阈值r1＝20.0，向量阈值r2＝1.0，对于当前研究的数据点p，在另一个断骨模型中遍历包围盒内数据点p′，同时满足以下三个条件则视该点p′为点p的匹配点：
①坐标的平方距离小于阈值，即(x1-x2)2+(y1-y2)2+(z1-z2)2<r1
②(xp+xp′)2+(yp+yp′)2+(zp+zp′)2<r2，
(xp,yp,zp)(xp′,yp′,zp′)分别是点p和p′处的法向量
③点p和p′处的法向量夹角为γ，cosγ>0。
若对点p，能找到至少5个p′满足上述条件，则将p加入至断面点集。(5) 分别对断骨模型1和断骨模型2中的包围盒内数据点进行上述遍历，最终提取得到两个断面点集P1和P2。点集搜索的结果如图13所示。
断骨模型预配准
在前面的步骤中两个断骨模型的主轴线已经对准(全部和x轴重合)，断面的点集已经提取完成，本步骤中利用断面的点集进行预配准，目的是使两个断骨模型的断面大致吻合。
通过PCA算法分别计算两个断面点集的主方向，以x轴增大的方向为正方向，将两个主方向向量都调整为正方向。
提取两个主方向向量在y-o-z平面上的投影向量v1v2，固定断骨2，将断骨1 以主轴线为轴进行旋转，直至v1v2平行，此时两个断面在主轴线方向上可以做到大致吻合。
缩短断骨之间的间距，便于后续的精配准，至此，预配准过程结束。预配准的结果如图14所示。
断骨模型精配准
预配准结束后，利用ICP算法对两组断骨模型点集进行精配准操作，进一步提高其准确性。
将两部分点集分别记为U与P。
计算最近点，即对于集合U中的每一个点，在集合P中都找出距该点最近的对应点，设集合P中由这些对应点组成的新点集为Q＝{qi,i＝0,1,2,...,n}。
采用最小均方根法，计算点集U与Q之间的配准，使得到配准变换矩阵R， T，其中R是3×3的旋转矩阵，T是3×1的平移矩阵。
计算坐标变换，即对于集合U，用配准变换矩阵R，T进行坐标变换，得到新的点集U1，即U1＝RU+T。
计算U1与Q之间的均方根误差，如小于预设的极限值ε，则结束，否则，以点集U1替换U，重复上述步骤。最终精配准的结果如图15所示。
虚拟钢板预弯
在本步骤，主要利用对断骨断裂部位曲面数据的相应操作来模拟出钢板的大概形态。
用户在拼接好的断骨断面附近进行点选(如图16)，确定钢板模型的大概形状以及大小。
记录下点选的三角平面值，选中范围内的所有表面三角面片。
计算每个三角面片的法向量值，并记录下来。
将每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位。
所得到的加厚部分即为模拟的钢板模型三维数据(如图17)，可将其作为结果导出输出。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。
Claims (10)
1.一种全自动骨折钢板模型个性化重构方法，其特征在于包括如下步骤：
断骨图像分离步骤：
对采集的断骨CT图像序列进行至少包括二值化的预处理；提取二值化后CT图像序列中的每幅图像的断骨轮廓；
根据该断骨轮廓，应用snake模型完成CT图像中的断骨一和断骨二图像区域分离；分别得到断骨一的图片序列集和断骨二的图片序列集；
断骨三维重构步骤：
以二值化后的CT图像序列，三维空间中的相邻像素点为顶点，构建多个正方体体元；
根据像素图像中的灰度值，将图像中的像素点划分为实点和虚点，将所述正方体体元中既包含实点又包含虚点的体元定义为边界体元；取边界体元中连接实点和虚点的边的中点作为等值面三角形的顶点，以此方法在所述的边界体元中构建等值面三角形；
提取两个断骨的图像序列集中的全部等值面三角形，分别组成断骨一和断骨二的三维模型；
断骨轴线提取步骤：
生成所述断骨一和断骨二三维模型的初始样本矩阵Z；对样本矩阵进行中心化，得到矩阵X；计算该矩阵X的协方差矩阵C和协方差矩阵C的特征值，选取最大特征值对应的特征向量作为所述断骨一和断骨二的断骨轴线；
断骨截面分割步骤：
在所述的断骨一和断骨二的三维模型上提取特征点；将特征点按长度缓冲值以及网格模型划分为首尾两部分；分别形成断骨一和断骨二首尾特征点的共4个特征点集合；
利用所述的4个特征点集合，建立包围盒；通过分析比较包围盒三轴方差值，确定表示断骨断面所在的2个包围盒，并进行包围盒对准；对准后，进行点集搜索，得到两个断骨模型截面的点集P1和P2；
断骨配准步骤：
该步骤包括以通过坐标轴旋转，拉近所述断骨一和断骨二距离，以及完成截面大致对齐的预配准步骤和基于ICP算法对两断骨模型点集进行精确配准的精确配准步骤；
虚拟钢板预弯步骤：
在拼接好断骨断面附近进行点选，确定钢板模型尺寸和形状；记录点选三角平面值，选中范围内所有表面三角面片，计算每个三角面片的法向量值，每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位，最终得到模拟钢板三维数据。
2.根据权利要求1所述的全自动骨折钢板模型个性化重构方法，其特征还在于所述的snake模型完成CT图像中的断骨一和断骨二图像区域自动分离具体过程如下：
分别选取CT图像集的第一张和最后一张的轮廓线作为正向和反向分割的初始曲线；
在轮廓线的控制点上定义snake能量函数
其中，snake能量函数中前两项：
为内部力，用来控制轮廓线的弹性形变；第三项
为定义的曲线的外部力，
表示v点处的灰度梯度，用来控制曲线的位置与局部特征吻合；
设置每张图像中snake算法的迭代次数，迭代过程中轮廓的控制点在局部范围内搜索，使总能量Etotal趋近于极小值；
从CT图像集的首部至尾部迭代来分离出断骨一的序列图像集，再反向重新迭代分离出断骨二的序列图像集。
3.根据权利要求2所述的全自动骨折钢板模型个性化重构方法，其特征还在于分离断骨图像之后，还具有图像去噪步骤：
将重合的像素提取出来，存入像素点集S；
对S中的第i个像素点，分别获取相邻10张图像中该像素邻近位置中像素点的灰度值g，若g>140，则记录为有效灰度值，对全部10张图片中的有效灰度值求和得到Gi，定义像素点集S与该图像集的相关值
其中，n为S中像素点的数量；
分别计算断骨一和断骨二两个图像集中重合部分的相关值M1和M2，在相关值小的图像集中舍去重合的像素点集。
4.根据权利要求1所述的全自动骨折钢板模型个性化重构方法，其特征还在于：对断骨进行三维重构的过程如下：
以图像集的相邻像素点为顶点构建w×h×n个正方体体元，w：图像宽度，h：图像高度，n为图像数量；
将图像中灰度值大于阈值的像素点规定为实点，反之为虚点，顶点中既包含实点又包含虚点的体元为边界体元；
取边界体元中连接实点和虚点的边的中点作为等值面三角形的顶点，以此方法在边界体元中构建等值面三角形；
分别提取两个图像集中所有的等值面三角形构成两个断骨三维模型。
5.根据权利要求4所述的全自动骨折钢板模型个性化重构方法，其特征还在于构建得到断骨模型后，还具有模型降噪步骤：
随机选择模型中一个未遍历过的面片作为起始面片，用递归的方法对相邻面片进行蔓延；
若蔓延结束时，蔓延到的面片数量超过阈值T；
将蔓延到的面片全部加入已蔓延集合中，否则，加入到孤立面片集合中；
重复上述步骤直到所有面片均被蔓延，然后删除孤立面片集合中所有的面片及其顶点，去噪完成。
6.根据权利要求1所述的全自动骨折钢板模型个性化重构方法，其特征还在于断骨截面分割步骤中特征点提取过程下：
计算每个点的曲率α和每个面片的法向量
并计算
与主轴线l的夹角余弦值β，当满足以下两个条件中任意一个时，即提取该点为断面分割的特征点：该点的曲率α>0.9；该点所在的面片β>0.35。
7.根据权利要求6所述的全自动骨折钢板模型个性化重构方法，其特征还在于特征点提取后，还具有特征点去噪步骤：
移动两个断骨模型的空间位置，使断骨轴线与x轴平行，方便后续沿坐标轴划分网格；
利用网格法寻找噪音点，将整个模型点集划分为a×a×b个网格，沿主轴线方向进行b等分，将整个网格模型建立一个三维坐标系，从(0,0,0)到(a,a,b)每个网格都可以用一个三维坐标(xi,yi,zi)表示，并统计每个网格中点的数量ni；
设定网格中点数目的阈值t，t＝5，若0<ni<t，则将此网格i列为待定网格；
遍历所有网格，若网格j满足nj≥t，则计算此网格j与当前待定网格i的直线平方距离d2＝(xi-xj)2+(yi-yj)2+(zi-zj)2，遍历完成后找出数值最小的5个平方距离并求平均值
设定平均距离的阈值r，若
则网格i中的点为好点，否则，为噪音点；
重复上述步骤，在待定网格中找出所有噪音点并去除，去噪效果如图8所示。
8.根据权利要求7所述的全自动骨折钢板模型个性化重构方法，其特征还在于特征点划分步骤如下：
设定一个长度缓冲值m＝l/2.2166，l为断骨模型沿主轴线方向的长度，m的作用是为后续的点集划分过程提供一个弹性缓冲，并根据断骨长度的不同提供不同的缓冲值；
利用模型去噪过程中构造的a×a×b的网格模型，沿主轴线进行了b等分，每个单位长度含有a×a个网格，统计这a2个网格中点的总数目n，若n小于阈值r，则m-1，若n大于阈值r，则初始化缓冲值m；
上述过程结束之后，进行下一组a2个网格的检测，循环此过程直到m<0或者遍历完全部b个网格组，经过上述过程遍历到的点组成一个新的点集；
按上述步骤从0到b-1遍历一遍得到点集S1，再反向遍历一遍得到点集S2，此时两个点集可能会有重合部分，即l1+l2>l，l1 l2分别是两个点集的长度，l是断骨模型的长度；
将两个点集的重合部分进行了均分，均分方法：
l0＝(l1+l2-l)/2
l1′＝l-l2+l0
l2′＝l-l1+l0
l1′和l2′是首尾两个点集最终的长度；
建立包围盒，利用包围盒算法分别计算四个点集(每个断骨模型包含首尾两个点集)的OBB包围盒，并记录长轴、中轴和短轴的方向和长度；
截面包围盒自动选取，求两个包围盒三轴方差的公式为：
l1 l2 l3代表包围盒1的长中短轴的长度；
l1′l2′l3′代表包围盒2的长中短轴的长度；将模型一中的两个包围盒与模型二中的两个包围盒进行两两比较，共计算得到4个方差值，选取方差最小的一组，即为两个断面所在的包围盒，即为提取出来的两个包围盒。
9.根据权利要求1所述的全自动骨折钢板模型个性化重构方法，其特征还在于所述的包围盒对准过程如下：
判断两个断骨模型的贴近一端是否是断面，若不是则进行调整，使断面相互贴近；
提取两个包围盒的三轴在y-o-z平面上的投影方向；
在y-o-z平面上计算长轴和中轴的角分线方向向量
以x轴增大的方向为正方向，将短轴统一为正方向，提取短轴在y-o-z平面上的投影向量；
固定模型一，将模型二以1°为最小单位绕主轴线进行旋转遍历，设旋转的角度为α，则模型二从初始位置进行旋转的旋转矩阵为：
旋转后的向量为
V1′＝AV1
V2′＝AV2
求旋转变换后，两个包围盒对应轴之间的夹角余弦值；
D＝cos a+cos b
对断骨2进行360°的旋转遍历，D的值最大时对应的角度即为包围盒配准时断骨2所需旋转的角度。
10.根据权利要求1所述的全自动骨折钢板模型个性化重构方法，其特征还在于所述对准后，进行点集搜索，得到两个断骨模型截面的点集P1和P2，
对包围盒对准后的模型重新计算面片法向提取得到包含该点的所有面片的原始法向量，其中面片i的原始法向量表示为vi＝(xi,yi,zi)，将所有面片的原始法向量的长度化为1；
面片i的法向量更新为(x′i,y′i,z′i)；
设该数据点处的法向量vp＝(xp,yp,zp)；
n为包含该数据点的面片数量，yp，zp同理；遍历提取出来的所有数据点，计算该点处的法向量与断骨轴线的夹角β，若cosβ>0.35则直接将该点提取至断面点集，否则进行如下步骤；
设置坐标距离阈值r1＝20.0，向量阈值r2＝1.0，对于当前研究的数据点p，在另一个断骨模型中遍历包围盒内数据点p′，同时满足以下三个条件则视该点p′为点p的匹配点：
坐标的平方距离小于阈值，即(x1-x2)2+(y1-y2)2+(z1-z2)2<r1；(xp+xp′)2+(yp+yp′)2+(zp+zp′)2<r2，(xp,yp,zp)(xp′,yp′,zp′)分别是点p和p′处的法向量，点p和p′处的法向量夹角为γ，cosγ>0；
若对点p，能找到至少5个p′满足上述条件，则将p加入至断面点集，(5)分别对断骨模型1和断骨模型2中的包围盒内数据点进行上述遍历，最终提取得到两个断面/截面点集P1和P2。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710542904.0A CN107330281B (zh) | 2017-07-05 | 2017-07-05 | 全自动骨折钢板模型个性化重构方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710542904.0A CN107330281B (zh) | 2017-07-05 | 2017-07-05 | 全自动骨折钢板模型个性化重构方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107330281A CN107330281A (zh) | 2017-11-07 |
| CN107330281B true CN107330281B (zh) | 2020-09-29 |
Family
ID=60196224
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710542904.0A Active CN107330281B (zh) | 2017-07-05 | 2017-07-05 | 全自动骨折钢板模型个性化重构方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107330281B (zh) |
Families Citing this family (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108932723A (zh) * | 2018-03-26 | 2018-12-04 | 天津工业大学 | 一种基于曲面形态的三维Snake主动脉夹层分割方法 |
| CN109035311B (zh) * | 2018-07-11 | 2021-10-01 | 大连理工大学 | 一种弯骨骨折自动配准及内固定钢板预弯建模方法 |
| CN111000626B (zh) * | 2019-10-17 | 2021-06-22 | 中南大学湘雅医院 | 一种钛合金医用内固定板的塑形方法 |
| CN113129456B (zh) * | 2019-12-30 | 2023-07-25 | 百度在线网络技术（北京）有限公司 | 车辆三维模型变形方法、装置和电子设备 |
| CN111383353B (zh) * | 2020-04-01 | 2023-05-23 | 大连理工大学 | 基于高斯混合模型和轮廓描述子的断骨模型配准方法 |
| CN112967374B (zh) * | 2021-02-20 | 2022-12-09 | 济南大学 | 一种骨科手术钢板数字预弯模型获取方法及系统 |
| CN113345112B (zh) * | 2021-05-25 | 2023-06-13 | 上海大学 | 一种长骨骨折断面点云预处理及配准方法 |
| CN115099040B (zh) * | 2022-06-29 | 2023-09-29 | 中国人民解放军91439部队 | 一种板状结构网格模型破口面积计算方法及系统 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105139454A (zh) * | 2015-08-06 | 2015-12-09 | 北京工业大学 | 一种三维ct图像中肝脏三维感兴趣区域的自动提取方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9785858B2 (en) * | 2008-09-26 | 2017-10-10 | Siemens Healthcare Gmbh | Method and system for hierarchical parsing and semantic navigation of full body computed tomography data |
| US10593035B2 (en) * | 2015-03-18 | 2020-03-17 | University Of South Florida | Image-based automated measurement model to predict pelvic organ prolapse |
-
2017
- 2017-07-05 CN CN201710542904.0A patent/CN107330281B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105139454A (zh) * | 2015-08-06 | 2015-12-09 | 北京工业大学 | 一种三维ct图像中肝脏三维感兴趣区域的自动提取方法 |
Non-Patent Citations (2)
| Title |
|---|
| 《Virtual plate pre-bending for the long bone fracture based on axis pre-alignment》;Bin Liu 等;《Computerized Medical Imaging and Graphics》;20140630;第38卷(第4期);全文 * |
| 《基于断骨模型自动配准的完全性骨折钢板预湾》;刘斌 等;《光电子 激光》;20090731;第20卷(第7期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107330281A (zh) | 2017-11-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107330281B (zh) | 全自动骨折钢板模型个性化重构方法 | |
| EP2036037B1 (en) | Methods and systems for segmentation using boundary reparameterization | |
| CN113178009B (zh) | 一种利用点云分割和网格修补的室内三维重建方法 | |
| CN106485695B (zh) | 基于统计形状模型的医学图像Graph Cut分割方法 | |
| CN102592136B (zh) | 基于几何图像中中频信息的三维人脸识别方法 | |
| CN102722890B (zh) | 基于光流场模型的非刚性心脏图像分级配准方法 | |
| Bódis-Szomorú et al. | Superpixel meshes for fast edge-preserving surface reconstruction | |
| CN105279762B (zh) | 一种口腔软硬组织ct序列与三维网格模型配准方法 | |
| CN110688947B (zh) | 一种同步实现人脸三维点云特征点定位和人脸分割的方法 | |
| CN106910242A (zh) | 基于深度相机进行室内完整场景三维重建的方法及系统 | |
| WO2016082797A1 (zh) | 一种基于单幅图像的三维场景结构建模与注册方法 | |
| CN104915965A (zh) | 一种摄像机跟踪方法及装置 | |
| CN113112609A (zh) | 一种面向肺部活检支气管镜的导航方法和系统 | |
| Sinko et al. | 3D registration of the point cloud data using ICP algorithm in medical image analysis | |
| Mayunga et al. | Semi-automatic building extraction utilizing Quickbird imagery | |
| CN108876769B (zh) | 一种左心耳ct图像分割方法 | |
| CN111179321B (zh) | 一种基于模板匹配的点云配准方法 | |
| US20040109594A1 (en) | Method for automatic construction of 2D statistical shape model for the lung regions | |
| CN107316327B (zh) | 一种断骨模型配准方法 | |
| CN114170284B (zh) | 基于主动标志点投射辅助的多视图点云配准方法 | |
| CN109965979A (zh) | 一种无需标志点的稳健的神经导航自动注册方法 | |
| CN109035311B (zh) | 一种弯骨骨折自动配准及内固定钢板预弯建模方法 | |
| CN109064473B (zh) | 一种2.5d超声宽景图像分割方法 | |
| CN109872353B (zh) | 基于改进迭代最近点算法的白光数据与ct数据配准方法 | |
| CN106407902A (zh) | 一种基于几何差异的飞机目标识别方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |