一种低剂量X射线CT图像重建方法与流程

本发明涉及一种医学影像的图像处理技术，具体涉及一种基于统计建模与最大后验估计框架的低剂量X射线CT图像重建方法。

背景技术：

X射线CT扫描已经广泛应用于临床医学影像诊断，但是CT扫描过程中过高的X射线辐射剂量对人体存在潜在风险，容易造成辐射损伤、诱发恶性肿瘤等。在保证图像质量的前提下，最大限度降低X射线使用剂量已经成为医学CT成像领域研究的关键技术之一。

为了降低X射线辐射剂量，可以通过各种硬件技术及软件技术降低CT扫描中的X射线使用剂量。常见的技术包括降低管电流、降低X射线曝光时间以及减少投影数据的采集量(即稀疏角度CT扫描)等方法。

降低CT扫描中的管电流(Low-mA)和扫描时间可以直接减少使用X线的辐射剂量，但是其相应的成像数据中，随机成分的比率将大大增加，直接导致图像质量的严重退化，难以用于临床诊断。为了在保证成像质量的前提下，最大限度降低X射线辐射剂量，基于降低管电流和扫描时间的低剂量CT图像重建方法相继被提出，例如基于统计模型的迭代重建方法以及基于投影数据滤波的解析重建方法。其中，基于统计模型的迭代重建方法根据采集投影数据的统计特性以及成像系统构建CT图像重建模型，可以实现低剂量CT图像的优质重建；基于投影数据滤波的解析重建方法同样可以根据采集投影数据的统计特性以及成像系统进行数据滤波建模，再通过解析重建方法实现快速优质的低剂量CT图像重建。

基于统计模型的迭代重建方法需要对目标函数进行几十甚至上百次的反复迭代求解，这使得重建CT图像的时间代价往往非常庞大。重建相同像素大小的CT图像时，基于统计模型的迭代重建方法所需时间往往远超经典的解析重建方法，因此仍不能满足临床CT实时显像的需求。

而传统的基于投影数据滤波的解析重建方法并未充分利用数据内在统计信息，在投影数据恢复过程中往往导致图像原有细节信息的丢失，从而使得相应CT图像的分辨率下降。

针对现有技术不足，提供一种能够确保重建图像的分辨率的低剂量CT图像重建方法因此甚为必要。

技术实现要素：

本发明的目的在于提供一种所获重建CT图像分辨率更高，成像质量更好，细节信息保留更充分的低剂量X射线CT图像重建方法。

为达到上述目的，本发明采用的技术方案是：

步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据；

步骤S2:通过成像过程的统计规律构建生成投影数据的统计模型；

步骤S3:根据投影数据和弦图数据的结构特征与实际应用中的需求，构建数据先验的统计模型；

步骤S4:结合步骤S2与步骤S3，构建以投影数据信息为条件的弦图数据统计生成模型，并利用最大后验估计方法，构造弦图数据复原算法；

步骤S5:以步骤S1获得的投影数据为输入，应用步骤S4的弦图数据复原算法，获得复原弦图数据及其它统计变量；

若步骤S3中统计模型以重建的CT图像为其统计变量，那么步骤S5直接得到重建CT图像；

步骤S6:根据步骤S5所获的弦图数据进行CT图像重建，得到输出结果。

所述步骤S1中获取的CT设备的成像系统参数包括X射线入射光子强度I₀、系统电子噪声的方差

所述步骤S2中，通过成像过程的统计规律构建的投影数据统计生成模型为：

p＝I+ε,

其中，p为感受器上观测的原始投影数据，I为到达感受器的X射线光子强度，y为弦图数据，ε为系统电子噪声，它们的第i个分量分别代表第i个数据点上对应的数据；代表第i个数据点上电子噪声满足的分布，通常假设为一个以σ_ε为方差的正态分布，其形式为：

P{I_i|y_i}代表感第i个数据点上射线光子强度满足的条件分布，其概率密度函数为以下分布形式：

其中，表示以λ为均值的泊松分布，I₀为第i个数据点上的X射线入射光子强度。

所述步骤S2中，通过成像过程的统计规律构建的投影数据统计生成模型表达为如下的条件概率形式：

其中，P(I_i|y_i)与分别由公式(1)与(2)定义。

所述步骤S3中，根据投影数据和弦图数据的结构特征与实际应用中的需求，构建数据先验的统计模型；的具体形式应根据实际情况与对运算效率与运算精度的需求确定，可归纳为如下的表达形式：

其中，q是必要的辅助变量，其中，N是向量q的元素总数，L(q|0,b)是均值为0，尺度参数为b的拉普拉斯分布，P(b)是关于参数b的无信息先验，值恒为常数(即足够大范围内的均匀分布)，Ω_q是归一化常数，其值为Ω_q＝∫_f(y)＝qydy，f(y)是根据实际需要而确定的映射。

所述步骤S4中构建的统计生成模型为如下完整后验分布形式：

其中，P(I,p|y)与P(y,q,b)分别由公式(3)与(4)定义。

所述步骤S4中，根据最大后验估计方法，由统计模型转化得弦图数据复原模型为如下的优化：

s.t.f(y)＝q。

由于实际问题中f(·)往往为非线性映射，为计算方便起见，可通过变量替换q＝h(z)将(6)转化为如下便于求解的等价形式：

s.t.g(y)＝z,

其中，g(·)一般为线性映射，且h(z)满足h(g(y))＝f(y)。

所述步骤S5采用交替方向乘子法求解步骤S4中的弦图数据复原模型公式(7)，具体步骤包括：

S4.1)给出公式(7)的增广拉格朗日函数：

其中，Λ是乘子，μ是一个大于0的数；

S4.2)建立交替方向乘子法的迭代格式与终止条件：

迭代格式为:

Λ^k+1＝Λ^k+μ^k(g(y^k+1)-Z^k+1) (12)

μ^k+1＝ρμ^k (13)

其中，ρ是一个大于1的正数，一般设为ρ＝1.5，

迭代终止条件为：

S4.3)对问题(8)、(9)、(10)和(11)进行求解，给出迭代的具体算式；

S4.4)设置迭代的初始值设置为：y⁰＝lnI₀-lnρ，I⁰＝round(p)，z⁰＝g(y⁰)，Λ⁰＝0，其中，round(·)为取整函数，p为步骤S1所获得的投影数据；

S4.5)进行(8)-(13)的迭代运算，直到满足终止条件，若模型中以重建CT图像为其统计变量，直接获取最终CT图像重建结果；否则，需执行步骤S6，即以当前的y^k作为输出的弦图数据，并对该数据通过解析重建获取重建CT图像。

所述的(9)式即求解如下的问题：

其解I^k+1的任意分量满足

其中，式即：

因此，通过下面两步来求解问题(14)：

S9.1)将I^k中不满足(16)式第一个不等式的分量，以步长1不停下降，直到满足(16)式第一个不等式，作为I^k+1的对应分量；

S9.2)将I^k中不满足(16)式第二个不等式的分量，以步长1不停上升，直到满足(16)式第二个不等式，作为I^k+1的对应分量。

所述(10)式即求解如下的问题：

根据函数h(·)的形式，一般有对应的阈值算子形式解析解：

其中，是对应阈值算子。

所述的(11)式即求解如下的问题：

其解为：

所述的(8)式即求解如下的问题：

针对不同g(·)形式，用FISTA算法，直接调用求解。

所述的步骤S6，采用滤波反投影算法，对步骤S5输出的弦图数据进行CT图像重建。

本发明相比传统方法，更加充分准确的对投影数据及其弦图数据的统计信息进行统计建模，从而所获重建CT图像分辨率更高，成像质量更好，细节信息保留更充分。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1为本发明的流程图。

图2为实例1中仿真使用的假人体模数据，右下角红框中的图像为原图红框中图像提高对比度并放大四倍的结果。

图3为(a)实例1中使用的低剂量CT投影数据重建的CT图像结果及(b)利用实例1方法重建的CT图像结果，图像右下角红框中的图像为原图红框中图像提高对比度并放大四倍的结果。

图4为(a)实例1中使用的不含噪弦图数据图像；(b)低剂量情形下(含噪)弦图数据图像；(c)实例1方法估计的弦图数据图像；(d)(c)与(b)的差值绝对值图像(注所有弦图只展示左半部分，右半部分与左半部分相似)。

图5为弦图数据的近似分片平面先验展示。左上角为灰度图展示，右上角为红框处图像的放大图。可以看出，弦图数据可以较好地近似若干平面的拼合。

图6为CT图像数据的平滑先验展示。左上角为灰度图展示，右上角为红框处图像的放大图。可以看出，CT图像数据可以较好地近似若干水平平面的拼合。

具体实施方式

结合以下实例对本发明作进一步描述。

实施例1：

采用如图2所示的假人体模图像作为本发明的计算机仿真实验对象。体模图像大小设为512×512，模拟CT设备的X射线源到旋转中心和探测器的距离分别为1361.2mm和615.18mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.85mm，弦图数据的片段如图4所示。选择弦图的分片平面近似特征(如图5所示)为先验，本发明依次包括如下步骤:

步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据p；所获取的CT设备的成像系统参数包括X射线入射光子强度I₀、系统电子噪声的方差σ_ε。

步骤S2:通过成像过程的统计规律构建数据生成的统计模型；表达为如下的条件概率分布：

其中，p为感受器上观测的原始投影数据，I为到达感受器的X射线光子强度，ε为系统电子噪声，它们的第i个分量分别代表第i个数据点上对应的数据，P(I_i|y_i)与分别由(1)与(2)定义。

步骤S3:根据弦图近似分片平面的结构特征，如图5所示，可知投影数据的二阶差商场是稀疏的。故可以假设其二阶差商经过f(·)＝ln(|·|+∈)-ln(∈)变换后满足拉普拉斯分布，因此，我们可以把先验分布表达为下式：

其中∈是一个微小的数，取值为10^-16，D₂是二阶差商矩阵，q是必要的辅助变量，其中，N是向量q的元素总数，L(q|0,b)是均值为0，尺度参数为b的拉普拉斯分布，P(b)是关于参数b的无信息先验，值恒为常数(即足够大范围内的均匀分布)，Ω_q是归一化常数，其值为Ω_q＝∫_f(y)＝qydy，f(y)是根据实际需要而确定的映射。

步骤S4:结合步骤S2与步骤S3，构建如下后验分布统计模型：

其中，P(I,p|y)与P(y,q|b)分别由(3)与(4)定义。

根据最大后验估计方法，可得如下弦图数据复原优化问题：

R.t.ln(|D₂y|+∈)-ln(∈)＝q，

其中，N是向量q的元素总数。以上优化问题可等价转换为：

R.t.D₂y＝z，

可以通过如下的步骤对模型进行求解。

步骤S5:基于步骤S1所输入投影数据及所设模型参数，应用步骤S4的弦图复原模型构造弦图复原求解算法。

求解的过程利用如下的迭代格式:

Λ^k+1＝Λ^k+μ^k(D₂y^k+1-z^k+1) (12)

μ^k+1＝ρμ^k (13)

其中，ρ是一个大于1的正数，一般取值ρ＝1.5。以下给出迭代的细节：

A.(8)式即求解如下的问题：

可直接调用如下FISTA算法进行求解：

输入：

步1.设置α为I₀的最大元素值的10倍，设置最大迭代步数为T＝6,命a⁽¹⁾＝1,t＝1,

步2.当t≤T,重复迭代以下操作：

a)

b)

c)

d)t＝t+1

步3.输出

B.(9)式即求解如下的问题：

可通过以下对问题求解：

(S9.1)将I^k中不满足(16)式第一个不等式的分量，以步长1不停下降，直到满足(35)式第一个不等式，作为的对应分量。

(S9.2)将中不满足(16)式第二个不等式的分量，以步长1不停上升，直到满足(16)式第二个不等式，作为I^k+1的对应分量。

C.(10)式即求解如下的问题：

根据函数h(·)的形式，一般有对应的阈值算子形式的解析解：

阈值算子定义如下:

其中，c₁＝β-∈，

D.(11)式的解为：

我们设置迭代的初始值为：y⁰＝lnI₀-lnp，I⁰＝round(p)，z⁰＝D₂y⁰，Λ⁰＝0，其中，round(·)为取整运算，p为步骤S1所获得的投影数据。迭代终止条件为：

步骤S6:根据步骤S5得到的复原弦图数据进行CT图像重建，得到最终输出的CT图像。

实施例2.

采用如图2所示的假人体模图像(弦图数)作为本发明的计算机仿真实验对象。体模图像大小设为512×512，模拟CT设备的X射线源到旋转中心和探测器的距离分别为1361.2mm和615.18mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.85mm，弦图数据的片段如图4所示。选CT图像的平滑先验(如图6所示)为先验，本发明依次包括如下步骤:

步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据p；所获取的CT设备的成像系统参数包括X射线入射光子强度I₀、系统电子噪声的方差σ_ε。

步骤S2:通过成像过程的统计规律构建数据生成的统计模型；表达为如下的条件概率分布：

其中，p为感受器上观测的原始投影数据，I为到达感受器的X射线光子强度，ε为系统电子噪声，它们的第i个分量分别代表第i个数据点上对应的数据，P(I_i|y_i)与分别由(1)与(2)定义。

步骤S3:根据CT图像的平滑特征(可由若干水平平面近似拟合)，如图6所示，可知CT图像的一阶差商场是稀疏的，故可以假设其一阶差商经过f(·)＝ln(|·|+∈)-ln(∈)变换后是满足拉普拉斯分布的。因此，我们可以把先验分布表达为下式：

其中∈是一个微小的数，取值为10^-16，A是成像的系统矩阵，D₁是一阶差商矩阵。

步骤S4:结合步骤S2与步骤S3，构建完整的统计模型，如下的后验分布形式：

其中，P(I,p|y)与P(y,q,b)分别由(3)与(4)定义。

根据最大后验估计方法，由统计模型转化的弦图数据复原模型可表达为如下的优化问题：

s.t.ln(|D₁A^-1y|+∈)-ln(∈)＝q，

其中，N是向量q的元素总数。该优化问题可以等价转换为如下的问题：

s.t.D₁A^-1y＝z。

步骤S5:基于步骤S1所输入投影数据及其所设模型参数，对步骤S4中弦图复原模型进行求解。

求解过程采用如下迭代格式:

Λ^k+1＝Λ^k+μ^k(g(y^k+1)-z^k+1) (12)

μ^k+1＝ρμ^k (13)

其中，ρ是一个大于1的正数，一般取值ρ＝1.5。以下给出迭代细节：

A.(8)式即求解如下的问题：

可等价转化为：

s.t.Ax＝y。

该优化问题可以用ADMM算法来求解，其增广朗日函数为：

该问题可通过如下迭代来求解

Γ^(t+1)＝Γ^k+τ^k(Ax^(t+1)-y^(t+1)) (70)

τ^(t+1)＝ρτ^(t+1) (71)

其中，(68)即求解如下的问题：

该问题可直接调用如下FISTA算法进行求解：

输入：

步1.设置α为I₀的最大元素值的10倍，设置最大迭代步数为T＝6,命a⁽¹⁾＝1,l＝1,

步2.当l≤T,重复迭代以下操作：

a)

b)

c)

d)l＝l+1

步3.输出

输入：

步1.设置α为I₀的最大元素值的10倍，设置最大迭代步数为T＝6,命a⁽¹⁾＝1,l＝1,

步2.当l≤T,重复迭代以下操作：

a)

b)

c)

d)l＝l+1

e)步3.输出

该迭代算法以τ⁽⁰⁾＝μ^k,y⁽⁰⁾＝y^k和x⁽⁰⁾＝x^k为初始,输出y^k+1＝y^(t+1),x^k+1＝x^(t+1)。我们一般采用的输出为：x^k+1＝A^-1y^k+1，可避免其它迭代步中A^-1y^k+1的运算，将其作为最终输出的重建CT图像结果。

迭代终止条件可设为：

经验上，设置最大迭代步数为5即可保证实验效果。

B.(9)式即求解如下的问题：

可以通过下面两步来求解问题(14)

(S9.1)将I^k中不满足(16)式第一个不等式的分量，以步长1不停下降，直到满足(35)式第一个不等式，作为的对应分量。

(S9.2)将中不满足(16)式第二个不等式的分量，以步长1不停上升，直到满足(16)式第二个不等式，作为I^k+1的对应分量。

C.(10)式即求解如下的问题：

根据函数h(·)的形式，一般有对应阈值算子形式的解析解：

阈值算子定义与(18)式相同。

D.(11)式存在解析解：

我们设置迭代的初始值设置为：y⁰＝lnI₀-lnp，I⁰＝round(p)，z⁰＝g(y⁰)，Λ⁰＝0。其中，round(·)为取整运算，p是步骤S1所获得的投影数据。迭代终止条件为：

输出为x^k+1作为最终CT图像重建结果。