CN106780641B - 一种低剂量x射线ct图像重建方法 - Google Patents

一种低剂量x射线ct图像重建方法 Download PDF

Info

Publication number
CN106780641B
CN106780641B CN201611002595.XA CN201611002595A CN106780641B CN 106780641 B CN106780641 B CN 106780641B CN 201611002595 A CN201611002595 A CN 201611002595A CN 106780641 B CN106780641 B CN 106780641B
Authority
CN
China
Prior art keywords
statistical
image
ray
data
dose
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
Application number
CN201611002595.XA
Other languages
English (en)
Other versions
CN106780641A (zh
Inventor
谢琦
孟德宇
马建华
赵谦
徐宗本
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southern Medical University
Xian Jiaotong University
Original Assignee
Southern Medical University
Xian Jiaotong University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southern Medical University, Xian Jiaotong University filed Critical Southern Medical University
Priority to CN201611002595.XA priority Critical patent/CN106780641B/zh
Publication of CN106780641A publication Critical patent/CN106780641A/zh
Application granted granted Critical
Publication of CN106780641B publication Critical patent/CN106780641B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/008Specific post-processing after tomographic reconstruction, e.g. voxelisation, metal artifact correction
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10072Tomographic images
    • G06T2207/10081Computed x-ray tomography [CT]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/421Filtered back projection [FBP]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Apparatus For Radiation Diagnosis (AREA)

Abstract

一种低剂量X射线CT图像重建方法,首先获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据;通过成像过程的统计规律构建弦图数据的统计生成模型;根据投影数据和图像的结构特征与实际应用中的需求,构建弦图数据先验的统计模型;构建完整的统计模型,并根据最大后验估计方法,将模型转化为弦图数据复原模型;应用弦图数据复原模型,得到估计的弦图数据与其它统计变量;根据得到的弦图数据进行CT图像重建,得到输出CT图,本发明旨在建立基于弦图数据生成原理及其内在统计先验的高质量弦图数据复原模型与方法,并进而实现能够大幅减少图像噪声/伪影并恢复图像细节的CT图像优质重建。

Description

一种低剂量X射线CT图像重建方法
技术领域
本发明涉及一种医学影像的图像处理技术,具体涉及一种基于统计建模与最大后验估计框架的低剂量X射线CT图像重建方法。
背景技术
X射线CT扫描已经广泛应用于临床医学影像诊断,但是CT扫描过程中过高的X射线辐射剂量对人体存在潜在风险,容易造成辐射损伤、诱发恶性肿瘤等。在保证图像质量的前提下,最大限度降低X射线使用剂量已经成为医学CT成像领域研究的关键技术之一。
为了降低X射线辐射剂量,可以通过各种硬件技术及软件技术降低CT扫描中的X射线使用剂量。常见的技术包括降低管电流、降低X射线曝光时间以及减少投影数据的采集量(即稀疏角度CT扫描)等方法。
降低CT扫描中的管电流(Low-mA)和扫描时间可以直接减少使用X线的辐射剂量,但是其相应的成像数据中,随机成分的比率将大大增加,直接导致图像质量的严重退化,难以用于临床诊断。为了在保证成像质量的前提下,最大限度降低X射线辐射剂量,基于降低管电流和扫描时间的低剂量CT图像重建方法相继被提出,例如基于统计模型的迭代重建方法以及基于投影数据滤波的解析重建方法。其中,基于统计模型的迭代重建方法根据采集投影数据的统计特性以及成像系统构建CT图像重建模型,可以实现低剂量CT图像的优质重建;基于投影数据滤波的解析重建方法同样可以根据采集投影数据的统计特性以及成像系统进行数据滤波建模,再通过解析重建方法实现快速优质的低剂量CT图像重建。
基于统计模型的迭代重建方法需要对目标函数进行几十甚至上百次的反复迭代求解,这使得重建CT图像的时间代价往往非常庞大。重建相同像素大小的CT图像时,基于统计模型的迭代重建方法所需时间往往远超经典的解析重建方法,因此仍不能满足临床CT实时显像的需求。
而传统的基于投影数据滤波的解析重建方法并未充分利用数据内在统计信息,在投影数据恢复过程中往往导致图像原有细节信息的丢失,从而使得相应CT图像的分辨率下降。
针对现有技术不足,提供一种能够确保重建图像的分辨率的低剂量CT图像重建方法因此甚为必要。
发明内容
本发明的目的在于提供一种所获重建CT图像分辨率更高,成像质量更好,细节信息保留更充分的低剂量X射线CT图像重建方法。
为达到上述目的,本发明采用的技术方案是:
步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据;
步骤S2:通过成像过程的统计规律构建生成投影数据的统计模型;
步骤S3:根据投影数据和弦图数据的结构特征与实际应用中的需求,构建数据先验的统计模型;
步骤S4:结合步骤S2与步骤S3,构建以投影数据信息为条件的弦图数据统计生成模型,并利用最大后验估计方法,构造弦图数据复原算法;
步骤S5:以步骤S1获得的投影数据为输入,应用步骤S4的弦图数据复原算法,获得复原弦图数据及其它统计变量;
若步骤S3中统计模型以重建的CT图像为其统计变量,那么步骤S5直接得到重建CT图像;
步骤S6:根据步骤S5所获的弦图数据进行CT图像重建,得到输出结果。
所述步骤S1中获取的CT设备的成像系统参数包括X射线入射光子强度I0、系统电子噪声的方差
Figure GDA0002485112760000021
所述步骤S2中,通过成像过程的统计规律构建的投影数据统计生成模型为:
p=I+ε,
Figure GDA0002485112760000038
其中,p为感受器上观测的原始投影数据,I为到达感受器的X射线光子强度,y为弦图数据,ε为系统电子噪声,它们的第i个分量分别代表第i个数据点上对应的数据;
Figure GDA0002485112760000031
代表第i个数据点上电子噪声满足的分布,通常假设为一个以σε为方差的正态分布,其形式为:
Figure GDA0002485112760000032
P{Ii|yi}代表感第i个数据点上射线光子强度满足的条件分布,其概率密度函数为以下分布形式:
Figure GDA0002485112760000033
其中,
Figure GDA0002485112760000034
表示以λ为均值的泊松分布,I0为第i个数据点上的X射线入射光子强度。
所述步骤S2中,通过成像过程的统计规律构建的投影数据统计生成模型表达为如下的条件概率形式:
Figure GDA0002485112760000035
其中,P(Ii|yi)与
Figure GDA0002485112760000036
分别由公式(1)与(2)定义。
所述步骤S3中,根据投影数据和弦图数据的结构特征与实际应用中的需求,构建数据先验的统计模型;的具体形式应根据实际情况与对运算效率与运算精度的需求确定,可归纳为如下的表达形式:
Figure GDA0002485112760000037
其中,q是必要的辅助变量,其中,N是向量q的元素总数,L(q|0,b)是均值为0,尺度参数为b的拉普拉斯分布,P(b)是关于参数b的无信息先验,值恒为常数(即足够大范围内的均匀分布),Ωq是归一化常数,其值为Ωq=∫f(y)=qydy,f(y)是根据实际需要而确定的映射。
所述步骤S4中构建的统计生成模型为如下完整后验分布形式:
Figure GDA0002485112760000041
其中,P(I,p|y)与P(y,q,b)分别由公式(3)与(4)定义。
所述步骤S4中,根据最大后验估计方法,由统计模型转化得弦图数据复原模型为如下的优化:
Figure GDA0002485112760000042
s.t.f(y)=q。
由于实际问题中f(·)往往为非线性映射,为计算方便起见,可通过变量替换q=h(z)将(6)转化为如下便于求解的等价形式:
Figure GDA0002485112760000043
s.t.g(y)=z,
其中,g(·)一般为线性映射,且h(z)满足h(g(y))=f(y)。
所述步骤S5采用交替方向乘子法求解步骤S4中的弦图数据复原模型公式(7),具体步骤包括:
S4.1)给出公式(7)的增广拉格朗日函数:
Figure GDA0002485112760000044
其中,Λ是乘子,μ是一个大于0的数;
S4.2)建立交替方向乘子法的迭代格式与终止条件:
迭代格式为:
Figure GDA0002485112760000051
Figure GDA0002485112760000052
Figure GDA0002485112760000053
Figure GDA0002485112760000054
Λk+1=Λkk(g(yk+1)-zk+1) (12)
μk+1=ρμk (13)
其中,ρ是一个大于1的正数,一般设为ρ=1.5,
迭代终止条件为:
Figure GDA0002485112760000055
S4.3)对问题(8)、(9)、(10)和(11)进行求解,给出迭代的具体算式;
S4.4)设置迭代的初始值设置为:y0=ln I0-lnp,I0=round(p),z0=g(y0),
Figure GDA0002485112760000056
Λ0=0,其中,round(·)为取整函数,p为步骤S1所获得的投影数据;
S4.5)进行(8)-(13)的迭代运算,直到满足终止条件,若模型中以重建CT图像为其统计变量,直接获取最终CT图像重建结果;否则,需执行步骤S6,即以当前的yk作为输出的弦图数据,并对该数据通过解析重建获取重建CT图像。
所述的(9)式即求解如下的问题:
Figure GDA0002485112760000057
其解Ik+1的任意分量满足
Figure GDA0002485112760000058
其中,
Figure GDA0002485112760000061
(15)式即:
Figure GDA0002485112760000062
因此,通过下面两步来求解问题(14):
S9.1)将Ik中不满足(16)式第一个不等式的分量,以步长1不停下降,直到满足(16)式第一个不等式,作为Ik+1的对应分量;
S9.2)将Ik中不满足(16)式第二个不等式的分量,以步长1不停上升,直到满足(16)式第二个不等式,作为Ik+1的对应分量。
所述(10)式即求解如下的问题:
Figure GDA0002485112760000063
根据函数h(·)的形式,一般有对应的阈值算子形式解析解:
Figure GDA0002485112760000064
其中,是对应阈值算子。
所述的(11)式即求解如下的问题:
Figure GDA0002485112760000066
其解为:
Figure GDA0002485112760000067
所述的(8)式即求解如下的问题:
Figure GDA0002485112760000068
针对不同g(·)形式,用FISTA算法,直接调用求解。
所述的步骤S6,采用滤波反投影算法,对步骤S5输出的弦图数据进行CT图像重建。
本发明相比传统方法,更加充分准确的对投影数据及其弦图数据的统计信息进行统计建模,从而所获重建CT图像分辨率更高,成像质量更好,细节信息保留更充分。
附图说明
利用附图对本发明作进一步的说明,但附图中的内容不构成对本发明的任何限制。
图1为本发明的流程图。
图2为实例1中仿真使用的假人体模数据,右下角红框中的图像为原图红框中图像提高对比度并放大四倍的结果。
图3为(a)实例1中使用的低剂量CT投影数据重建的CT图像结果及(b)利用实例1方法重建的CT图像结果,图像右下角红框中的图像为原图红框中图像提高对比度并放大四倍的结果。
图4为(a)实例1中使用的不含噪弦图数据图像;(b)低剂量情形下(含噪)弦图数据图像;(c)实例1方法估计的弦图数据图像;(d)(c)与(b)的差值绝对值图像(注所有弦图只展示左半部分,右半部分与左半部分相似)。
图5为弦图数据的近似分片平面先验展示。左上角为灰度图展示,右上角为红框处图像的放大图。可以看出,弦图数据可以较好地近似若干平面的拼合。
图6为CT图像数据的平滑先验展示。左上角为灰度图展示,右上角为红框处图像的放大图。可以看出,CT图像数据可以较好地近似若干水平平面的拼合。
具体实施方式
结合以下实例对本发明作进一步描述。
实施例1:
采用如图2所示的假人体模图像作为本发明的计算机仿真实验对象。体模图像大小设为512×512,模拟CT设备的X射线源到旋转中心和探测器的距离分别为1361.2mm和615.18mm,旋转角在[0,2π]间采样值为1160,每个采样角对应672个探测器单元,探测器单元的大小为1.85mm,弦图数据的片段如图4所示。选择弦图的分片平面近似特征(如图5所示)为先验,本发明依次包括如下步骤:
步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据p;所获取的CT设备的成像系统参数包括X射线入射光子强度I0、系统电子噪声的方差σε
步骤S2:通过成像过程的统计规律构建数据生成的统计模型;表达为如下的条件概率分布:
Figure GDA0002485112760000081
其中,p为感受器上观测的原始投影数据,I为到达感受器的X射线光子强度,ε为系统电子噪声,它们的第i个分量分别代表第i个数据点上对应的数据,P(Ii|yi)与
Figure GDA0002485112760000084
分别由(1)与(2)定义。
步骤S3:根据弦图近似分片平面的结构特征,如图5所示,可知投影数据的二阶差商场是稀疏的。故可以假设其二阶差商经过f(·)=ln(|·|+∈)-ln(∈)变换后满足拉普拉斯分布,因此,我们可以把先验分布表达为下式:
Figure GDA0002485112760000082
其中∈是一个微小的数,取值为10-16,D2是二阶差商矩阵,q是必要的辅助变量,其中,N是向量q的元素总数,L(q|0,b)是均值为0,尺度参数为b的拉普拉斯分布,P(b)是关于参数b的无信息先验,值恒为常数(即足够大范围内的均匀分布),Ωq是归一化常数,其值为Ωq=∫f(y)=qydy,f(y)是根据实际需要而确定的映射。
步骤S4:结合步骤S2与步骤S3,构建如下后验分布统计模型:
Figure GDA0002485112760000083
其中,P(I,p|y)与P(y,q|b)分别由(3)与(4)定义。
根据最大后验估计方法,可得如下弦图数据复原优化问题:
Figure GDA0002485112760000091
s.t.ln(|D2y|+∈)-ln(∈)=q,
其中,N是向量q的元素总数。以上优化问题可等价转换为:
Figure GDA0002485112760000092
s.t.D2y=z,
可以通过如下的步骤对模型进行求解。
步骤S5:基于步骤S1所输入投影数据及所设模型参数,应用步骤S4的弦图复原模型构造弦图复原求解算法。
求解的过程利用如下的迭代格式:
Figure GDA0002485112760000093
Figure GDA0002485112760000094
Figure GDA0002485112760000095
Figure GDA0002485112760000096
Λk+1=Λkk(D2yk+1-zk+1) (12)
μk+1=ρμk (13)
其中,ρ是一个大于1的正数,一般取值ρ=1.5。以下给出迭代的细节:
A.(8)式即求解如下的问题:
Figure GDA0002485112760000097
可直接调用如下FISTA算法进行求解:
算法1:FISTA算法求解子问题(50)
输入:
Figure GDA0002485112760000101
步1.设置α为I0的最大元素值的10倍,设置最大迭代步数为T=6,命a(1)=1,t=1,
Figure GDA0002485112760000102
步2.当t≤T,重复迭代以下操作:
Figure GDA0002485112760000103
Figure GDA0002485112760000104
Figure GDA0002485112760000105
d)t=t+1
步3.输出
Figure GDA0002485112760000106
B.(9)式即求解如下的问题:
Figure GDA0002485112760000107
可通过以下对问题求解:
(S9.1)将Ik中不满足(16)式第一个不等式的分量,以步长1不停下降,直到满足(35)式第一个不等式,作为
Figure GDA0002485112760000108
的对应分量。
(S9.2)将
Figure GDA0002485112760000109
中不满足(16)式第二个不等式的分量,以步长1不停上升,直到满足(16)式第二个不等式,作为Ik+1的对应分量。
C.(10)式即求解如下的问题:
Figure GDA00024851127600001010
根据函数h(·)的形式,一般有对应的阈值算子形式的解析解:
Figure GDA00024851127600001011
阈值算子
Figure GDA0002485112760000111
定义如下:
Figure GDA0002485112760000112
其中,c1=β-∈,
Figure GDA0002485112760000113
D.(11)式的解为:
Figure GDA0002485112760000114
我们设置迭代的初始值为:y0=lnI0-lnp,I0=round(p),z0=D2y0
Figure GDA0002485112760000115
Λ0=0,其中,round(·)为取整运算,p为步骤S1所获得的投影数据。迭代终止条件为:
Figure GDA0002485112760000116
步骤S6:根据步骤S5得到的复原弦图数据进行CT图像重建,得到最终输出的CT图像。
实施例2.
采用如图2所示的假人体模图像(弦图数)作为本发明的计算机仿真实验对象。体模图像大小设为512×512,模拟CT设备的X射线源到旋转中心和探测器的距离分别为1361.2mm和615.18mm,旋转角在[0,2π]间采样值为1160,每个采样角对应672个探测器单元,探测器单元的大小为1.85mm,弦图数据的片段如图4所示。选CT图像的平滑先验(如图6所示)为先验,本发明依次包括如下步骤:
步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据p;所获取的CT设备的成像系统参数包括X射线入射光子强度I0、系统电子噪声的方差σε
步骤S2:通过成像过程的统计规律构建数据生成的统计模型;表达为如下的条件概率分布:
Figure GDA0002485112760000121
其中,p为感受器上观测的原始投影数据,I为到达感受器的X射线光子强度,ε为系统电子噪声,它们的第i个分量分别代表第i个数据点上对应的数据,P(Ii|yi)与
Figure GDA0002485112760000122
分别由(1)与(2)定义。
步骤S3:根据CT图像的平滑特征(可由若干水平平面近似拟合),如图6所示,可知CT图像的一阶差商场是稀疏的,故可以假设其一阶差商经过f(·)=ln(|·|+∈)-ln(∈)变换后是满足拉普拉斯分布的。因此,我们可以把先验分布表达为下式:
Figure GDA0002485112760000123
其中∈是一个微小的数,取值为10-16,A是成像的系统矩阵,D1是一阶差商矩阵。
步骤S4:结合步骤S2与步骤S3,构建完整的统计模型,如下的后验分布形式:
Figure GDA0002485112760000124
其中,P(I,p|y)与P(y,q,b)分别由(3)与(4)定义。
根据最大后验估计方法,由统计模型转化的弦图数据复原模型可表达为如下的优化问题:
Figure GDA0002485112760000125
s.t.ln(|D1A-1y|+∈)-ln(∈)=q,
其中,N是向量q的元素总数。该优化问题可以等价转换为如下的问题:
Figure GDA0002485112760000126
s.t.D1A-1y=z。
步骤S5:基于步骤S1所输入投影数据及其所设模型参数,对步骤S4中弦图复原模型进行求解。
求解过程采用如下迭代格式:
Figure GDA0002485112760000131
Figure GDA0002485112760000132
Figure GDA0002485112760000133
Figure GDA0002485112760000134
Λk+1=Λkk(g(yk+1)-zk+1) (12)
μk+1=ρμk (13)
其中,ρ是一个大于1的正数,一般取值ρ=1.5。以下给出迭代细节:
A.(8)式即求解如下的问题:
Figure GDA0002485112760000135
可等价转化为:
Figure GDA0002485112760000136
s.t.Ax=y。
该优化问题可以用ADMM算法来求解,其增广朗日函数为:
Figure GDA0002485112760000137
该问题可通过如下迭代来求解
Figure GDA0002485112760000138
Figure GDA0002485112760000139
Γ(t+1)=Γkk(Ax(t+1)-y(t+1)) (70)
τ(t+1)=ρτ(t+1) (71)
其中,(68)即求解如下的问题:
Figure GDA0002485112760000141
该问题可直接调用如下FISTA算法进行求解:
算法2:FISTA算法求解子问题(70):
输入:
Figure GDA0002485112760000142
步1.设置α为I0的最大元素值的10倍,设置最大迭代步数为T=6,命a(1)=1,l=1,
Figure GDA0002485112760000143
步2.当l≤T,重复迭代以下操作:
a)
Figure GDA0002485112760000144
b)
Figure GDA0002485112760000145
c)
Figure GDA0002485112760000146
d)l=l+1
步3.输出
Figure GDA0002485112760000147
(69)即求解如下的问题:
Figure GDA0002485112760000148
亦可调用FISTA算法进行求解:
算法3:FISTA算法求解子问题(73)
输入:
Figure GDA0002485112760000149
步1.设置α为I0的最大元素值的10倍,设置最大迭代步数为T=6,命a(1)=1,l=1,
Figure GDA00024851127600001410
步2.当l≤T,重复迭代以下操作:
a)
Figure GDA0002485112760000151
b)
Figure GDA0002485112760000152
c)
Figure GDA0002485112760000153
d)l=l+1
e)步3.输出
Figure GDA0002485112760000154
该迭代算法以τ(0)=μk,y(0)=yk和x(0)=xk为初始,输出yk+1=y(t+1),xk+1=x(t+1)。我们一般采用的输出为:xk+1=A-1yk+1,可避免其它迭代步中A-1yk+1的运算,将其作为最终输出的重建CT图像结果。
迭代终止条件可设为:
Figure GDA0002485112760000155
经验上,设置最大迭代步数为5即可保证实验效果。
B.(9)式即求解如下的问题:
Figure GDA0002485112760000156
可以通过下面两步来求解问题(14)
(S9.1)将Ik中不满足(16)式第一个不等式的分量,以步长1不停下降,直到满足(35)式第一个不等式,作为
Figure GDA0002485112760000157
的对应分量。
(S9.2)将
Figure GDA0002485112760000158
中不满足(16)式第二个不等式的分量,以步长1不停上升,直到满足(16)式第二个不等式,作为Ik+1的对应分量。
C.(10)式即求解如下的问题:
Figure GDA0002485112760000161
根据函数h(·)的形式,一般有对应阈值算子形式的解析解:
Figure GDA0002485112760000162
阈值算子
Figure GDA0002485112760000163
定义与(18)式相同。
D.(11)式存在解析解:
Figure GDA0002485112760000164
我们设置迭代的初始值设置为:y0=lnI0-lnp,I0=round(p),z0=g(y0),
Figure GDA0002485112760000165
Λ0=0。其中,round(·)为取整运算,p是步骤S1所获得的投影数据。迭代终止条件为:
Figure GDA0002485112760000166
输出为xk+1作为最终CT图像重建结果。

Claims (10)

1.一种低剂量X射线CT图像重建方法,其特征在于包括如下步骤:
步骤S1:获取CT设备的成像系统参数和低剂量CT扫描协议下的投影数据;
步骤S2:通过成像过程的统计规律构建生成投影数据的统计模型;所述步骤S2中,通过成像过程的统计规律构建的投影数据统计生成模型为:
p=I+ε,
s.t.Ii~P(Ii|yi),
Figure FDA0002485112750000011
其中,p为感受器上观测的原始投影数据,I为到达感受器的X射线光子强度,y为弦图数据,ε为系统电子噪声,它们的第i个分量分别代表第i个数据点上对应的数据;
Figure FDA0002485112750000012
代表第i个数据点上电子噪声满足的分布,假设为一个以σε为方差的正态分布,其形式为:
Figure FDA0002485112750000013
P(Ii|yi)代表第i个数据点上射线光子强度满足的条件分布,其概率密度函数为以下分布形式:
Figure FDA0002485112750000014
其中,
Figure FDA0002485112750000015
表示以λ为均值的泊松分布,I0为第i个数据点上的X射线入射光子强度;
所述步骤S2中,通过成像过程的统计规律构建的投影数据统计生成模型表达为如下的条件概率形式:
Figure FDA0002485112750000016
其中,P(Ii|yi)与
Figure FDA0002485112750000017
分别由公式(1)与(2)定义;
步骤S3:根据投影数据和弦图数据的结构特征与实际应用中的需求,构建数据先验的统计模型;所述步骤S3中,根据投影数据和弦图数据的结构特征与实际应用中的需求,构建数据先验的统计模型的具体形式应根据实际情况与对运算效率与运算精度的需求确定,可归纳为如下的表达形式:
Figure FDA0002485112750000021
其中,q是必要的辅助变量,其中,N是向量q的元素总数,L(q|0,b)是均值为0,尺度参数为b的拉普拉斯分布,P(b)是关于参数b的无信息先验,值恒为常数(即足够大范围内的均匀分布),Ωq是归一化常数,其值为Ωq=∫f(y)=qydy,f(y)是根据实际需要而确定的映射;
步骤S4:结合步骤S2与步骤S3,构建以投影数据信息为条件的弦图数据统计生成模型,并利用最大后验估计方法,构造弦图数据复原算法;
步骤S5:以步骤S1获得的投影数据为输入,应用步骤S4的弦图数据复原算法,获得复原弦图数据及其它统计变量;
若步骤S3中统计模型以重建的CT图像为其统计变量,那么步骤S5直接得到重建CT图像;
步骤S6:根据步骤S5所获的弦图数据进行CT图像重建,得到输出结果。
2.根据权利要求1所述的低剂量X射线CT图像重建方法,其特征在于:所述步骤S1中获取的CT设备的成像系统参数包括X射线入射光子强度I0、系统电子噪声的方差
Figure FDA0002485112750000023
3.根据权利要求1至2任一项所述的低剂量X射线CT图像重建方法,其特征在于:所述步骤S4中构建的统计生成模型为如下完整后验分布形式:
Figure FDA0002485112750000022
其中,P(I,p|y)与P(y,q,b)分别由公式(3)与(4)定义。
4.根据权利要求3所述的低剂量X射线CT图像重建方法,其特征在于:
所述步骤S4中,根据最大后验估计方法,由统计模型转化得弦图数据复原模型为如下的优化:
Figure FDA0002485112750000031
s.t.f(y)=q
由于实际问题中f(·)往往为非线性映射,为计算方便起见,可通过变量替换q=h(z)将(6)转化为如下便于求解的等价形式:
Figure FDA0002485112750000032
s.t.g(y)=z,
其中,g(·)一般为线性映射,且h(z)满足h(g(y))=f(y)。
5.根据权利要求4所述的低剂量X射线CT图像重建方法,其特征在于:所述步骤S5采用交替方向乘子法求解步骤S4中的弦图数据复原模型公式(7),具体步骤包括:
S4.1)给出公式(7)的增广拉格朗日函数:
Figure FDA0002485112750000033
其中,Λ是乘子,μ是一个大于0的数;
S4.2)建立交替方向乘子法的迭代格式与终止条件:
迭代格式为:
Figure FDA0002485112750000034
Figure FDA0002485112750000035
Figure FDA0002485112750000036
Figure FDA0002485112750000037
Λk+1=Λkk(g(yk+1)-zk+1) (12)
μk+1=ρμk (13)
其中,ρ是一个大于1的正数,一般设为ρ=1.5,
迭代终止条件为:
Figure FDA0002485112750000041
S4.3)对问题(8)、(9)、(10)和(11)进行求解,给出迭代的具体算式;
S4.4)设置迭代的初始值设置为:y0=ln I0-ln p,I0=round(p),z0=g(y0),
Figure FDA0002485112750000042
Λ0=0,其中,round(·)为取整函数,p为步骤S1所获得的投影数据;
S4.5)进行(8)-(13)的迭代运算,直到满足终止条件,若模型中以重建CT图像为其统计变量,直接获取最终CT图像重建结果;否则,需执行步骤S6,即以当前的yk作为输出的弦图数据,并对该数据通过解析重建获取重建CT图像。
6.根据权利要求5所述的低剂量X射线CT图像重建方法,其特征在于:所述的(9)式即求解如下的问题:
Figure FDA0002485112750000043
其解Ik+1的任意分量满足
Figure FDA0002485112750000044
其中,
Figure FDA0002485112750000045
(15)式即:
Figure FDA0002485112750000046
因此,通过下面两步来求解问题(14)
S9.1)将Ik中不满足(16)式第一个不等式的分量,以步长1不停下降,直到满足(16)式第一个不等式,作为Ik+1的对应分量;
S9.2)将Ik中不满足(16)式第二个不等式的分量,以步长1不停上升,直到满足(16)式第二个不等式,作为Ik+1的对应分量。
7.根据权利要求5所述的低剂量X射线CT图像重建方法,其特征在于:所述(10)式即求解如下的问题:
Figure FDA0002485112750000051
根据函数h(·)的形式,一般有对应的阈值算子形式解析解:
Figure FDA0002485112750000052
其中,
Figure FDA0002485112750000053
是对应阈值算子。
8.根据权利要求5所述的低剂量X射线CT图像重建方法,其特征在于:所述的(11)式即求解如下的问题:
Figure FDA0002485112750000054
其解为:
Figure FDA0002485112750000055
9.根据权利要求5所述的低剂量X射线CT图像重建方法,其特征在于:所述的(8)式即求解如下的问题:
Figure FDA0002485112750000056
针对不同g(·)形式,用FISTA算法,直接调用求解。
10.根据权利要求1所述的低剂量X射线CT图像重建方法,其特征在于:所述的步骤S6,采用滤波反投影算法,对步骤S5输出的弦图数据进行CT图像重建。
CN201611002595.XA 2016-11-14 2016-11-14 一种低剂量x射线ct图像重建方法 Active CN106780641B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201611002595.XA CN106780641B (zh) 2016-11-14 2016-11-14 一种低剂量x射线ct图像重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201611002595.XA CN106780641B (zh) 2016-11-14 2016-11-14 一种低剂量x射线ct图像重建方法

Publications (2)

Publication Number Publication Date
CN106780641A CN106780641A (zh) 2017-05-31
CN106780641B true CN106780641B (zh) 2020-07-28

Family

ID=58968329

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201611002595.XA Active CN106780641B (zh) 2016-11-14 2016-11-14 一种低剂量x射线ct图像重建方法

Country Status (1)

Country Link
CN (1) CN106780641B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107945241B (zh) * 2017-10-25 2021-01-08 首都师范大学 一种基于边界信息扩散的x射线cl图像重建算法
WO2019128660A1 (zh) * 2017-12-29 2019-07-04 清华大学 训练神经网络的方法和设备、图像处理方法和设备以及存储介质
CN109035169B (zh) * 2018-07-19 2020-06-12 西安交通大学 一种无监督/半监督ct图像重建深度网络训练方法
CN110136218B (zh) * 2019-03-28 2022-10-28 中国人民解放军战略支援部队信息工程大学 基于噪声生成机制与数据驱动紧框架的ct投影去噪重建方法及装置
CN110084864B (zh) * 2019-04-09 2023-04-28 南京航空航天大学 一种基于能谱ct的电子密度图像重建方法
CN115115738B (zh) * 2022-08-29 2022-11-15 威海市博华医疗设备有限公司 一种对肺癌影像成像的修正方法及装置
CN116071450B (zh) * 2023-02-07 2024-02-13 深圳扬奇医芯智能科技有限公司 鲁棒的低剂量ct成像算法、装置
CN117274080B (zh) * 2023-09-13 2024-04-12 西安交通大学 一种低剂量ct弦图恢复方法及相关装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1632830A (zh) * 2003-12-22 2005-06-29 中国科学院自动化研究所 脑缺血病灶区的自动分割方法
CN102509322A (zh) * 2011-11-11 2012-06-20 刘华锋 一种基于卡尔曼滤波的pet图像重建方法
CN103413280A (zh) * 2013-08-26 2013-11-27 南方医科大学 一种低剂量x射线ct图像重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1632830A (zh) * 2003-12-22 2005-06-29 中国科学院自动化研究所 脑缺血病灶区的自动分割方法
CN102509322A (zh) * 2011-11-11 2012-06-20 刘华锋 一种基于卡尔曼滤波的pet图像重建方法
CN103413280A (zh) * 2013-08-26 2013-11-27 南方医科大学 一种低剂量x射线ct图像重建方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Nonlinear Sinogram Smoothing for Low-Dose X-Ray CT;Tianfang Li;《IEEE Transactions on nuclear science》;20041130;第2505-2513页 *
Penalizaed Weighted Alpha-Divergence Approach to Sinogram Restoration for Low-Dose X-ray Computed Tomography;Zhaoying Bian.etc.;《2012 IEEE Nuclear Science Symposium and Medical Imaging Conference Record》;20121231;第3675-3678页 *

Also Published As

Publication number Publication date
CN106780641A (zh) 2017-05-31

Similar Documents

Publication Publication Date Title
CN106780641B (zh) 一种低剂量x射线ct图像重建方法
Hu et al. Artifact correction in low‐dose dental CT imaging using Wasserstein generative adversarial networks
Heinrich et al. Residual U-net convolutional neural network architecture for low-dose CT denoising
CN110998602B (zh) 使用深度学习方法对3d牙颌面结构的分类和3d建模
Huang et al. CaGAN: A cycle-consistent generative adversarial network with attention for low-dose CT imaging
CN108898642B (zh) 一种基于卷积神经网络的稀疏角度ct成像方法
RU2705014C1 (ru) Способы и системы для обработки изображений
Kadimesetty et al. Convolutional neural network-based robust denoising of low-dose computed tomography perfusion maps
JP2020099662A (ja) X線ctシステム及び方法
CN112348936B (zh) 一种基于深度学习的低剂量锥束ct图像重建方法
Ding et al. Low-dose CT with deep learning regularization via proximal forward–backward splitting
Zhang et al. Low-dose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted least-squares
CN112770838A (zh) 使用自关注深度学习进行图像增强的系统和方法
Li et al. An efficient iterative cerebral perfusion CT reconstruction via low-rank tensor decomposition with spatial–temporal total variation regularization
CN104933744B (zh) Ct图像重建方法和系统
Bao et al. Few‐view CT reconstruction with group‐sparsity regularization
Zhang et al. Limited angle CT reconstruction by simultaneous spatial and Radon domain regularization based on TV and data-driven tight frame
Wang et al. An end-to-end deep network for reconstructing CT images directly from sparse sinograms
Yao et al. Multi-energy computed tomography reconstruction using a nonlocal spectral similarity model
Zeng et al. Full-spectrum-knowledge-aware tensor model for energy-resolved CT iterative reconstruction
CN113870138A (zh) 基于三维U-net的低剂量CT图像去噪方法及系统
Xu et al. Dictionary learning based low-dose X-ray CT reconstruction
Karimi et al. Reducing streak artifacts in computed tomography via sparse representation in coupled dictionaries
Dutta et al. Quantum denoising-based super-resolution algorithm applied to dental tomography images
Bousse et al. Systematic review on learning-based spectral CT

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant