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

技术特征：

1.一种低剂量X射线CT图像重建方法，其特征在于包括如下步骤：

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

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

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

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

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

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

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

2.根据权利要求1所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S1中获取的CT设备的成像系统参数包括X射线入射光子强度I₀、系统电子噪声的方差

3.根据权利要求1所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S2中，通过成像过程的统计规律构建的投影数据统计生成模型为：

p＝I+ε,

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

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

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

4.根据权利要求3所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S2中，通过成像过程的统计规律构建的投影数据统计生成模型表达为如下的条件概率形式：

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

5.根据权利要求1所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S3中，根据投影数据和弦图数据的结构特征与实际应用中的需求，构建数据先验的统计模型；的具体形式应根据实际情况与对运算效率与运算精度的需求确定，可归纳为如下的表达形式：

$(y,q,b)~P(y,q,b)=P\left(y\right|q\left)L\right(q|0,b)P\left(b\right)=\frac{\mathrm{\δ}\left(f\right(y)=q)}{{\mathrm{\Ω}}_q}.\frac{1}{b^N}{e}^{-\frac{\left|\right|q|{|}_1}{b}}P\left(b\right),---\left(4\right)$

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

6.根据权利要求1至5任一项所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S4中构建的统计生成模型为如下完整后验分布形式：

$P(y,I,q,b|p)=\frac{P(I,p|y)P(y,q,b)}{P\left(p\right)},---\left(5\right)$

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

7.根据权利要求6所述的低剂量X射线CT图像重建方法，其特征在于：

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

$\begin{array}{c}\underset{y,I,q,b}{min}\underset{i}{\Σ}(\frac{{(p_i-I_i)}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}-I_i(ln{I}_0^i)+I_i{y}_i+ln{I}_i!+I_0{e}^{-y_i})+\frac{1}{b}\left|\right|q|{|}_1+Nln\left(b\right)\\ \begin{array}{cc}s.t.& f\left(y\right)=q\end{array}\end{array}---\left(6\right)$

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

$\begin{array}{c}\underset{y,I,z,b}{\min }\underset{i}{\Σ}(\frac{{(p_i-I_i)}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}-I_i(\ln I_0^i)+I_i{y}_i+\ln I_i!+I_0{e}^{-y_i})+\frac{1}{b}\left|\right|h\left(z\right)|{|}_1+N\ln \left(b\right)\\ \begin{array}{cc}s.t.& g\left(y\right)=z,\end{array}\end{array}---\left(7\right)$

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

8.根据权利要求7所述的低剂量X射线CT图像重建方法，其特征在于：所述步骤S5采用交替方向乘子法求解步骤S4中的弦图数据复原模型公式(7)，具体步骤包括：

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

$\begin{array}{c}{L}_{\mathrm{\μ}}(y,I,z,b,\mathrm{\Λ})\\ =\underset{i}{\Σ}(\frac{(p_i-I_i)}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}-I_i(\ln I_0^i)+I_i{y}_i+\ln I_i!+I_0{e}^{-y_i})+\frac{1}{b}\left|\right|h\left(z\right)|{|}_1\\ +N\ln \left(b\right)+{\mathrm{\Λ}}^T\left(g\right(y)-z)+\frac{\mathrm{\μ}}{2}\left|\right|g\left(y\right)-z|{|}_2^2,\end{array}$

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

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

迭代格式为:

$y^{k+1}=\arg \underset{y}{min}{L}_{{\mathrm{\μ}}^k}(y,I^k,z^k,b^k,{\mathrm{\Λ}}^k)---\left(8\right)$

$I^{k+1}=\arg \underset{I}{min}{L}_{{\mathrm{\μ}}^k}(y^{k+1},I,z^k,b^k,{\mathrm{\Λ}}^k)---\left(9\right)$

$z^{k+1}=\arg \underset{z}{min}{L}_{{\mathrm{\μ}}^k}(y^{k+1},I^{k+1},z,b^k,{\mathrm{\Λ}}^k)---\left(10\right)$

$b^{k+1}=\arg \underset{b}{min}{L}_{{\mathrm{\μ}}^k}(y^{k+1},I^{k+1},z^{k+1},b,{\mathrm{\Λ}}^k)---\left(11\right)$

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

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

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

迭代终止条件为：

$\frac{\left|\right|g\left(y^{k+1}\right)-z^{k+1}|{|}_2^2}{\left|\right|g\left(y^{k+1}\right)|{|}_2^2}\≤{10}^{-4};$

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

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

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

9.根据权利要求8所述的低剂量X射线CT图像重建方法，其特征在于：所述的(9)式即求解如下的问题：

$\underset{I}{min}\underset{i}{\Σ}(\frac{{(p_i-I_i)}^2}{2{\mathrm{\σ}}_{\mathrm{\ϵ}}^2}-I_i(ln{I}_0^i)+I_i{y}_i^{k+1}+ln{I}_i!),---\left(14\right)$

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

$\left\{\begin{array}{c}F(I_i^{k+1}-1)\≥F\left(I_i^{k+1}\right)\\ F(I_i^{k+1}+1)\≥F\left(I_i^{k+1}\right)\end{array}\right.,---\left(15\right)$

其中，(15)式即：

$\left\{\begin{array}{c}\frac{(2p_i-2I_i+1)}{2{\mathrm{\σ}}^2}-(y_i^{k+1}-\ln I_0^i)-\ln I_i\≥0\\ (y_i^k-\ln I_0^i)+\ln (I_i+1)-\frac{(2p_i-2I_i-1)}{2{\mathrm{\σ}}^2}\≥0\end{array}\right.,---\left(16\right)$

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

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

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

10.根据权利要求8所述的低剂量X射线CT图像重建方法，其特征在于：所述(10)式即求解如下的问题：

$\underset{z}{min}\frac{1}{{\mathrm{b\μ}}^k}\left|\right|h\left(z\right)|{|}_1+\frac{1}{2}||g\left(y^{k+1}\right)+{\mathrm{\μ}}^{-1}{\mathrm{\Λ}}^k-z|{|}_2^2,---\left(17\right)$

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

$z^{k+1}=S_{\frac{1}{{\mathrm{b\μ}}^k}}\left(g\right(y^{k+1})+{\mathrm{\μ}}^{-1}{\mathrm{\Λ}}^k),---\left(18\right)$

其中，是对应阈值算子。

11.根据权利要求8所述的低剂量X射线CT图像重建方法，其特征在于：所述的(11)式即求解如下的问题：

$\underset{b}{min}\frac{1}{b}\left|\right|h\left(z^{k+1}\right)|{|}_1+Nln\left(b\right),$

其解为：

$b^{k+1}=\frac{\left|\right|h\left(z^{k+1}\right)|{|}_1}{N}---\left(19\right)$

12.根据权利要求8所述的低剂量X射线CT图像重建方法，其特征在于：所述的(8)式即求解如下的问题：

$\underset{y}{min}\underset{i}{\Σ}(I_i^k{y}_i+I_0{e}^{-y_i})+\frac{\mathrm{\μ}}{2}\left|\right|g\left(y^{k+1}\right)+{\mathrm{\μ}}^{-1}{\mathrm{\Λ}}^k-z^k|{|}_2^2$

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

13.根据权利要求1所述的低剂量X射线CT图像重建方法，其特征在于：所述的步骤S6，采用滤波反投影算法，对步骤S5输出的弦图数据进行CT图像重建。