CN106127825B - 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 - Google Patents
一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 Download PDFInfo
- Publication number
- CN106127825B CN106127825B CN201610428470.7A CN201610428470A CN106127825B CN 106127825 B CN106127825 B CN 106127825B CN 201610428470 A CN201610428470 A CN 201610428470A CN 106127825 B CN106127825 B CN 106127825B
- Authority
- CN
- China
- Prior art keywords
- image
- projection data
- dose
- reconstruction
- ray
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 38
- 230000009466 transformation Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000005457 optimization Methods 0.000 claims description 7
- 238000012937 correction Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 description 8
- 238000011084 recovery Methods 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000013170 computed tomography imaging Methods 0.000 description 4
- 238000003759 clinical diagnosis Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000002591 computed tomography Methods 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000002401 inhibitory effect Effects 0.000 description 2
- 238000012805 post-processing Methods 0.000 description 2
- 238000004088 simulation Methods 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 208000026350 Inborn Genetic disease Diseases 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000011281 clinical therapy Methods 0.000 description 1
- 238000002247 constant time method Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000002059 diagnostic imaging Methods 0.000 description 1
- 230000010339 dilation Effects 0.000 description 1
- 208000016361 genetic disease Diseases 0.000 description 1
- 231100000405 induce cancer Toxicity 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 208000032839 leukemia Diseases 0.000 description 1
- 239000002184 metal Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000013179 statistical model Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T11/00—2D [Two Dimensional] image generation
- G06T11/003—Reconstruction from projections, e.g. tomography
- G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
一种基于广义惩罚加权最小二乘的X射线CT图像重建方法,包括如下步骤:(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据qe;(2)对(1)中获取的投影数据qe进行系统校正和对数变换后得到投影数据y;对y进行逐个数据点上的方差估计,获得每个点上的方差集合(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型;(4)对(3)中的低剂量CT重建模型进行求解,建立算法的全局收敛性以收敛点作为解,再进行图像重建,得到最终的CT重建图像。本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。
Description
技术领域
本发明涉及医学影像的计算机断层成像技术领域,特别是涉及一种基于广义惩罚加权最小二乘的X射线CT图像重建方法。
背景技术
X射线计算机断层成像技术是现代医学影像学的杰出代表,已经广泛应用于临床诊断和治疗。然而偏高的X射线照射剂量可诱发癌症、白血病或者其他遗传性疾病。如何最大限度地降低X射线使用剂量已经成为医学CT成像领域研究的关键技术。
当前,临床上采用的低剂量CT扫描方式虽然可以在一定程度上减少X射线照射剂量,但低剂量的投影数据将受到量子噪声和电子噪声污染,而采用滤波反投影算法重建的CT图像质量严重退化,直接影响临床诊断的准确性。
对于这一问题有两种传统的后处理解决方法:一是直接对低剂量CT图像进行滤波,以减少图像的噪声和伪影,属于图像后处理技术;二是根据投影数据满足的统计学规律,完成基于统计的CT图像迭代重建。第一种的处理技术,直接减少图像的噪声和伪影。图像后处理技术因其简单且易于操作已在低剂量CT图像恢复中得到广泛的研究和应用。由于低剂量CT图像中噪声和伪影分布的复杂性,使得高精度的滤波器设计困难极大。第二种通过CT系统建模,构建图像重建模型,通过优化目标函数实现图像重建。相对于经典的FBP算法,迭代重建算法通过系统建模(系统光学模型和系统统计模型)对CT成像几何、X射线的能谱特性、射束硬化效应、散射和噪声特性进行准确描述,而且易于加入先验信息约束,因此特别适合低剂量CT图像优质重建。统计迭代重建在抑制图像噪声和伪影以及提高空间分辨率等方面都有上佳表现,是当前低剂量CT成像领域的研究热点。然而,由于统计迭代重建需要反复进行投影与反投影运算,且CT图像数据量庞大,导致CT图像重建速度特别慢,难以满足临床中实时交互的需求。
不同于上述两种方法的另一种策略是,根据投影数据的噪声统计特性建立数据恢复模型,得到更接近于理想值的投影数据后采用FBP算法重建出CT图像。该类方法不仅可以避免统计迭代重建方法在计算速度上的劣势,而且考虑了投影数据的噪声统计特性,去除噪声的同时可以较好地保持图像的空间分辨率。诸多研究中,一般都是将投影数据的噪声看作高斯或者泊松噪声来处理,大量的临床低剂量CT实验和理论分析,提出投影数据(对数变换后)的噪声近似满足均值与方差呈非线性关系的高斯分布,基于这一重要特点,多种低剂量CT投影数据恢复的惩罚加权最小二乘(penalized weighted least-squares,简记为PWLS)方法相继被提出,例如基于K—L变换的PWLS方法,多尺度PWLS方法等。这种处理方法相比传统的滤波反投影方法能够在一定程度上改善CT重建图像的质量,但重建的图像在对比度和信噪比等方面仍然存在一定的不足,从而影响临床诊断结果。
因此,针对现有技术不足,提供一种可以有效地抑制低剂量CT图像中的噪声和伪影,同时较好地保持图像的空间分辨率的基于广义惩罚加权最小二乘(generalizedpenalized weighted least-squares,简记为GPWLS)的X射线CT图像重建方法以克服现有技术不足甚为必要。
发明内容
本发明的目的在于避免现有技术的不足之处而提供一种基于广义惩罚加权最小二乘的X射线CT图像重建方法,该方法可以有效地抑制低剂量CT图像中的噪声和伪影,同时较好地保持图像的空间分辨率。
本发明的上述目的通过如下技术手段实现。
提供一种基于广义惩罚加权最小二乘的X射线CT图像重建方法,包括如下步骤:
(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据qe;
(2)对(1)中获取的投影数据qe进行系统校正和对数变换后得到投影数据y,y={y1,y2,...,yM},i=1,2,……,M,yi是投影数据y的第i个分量,M是正整数,M是投影数据y的分量的个数;对y进行逐个数据点上的方差估计,获得每个点上的方差集合
(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型;
(4)对(3)中的低剂量CT重建模型进行求解,建立算法的全局收敛性以收敛点作为解,再进行图像重建,得到最终的CT重建图像。
优选的,上述步骤(1)中获取CT设备的系统参数包括X射线入射光子强度I0和系统电子噪声的方差
优选的,上述步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。
优选的,上述步骤(3)中建立的低剂量CT重建模型为:
其中,“T”表示转置运算,I为单位矩阵,
β1、β2为两个超参数,β1>0,β2>0,
R(q)是惩罚项,R(q)的表达式为:
其中,Ni为投影数据中第i个像素的一阶邻域,qi为投影数据中第i个像素值,qm为投影数据中第i个像素的邻域的像素值,ωim表示qi的一阶邻域内的qm所占的权重,并且ωim在水平方向取值为1,在垂直方向取值为0.25。
优选的,上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。
优选的,上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解,具体是:
引入变量p,得到
其中,
根据式(Ⅳ),将式(Ⅰ)写成如下形式:
其中,J(p,q)是关于p和q的函数;
(4.1)给定初始点q0=y,令k=0,k是迭代次数,k+1表示第k+1次迭代;
(4.2)使用交替方向法依次求pk+1和qk+1,
其中,SG(qk)表示式(Ⅴ)中对p的优化解,SR(pk+1)表示式(Ⅴ)中对q的优化解;
(4.3)判断是否满足迭代终止条件:如果满足迭代终止条件,则停止迭代,进入步骤(4.5),否则进入步骤(4.4);
(4.4)令k=k+1,进入步骤(4.2);
(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数,对Ⅵ(a)式求导并令其为零,得到Ⅵ(a)的解为
对于式Ⅵ(b),使用Gauss-Seidel算法求解,其对应的迭代格式为其中 表示像素qi左上方的邻域,表示像素qi右上方的邻域,Ni表示像素的一阶邻域。
优选的,上述步骤(4)中,具体通过滤波反投影方法进行图像重建,得到最终的CT重建图像。
本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法,包括如下步骤:(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据qe;(2)对(1)中获取的投影数据qe进行系统校正和对数变换后得到投影数据y,y={y1,y2,...,yM},i=1,2,L M,yi是投影数据y的第i个元素,M是正整数,M是投影数据y的分量的个数;对y进行逐个数据点上的方差估计,获得每个点上的方差集合(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型;(4)对(3)中的低剂量CT重建模型进行求解,建立算法的全局收敛性以收敛点作为解,再进行图像重建,得到最终的CT重建图像。该基于广义惩罚加权最小二乘的X射线CT图像重建方法,构造广义加权最小二乘保真项,并且将惩罚先验信息引入投影数据的恢复过程中,从而达到抑制噪声的目的,最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。
附图说明
利用附图对本发明作进一步的说明,但附图中的内容不构成对本发明的任何限制。
图1是本发明一种基于广义惩罚加权最小二乘的X射线CT图像重建方法的流程图。
图2为本发明实施例2中用到的数值体模,其中(a)为Clock体模,(b)为Shepp-Logan体模。
图3为数值Clock体模的重建结果,其中(a)为FBP算法(Ramp窗)的重建结果;(b)为FBP算法(Hanning窗)的重建结果;(c)为PWLS方法的重建结果;(d)为本发明方法的重建结果。
图4为数值Shepp-Logan体模的重建结果,其中(a)为FBP算法(Ramp窗)的重建结果;(b)为FBP算法(Hanning窗)的重建结果;(c)为PWLS方法的重建结果;(d)为本发明方法的重建结果。
图5为对应图4中不同方法重建图像的水平剖面图,其中(a)为FBP算法(Ramp窗)重建图像的水平剖面图;(b)为FBP算法(Hanning窗)重建图像的水平剖面图;(c)为PWLS方法重建图像的水平剖面图;(d)为GPWLS方法重建图像的水平剖面图。
图6为对应图4中不同方法重建图像感兴趣区域的局部放大图,其中(a)为FBP算法(Ramp窗)重建图像感兴趣区域的局部放大图;(b)为FBP算法Hanning窗)重建图像感兴趣区域的局部放大图;(c)为PWLS方法重建图像感兴趣区域的局部放大图;(d)为GPWLS方法重建图像感兴趣区域的局部放大图。
具体实施方式
结合以下实施例对本发明作进一步描述。
实施例1。
一种基于广义惩罚加权最小二乘的X射线CT图像重建方法,如图1所示,包括如下步骤。
(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据qe;
(2)对(1)中获取的投影数据qe进行系统校正和对数变换后得到投影数据y,y={y1,y2,...,yM},i=1,2,……,M,yi是投影数据y的第i个分量,M是正整数,M是投影数据y的分量的个数;对y进行逐个数据点上的方差估计,获得每个点上的方差集合
(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型,简称构建GPWLS模型;
(4)对(3)中的低剂量CT重建模型进行求解,建立算法的全局收敛性以收敛点作为解,再通过滤波反投影方法进行图像重建,得到最终的CT重建图像。
具体的,上述步骤(1)中获取CT设备的系统参数包括X射线入射光子强度I0和系统电子噪声的方差
具体的,上述步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。
具体的,上述步骤(3)中建立的低剂量CT重建模型为:
其中,“T”表示转置运算,I为单位矩阵,
β1、β2为两个超参数,β1>0,β2>0,
R(q)是惩罚项,R(q)的表达式为:
其中,Ni为投影数据中第i个像素的一阶邻域,qi为投影数据中第i个像素值,qm为投影数据中第i个像素的邻域的像素值,ωim表示qi的一阶邻域内的qm所占的权重,并且ωim在水平方向取值为1,在垂直方向取值为0.25。
具体的,上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。
步骤(4)是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解,具体是:
引入变量p,得到
其中,
根据式(Ⅳ),将式(Ⅰ)写成如下形式:
其中,J(p,q)是关于p和q的函数;
(4.1)给定初始点q0=y,令k=0,k是迭代次数,k+1表示第k+1次迭代;
(4.2)使用交替方向法依次求pk+1和qk+1,
其中,SG(qk)表示式(Ⅴ)中对p的优化解,SR(pk+1)表示式(Ⅴ)中对q的优化解;
(4.3)判断是否满足迭代终止条件:如果满足迭代终止条件,则停止迭代,进入步骤(4.5),否则进入步骤(4.4);
(4.4)令k=k+1,进入步骤(4.2);
(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数,对Ⅵ(a)式求导并令其为零,得到Ⅵ(a)的解为
对于式Ⅵ(b),使用Gauss-Seidel算法求解,其对应的迭代格式为其中 表示像素qi左上方的邻域,表示像素qi右上方的邻域,Ni表示像素的一阶邻域。
下面给出GPWLS方法的全局收敛性分析。根据最优化理论,非扩张渐近正则算子Γ以及Γ的不动点解集非空,则对任意初始点q0都存在序列{Γ(qk)}弱收敛到Γ的不动点。因此,要分析GPWLS方法的收敛性,首先应该证明是Γ非扩张渐近正则算子,以及目标函数的解集非空,下面给出几个相关的概念。
定义1如果对任意的x1,x2∈RN都有||Ρ(x1)-Ρ(x2)||≤||x1-x2||,则Ρ称为非扩张算子。如果存在非扩张算子A以及κ∈(0,1)使得Ρ=(1-κ)I+κA,则称是κ—平均非扩张算子。
定义2设C是Banach空间X上的一个闭凸集,如果对任意的x∈C,F:C→C满足:
则称F:C→C是渐近正则的。
定义3设函数如果至少存在一点使得并且对任意x∈X都有则称在X上是适定的。如果对X中满足||xk||→∞的任意子列都则称在X上是强制的。
在给出GPWLS方法的收敛性定理之前,首先给出如下几个引理。
引理1设是凸的下半连续的函数以及β>0,定义如下:对任意x,定义则S为—平均非扩张算子。
引理2公式Γ(·)=SR(SG(·))定义的算子Γ是非扩张的。
证明:由于(Σ-1+β1I)是对称正定矩阵且其特征值大于β1,则对任意的u,v有:
||β1(Σ-1+β1I)-1(u-v)||2≤||u-v||2;
因此,
故引理得证。
引理3公式qk+1=Γ(qk)产生的点列{qk}满足
证明:在点(pk,qk)第一个分量pk+1处对目标函数进行泰勒展开可以得到:
由于,
pk+1=arg minp J(p,qk),
则,
进一步,
由于Σ-1+β1I对称正定且其特征值大于β1,则:
由于SR(q)为非扩张算子,则:
可以得到:
因此收敛,则
故引理得证。
引理4对任意初始值q0,设{qk}是由qk+1=Γ(qk)产生的序列,则Γ为渐近正则算子。
证明:由Γ的定义,则:
故引理得证。
引理5由公式(Ⅴ)定义的函数J(p,q)是强制的。
证明:由R(q)的定义,我们可以得到:
其中,Dh,Dv分别为水平和垂直方向的差分矩阵,注意到:
其中,进一步,
注意到上述矩阵是满秩,因此,当时,J(p,q)→∞,
故引理得证。
引理6公式Γ(·)=SR(SG(·))定义的Γ的不动点解集非空。
证明:由于J是强制的,因此J的极小值构成的集合非空,设(p*,q*)是J(x,u)的一个极小值即,
因此,我们可以得到
则q*=SR(SG(q*))=Γ(q*),即q*是Γ的一个不动点。故引理得证。
由于目标函数J是严格凸的,因此,Γ的不动点即是J的全局最优点。根据最优化理论序列{qk}收敛到J的不动点,即极小值点。下面给出GPWLS算法的全局收敛。
定理1给定初始值q0∈Rn,设{qk}是由qk+1=Γ(qk)产生的序列,则qk收敛到J的不动点即极小值点。
本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法,验证了算法的全局收敛性。本发明构造广义加权最小二乘保真项,并且将惩罚先验信息引入投影数据的恢复过程中,从而达到抑制噪声的目的,最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明,本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。
实施例2。
采用图1所示的Shepp-Logan数字体模图像和Clock数字体模图像作为本发明的计算机仿真实验对象。体模图像大小设定位512×512,模拟CT机的X射线源到旋转中心和探测器的距离分别为570mm和1040mm,旋转角在[0,2π]间采样值为1160,每个采样角对应672个探测器单元,探测器单元的大小为1.407mm。通过CT系统仿真生成大小为1160×672的投影数据,其中X射线的入射光子强度I0为5.0×104。在实际的CT数据采集中,投影数据和系统参数即入射光子强度I0和系统电子噪声的方差均可以直接获取。
对步骤1中模拟生成的CT投影数据和系统参数和I0,通过步骤2系统校正和对数变换后得到的投影数据y,以及方差采用如下公式:
在本实例中,I0为第i个数据点的X射线入射光子强度,即I0为5.0×104。为系统电子噪声的方差。式(Ⅱ)中Ni是投影数据中第i个像素的一阶邻域,ωim在水平方向取值为1,在垂直方向取值为0.25。不同于图像域的二阶邻域相似性,投影数据的相似性仅存在于相邻探测器单元和相邻投影角度这两个维度之间,而且前者的相似概率大于后者,因此水平与垂直方向使用不同的权重。根据步骤3利用本发明所示模型(Ⅰ)对待估计投影数据的理想值q进行求解。
本发明使用(Ⅰ)的广义加权最小二乘代替加权最小二乘。注意到,当β1充分大时,广义加权最小二乘就等价于加权最小二乘。进一步,引入一个新的变量p,可以得到:
根据(Ⅳ)式,可以将式(Ⅰ)写成如下形式:
求解式(Ⅴ)中的目标函数,本发明使用交替方向法依次求p和q,
SG(qk)表示式子Ⅴ中对p的优化解,SR(pk+1)表示式子Ⅴ中q的优化解。
给定初始点p0,产生迭代序列:p0,p1,q1,p2,q2,…,pk,qk,…。
满足如下的关系:
qk+1=SR(SG(qk))……(Ⅶ);
可以转化为:
qk+1=T(qk)……(Ⅷ);
由于Ⅵ(a)式中的目标函数为光滑的二次凸函数,对Ⅵ(a)式求导并令其为零,可以得到Ⅵ(a)的解为对于Ⅵ(b)式,使用Gauss-Seidel算法求解,其对应的迭代格式为其中 表示像素qi左上方的邻域,表示像素qi右上方的邻域,Ni表示像素的一阶邻域。
为了验证本发明所示方法的效果,对系统校正和对数变换后的投影数据建立广义惩罚加权最小二乘(GPWLS)模型,然后对估计的投影数据q采用传统的滤波反投影(FBP)算法进行重建,得到重建图像。为了对比本发明技术方案的效果,将本发明方法与Ramp窗FBP算法、Hanning窗FBP以及(PWLS)重建图像的影像进行比较,图3给出了Clock体模不同方法的重建结果。图3(a)为Ramp窗FBP算法重建的结果;图3(b)为Hanning窗FBP算法重建的结果;(c)为PWLS方法重建的结果;(d)为GPWLS方法重建的结果。从视觉效果评价来看,GPWLS方法在抑制图像噪声和条形伪影方面具有更佳的表现。为了进一步定量分析GPWLS方法,计算重建图像的信噪比(SNR),均方误差(NMSE)和以及重建图像中两个感兴趣区域(如图2(a)中所示)的对比度(CNR)。SNR,NMSE以及CNR的计算公式如下:
其中,μ(i)表示图像μ的第i个像素,表示μ的均值,μture(i)表示真实的体模图像μture的第i个像素,表示感兴趣(ROI)区域内图像μ的均值,表示感背景(Background)区域内图像μ的均值,表示ROI区域内图像μ的方差,表示Background区域内图像μ的方差。本发明中选取的ROI区域大小为20×20。由表1的信噪比和均方误差分析可以看出,与其他三种方法相比本发明在提高重建图像的信噪比,降低重建图像均方误差方面均有更佳表现。此外,由表2可以看出,在高对比度区域ROI1,四种方法重建图像的对比度相差不大。然而,在低对比度区域ROI2,GPWLS方法重建的图像具有最高的对比度。
表1图3中重建图像的信噪比和均方误差
表2图3中重建图像的对比度
图4为低剂量的Shepp-Logan投影数据由不同方法重建的结果。图4(a)为Ramp窗FBP算法重建的结果;图4(b)为Hanning窗FBP算法重建的结果;4(c)为PWLS方法重建的结果;4(d)为GPWLS方法重建的结果。从图4可以看出PWLS方法和GPWLS两种方法较好地抑制了图像中的条形伪影。从视觉评价来看,图4(d)中CT图像的噪声水平比图4(c)中CT图像的噪声水平低。
图5是重建结果的水平剖面图,从图中可以看出通过本发明重建的图像更接近真实的体模图像。为了更清晰的比较各种方法对噪声和伪影的抑制效果,图6给出了对应图4中重建结果的感兴趣区域的局部放大图。从图6可以看出,本发明在噪声去除、条形伪影抑制以及结构信息保持方面都比其他方法具有更佳的表现。
本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法,构造广义加权最小二乘保真项,并且将惩罚先验信息引入投影数据的恢复过程中,从而达到抑制噪声的目的,最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明,本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。
最后应当说明的是,以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制,尽管参照较佳实施例对本发明作了详细说明,本领域的普通技术人员应当理解,可以对本发明的技术方案进行修改或者等同替换,而不脱离本发明技术方案的实质和范围。
Claims (6)
1.一种基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于,包括如下步骤:
(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据qe;
(2)对(1)中获取的投影数据qe进行系统校正和对数变换后得到投影数据y,y={y1,y2,...,yM},i=1,2,……,M,yi是投影数据y的第i个分量,M是正整数,M是投影数据y的分量的个数;对y进行逐个数据点上的方差估计,获得每个点上的方差集合
(3)根据(2)中得到的方差建立广义惩罚加权最小二乘的低剂量CT重建模型;
(4)对(3)中的低剂量CT重建模型进行求解,建立算法的全局收敛性以收敛点作为解,再进行图像重建,得到最终的CT重建图像;
所述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解,具体是:
引入变量p,得到
其中,β1为超参数,β1>0,I为单位矩阵,
根据式(Ⅳ),将式(Ⅰ)写成如下形式:
其中,J(p,q)是关于p和q的函数,
β2为超参数,β2>0,R(q)是惩罚项;
(4.1)给定初始点q0=y,令k=0,k是迭代次数,k+1表示第k+1次迭代;
(4.2)使用交替方向法依次求pk+1和qk+1,
其中,SG(qk)表示式(Ⅴ)中对p的优化解,SR(pk+1)表示式(Ⅴ)中对q的优化解;q为投影数据中像素值;
(4.3)判断是否满足迭代终止条件:如果满足迭代终止条件,则停止迭代,进入步骤(4.5),否则进入步骤(4.4);
(4.4)令k=k+1,进入步骤(4.2);
(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数,对Ⅵ(a)式求导并令其为零,得到Ⅵ(a)的解为
2.根据权利要求1所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于,步骤(1)中获取CT设备的系统参数包括X射线入射光子强度I0和系统电子噪声的方差
3.根据权利要求2所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于,步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。
4.根据权利要求1至3任意一项所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于:
所述步骤(3)中建立的低剂量CT重建模型为:
其中,“T”表示转置运算,I为单位矩阵,
R(q)的表达式为:
其中,Ni为投影数据中第i个像素的一阶邻域,qi为投影数据中第i个像素值,qm为投影数据中第i个像素的邻域的像素值,ωim表示qi的一阶邻域内的qm所占的权重,并且ωim在水平方向取值为1,在垂直方向取值为0.25。
5.根据权利要求4所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于,所述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。
6.根据权利要求5所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法,其特征在于,所述步骤(4)中,具体通过滤波反投影方法进行图像重建,得到最终的CT重建图像。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610428470.7A CN106127825B (zh) | 2016-06-15 | 2016-06-15 | 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610428470.7A CN106127825B (zh) | 2016-06-15 | 2016-06-15 | 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106127825A CN106127825A (zh) | 2016-11-16 |
CN106127825B true CN106127825B (zh) | 2019-12-03 |
Family
ID=57470374
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610428470.7A Active CN106127825B (zh) | 2016-06-15 | 2016-06-15 | 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106127825B (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106683146B (zh) * | 2017-01-11 | 2021-01-15 | 上海联影医疗科技股份有限公司 | 一种图像重建方法和图像重建算法的参数确定方法 |
US10327727B2 (en) * | 2017-05-12 | 2019-06-25 | Varian Medical Systems, Inc. | Automatic estimating and reducing scattering in computed tomography scans |
CN111968192A (zh) * | 2020-06-29 | 2020-11-20 | 深圳先进技术研究院 | 一种ct图像的构建方法、ct设备以及存储介质 |
WO2022000192A1 (zh) * | 2020-06-29 | 2022-01-06 | 深圳先进技术研究院 | 一种ct图像的构建方法、ct设备以及存储介质 |
CN112489153B (zh) * | 2020-11-26 | 2024-06-21 | 深圳先进技术研究院 | 一种图像重建方法及应用 |
WO2022109928A1 (zh) * | 2020-11-26 | 2022-06-02 | 深圳先进技术研究院 | 一种图像重建方法及应用 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103810734A (zh) * | 2014-02-28 | 2014-05-21 | 南方医科大学 | 一种低剂量x射线ct投影数据恢复方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US9036885B2 (en) * | 2012-10-28 | 2015-05-19 | Technion Research & Development Foundation Limited | Image reconstruction in computed tomography |
-
2016
- 2016-06-15 CN CN201610428470.7A patent/CN106127825B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103810734A (zh) * | 2014-02-28 | 2014-05-21 | 南方医科大学 | 一种低剂量x射线ct投影数据恢复方法 |
Non-Patent Citations (2)
Title |
---|
低剂量 CT 成像与稀疏角度重建研究;刘文磊;《中国优秀硕士学位论文全文数据库(信息科技辑)》;20140315(第3期);正文第27页 * |
基于变分正则化的低剂量CT成像方法研究;牛善洲;《中国博士学位论文全文数据库(医药卫生科技辑)》;20160315(第3期);正文第29-30页,第35-36页,第46页,第52-53页 * |
Also Published As
Publication number | Publication date |
---|---|
CN106127825A (zh) | 2016-11-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN106127825B (zh) | 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 | |
Lu et al. | Iterative reconstruction of low-dose CT based on differential sparse | |
Hu et al. | Artifact correction in low‐dose dental CT imaging using Wasserstein generative adversarial networks | |
CN103810734B (zh) | 一种低剂量x射线ct投影数据恢复方法 | |
CN104751429B (zh) | 一种基于字典学习的低剂量能谱ct图像处理方法 | |
CN103810733B (zh) | 一种稀疏角度x射线ct图像的统计迭代重建方法 | |
CN105321155A (zh) | 一种cbct图像环形伪影消除方法 | |
CN101980302A (zh) | 投影数据恢复导引的非局部平均低剂量ct重建方法 | |
CN103810735A (zh) | 一种低剂量x射线ct图像统计迭代重建方法 | |
CN103413280A (zh) | 一种低剂量x射线ct图像重建方法 | |
Fang et al. | Removing ring artefacts for photon-counting detectors using neural networks in different domains | |
CN114187235B (zh) | 一种对伪影不敏感的医学图像的形变场提取方法和配准方法和装置 | |
CN105844678A (zh) | 基于全广义变分正则化的低剂量x射线ct图像重建方法 | |
Liu et al. | DFSNE-Net: Deviant feature sensitive noise estimate network for low-dose CT denoising | |
CN111080736B (zh) | 一种基于稀疏变换的低剂量ct图像重建方法 | |
Wu et al. | Data-iterative optimization score model for stable ultra-sparse-view CT reconstruction | |
Kaviani et al. | Image reconstruction using UNET-transformer network for fast and low-dose PET scans | |
CN109658464B (zh) | 基于加权核范数极小的稀疏角ct图像重建方法 | |
CN112116677B (zh) | 一种基于低维流形先验的低剂量ct重建方法 | |
Lee et al. | Four-dimensional CBCT reconstruction based on a residual convolutional neural network for improving image quality | |
Kamasak | Clustering dynamic PET images on the Gaussian distributed sinogram domain | |
Feng et al. | Multi-Dimensional Spatial Attention Residual U-Net (Msaru-Net) for Low-Dose Lung Ct Image Restoration | |
Liang et al. | A model-based deep learning reconstruction for X-Ray CT | |
Song et al. | Ultrasonic image processing based on fusion super-resolution reconstruction of familiar models | |
CN112489153B (zh) | 一种图像重建方法及应用 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
CB02 | Change of applicant information | ||
CB02 | Change of applicant information |
Address after: 341000 Ganzhou economic and Technological Development Zone, Jiangxi Applicant after: Gannan Normal University Address before: 341000 Ganzhou economic and Technological Development Zone, Jiangxi, GanNan Normal University Applicant before: Gan Nan Normal College |
|
GR01 | Patent grant | ||
GR01 | Patent grant |