CN105701847A - 一种改进权系数矩阵的代数重建方法 - Google Patents

Abstract

本发明涉及一种改进权系数矩阵的代数重建方法，属于CT投影数据重建方法技术领域。该方法主要为解决Ray-Box Intersection权系数矩阵算法存在的逻辑漏洞以及基于双线性插值的SART+TVM重建图像细节模糊的问题，通过对Ray-Box Intersection权系数矩阵算法的改进，缩短了重建时间，保持了图像细节，提高了信噪比；同时利用TVM算法作为约束项，加快了迭代重建收敛速度，减少了重建图像的伪影和噪声，有效提高了重建图像质量。

Description

一种改进权系数矩阵的代数重建方法

技术领域

本发明属于CT投影数据重建方法技术领域，涉及一种基于改进的Ray-BoxIntersection权系数矩阵的联合代数重建方法。

背景技术

计算机层析成像(ComputedTomography,CT)技术是一种重要的无损检测技术，其核心是图像重建技术，即利用多个采样视角下的投影数据重建出物体内、外部结构特征的二维或三维图像。CT重建算法主要包括解析重建算法和迭代重建算法。解析重建需要在180°或360°范围内均匀采样足够的投影数据，但在工业CT***实际应用中，常出现由于扫描工件太长使得采样角度有限，或者X射线无法穿透高密度区域导致采样计数不足，或者减少采样密度以提高扫描速度等情况，这些情况均属于不完备投影数据重建问题，无法满足解析重建的条件，需采用迭代重建算法。

近几年，解决不完备投影数据重建问题，一般采用基于图像全变差的迭代重建算法，尤其是基于图像全变差最小化的联合代数重建算法(SimultaneousAlgebraicReconstructionTechniqueBasedTotalVariationMinimization,SART+TVM)，该算法利用CT图像梯度变换的稀疏性，并与SART算法相结合，在每次迭代过程中，先利用SART算法获得初步的重建图像，然后利用TVM算法沿图像梯度方向调整图像全变差。SART+TVM算法重建图像质量有多个影响因素，包括：SART迭代次数、SART松弛因子、权系数矩阵、TVM迭代次数、TVM调节因子等，而权系数矩阵是决定SART+TVM算法重建图像质量的关键因素之一。

双线性插值算法是当前最常用的求解权系数矩阵的方法，其利用相邻像素点进行比例插值求解当前点的像素值，能得到比较完整的结构信息，但图像边缘易模糊，相比之下，Ray-BoxIntersection算法在原理上更接近真实投影过程，但因其存在逻辑漏洞，例如：射线与某像素块不相交，却判断为相交，或权重值大于射线穿过像素块的实际长度值，导致重建图像质量差，伪影多。因此，需要对Ray-BoxIntersection算法进行改进，弥补其自身存在的逻辑漏洞，提高重建图像质量，尤其是重建图像细节要突出。

发明内容

有鉴于此，本发明的目的在于提供一种改进权系数矩阵的代数重建方法，该方法针对Ray-BoxIntersection权系数矩阵算法存在的逻辑漏洞以及基于双线性插值的SART+TVM重建图像细节模糊的问题，对Ray-BoxIntersection权系数矩阵算法进行改进，从而缩短重建时间，保持图像细节，提高信噪比；同时利用TVM算法作为约束项，加快了迭代重建收敛速度，减少了重建图像的伪影和噪声，有效提高了重建图像质量。

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

一种改进权系数矩阵的代数重建方法，包括以下步骤：

S1：获得投影数据p_i，初始化CT扫描参数，i＝0,1,2,...,N，N为投影视角总数；

S2：将待重建图像x_j赋初值，j＝0,1,2,...,M-1，M为图像像素总数，k为SART迭代次数；

S3：依据射线驱动方式，计算该投影方向下的权系数矩阵A＝{a_ij}，a_ij为第i条射线对第j个像素块的加权值；

S4：正投影，获取第i条射线的模拟投影值

S5：依据射线的实测投影值p_i、模拟投影值和权系数矩阵A，求出第i条射线的修正值 ${\mathrm{\Δ}}_i=(p_i-{\tilde{P}}_i)/\underset{m=1}{\overset{M}{\Σ}}{a}_{im}{x}_m^{\left(k\right)};$

S6：i＝i+1，重复步骤S3-S5，直至完成该投影方向下的所有射线的修正值时，根据以下SART迭代公式进行反投影更新图像得到

$x_j^{(k+1)}=x_j^k+{\mathrm{\λ}}^k\frac{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}\[a_{ij}{\mathrm{\Δ}}_i\]}{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}{a}_{ij}}$

式中和分别为第k+1和第k次子迭代过程中第j个像素块的像素值，I_θ为角度θ下的所有射线的集合，λ^k为第k次子迭代过程中松弛因子；

S7：重复步骤S2-S6，直至完成所有投影角度的修正；

S8：对更新后图像进行梯度下降法调整图像全变差，图像全变差公式为：

$\frac{u_{n+1}-u_n}{\mathrm{\Δ}t}=\▿\·\frac{\▿u_n}{\left|\right|\▿u_n\left|\right|}-\mathrm{\α}(u_n-x_j^{(k+1)})$

1≤n＜N_TVM,N_TVM为TVM的迭代次数，u_n为全变差最小化图像，▽u_n代表图像u_n的梯度，||·||代表L1范数算子，α为TVM调解因子；

S9：令判断是否达到收敛条件ε为大于0的极小值，是则转向步骤S11，不是则转向步骤S10；

S10：重设迭代图像重复步骤S3-S9；

S11：退出循环，得到重建图像。

进一步，所述步骤S3具体包括：

S31：将射线束等角度离散化，利用射线与探测器一一对应的关系，可得到每条射线的位置信息；

S32：利用射线和图像区域的代数关系，得到射线在图像区域内的起点(x₁,y₁)和终点(x₂,y₂)；

S33：求取射线在重建区域内经过的所有像素块的左下角坐标集合{(Indx,Indy)}；

S34：计算射线与所经过的像素块的交点；

S35：计算同一像素块上两交点之间的距离，并将其当作该像素块对于该条射线的权重，特殊地，当射线水平或者垂直于像素块时，权重为1；

S36：重复步骤S32-S35，直至遍历所有射线。

本发明的有益效果在于：本发明提供的方法缩短了重建时间，保持了图像细节，提高了信噪比；同时利用TVM算法作为约束项，加快了迭代重建收敛速度，减少了重建图像的伪影和噪声，有效提高了重建图像质量。

附图说明

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

图1为SART+TVM重建算法流程图；

图2为改进前算法流程图；

图3为改进后算法流程图；

图4为射线与重建区域示意图。

具体实施方式

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

图1为SART+TVM重建算法流程图，图2为改进前算法流程图，图3为改进后算法流程图，图4为射线与重建区域示意图，如图所示，本发明主要针对不完备投影数据重建问题，采用SART+TVM重建算法，即在SART迭代重建算法的基础上引入图像全变差最小化约束。

本方法具体包括以下步骤：

S1：获得投影数据p_i，初始化CT扫描参数，i＝0,1,2,...,N，N为投影视角总数；

S2：将待重建图像x_j赋初值，j＝0,1,2,...,M-1，M为图像像素总数，k为SART迭代次数；

S3：依据射线驱动方式，计算该投影方向下的权系数矩阵A＝{a_ij}，a_ij为第i条射线对第j个像素块的加权值；

S4：正投影，获取第i条射线的模拟投影值

S5：依据射线的实测投影值p_i、模拟投影值和权系数矩阵A，求出第i条射线的修正值 ${\mathrm{\Δ}}_i=(p_i-{\tilde{P}}_i)/\underset{m=1}{\overset{M}{\Σ}}{a}_{im}{x}_m^{\left(k\right)};$

S6：i＝i+1，重复步骤S3-S5，直至完成该投影方向下的所有射线的修正值时，根据以下SART迭代公式进行反投影更新图像得到

$x_j^{(k+1)}=x_j^k+{\mathrm{\λ}}^k\frac{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}\[a_{ij}{\mathrm{\Δ}}_i\]}{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}{a}_{ij}}$

式中和分别为第k+1和第k次子迭代过程中第j个像素块的像素值，I_θ为角度θ下的所有射线的集合，λ^k为第k次子迭代过程中松弛因子；

S7：重复步骤S2-S6，直至完成所有投影角度的修正；

S8：对更新后图像进行梯度下降法调整图像全变差，图像全变差公式为：

$\frac{u_{n+1}-u_n}{\mathrm{\Δ}t}=\▿\·\frac{\▿u_n}{\left|\right|\▿u_n\left|\right|}-\mathrm{\α}(u_n-x_j^{(k+1)})$

1≤n＜N_TVM,N_TVM为TVM的迭代次数，u_n为全变差最小化图像，▽u_n代表图像u_n的梯度，||·||代表L1范数算子，α为TVM调解因子；

S9：令判断是否达到收敛条件ε为大于0的极小值，是则转向步骤S11，不是则转向步骤S10；

S10：重设迭代图像重复步骤S3-S9；

S11：退出循环，得到重建图像。

上面步骤S3中权系数矩阵算法采用改进的Ray-BoxIntersection算法，其具体实施步骤为：

S31：将射线束等角度离散化，利用射线与探测器一一对应的关系，可得到每条射线的位置信息；

S32：利用射线和图像区域的代数关系，得到射线在图像区域内的起点(x₁,y₁)和终点(x₂,y₂)；

S33：求取射线在重建区域内经过的所有像素块的左下角坐标集合{(Indx,Indy)}；

S34：计算射线与所经过的像素块的交点；

S35：计算同一像素块上两交点之间的距离，并将其当作该像素块对于该条射线的权重，特殊地，当射线水平或者垂直于像素块时，权重为1；

S36：重复步骤S32-S35，直至遍历所有射线。

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

Claims (2)

1.一种改进权系数矩阵的代数重建方法，其特征在于：包括以下步骤：

S1：获得投影数据p_i，初始化CT扫描参数，i＝0,1,2,...,N，N为投影视角总数；

S2：将待重建图像x_j赋初值，j＝0,1,2,...,M-1，M为图像像素总数，k为SART迭代次数；

S3：依据射线驱动方式，计算该投影方向下的权系数矩阵A＝{a_ij}，a_ij为第i条射线对第j个像素块的加权值；

S4：正投影，获取第i条射线的模拟投影值

S5：依据射线的实测投影值p_i、模拟投影值和权系数矩阵A，求出第i条射线的修正值 ${\mathrm{\Δ}}_i=(p_i-\widetilde{P_i})/\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{a}_{\mathrm{im}}{x}_m^{\left(k\right)};$

S6：i＝i+1，重复步骤S3-S5，直至完成该投影方向下的所有射线的修正值时，根据以下SART迭代公式进行反投影更新图像得到

$x_j^{(k+1)}=x_j^k+{\mathrm{\λ}}^k\frac{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}\[a_{ij}{\mathrm{\Δ}}_i\]}{\underset{i\⋐I_{\mathrm{\θ}}}{\Σ}{a}_{ij}}$

式中和分别为第k+1和第k次子迭代过程中第j个像素块的像素值，I_θ为角度θ下的所有射线的集合，λ^k为第k次子迭代过程中松弛因子；

S7：重复步骤S2-S6，直至完成所有投影角度的修正；

S8：对更新后图像进行梯度下降法调整图像全变差，图像全变差公式为：

$\frac{u_{n+1}-u_n}{\mathrm{\Δ}t}=\▿\·\frac{\▿u_n}{\left|\right|\▿u_n\left|\right|}-\mathrm{\α}(u_n-x_j^{(k+1)})$

1≤n＜N_TVM,N_TVM为TVM的迭代次数，u_n为全变差最小化图像，代表图像u_n的梯度，||·||代表L1范数算子，α为TVM调解因子；

S9：令判断是否达到收敛条件ε为大于0的极小值，是则转向步骤S11，不是则转向步骤S10；

S10：重设迭代图像重复步骤S3-S9；

S11：退出循环，得到重建图像。

2.根据权利要求1所述的一种改进权系数矩阵的代数重建方法，其特征在于：所述步骤S3具体包括：

S31：将射线束等角度离散化，利用射线与探测器一一对应的关系，可得到每条射线的位置信息；

S32：利用射线和图像区域的代数关系，得到射线在图像区域内的起点(x₁,y₁)和终点(x₂,y₂)；

S33：求取射线在重建区域内经过的所有像素块的左下角坐标集合{(Indx,Indy)}；

S34：计算射线与所经过的像素块的交点；

S35：计算同一像素块上两交点之间的距离，并将其当作该像素块对于该条射线的权重，特殊地，当射线水平或者垂直于像素块时，权重为1；

S36：重复步骤S32-S35，直至遍历所有射线。