CN105205842B - 一种x射线成像系统中的变电流投影融合方法 - Google Patents

Abstract

本发明公开了一种X射线成像系统中的变电流投影融合方法，克服了现有技术中，已有变电流投影融合技术中，人工选择融合阈值的缺陷问题。该发明具体步骤如下：（1）采集成像物体在各投影角度下不同电流的投影数据；（2）构造不同电流投影的灰度值‑有效边缘梯度序列；（3）使用动态时间弯曲方法求解最佳融合阈值；（4）计算缩放因子；（5）投影图像融合；（6）CT图像重建。该发明方法能够避免人工选择融合阈值的主观性，有效扩展探测器的动态范围。探测器动态范围得到扩展后，实验对于探测器的成像需求得到降低，从而减少硬件成本。基于融合后投影的CT图像SNR优于传统人工选择固定阈值的融合方法。

Description

一种X射线成像系统中的变电流投影融合方法

技术领域

该发明涉及X射线成像技术，特别是涉及一种X射线成像系统中的变电流投影融合方法。

背景技术

工业计算机断层成像(Computed Tomography,CT)利用X射线的透射能力对物体成像，通过采集物体在不同角度下的投影数据，基于一定的重建算法获取物体内部的二维或三维CT图像。为了获得优质的成像效果，CT数据采集需要根据物体设置合适的扫描参数。然而，对于一些衰减系数差异明显的多材料物体，传统单一扫描参数的成像模式容易造成探测器曝光不足或曝光过度的现象，由此带来的CT图像伪影或信息缺失极大地影响了CT图像质量，变扫描参数成像方法是解决这类问题的常用方法。

经过对现有技术的文献检索发现，变扫描参数成像方法可以分为变电压和变电流两类。变电压成像方法根据扫描电压决定X射线穿透物体厚度的原理，2009年，杨莹等人在《光子学报》上发表文章“基于小波变换的X射线图像超动态范围重建”，提出了利用小波变换对不同电压下采集的投影进行融合的方法，该方法首先对投影进行多分辨率分解，然后基于不同策略分别计算高、低频系数，最后经过小波逆变换融合出高动态范围的投影。基于小波变换的投影融合方法在求解重构系数时易受伪边缘的影响，融合后的投影可能出现灰阶混乱。为了避免该问题，2013年，Chen等人在《Optik》上发表文章“Multi-voltage imagestack reconstruction in X-ray digital imaging”，提出了一种灰度加权模型，对不同电压下采集的投影序列，首先选择对比度最好的区域作为最佳灰度带，提取待融合的子图像，再通过设置合适的加权系数完成子图像的融合。这种方法没有考虑不同电压的能谱变化，加权系数的求解精度有待提高。2014年，Liu等人在《Journal of X-ray science andtechnology》上发表文章“Multi-energy image sequence fusion based on variableenergy X-ray imaging”，文章所提方法在计算加权系数前，首先通过射束硬化校正保证投影灰度值和物体厚度之间服从线性关系,但是硬化校正要求阶梯状工件和被扫描物体是同一材料的，多材料物体的投影融合问题有待进一步研究。

变电流成像方法基于改变电流只会使光子数量发生变化，不会影响X射线能谱特性的原理，研究该类方法对多材料物体成像具有较高的实用价值。2010年，Kramer等人在《Measurement Science and Technology》上发表文章“Multi-energy image stackfusion in computed tomography”，利用人工设置的融合阈值判断未超出探测器动态范围的投影区域，并计算缩放因子完成投影融合，在对铝和有机玻璃组成的楔形工件成像中获得了更加准确的几何测量结果，但是人工选择投影融合阈值存在主观性。

发明内容

本发明克服了已有变电流投影融合技术中，人工选择融合阈值的缺陷，提供了一种X射线成像系统中的变电流投影融合方法。

本发明的技术解决方案是，提供一种具有以下步骤的X射线成像系统中的变电流投影融合方法：一种X射线成像系统中的变电流投影融合方法，具体步骤如下：(1)采集成像物体在各投影角度下不同电流的投影数据；(2)构造不同电流投影的灰度值-有效边缘梯度序列；(3)使用动态时间弯曲(Dynamic Time Wraping，DTW)方法求解最佳融合阈值；(4)计算缩放因子；(5)投影图像融合；(6)CT图像重建。

所述步骤(1)包括设置合适的扫描电压和电流，采集不同角度下的变电流投影图像，在探测器没有校正的情况下，大、小电流对应的投影分别处于曝光过度和未曝光过度的状态，并且不同电流的投影图像之间存在均未曝光过度的重叠区域，该重叠区域至少为物体投影区域的10％。

所述步骤(2)根据当前投影的灰度值，分别在相同角度下大、小电流的投影中首先确定有效边缘点的位置，然后计算其梯度，每一个灰度值Val都对应了一个有效边缘梯度g_eff，因此可以构成两组灰度值-有效边缘梯度序列。

所述步骤(3)首先将步骤(2)中构造的不同电流序列进行预处理，并构造具有梯度方向的新序列X'和Y'，然后以DTW距离最短为目标计算新序列X'和Y'间的最佳相似匹配路径Q^*，最后根据匹配中疑似突变点附近的子序列X_sub和Y_sub计算相关系数ρ，相关系数最小的点记为最佳融合阈值Th所在的点。

所述步骤(4)根据步骤(3)确定的最佳融合阈值判断大电流投影中未超出探测器动态范围的区域，并计算不同电流投影在相同位置区域的灰度值比例关系,以此作为缩放因子

所述步骤(5)以小电流投影图像为基底，使用步骤(4)计算缩放因子s，将大电流投影中未超出探测器动态范围的区域经过比例缩放后加以替换，融合后的投影

所述步骤(6)针对步骤(5)中融合的投影，使用同一种解析类或迭代类图像重建算法重建CT图像，并以信噪比(Signal to Noise Ratio,SNR)评价融合效果。

与现有技术相比，本发明X射线成像系统中的变电流投影融合方法具有以下优点：该发明方法能够避免人工选择融合阈值的主观性，有效扩展探测器的动态范围。探测器动态范围得到扩展后，实验对于探测器的成像需求得到降低，从而减少硬件成本。基于融合后投影的CT图像SNR优于传统人工选择固定阈值的融合方法。

附图说明

图1是本发明X射线成像系统中的变电流投影融合方法实验过程中实际工业零件位置关系示意图；

图2是本发明X射线成像系统中的变电流投影融合方法的流程图；

图3是本发明X射线成像系统中的变电流投影融合方法在大电流下的实际CT图像；

图4是本发明X射线成像系统中的变电流投影融合方法在不同阈值下融合投影实际CT图像比较：(a)小电流；(b)最小融合阈值；(c)平均融合阈值；(d)本专利方法；

图5是本发明X射线成像系统中的变电流投影融合方法中实际变电流投影融合CT图像的SNR比较图。

具体实施方式

下面结合附图和具体实施方式对本发明X射线成像系统中的变电流投影融合方法作进一步说明：本发明设计一套易于操作且精度较高的变电流投影融合方法。本发明是通过以下技术方案实现的，具体步骤如下：

本发明的技术解决方案是，提供一种具有以下步骤的X射线成像系统中的变电流投影融合方法：具体步骤如下：(1)采集成像物体在各投影角度下不同电流的投影数据；(2)构造不同电流投影的灰度值-有效边缘梯度序列；(3)使用动态时间弯曲(Dynamic TimeWraping，DTW)方法求解最佳融合阈值；(4)计算缩放因子；(5)投影图像融合；(6)CT图像重建。

所述步骤(1)包括设置合适的扫描电压和电流，采集不同角度下的变电流投影图像，在探测器没有校正的情况下，大、小电流对应的投影分别处于曝光过度和未曝光过度的状态，并且不同电流投影的图像之间存在均未曝光过度的重叠区域，该重叠区域至少为物体投影区域的10％。

所述步骤(2)根据当前投影的灰度值，分别在相同角度下大、小电流的投影中首先确定有效边缘点的位置，然后计算其梯度，每一个灰度值Val都对应了一个有效边缘梯度g_eff，因此可以构成两组灰度值-有效边缘梯度序列。

所述步骤(3)首先将步骤(2)中构造的不同电流序列进行预处理，并构造具有梯度方向的新序列X'和Y'，然后以DTW距离最短为目标计算新序列X'和Y'间的最佳相似匹配路径Q^*，最后根据匹配中疑似突变点附近的子序列X_sub和Y_sub计算相关系数ρ，相关系数最小的点记为最佳融合阈值Th所在的点。

所述步骤(4)根据步骤(3)确定的最佳融合阈值判断大电流投影中未超出探测器动态范围的区域，并计算不同电流投影在相同位置区域的灰度值比例关系,以此作为缩放因子

所述步骤(5)以小电流投影图像为基底，使用步骤(4)计算缩放因子s，将大电流投影中未超出探测器动态范围的区域经过比例缩放后加以替换，融合后的投影

所述步骤(6)针对步骤(5)中融合的投影，使用同一种解析类或迭代类图像重建算法重建CT图像，并以信噪比SNR评价融合效果。

步骤(1)采集实际变电流投影数据的具体方法为，首先确定一个可以穿透被成像物体的扫描电压，然后分别在大电流和小电流下采集不同角度下的投影数据。小电流投影应保证未曝光过度的状态，具体判断方法为，在未做校正的情况下，探测器实际最大输出值应小于最大输出的理论值；大电流投影已经达到曝光过度的状态，具体判断方法为，在未做校正的情况下，探测器实际最大输出值等于或接近最大输出的理论值，部分物体区域的投影灰度值等于背景的灰度值。此外，大电流投影和小电流投影之间应该保证一定的均未曝光过度的重叠区域，该重叠区域至少为物体投影区域的10％；

步骤(2)构造不同电流投影的灰度值-有效边缘梯度序列的具体方法为，如果探测器上相邻两个像素的灰度值分别大于和小于Val(Val介于0和当前投影的最大灰度值之间)，那么就将这组点对记为灰度值Val下的有效边缘。假设投影中坐标(u,v)的灰度值为p(u,v)，其梯度向量g_edge(u,v)可表示为g_edge(u,v)＝(g_up(u,v),g_down(u,v),g_left(u,v),g_right(u,v))^T，公式(1)用于计算上下左右4个方向的梯度

\*MERGEFORMAT(1)

全部像素点的梯度矩阵记为G＝(g_edge(1,1),g_edge(1,2),...,g_edge(U,V))，其中U和V分别表示投影的高度和宽度，整幅投影的有效边缘梯度g_eff定义为

\*MERGEFORMAT(2)

对某个角度下两幅不同电流的投影，采样从0到当前投影的最大灰度值之间的不同灰度值Val，每个灰度值Val都有一个有效边缘梯度g_eff与之对应，由此可以得到两幅不同电流投影的两组灰度值-有效边缘梯度序列。

步骤(3)求解最佳融合阈值的具体方法为，为了解决电流大小不同导致步骤(2)构造的灰度值-有效边缘梯度序列的灰度值和有效边缘梯度在幅度上的伸缩问题，需要对原始的灰度值-有效边缘梯度序列进行预处理。首先根据采样得到的两种电流序列的有效边缘梯度峰值点位置关系，选择小电流序列中与大电流序列形状相似的前一段；然后对新选择的小电流序列进行重新采样，采样步长可以根据实际需求确定；最后分别针对灰度值和有效边缘梯度，对原始的大电流序列和新选择的小电流序列做归一化处理，得到大、小电流下经过归一化处理后的灰度值-有效边缘梯度序列。两个序列的灰度值范围均为[0，1]，若对不同电流投影均以当前投影最大灰度值相同的百分比为步长进行采样，则可忽略因灰度值差异带来的伸缩，只考虑有效边缘梯度的差异。将大、小电流下长度分别为n和m的灰度值-有效边缘梯度序列记为X和Y(若DTW相似匹配过程中，序列X和Y后方补0数量相等，则有n＝m)，序列中任意两点x_i和y_j表示对不同电流投影，采样归一化处理后的灰度值，对应的有效边缘梯度g_eff。记x_i和y_j之间的对应关系矩阵为M，M_ij表示x_i和y_j的距离

d(x_i,y_j)＝|x_i-y_j|.\*MERGEFORMAT(3)

两个序列之间的相似匹配关系用路径Q＝(q₁,q₂,...,q_K)表示，路径的最后一个点K∈[max(m,n),m+n-1]，路径的第k个元素q_k＝(i,j)表示序列X中第i个元素x_i和序列Y中第j个元素y_j互相匹配。将两点间的距离d(x_i,y_j)简记为d_k，最佳路径Q^*具有最小弯曲代价

\*MERGEFORMAT(4)

路径Q^*可以利用基于累积增量矩阵的动态规划方法计算，累积增量矩阵D定义为

\*MERGEFORMAT(5)

在获得累积增量矩阵D后，确定最优弯曲路径的具体方法为，从弯曲路径的终点开始依次回溯到起点。对于累积增量矩阵D中的某个元素D(i,j)，分别比较三个相邻点D(i-1,j)，D(i,j-1)和D(i-1,j-1)，选择其中的最小点作为下一步回溯点。

在使用DTW方法对序列进行相似匹配时可能存在以下问题，如果采样的间隔太小，那么在每个采样间隔内的点数很少，计算的投影有效边缘梯度平均效果不明显，结果容易受到噪声影响。为此，本专利在原始距离的基础上引入序列的增量方向，构造了一种带有梯度方向的灰度值-有效边缘梯度序列。大电流下的新序列可以通过公式(6)计算

\*MERGEFORMAT(6)

同理也可以得到小电流下的新序列

考虑到噪声等因素的影响，在大电流有效梯度序列和小电流有效梯度序列的相似匹配点对中，可能出现一个点对应多个点或是多个点对应一个点的情况，将所有符合这一条件的点记为探测器超动态范围的疑似突变点，并且使用公式(7)计算疑似突变点附近长度为L_new的子序列X_sub和Y_sub的相关系数ρ，和表示子序列的均值，子序列的长度取决于疑似突变点的对应情况，之后根据互相关系数的差异判断真正的变电流投影融合阈值，相关系数ρ最小的点记为最佳融合阈值Th所在的点

\*MERGEFORMAT(7)

步骤(4)计算缩放因子的具体方法为，假设步骤(3)求出的最佳融合阈值为Th，p₁和p₂分别表示同一角度下大、小电流的两幅投影，缩放因子s可通过一定条件下不同电流投影中相同区域的灰度值求解

\*MERGEFORMAT(8)

步骤(5)投影图像融合的具体方法为，使用步骤(4)求到的缩放因子s，融合后的投影p_fusion(u,v)可以用公式(9)计算

\*MERGEFORMAT(9)

步骤(6)图像重建的具体方法为，本发明专利对重建算法并没有明确限制，可以是解析类算法也可以是迭代类的重建算法，后续如何由投影数据获得三维CT图像为本领域公知，在此处不再赘述。需要强调的是，本发明专利所述方法在比较融合前后的CT图像质量时，均针对同一种重建算法，评价指标SNR用公式(10)计算

\*MERGEFORMAT(10)

其中，Ω_ROI表示CT图像中的均匀区域。

具体实验过程：为了进一步验证所提投影融合方法的有效性，本发明针对实际工业零件开展了实验。该工业零件由7个圆柱组成，外层圆柱为低衰减的铝材料，内层圆柱为高衰减的钢材料。其内部结构的相对位置关系如附图1所示，圆柱直径的标称值为6mm。

实验中使用的CT系统包括了德国YXLon公司生产的FXE 225.48射线源和美国Varian公司生产的4030E探测器。扫描电压为180kV，电流包括100μA和300μA两种，射线源到旋转轴距离为177.66mm，射线源到探测器距离为1033.68mm，采集投影时在360°范围内每隔0.5°采集一张投影，共720张投影，图像重建选择了解析类的FDK算法，重建体素尺寸为0.022mm³。

按照本专利所提附图2的投影融合流程，构造灰度值-有效边缘梯度序列时采样步长均选为当前投影最大灰度值的0.5％，本专利方法可以自动判断不同角度下变电流投影融合的最佳阈值。在300μA的大电流条件下探测器已经曝光过度，重建CT图像如附图3所示，重建CT图像中外层低衰减的铝材料信息缺失。将本专利方法的融合效果与未融合的小电流CT图像，以及传统人工选择固定融合阈值的方法作比较。在人工固定融合阈值的选择上，使用了本专利方法计算的不同角度下融合阈值的最小值和均值。两组对比实验各角度均采用固定不变的融合阈值，实际数据受到探测器本底噪声、响应不一致、散射等因素的影响，不同角度投影使用相同融合阈值带来的影响十分明显。从附图4的比较中可以看到，一方面，附图4(a)的小电流CT图像均匀性更差；另一方面，不恰当的融合阈值导致投影融合时缩放因子计算错误，进而导致附图4(b)和附图4(c)固定阈值融合后的CT图像在融合边界处出现跳跃。使用本专利所提融合方法得到的如附图4(d)所示的融合结果，既能使CT图像均匀性得到改善，又能避免融合后的CT图像在融合边界处出现跳跃的现象。

在比较不同材料圆柱区域的CT图像SNR时，选择了各圆柱区域中心附近170×170个像素的正方形区域。从附图5的比较中可以发现，融合后的CT图像SNR明显优于融合前的小电流CT图像，并且本专利方法优于两种固定融合阈值CT图像的SNR。

Claims (4)

1.一种X射线成像系统中的变电流投影融合方法，其特征在于：具体步骤如下：

(1)采集成像物体在各投影角度下不同电流的投影数据；设置合适的扫描电压和电流，采集不同角度下变电流的投影图像，在探测器没有校正的情况下，大、小电流对应的投影分别处于曝光过度和未曝光过度的状态，并且不同电流的投影图像之间存在均未曝光过度的重叠区域，该重叠区域至少为物体投影区域的10％；

(2)构造不同电流投影的灰度值-有效边缘梯度序列；根据当前投影的灰度值，分别在相同角度下大、小电流的投影中首先确定有效边缘点的位置，然后计算有效边缘点梯度，每一个灰度值Val都对应了一个有效边缘梯度g_eff，因此构成两组灰度值-有效边缘梯度序列；

(3)使用动态时间弯曲方法求解最佳融合阈值；首先将步骤(2)中构造的不同电流投影的灰度值-有效边缘梯度序列进行预处理，并构造具有梯度方向的新序列X'和Y'，然后以动态时间弯曲距离最短为目标计算新序列X'和Y'间的最佳相似匹配路径Q^*，最后根据匹配中疑似突变点附近的子序列X_sub和Y_sub计算相关系数ρ，相关系数最小的点记为最佳融合阈值Th所在的点；

(4)计算缩放因子；

(5)投影图像融合；

(6)CT图像重建。

2.根据权利要求1所述的X射线成像系统中的变电流投影融合方法，其特征在于：所述步骤(4)根据步骤(3)确定的最佳融合阈值判断大电流投影中未超出探测器动态范围的区域，并计算不同电流投影在相同位置区域的灰度值比例关系,以此作为缩放因子其中p₁(u,v)和p₂(u,v)分别表示坐标位置(u,v)处大、小电流投影的灰度值,Th为最佳融合阈值。

3.根据权利要求1所述的X射线成像系统中的变电流投影融合方法，其特征在于：所述步骤(5)以小电流投影图像为基底，使用步骤(4)计算缩放因子s，将大电流投影中未超出探测器动态范围的区域经过比例缩放后加以替换，融合后的投影其中p₁(u,v)和p₂(u,v)分别表示坐标位置(u,v)处大、小电流投影的灰度值。

4.根据权利要求1所述的X射线成像系统中的变电流投影融合方法，其特征在于：所述步骤(6)针对步骤(5)中融合的投影，使用同一种解析类或迭代类图像重建算法重建CT图像，并以信噪比评价融合效果。