CN101566590B - 面阵探测器射线数字成像中的散射强度分布获取方法 - Google Patents

Abstract

本发明公开了一种面阵探测器射线数字成像中的散射强度分布获取方法，分别采集1幅含滤波板的投影图像、1幅含滤波板和被检测物体的投影图像、1幅含射束衰减网格和被检测物体的投影图像以及1幅含射束衰减网格的投影图像，依次计算射束衰减网格中各小球投影中心的位置和射束衰减网格中各小球投影中心处的散射强度值，根据射束衰减网格中各小球投影中心处的散射强度值，采用Catmull-Clark细分曲面逼近的方法进行散射场图像拟合，得到散射强度分布。本发明不要求网格完全挡住透射的射线，而是利用网格引起射线的部分衰减来计算散射的强度分布，可应用于低、中、高各能量段的射线DR/CT扫描系统，具有良好的适用性与通用性。

Description

面阵探测器射线数字成像中的散射强度分布获取方法

技术领域

本发明属于射线无损检测领域，涉及一种面阵探测器射线数字成像中的散射强度分布获取方法。

背景技术

面阵探测器主要包括平板探测器和图像增强器两种，目前已广泛应用于射线DR(Digital Radiography)和CT(Computed Tomography)成像系统中，其中DR成像是CT成像的基础。

在成像过程中，由于X射线与被检测物体的相互作用，将产生以康普顿效应为主的各种散射，而散射是影响DR图像质量的重要因素，主要表现为降低图像的清晰度和对比度。散射的强度分布与很多因素有关，包括射线能量、物体特性(大小、形状、材质)、扫描位置关系等。散射一直是DR成像研究的重要问题，而对于二维开放成像的面阵探测器，散射的抑制与校正远比一维线阵探测器复杂。

散射的抑制与校正技术总体上可分为两类：一类为减少到达探测器的散射射线的技术，包括准直器、空气隙、滤线器、扫描狭缝技术、不透射线的铅条或铅板技术等，另一类为数字技术(即数字图像后处理技术)，包括高通滤波、散射卷积、散射反卷积等。近年来很多研究人员采用蒙特卡罗方法(Monte Carlo Method，MCM)进行散射仿真，通过专门的程序计算得到特定条件下的散射强度分布，进而分析散射的表现并进行相应的校正。采用蒙特卡罗程序仿真散射可以得到比较准确的理论结果，但这种方法每生成一幅投影图像都要跟踪数百万甚至上亿个光子，计算耗时数十甚至上百小时，实际应用的时间难以接受，因此多限于理论研究。文献Ruola Ning，Xiangyang Tang，D.L.Conover，X-Ray scatter suppression algorithm for cone beam volume CT，Proc.SPIEvol.4682，774-781，2002中使用Beam Stop Array(BSA)进行锥束CT散射伪影校正，该方法要求BSA中的铅球(或者铅盘、铅柱)完全挡住透射的X射线，对较低扫描电压的医学CT有较好的效果，但对扫描电压较高的工业CT，该方法所要求的“完全挡住射线”可能导致铅球尺寸过大而难以实施。另一方面，BSA本身也是散射源，如何进一步降低加入BSA后的影响，仍是需要研究的问题。

发明内容

为了克服现有技术的适用性与通用性的不足，本发明提供一种基于射束衰减网格(Beam Attenuation Array，BAA)的散射强度分布获取方法，该方法不要求网格完全挡住透射的射线，而是利用网格引起射线的部分衰减来计算散射的强度分布，可应用于低、中、高各能量段的射线DR/CT扫描系统。

本发明解决其技术问题所采用的技术方案包括以下步骤：

(1)在射线源与工作台之间尽量靠近工作台的位置放置滤波板，该位置以滤波板不与工作台和被检测物体发生碰撞干涉为准，通过面阵探测器采集1幅含滤波板的投影图像，该投影图像为投影图像1，本次扫描为扫描1；

(2)保持扫描1中的滤波板不动，将被检测物体置于工作台上，以扫描1的扫描参数采集1幅含滤波板和被检测物体的投影图像，该投影图像为投影图像2，本次扫描为扫描2；

(3)在滤波板放置位置将滤波板替换成射束衰减网格，被检测物体不动，以扫描1的扫描参数采集1幅含射束衰减网格和被检测物体的投影图像，该投影图像为投影图像3，本次扫描为扫描3；

(4)保持扫描3中的射束衰减网格不动，移除被检测物体，以扫描1的扫描参数采集1幅含射束衰减网格的投影图像，该投影图像为投影图像4，本次扫描为扫描4；

(5)根据投影图像4计算射束衰减网格中各小球投影中心的位置；

(6)计算射束衰减网格中各小球投影中心处的散射强度值；

(7)根据射束衰减网格中各小球投影中心处的散射强度值，采用Catmull-Clark细分曲面逼近的方法进行散射场图像拟合，得到散射强度分布，即散射强度分布是以散射场图像的形式表示。

在上述方法第1-2步中，用到的滤波板与射束衰减网格中用于固定小球的基板的材质与尺寸完全相同。

在上述方法第3-4步中，用到的射束衰减网格由基板和在其上均匀分布的小球阵列组成。对射束衰减网格的一般要求是：在满足小球对被检测物体的投影灰度值有一定影响量(投影图像3中所有小球投影区域的最小灰度值比投影图像2中该位置的灰度值减小20-80)的前提下，基板上的小球直径越小越好，并使小球球心的间距为小球直径的1.5-3.5倍；射束衰减网格基板材质与小球材质相同，或基板材质密度是小球材质密度的1-1.5倍；射束衰减网格基板厚度为小球直径的0.2-0.8倍，基板形状一般为矩形且不小于面阵探测器的成像窗口；小球被粘贴于基板的一个侧面上，所用粘贴剂的密度与小球材质密度相比越小越好，该侧面在扫描3中朝向射线源。对于医学DR/CT成像系统，典型的射束衰减网格由有机玻璃基板上粘贴均匀分布的塑料小球阵列组成；对于工业DR/CT成像系统，典型的射束衰减网格由铜基板上粘贴均匀分布的钢珠阵列组成。

在上述方法第1-4步中，从降低数字图像噪声和提高算法精度的角度考虑，可分别采集多幅投影图像，然后分别将各次扫描得到的多幅投影图像按对应象素进行灰度叠加并平均，得到1幅投影图像。

在上述方法第5步中，根据投影图像4计算射束衰减网格中各小球投影中心的位置的步骤如下：

1)将投影图像4复制一幅；

2)查找复制图像的小球阵列投影区域中的最小灰度；

3)以该最小灰度的110％-130％为阈值，对复制图像进行二值化，对小于阈值的象素视为位于小球投影区域，将其灰度设为1，其余象素灰度设为0，这样就得到与小球个数相同且灰度为1的投影区域，各小球投影中心必然位于这些投影区域之内；

4)生成一个位置-灰度链表，用于存储各小球投影中心的位置及相应灰度；

5)对二值化图像进行逐行逐象素扫描，若当前象素灰度为1，则进行如下处理(处理完后从下一个象素开始继续进行扫描)：以该象素为种子象素，设置一个比较位，该比较位初始值为种子象素对应于投影图像4中的位置及相应灰度；从种子象素开始用0值在二值化图像中进行区域填充，每填充一个象素就将投影图像4中该象素位置的灰度与比较位的灰度进行比较，若小于比较位的灰度，则用投影图像4中该象素的灰度与位置替换比较位的灰度与位置；若本次区域填充结束，则将比较位的位置与灰度加入位置-灰度链表尾。

在上述方法第6步中，计算射束衰减网格中各小球投影中心处的散射强度值的方法是相同的，以计算一个小球投影中心处的散射强度值为例，进行其计算方法的分析和推导。由于扫描1至扫描4的扫描参数完全相同，且射束衰减网格基板与滤滤板完全相同，因此可以认为射束衰减网格基板与滤滤板对投影图像灰度的贡献相同，从而认为其是探测器增益的一部分而忽略不计；考虑到小球直径较小，相应的散射值也较小，而且基板材质密度大于或等于小球材质密度，因此可以认为基板全部吸收了小球阵列所产生的散射射线。另一方面，由于小球阵列较被检测物体更远离面阵探测器，因此即使有少量透过基板的由小球阵列产生的散射射线，也会由于被检测物体对其的进一步吸收和空气隙效应，难以到达面阵探测器。

对扫描1，设象素G为某小球的投影中心，D为小球直径，象素G接收到的初始入射射线强度为I₀；对扫描4，象素G接收到的透射射线强度为I₂，根据Beer定律，有

$I_2=I_0{e}^{{-u}_b\left(E_{a0}\right)D}---\left(1\right)$

对扫描2，设初始入射射线强度为I₀，I₀透过被检测物体被象素G接收到的透射射线强度为I₁，象素G接收到的所有散射射线强度为S₁，则象素G接收到的总射线强度C₁＝I₁+S₁。以被检测物体为研究对象，有

$I_1=C_1-S_1$

$={\∫}_{E_{\min }}^{E_{\max }}{N}_0\left(E\right)e^{-{\∫}_0^L{u}_p\left(E\right)\mathrm{dL}}\mathrm{dE}$

$={\∫}_{E_{\min }}^{E_{\max }}{N}_0\left(E\right)e^{-{\mathrm{\μ}}_p\left(E_{a0}\right)L}\mathrm{dE}---\left(2\right)$

$=e^{-{\mathrm{\μ}}_p\left(E_{a0}\right)L}{\∫}_{E_{\min }}^{E_{\max }}{N}_0\left(E\right)\mathrm{dE}$

$=I_0{e}^{{-\mathrm{\μ}}_p\left(E_{a0}\right)L}$

对扫描3，设初始入射射线强度仍为I₀，I₀透过小球后的透射射线强度为I₂，I₂透过被检测物体被象素G接收的透射强度为I₃，象素G接收到的所有散射射线强度为S₂，象素G接收到的总射线强度C₂＝I₃+S₂。以被检测物体为研究对象，有

$I_3=C_2-S_2$

$={\∫}_{E_{\min }}^{E_{\max }}{N}_2\left(E\right)e^{{-\mathrm{\μ}}_p\left(E_{a2}\right)L}\mathrm{dE}---\left(3\right)$

$=e^{{-\mathrm{\μ}}_p\left(E_{a2}\right)L}{\∫}_{E_{\min }}^{E_{\max }}{N}_2\left(E\right)\mathrm{dE}$

$=I_2{e}^{-{\mathrm{\μ}}_p\left(E_{a2}\right)L}$

以上各式中，L为射线在象素G位置穿越被检测物体的长度，μ_b(E_a0)为小球在I₀的等效能量E_a0下的线性衰减系数，μ_p(E_a0)为被检测物体在I₀的等效能量E_a0下的线性衰减系数，μ_p(E_a2)为被检测物体在I₂的等效能量E_a2下的线性衰减系数，N₀(E)为I₀光子数量分布，N₂(E)为I₂光子数量分布。

式(2)除以式(3)，得

$\frac{C_2-S_2}{C_1-S_1}=\frac{I_2{e}^{-{\mathrm{\μ}}_p\left(E_{a2}\right)L}}{I_0{e}^{-{\mathrm{\μ}}_p\left(E_{a0}\right)L}}---\left(4\right)$

将式(1)代入上式，有

$\frac{C_2-S_2}{C_1-S_1}=\frac{e^{-u_b\left(E_{a0}\right)D}{e}^{-{\mathrm{\μ}}_p\left(E_{a2}\right)L}}{e^{-{\mathrm{\μ}}_p\left(E_{a0}\right)L}}=e^{-u_b\left(E_{a0}\right)D}{e}^{[{\mathrm{\μ}}_p\left(E_{a0}\right)-{\mathrm{\μ}}_p\left(E_{a2}\right)]L}---\left(5\right)$

在扫描3中，由于小球阵列的存在，使得作用于被检测物体的总射线强度比扫描2时略有减少，相应的散射射线强度也会略有减少，即S₂略小于S₁。另外，由于I₀穿越小球成为I₂的过程存在射束硬化效应，造成N₂(E)的平均能量略大于N₀(E)的平均能量，即μ_p(E_a2)略小于μ_p(E_a0)。设S₂＝αS₁(0＜α＜1)， $\mathrm{\β}=e^{-u_b\left(E_{a0}\right)D}$ (0＜β＜1)， $\mathrm{\γ}=e^{[{\mathrm{\μ}}_p\left(E_{a0}\right)-{\mathrm{\μ}}_p\left(E_{a2}\right)]L}$ (γ＞1)，则由式(5)可得扫描2中小球投影中心处的散射强度值为

$S_1=\frac{C_2-\mathrm{\β\γ}{C}_1}{\mathrm{\α}-\mathrm{\β\γ}}---\left(6\right)$

关于α，可以进一步作如下分析：小球阵列在投影范围内是近似均匀分布的，因此可以认为小球阵列对射线强度的减弱是近似均匀的，即不管是否有小球阵列，在小射束锥角内射线强度都是均匀分布的，只是大小略有不同。小球阵列对总射线强度的影响可以从投影图像1和投影图像4的比较中得到，本发明将α近似取为

关于β，由式(1)可得

$\mathrm{\β}=e^{-u_b\left(E_{a0}\right)D}=I_2/I_0---\left(8\right)$

关于γ，分析可知L增大会使μ_p(E_a0)-μ_p(E_a2)减小，即μ_p(E_a0)-μ_p(E_a2)是L的函数，而L是该位置灰度的函数，因此可以认为γ是小球投影中心灰度的函数。本发明将γ近似取为

$\mathrm{\γ}=e^{\frac{\mathrm{Gray}-\mathrm{MinGray}+c}{{\mathrm{MinGray}}^2}}---\left(9\right)$

其中Gray是投影图像3中当前计算的小球投影中心灰度，MinGray是该幅投影所有小球投影中心的最小灰度，c是常数，一般取为(0.1～0.3)MinGray。

一般可认为面阵探测器象素的输出灰度值与其所接收到的射线强度成正比，因此上述各式中象素G接收到的总射线强度可直接由该象素的灰度代替。上述分析与计算是对一个小球投影中心象素G进行的，对其余小球投影中心象素进行完全相似的计算，可以得到覆盖整个成像区域的所有小球投影中心处的散射值。

在上述方法第7步中，采用Catmull-Clark细分曲面逼近的方法进行散射场图像拟合时，由于该细分算法进行多层细分后，控制网格顶点的数量和密度会急剧增加，很容易就超过投影图像本身像素点的数量和密度，此时只需要对控制网格顶点的X和Y坐标进行四舍五入取整作为图像像素位置，顶点的相应Z坐标值就是该像素位置的散射值。Catmull-Clark细分层数主要取决于投影图像的分辨率，一般情况下细分4-7层即可。

上述的面阵探测器射线数字成像中的散射强度分布获取方法，可应用于医学和工业DR/CT成像系统散射校正中的获取散射强度分布步骤，得到散射场图像后散射校正的基本方法是：原始投影图像的各象素灰度减去与其对应的散射场图像的各象素灰度。

本发明的有益效果是：提出了一种新颖的面阵探测器射线数字成像中的散射强度分布获取方法，通过一定的扫描和计算可以得到比较准确的散射场图像，不要求网格完全挡住透射的射线，而是利用网格引起射线的部分衰减来计算散射的强度分布。该方法可应用于医学和工业领域的低、中、高各能量段的射线DR/CT扫描系统，具有良好的适用性与通用性。

下面结合附图和实施例对本发明进一步说明。

附图说明

图1为本发明获取面阵探测器射线数字成像中的散射强度分布流程；

图2为扫描1示意图；

图3为扫描2示意图；

图4为扫描3示意图；

图5为扫描4示意图；

图6为滤波板示意图；

图7为射束衰减网格示意图；

图8为计算一个小球投影中心散射强度中扫描2示意图；

图9为计算一个小球投影中心散射强度中扫描3示意图；

图10为投影图像2与获取的散射场图像的线性灰度比较。

具体实施方式

对一个直径为38mm，材料为铁的圆柱零件进行DR扫描，X射线源采用德国YXLON Y.TU 450-D02，面阵探测器采用美国Varian PaxScan2520，应用本发明方法获取散射强度分布，执行以下步骤：

(1)在射线源与工作台之间尽量靠近工作台的位置(不能与工作台和被检测物体发生碰撞干涉)放置滤波板，滤波板材质为铜，厚度1mm，扫描参数为电压350kV、电流0.22mA、探测器输出3幅/秒，通过面阵探测器采集1幅含滤波板的投影图像，该投影图像为投影图像1，本次扫描为扫描1；

(2)保持扫描1中的滤波板不动，将被检测物体置于工作台上，以扫描1的扫描参数采集1幅含滤波板+被检测物体的投影图像，该投影图像为投影图像2，本次扫描为扫描2；

(3)在滤波板放置位置将滤波板替换成射束衰减网格，射束衰减网格的基板材质和尺寸与滤波板完全相同，小球阵列为3mm直径的钢珠，间距5mm，被检测物体不动，以扫描1的扫描参数采集1幅含射束衰减网格+被检测物体的投影图像，该投影图像为投影图像3，本次扫描为扫描3；

(4)保持扫描3中的射束衰减网格不动，移除被检测物体，以扫描1的扫描参数采集1幅含射束衰减网格的投影图像，该投影图像为投影图像4，本次扫描为扫描4；

(5)根据投影图像4计算射束衰减网格中各小球投影中心的位置，具体步骤如下：

1)将投影图像4复制一幅；

2)查找复制图像的小球阵列投影区域中的最小灰度；

3)以该最小灰度的130％为阈值，对复制图像进行二值化，对小于阈值的象素视为位于小球投影区域，将其灰度设为1，其余象素灰度设为0，这样就得到与小球个数相同且灰度为1的投影区域，各小球投影中心必然位于这些投影区域之内；

4)生成一个位置-灰度链表，用于存储各小球投影中心的位置及相应灰度；

5)对二值化图像进行逐行逐象素扫描，若当前象素灰度为1，则进行如下处理：以该象素为种子象素，设置一个比较位，该比较位初始值为种子象素对应于投影图像4中的位置及相应灰度；从种子象素开始用0值在二值化图像中进行区域填充，每填充一个象素就将投影图像4中该象素位置的灰度与比较位的灰度进行比较，若小于比较位的灰度，则用投影图像4中该象素的灰度与位置替换比较位的灰度与位置；若本次区域填充结束，则将比较位的位置与灰度加入位置-灰度链表尾。处理完后从下一个象素开始继续进行扫描。

(6)计算射束衰减网格中各小球投影中心处的散射强度值，以步骤(5)中计算得到的一个小球投影中心(798，536)为例，投影图像2中该位置的象素灰度C₁＝1739，投影图像3中该位置的象素灰度C₂＝1289，由式(7)计算得到α＝0.91，由式(8)计算得到β＝0.70，由式(9)计算得到γ＝1.01，代入式(6)计算得到 $S_1=\frac{C_2-\mathrm{\β\γ}{C}_1}{\mathrm{\α}-\mathrm{\β\γ}}=293.24$ 依次计算所有小球投影中心的散射强度值。

(7)根据射束衰减网格中各小球投影中心处的散射强度值，采用Catmull-Clark细分曲面逼近的方法进行散射场图像拟合，细分6层后，对控制网格顶点的X和Y坐标进行四舍五入取整作为图像像素位置，顶点的相应Z坐标值就是该像素位置的散射值，从而得到散射强度分布，即散射强度分布是以散射场图像的形式表示。

从图10给出的投影图像2与获取的散射场图像的线性灰度比较可以看出，本发明获取的散射场图像灰度变化平滑，且与投影图像2的灰度存在一定的对应关系，这与理论分析的散射强度分布情况基本相符，表明本发明方法的可行性和有效性。

Claims (6)

1.面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于包括下述步骤：

(1)在射线源与工作台之间尽量靠近工作台的位置放置滤波板，该位置以滤波板不与工作台和被检测物体发生碰撞干涉为准，通过面阵探测器采集1幅含滤波板的投影图像，该投影图像为投影图像1，本次扫描为扫描1；

(2)保持扫描1中的滤波板不动，将被检测物体置于工作台上，以扫描1的扫描参数采集1幅含滤波板和被检测物体的投影图像，该投影图像为投影图像2，本次扫描为扫描2；

(3)在滤波板放置位置将滤波板替换成射束衰减网格，被检测物体不动，以扫描1的扫描参数采集1幅含射束衰减网格和被检测物体的投影图像，该投影图像为投影图像3，本次扫描为扫描3；

(4)保持扫描3中的射束衰减网格不动，移除被检测物体，以扫描1的扫描参数采集1幅含射束衰减网格的投影图像，该投影图像为投影图像4，本次扫描为扫描4；

(5)根据投影图像4计算射束衰减网格中各小球投影中心的位置，包括以下步骤：

1)将投影图像4复制一幅；

2)查找复制图像的小球阵列投影区域中的最小灰度；

3)以该最小灰度的110％-130％为阈值，对复制图像进行二值化，对小于阈值的象素视为位于小球投影区域，将其灰度设为1，其余象素灰度设为0，得到与小球个数相同且灰度为1的投影区域，各小球投影中心必然位于这些投影区域内；

4)生成一个位置-灰度链表，用于存储各小球投影中心的位置及相应灰度；

5)对二值化图像进行逐行逐象素扫描，若当前象素灰度为1，则进行如下处理：以该象素为种子象素，设置一个比较位，该比较位初始值为种子象素对应于投影图像4中的位置及相应灰度；从种子象素开始用0值在二值化图像中进行区域填充，每填充一个象素就将投影图像4中该象素位置的灰度与比较位的灰度进行比较，若小于比较位的灰度，则用投影图像4中该象素的灰度与位置替换比较位的灰度与位置；若本次区域填充结束，则将比较位的位置与灰度加入位置-灰度链表尾；处理完后从下一个象素开始继续进行扫描；

(6)根据

计算得到该象素接收到的所有散射射线强度

依次计算射束衰减网格中所有小球投影中心的散射强度值，其中，I₀为象素G接收到的初始入射射线强度，I₂为扫描1时象素G接收到的透射射线强度，Gray是投影图像3中当前计算的小球投影中心灰度，MinGray是该幅投影所有小球投影中心的最小灰度，c是常数，取为(0.1～0.3)MinGray；

(7)根据射束衰减网格中各小球投影中心处的散射强度值，采用Catmull-Clark细分曲面逼近的方法进行散射场图像拟合，得到散射强度分布。

2.根据权利要求1所述的面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于：所述的滤波板与射束衰减网格中用于固定小球的基板的材质与尺寸完全相同。

3.根据权利要求1所述的面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于：所述的射束衰减网格由基板和在其上均匀分布的小球阵列组成，在满足投影图像3中所有小球投影区域的最小灰度值比投影图像2中该位置的灰度值减小20-80的前提下，基板上的小球直径越小越好，并使小球球心的间距为小球直径的1.5-3.5倍；射束衰减网格基板材质与小球材质相同，或基板材质密度是小球材质密度的1-1.5倍；射束衰减网格基板厚度为小球直径的0.2-0.8倍，基板形状为矩形且不小于面阵探测器的成像窗口；小球被粘贴于基板的一个侧面上，所用粘贴剂的密度与小球材质密度相比越小越好，该侧面在扫描3中朝向射线源。

4.根据权利要求1所述的面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于：所述的射束衰减网格在医学直接数字化照相和计算机射线照相成像系统中由有机玻璃基板上粘贴均匀分布的塑料小球阵列组成；所述的射束衰减网格在工业直接数字化照相和计算机射线照相成像系统中由铜基板上粘贴均匀分布的钢珠阵列组成。

5.根据权利要求1所述的面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于：所述的步骤(1)至步骤(4)中，分别采集多幅投影图像，然后分别将各次扫描得到的多幅投影图像按对应象素进行灰度叠加并平均，得到1幅投影图像。

6.根据权利要求1所述的面阵探测器射线数字成像中的散射强度分布获取方法，其特征在于：所述的步骤(7)中，Catmull-Clark细分层数为4-7层，此时对控制网格顶点的X和Y坐标进行四舍五入取整作为图像像素位置，顶点的相应Z坐标值就是该像素位置的散射值。

Images