CN109813259B

CN109813259B - 高动态x射线成像方法、存储介质和装置

Info

Publication number: CN109813259B
Application number: CN201910048057.1A
Authority: CN
Inventors: 郭鹏程; 陈平; 杨硕; 魏交统; 刘宾; 苏新彦; 王鉴
Original assignee: North University of China
Current assignee: North University of China
Priority date: 2019-01-18
Filing date: 2019-01-18
Publication date: 2021-06-15
Anticipated expiration: 2039-01-18
Also published as: CN109813259A

Abstract

本发明提供一种高动态X射线成像方法、存储介质和装置，该方法包括：步骤11：采用预设低能量进行全角度扫描，获得第一图像集；步骤12：提取每个第一图像的工件投影区域和边缘轮廓信息；计算每个第一图像的灰度值，将最大或最小灰度值对应的扫描角度记为θ1；步骤14：固定θ1角度逐步调整X射线能量对被测物进行成像，获得第二图像集，定位目标第二图像，将目标第二图像的平均灰度值记为H1、管电压记为V1、管电流记为I1；步骤15：从θ1开始全角度扫描，获得第三图像集；步骤16：基于非线性扩散方程将第三图像集恢复为完整的CT二维图像集。基于本发明的方法，相比现有的递变能量CT成像技术，能极大减小曝光次数，延长射线管的使用寿命。

Description

高动态X射线成像方法、存储介质和装置

技术领域

本发明涉及仪器领域，特别涉及一种高动态X射线成像方法、存储介质和装置。

背景技术

X射线CT成像是通过探测器获取每个投影角度上的射线携带工件的衰减信息，通过CT重建获取工件的3D结构。为了保证得到较高的图像质量，要求成像系统的射线能量、探测器的动态范围和工件等效衰减厚度三者匹配。但是对于一些复杂结构件(如：外形结构复杂的支撑件、内部空洞较多的铸造件、内部结构复杂的仪表等)，在CT旋转角度上，等效厚度差异较大，结构越复杂，差异越明显。如此，按照传统的固定能量的成像方式，势必在某些投影角度上，能量、动态范围和等效厚度，三者不匹配，极易出现厚度厚的地方剂量不足，噪声水平高；薄的地方，射线能量过高，散射现象严重，信息缺失，CT重建质量差，无法进行相关的CT检测工作。

对于复杂异形工件的CT成像，早期基于低动态成像系统，提出了递变能量CT成像技术，通过获取每个投影角度下递变能量投影序列，并通过融合扩大成像系统动态范围，有效克服了因成像探测器动态范围不足而引起的投影信息大量缺失的问题。随着更高动态范围的探测器出现，容纳的复杂结构件的等效厚度差异性得到极大地提高。但是采用每个投影角度下的递变能量成像，通过融合实现高动态CT，由于每个投影角度下采集投影序列过多，一方面影响工件内部投影灰度与射线能量的对应性，另一方曝光次数多，影响射线管使用寿命。

发明内容

有鉴于此，本发明提供一种高动态X射线成像方法、存储介质和装置，以解决递变能量CT成像技术存在曝光过度的问题。

本发明提供一种高动态X射线成像方法，包括：

步骤11：采用预设低能量的X射线对被测物进行预设角度间隔△θ的全角度扫描，获得第一图像集；

步骤12：提取每个第一图像的工件投影区域和边缘轮廓信息；

步骤13：计算每个第一图像工件投影区域的像素灰度值，将最大像素灰度值或最小像素灰度值对应的第一图像的扫描角度记为θ1；

步骤14：固定θ1角度对被测物进行X射线能量递变成像，获得第二图像集，选择第二图像集中清晰度指标最大值对应的图像为目标第二图像，将目标第二图像中的平均灰度值记为H1，目标第二图像对应的X射线能量的管电压记为V1、管电流记为I1；

步骤15：从θ1开始对被测物进行△θ间隔的全角度扫描，获得第三图像集；每个角度成像时X射线能量的管电流为I1，通过调整当前管电压V，使得当前角度对应的第三图像的平均灰度值H与H1的偏差值满足预设条件；

步骤16：基于非线性扩散方程恢复第三图像集中因散射而损失掉的结构信息，获得内部结构完整的CT二维图像集，非线性扩散方程的初始条件为边缘轮廓信息。

本发明还提供一种非瞬时计算机可读存储介质，非瞬时计算机可读存储介质存储指令，指令在由处理器执行时使得处理器执行上述的高动态X射线成像方法中的步骤。

本发明还提供一种高动态X射线成像装置，包括存储器、处理器以及存储在存储器中并可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述的高动态X射线成像方法中的步骤。

本发明的高动态X射线成像方法，通过步骤15灰度变化自适应调整每个角度下的X射线能量的管电压V，相比现有的递变能量CT成像技术能极大减小曝光次数，延长射线管的使用寿命。

此外，本发明还通过定位视角θ1(最小厚度视角或最大厚度视角)确定起始能量(V1和I1)，减小步骤15的能量调整难度和步骤16非线性扩散方程的计算难度，提高本发明方法的执行效率和CT图像的构建精度。

附图说明

图1为本发明高动态X射线成像方法的流程图；

图2为本发明高动态X射线成像装置的结构图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明进行详细描述。

如图1所示，本发明的高动态X射线成像方法，包括：

步骤11(S11)：采用预设低能量的X射线对被测物进行预设角度间隔△θ的全角度扫描，获得第一图像集；

预设低能量的获取的第一图像要求：可明显分辨工件投影区域和背景区域。

当CT的X射线工作管电压范围是：0～450kV，管电流为2mA，考虑到图像质量要求，预设低能量可设为110VA至130VA中的任一值，例如，管电压设为60kV，管电流设为2mA。

步骤12(S12)：提取每个第一图像的工件投影区域和边缘轮廓信息；

步骤11输出的第一图像，背景区域和工件投影区域具有灰度差异，工件投影边界清晰，通过对第一图像进行二值化处理，可提取或标记每个第一图像的工件投影区域和边缘轮廓信息。

步骤13(S13)：计算每个第一图像工件投影区域的像素灰度值，将最大像素灰度值或最小像素灰度值对应的第一图像的扫描角度记为θ1；

在完成轮廓提取的同时，对每个第一图像(不同扫描角度)工件投影区域的像素灰度值进行计算统计。

依据X射线成像原理，放大比固定且相同透照能量下，等效厚度越薄，穿透工件到达探测器的光子数越多，对应的图像灰度越高。计算各投影角度下工件投影区域的像素灰度值H(i,j)，i和j为像素坐标，将H(i,j)中的最大值或最小值记为H_m，如下式：

或

w1与h1分别表示工件投影区域像素宽度与像素高度。

当H_m为像素灰度最大值时，θ1为最小厚度视角，当H_m为像素灰度最小值时，θ1为最大厚度视角，最大厚度视角和最小厚度视角均可作为步骤15的起始视角，确保步骤15全角度扫描时，管电压单向变化满足预设条件，单向递增或单向递减，降低步骤15中管电压自适应调整难度。

步骤14(S14)：固定θ1角度对被测物进行X射线能量递变成像，获得第二图像集，选择第二图像集中清晰度指标最大值对应的图像为目标第二图像，将目标第二图像中的平均灰度值记为H1，目标第二图像对应的X射线能量的管电压记为V1、管电流记为I1；

步骤14是为了选出θ1角度下最清晰的成像作为目标第二图像(或称为标准图像)，其中最清晰的筛选标准是：清晰度指标，清晰度指标可参考现有算法计算，例如：Brenner梯度函数、Tenengrad梯度函数、Laplacian梯度函数、SMD(灰度方差)函数等。为了便于计算，清晰度指标也可以设定为：图像的平均灰度值或图像中工件投影区域的平均灰度值。

式2给出了平均灰度值H1的计算方法：

或

H1即可以设定为整幅图像的平均灰度，也可以设定为一幅图像中工件投影区域的平均灰度；式中，w2与h2分别表示整幅图像的像素宽度与像素高度，w1与h1分别表示工件投影区域的像素宽度与像素高度。

步骤15(S15)：从θ1开始对被测物进行△θ间隔的全角度扫描，获得第三图像集；每个角度成像时X射线能量的管电流为I1，通过调整当前管电压V，使得当前角度对应的第三图像的平均灰度值H与H1的偏差值满足预设条件；

步骤15的自适应变能量扫描与步骤11的低能量预扫描采取相同的采样间隔△θ，确保各角度高低能量投影图像完全匹配。

依据步骤14，设定最小厚度视角或最大厚度视角为步骤15的初始角度θ1，V1为初始电压。每个角度成像可包含以下步骤：

步骤151：令当前管电压V＝V1、管电流为I1，获取当前角度的第三图像，计算第三图像的平均灰度值H；

步骤152：判断H与H1的偏差是否满足预设条件，如果是，则调整当前角度值，返回步骤151，如果否，则执行步骤153；

步骤153：根据H与H1的比值，调整当前管电压V，并更新当前角度的第三图像和对应的H，返回步骤152。

H与H1的定义相同，H和H1为整幅图像的平均灰度，或H和H1为一幅图像中工件投影区域的平均灰度。

步骤151-153是基于成像灰度-电压关系，每个角度成像时，可先采用V1成像获得第三图像，通过第三图像的平均灰度H与H1的差异，计算该角度下对应的电压量V，再采用V成像更新第三图像，并检测更新后的第三图像的H是否符合预设条件，符合则进行下一个角度成像，否则继续调整V，直至全角度扫描结束。通过H调整V的原理说明如下：

X射线的总强度T可以由以下参数来表示:管电流I(mA)，管电压V(kV)及靶材料原子序数Z，常数K_i＝1.1～1.4×10^-6：

T＝K_iZIV² (3)

而H与T成线性关系，

H＝αT＝αK_iZIV² (4)

其中α为比例系数，一般取8.2，根据公式(4)可得，θ1下管电压V1与图像平均灰度H1间的关系为

当扫描成像的过程中电流恒定不变，则X射线能量只与管电压V有关。使用V1对下一角度成像得到即时平均灰度H，依据其与H1的差异调整电压V，为方便调整，采用下取整舍去误差。

即

为确保能量自适应成像信息的有效性，会对调整电压V后所得图像的H与H1比较，若变化量的绝对值超过H1的预设值(例如6％)，则视为电压调节不到位，会继续调整，直至H达到标准才会进去下一扫描角度成像。

预设条件的设置可确保每个第三图像均具有较高图像质量，最大限度减小步骤16所需的恢复信息，并确保CT二维图像的质量。

步骤16(S16)：基于非线性扩散方程恢复第三图像集中因散射而损失掉的结构信息，获得内部结构完整的CT二维图像集，非线性扩散方程的初始条件为边缘轮廓信息。

待步骤15扫描采集结束，依据预扫描过程中获取每个角度的边缘轮廓信息，并基于非线性扩散方程建立投影散射恢复的数学模型，将步骤12获得的边缘轮廓信息作为方程初始条件，对于步骤15获得的每个角度的第三图像(高能投影图像)中因散射而损失掉的结构信息进行恢复，即对非线性扩散方程加以约束求得该扩散方程的唯一解，得到待恢复区域的图像灰度，最后获得轮廓细节清晰、内部结构完整的CT二维图像。

常用的非线性扩散方程为热传导方程，经典的热传导方程初值问题为

其中div为散度算子，

为梯度算子，初始条件为零时刻的u。将u替换为H，添加边缘函数

使热传导方程具有了各向异性的扩散方式，即PM方程：

其中

表示图像灰度值的梯度，

作为图像的幅值表示了图像灰度值的变化幅度，扩散函数g取单调递减函数。采用有限差分法求解方程，加入第一图像的边缘轮廓信息作为初始条件H₀，t作为尺度参数，为求最终投影恢复图像H_t(CT二维图像)，将PM方程离散变形可得恢复模型：

即

方程离散化处理，对于图像上各个像素(i,j)求偏导得：

其中

Δ(H₀)_i,j＝(H_xx)_i,j+(H_yy)_i,j＝((H_x)_x)_i,j+((H_y)_y)_i,j (13)

为保持数值的稳定性，且假设图像水平与竖直方向上的差分步长一致，采用一种特殊的差分格式

经过离散运算可求得H_t，即最终得到的投影散射恢复图像(CT二维图像)。

在步骤16之后还包括：

步骤17：基于CT二维图像集构建被测物的CT三维图像。

采用常用的解析重建或迭代重建可实现复杂结构件的CT三维图像。在实际应用中，考虑到重建效率，可采用三维FDK解析重建算法进行重建。

在极坐标系下，FDK算法的表达式为：

其中β表示旋转角度，D表示光源到旋转中心的距离，权函数D/(D-S)中S为从重建点到虚拟探测器的距离，将虚拟探测器设置在旋转中心，可使算法分析和实现得到简化。由极坐标定义，对于某重建切层，

l是探测器上的线性坐标，

旋转过程中坐标系绕

轴旋转，转轴方向不做滤波，

表示锥形束投影数据，

是锥形束射线入射角的余弦函数，h(l)为一维斜坡滤波器的卷积核，常用R-L或S-L滤波函数。

如图2所示，包括：

低能量预扫描模块：采用预设低能量的X射线对被测物进行预设角度间隔△θ的全角度扫描，获得第一图像集；

起始角度定位模块：提取每个第一图像的工件投影区域和边缘轮廓信息；计算每个第一图像工件投影区域的像素灰度值，将最大像素灰度值或最小像素灰度值对应的第一图像的扫描角度记为θ1；

X射线高能量定位模块：固定θ1角度对被测物进行X射线能量递变成像，获得第二图像集，选择第二图像集中清晰度指标最大值对应的图像为目标第二图像，将目标第二图像的平均灰度值记为H1，目标第二图像对应的X射线能量的管电压记为V1、管电流记为I1；

高能量扫描模块：从θ1开始对被测物进行△θ间隔的全角度扫描，获得第三图像集；每个角度成像时X射线能量的管电流为I1，通过调整当前管电压V，使得当前角度对应的第三图像的平均灰度值H与H1的偏差值满足预设条件；

图像恢复模块：基于非线性扩散方程恢复第三图像集中因散射而损失掉的结构信息，获得内部结构完整的CT二维图像集，非线性扩散方程的初始条件为边缘轮廓信息。

图像恢复模块之后还包括：

三维重建模块：基于CT二维图像集重建被测物的CT三维图像，或采用FDK算法将CT二维图像集重建为被测物的CT三维图像。

可选地，110VA≤预设低能量≤130VA。

提取每个第一图像的工件投影区域和边缘轮廓信息，提取方法可选用二值化处理方法。

为了便于计算，清晰度指标可选用：图像的平均灰度值或图像中工件投影区域的平均灰度值。

高能量扫描模块中又包括：

初始成像模块：令当前管电压V＝V1、管电流为I1，获取当前角度的第三图像，计算第三图像的平均灰度值H；

判断模块：判断H与H1的偏差是否满足预设条件，如果是，则调整当前角度值，返回初始成像模块，如果否，则执行调整模块；

调整模块：根据H与H1的比值，调整当前管电压V，并更新当前角度的第三图像和对应的H，返回判断模块。

非线性扩散方程可设为：

或者设为：

假设第三图像水平与竖直方向的差分一致，则非线性扩散方程可简化为：

H1也可以设定为目标第二图像工件投影区域的平均灰度值，相应地H也设定为第三图像工件投影区域的平均灰度值。

需要说明的是，本发明的高动态X射线成像装置的实施例，与高动态X射线成像方法的实施例原理相同，相关之处可以互相参照。

以上所述仅为本发明的较佳实施例而已，并不用以限定本发明的包含范围，凡在本发明技术方案的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种高动态X射线成像方法，其特征在于，所述方法包括：

步骤13：计算每个第一图像所述工件投影区域的像素灰度值，将最大像素灰度值或最小像素灰度值对应的第一图像的扫描角度记为θ1；

步骤14：固定所述θ1角度对所述被测物进行X射线能量递变成像，获得第二图像集，选择所述第二图像集中清晰度指标最大值对应的图像为目标第二图像，将所述目标第二图像的平均灰度值记为H1，所述目标第二图像对应的所述X射线能量的管电压记为V1、管电流记为I1；

步骤15：从所述θ1开始对所述被测物进行所述△θ间隔的全角度扫描，获得第三图像集；每个角度成像时所述X射线能量的管电流为I1，通过调整当前管电压V，使得当前角度对应的第三图像的平均灰度值H与所述H1的偏差值满足预设条件；

步骤16：基于非线性扩散方程恢复所述第三图像集中因散射而损失掉的结构信息，获得内部结构完整的CT二维图像集，所述非线性扩散方程的初始条件为所述边缘轮廓信息。

2.根据权利要求1所述的方法，其特征在于，所述步骤16之后还包括：

步骤17：基于所述CT二维图像集重建所述被测物的CT三维图像。

3.根据权利要求1所述的方法，其特征在于，110VA≤所述预设低能量≤130VA。

4.根据权利要求1所述的方法，其特征在于，步骤12包括：基于二值化处理，提取每个第一图像的工件投影区域和边缘轮廓信息。

5.根据权利要求1所述的方法，其特征在于，所述清晰度指标为：图像的平均灰度值或图像中所述工件投影区域的平均灰度值。

6.根据权利要求1所述的方法，其特征在于，所述每个角度成像时所述X射线能量的管电流为I1，通过调整当前管电压V，使得当前角度对应的第三图像的平均灰度值H与所述H1的偏差值满足预设条件包括：

步骤151：令当前管电压V＝V1、管电流为I1，获取当前角度的第三图像，计算所述第三图像的平均灰度值H；

步骤152：判断所述H与所述H1的偏差是否满足预设条件，如果是，则调整所述当前角度值，返回步骤151，如果否，则执行步骤153；

步骤153：根据所述H与所述H1的比值，调整当前管电压V，并更新所述当前角度的第三图像和对应的所述H，返回步骤152。

7.根据权利要求1所述的方法，其特征在于，所述非线性扩散方程包括：

其中

表示图像灰度值的梯度，

作为图像的幅值表示了图像灰度值的变化幅度，扩散函数g取单调递减函数。

8.根据权利要求1所述的方法，其特征在于，所述非线性扩散方程包括：

其中，第一图像的边缘轮廓信息作为初始条件H₀，t作为尺度参数，H_t为最终投影恢复图像，

表示图像灰度值的梯度，

9.根据权利要求8所述的方法，其特征在于，令所述第三图像水平与竖直方向的差分一致，则所述非线性扩散方程为：

(H_x)_i,j表示第三图像水平方向差分，(H_y)_i,j表示第三图像竖直方向差分，H_ij表示第三图像在(i，j)像素位置的图像灰度值。

10.根据权利要求2所述的方法，其特征在于，所述步骤17包括：采用FDK算法将所述CT二维图像集重建为所述被测物的CT三维图像。

11.根据权利要求1-10任一所述的方法，其特征在于，所述H1为所述目标第二图像所述工件投影区域的平均灰度值，所述H为所述第三图像所述工件投影区域的平均灰度值。

12.一种非瞬时计算机可读存储介质，所述非瞬时计算机可读存储介质存储指令，其特征在于，所述指令在由处理器执行时使得所述处理器执行如权利要求1至11中任一所述的高动态X射线成像方法中的步骤。

13.一种高动态X射线成像装置，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1至11中任一所述的高动态X射线成像方法中的步骤。