一种点云重构CT图三维复现的方法

一种点云重构ct图三维复现的方法
技术领域
1.本发明涉及医学影像处理技术领域，具体为一种点云重构ct图三维复现的方法。


背景技术：

2.在医学上，ct图作为常用的检测手段被广泛应用，但在长期使用中也显露出了其在直观可视化方面的短板，因此需要采取一些手段来解决该问题，而最有效的方法就是将其复原为立体模型，称为三维复现。


技术实现要素：

3.本发明提供一种点云重构ct图三维复现的方法，解决了ct图在三维重建后精度较差，可视化困难的问题。
4.为实现上述目的，本发明提供如下技术方案：一种点云重构ct图三维复现的方法，包括：
5.s1、获得图像，并采用最大类间方差法对图像进行二值化处理；
6.s2、使用高斯滤波器与图像进行卷积，以平滑图像，减少图像中较为明显的噪声影响，使提取的轮廓更加精确，最后得到降噪处理的图像；
7.s3、得到降噪处理的图像后，根据公式，
[0008][0009]
计算出图像梯度强度g和方向θ,并采用nms继续处理图像，最终完成了图像边缘提取；
[0010]
s4、通过计算得到每张图像之间的z坐标差值δz，以z＝0为初始值，按图像层级顺序为每一组数据额外赋z值，得到一组包含x,y,z位置的新的数据集，利用平面相交分割图像，并采用插值的方式将缺失层补齐，利用散点pchip方法对插值拟合，最终得到散点集；
[0011]
s5、最后根据所述散点集中各个散点的空间坐标，选择以狄洛尼三角网和二维平面的alpha shapes方法为基础，将其结合并拓展到三维平面，构造出本身的形状。
[0012]
优选的，所述s1中，所述最大类间方差法是将灰度图的长宽记为h和w，然后选取一个阈值t，将小于该阈值的灰度归类为前景灰度sum
front
，大于等于该阈值的归类为背景灰度sum
back
，前景灰度和背景灰度的总量应当等于h
×
w，
[0013]
即有式：
[0014]
len
front
+len
back
＝h
×w[0015]
则式：
[0016][0017]
其中m为全局平均灰度，front
p
和back
p
为像素被分为前景灰度和背景灰度的概率，然后计算出前景灰度和背景灰度的平均值front
average
和back
average
带入公式：
[0018]
v＝front
p
*back
p
*(front
average-back
average
)2[0019]
计算出类间方差v，能够使类间方差最大化的灰度k即为所求阈值，再将图像各像素点灰度值和阈值k代入式m＝front
p
×
front
average
+back
p
×
back
average
中，即可完成图像二值化，图像二值化后，得到包含n个二值矩阵的元组m(h,w,n)，其中h,w为图像长宽，n为层数。
[0020]
优选的，所述s3中，所述nms是将采集到的所有可能边缘点与其正负梯度上的其他点进行比较，其中像素梯度强度最大的点保留下来，作为最有可能的边缘点，其余点被抑制舍去，然后便得到了包含n个二值矩阵的二值边缘元组b(x,y,n)，其中x,y为边缘坐标，n为层数，完成点云数据集的初步获取。
[0021]
优选的，所述s4中，所述数据集设为location(x,y,z)，并将其看成n个z平面，其中n为大于0的自然数，首先对其进行切割，用m个垂直于原数据集平面的紧凑x平面与其相交，得到m组散点location_new(x,y,z)，然后使用n组y平面与m组x平面相交，得到散点集location_last(x,y,z)，其中包括m
×
n个散点数据。
[0022]
优选的，所述s4中，在进行所述pchip方法时，根据三次样条插值函数的性质，所述平面x单调增加，选择逐两层拟合，保证三维复现的模型相似度。
[0023]
优选的，所述s5中，所述狄洛尼三角网采用逐点插入lawson算法形成，建立超四面体将点集location_last(x,y,z)包围在四面体内部，然后遍历locatio n_last(x,y,z)，将其中的数据点标号并逐个插入，同时以对角线交换的方式保证所形成的三角网具有最大化最小角特性，使得狄洛尼三角网规则化。
[0024]
优选的，随后将符合条件的三角形作为下一步插点时的边界三角形，继续插入点集中的数据点，循环执行以上步骤，直至所有点全部存在于狄洛尼三角网中，得到一个m
×
4的矩阵delaunay(a,b,c,d)，其中m为四面体个数，a,b,c,d为四面体四个端点的序号。
[0025]
优选的，根据所述矩阵delaunay(a,b,c,d)，提取其中每个四面体的面，保存到n
×
3矩阵triangles(a,b,c)内，其中n为面的个数，a,b,c为形成面三个端点的序号，得到triangles后进行判断，将其中所有重复出现的面舍去，保留仅出现过一次的面，构成这些面的点即为所需的边界点，将这些散点放置在新的点集中，记为location_board(x,y,z)。
[0026]
优选的，设置固定判别半径r，作为后续操作的重要参数之一，根据设置的半径r，遍历location_board(x,y,z)，逐次单独将其中的每两个点提取，并过所述两点绘制半径为r的球体，最终得到的所有线段将形成连接的闭合面，由这些闭合面构成的三维体即为点云重构三维复现所求立体。
[0027]
与现有技术相比，本发明的有益效果：
[0028]
1、本发明中，采用otsu法将图像二值化处理，图像二值化意将非二值图像经过计算变成二值图像，它是进行图像分割最简单的一种方法，即在后续图像处理技术时，可以简化后期的处理，提高处理速度，舍弃一部分无用信息。
[0029]
2、本发明中，使用高斯滤波器与图像进行卷积，以平滑图像，减少图像中较为明显的噪声影响，使提取的轮廓更加精确。
[0030]
3、本发明中，采用nms处理图像，能够获取到单像素边缘，用这些点云进行数据处理可以使得建立的模型更加准确。
[0031]
4、本发明中，相较于其他拟合方法spline(三次样条数据插值)和makima(修正akima分段三次hermite插值)，本发明采用更适合非平滑点的pchip(分段三次hermite插值多项式)方法，该方法在不平滑数据上的表现优于makima和spline，其无过冲且震荡较少，能够较大限度的保留数据的形状和相应单调性。
[0032]
5、本发明中，相较于常用的凸包算法，本发明采用以狄洛尼三角网和二维平面的alpha shapes方法为基础，将其结合并拓展到三维平面，能够有效避免边界线判断错误等问题，优化模型精度。
附图说明
[0033]
附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。
[0034]
在附图中：
[0035]
图1是本发明二值化提取效果图；
[0036]
图2是本发明边缘提取效果图；
[0037]
图3是本发明平面相交切割示意图；
[0038]
图4是本发明平面相交后得到的点云图；
[0039]
图5是本发明插值拟合后局部点云图；
[0040]
图6是本发明五点形成的狄洛尼三角网示意图；
[0041]
图7是本发明提取边界点后的狄洛尼三角网示意图；
[0042]
图8是本发明鼻腔后部三维重现效果图。
具体实施方式
[0043]
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0044]
实施例：一种点云重构ct图三维复现的方法，包括：
[0045]
(1)图像数据处理：
[0046]
为了实现三维点云重构，需要对取得的图像进行数据化，得到合适的参数和数据后，才可以进行三维复现。故第一步采用二值化的方式处理图像。图像二值化意将非二值图像经过计算变成二值图像，它是进行图像分割最简单的一种方法，即在后续图像处理技术时，可以简化后期的处理，提高处理速度，舍弃一部分无用信息。
[0047]
对图像进行二值化处理，得到图像后，对其进行二值化处理，常见的二值化一般算
法为其中，t为二值化阈值，即将大于某个预定值的灰度值设定为1，小于某个预定制的灰度值设定为0，就起到了分离图像信息的作用。二值化阈值的设定分为固定阈值和自适应阈值，对于ct图像，固定阈值可能会导致应用在不同扫描部位上的效果方差较大，导致精度差距变大；
[0048]
本发明采用otsu法，即最大类间方差法来进行图像二值化，如图1所示，首先将灰度图的长宽记为h和w，然后选取一个阈值t，将小于该阈值的灰度归类为前景灰度sum
front
，大于等于该阈值的归类为背景灰度sum
back
，前景灰度和背景灰度的总量应当等于h
×
w，
[0049]
即有式：
[0050]
len
front
+len
back
＝h
×w[0051]
则式：
[0052][0053]
其中m为全局平均灰度，front
p
和back
p
为像素被分为前景灰度和背景灰度的概率，然后计算出前景灰度和背景灰度的平均值front
average
和back
average
带入公式：
[0054]
v＝front
p
*back
p
*(front
average-back
average
)2[0055]
计算出类间方差v，能够使类间方差最大化的灰度k即为所求阈值，再将图像各像素点灰度值和阈值k代入式m＝front
p
×
front
average
+back
p
×
back
average
中，即可完成图像二值化。
[0056]
(2)边缘提取：
[0057]
图像二值化后，得到包含n个二值矩阵的元组m(h,w,n)，其中h,w为图像长宽，n为层数，而在该元组内，只有其中的边缘部分是三维复现时必须的，所以需要将这些边缘提取出来。
[0058]
利用高斯平滑滤波，如图2所示，使用高斯滤波器与图像进行卷积，以平滑图像，减少图像中较为明显的噪声影响，使提取的轮廓更加精确。下式为高斯滤波核的生成公式。
[0059][0060]
若图像中一块3
×
3的区域为a，需要滤波的像素点为e，经过高斯滤波后e点亮度值为
[0061][0062]
其中，*为卷积运算，h为高斯卷积核；
[0063]
得到降噪处理的图像后，根据下式计算图像梯度强度g和方向θ
[0064][0065]
至此，图像边缘的大致提取完成，但目前得到的边缘数据仍不是最精准的，对于三维复现，如果能够获取到单像素边缘，用这些点云进行数据处理可以使得建立的模型更加准确，所以选择采用nms(非极大值抑制)继续处理图像，将采集到的所有可能边缘点与其正负梯度上的其他点进行比较，其中像素梯度强度最大的点保留下来，作为最有可能的边缘点，其余点被抑制舍去，然后便得到了包含n个二值矩阵的二值边缘元组b(x,y,n)，其中x,y为边缘坐标，n为层数，完成了点云数据集的初步获取。
[0066]
(3)函数插值：
[0067]
包含边缘信息的数据对应了三维复现时的x，y位置信息，但这些位置信息并不包含z轴位置信息，那么就需要将这些数据对应还原到空间坐标中。并且由于每张图像z轴差距较大，则需要为其增添数个过度层；
[0068]
利用平面相交分隔图像，如图3所示，根据每一台ct扫描仪的基本参数，再通过计算可以得到每张图像之间的z坐标差值δz，以z＝0为初始值，按图像层级顺序为每一组数据额外赋z值，得到一组包含x,y,z位置的新的数据集location(x,y,z)；但这些数据的每两个z值之间仍有较大空隙，对实施三维复现的精确度有很大的影响，所以考虑采用插值的方式将缺失层补齐；对于已经取得的数据集location(x,y,z)，可以将其看成n个z平面，首先对其进行切割，用m个垂直于原数据集平面的紧凑x平面与其相交，即可得到m组散点locat ion_new(x,y,z)；
[0069]
如图4所示，以其中一组散点为例，在进行拟合之前，由于函数的特性，需要尽可能保证其中任意两点横坐标不相同，否则有可能会导致拟合函数失败，故在保存点时，需要遍历当前组内散点，若有相同点，则将该点移动δx(δx
→
0)后再加入该散点组。
[0070]
利用散点pchip插值拟合，如图5所示，当得到无重复的散点组后，开始着手函数插值，首先需要将这些散点拟合为一个函数，可选的拟合方法有spline(三次样条数据插值)，pchip(分段三次hermite插值多项式)和makima(修正akima分段三次hermite插值)。由于边缘提取后的轮廓散点并不平滑，所以采用更适合非平滑点的pchip方法，该方法在不平滑数
据上的表现优于makima和spline，其无过冲且震荡较少，能够较大限度的保留数据的形状和相应单调性。而根据三次样条插值函数的性质，x应当是严格单调增加，但由图像分解成的边缘散点并无此类限制，故为了防止散点集内出现类似点导致无法拟合的情况，选择逐两层拟合，尽可能的保证三维复现的模型相似度。
[0071]
同时在插值拟合后，对于每个不同的x平面，都有一个，然后使用n组y平面与m组x平面相交，得到散点集location_last(x,y,z)，其中包括m
×
n个散点数据。
[0072]
(4)点云重构：
[0073]
在得到最终的散点数据后，我们需要根据这些散点的空间坐标，构造出本身的形状，也就是点云重构，选择以狄洛尼三角网和二维平面的alpha shapes方法为基础，将其结合并拓展到三维平面，相较于常用的凸包算法，能够有效避免边界线判断错误等问题，优化模型精度。
[0074]
狄洛尼三角剖分作为被广泛使用的散点构面方法，其生成的狄洛尼三角网最外层边界将会形成一个凸边形外壳，在点云重构中有助于将其边界点提取并进行三维面构造。
[0075]
而本文采用逐点插入lawson算法形成狄洛尼三角网，如图6所示，首先建立超四面体t1将点集location_last(x,y,z)包围在四面体t1内部，然后遍历loca tion_last(x,y,z)，将其中的数据点标号并逐个插入。
[0076]
以p0点为例，连接p0和四面体的四个顶点a1，b1，c1，d1，生成六个新三角形对这六个三角形依次进行空外接圆检验，保证其中任意四点不可共圆，即狄洛尼三角网中任意三角形外接圆区域内不存在任何其他点。同时以对角线交换的方式保证所形成的三角网具有最大化最小角特性，使得狄洛尼三角网更接近于规则化。随后将符合条件的三角形作为下一步插点时的边界三角形继续插入点集中的数据点，循环执行上述步骤，直至所有点全部存在于狄洛尼三角网中，得到一个m
×
4的矩阵delaunay(a,b,c,d)，其中m为四面体个数，a,b,c,d为四面体四个端点的序号。
[0077]
对边界点提取，如图7所示，根据矩阵delaunay(a,b,c,d)，提取其中每个四面体的面，保存到n
×
3矩阵triangles(a,b,c)内，其中n为面的个数，a,b,c为形成面三个端点的序号。得到triangles后进行判断，将其中所有重复出现的面舍去，保留仅出现过一次的面，构成这些面的点即为所需的边界点，将这些散点放置在新的点集中，记为location_board(x,y,z)。
[0078]
最后利用alpha shapes方法，如图8所示，设置固定判别半径r，该值将作为后续操作的重要参数之一，其值决定了构成模型的精度，越小模型精度越高；得到半径r后，遍历location_board(x,y,z)，逐次单独将其中的每两个点提取，设提取两点分别为q1,，q2，并过该两点绘制半径为r的球体(排除两点间距离为2r的情况)，若其中任意一个球体中没有除q1，q2外的其他数据点，则连接q1，q2，若两点间距离大于2r，则跳过该两点。
[0079]
最终得到的所有线段将形成连接的闭合面，由这些闭合面构成的三维体即为点云重构三维复现所求立体。
[0080]
最后应说明的是：以上所述仅为本发明的优选实例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保
护范围之内。