计算机应用   2016, Vol. 36 Issue (10): 2916-2921  DOI: 10.11772/j.issn.1001-9081.2016.10.2916
0

引用本文 

何琳, 张权, 上官宏, 张文, 张鹏程, 刘祎, 桂志国. 自适应加权全变分的低剂量CT统计迭代算法[J]. 计算机应用, 2016, 36(10): 2916-2921.DOI: 10.11772/j.issn.1001-9081.2016.10.2916.
HE Lin, ZHANG Quan, SHANGGUAN Hong, ZHANG Wen, ZHANG Pengcheng, LIU Yi, GUI Zhiguo. Statistical iterative algorithm based on adaptive weighted total variation for low-dose CT[J]. JOURNAL OF COMPUTER APPLICATIONS, 2016, 36(10): 2916-2921. DOI: 10.11772/j.issn.1001-9081.2016.10.2916.

基金项目

国家自然科学基金资助项目(61271357);国家重大科学仪器设备开发专项(2014YQ24044508);山西省自然科学基金资助项目(2015011046);中北大学2013年校科学基金资助项目

通信作者

桂志国(1972—),男,天津蓟县人,教授,博士,主要研究方向:信号与信息处理、图像处理和识别、图像重建,E-mail:gzgtg@163.com

作者简介

何琳(1991—),女,山西运城人,硕士研究生,主要研究方向:图像处理与重建;
张权(1974—),男,山西大同人,副教授,博士,主要研究方向:图像处理、科学可视化;
上官宏(1988—),女,山西临汾人,博士研究生,主要研究方向:图像处理、医学图像重建;
张文(1992—),男,湖北孝感人,硕士研究生,主要研究方向:图像处理与重建;
张鹏程(1984—),男,内蒙古巴彦淖尔人,讲师,博士,主要研究方向:剂量计算、方案优化;
刘祎(1987—),女,河南睢县人,讲师,博士,主要研究方向:图像处理、医学图像重建

文章历史

收稿日期:2016-03-16
修回日期:2016-06-14
自适应加权全变分的低剂量CT统计迭代算法
何琳1, 张权1, 上官宏1, 张文1, 张鹏程1, 刘祎1, 桂志国1,2    
1. 电子测试技术国家重点实验室(中北大学), 太原 030051 ;
2. 仪器科学与动态测试教育部重点实验室(中北大学), 太原 030051
摘要: 针对低剂量计算机断层扫描(LDCT)重建图像时出现条形伪影和脉冲噪声的现象,提出一种自适应加权全变分的LDCT统计迭代重建算法。该算法克服了传统全变分(TV)算法在去除条形伪影的同时引入阶梯效应的缺点,把基于加权方差的加权因子与TV模型相结合提出自适应加权全变分模型,然后再把新模型应用到惩罚加权最小二乘(PWLS)重建算法中,这样就可以对图像的不同区域进行不同强度的去噪,从而取得噪声抑制和边缘保持的良好效果。采用Shepp-Logan模型和数字骨盆体模来验证算法的有效性,实验结果表明,所提算法的归一化均方距离和归一化平均绝对距离均比滤波反投影(FBP)、PWLS、惩罚加权最小二乘的中值先验(PWLS-MP)以及惩罚加权最小二乘的全变分(PWLS-TV)算法的值小,且可分别获得40.91 dB和42.25 dB的峰值信噪比。实验结果表明,该算法重建出的图像在有效去除条形伪影的同时对图像的边缘和细节起到很好的保护作用。
关键词: 低剂量计算机断层扫描    统计迭代重建    惩罚加权最小二乘    全变分    加权方差    
Statistical iterative algorithm based on adaptive weighted total variation for low-dose CT
HE Lin1, ZHANG Quan1, SHANGGUAN Hong1, ZHANG Wen1, ZHANG Pengcheng1, LIU Yi1, GUI Zhiguo1,2     
1. National Key Laboratory for Electronic Measurement Technology (North University of China), Taiyuan Shanxi 030051, China ;
2. Key Laboratory of Instrumentation Science and Dynamic Measurement, Ministry of Education (North University of China), Taiyuan Shanxi 030051, China
Abstract: Concerning the streak artifacts and impulse noise of the Low-Dose Computed Tomography (LDCT) reconstructed images, a statistical iterative reconstruction method based on adaptive weighted Total Variation (TV) for LDCT was presented. Considering the shortage that traditional TV may bring staircase effect while suppressing streak artifacts, an adaptive weighted TV model that combined the weighting factor based on weighted variation and TV model was proposed. Then, the new model was applied to the Penalized Weighted Least Square (PWLS). Different areas of the image were processed with different de-noising intensities, so as to achieve a good effect of noise suppression and edge preservation. The Shepp-Logan model and the digital pelvis phantom were used to test the effectiveness of the proposed algorithm. Experimental results show that the proposed method has smaller Normalized Mean Square Distance (NMSD) and Normal Average Absolute Distance (NAAD) in the two experiment images, compared with the Filtered Back Projection (FBP), PWLS, PWLS-Median Prior (PWLS-MP) and PWLS-TV algorithms. Meanwhile, the proposed method get Peak Signal-To-Noise Ratio (PSNR) of 40.91 dB and 42.25 dB respectively. Experimental results show that the proposed algorithm can well preserve image details and edges, while eliminating streak artifacts effectively.
Key words: Low-Dose Computed Tomography (LDCT)    statistical iterative reconstruction    Penalized Weighted Least Square (PWLS)    Total Variation (TV)    weighted variation    
0 引言

计算机断层扫描(Computed Tomography,CT)技术在临床医学应用方面取得令人满意的成绩,但高剂量的CT辐射会对人体健康带来极大的隐患,因此,近年来CT学者往往在降低CT辐射剂量的同时尽可能多地重建出低噪声、高分辨率的CT图像。低剂量CT(Low-Dose CT,LDCT)重建的方法主要有投影域降噪算法、后处理算法和统计迭代重建算法,本文主要采用统计迭代重建算法。

统计迭代重建算法也可称为图像域直接重建算法,其优化的变量为重建图像本身,同时由于其充分研究并利用了图像数据噪声的统计数学模型,因此可以得到高质量的重建图像。近年来,国内外的许多学者针对LDCT的统计迭代重建算法进行了相关的研究。Lui等[1]提出了一种新的噪声补偿 CT重建方法,根据局部噪声的统计特性对无噪图像进行无参估计,从而进一步处理LDCT图像中含有的非稳态噪声,最终获得高信噪比的重建图像。Rust等[2]使用非线性高斯滤波方法对重建图像进行滤波降噪,从而得到较好的降噪效果且能够很好地保留图像的边缘和细节。Chen等[3]提出一种基于局部相似性的混合先验模型,通过邻域像素间的高斯加权距离来定义先验信息,将该先验模型用于最大后验概率重建算法中可以获得优质的图像。Zhang等[4-5]把基于非局部均值和自适应非局部均值的先验模型用于LDCT图像的统计迭代重建算法中,相比之下,自适应非局部先验模型在去除条形伪影的同时能保留图像更多的边缘细节信息。Wang等[6]使用惩罚加权最小二乘(Penalized Weighted Least Square,PWLS)重建算法进行图像重建,针对图像空间像素间的几何相关性,采用马尔可夫随机场吉布斯函数进行数学建模,采用逐次超松弛迭代算法进行求解,取得了比较满意的图像效果。Xu等[7]把字典学习引入到LDCT重建算法中,通过交替迭代最小化方法最优化目标函数,这样不仅能够很好地去除图像噪声,还保留了图像更多的细节纹理信息结构。Tian等[8]将全变分(Total Variation,TV)能量最小化函数引入到图像空间中,由于TV正则项具有良好的边缘保持特性,因此LDCT重建算法获得了较好的重建图像。Zhu等[9]针对图像中含有的泊松噪声提出一种改进的TV模型,与传统的TV正则化模型相比,它不仅能去除泊松噪声,还能获得与原始图像更为接近的重建图像。

Chao等[10]把方差引入到图像重建中,形成一种新的基于图像灰度方差的扩散系数,提出一种新的自适应扩散模型,它克服了原始扩散模型容易产生阶梯效应和模糊图像边缘的缺点。Lou等[11]使用各向同性和各向异性扩散的加权TV正则化算法来进行图像降噪,相对于各向同性扩散模型,各向异性扩散模型能更好地保持图像的边缘信息。TV先验模型常用于统计迭代重建算法中,特别是在投影数据缺失的情况下,能够取得十分理想的重建效果[12-13]。TV正则化方法在去除图像噪声的同时能较好地保持图像的边缘细节信息,但是容易使图像的平滑区域产生阶梯效应。本文提出一种自适应加权全变分(Adaptive Weighted TV,AWTV)的低剂量CT统计迭代重建算法,首先把基于加权方差的扩散系数引入到传统TV模型中,然后把改进的自适应TV模型与PWLS重建算法相结合构成新的目标函数,再使用交替迭代重建方法对目标函数进行分步求解。

1 噪声模型

精确的投影数据和噪声建模是LDCT重建过程中至关重要的一个环节。LDCT图像的噪声表现在图像上往往是一些条形伪影和脉冲噪声,使得图像被噪声严重污染而无法看清细节信息。根据投影数据噪声模型的特点,Li等[14]认为经系统校准和对数变换后的LDCT投影数据近似地服从非平稳高斯分布,均值和方差之间呈非线性递增关系,其统计模型描述如下:

$\mathbf{\delta }_{i}^{2}={{\varepsilon }_{i}}\text{exp}\left( \frac{{{\mathbf{P}}_{i}}}{\eta } \right)$ (1)

其中:Piδi2分别表示第i个探测器上测得的投影数据的均值和方差;εi是自适应于不同探测器的参数;η表示用于描述扫描系统校准过程中的参数。

2 自适应加权全变分的LDCT统计迭代算法 2.1 惩罚加权最小二乘法

Fessler[15]首次提出加权最小二乘(Weighted Least Square,WLS)算法,其在图像域目标函数为:

$\Phi \left( \mathbf{u} \right)={{\left( \mathbf{y}-\mathbf{Gu} \right)}^{^{\operatorname{T}}}}{{\Sigma }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)$ (2)

其中:G表示系统矩阵; y是经系统校准和对数变换后的投影值;u表示待处理的图像;T表示转置运算符;Σ表示对角元素为δi2的对角矩阵,即第i个探测元上y值的方差。

WLS的目标函数是一个病态问题,最小化其目标函数往往得不到理想的结果。为解决这一问题,在WLS的目标函数中加入有平滑约束作用的惩罚项或正则项,从而得到惩罚加权最小二乘重建算法。PWLS目标函数可以改写为:

$\Phi \left( \mathbf{u} \right)={{\left( \mathbf{y}-\mathbf{Gu} \right)}^{\operatorname{T}}}{{\sum }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)+\beta R\left( \mathbf{u} \right)$ (3)

其中:(y-Gu)TΣ-1(y-Gu)表示WLS的目标函数;R(u)表示惩罚项; β表示系统的平滑参数。

2.2 自适应加权全变分模型 2.2.1 全变分模型

1992年,Rudin等[16]首次提出一种经典的图像降噪模型——TV模型,即为极小化如下的能量泛函:

$\underset{\mathbf{u}}{\mathop{\text{arg min}}}\,\left\{ {{\int_{\Omega }{\left| D\mathbf{l} \right|}}_{1}}dxdy+\frac{\lambda }{2}\int_{\Omega }{\left| \mathbf{u}-\mathbf{l} \right|_{2}^{2}dxdy} \right\}$ (4)

其中:Ω|Dl|表示l的全变分;Ω∈R是一个开区域。

图像降噪时,往往将式(4)改写为:

$\underset{\mathbf{u}}{\mathop{\text{arg min}}}\,\left\{ \int_{\Omega }{{{\left\| \nabla \mathbf{l} \right\|}_{1}}dxdy}+\frac{\lambda }{2}\int_{\Omega }{\left\| \mathbf{u}-\mathbf{l} \right\|_{2}^{2}\text{d}x\text{d}y} \right\}$ (5)

其中:${{\int_{\Omega }{\left\| \nabla \mathbf{l} \right\|}}_{1}}\text{d}x\text{d}y$表示l的全变分,称为图像l的先验项或正则项;Ωu-l22 dxdy表示数据保真项;λ为正则化参数,它在保持图像边缘信息和去除图像噪声的平衡方面起着十分重要的作用。

2.2.2 自适应加权全变分模型

TV 模型的解在有界变差空间中,所以它允许不连续解的存在,这样就能够在降噪的同时保持图像的边缘信息。但TV算法在强噪声的情况下,往往会使图像的平滑区域产生阶梯效应,并且不能很好地保护图像的纹理信息。因此,在原TV模型的基础上引入一个基于加权方差的加权因子,从而得到一个自适应加权TV模型,该先验模型能够根据图像噪声水平和局部结构特征来自适应地调节平滑的程度,这样不仅可以有效地去除TV模型产生的阶梯效应,还能很好地保护图像的边缘和细节纹理信息。

利用幅值相似度模糊函数m和空间邻近度模糊函数d来定义加权方差(Weighted Variance,WV)[17]

${{\mathbf{m}}_{i,j}}=\text{exp}\left( -\frac{{{\left( {{\mathbf{m}}_{i,j}}-{{\mathbf{m}}_{x,y}} \right)}^{6}}}{{{B}^{6}}} \right)$ (6)
${{\mathbf{d}}_{i,j}}=\text{exp}\left( -\frac{{{\left( i-x \right)}^{2}}+{{\left( j-y \right)}^{2}}}{{{D}^{2}}} \right)$ (7)

加权系数s表示如下:

$\mathbf{s}\left( i,j \right)={{\mathbf{m}}_{i,j}}\times {{\mathbf{d}}_{i,j}}$ (8)

因此,在图像3×3邻域内,基于加权方差的加权因子可表示为:

$\begin{align} & \mathbf{\omega }\left( \mathbf{l}\left( x,y \right) \right)=\mathbf{g}\left( \nabla \mathbf{l}_{t}^{i}\left( x,y \right)\cdot \mathbf{\sigma }_{t,W,N}^{2}\left( x,y \right) \right)= \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{1}{1+{{\left( \frac{\nabla \mathbf{l}_{t}^{i}\left( x,y \right)\cdot \mathbf{\sigma }_{t,W,N}^{2}\left( x,y \right)}{K} \right)}^{2}}} \\ \end{align}$ (9)
$\mathbf{\sigma }_{t,W,N}^{2}\left( x,y \right)=1+\frac{\mathbf{\sigma }_{t,W}^{2}\left( x,y \right)-\text{min}\mathbf{\sigma }_{t,W}^{2}}{Max\mathbf{\sigma }_{t,W}^{2}-\text{min}\mathbf{\sigma }_{t,W}^{2}}\times 254$ (10)
$\mathbf{\sigma }_{t,W}^{2}\left( x,y \right)=\frac{\sum{\mathbf{s}\left( i,j \right)\times \mathbf{\sigma }_{t}^{2}\left( x,y \right)}}{\sum{\mathbf{s}\left( i,j \right)}}$ (11)
$\mathbf{\sigma }_{t}^{2}\left( x,y \right)={{\left[ {{\mathbf{l}}_{t}}\left( x+i,y+j \right)-\overline{{{\mathbf{l}}_{\operatorname{t}}}}\left( i,j \right) \right]}^{2}}$ (12)

其中:x、y表示当前像素点;i、j表示像素的邻域;K表示控制扩散强度的系数;σt2(x,y)表示灰度方差;σt,W2(x,y)表示加权方差;σt,W,N2(x,y)表示归一化的加权方差;ω(l(x,y))表示基于加权方差的加权因子;lt(i,j)表示3×3窗口中各像素点的均值;B为控制幅值相似度的参数,D为控制空间邻近度的参数,其值分别为7和3;表示梯度。

自适应加权全变分模型的表达式可以改写为:

$\begin{align} & \underset{\mathbf{u}}{\mathop{\text{arg min}}}\,\left\{ \int_{\Omega }{\mathbf{\omega }\left( \mathbf{l}\left( x,y \right) \right){{\left\| \nabla \mathbf{l} \right\|}_{1}}\text{d}x\text{d}y}\text{+} \right. \\ & \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{\lambda }{2}\int_{\Omega }{\left\| \mathbf{u}-\mathbf{l} \right\|_{2}^{2}\text{d}x\text{d}y} \right\} \\ \end{align}$ (13)

ω(l(x,y))${{\left\| \nabla \mathbf{l} \right\|}_{1}}$都是基于图像的梯度来计算,因此ω(l(x,y))${{\left\| \nabla \mathbf{l} \right\|}_{1}}$是有效的,从而有效地验证了式(13)的正确性。

2.3 自适应加权全变分的低剂量CT统计迭代算法

通过使用RAWTV来表示自适应加权全变分的模型,则自适应加权全变分的LDCT重建算法的新目标函数可写为如下的形式:

$\Phi \left( \mathbf{u} \right)={{\left( \mathbf{y}-\mathbf{Gu} \right)}^{\operatorname{T}}}{{\Sigma }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)+\beta {{R}_{\text{AWTV}}}\left( \mathbf{u} \right)$ (14)

为了更好地求解式(14),引入一个和u有相同含义的中间变量f,则可转化为如下的约束优化问题:

$\left\{ \begin{align} & \underset{\mathbf{u}}{\mathop{\arg \min }}\,{{\left( \mathbf{y}-\mathbf{Gu} \right)}^{\operatorname{T}}}{{\sum }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)+\beta {{R}_{\text{AWTV}}}\left( \mathbf{f} \right) \\ & \left\| \mathbf{u}-\mathbf{f} \right\|_{2}^{2}=0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ \end{align} \right.$ (15)

利用Lagrange定理,将式(15)转化为如下的无约束优化问题[18]:

$\begin{align} & \underset{\mathbf{u}}{\mathop{\arg \min }}\,{{\left( \mathbf{y}-\mathbf{Gu} \right)}^{\text{T}}}{{\sum }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)+ \\ & \ \ \ \ \ \ \ \ \ \ \ \ \beta {{R}_{\text{AWTV}}}\left( \mathbf{f} \right)+\frac{\lambda }{2}\left\| \mathbf{u}-\mathbf{f} \right\|_{2}^{2} \\ \end{align}$ (16)

对于式(16),采用交替方向法进行优化,即转化为以下的两个问题进行求解:

问题1

$\mathbf{u}=\underset{\mathbf{u}}{\mathop{\text{arg min}}}\,{{\left( \mathbf{y}-\mathbf{Gu} \right)}^{\text{T}}}{{\sum }^{-1}}\left( \mathbf{y}-\mathbf{Gu} \right)+\frac{\lambda }{2}\left\| \mathbf{u}-\mathbf{f} \right\|_{2}^{2}$ (17)

问题2

$\mathbf{f}=\underset{\mathbf{f}}{\mathop{\text{arg min}}}\,\beta {{R}_{\text{AWTV}}}\left( \mathbf{f} \right)+\frac{\lambda }{2}{{\left\| \mathbf{u}-\mathbf{f} \right\|}^{2}}$ (18)

首先,利用可分离抛物面替代算法[19]求解问题1,可得:

$\begin{align} & \mathbf{u}_{\mathbf{j}}^{\mathbf{k+1}}=\mathbf{u}_{\mathbf{j}}^{\mathbf{k}}-\frac{\sum\limits_{i=1}^{M}{\left( \left( \frac{1}{\delta _{i}^{2}} \right){{G}_{ij}}\left( {{\left[ G{{u}^{k}} \right]}_{i}}-{{y}_{i}} \right) \right)}}{\sum\limits_{i=1}^{M}{\left( \left( \frac{1}{\delta _{i}^{2}} \right){{G}_{ij}}\sum\limits_{t=1}^{N}{{{G}_{it}}} \right)+\beta }} \\ & \ \ \ \ \ \ \ \ \ \text{+}\frac{\beta \left( u_{j}^{k}-f_{j}^{k} \right)}{\sum\limits_{i=1}^{M}{\left( \left( \frac{1}{\delta _{i}^{2}} \right){{G}_{ij}}\sum\limits_{t=1}^{N}{{{G}_{it}}} \right)+\beta }} \\ \end{align}$ (19)

其次,利用梯度下降流法来解决问题2:1)利用变分法计算f所满足的欧拉-拉格朗日方程:

$\text{div}\left( \frac{\mathbf{\omega }\left( \mathbf{f}\left( x \right) \right)\cdot \nabla \mathbf{f}}{{{\left\| \nabla \mathbf{f} \right\|}_{1}}} \right)-\lambda \left( \mathbf{u}-\mathbf{f} \right)=0$ (20)

2) 运用梯度下降流得到所满足的偏微分方程:

$\frac{\partial \mathbf{f}}{\partial {{\mathbf{f}}_{i,j}}}=\text{div}\left( \frac{\mathbf{\omega }\left( \mathbf{f}\left( x \right) \right)\cdot \nabla \mathbf{f}}{{{\left\| \nabla \mathbf{f} \right\|}_{1}}} \right)-\lambda \left( \mathbf{u}-\mathbf{f} \right)$ (21)

3) 采用数值计算方法求解此偏微分方程,则有:

$\begin{align} & \frac{\partial \mathbf{f}}{\partial {{\mathbf{f}}_{i,j}}}\approx \frac{\mathbf{\omega }\left( \mathbf{f}\left( i,j \right) \right)\cdot ({{\mathbf{f}}_{\operatorname{i},j}}-{{\mathbf{f}}_{\operatorname{i}-1,j}})}{\sqrt{\psi +{{({{\mathbf{f}}_{i,j}}-{{\mathbf{f}}_{i-1,j}})}^{2}}+{{({{\mathbf{f}}_{s,t}}-{{\mathbf{f}}_{i,j-1}})}^{2}}}} \\ & \ \ \ \ \ \ +\frac{\mathbf{\omega }\left( \mathbf{f}\left( i,j \right) \right)\cdot ({{\mathbf{f}}_{i,j}}-{{\mathbf{f}}_{i,j-1}})}{\sqrt{\psi +{{({{\mathbf{f}}_{i,j}}-{{\mathbf{f}}_{i-1,j}})}^{2}}+{{({{\mathbf{f}}_{s,t}}-{{\mathbf{f}}_{i,j-1}})}^{2}}}} \\ & \ \ \ \ \ -\frac{\mathbf{\omega }\left( \mathbf{f}\left( i,j \right) \right)\cdot ({{\mathbf{f}}_{i+1,j}}-{{\mathbf{f}}_{i,j}})}{\sqrt{\psi +{{({{\mathbf{f}}_{i+1,j}}-{{\mathbf{f}}_{i,j}})}^{2}}+{{({{\mathbf{f}}_{i+1,j}}-{{\mathbf{f}}_{i+1,j-1}})}^{2}}}} \\ & \ \ \ \ -\frac{\mathbf{\omega }\left( \mathbf{f}\left( i,j \right) \right)\cdot ({{\mathbf{f}}_{i,j+1}}-{{\mathbf{f}}_{i,j}})}{\sqrt{\psi +{{({{\mathbf{f}}_{i,j+1}}-{{\mathbf{f}}_{i,j}})}^{2}}+{{({{\mathbf{f}}_{i,j+1}}-{{\mathbf{f}}_{i-1,j+1}})}^{2}}}} \\ \end{align}$ (22)

其中:ψ是一个非常小的正参数,ψ=10-8

2.4 算法步骤

重建算法步骤概括如下:

步骤1 利用滤波反投影(Filtered Back Projection,FBP)重建算法,得到初始化图像f0;

步骤2 接着把初始化图像f0代入式(19),计算出迭代重建图像u;

步骤3 利用式(6)~(12)计算出基于加权方差的加权因子ω(f(x,y));

步骤4 根据式(22)采用自适应加权全变分模型计算出中间变量f;

步骤5 把求得的f代入式(19),更新重建图像u

重复步骤2~5,经过反复实验,选取视觉效果最好的重建结果作为最终的重建图像。

3 实验结果与分析

本文采用Shepp-Logan模型和数字骨盆模型来验证自适应加权全变分的低剂量CT统计迭代算法的有效性和可行性,将PWLS-AWTV算法分别与惩罚加权最小二乘(PWLS)算法[6]、惩罚加权最小二乘中值先验(PWLS-Median Prior,PWLS-MP)算法[20]以及惩罚加权最小二乘全变分(PWLS-TV)算法[12]进行实验结果的比较。图 1为两种模型数据。

图 1 两种模型数据(256mm×256mm)

通过向原始投影数据添加式(1)中的高斯噪声来进行实验的仿真模拟,重建算法在不同的噪声情况下会有不一样的处理结果。因此,本文研究η为22000,εi为200、300以及400时三种不同噪声情况下AWTV算法的重建图像。

为了更加清楚地看清AWTV重建图像的细节边缘区域,选取图像的三个感兴趣区(Region Of Interest,ROI)进行局部放大,图 2中标注了Shepp-Logan模型的三个ROI。

图 2 Shepp-Logan模型的三个ROI

图 3为不同噪声方差下的AWTV重建图像及其ROI。由图 3(b)可知,随着噪声水平的增加,AWTV算法会模糊图像的细节和有价值的边缘,图像质量变差。当εi为200,η为22000时重建图像的质量最好,图像在滤除噪声的同时能够很好地保护图像的边缘细节信息。

图 3 不同噪声方差下的AWTV重建图像及其ROI

下面分析εi为200,η为22000的噪声水平下,本文算法和其他三种对比重建算法的图像对比结果。图 4为Shepp-Logan模型各重建算法对比结果及其三个ROI对应的细节局部放大图。

图 4 Sheep-Logan模型各重建算法图像及其ROI

图 4可看出,PWLS算法可以去除LDCT图像中含有的部分条形伪影。PWLS-MP重建图像的质量较PWLS算法有较明显的改善,但边缘细节区域仍有条形伪影。PWLS-TV重建算法虽然抑制了条形伪影的产生,但PWLS-TV模型在降噪的同时导致虚假边缘的产生,从而引起阶梯效应。PWLS-AWTV重建算法不仅去除了PWLS-TV算法引入的阶梯伪影,还较好地保持图像的纹理和边缘信息,从视觉效果上来看,PWLS-AWTV算法明显优于其他三种重建算法。

图 5为数字骨盆模型的三个ROI,图 6为数字骨盆模型各重建算法对比结果图及其对应的三个ROI细节对比图。

图 5 数字骨盆模型的三个ROI
图 6 数字骨盆模型各重建算法对比结果及其ROI

图 6可知,PWLS算法虽然改善了重建图像中的条形伪影现象,但图像的细节边缘信息仍不清晰。PWLS-MP算法虽然去除了明显的条形伪影,但图像中仍有少量区域含有条形伪影,这样会影响医生作出精确的诊疗判断。PWLS-TV重建算法在消除条形伪影的同时带来阶梯效应,在图像的边缘细节区域阶梯伪影尤为严重。PWLS-AWTV算法在平滑图像的同时可以较好地保持图像的纹理和边缘信息。因此,本文算法在降噪能力、伪影抑制以及分辨率保持方面都优于其他的几种对比重建算法。

为了进一步更好地验证本文算法的优越性,图 7描绘了两种模型的原始图像和各重建算法在第125列上的侧面轮廓线的对比图。本文算法无论是在边缘区域还是在背景区域都更接近原始图像的侧面轮廓,并且图像噪声波动最小,同时还能够取得噪声抑制与边缘保持之间良好的平衡效果。

图 7 两种模型第125列侧面轮廓线的对比

为了更清楚地描述不同算法的各种质量参数,本文采用归一化均方距离(Normalized Mean Square Distance,NMSD)、归一化平均绝对距离(Normal Average Absolute Distance,NAAD)和峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)等对各重建算法进行定量的描述:

$NMSD=\sqrt{\frac{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{{{\left( \mathbf{r}\left( i,j \right)-\mathbf{q}\left( i,j \right) \right)}^{2}}}}}{\sum\limits_{i=1}^{M}{\sum\limits_{i=1}^{M}{{{\left( \mathbf{r}\left( i,j \right)-\mathbf{\bar{r}}\left( i,j \right) \right)}^{2}}}}}}$ (23)
$NAAD=\frac{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\left| \mathbf{r}\left( i,j \right)-\mathbf{q}\left( i,j \right) \right|}}}{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\left| \mathbf{r}\left( i,j \right) \right|}}}$ (24)
$\operatorname{PSNR}\text{=}10\text{lg}\left( \frac{{{255}^{2}}}{MAE} \right)$ (25)
$MAE=\frac{{{\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{\left[ \left( \mathbf{r}\left( i,j \right)-\mathbf{q}\left( i,j \right) \right) \right]}}}^{2}}}{\left( M\times N \right)}$ (26)

其中:M×N表示图像的大小;rq分别表示原始图像和重建图像;表示原始图像的均值。

表 1为Shepp-Logan模型和数字骨盆模型采用PWLS、PWLS-MP、PWLS-TV以及PWLS-AWTV算法重建的图像相对于原始图像的NMSD、NAAD、PSNR以及重建时间等质量参数的对比结果。由表 1可知,本文算法的NMSD和NAAD值都比其他的对比重建算法小,PSNR比其他算法大。表明本文算法重建出的图像与原始图像质量相当,相似度很高,因此本文算法具有很高的可行性。

表 1 Shepp-Logan和数字骨盆的模型各算法的客观评价参数
4 结语

本文提出一种自适应加权全变分的低剂量CT统计迭代算法,该算法先利用FBP算法获得LDCT重建图像,然后利用该图像作为初始化图像进行PWLS迭代重建,接着使用自适应加权全变分模型对得到的图像进行降噪处理,使用交替迭代重建的方法得到最终的重建图像。通过与其他重建算法进行仿真对比分析,无论是在主观的视觉效果方面,还是在客观的质量评价方面,本文算法不仅可以有效地去除条形伪影等噪声,而且保留了图像更多的细节信息,使得重建图像更加接近原始图像。本文方法对LDCT重建图像有不错的效果,接下来要对统计迭代重建算法作进一步的研究。

参考文献
[1] LUI D, CAMERON A, MODHAFAR A, et al. Low-dose computed tomography via spatially adaptive Monte-Carlo reconstruction[J]. Computerized Medical Imaging and Graphics, 2013, 37 (7/8) : 438-449. (0)
[2] RUST G F, AURICH V, REISER M. Noise dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20 mAs with nonlinear Gaussian filter chains [C]//Proceedings of SPIE 4683. Bellingham, WA: SPIE Press, 2002: 186-197. (0)
[3] CHEN Y, LI Y, YU W, et al. Joint-map tomographic reconstruction with patch similarity based mixture prior model[J]. Multiscale Modeling and Simulation, 2011, 9 (4) : 1399-1419. doi: 10.1137/100814184 (0)
[4] ZHANG H, MA J H, WANG J, et al. Statistical image reconstruction for low-dose CT using nonlocal means-based regularization[J]. Computerized Medical Imaging and Graphics, 2014, 38 : 423-435. doi: 10.1016/j.compmedimag.2014.05.002 (0)
[5] ZHANG H, MA J H, WANG J, et al. Statistical image reconstruction for low-dose CT using nonlocal means-based regularization, Part II: an adaptive approach[J]. Computerized Medical Imaging and Graphics, 2015, 43 : 26-35. doi: 10.1016/j.compmedimag.2015.02.008 (0)
[6] WANG J, LI T F, LU H B, et al. Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography[J]. IEEE Transactions on Medical Imaging, 2006, 25 (10) : 1272-1283. doi: 10.1109/TMI.2006.882141 (0)
[7] XU Q, YU H Y, MOU X Q, et al. Low-dose X-ray CT reconstruction via dictionary learning[J]. IEEE Transactions on Medical Imaging, 2012, 31 (9) : 1682-1697. doi: 10.1109/TMI.2012.2195669 (0)
[8] TIAN Z, JIA X, YUAN K H, et al. Low-dose CT reconstruction via edge-preserving total variation regularization[J]. Physics in Medicine and Biology, 2011, 56 (18) : 5949-5967. doi: 10.1088/0031-9155/56/18/011 (0)
[9] ZHU Y, ZHAO M L, ZHAO Y S, et al. Noise reduction with low-dose CT data based on a modified ROF model[J]. Optics Express, 2012, 20 (16) : 17987-18004. doi: 10.1364/OE.20.017987 (0)
[10] CHAO S M, TSAI D M. An improved anisotropic diffusion model for detail- and edge-preserving smoothing[J]. Pattern Recognition Letters, 2010, 31 (13) : 2012-2023. doi: 10.1016/j.patrec.2010.06.004 (0)
[11] LOU Y F, ZENG T Y, OSHER S, et al. A weighted difference of anisotropic and isotropic total variation model for image processing[J]. SIAM Journal on Imaging Sciences, 2015, 8 (3) : 1798-1823. doi: 10.1137/14098435X (0)
[12] TANG J, NETT B E, CHEN G H. Performance comparison between total variation-based compressed sensing and statistical iterative reconstruction algorithms[J]. Physics in Medicine & Biology, 2009, 54 (19) : 5781-5804. (0)
[13] XU Q, MOU X, WANG G, et al. Statistical interior tomography[J]. IEEE Transactions on Medical Imaging, 2011, 30 (5) : 1116-1128. doi: 10.1109/TMI.2011.2106161 (0)
[14] LI T F, LI X, WANG J, et al. Nonlinear sinogram smoothing for low-dose X-ray CT[J]. IEEE Transactions on Nuclear Science, 2004, 51 (5) : 2505-2513. doi: 10.1109/TNS.2004.834824 (0)
[15] FESSLER J A. Penalized weighted least-squares image reconstruction for positron emission tomography[J]. IEEE Transactions on Medical Imaging, 1994, 13 (2) : 290-300. doi: 10.1109/42.293921 (0)
[16] RUDIN L, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D — Nonlinear Phenomena, 1992, 60 (1/2/3/4) : 259-268. (0)
[17] KIANI M, SEID H P. State estimation of nonlinear dynamic systems using weighted variance-based adaptive particle swarm optimization[J]. Applied Soft Computing, 2015, 34 (C) : 1-17. (0)
[18] LU X Q, SUN Y, YUAN Y, et al. Image reconstruction by an alternating minimization[J]. Neuro Computing, 2011, 74 (5) : 661-670. (0)
[19] ELBKRI I, FESSLER J. Statistical image reconstruction for polyene-rgetic X-ray computed tomography[J]. IEEE Transactions on Medical Imaging, 2002, 21 : 89-99. doi: 10.1109/42.993128 (0)
[20] HSIAO T, RANGARAJAN A, GINDI G. A new convex edge-preserving median prior with applications to tomography[J]. IEEE Transactions on Medical Imaging, 2003, 22 (5) : 580-585. doi: 10.1109/TMI.2003.812249 (0)
[21] 钱姗姗, 黄静, 马建华, 等. 基于投影数据非单调性全变分恢复的低剂量CT重建[J]. 电子学报, 2011, 39 (7) : 1702-1707. ( QIAN S S, HUANG J, MA J H, et al. Nonmonotone total variation minimization based projection restoration for low-dose CT reconstruction[J]. Acta Electronica Sinica, 2011, 39 (7) : 1702-1707. ) (0)
[22] 王丽艳, 韦志辉. 低剂量CT 的线性Bregman迭代重建算法[J]. 电子与信息学报, 2013, 35 (10) : 2418-2424. ( WANG L Y, WEI Z H. Linearized bregman iterations for low-dose CT reconstruction[J]. Journal of Electronics and Information, 2013, 35 (10) : 2418-2424. ) (0)