A Correction Method for Hardening Artifacts in CT Images Based on Integral Invariance
-
摘要:
常规CT系统中，X射线是多能的，导致CT图像中出现硬化伪影。对于多材料成像对象，在硬化校正时通常需要能谱或材料等先验信息。对此，本文提出一种适用于单电压条件下多材料物体CT成像的硬化伪影校正方法。该方法将X射线透射率图像看作是多个单能透射率图像的非负加权组合，以高斯分布的极大似然函数为目标函数，并结合不同角度下投影积分不变性作为约束条件，构建X射线透射率分解模型。通过求解该模型，可以获得窄能谱投影信息，进而重建出不同能量下的CT图像。实验结果表明，相比于直接重建的CT图像，本文方法重建得到的图像硬化伪影明显减弱。
Abstract:In conventional computed tomography (CT) systems, X-ray beams used for imaging are polyenergetic. These beams cause hardening artifacts in the reconstructed images owing to their polychromatic nature. For a multi-material imaging object, prior information about the X-ray spectrum and composition of the imaging object is required for correcting the beam hardening. To avoid relying on such prior information, this study proposes a beam-hardening artifact correction method for CT imaging of multi-material objects under single voltage. In this method, the X-ray image is considered as a non-negative weighted sum of multiple single-energy X-ray images. The decomposition model for the X-ray transmission images is constructed using the maximum likelihood function of a Gaussian distribution as the objective function under the constraint of the projection integral invariance at different angles. To solve this model, the projection of the spatial terms on the X-ray path is obtained. Thus, CT images for X-ray beams with different energies can be reconstructed. For imaging experiments with dual-material and multi-material models, compared with directly reconstructed CT images, the hardening artifacts in the reconstructed images were significantly reduced using this method, indicating its effectiveness.
-
X射线计算机断层扫描技术（computed tomography，CT）作为一种重要的成像技术，在无损检测和材料分析领域起重要作用。在常规CT系统中，X射线是连续能谱的。当X射线穿过物体时，低能光子穿透能力弱，高能光子穿透能力强，低能射线相较于高能射线更容易被衰减，通过物体后的平均能量比通过物体前的平均能量高，会出现“射束硬化”现象[1]。
射束硬化伪影校正是CT研究中的热点，实现方法通常可以分为两大类：硬件方法和软件方法[2]。硬件方法是指通过改进CT成像设备来减少射束硬化效应，其中典型的方法有滤波片校正、光子计数探测器等。滤波片校正法是对多能射线束进行滤波，可以在一定程度上抑制硬化伪影现象[3-4]。光子计数探测器可以通过合并一个较低阈值以上的光子，减少了较低能量光子的影响，从而降低射束硬化效应[5]。
软件校正方法主要有深度学习法、迭代法、线性化法和多能成像法等[2]。深度学习法对大规模样本进行训练，提取射束硬化伪影的特征并自动校正[6-8]，但需要标注大量数据，存在过拟合和泛化性差的问题。迭代法是一种不断用变量的旧值递推新值的过程。Brabant等[9]提出在迭代重建过程时，对射束硬化效应进行建模，将其融合到正投影过程中，达到消除伪影的目的。线性化法是通过线性变换对非线性数据进行校正，从而实现硬化伪影校正。通常可以利用多项式拟合将多能连续谱投影数据转化为单能投影数据，但该方法需要利用特定模体投影数据或其它先验信息进行估计，且难以适用于多材料构成的成像对象[10-12]。Ingacheva等[13]将投影数据进行单参数线性化，利用各角度下投影之和相同，通过求解参数最优解，对投影数据校正后进行重建。多能成像法是通过对多种不同能谱下的投影数据分解，获取基投影图像，进而合成出近似单能图像。对于物质成分已知且彼此不混合的被测物，罗婷等[14]把被检测物体的已知物质作为基材料，对CT投影数据采集过程进行数学建模，对基材料图像迭代求解，得到伪单能图像，可以显著地校正硬化伪影。Wei等[15-16]利用盲分解法，将能谱、材料信息作为未知变量，对多电压投影进行分解，得到基图像和能量，但该方法需要多个电压数据，不适用于单电压CT成像。Sajja等[17]改进了双能锥束CT技术，改善软组织对比度，减弱射束硬化伪影和金属伪影，该方法同样需要至少两个电压下的投影数据。
尽管多材料物体的射束硬化校正方法很多，但除深度学习法外，大多需要能谱信息，材料组分等作为先验信息。在这些先验信息未知条件下，则需要获取多电压投影。对此，本文提出一种单电压CT成像射束硬化伪影校正方法，适用于多材料物体成像，且不需要能谱，材料组分等先验信息。
1. 理论与方法
实际中X射线是连续能谱，然而传统CT重建算法是基于单能假设，导致重建图像中会出现硬化伪影。利用窄能谱X射线图像进行重建，可以减弱硬化伪影的影响。连续谱X射线图像可以看作是窄能谱X射线图像的线性组合，那么可以构建多能X射线透射率[18]分解模型。对X射线透射率模型进行求解，可以获取窄能谱X射线图像，从而实现硬化伪影校正。
1.1 多能X射线透射率分解模型
在实际应用中，X射线源发出的是连续能谱X射线，衰减系数不仅与位置相关，还与射线能量有关。根据Beer定理[19]
$$ I = {I_0}\int_0^{{E_{\max }}} {S(E){e^{ - \int_L {u(x,E)dx} }}dE} {\text{ }}, $$ (1) 式中：I为探测器接收到的射线强度，I0为射线的初始强度，Emax表示成像电压下X射线最高能量，S(E)表示归一化等效能谱，与X射线能谱、探测器效率等因素有关，u(x,E)表示路径L上x处，能量为E时X射线的衰减系数。S(E)满足
$$ 1 = \int_0^{{E_{\max }}} {S(E)dE{\text{ }}} . $$ (2) 对（2）进行离散化，第r个窄能谱段能量近似为Er，可以得到
$$ \frac{I}{{{I_0}}} = \sum\limits_{r = 1}^R {{s_r}{e^{ - \sum\limits_{k = 1}^K {{u_{rk}}{d_k}} }}} {\text{ }}, $$ (3) 其中，I/I0为X射线透射率，sr为第r个窄能谱段X射线透射率对应的加权系数，k表示第k种基材料或基效应。urk为第k种基材料在能量Er时的衰减系数，dk为X射线穿过路径上第k种基材料（效应）的等效长度。其中
$$ 1 = \sum\limits_{r = 1}^R {{s_r}{\text{ }}} , $$ (4) 将投影图像中第m个像素对应的透射率记为fm，m=1,2,
$\cdots $ ,M。由式（4），将所有像素的透射率写成矩阵形式为$$ F = S{e^{ - UD}}{\text{ }}, $$ (5) 其中，F=(f1,f2,
$\cdots $ ,fM)，为投影图像各像素对应的X射线透射率构成的图像。S=(s1,s2,$\cdots $ ,sR)，为构成X射线透射率图像的窄能谱图像加权系数。U=(urk)RK，U的一列对应一材料在各窄谱段的衰减系数，可以通过查询材料衰减系数数据库，得到各基材料在各能量段的衰减系数。D=(dkm)KM，D的一行是同种材料在各像素上的衰减长度。UD的每一行对应一个窄能谱投影。根据S，D的物理意义，S和D为非负矩阵，即S≥0，D≥0。由于各探元获得的测量值相互独立，所以考虑根据各探元对应的fm所服从的分布建立极大似然估计模型，对式（5）中的能谱权重矩阵S，各材料在各个像素中的衰减长度D进行估计。
X射线的成像过程中存在量子噪声、光子噪声和电子噪声。量子噪声是由于检测有限数量的X射线量子时的波动而出现的，是CT图像中的统计噪声产生的主要原因[20]。光子噪声是由于光子的随机到达引起的噪声，通常与光的粒子性相关。电子电路接收模拟信号的过程可能会受到一些噪声的影响，这些噪声被称为电子噪声。3种噪声都被认为是CT图像中的噪声，通常被认为服从加性高斯分布[20]，即
$$ {f}_{m}\text{ }～\text{ }N({\tau }_{m},{\sigma }_{m}^{2})\text{ }, $$ (6) 式中，τm为fm的期望，σm为fm的标准差。将式（6）化为标准正态分布，可以得到
$$ \frac{{f}_{m}-{\tau }_{m}}{{\sigma }_{m}}～\text{ }N(0,\text{ }1)\text{ }, $$ (7) 概率分布函数为
$$ \varphi \left( {\frac{{{f_m} - {\tau _m}}}{{{\sigma _m}}}} \right) = \frac{1}{{\sqrt {2\pi } {\sigma _m}}}{e^{ - \tfrac{{{{({f_m} - {\tau _m})}^2}}}{{2\sigma _m^2}}}}{\text{ }}, $$ (8) 在fm已知的条件下，τm为第m探测单元上透射率的期望值[21]，即
$$ {\tau _m} = \sum\limits_{r = 1}^R {{s_r}} {e^{ - \sum\limits_{k = 1}^K {{u_{rk}}{d_{km}}} }}{\text{ }}, $$ (9) 用ωm对式（8）中σm进行估计，ωm其中为各角度下每点的局部方差，计算过程如式（10）所示：
$$ {\omega _m} = \frac{{\displaystyle\sum\limits_{m \in \Theta } {{{(f(m) - \hat f(m))}^2}} }}{N}{\text{ }}, $$ (10) 其中，f (m)代表X射线透射率图像中第m个像素的f值，
$ \hat{f}(m) $ 代表X射线透射率图像中以第m个像素点为中心，以N为窗口宽度的平均f值，Θ为像素点m所构成的集合。因此，建立高斯对数似然函数$$ \begin{split}L({s_r},{d_{km}}) =\;& - \frac{1}{2}m\ln (2\pi ) - \\ \;& \frac{1}{2}\sum\limits_{m = 1}^M {{{\left( {\left( {{f_m} - \sum\limits_{r = 1}^R {{s_r}{e^{ - \sum\limits_{k = 1}^K {{u_{rk}}{d_{km}}} }}} } \right)/{\omega _m}} \right)}^2}} . \end{split}$$ (11) 那么，式（11）可以转化为优化问题：
$$ \min \sum\limits_{m = 1}^M {{{\left( {{f_m} - \sum\limits_{r = 1}^R {{s_r}{e^{ - \sum\limits_{k = 1}^K {{u_{rk}}{d_{km}}} }}} } \right)}^2}} /\omega _m^2{\text{ }}. $$ (12) 由于各像素点具有异方差性，可以使用加权最小二乘估计思想[22]对式（12）进行模型转化，进而对sr，dkm进行估计。结合式(4)，得到优化模型如下：
$$ \begin{gathered} \min \left\| {\frac{{F - S{e^{ - UD}}}}{\Omega }} \right\|_F^2 \\ s.t. S \geqslant 0,{\text{ }}D \geqslant 0,{\text{ }}1 = \sum\limits_{r = 1}^R {{s_r}} {\text{ }}. \\ \end{gathered} $$ (13) 其中，Ω=(ω1,ω2 ……,ωM)，为X射线透射率图像各角度下每点的局部方差矩阵。在单电压条件下，式（13）中的解具有不确定性，所以引入其它限制条件。
1.2 基于积分不变性的多能X射线透射率分解模型
CT 成像扫描中物体本身是不变的，各材料的分布也是不变的。以二维CT成像为例，平行束X射线的投影过程如图1所示，图中阴影区域H表示一个单材料物体。以s表示探测器的坐标位置，射线路径可以表示为L(s)。射线沿着路径L(s)穿过物体时，穿过物体的长度记为d(s)，则d(s)对s的积分表示阴影区域H的面积S面积，即
$$ {S}_{面积}={\displaystyle {\int }_{-\infty }^{+\infty }d\left(s\right)\text{ }}\text{d}s. $$ (14) 这一积分值与角度无关，由物体本身决定。对于单能X射线，如果衰减系数为u，则沿着路径L(s)的投影为ud(s)，从而
$$ u{S}_{面积}={\displaystyle {\int }_{-\infty }^{+\infty }ud\left(s\right)\text{ }}\text{d}s, $$ (15) 即uS面积与旋转角度无关，此为投影积分不变性[13,22-23]。同理，在多材料时，各材料对应的d(s)，其积分为各材料的面积，其值仅由物体本身决定。在不同角度下，是不变的。
假设第j个扫描角度下X射线投影像素数为Mj（j=1, 2,……, J），将所有扫描角度下投影图像像素数之和记为M，则
$$ M = \sum\limits_{j = 1}^J {{M_j}{\text{ }}{\text{.}}} $$ (16) 在平行束X射线成像下，对式（14），在离散情况下，积分变为求和，从而第k种材料在第j个角度下的面积可表示为
$$ {T_{kj}} = \sum\limits_{i = 1}^{{M_j}} {{d_{k\left( {i + \sum\limits_{g = 1}^{j - 1} {{M_g}} } \right)}}} {\text{ }}{\text{.}} $$ (17) 在不同角度下，该值均相等，因此有
$$ {T_{kj}} = {T_{k(j + 1)}},j = 1,2, \cdots ,J - 1,k = 1,2, \cdots ,K{\text{ }}. $$ (18) 在实际成像中，当X射线扇束角或锥角较小时，可近似看作平行束，用式(16)近似[13]。
由式（13）和式（16）及变量固有的非负性，可得X射线透射率分解模型如下：
$$ \begin{gathered} \min \left\| {\frac{{F - S{e^{ - UD}}}}{\Omega }} \right\|_F^2 \\ s.t. \;\;{\text{ }}S \geq 0{\text{ }}; \\ \qquad D \geq 0{\text{ ;}} \\ \qquad 1 = \sum\limits_{r = 1}^R {{s_r}} {\text{ ;}} \\ \qquad {T_{kj}} = {T_{k(j + 1)}},j = 1,2, \cdots , J - 1, \\ \qquad k = 1,2, \cdots ,K{\text{.}} \end{gathered} $$ (19) 由于X射线透射率分解模型是一个具有非负约束的盲分离模型，仿照非负矩阵问题算法的求解方法，利用KKT条件对式（19）求解算法进行推导。在该模型中，将能谱段等间距划分，每个能量段以中间能量作参考，确定U值[16]。不考虑能谱归一化约束和积分不变性约束，可得S，D的迭代求解公式分别为
$$ S = S \odot \dfrac{{\left( {\dfrac{{S\left( {{e^{ - UD}}} \right)}}{\Omega }} \right){{\left( {{e^{ - UD}}} \right)}^T}}}{{\dfrac{F}{{{\Omega ^2}}}{{\left( {{e^{ - UD}}} \right)}^T}}}, $$ (20) $$ D = D \odot \dfrac{{{U^T}\left( {\left( {{S^T}\left( {\dfrac{{S\left( {{e^{ - UD}}} \right)}}{{{\Omega ^2}}}} \right)} \right) \odot \left( {{e^{ - UD}}} \right)} \right)}}{{{U^T}\left( {\left( {{S^T}\dfrac{F}{{{\Omega ^2}}}} \right) \odot \left( {{e^{ - UD}}} \right)} \right)}}{\text{ }}. $$ (21) 在每次对S进行更新后，对其进行归一化，确保其满足式(4)约束。然后对D的值重新进行强制分配，使其尽量满足式（21）约束，具体为
$$ {d_{k\left( {i + \sum\limits_{g = 1}^{j - 1} {{M_g}} } \right)}} = {d_{k\left( {i + \sum\limits_{g = 1}^{j - 1} {{M_g}} } \right)}}\left( {\frac{{\displaystyle\frac{1}{J}\displaystyle\sum\limits_{m = 1}^M {{d_{km}}} }}{{\displaystyle\sum\limits_{i = 1}^{{M_j}} {{d_{k\left( {i + \sum\limits_{g = 1}^{j - 1} {{M_g}} } \right)}}} }}} \right){\text{ }}. $$ (22) 由于X射线能谱的具体值是未知的，在某一能量段下，中间能量的强度比两端能量的强度相对较大，因此用二次函数作为能谱的近似。假定能谱等间隔划分为R段，划分能量为
$$ 0 = {E_0} < {E_1} < \cdots < {E_R} = {E_{\max }}{\text{ }}, $$ (23) 则sr初始值取值为
$$ {s_r} = \frac{{\displaystyle\int_{{E_{r - 1}}}^{{E_r}} {E\left( {{E_{\max }} - E} \right)} {\text{ }}dE}}{{\displaystyle\int_0^{{E_{\max }}} {E\left( {{E_{\max }} - E} \right)} {\text{ }}dE}}{\text{ }}, $$ (24) 对dkm初始值取值为
$$ {d_{km}} = \frac{{{f_m}}}{{K{u_{\bar rk}}}},\bar r = \left[ {0.5 + \frac{R}{2}} \right]{\text{ }}, $$ (25) 其中“[ ]”表示取整。
算法流程图如图2。
2. 实验结果
为了验证方法对射束硬化伪影的校正作用，分别使用双材料模体（钛合金与钢）、多材料模体（铝、氧化铝及氮化铝）进行实验验证。第1个模体外部为钛合金材料，内部圆形区域为钢材料。钛合金部分的横截面为圆形的一部分，直径约为10.0 mm，钢部分直径约为2.5 mm。第2个模体铝、氧化铝及氮化铝3种材料横截面均为圆形，直径约为5.0 mm。CT成像系统为中北大学智能探测技术与装备山西省重点实验室的YXLON FF20 CT系统，采用FBP算法进行重建。在单电压条件下，对于多材料物体的硬化伪影校正，文献[9]中的方法不需要能谱和材料信息，使用迭代重建算法，对硬化伪影进行校正；文献[13]中利用投影积分不变性对投影数据进行单参数线性化，校正后进行重建。所以选择文献[9]与文献[13]中方法重建结果与本算法重建结果、直接重建结果进行对比。
2.1 钛合金与钢实验
样品一实验中，管电压为160 kV，管电流为40 μA。射线源与物体中心的距离为200.000 mm，射线源与探测器的距离为779.577 mm，采样间隔为1°，采样角度为360°。在上述距离条件下，X射线扇束角最大值为1.724°，可将扇形束投影数据近似看成是平行束投影数据。将能量范围等间隔划分为16个能量区间，记R=16，每个能量区间10 keV。选择光电效应和康普顿效应为基效应，记K=2。本实验中选择窗口宽度经验值为5，计算X射线透射率图像各角度下每点的局部方差，记N=5。分解后得到投影图像进行重建，重建图像大小为256×256，像素大小为0.045 mm×0.045 mm。为了方便进行对比，所以选择与直接重建图像灰度值相近的能量区间分解投影的重建图进行比较，此处选择第十能量区间下的投影重建图。对于文献[9]中的算法，参考文章中选取参数α=0.004，β=3。
直接重建图像如图3（a）所示，积分不变性线性化校正方法重建图像如图3（b）所示，迭代校正方法重建图像如图3（c）所示，本文方法重建图像如图3（d）所示。不难发现，直接重建图像中，硬化伪影明显，尤其是钛合金区域，图像边缘像素灰度值高于中间部分。如图3所示，对比重建图像的同一列灰度值变化，本文方法硬化伪影校正后的重建图像，钛合金区域灰度值差异明显变小，优于直接重建图像和对比方法校正后的重建图像。
为了更加直观的观察伪影情况，对图像的伪彩色图进行对比。从图4（a）可知，直接重建的图像存在较为严重的硬化伪影，钛合金部分明显中间为黄绿色，边缘为红色，灰度值不均匀；从图4（b）观察，边缘亮度为蓝色，数值小于中间亮度，说明存在过校正现象；图4（c）观察发现，钛合金区域仍存在硬化伪影，校正效果不佳，表明迭代校正算法中，模型的参数选取对材料组分依赖性较强，文献中参数不适用于钛合金材料的伪影校正；相比之下，本文方法所得图像中间和边缘数值接近，图像更为均匀。
![]()
图 4 重建图像一列（图3（a）红线位置）灰度值变化Figure 4. Change in the grayscale values of a column [red line position in Figure 2 (a)] in the reconstruction result射束硬化伪影会导致重建图像的中间区域灰度值偏低，边缘区域灰度值偏高，极差值相对较大。极差均值比（extreme mean ratio，EMR）可以作为量化伪影的严重程度的指标，计算方法如式（26）所示，
$$ EMR = {G_{region}}/{\mu _{region}}{\text{ }}. $$ (26) 式中：Gregion为区域极差，μregion为区域均值。EMR越小，表示图像均匀性越好，像素值一致性越高。所以此处采用EMR对重建结果进行评价。表1展示了模体重建图像中钛合金和钢的极差、均值和EMR。
表 1 不同方法下极差均值比Table 1. Mean range ratio for different methods直接重建 积分不变性线性化校正方法重建 迭代校正方法重建 本文方法重建 钛合金 钢 钛合金 钢 钛合金 钢 钛合金 钢 极差Gregion 0.161 0.169 0.309 0.355 0.236 0.156 0.118 0.149 均值μregion 0.183 0.259 0.656 0.654 0.270 0.264 0.157 0.245 EMR 0.881 0.655 2.119 0.543 0.875 0.591 0.750 0.607 本方法校正后，钛合金区域和钢的EMR值都低于直接重建图像的EMR值。尽管积分不变性线性化方法与迭代校正方法中对钢区域的射束硬化形成校正作用相对较好，但钛合金区域的EMR值高于本方法重建得到的钛合金区域的EMR值，说明图像的均匀性不佳。灰度值对比曲线、伪彩色图以及EMR值均表明，本方法对重建图像的射束硬化伪影现象具有一定抑制作用。
2.2 铝、氧化铝及氮化铝实验
样品二实验中，管电压为80 kV，管电流为40 μA。射线源与物体中心的距离为400.000 mm，射线源与探测器的距离为780.577 mm，采样间隔为1°，采样角度为360°。在上述距离条件下，X射线扇束角最大值为0.755°，可将扇形束投影数据近似看成是平行束投影数据。将能量范围等间隔划分为8个能量区间，记R=10，每个能量区间为10 keV。选择光电效应和康普顿效应为基效应，记K=2。本实验中选择窗口宽度经验值为5，计算X射线透射率图像各角度下每点的局部方差，记N=5。分解后得到投影图像，进行重建，重建图像大小为256×256，重建像素大小为
0.0635 mm×0.0635 mm。对于文献[9]中的算法，参考文章中选取参数α=0.004，β=3。此处选择第五能量区间下的投影重建图与直接重建结果、积分不变性线性化校正[13]重建结果、迭代校正[9]重建结果进行对比（图5）。其中，左侧为铝，右侧上方为氮化铝，右侧下方为氧化铝。![]()
图 7 重建图像一列（图6（a）红线位置）灰度值变化Figure 7. Change in the grayscale values of a column [red line position in Figure 5 (a)] in the reconstruction result对比重建图像的同一条线上的灰度值变化，不难发现，本文方法硬化伪影校正后的重建图像，各材料的硬化伪影得到了很好的校正，边缘亮度，优于直接重建图像和积分不变性线性化方法与迭代校正方法校正后的重建图像。
为了观察伪影情况，采用图像的伪彩色图进行对比，如图8所示。在直接重建图像中，可以观察到三种材料的边缘值均大于内部值，存在明显的硬化伪影。对比方法重建图像，均匀性不佳。相比之下，本文方法得到的图像同组分区域颜色相近，图像更为均匀，表明本方法对重建图像的射束硬化伪影现象具有抑制作用。
![]()
图 8 重建图像一列（图6（a）绿线位置）灰度值变化Figure 8. Change in the grayscale values of a column [green line position in Figure 5 (a)] in the reconstruction result通过计算EMR值，不难发现，对于多材料模体的硬化伪影校正，积分不变性线性化方法校正效果不佳，迭代校正方法对铝和氧化铝材料的硬化伪影可以起到校正效果，但对氮化铝材料不能进行硬化伪影校正。相比之下，本文方法的射束硬化形成校正作用明显，得到的图像均匀性更好。灰度值对比曲线与EMR值均表明，本方法对多材料模体重建图像的射束硬化伪影现象具有一定抑制作用。
表 2 不同方法下极差均值比Table 2. Mean range ratio under different methods直接重建 积分不变性线性化校正方法重建 迭代校正方法重建 本文方法重建 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 极差Gregion 0.018 0.038 0.018 0.043 0.059 0.026 0.017 0.038 0.021 0.015 0.031 0.009 均值μregion 0.021 0.028 0.018 0.028 0.043 0.024 0.021 0.029 0.018 0.021 0.024 0.015 EMR 0.870 1.370 1.016 1.507 1.351 1.081 0.780 1.309 1.163 0.734 1.277 0.629 3. 结论
针对单电压条件下多材料物体CT成像中存在硬化伪影的问题，本文提出了一种基于投影积分不变性约束的单电压CT图像硬化伪影校正方法。该方法以投影积分不变性作为约束条件，构建了X射线透射率分解模型，求解模型以获得窄能谱投影，进而重建得到硬化伪影校正的图像。实验结果中，该方法所得CT图像相比于直接重建、对比方法重建的CT图像，硬化伪影明显减弱。理论上，单电压下，模型（13）必然是多解的，求解所得窄能谱投影未必是真正的窄能谱投影，并且其解会受到始值影响。利用投影积分不变性，约束了模型（13）解的结构，从而一定程度上达到硬化伪影校正的目的。想要实现真正的窄能谱成像，仍然需要多个电压下的投影数据。此外，在实际成像中，对扇束和锥束可以近似转化为平行束的角度阈值界定，需要进一步研究。
-
![]()
图 4 重建图像一列（图3（a）红线位置）灰度值变化
Figure 4. Change in the grayscale values of a column [red line position in Figure 2 (a)] in the reconstruction result
![]()
图 7 重建图像一列（图6（a）红线位置）灰度值变化
Figure 7. Change in the grayscale values of a column [red line position in Figure 5 (a)] in the reconstruction result
![]()
图 8 重建图像一列（图6（a）绿线位置）灰度值变化
Figure 8. Change in the grayscale values of a column [green line position in Figure 5 (a)] in the reconstruction result
表 1 不同方法下极差均值比
Table 1 Mean range ratio for different methods
直接重建 积分不变性线性化校正方法重建 迭代校正方法重建 本文方法重建 钛合金 钢 钛合金 钢 钛合金 钢 钛合金 钢 极差Gregion 0.161 0.169 0.309 0.355 0.236 0.156 0.118 0.149 均值μregion 0.183 0.259 0.656 0.654 0.270 0.264 0.157 0.245 EMR 0.881 0.655 2.119 0.543 0.875 0.591 0.750 0.607 
下载: 导出CSV
表 2 不同方法下极差均值比
Table 2 Mean range ratio under different methods
直接重建 积分不变性线性化校正方法重建 迭代校正方法重建 本文方法重建 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 铝 氧化铝 氮化铝 极差Gregion 0.018 0.038 0.018 0.043 0.059 0.026 0.017 0.038 0.021 0.015 0.031 0.009 均值μregion 0.021 0.028 0.018 0.028 0.043 0.024 0.021 0.029 0.018 0.021 0.024 0.015 EMR 0.870 1.370 1.016 1.507 1.351 1.081 0.780 1.309 1.163 0.734 1.277 0.629
下载: 导出CSV
-
[1] 张秀琰, 陈明, 郑永果. 一种X射线射束硬化校正方法[J]. 软件, 2020, 41(2): 256-259. ZHANG X Y, CHEN M, ZHENG Y G, A Correction Method of X-ray Beam Hardening[J]. Computer engineering & Software, 2020, 41(2): 256-259. (in Chinese).
[2] 张俊, 李磊, 张峰, 等. X射线CT射束硬化校正方法综述[J]. CT理论与应用研究, 2013, 22(1): 195-204. ZHANG J, LI L, ZHANG F, et al. Review of the methods for beam hardening correction in X-ray computed tomography[J]. CT Theory and Applications, 2013, 22(1): 195-204. (in Chinese).
[3] 陈世杰, 马巍, 李国玉. 冻土CT图像中环形伪影和射线束硬化伪影的成因及校正方法[J]. 冰川冻土, 2020, 42(4): 1407-1416. CHEN S J, MA W, LI G Y. Ring artifact and beam hardening artifact in CT image of frozen soil: causes and correction methods[J]. Journal of Glaciology and Geocryology, 2020, 42(4): 1407-1416. (in Chinese).
[4] HOMOLKA P, FIGL M. Equivalent thicknesses of beam hardening filters consisting of aluminium, copper, Al/Cu and Al/Gold combinations and plumbiferous acrylic for 40 to 150 kVp diagnostic spectra[J]. Journal Of Radiological Protection, 2018, 38(4): 1269-1283. DOI: 10.1088/1361-6498/aadaf4.
[5] MCGUIRE A M, SMITH C D, CHAMBERLIN J H, et al. Reduction of beam hardening artifact in photon-counting computed tomography: Using low-energy threshold polyenergetic reconstruction[J]. Journal of Cardiovascular Computed Tomography, 2023, 17(5): 356-357. DOI: 10.1016/j.jcct.2023.08.007.
[6] 綦振国, 杨晨菲. 卷积神经网络在射线检测中的应用浅析[J]. 无损探伤, 2023, 47(4): 11-13+41. QI Z G, YANG C F. Application of convolutional neural network in radiographic testing[J]. Nondestructive Testing Technology, 2023, 47(4): 11-13+41. (in Chinese).
[7] 史再峰, 谢向桅, 曹清洁, 等. 基于卷积神经网络的能谱CT射束硬化伪影抑制方法[J]. 南开大学学报(自然科学版), 2021, 54(2): 74-79,84. SHI Z F, XIE X Z, CAO Q J, et al. Beam hardening artifacts reduction in spectral CT based on convolutional neural network[J]. Acta Scientiarum Naturalium Universitatis Nankaiensis, 2021, 54(2): 74-79,84. (in Chinese).
[8] KALARE K, BAJPAI M, SARKAR S. et al. Deep neural network for beam hardening artifacts removal in image reconstruction[J]. Applied Intelligence, 2022, 52: 6037-6056. DOI: 10.1007/s10489-021-02604-y.
[9] BRABANT L, PAUWELS E, DIERICK, M, et al. A novel beam hardening correction method requiring no prior knowledge, incorporated in an iterative reconstruction algorithm[J]. NDT& E International, 2012, 51: 68-73.
[10] XIU G Y, YUAN C Y, CHEN X H, et al. An innovative beam hardening correction method for computed tomography systems[J]. Traitement du Signal, 2019, 36(6): 515-520. DOI: 10.18280/ts.360606.
[11] 杨柳, 宋勇, 李廷, 等. 基于蒙卡模拟能谱划分的X射线图像硬化校正方法[J]. 中国医学物理学杂志, 2024, 41(9): 1145-1151. DOI: 10.3969/j.issn.1005-202X.2024.09.013. YANG L, SONG Y, LI T, et al. Beam hardening correction method based on Monte Carlo energy spectrum partitioning for X-ray image[J]. Chinese Journal of Medical Physics, 2024, 41(9): 1145-1151. DOI: 10.3969/j.issn.1005-202X.2024.09.013. (in Chinese).
[12] 王星艺, 孙跃文, 曾天辰, 等. 文物CT图像伪影校正方法研究[J]. 文物保护与考古科学, 2022, 34(6): 115-120. WANG X Y, SUN Y W, ZENG T C, et al. Research on the artifact correction of CT images of cultural relics[J]. Sciences of Conservation and Archaeology, 2022, 34(6): 115-120. (in Chinese).
[13] INGACHEVA, ANASTASIA, CHUKALINA, et al. Polychromatic CT Data Improvement with One-Parameter Power Correction[J]. Mathematical Problems in Engineering, 2019, 2019: 1405365. DOI: 10.1155/2019/1405365.
[14] 罗婷, 赵星, 赵云松, 等. X射线CT正交基材料分解成像方法及其在校正金属伪影中的应用[J]. 光学学报, 2024, 44(8): 45-59. LUO T, ZHAO X, ZHAO Y S, et al. Orthogonal multi-material decomposition for X-Ray CT and application in metal artifact correction[J]. Acta Optica Sinica, 2024, 44(8): 45-59. (in Chinese).
[15] WEI J T, CHEN P, HAN Y, et al. A multi-energy computed tomography method without image segmentation or prior knowledge of X-ray spectra or materials[J]. Heliyon, 2022, 8(11): e11584. DOI: 10.1016/j.heliyon.2022.e11584.
[16] WEI J T, CHEN P, HAN Y, et al. Blind separation model of multi-voltage projections for the hardening artifact correction in computed tomography[J]. Biomedical Signal Processing and Control, 2021, 64: 102236. DOI: 10.1016/j.bspc.2020.102236.
[17] SAJJA S, LEE Y, ERIKSSON M, et al. Technical principles of dual-energy cone beam computed tomography and clinical applications for radiation therapy[J]. Advances in Radiation Oncology, 2019, 5(1): 1-16.
[18] 王毅, 李勤, 江孝国, 等. 衰减透射法测量高能X射线能谱[J]. 强激光与粒子束, 2012, 24(6): 1466-1470. WANG Y, LI Q, JIANG X G, et al. 衰减透射法测量高能X射线能谱[J]. High Power Laser and Particle Beams, 2012, 24: 1466-1470. (in Chinese). [19] 王黎明, 赵英亮, 韩焱. 一种通用CT图像硬化校正算法研究[J]. 测试技术学报, 2005(1): 52-56. DOI: 10.3969/j.issn.1671-7449.2005.01.013. WANG L M, ZHAO Y L, HAN Y. A universal method of CT image hardening correction[J]. Journal of Test and Measurement Technology, 2005(1): 52-56. DOI: 10.3969/j.issn.1671-7449.2005.01.013. (in Chinese).
[20] DIWAKAR M, KUMAR M. A review on CT image noise and its denoising[J]. Biomedical Signal Processing and Control, 2018, 42: 73-88. DOI: 10.1016/j.bspc.2018.01.010.
[21] 林少春, 马建华, 陈武凡. 基于有序子集的惩罚加权最小二乘低剂量CT优质重建算法研究[J]. 中国组织工程研究与临床康复, 2008(4): 679-682. LIN S C, MA J H, CHEN W F. High quality reconstruction of low-dose computed tomography based on ordered subsets penalized weighted least-squares[J]. Journal of Clinical Rehabilitative Tissue Engineering Research, 2008(4): 679-682. (in Chinese).
[22] WEI Y, YU H, WANG G. Integral invariants for computed tomography[J]. IEEE Signal Processing Letters, 2006, 13(9): 549-552. DOI: 10.1109/LSP.2006.874452.
[23] TANG S, HUANG K, CHENG Y, et al. Optimization based beam-hardening correction in CT under data integral invariant constraint[J]. Physics in Medicine & Biology, 2018, 63(13): 135015.
计量
- 文章访问数: 2
- HTML全文浏览量: 5
- PDF下载量: 6
