CN103810735A - 一种低剂量x射线ct图像统计迭代重建方法 - Google Patents

Abstract

本发明公开一种低剂量X射线CT图像统计迭代重建方法，包括步骤：对CT设备的低剂量X射线CT图像的投影数据y_raw进行图像重建，得到迭代初始图像μ^init；对投影数据y_raw进行数据恢复处理，得到恢复后的投影数据y_restored，对恢复后的投影数据y_restored进行图像重建，得到参考图像μ^ref；以参考图像μ^ref和迭代初始图像_μ ^init，按R(μ^init)=φ(μ^init-SRNLM(μ^init))构建边缘保持先验R(μ^init)，其中，φ()为能量势函数，SRNLM(μ^init)为由参考图像μ^ref引导的非局部均值滤波；根据迭代初始图像μ^init边缘保持先验R(μ^init)，利用统计迭代公式进行迭代计算，得到迭代重建图像μ^iter；当对重建图像μ^iter进行迭代计算的结果满足迭代终止条件时，迭代终止，获得低剂量X射线CT图像最终的重建图像。本发明能够有效地去除图像噪声并抑制条形伪影，同时较好地保持图像细节信息。

Description

一种低剂量X射线CT图像统计迭代重建方法

技术领域

本发明涉及医学影像的计算机处理技术，尤其是涉及一种低剂量X射线CT图像统计迭代重建方法。

背景技术

X射线CT扫描已经广泛应用于临床医学影像诊断，但是CT扫描过程中过高的X射线辐射剂量会存在致癌风险。为了降低对使用者的损害，如何最大限度地降低X射线使用剂量已经成为医学CT成像领域研究的关键技术之一。

为了降低X射线辐射剂量，现有技术中使用的最简便的途径就是降低CT扫描过程中的管电流和扫描时间。在此基础上，现有技术针对低剂量X射线CT图像重建主要包括滤波反投影方法及统计迭代重建方法。

1、滤波反投影方法（Filtered Back-Projection，FBP）。

针对低剂量X射线CT图像，目前临床设备中常用的图像重建技术仍为滤波反投影方法。滤波反投影方法是通过对采集的投影数据进行傅立叶变换并采用一定的滤波函数（一般为斜坡滤波）进行滤波，然后对处理后的CT图像逐个角度进行反投影，得到重建图像。

但是，由于降低了管电流和扫描时间，低剂量CT投影数据中含有大量的噪声，基于滤波反投影方法重建的图像质量存在严重的退化现象，图像中存在大量的噪声和条形伪影，难以满足临床诊断需要。

2、基于统计模型的迭代重建方法，也称为统计迭代重建方法。

为了在保证图像质量的前提下大幅降低X射线辐射剂量，诸多基于降低管电流和扫描时间的低剂量CT图像重建方法相继提出。其中，基于统计模型的迭代重建方法，也称为统计迭代重建方法，通过对采集的投影数据的噪声以及成像系统进行图像重建模型构建，可以实现低剂量CT图像优质重建。为了保证统计迭代重建方法的鲁棒性，通常会引入一定的先验信息进行约束。

但是，现有技术使用基于边缘保持先验的统计迭代重建方法，由于低剂量CT重建图像中条形伪影不同于一般的图像噪声，常常把条形伪影作为一种图像结构进行保持，从而使得低剂量CT重建图像中条形伪影依然存在，不能起到抑制条形伪影的作用。因此，现有技术中基于一般边缘保持先验的统计迭代重建方法，虽然可以有效地去除图像噪声，却不能抑制图像中的条形伪影。

发明内容

针对现有技术的不足，本发明提供一种基于改进的边缘保持先验的低剂量X射线CT图像统计迭代重建方法，能够针对降低管电流和扫描时间的低剂量CT图像进行重建，能够有效地去除图像噪声并抑制条形伪影，同时较好地保持图像细节信息。

本发明采用如下技术方案实现：一种低剂量X射线CT图像统计迭代重建方法，其包括步骤：

获取CT设备的系统参数和低剂量X射线CT图像的投影数据y_raw，对投影数据y_raw进行图像重建，得到迭代初始图像μ^init；

对投影数据y_raw进行数据恢复处理，得到恢复后的投影数据y_restored，对恢复后的投影数据y_restored进行图像重建，得到参考图像μ^ref；

以参考图像μ^ref和迭代初始图像μ^init，按R(μ^init)=φ(μ^init-SRNLM(μ^init))构建边缘保持先验R(μ^init)，其中，φ()为能量势函数，SRNLM(μ^init)为由参考图像μ^ref引导的非局部均值滤波；

根据迭代初始图像μ^init边缘保持先验R(μ^init)，利用统计迭代公式进行迭代计算，得到迭代重建图像μ^iter；

当对重建图像μ^iter进行迭代计算的结果满足迭代终止条件时，迭代终止，获得低剂量X射线CT图像最终的重建图像。

其中，对投影数据y_raw进行数据恢复处理的步骤，是通过迭代公式

进行迭代求解，得到恢复后的投影数据y_restored；

其中，y表示待恢复的投影数据；Σ为一对角矩阵，其对角线元素为对应投影数据y_raw的方差R(y)为先验约束项；λ>0为正则化调节参数。

其中，投影数据y_raw的方差

通过投影数据y_raw的逐个数据点估计得到，估计公式为

其中I₀为投影数据y_raw的第i个数据点的X射线入射光子强度，

为CT设备的系统电子噪声的方差。

其中，所述迭代公式使用梯度下降法或者高斯-塞德尔优化算法进行迭代求解获得恢复后的投影数据y_restored。

其中，对投影数据y_raw及对恢复后的投影数据y_restored均是采用滤波反投影方法进行图像重建。

其中，由参考图像μ^ref引导的非局部均值滤波SRNLM(μ^init)的计算公式为：

$\mathrm{SRNLM}\left({\mathrm{\μ}}_i^{\mathrm{init}}\right)=\frac{1}{Z\left(i\right)}{\mathrm{\Σ}}_{j\∈N_i}\exp \{-{\left|\right|{\mathrm{\μ}}^{\mathrm{init}}\left(V_i\right)-{\mathrm{\μ}}^{\mathrm{ref}}\left(V_j\right)\left|\right|}_{2,a}^2/h^2\}{\mathrm{\μ}}_j^{\mathrm{init}}$

其中，

为归一化因子，V_i和V_j分别为图像μ^init中像素点i处和图像μ^ref中像素点j处的两个局部相似窗，μ^init(V_i)和μ^ref(V_j)分别为迭代初始图像μ^init在局部相似窗V_i内的图像矩阵块和参考图像μ^ref在局部相似窗V_j内的图像矩阵块，N_i为以像素点i为中心的搜索窗，符号||||₂表示计算两个相似窗内的图像矩阵块之间的欧氏距离，参数h是一个平滑调节因子。

其中，统计迭代公式为μ^iter=μ^init-α(A^T(Σ^-1(Aμ^init-y_raw)))-βR′(μ^init)，其中，A为CT设备的成像系统矩阵或系统投影矩阵，A^T表示A的转置；Σ为一对角矩阵，对角元素为对应投影数据y_raw的方差

；参数β为一非负实数，用于刻画先验约束的强度；参数α计算公式为

G=A^T(Σ^-1(Aμ^iter-y_raw))；R′(μ^init)计算公式为R′(μ^init)=psign(μ^init-SRNLM(μ^init))|μ^init-SRNLM(μ^init)|^p-1，参数p取值为区间[1,2]中的实数。

其中，所述的迭代终止条件为相邻两次对重建图像μ^iter进行迭代计算的结果的相对均方误差小于设定的阈值δ，即δ为一正实数。

其中，阈值δ=0.001。

与现有技术相比，本发明具有如下有益效果：

本发明通过引入参考图像引导的非局部均值先验信息，进行约束低剂量CT统计迭代重建，有效解决一般的边缘保持先验，如非局部均值先验，对图像中条形伪影的结构保持作用，可以在去除图像噪声的同时有效地抑制图像中的条形伪影，最终实现低剂量CT图像的优质重建。

附图说明

图1是本发明一个优选实施例的流程示意图。

图2是本发明实施例中仿真所使用的XCAT体模图像示意图。

图3是对图2使用本发明获得的CT重建图像。

图4是对图2直接通过FBP法获得的CT重建图像。

图5是对图2使用现有技术中基于非局部均值边缘保持先验的统计迭代重建方法获得的CT重建图像。

具体实施方式

本发明提供一种基于改进的边缘保持先验的低剂量X射线CT图像统计迭代重建方法，能够针对降低管电流和扫描时间的低剂量CT图像进行重建，能够有效地去除图像噪声并抑制条形伪影，同时较好地保持图像细节信息。

如图1所示，本发明的一个优选实施例包括如下实现步骤：

步骤S1.获取CT设备的系统参数和低剂量X射线CT图像的投影数据y_raw。

其中，CT设备的系统参数主要包括X射线入射光子强度I₀，系统电子噪声的方差

等。

步骤S2.对步骤S1中获取的投影数据y_raw进行数据恢复处理，得到恢复后的投影数据y_restored。

对投影数据y_raw进行数据恢复处理的数据恢复模型为

$y_{\mathrm{restored}}=\underset{y\≥0}{\arg \min }\{{(y_{\mathrm{raw}}-y)}^T{\mathrm{\Σ}}^{-1}(y_{\mathrm{raw}}-y)+\mathrm{\λR}\left(y\right)\}---\left(I\right)$

其中：y表示待求解的恢复后的投影数据；Σ为一对角矩阵，其对角线元素为对应投影数据y_raw的方差

；R(y)为先验约束项，可以选择二次平滑先验约束形式；λ>0为正则化调节参数，用于刻画先验约束的强度。

具体的，投影数据y_raw的方差

通过投影数据y_raw的逐个数据点估计得到，可以通过公式

估计得到，其中I₀为第i个投影数据y_raw的数据点的X射线入射光子强度，

为CT设备的系统电子噪声的方差。当然，方差

也可以通过局部邻域方差估计等其他方式得到。

具体的，式（I）中数据恢复模型的求解，可以通过梯度下降法或者高斯-塞德尔优化算法（Gauss-Seidel method）进行迭代求解，得到恢复后的投影数据y_restored。

步骤S3.对步骤S2中得到恢复后的投影数据y_restored采用解析重建方法进行图像重建，得到重建后的低剂量CT图像，作为统计迭代重建中先验构建的参考图像μ^ref。

其中，步骤S3中采用的解析重建方法可以为现有技术中的滤波反投影方法（Filtered Back-Projection，FBP）。

步骤S4.对步骤S1中获取的投影数据y_raw采用解析重建方法进行图像重建，得到重建后的低剂量X射线CT图像，作为统计迭代重建步中的迭代初始图像μ^init。

其中，步骤S4中采用的解析重建方法也可以为现有技术中的滤波反投影方法（Filtered Back-Projection，FBP）。

步骤S5.根据步骤S3得到的参考图像μ^ref和步骤S4得到的迭代初始图像μ^init，构建本发明改进的边缘保持先验。

上述改进的边缘保持先验，其先验构建形式为：

R(μ^init)=φ(μ^init-SRNLM(μ^init))  （II）其中，φ()为能量势函数，SRNLM(μ^init)为由参考图像μ^ref引导的非局部均值滤波。

具体的，能量势函数φ()可以为p范数形式，如：当p=1时，即为1范数形式，当p=2时，则为2范数形式。

具体的，由参考图像μ^ref引导的非局部均值滤波SRNLM(μ^init)，其表达式为：

$\mathrm{SRNLM}\left({\mathrm{\μ}}_i^{\mathrm{init}}\right)=\frac{1}{Z\left(i\right)}{\mathrm{\Σ}}_{j\∈N_i}\exp \{-{\left|\right|{\mathrm{\μ}}^{\mathrm{init}}\left(V_i\right)-{\mathrm{\μ}}^{\mathrm{ref}}\left(V_j\right)\left|\right|}_{2,a}^2/h^2\}{\mathrm{\μ}}_j^{\mathrm{init}}---\left(\mathrm{III}\right)$

其中，

为归一化因子，V_i和V_j分别为图像μ^init中像素点i处和图像μ^ref中像素点j处的两个局部相似窗，μ^init(V_i)和μ^ref(V_j)分别为图像μ^init在局部相似窗V_i内的图像矩阵块和图像μ^ref在局部相似窗V_j内的图像矩阵块。N_i为以像素点i为中心的搜索窗。符号||||₂表示计算两个相似窗内的图像矩阵块之间的欧氏距离；参数h是一个平滑调节因子，可以根据图像噪声水平来决定。

步骤S6.根据步骤S4中得到的迭代初始图像μ^init和步骤S5中构建的改进的边缘保持先验R(μ^init)，应用统计迭代公式进行迭代计算，得到迭代重建图像μ^iter。

统计迭代公式为μ^iter=μ^init-α(A^T(Σ^-1(Aμ^init-y_raw)))-βR′(μ^init)，其中:

A为CT设备的成像系统矩阵，也称为系统投影矩阵，由具体的CT成像系统决定，可以采用不同的方法来计算，比如：基于面积加权的系统矩阵计算，基于体素加权的系统矩阵计算；A^T表示A的转置。

Σ为一对角矩阵，对角元素为对应投影数据y_raw的方差

。

参数β为一非负实数，用于刻画先验约束的强度。

参数α计算公式为其中，G=A^T(Σ^-1(Aμ^iter-y_raw))。

R′(μ^init)计算公式为R′(μ^init)=psign(μ^init-SRNLM(μ^init))|μ^init-SRNLM(μ^init)|^p-1,其中参数p取值即步骤S5中所述的区间[1,2]中的一个实数。

步骤S7.判断迭代重建图像μ^iter是否满足迭代终止条件。若是，转入步骤S8若否，以迭代重建图像μ^iter更新迭代初始图像μ^init，重复步骤S6-S7。

所述的迭代终止条件为相邻两次迭代重建结果的相对均方误差小于设定的阈值δ，即

δ为一正实数。

比如，在一个实施例中，迭代终止条件中设定的阈值δ取值为0.001。

步骤S8.迭代终止，该迭代重建图像即是最终重建图像。

为了对比本发明所示方法的效果，采用图2所示的XCAT数字体模图像作为本发明的计算机仿真实验对象，该XCAT数字体模图像像素矩阵大小设为512×512，模拟CT机的X射线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间，采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.407mm。

通过CT系统仿真生成大小为1160×672的低剂量CT投影数据y_raw，其中X射线的入射光子强度I₀为1.0×10⁵，系统电子噪声的方差

为11.0。对模拟生成的同组CT图像的投影数据，采用本发明前述基于非局部均值边缘保持先验的统计迭代重建方法进行重建，得到CT重建图像，如图3所示。

图4给出了对图2直接采用FBP重建方法得到的CT重建图像。比较图4和图3中的重建结果，可以看出，本发明的统计迭代重建方法可以有效地减少重建图像中的噪声和伪影，且明显改善了CT重建图像的质量。

图5给出了对图2采用现有技术中基于非局部均值边缘保持先验的统计迭代重建方法得到的CT重建图像。将图3和图5进行比较可见，本发明方法能够在去除图像噪声的同时，更好地减少图像中的伪影。

综上，本发明通过引入参考图像引导的非局部均值先验信息，进行约束低剂量CT统计迭代重建，有效解决一般的边缘保持先验，如非局部均值先验，对图像中条形伪影的结构保持作用，可以在去除图像噪声的同时有效地抑制图像中的条形伪影，最终实现低剂量CT图像的优质重建。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (9)

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

获取CT设备的系统参数和低剂量X射线CT图像的投影数据y_raw，对投影数据y_raw进行图像重建，得到迭代初始图像μ^init；

对投影数据y_raw进行数据恢复处理，得到恢复后的投影数据y_restored，对恢复后的投影数据y_restored进行图像重建，得到参考图像μ^ref；

以参考图像μ^ref和迭代初始图像μ^init，按R(μ^init)=φ(μ^init-SRNLM(μ^init))构建边缘保持先验R(μ^init)，其中，φ()为能量势函数，SRNLM(μ^init)为由参考图像μ^ref引导的非局部均值滤波；

根据迭代初始图像μ^init边缘保持先验R(μ^init)，利用统计迭代公式进行迭代计算，得到迭代重建图像μ^iter；

当对重建图像μ^iter进行迭代计算的结果满足迭代终止条件时，迭代终止，获得低剂量X射线CT图像最终的重建图像。

2.根据权利要求1所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，对投影数据y_raw进行数据恢复处理的步骤，是通过迭代公式

进行迭代求解，得到恢复后的投影数据y_restored；

其中，y表示待恢复的投影数据；Σ为一对角矩阵，其对角线元素为对应投影数据y_raw的方差

；R(y)为先验约束项；λ>0为正则化调节参数。

3.根据权利要求2所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，投影数据y_raw的方差

通过投影数据y_raw的逐个数据点估计得到，估计公式为

其中I₀为投影数据y_raw的第i个数据点的X射线入射光子强度，为CT设备的系统电子噪声的方差。

4.根据权利要求2所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，所述迭代公式使用梯度下降法或者高斯-塞德尔优化算法进行迭代求解获得恢复后的投影数据y_restored。

5.根据权利要求1所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，对投影数据y_raw及对恢复后的投影数据y_restored均是采用滤波反投影方法进行图像重建。

6.根据权利要求1所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，由参考图像μ^ref引导的非局部均值滤波SRNLM(μ^init)的计算公式为：

$\mathrm{SRNLM}\left({\mathrm{\μ}}_i^{\mathrm{init}}\right)=\frac{1}{Z\left(i\right)}{\mathrm{\Σ}}_{j\∈N_i}\exp \{-{\left|\right|{\mathrm{\μ}}^{\mathrm{init}}\left(V_i\right)-{\mathrm{\μ}}^{\mathrm{ref}}\left(V_j\right)\left|\right|}_{2,a}^2/h^2\}{\mathrm{\μ}}_j^{\mathrm{init}}$

其中，

为归一化因子，V_i和V_j分别为μ^init图像中像素点i和j处的两个局部相似窗，μ^init(V_i)和μ^ref(V_j)分别为迭代初始图像μi^nit在局部相似窗V_i内的图像矩阵块和参考图像μ^ref在局部相似窗V_j内的图像矩阵块，N_i为以像素点i为中心的搜索窗，符号||||₂表示计算两个相似窗内的图像矩阵块之间的欧氏距离，参数h是一个平滑调节因子。

7.根据权利要求1所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，统计迭代公式为μ^iter=μ^init-α(A^T(Σ^-1(Aμ^init-y_raw)))-βR′(μ^init)，其中，A为CT设备的成像系统矩阵或系统投影矩阵，A^T表示A的转置；Σ为一对角矩阵，对角元素为对应投影数据y_raw的方差

参数β为一非负实数，用于刻画先验约束的强度；参数α计算公式为

G=A^T(Σ^-1(Aμ^iter-y_raw))；R′(μ^init)计算公式为R′(μ^init)=psign(μ^init-SRNLM(μ^init))|μi^nit-SRNLM(μ^init)|^p-1，参数p取值为区间[1,2]中的实数。

8.根据权利要求1-7任何一项所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，所述的迭代终止条件为相邻两次对重建图像μ^iter进行迭代计算的结果的相对均方误差小于设定的阈值δ，即

δ为一正实数。

9.根据权利要求8所述一种低剂量X射线CT图像统计迭代重建方法，其特征在于，阈值δ=0.001。

Images