CN109118439B - 基于线积分的锥束ct去模糊方法 - Google Patents

Abstract

本发明公开了一种基于线积分的锥束CT去模糊方法。该方法首先对原始投影施行对数变换得到线积分数据；然后在线积分域建立基于系统点扩散函数的反卷积模型，包括保证数据完整性条件的数据保真项和控制噪声强度的正则化项，卷积核来自成像系统点扩散函数，并使用旋转对称的单参数高斯函数作为模型；最后使用反卷积处理后的线积分数据进行三维解析重建得到高分辨CT图像。本发明使用基于线积分的二维反卷积以解决模糊效应问题，快速而有效地提升了锥束CT图像的空间分辨率。

Description

基于线积分的锥束CT去模糊方法

技术领域

本发明属于医学成像技术领域，尤其涉及一种基于线积分的锥束CT去模糊方法。

背景技术

近年来，锥束CT由于采集速度快、图像各向同性好等特点，在微结构成像中发挥了重要作用。但是有限尺寸的X射线光源、平板探测器的分辨率和闪烁体的光子扩散效应等物理因素导致锥束CT图像存在的模糊效应，严重影响了图像空间分辨率，降低微结构的空间辨识度。

现有投影域去模糊算法难以完整而准确的模拟出整个系统的模糊效应，构建出的物理模型较为复杂，不易求解。与之相比，图像域去模糊算法数学模型简洁，并且卷积核容易由诸如实验直接测量法、模型求解的方法得到。对于三维重建，图像域去模糊方法进行三维卷积，计算复杂度较高，并且此类方法属于图像后处理方法，处理前需要先进行三维CT重建，难以直接并入到数据采集流程中。

发明内容

本发明所要解决现有技术的问题在于：1.数学模型复杂，2.对于三维数据，计算复杂度高，3.难以加入现有数据采集流程。针对现有技术的不足，本发明提出一种基于线积分的锥束CT去模糊方法，该方法本质上等价于图像域方法，而实际在线积分域施行的是二维反卷积，因此数学模型简洁且计算复杂度低。并且本发明不属于CT图像后处理方法，所施行的去模糊处理作用于三维重建前，因此能直接并入到已有图像采集流程以得到高分辨原始数据，增加数据采集的精度。

本发明的目的是通过以下技术方案来实现的：一种基于线积分的锥束CT去模糊方法，包括以下步骤：

(1)使用锥束CT成像系统获得原始投影数据p₀；

(2)对原始投影数据p₀施行对数变换得到二维线积分测量图像L′；

(3)在线积分域建立去模糊模型：引入全变分项，并对二维线积分测量图像L′进行去模糊处理，去模糊模型的目标函数如下所示：

其中L是需要求取的二维线积分修正图像，B₂代表二维卷积核矩阵，B₂与二维线积分修正图像L的乘积描述了锥束CT成像系统带来的系统模糊效应；二维线积分修正图像L的实际物理意义为沿着X射线的衰减值的累和，因此L受到非负约束；γ为正则化项因子，控制正则化项强度；||■||₂代表求取二范数，||■||TV代表求取全变分，

为置信项；

(4)求解线积分域去模糊模型的目标函数，得到二维线积分修正图像L；

(5)对所有角度的二维线积分测量图像L′均按照步骤(3)、(4)进行同样的处理，对求得的所有角度的二维线积分修正图像L进行三维解析重建，最终实现三维锥束CT图像的去模糊。

进一步地，所述步骤3中，二维卷积核矩阵B₂来自于锥束CT成像系统的点扩散函数PSF_σ，使用带单参数的二维高斯函数对所述点扩散函数进行模型化。

进一步地，所述步骤3中，全变分||L||_TV的形式如下：

||L||_TV＝∑_i，j W_i，j|L_i-L_j|＝||D₂·L||₁， (2)

其中W_i，j为权重矩阵，i和j分别表示行和列，根据像素位置将权重值设为0或1；

D₂是二维离散梯度算子

的矩阵形式，

的具体定义为：

式中，

表示矩阵

第i行第j列的元素，L_i，j表示图像L第i行第j列的元素。

进一步地，所述步骤4中，使用超松弛Chambolle-Pock方法进行目标函数的最小化，按照如下的迭代关系来对二维线积分修正图像L进行更新：

其中，k代表第k次迭代，L^k为第k次迭代的优化变量，

为第k次迭代的中间变量，γ是正则化项因子常数，p是超松弛参数，参数t和s是超松弛Chambolle-Pock方法的步长参数；T代表矩阵的转置运算；M是一个用作掩模的对角矩阵，用以解决卷积过程中的周期性边界问题，它给图像的内部像素分配权重为1，并对边界像素分配权重为0；P_γB(■)的具体形式为：

P_γB(z₂)＝min(γ，z₂) (5)

函数P₊(■)的具体形式为：

P₊(x)＝max(x，0) (6)B₂为卷积核矩阵，因此B₂·L和二维卷积

表示同一种线性变换，PSF_σ为锥束CT成像系统的点扩散函数，

表示二维卷积运算；根据卷积定理有如下关系：

B₂·L＝IFFT(FFT(PSF_σ).*FFT(I₂)) (7)

其中I₂是二维线积分修正图像L的空间域形式，“.*”表示分量相乘，FFT表示快速傅立叶变换，IFFT表示快速傅立叶逆变换；

z与

的相乘用PSF_σ和I_z表示，其中z是I_z的向量形式，如下定义：

式中Conj(■)表示求取每一个元素的共轭。

进一步地，所述步骤5中，所述三维解析重建采用标准FDK重建算法。

本发明的有益效果是：本发明在原有使用系统点扩散函数进行图像去模糊的理论下，使用基于线积分的二维反卷积实现锥束CT图像的去模糊，快速而有效地提升了锥束CT图像的空间分辨率，并且本方法易于并入已有锥束CT数据采集流程。

附图说明

图1为本发明基于线积分的锥束CT去模糊方法的整体流程图；

图2为本发明在数字模体数据上的实施结果，(a)高分辨参考图；(b)经过常规步骤重建得到的结果；(c)使用本发明的去模糊算法得到的结果；显示窗[0.0150.025]mm^-1；在(a)中使用白色虚线框标注的线对部分在图像下方进行了放大显示；

图3为本发明在物理模体数据上的实施结果，(a)高分辨参考图；(b)经过常规步骤重建得到的结果；(c)使用本发明的去模糊算法得到的结果；显示窗[0.0150.055]mm^-1；

图4为本发明在真实小鼠数据上的实施结果，(a1)-(a4)高分辨参考图；(b1)-(b4)经过常规步骤重建得到的结果；(c1)-(c4)使用本发明的去模糊算法得到的结果；(a1)中白色虚线框标注的部分在(a2)，(b2)，(c2)中被放大显示；(a3)中白色虚线框标注的部分在(a4)，(b4)，(c4)中被放大显示；第一、二行显示窗为[0.030.07]mm^-1，第三、四行显示窗为[0.0150.055]mm^-1。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

如图1所示，本发明提出的一种基于线积分的锥束CT去模糊方法，包括以下步骤：

(1)使用锥束CT成像系统获得原始投影数据p₀；

由锥束CT系统采集物理模体和小鼠数据(模拟仿真的数字模体数据不需要获得原始投影)。

(2)对原始投影数据p₀施行对数变换得到二维线积分测量图像L′；

对小鼠和物理模体数据的所有角度的投影施行对数变换得到二维线积分测量图像L′，模拟仿真的数字模体数据可通过模拟正投过程直接得到L′；

(3)在线积分域建立去模糊模型：引入全变分项，并对二维线积分图像L′进行去模糊处理；去模糊模型的目标函数如下所示：

其中L是需要求取的二维线积分修正图像，B₂代表二维卷积核矩阵，B₂与二维线积分修正图像L的乘积描述了锥束CT成像系统带来的系统模糊效应；二维线积分修正图像L的实际物理意义为沿着X射线的衰减值的累和，因此L受到非负约束；γ为正则化项因子，控制正则化项强度；||■||₂代表求取二范数，||■||_TV代表求取全变分，

为置信项；在线积分域进行去模糊处理解决的是系统模糊效应，避免了投影域模拟多个模糊因素及其求解的复杂性；本方法在线积分域建立去模糊模型，实际施行的是二维反卷积，相对于图像域去模糊方法，线积分域去模糊方法在计算效率方面具有较大优势；

(4)求解线积分域去模糊模型的目标函数，得到二维线积分修正图像L；

(4.1)全变分||L||_TV的形式如下：

||L||_TV＝∑_i，jW_i，j|L_i-L_j|＝||D₂·L||₁， (2)

其中W_i，j为权重矩阵，i和j分别表示行和列，根据像素位置将权重值设为0或1；

D₂是二维离散梯度算子

的矩阵形式，

的具体定义为：

式中，

表示矩阵

第i行第j列的元素，L_i，j表示图像L第i行第j列的元素。

(4.2)使用超松弛Chambolle-Pock方法进行目标函数的最小化，具体包括以下步骤：

(4.2.1)结合Chambolle-Pock方法和所要解决的去模糊问题，设计矩阵K以及函数F、G、F₁和F₂：

F(KL)＝F₁(B₂L)+F₂(D₂L)， (6)

其中，γ是正则化项常数；此步骤中，K矩阵的设计包含了图像模糊和图像求导两个运算过程，设计函数G(L)保证了优化变量的非负性，设计函数F₁(■)与F₂(■)满足了所解去模糊问题的目标函数中二范数和全变分项的要求；

(4.2.2)将(4.2.1)中设计的式(4)、(5)、(6)和(7)代入Chambolle-Pock方法进行推导，可得到如下求取F^*和G的近似映射prox(■)的公式：

prox_tG(L)＝P₊(L)＝max(0，L) (8)

其中，参数t和s是Chambolle-Pock方法的步长参数，本方法设置步长参数t和s满足st||K||²≤1关系以保证求解算法的快速收敛性；

(4.2.3)Chambolle-Pock算法使用如下迭代公式来对修正图像进行更新：

其中，k代表了第k次迭代，L^k为第k次迭代的优化变量，

为第k次迭代的中间变量，γ是正则化项因子常数，p是超松弛参数，T代表矩阵的转置运算。P_γB表示投影到球内区域γB＝{x| ||x||_∞≤γ}的函数，如下式：

P_γB(z₂)＝min(γ，z₂) (12)

函数P₊(■)具体形式为：

P₊(x)＝max(x，0) (13)

B₂为卷积核矩阵，因此B₂·L和二维卷积

表示同一种线性变换，PSF_σ为锥束CT成像系统的点扩散函数，

表示二维卷积运算。根据卷积定理有如下关系：

B₂·L＝IFFT(FFT(PSF_σ).*FFT(I₂)) (14)

其中I₂是二维线积分修正图像L的空间域形式，“.*”表示分量相乘，FFT表示快速傅立叶变换，IFFT表示快速傅立叶逆变换。对于锥束CT系统中的模糊效应问题，卷积过程与卷积定理所涉及到的周期性边界条件会导致图像产生“环绕”效应(图像底部的部分会模糊到图像的顶部)。本方法在式(11)中设计了一个用作掩模的对角矩阵M以解决卷积过程中的周期性边界问题，它给图像的内部像素分配权重为1，并对边界像素分配权重为0；

z与

的相乘用PSF_σ和I_z表示，其中z是I_z的向量形式，如下定义：

式中Conj(■)表示求取每一个元素的共轭。

(5)对所有角度的线积分数据都按照步骤(3)、(4)施行同样处理，对求得的线积分修正数据进行标准FDK三维重建，最终实现三维锥束CT图像的去模糊。

实施例

1.在三个实施例中，设置步长参数满足：st||K||²＝1。

2.在三个实施例中，设置超松弛参数：p＝1.7。

3.数字模体数据实施例中，正则化项参数γ＝2×10^-4；物理模体数据实施例中，正则化项参数γ＝3.5×10^-4；小鼠数据实施例中，正则化项参数γ＝2×10^-4。

4.上述技术方案中，直接用FDK算法得到的原始锥束CT图像如图2、3中的(b)和图4中的(b1)-(b4)，可以发现图中的结构边缘受到模糊效应影响，空间分辨率较低。

5.上述技术方案中，所述数字模体数据，使用本发明中的去模糊算法得到的结果，相比原始图像，空间分辨率提升21.4％，本发明的计算时间为常规三维图像域去模糊方法的4.6％。

6.上述技术方案中，所述物理模体数据，使用本发明中的去模糊算法得到的结果，相比原始图像，空间分辨率提升25.9％，本发明的计算时间为常规三维图像域去模糊方法的4.4％。

7.上述技术方案中，所述实验小鼠数据，使用本发明中的去模糊算法得到的结果，相比原始图像，空间分辨率提升12.5％，本发明的计算时间为常规三维图像域去模糊方法的4.9％。

以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种基于线积分的锥束CT去模糊方法，其特征在于，包括以下步骤：

(1)使用锥束CT成像系统获得原始投影数据p₀；

(2)对原始投影数据p₀施行对数变换得到二维线积分测量图像L′；

(3)在线积分域建立去模糊模型：引入全变分项，并对二维线积分测量图像L′进行去模糊处理，去模糊模型的目标函数如下所示：

其中L是需要求取的二维线积分修正图像，B₂代表二维卷积核矩阵，B₂与二维线积分修正图像L的乘积描述了锥束CT成像系统带来的系统模糊效应；二维线积分修正图像L的实际物理意义为沿着X射线的衰减值的累和，因此L受到非负约束；γ为正则化项因子，控制正则化项强度；‖■‖₂代表求取二范数，‖■‖_TV代表求取全变分，

为置信项；全变分‖L‖_TV的形式如下：

‖L‖_TV＝∑_i,jW_i,j|L_i-L_j|＝‖D₂·L‖₁, (2)

其中W_i,j为权重矩阵，i和j分别表示行和列，根据像素位置将权重值设为0或1；

D₂是二维离散梯度算子

的矩阵形式，

的具体定义为：

式中，

表示矩阵

第i行第j列的元素，L_i,j表示图像L第i行第j列的元素；

使用超松弛Chambolle-Pock方法进行目标函数的最小化，按照如下的迭代关系来对二维线积分修正图像L进行更新：

其中，k代表第k次迭代，L^k为第k次迭代的优化变量，

为第k次迭代的中间变量，γ是正则化项因子常数，p是超松弛参数，参数t和s是超松弛Chambolle-Pock方法的步长参数；T代表矩阵的转置运算；M是一个用作掩模的对角矩阵，用以解决卷积过程中的周期性边界问题，它给图像的内部像素分配权重为1，并对边界像素分配权重为0；P_γB(■)的具体形式为：

p_γB(z₂)＝min(γ,z₂) (5)

函数P₊(■)的具体形式为：

P₊(x)＝max(x,0) (6)

B₂为二维卷积核矩阵，因此B₂·L和二维卷积

表示同一种线性变换，PSF_σ为锥束CT成像系统的点扩散函数，

表示二维卷积运算；根据卷积定理有如下关系：

B₂·L＝IFFT(FFT(PSF_σ).*FFT(I₂)) (7)

其中I₂是二维线积分修正图像L的空间域形式，“.*”表示分量相乘，FFT表示快速傅立叶变换，IFFT表示快速傅立叶逆变换；

z与

的相乘用PSF_σ和I_z表示，其中z是I_z的向量形式，如下定义：

式中Conj(■)表示求取每一个元素的共轭；

(4)求解线积分域去模糊模型的目标函数，得到二维线积分修正图像L；

(5)对所有角度的二维线积分测量图像L′均按照步骤(3)、(4)进行同样的处理，对求得的所有角度的二维线积分修正图像L进行三维解析重建，最终实现三维锥束CT图像的去模糊。

2.根据权利要求1所述的一种基于线积分的锥束CT去模糊方法，其特征在于，步骤(3)中，二维卷积核矩阵B₂来自于锥束CT成像系统的点扩散函数PSF_σ，使用带单参数的二维高斯函数对所述点扩散函数进行模型化。

3.根据权利要求1所述的一种基于线积分的锥束CT去模糊方法，其特征在于，步骤(5)中，三维解析重建采用标准FDK重建算法。

Images