基于锥束CT的工业部件三维图像表面平滑方法与流程

本发明属于图像处理技术领域，涉及一种基于锥束CT的工业部件三维图像表面平滑方法。

背景技术：

工业部件的三维图像在获取和传输的过程中会受到各种各样噪声的干扰，而且三维图像的表面必须通过分割或者拟合的方法来得到，这都会造成一些伪影，比如洞或者不连续的部分。这会严重影响三维图像的视觉效果，也会给后续的图像分析和理解带来一定的困难。为了抑制噪声、改善图像质量，需要对图像进行表面平滑处理，而平滑质量的好坏将会直接影响后续的处理。所以，三维图像的表面平滑技术的研究具有较大的实际意义。

现有技术中，基于空域和变换域的图像平滑方法算法简单，计算复杂度不高，执行速度快，但是它们也有很多不足，最明显的一点就是对图像的纹理和细节保持不够理想，它们在去噪的同时模糊了边缘、纹理和细节等重要的信息。基于偏微分方程的图像平滑方法去噪效果明显，而且对边缘细节的保护比较好，但是因为它们平滑的前提是假设图像是分片常数的，所以它们在图像比较光滑的部分会产生阶梯效应。2005年Buades提出非局部均值算法，它把图像中出现很多周期性的冗余信息加以利用，在全局中搜索灰度相似块并充分利用像素之间的相关性，它不仅能够去除噪声，还能更好的保护边缘特征。2008年董斌提出基于水平集的非局部表面恢复算法，它利用隐式的水平集方法来表示表面，这样不仅数值简单，还能让拓扑结构变化灵活多变，它能使三维图像的表面在边界得到保护的同时减少了大量的噪声，而且表面的有些地方还得到了重造，让整个表面更加的平滑。

技术实现要素：

有鉴于此，本发明的目的在于提供一种基于锥束CT的工业部件三维图像表面平滑方法，该方法在扫描前将射线源和面阵探测器按照一定的放大比分别放置在待扫描工业部件的两侧，扫描开始后，将射线源对准面阵探测器的中心，射线源和探测器绕工业部件的中心旋转一周来获得完整的投影数据，然后将获得的投影数据传送给控制与图像处理系统并利用FDK重建算法进行三维图像重建；得到工业部件的三维图像后使用3D-RCV算法来获得它的表面，最后利用基于水平集的非局部表面恢复算法来对工业部件三维图像的表面进行平滑。经过该方法处理后，工业部件三维图像的表面在边界得到保护的同时减少了大量的噪声，而且在表面的有些地方还进行了重造，让表面更加的平滑。

为达到上述目的，本发明提供如下技术方案：

一种基于锥束CT的工业部件三维图像表面平滑方法，该方法包括以下步骤：

S1：检测装置安装：该检测装置包括射线源(1)、面阵探测器(2)以及控制及图像处理系统(5)，所述射线源(1)、面阵探测器(2)的信号线路与控制与图像处理系统(5)相连，射线源的圆形轨道(3)的中心与工业部件(4)某一横截面的圆心重合，射线源(1)、面阵探测器(2)分别固定在各自的圆形轨道上，使得射线源(1)产生的锥形射线束能够覆盖工业部件(4)的全体区域；

S2：射线源和探测器绕工业部件的中心旋转一周来获得完整的投影数据：在控制与图像处理系统(5)的控制下，首先将射线源(1)对准面阵探测器(2)的中心，射线源(1)和面阵探测器(2)绕工业部件的中心旋转一周来获得完整的投影数据，然后传送到控制与图像处理系统(5)中存储；

S3：工业部件的三维图像重建：根据得到的投影数据重建工业部件的三维图像；

S4：利用3D-RCV算法来得到工业部件三维图像的表面；

S5：使用基于水平集的非局部表面恢复算法来对工业部件的三维图像表面进行平滑；

S6：显示三维图像表面的平滑结果。

进一步，在步骤S3中，所述的根据得到的投影数据重建工业部件的三维图像具体包括：S31：对投影数据进行加权；S32：结合斜坡滤波器对加权后的投影数据进行一维滤波；S33：对步骤S32的滤波结果进行三维加权反投影；

具体算法如下：

1)面阵探测器实际采集到的投影数据为p(β,a,b)，则投影数据加权后为：

其中，β表示中心射线与y轴的夹角，a,b分别代表面阵探测器上的水平位置和竖直位置坐标，R表示射线源到旋转中心的距离；

2)对加权后的投影数据p′(β,a,b)进行一维滤波

其中“*”代表卷积算子，q为积分变量；

3)将卷积后的投影数据转化为虚拟探测器上的投影数据则待重建点(x,y,z)处的重建公式为：

其中(x,y,z)表示待重建点的三维坐标，f(x,y,z)表示在(x,y,z)处的图像灰度，U(x,y,β)＝R+xcosβ+ysinβ，表示对应射线在虚拟探测器上的水平位置坐标，表示对应射线在虚拟探测器上的竖直位置坐标。

进一步，在步骤S4中，利用3D-RCV算法来得到工业部件三维图像的表面通过求解以下偏微分方程得到：

其中，I(y)是点x的邻域N(x)内点q(x,y,z)的灰度值，x,y,z分别代表笛卡尔坐标系中的点对应的三个坐标分量的值，为水平集函数，为定义的初始轮廓，为函数的梯度，H_ε为Heaviside函数的正则化形式，δ_ε为一维Dirac测度的正则化形式，ε为正常数，div(·)表示求散度，c_i表示图像被划分为目标即轮廓面内部和背景即轮廓面外部两个同质区域时的平均灰度，μ,v≥0，λ₁,λ₂＞0是各项的权值系数，t为引入的人工变量，m是正则化常数，σ是高斯函数的标准差，r是邻域N(x)的半径。

进一步，在步骤S5中，所述使用基于水平集的非局部表面恢复算法来对工业部件的三维图像表面进行平滑具体包括以下步骤：

S51：采用符号距离函数φ来表示三维图像的表面：

其中，t表示图像的像素点，它的坐标为(x,y,z),Σ表示物体边界的内部，S表示区域Σ的边界；

S52：在φ的0水平集的一个狭窄的带上进行计算，把狭窄带记作Σ_η，η是狭窄带的宽度；

S53：选择权重ω(x,y)和相似函数D(x,y)：

其中：b₁和b₂是两个正则化参数，N_η是x在Σ_η中的邻域，φ[x]是φ中心在x处的3D块；

S54：对能量函数J(u)的梯度下降流使用有限差分法得到离散的迭代格式如下：

ω_jl可以由S53的式子计算得出，dt令为：

j遍历完Σ_η中所有的网格点；

S55：给定停止迭代步数k。

本发明的有益效果在于：本发明提供的方法能够很好的克服现有技术中存在的问题，使得工业部件三维图像的表面在边界得到保护的同时减少了大量的噪声，而且在表面的有些地方还进行了重造，让表面更加的平滑。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明的工业部件扫描结构示意图；

图2为本发明的圆周轨迹锥束CT的FDK重建算法几何结构示意图；

图3为本发明的基于锥束CT的工业部件三维图像表面平滑方法流程图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明的工业部件扫描结构示意图，如图所示，将射线源(1)、面阵探测器(2)分别固定在各自的圆形轨道上，以起始射线源(1)到工业部件(4)中心轴的垂足为坐标原点O建立空间直角坐标系O-xyz，x轴为原点与射线源(1)的连线并且正方向为从原点指向射线源(1)，y轴为沿工业部件(4)横向并垂直于x轴的坐标轴(图2)，z轴为与工业部件(4)中心轴线重合的坐标轴并且以工业部件一侧作为正方向(图1)。以坐标原点O为旋转中心，(x,y,z)表示被重建点坐标，在过坐标原点O并且垂直于中心射线的位置处引入虚拟探测器，射线源(1)位于S，SO′表示锥形射线的中心射线，SK表示经过被重建点的一条射线，K′表示K点在探测器中心层(z＝0)的投影，被重建点在直线SK′上的投影为M，中心射线与y轴成β角，κ为射线SK的锥角(图2为本发明的圆周轨迹锥束CT的FDK重建算法几何结构示意图)。

图3为本发明的基于锥束CT的工业部件三维图像表面平滑方法流程图，如图所示，基于锥束CT的工业部件三维图像表面平滑方法包括以下步骤：

S1.检测装置安装：锥束CT扫描装置，包括射线源(1)、面阵探测器(2)以及控制及图像处理系统(5)，所述射线源(1)、面阵探测器(2)的信号线路与控制及图像处理系统(5)相连，射线源的圆形轨道(3)的中心与工业部件(4)某一横截面的圆心重合，射线源(1)、面阵探测器(2)分别固定在各自的圆形轨道上，使得射线源(1)产生的锥形射线束能够覆盖工业部件(4)的全体区域；

S2.扫描：在控制与图像处理系统(5)的控制下，首先将射线源(1)对准面阵探测器(2)的中心，射线源(1)和面阵探测器(2)绕工业部件的中心旋转一周来获得完整的投影数据，然后传送到控制与图像处理系统(5)中存储；

S3.工业部件的三维图像的重建：根据得到的投影数据重建工业部件的三维图像，主要包含三个步骤：S31.对投影数据进行加权；S32.结合斜坡滤波器对加权后的投影数据进行一维滤波；S33.对步骤S32的滤波结果进行三维加权反投影。

本发明采用的重建算法如下：

(1)面阵探测器实际采集到的投影数据为p(β,a,b)，则投影数据加权后为：

其中，β表示中心射线与y轴的夹角，a,b分别代表面阵探测器上的水平位置和竖直位置坐标，R表示射线源到旋转中心的距离。

(2)对加权后的投影数据p′(β,a,b)进行一维滤波

其中“*”代表卷积算子，q为积分变量。

(3)将卷积后的投影数据转化为虚拟探测器上的投影数据则待重建点(x,y,z)处的重建公式为：

其中(x,y,z)表示待重建点的三维坐标，f(x,y,z)表示在(x,y,z)处的图像灰度，U(x,y,β)＝R+xcosβ+ysinβ，表示对应射线在虚拟探测器上的水平位置坐标，表示对应射线在虚拟探测器上的竖直位置坐标。

S4.利用3D-RCV算法来得到工业部件三维图像的表面；

通过求解以下偏微分方程得到：

其中，I(y)是点x的邻域N(x)内点q(x,y,z)的灰度值，x,y,z分别代表笛卡尔坐标系中的点对应的三个坐标分量的值，为水平集函数，为定义的初始轮廓，为函数的梯度，H_ε为Heaviside函数的正则化形式，δ_ε为一维Dirac测度的正则化形式，ε为正常数，div(·)表示求散度，c_i表示图像被划分为目标即轮廓面内部和背景即轮廓面外部两个同质区域时的平均灰度，μ,v≥0，λ₁,λ₂＞0是各项的权值系数，t为引入的人工变量，m是正则化常数，σ是高斯函数的标准差，r是邻域N(x)的半径；

S5.使用基于水平集的非局部表面恢复算法来对表面进行平滑。

本发明采用的平滑算法如下：

S51：采用符号距离函数φ来表示三维图像的表面：

其中，t表示图像的像素点，它的坐标为(x,y,z),Σ表示物体边界的内部，S表示区域Σ的边界；

S52：在φ的0水平集的一个狭窄的带上进行计算，把狭窄带记作∑_η，η是狭窄带的宽度；

S53：选择权重ω(x,y)和相似函数D(x,y)：

其中：b₁和b₂是两个正则化参数，N_η是x在∑_η中的邻域，φ[x]是φ中心在x处的3D块；

S54：对能量函数J(u)的梯度下降流使用有限差分法得到离散的迭代格式如下：

ω_jl可以由S53的式子计算得出，dt令为：

j遍历完∑_η中所有的网格点；

S55：给定停止迭代步数k。

S6.显示三维图像表面的平滑结果。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其做出各种各样的改变，而不偏离本发明权利要求书所限定的范围。