CN101980302A - 投影数据恢复导引的非局部平均低剂量ct重建方法 - Google Patents

Abstract

一种投影数据恢复导引的非局部平均低剂量CT重建方法，包括以下步骤：(1)利用CT成像设备采集低剂量投影数据；(2)对低剂量投影数据进行变换；(3)对步骤(1)获取的低剂量投影数据进行传统的FBP重建；(4)对变换后的低剂量投影数据进行滤波；(5)对步骤(4)滤波后的低剂量投影数据进行传统的FBP重建；(6)对步骤(5)获取的标准剂量图像计算权值矩阵；(7)利用步骤(6)获取的权值矩阵对步骤(3)获取的低剂量重建图像进行加权平均滤波，得到恢复后的低剂量图像。该方法在使用低剂量放射的同时，仍能保证产生高质量的CT重建图像，本发明方法具有良好的鲁棒性，在噪声消除和伪影抑制两方面均有上佳表现。

Description

投影数据恢复导引的非局部平均低剂量CT重建方法

技术领域

本发明涉及一种医学影像的图像重建方法，具体来说涉及一种投影数据恢复导引的非局部平均低剂量CT重建方法。

背景技术

近些年来，病人接受CT检查时所接受的放射剂量及其相关的癌症风险，正越来越受到关注。随着多排CT，双源CT在临床中的广泛应用，新的CT设备的使用造成更大的X射线使用剂量，使得人们越来越关注于CT剂量对人体的潜在危害以及如何在保证图像质量的前提下最大限度地降低X射线使用剂量。

传统的射线剂量降低方法主要包括降低X射线球管电压，球管电流，以及缩短暴光时间等。但是，采用上述方法重建的图像质量将严重退化，难以满足临床诊断需求。当前，除了先进的硬件方法外，解决低剂量CT图像重建的另一有效途径是优良的重建方法。例如，利用投影数据噪声的统计特性对实际的低剂量成像系统进行建模，采用迭代算法对测量数据进行重建。然而，由于CT投影数据量庞大，迭代计算耗时长，难以满足临床中实时交互的要求，其方法多限于理论探讨。另外，直接对低剂量图像进行滤波的后处理方法也只能从图像自身获取信息，其处理精度强依赖于图像数据的噪声和伪影特性。然而低剂量图像的噪声和伪影特性极不稳定且难以精确刻画，也使得多数基于后处理的低剂量CT图像滤波方法很难获得临床应用。另外，基于投影数据滤波的低剂量重建方法可以有效抑制重建图像中的噪声和伪影，然而其重建图像的分辨率相对于未经处理的投影数据直接重建图像有所下降，这也是该类方法实际应用中所面对的最大问题。

发明内容

本发明的目的在于提供一种投影数据恢复导引的非局部平均低剂量CT重建方法，该方法在使用低剂量放射的同时，仍能保证产生高质量的CT重建图像，本发明方法具有良好的鲁棒性，在噪声消除和伪影抑制两方面均有上佳表现。

本发明的目的可通过以下的技术措施来实现：

一种投影数据恢复导引的非局部平均低剂量CT重建方法，其特征在于包括以下步骤：

(1)利用CT成像设备采集低剂量投影数据；

(2)对步骤(1)获取的低剂量投影数据进行变换；

(3)对步骤(1)获取的低剂量投影数据进行传统的FBP重建；

(4)对步骤(2)变换后的低剂量投影数据进行滤波；

(5)对步骤(4)滤波后的低剂量投影数据进行传统的FBP重建；

(6)对步骤(5)获取的标准剂量图像计算权值矩阵；

(7)利用步骤(6)获取的权值矩阵对步骤(3)获取的低剂量重建图像进行加权平均滤波，得到恢复后的低剂量图像。

所述步骤(2)中的数据变换为Anscombe变换。

在进行所述步骤(4)前，将所述低剂量投影数据中满足Poisson分布的低剂量投影数据转化成Gaussian分布后，再进行步骤(4)。

所述步骤(4)中的滤波采用BM3D滤波。

在进行所述步骤(5)前，将所述滤波后的低剂量投影数据进行Anscombe逆变换后，再进行步骤(5)。

所述步骤(6)中的计算权值矩阵的过程为：

(a)选择一个包含图像中几何信息的大方形邻域，在大方形邻域内选择两个大小相同的小方形邻域；通过相似性测度计算两个小方形邻域的欧氏距离。

(b)在选定的小方形邻域内进行两像素间的灰度值比较的同时，利用两像素间相似性来获得刻画两像素关系中的权值量。

所述步骤(a)中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数。

所述步骤(b)中的权值量定义为

g_kj＝exp(-D(k，j)/h²)，

λ(v_k)＝{λ(l)：l∈v_k}，λ(v_j)＝{λ(l)：l∈v_j}，N_j表示像素点j处能够包含图像几何结构信息的大方形邻域，所述大方形领域以像素点j为中心；ω_kj为三j维权值矩阵中一个值，其中x，y轴坐标为像素点k在原图像中的坐标，z轴为以像素点k为中心的方形邻域中像素点j的位置，通过图像域内两像素点间的欧几里德距离的反比例函数确定；V_k和V_j为以像素点k和像素点j为中心的小方形邻域，所述小方形领域分别以像素点k和像素点j为中心；λ(V_k)和λ(V_j)为此两个小方形邻域中所有像素灰度值数组||·||代表此两个像素点所在区域的加权欧几里德距离；参数h为计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。

所述步骤(7)中的加权平均滤波定义为

其中N_k表示像素点k处的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。

本发明的CT重建方法相比现有技术具有以下有益效果：

1、本方法有效地结合CT投影域数据恢复及图像域数据恢复的优势，结合NL-means滤波权值设计的特点，将两者巧妙地联系融合，以实现低剂量CT图像的优质重建；

2、本文方法较好的保持的图像边缘，并且抑制匀质区域的噪声，同时具有较高的信噪比。

附图说明

图1是本发明投影数据恢复导引的非局部平均低剂量CT重建方法的流程图；

图2(a)是修定的Shepp-Logan体模图像数据；

图2(b)是图2(a)中方框所示的高衰减区域的局部放大示意图；

图3(a)是Shepp-Logan体模低剂量数据采用FBP算法直接重建图像；

图3(b)是Shepp-Logan体模低剂量数据采用双边滤波后的FBP重建图像；

图3(c)是Shepp-Logan体模低剂量数据采用BM3D滤波后的FBP重建图像；

图3(d)是Shepp-Logan体模采用本文方法非局部平均低剂量CT重建方法重建图像；

图4(a)至图4(d)分别是图3(a)至图3(d)的重建图像中方框的局部放大图像；

图5是图3(a)至图3(d)的不同方法所得的重建图像中线[150 250]区间剖面图；

图6是图3(a)至图3(d)的不同方法的剂量条件下不同重建方法获得的均方误差统计箱线图；

图7(a)是仿真的临床10mAs投影数据采用直接FBP重建的腹部图像；

图7(b)是仿真的临床10mAs投影数据采用投影域双边滤波后FBP重建的腹部图像；

图7(c)是仿真的临床10mAs投影数据采用投影域BM3D滤波后FBP重建的腹部图像；

图7(d)是仿真的临床10mAs投影数据采用本发明方法获得的重建的腹部图像；

图8(a)至图8(d)分别是图7(a)至图7(d)的重建图像中方框的局部放大图像；

图9(a)是临床10mAs投影数据采用直接FBP重建的肝部图像；

图9(b)是临床10mAs投影数据采用投影域双边滤波后FBP重建的肝部图像；

图9(c)是临床10mAs投影数据采用投影域BM3D滤波后FBP重建的肝部图像；

图9(d)是临床10mAs投影数据采本发明方法获得的重建的肝部图像；

图10(a)至图10(d)分别是图9(a)至图9(d)的重建图像中方框的局部放大图像。

具体实施方式

本发明的投影数据恢复导引的非局部平均低剂量CT重建方法的具体实施步骤如图1所示，具体如下：

1、利用CT设备采集低剂量投影数据，射线剂量为标准剂量的1/10至1/20。

2、对采集的低剂量投影数据进行变换，将满足Poisson泊松分布的低剂量投影数据转化成Gaussian高斯分布，具体为：假设服从Poisson分布的变量x均值为m，方差为v，经过Anscombe安斯库姆变换：

则得到服从Gaussian分布且方差近似为1的变量。

3、对获取的低剂量投影数据直接进行传统的FBP(Filter back-projection，FBP，滤波反投影)重建。

4、变换后的低剂量投影数据进行BM3D滤波(Block-Matching and 3D filtering，BM3D)。

5、对滤波后的低剂量投影数据进行传统的FBP重建。

6、对运用滤波后低剂量投影数据FBP重建的图像计算权值矩阵，权值矩阵为一三维矩阵，即对滤波后重建图像的每一像素点k，有一p＝N₀维列向量，其中N₀为以k为中心一较大方形邻域内的像素个数。列向量中的每一点ω_kj为像素点j与像素点k的相似性测度，j∈N₀。

ω_kj的具体计算方法为：

g_kj＝exp(-D(k，j)/h²)，

λ(v_k)＝{λ(l)：l∈v_k}，λ(v_j)＝{λ(l)：l∈v_j}，N_j表示像素点j处能够包含图像几何结构信息的大放邻域(称之为搜索邻域)，几何信息包括边缘、拐角和纹理等；V_k和V_j则设定为以像素点k和像素点j为中心的小方形邻域(相似邻域)；λ(V_k)和λ(V_j)为此两个邻域中所有像素灰度值数组；||·||代表此两个像素点所在区域的加权欧几里德距离；参数h用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系，即ω_kj通过图像域内两像点间的欧几里德距离的反比例函数确定。本方法实验中，综合处理效果和计算复杂度的考虑，相似邻域窗和搜索邻域窗分别设为7×7和11×11。参数h的选取与图像噪声的方差相关，可记为h＝F(σ)，σ为图像噪声方差。本方法实验中，h为人工设定。

7、利用获取的权值矩阵对未经处理的低剂量投影数据直接FBP重建图像进行加权平均滤波，得到恢复后的低剂量图像。其中的加权平均滤波定义为N_k表示像素点k处的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。

下面对具体数据采用本方法重建图像来说明本发明方法的效果。

首先采用如图2(a)所示的修定的shepp-logan体模图像以及该图像的局部细节图像图2(b)作为本发明的模拟实验对象。体模大小设为512×512，射源到旋转中心和探测器的距离分别为541mm和949.975mm，旋转角在[0，2π]间采样值为984，探测器单元为1.0239mm。通过转换概率矩阵A得到投影数据(sinogram)，然后通过调整光子总数值模拟生成低剂量投影数据，成像几何由Fessler等人提供的ASPIRE软件系统生成。重建过程中使用传统扇形束FBP算法，滤波过程采用Hanning窗进行，其截止频率设为奈奎斯特频率的80％。

图3(a)至图3(d)描述了在相比标准剂量图像，未经任何处理的直接FBP低剂量图像沿体模长轴方向会产生严重的条形伪影，本文方法和其他方法的效果对比。图3(a)为低剂量投影数据直接进行传统FBP的图像，图3(b)为对变换后的低剂量投影数据进行双边滤波后将传统FBP图像，其中体模长轴方向噪声有所减弱，但依然明显。图3(c)为利用BM3D对变换后的低剂量投影数据进行滤波后FBP重建获得的图像；图3(d)为本发明方法恢复的低剂量图像。可以对比看出本发明方法保持了良好的边缘并且在计量衰减大的区域有更佳的抑制噪声的效果。图4(a)至图4(d)分别是图3(a)至图3(d)的重建图像中方框的局部放大图像。

从图5描述了对低剂量投影数据重建图像上述图2(a)至图2(d)4种方法结果的水平剖面图，可以看出，本文方法较好的保持的图像边缘，并且抑制匀质区域的噪声。

对图3(a)、3(b)、3(c)、3(d)所示的采用不同方法重建后图像中方框内的匀质区域计算相应信噪比SNR(Sino-to-noise ratio)，如下表1所示。另外，不同方法获得的重建图像的归一化均方误差值NMSE(Normalized Mean Square Error)，如下表2所示。

表2

上面两个表可以定量说明本发明方法在抑制噪声和保持边缘两方面的上佳表现。

为了评估本文方法在不同(低)剂量条件下的重建图像质量，分析10组不同的低剂量投影数据(入射光子数在[1e4 5.5e4]区间)，其SNR如下表3所示。

表3

实验中分别采用直接的FBP重建，基于投影数据双边滤波的FBP重建，基于投影数据BM3D滤波的重建和本文方法对每一组数据实现低剂量CT重建，其重建结果分别与理想图像进行均方误差测量。为保持实验的公平性，10组数据的重建过程中不同方法设所计的参数均保持不变。运用Matlab中的统计工具包对获得均方误差进行统计分析以生成箱线图，如图6所示。

图7(a)至图7(d)中实验通过采用一组临床腹部低剂量投影数据来进一步验证本文方法的有效性。所有滤波图像均在同一窗宽窗位下显示，显示窗大小为[0.000765，0.001120]。数据由商业单排CT扫描仪采集，射源到旋转中心和探测器的距离分别为475.00mm和845.24mm，旋转角在间采样值为480，探测器单元为0.92mm，共384个。球管电压和球管电流分别为120kv和120mAs，光子计数为2.5×10⁶。通过降低光子计数为2.5×10⁵生成低剂量的投影数据(约10mAs)。图7(a)为低剂量投影数据直接进行传统FBP重建，图7(b)为采用双边滤波后的FBP算法重建图像，图7(c)为采用BM3D滤波后的FBP算法重建图像，图7(d)本发明方法恢复的低剂量图像。可以看出，由于X射线剂量较低，使用FBP直接重建的图像存在较多的噪声和伪影，尤其在高衰减区更为严重。相比而言，本文方法的重建结果较好的抑制了图像中的噪声。

图8(a)至图8(d)显示了图7(a)至图7(d)中方框内的图像，低剂量导致高衰减区出现严重的放射状伪影，辅助说明本发明方法的效果。对图7(a)至图7(d)选定的匀质区域计算信噪比，如下表4所示，可验证本发明方法的有效性：

表4

图9(a)至图9(d)为临床肝部低剂量投影数据不同重建方法重建的实验结果对比图。所有滤波图像均在同一窗宽窗位下显示，显示窗大小为[0.008930.01507]，搜索窗为11×11，相似窗为7×7，h为0.0002。图10(a)至图10(d)为图9(a)至图9(d)方框内区域的局部放大图，辅助说明重建图像的效果。

本发明的实施方式不限于此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (10)

1.一种投影数据恢复导引的非局部平均低剂量CT重建方法，其特征在于包括以下步骤：

(1)利用CT成像设备采集低剂量投影数据；

(2)对步骤(1)获取的低剂量投影数据进行变换；

(3)对步骤(1)获取的低剂量投影数据进行传统的FBP重建；

(4)对步骤(2)变换后的低剂量投影数据进行滤波；

(5)对步骤(4)滤波后的低剂量投影数据进行传统的FBP重建；

(6)对步骤(5)获取的标准剂量图像计算权值矩阵；

(7)利用步骤(6)获取的权值矩阵对步骤(3)获取的低剂量重建图像进行加权平均滤波，得到恢复后的低剂量图像。

2.根据权利要求1所述的CT重建方法，其特征在于：所述步骤(2)中的数据变换为安斯库姆变换。

3.根据权利要求1所述的CT重建方法，其特征在于：在进行所述步骤(4)前，将所述低剂量投影数据中满足Poisson分布的低剂量投影数据转化成Gaussian分布后，再进行步骤(4)。

4.根据权利要求1所述的CT重建方法，其特征在于：所述步骤(4)中的滤波采用BM3D滤波。

5.根据权利要求1所述的CT重建方法，其特征在于：在进行所述步骤(5)前，将所述滤波后的低剂量投影数据进行Anscombe逆变换后，再进行步骤(5)。

6.根据权利要求1所述的CT重建方法，其特征在于：所述步骤(6)中的计算权值矩阵的过程为：

(a)选择一个包含图像中几何信息的大方形邻域，在大方形邻域内选择两个大小相同的小方形邻域，通过相似性测度计算两个小方形邻域的欧氏距离；

(b)在选定的小方形邻域内进行两像素间的灰度值比较的同时，利用两像素间相似性来获得刻画两像素关系中的权值量。

7.根据权利要求6所述的CT重建方法，其特征在于：所述步骤(a)中的几何信息包括边缘、拐角和纹理。

8.根据权利要求6所述的CT重建方法，其特征在于：所述步骤(a)中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数。

9.根据权利要求6所述的CT重建方法，其特征在于：所述步骤(b)中的权值量定义为

g_kj＝exp(-D(k，j)/h²)，

λ(v_k)＝{λ(l)：l∈v_k}，λ(v_j)＝{λ(l)：l∈v_j}，N_j表示像素点j处能够包含图像几何结构信息的大方形邻域，所述大方形领域以像素点j为中心；ω_kj为三j维权值矩阵中一个值，其中x，y轴坐标为像素点k在原图像中的坐标，z轴为以像素点k为中心的方形邻域中像素点j的位置，通过图像域内两像素点间的欧几里德距离的反比例函数确定；V_k和V_j为以像素点k和像素点j为中心的小方形邻域，所述小方形领域分别以像素点k和像素点j为中心；λ(V_k)和λ(V_j)为此两个小方形邻域中所有像素灰度值数组；||·||代表此两个像素点所在区域的加权欧几里德距离；参数h为计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。

10.根据权利要求1所述的CT重建方法，其特征在于：所述步骤(7)中的加权平均滤波定义为

其中N_k表示像素点k处的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。

Images