专利名称：投影数据恢复导引的非局部平均低剂量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)获取的低剂量重建图像进行加权平 均滤波，得到恢复后的低剂量图像。
4
所述步骤(2)中的数据变换为Anscombe变换。在进行所述步骤(4)前，将所述低剂量投影数据中满足Poisson分布的低剂量投 影数据转化成Gaussian分布后，再进行步骤(4)。所述步骤(4)中的滤波采用BM3D滤波。在进行所述步骤(5)前，将所述滤波后的低剂量投影数据进行Anscombe逆变换 后，再进行步骤(5)。所述步骤(6)中的计算权值矩阵的过程为(a)选择一个包含图像中几何信息的大方形邻域，在大方形邻域内选择两个大小 相同的小方形邻域；通过相似性测度计算两个小方形邻域的欧氏距离。(b)在选定的小方形邻域内进行两像素间的灰度值比较的同时，利用两像素间相 似性来获得刻画两像素关系中的权值量。所述步骤(a)中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧 几里德距离的反比例函数。所述步骤(b)中的权值量定义为 =叫^^ ’ gk. = θχρ (-D (k，j) /h2)，
义(V7)厄，λ (Vk) = {λ (1) :1 e vk}，λ (Vj) = {λ (1) :1 e 乃}， 表示像 素点j处能够包含图像几何结构信息的大方形邻域，所述大方形领域以像素点j为中心； "kJ为三j维权值矩阵中一个值，其中χ，y轴坐标为像素点k在原图像中的坐标，ζ轴为以 像素点k为中心的方形邻域中像素点j的位置，通过图像域内两像素点间的欧几里德距离 的反比例函数确定；Vk和、为以像素点k和像素点j为中心的小方形邻域，所述小方形领 域分别以像素点k和像素点j为中心；λ (Vk)和λ (Vj)为此两个小方形邻域中所有像素灰 度值数组I I ·| I代表此两个像素点所在区域的加权欧几里德距离；参数h为计算像素点间 权值的指数函数同邻域相似性测度的反比例衰减关系。所述步骤(7)中的加权平均滤波定义为= Σ’其中Nk表示像素点k处
的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。本发明的CT重建方法相比现有技术具有以下有益效果1、本方法有效地结合CT投影域数据恢复及图像域数据恢复的优势，结合 NL-means滤波权值设计的特点，将两者巧妙地联系融合，以实现低剂量CT图像的优质重 建；2、本文方法较好的保持的图像边缘，并且抑制勻质区域的噪声，同时具有较高的 信噪比。
图1是本发明投影数据恢复导引的非局部平均低剂量CT重建方法的流程图；图2 (a)是修定的Sh印p-Logan体模图像数据；图2(b)是图2(a)中方框所示的高衰减区域的局部放大示意图；图3 (a)是Sh印p-Logan体模低剂量数据采用FBP算法直接重建图像；图3 (b)是Sh印p-Logan体模低剂量数据采用双边滤波后的FBP重建图像；
图3 (c)是Sh印p-Logan体模低剂量数据采用BM3D滤波后的FBP重建图像；图3 (d)是Sh印p-Logan体模采用本文方法非局部平均低剂量CT重建方法重建图 像；图4(a)至图4(d)分别是图3(a)至图3(d)的重建图像中方框的局部放大图像；图5是图3(a)至图3(d)的不同方法所得的重建图像中线[150 250]区间剖面 图；图6是图3(a)至图3(d)的不同方法的剂量条件下不同重建方法获得的均方误差 统计箱线图；图7 (a)是仿真的临床IOmAs投影数据采用直接FBP重建的腹部图像；图7 (b)是仿真的临床IOmAs投影数据采用投影域双边滤波后FBP重建的腹部图 像；图7 (c)是仿真的临床IOmAs投影数据采用投影域BM3D滤波后FBP重建的腹部图 像；图7(d)是仿真的临床IOmAs投影数据采用本发明方法获得的重建的腹部图像；图8(a)至图8(d)分别是图7(a)至图7(d)的重建图像中方框的局部放大图像；图9 (a)是临床IOmAs投影数据采用直接FBP重建的肝部图像；图9 (b)是临床IOmAs投影数据采用投影域双边滤波后FBP重建的肝部图像；图9 (c)是临床IOmAs投影数据采用投影域BM3D滤波后FBP重建的肝部图像；图9 (d)是临床IOmAs投影数据采本发明方法获得的重建的肝部图像；图10(a)至图10(d)分别是图9 (a)至图9(d)的重建图像中方框的局部放大图像。
具体实施例方式本发明的投影数据恢复导引的非局部平均低剂量CT重建方法的具体实施步骤如 图1所示，具体如下1、利用CT设备采集低剂量投影数据，射线剂量为标准剂量的1/10至1/20。2、对采集的低剂量投影数据进行变换，将满足Poisson泊松分布的低剂量投影数 据转化成Gaussian高斯分布，具体为假设服从Poisson分布的变量χ均值为m，方差为V，
经过Anscombe安斯库姆变换d χ则得到服从Gaussian分布且方差近似为1的变量。3、对获取的低剂量投影数据直接进行传统的FBP(Filter back-projection，FBP， 滤波反投影)重建。4、变换后的低剂量投影数据进行BM3D滤波(Block-Matching and 3D filtering, BM3D)。5、对滤波后的低剂量投影数据进行传统的FBP重建。6、对运用滤波后低剂量投影数据FBP重建的图像计算权值矩阵，权值矩阵为一三 维矩阵，即对滤波后重建图像的每一像素点k，有一 ρ = N0维列向量，其中Ntl为以k为中心 一较大方形邻域内的像素个数。列向量中的每一点《kj为像素点j与像素点k的相似性测
度，j e N00
Qkj的具体计算方法为 轻仏知‘gk. = eXp(-D(k, j)/h2)，
D(k,j) =IIA(Vi)-^(V7)H22, λ (Vk) = {λ ⑴l G vk}, λ (Vj) = {λ ⑴1 G Vj}，Nj 表示
像素点j处能够包含图像几何结构信息的大放邻域(称之为搜索邻域)，几何信息包括边 缘、拐角和纹理等；Vk和、则设定为以像素点k和像素点j为中心的小方形邻域(相似邻 域)；λ (Vk)和λ (Vj)为此两个邻域中所有像素灰度值数组；ι ι · 11代表此两个像素点所 在区域的加权欧几里德距离；参数h用于计算像素点间权值的指数函数同邻域相似性测度 的反比例衰减关系，即通过图像域内两像点间的欧几里德距离的反比例函数确定。本 方法实验中，综合处理效果和计算复杂度的考虑，相似邻域窗和搜索邻域窗分别设为7X7 和11X11。参数h的选取与图像噪声的方差相关，可记为h = F(o)，σ为图像噪声方差。 本方法实验中，h为人工设定。7、利用获取的权值矩阵对未经处理的低剂量投影数据直接FBP重建图像进行加 权平均滤波，得到恢复后的低剂量图像。其中的加权平均滤波定义为ξ Α/ΦΟ ’Nk
表示像素点k处的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。下面对具体数据采用本方法重建图像来说明本发明方法的效果。首先采用如图2(a)所示的修定的sh印p-logan体模图像以及该图像的局部细节 图像图2(b)作为本发明的模拟实验对象。体模大小设为512X512，射源到旋转中心和探 测器的距离分别为541mm和949. 975mm，旋转角在W，2ji]间采样值为984，探测器单元为 1.0239mm。通过转换概率矩阵A得到投影数据(sinogram)，然后通过调整光子总数值模拟 生成低剂量投影数据，成像几何由Fessler等人提供的ASPIRE软件系统生成。重建过程中 使用传统扇形束FBP算法，滤波过程采用Harming窗进行，其截止频率设为奈奎斯特频率的 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(Sin0-t0-n0ise ratio)，如下表1所示。另外，不同方法获得的重建 图像的归一化均方误差值NMSE (Normalized Mean Square Error)，如下表2所示。
权利要求
一种投影数据恢复导引的非局部平均低剂量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.根据权利要求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)中的权值量定义为’ =^|>,、= exp(-D(k, j)/h2)，D{Kj)=U{yk)-^Vj)f2, λ (vk) = {λ (1)1 e vj , λ (Vj) = {λ (1) 1 e Vj}，Nj表示像素点j处能够包含图像几何结构信息的大方 形邻域，所述大方形领域以像素点j为中心；《kj为三j维权值矩阵中一个值，其中χ，y轴 坐标为像素点k在原图像中的坐标，ζ轴为以像素点k为中心的方形邻域中像素点j的位 置，通过图像域内两像素点间的欧几里德距离的反比例函数确定；Vk和Vj为以像素点k和 像素点j为中心的小方形邻域，所述小方形领域分别以像素点k和像素点j为中心；λ (Vk) 和λ (Vj)为此两个小方形邻域中所有像素灰度值数组；N · I I代表此两个像素点所在区域 的加权欧几里德距离；参数h为计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。
10.根据权利要求1所述的CT重建方法，其特征在于所述步骤(7)中的加权平均滤 波定义为= Σ’其中Nk表示像素点k处的大方形邻域，u为待恢复低剂量图像，f(k)为恢复后图像在像素点k的值。
全文摘要
一种投影数据恢复导引的非局部平均低剂量CT重建方法，包括以下步骤(1)利用CT成像设备采集低剂量投影数据；(2)对低剂量投影数据进行变换；(3)对步骤(1)获取的低剂量投影数据进行传统的FBP重建；(4)对变换后的低剂量投影数据进行滤波；(5)对步骤(4)滤波后的低剂量投影数据进行传统的FBP重建；(6)对步骤(5)获取的标准剂量图像计算权值矩阵；(7)利用步骤(6)获取的权值矩阵对步骤(3)获取的低剂量重建图像进行加权平均滤波，得到恢复后的低剂量图像。该方法在使用低剂量放射的同时，仍能保证产生高质量的CT重建图像，本发明方法具有良好的鲁棒性，在噪声消除和伪影抑制两方面均有上佳表现。
文档编号A61B6/03GK101980302SQ20101051753
公开日2011年2月23日 申请日期2010年10月22日 优先权日2010年10月22日
发明者刘楠, 陈武凡, 马建华, 黄静 申请人:南方医科大学