CN102663790A - 一种稀疏角度ct图像的重建方法 - Google Patents

Abstract

本发明涉及一种图像处理方法，具体涉及一种稀疏角度CT图像的重建方法，该方法由以下步骤组成：1）获取CT机的系统参数和一周期内的稀疏角度投影数据，并将所获得的稀疏角度投影数据按采集的时间顺序等分成若干组；2）以式（Ⅰ）为重建模型采用辅助函数法求解得到的式（Ⅱ）所示的迭代运算方法进行循环迭代运算，当循环次数达到预设的次数时即停止迭代运算，并将当前所得到的迭代运算结果即作为重建图像。本发明重建方法所得到的重建图像不仅组织结构明显清晰，而且整幅图的对比度显著提升。

Description

一种稀疏角度CT图像的重建方法

技术领域

本发明涉及图像数据处理，具体涉及一种稀疏角度CT图像的重建方法，重建后的CT图像可用于临床诊断和图像引导放射治疗。

背景技术

计算机断层成像技术，即CT（Computed Tomography）是获取病人断层结构信息的必不可少的成像手段。CT图像重建是指，通过探测器采集到的病人的投影数据，利用某种重建方法所得到病人的断层信息。随着平板探测器的飞速发展，锥束CT（Cone Beam CT，CBCT）广泛应用于临床诊断（如胸部成像）和肿瘤放射治疗中的图像引导放射治疗。与传统的乳房X线机照射相比，CBCT重建出的胸部体积信息，其组织检测精度大为提高。在图像引导放射治疗中，CBCT可以快速重建出病人在某个摆放位置下的高空间分辨力3D图像。

高剂量的X射线照射人体，CT可以重建出高质量的医学图像。但是，病人也会随之接收更多的X射线照射，这就增加了病人的诱发癌症和基因缺损以及其他疾病的概率，尤其是在图像引导放射治疗中。CBCT的广泛应用，更是以增大的X射线使用剂量为代价获取高质量医学图像。因此，现在医院所用的CT在保证图像质量的前提下还存在以下缺陷：1、要获取几百至近千个CT投影数据，利用滤波反投影（filtered back-projection，FBP）重建方法，才能得到满意的CT图像，因此难免对病人的照射剂量过大；2.获取近千个CT投影数据的扫描时间内，病人有可能会发生身体运动，对后续的成像质量也会产生一定的影响，因此对病人的扫描时间有待进一步缩短。

在现有技术中，解决上述问题1的一种方法是，采用降低X射线球管管电流（mAs）来达到降低照射剂量的目的。但是，mAs的降低，投影数据中的噪声含量也随之增加，这个时候利用FBP方法重建出的图像中含有大量的噪声，再经过图像域的一般滤波方法处理依然不能得到较好的图像质量，甚至会滤掉图像中的细节信息。解决上述问题2的一种方法是，在维持mAs正常的情况下，通过增大采样间隔，降低投影数据的数量，达到总体剂量的下降。但是由于投影数量大为减少，以至于用传统的FBP方法重建出来的图像含有严重的条状伪影。

公开号为CN 102314698A专利申请公开了一种“基于阿尔法散度约束的全变分最小化剂量CT重建方法”，该方法先利用射线源旋转一周所采集984个低剂量投影数据通过FBP算法得到初值图像，然后利用改进的EM（expectation maximization）算法进行保真重建，最后，进行全变分恢复，并与改进的EM算法进行交替迭代，直至满足相邻两次的全变分恢复的图像差小于某一设定的阈值。虽然上述专利申请所述的方法能够在降低mAs情况下通过采集大量低剂量投影数据重建出优质CT图像，但是仍存在以下方面的缺陷：1、虽然降低了管电流，但是需要获取近千个投影数据（984个），因此扫描时间比较长，不可避免会增加病人的运动带来对后期重建图像的质量影响；2、公开号为CN 102314698A的专利申请所述的方法不能适用于稀疏投影CT图像重建，具体表现为重建图像比较模糊。第一个缺点显而易见，容易理解，而对于此缺点，我们将具体分析为何公开号为CN 102314698A的专利申请所述的方法不适用于稀疏投影CT图像重建：

在公开号为CN 102314698A的专利申请的重建模式（降低管电流）下，由于采集了近千个低剂量投影数据，因此利用FBP重建方法得到的图像，虽然还有大量的噪声，但是图像本身并不模糊，而此时的全变分正好发挥了它的去噪优势，因此得到了高质量图像。而在采集少量正常管电流下的投影数据的情况，由于投影数较少，因此利用FBP重建出的图像本身，不仅含有大量的噪声，而且还含有严重的条状伪影和结构模糊，再利用该专利中的全变分图像恢复方法，此时的全变分图像恢复方法将承担消除噪声、伪影和模糊的艰巨任务，对于去噪声，它表现依然优秀，但是对于去模糊而言，效果却很差，因为，对含有模糊的图像的全变差去模糊，必须要有相应的卷积矩阵进行处理才能去模糊，而卷积矩阵通常是不容易得到的，因此只能“盲”去模糊，从而结果不佳。加之，少量投影下的改进的EM重建的每次迭代后的图像都是含有大量模糊成分的，因此公开号为CN 102314698A的专利申请所述的方法重建稀疏角度CT图像效果欠佳。

发明内容

本发明所要解决的技术问题是提供一种适用于稀疏角度CT图像的重建方法，该方法可显著提高重建图像质量。

本发明解决的上述问题的技术解决方案是：

一种稀疏角度CT图像的重建方法，该方法由以下步骤组成：

1）获取CT机的系统参数和一周期内的稀疏角度投影数据，并将所获得的稀疏角度投影数据按采集的时间顺序等分成若干组；

2）以下式（Ⅰ）为重建模型采用辅助函数法求解得到的下式（Ⅱ）所示的迭代运算方法按下述步骤进行重建：

D_α(y,Wx)+β||x||_TV    （Ⅰ）

式（Ⅰ）中，

式中，y表示步骤1）所述的一周期内的稀疏角度投影数据，x表示待重建图像，W表示CT系统矩阵，i表示CT系统中探测元的编号，j表示待重建图像的像素点编号，m表示每一组投影数据所对应的探测元的数量；α表示调整y和Wx之间差别的参数，α∈(-∞,+∞)； ${\left|\right|x\left|\right|}_{\mathrm{TV}}=\underset{n=1,u=1}{\overset{n_{\max },u_{\max }}{\mathrm{\Σ}}}{\▿x}_{n,u}=\underset{n=1,u=1}{\overset{n_{\max },u_{\max }}{\mathrm{\Σ}}}\sqrt{{(x_{n,u}-x_{n-1,u})}^2+{(x_{n,u}-x_{n,u-1})}^2+\mathrm{\ϵ}},$ 式中，x表示待重建图像，n、u分别表示待重建图像像素点所在的行数和列数，ε为大于0小于10^-5的常数；β为正则化系数；

$x_j^{t+1}=x_j^t*{\left\{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{W_{\mathrm{ij}}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{ij}}+\mathrm{\α\β}\▿{\left|\right|x_j\left|\right|}_{\mathrm{TV}|x_j=x_j^t}}{\left[\frac{y_i}{{\left(\mathrm{Wx}\right)}_i}\right]}^{\mathrm{\α}}\right\}}^{1/\mathrm{\α}}---\left(\mathrm{II}\right)$

式（Ⅱ）中，x、i、j、m、y、W、α和β的定义与式（Ⅰ）相同，t为自然数，表示迭代运算的次数，N表示构成待重建图像x的像素点个数，

表示全变差项在x_j处的导数；

2.1）将所述的一个周期内的稀疏角度投影数据进行滤波反投影（FilteredBackprojection，FBP）重建，得到初始化待重建图像数据[即式（Ⅱ）中x_j ^t的t=0]；

2.2）将步骤1）所得到的第一组稀疏角度投影数据和步骤2.1）所得到的初始化待重建图像数据代入下式（Ⅱ）进行迭代运算，再将所得到的迭代运算结果和第二组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，然后以同样的方法依次进行下一次迭代运算，直至完成步骤1）所得到的每一组稀疏角度投影数据的迭代运算后便进行第二个循环的迭代运算；在第二个循环的开始，则将第一个循环最后一次迭代运算结果和第一组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，直至完成步骤1）所得到的每一组稀疏角度投影数据的迭代运算后便进行第三个循环的迭代运算；依次不断地循环迭代运算，当循环次数达到预设的次数时即停止迭代运算，并将当前所得到的迭代运算结果即作为重建图像。

上述方案中，所述的调整y和Wx之间差别的参数α和正则化系数β的取值可通过有限次实验获得，本发明人通过大量实验所确定的最优值为，α=-1，β=0.001。

由于本发明所述方法重建的是稀疏角度投影数据下的CT图像，如果要像公开号为CN102314698A的专利申请方法直接采用全变差去模糊，那么必须要有相应的卷积矩阵进行处理才能去模糊，而卷积矩阵通常是不容易得到的。为了解决上述技术问题，本发明采用构造辅助函数的方法将对式（Ⅰ）所示的CT重建模型的求解转化为对所构建的辅助函数进行求解，得到如式（Ⅱ）所示的求解结果，并利用式（Ⅱ）进行迭代运算达到了进一步去模糊的目的。上述采用辅助函数的方法进行求解的过程如下所述：

设上述（Ⅰ）的CT重建模型表示为f(x)，即f(x)＝D_α(y,Wx)+β||x||_TV，在求解f(x)时，利用优化变换的思想，将问题转换成易于求解的优化问题，可以选择一个用来逼近CT重建模型极小值的辅助函数F(x,x^t)，通过不断的迭代计算来得到f(x)的极小值，但是辅助函数要满足如下两个条件：

F(x^t,x^t)＝f(x^t)，且F(x,x^t)≥f(x)；

x^t+1＝min_x F(x,x^t)时，f(x^t+1)＜f(x^t)，即最小化F(x,x^t)的过程也是在最小化f(x)。

因此，我们定义辅助函数为：

F(x,x^t)＝F_α(x,x^t)+β||x||_TV  （Ⅲ）

式（Ⅲ）中， $F_{\mathrm{\α}}(x,x^t)=\frac{1}{\mathrm{\α}(1-\mathrm{\α})}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[{\mathrm{\αy}}_i+{(1-\mathrm{\α})\left(\mathrm{Wx}\right)}_i-y_i^{\mathrm{\α}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_{\mathrm{ij}}{\left(\frac{W_{\mathrm{ij}}{x}_j}{{\mathrm{\λ}}_{\mathrm{ij}}}\right)]}^{1-\mathrm{\α}}$ 式中，y、x、W、i、j、m和α与的定义与式（Ⅰ）相同，

x^t表示待重建图像x迭代了t次后的图像；||x||_TV为全变差项，其定义也与式（Ⅰ）相同。

容易证明，上述辅助函数满足构建条件，F(x^t,x^t)＝f(x^t)，且F(x,x^t)≥f(x)。

对（Ⅲ）进行求导，并令导数为0，整理就可以得到上述（Ⅱ）式。

为了便于理解本发明所述重建方法的有益效果，以下将对式（Ⅱ）的去模糊作用进行描述。假设迭代t次后的重建图像为x^t，并且图像中存在模糊，利用公式（Ⅱ）进行迭代时，分两种情况：当

以及周围相邻像素点组成的区域为均质区域时，这个时候

点处的全变差导数

趋近为0，公式（Ⅱ）将退化为α散度重建，使得均质区域更加光滑，便实现去噪功能；当

作为模糊边缘附近点时，那么它和周围临近点的像素值存在差异，使得

点处的全变差导数

将不趋近为0，这个时候根据

和周围像素点值不同，

将会有正负值之分，这就会使得公式（Ⅱ）中与

乘积的部分较没有

这一项时变大或者变小，这种变大或者变小的趋势，通过（Ⅱ）式等式右边的计算，正是让这一点的更新值更快更准的朝着

的理想值逼近，进而实现了图像的去模糊。

本发明所述重建方法的有益效果如下所述：

由本发明所述方法重建的是稀疏角度投影数据下的CT图像，其数据量相对较少，α散度重建后的图像边缘十分模糊，因此即使采用了全变差项去模糊，其效果仍不理想。本发明首先将一周期内的稀疏角度投影数据进行分组，以加速低频量和高频量的收敛速度，改善重建图的模糊程度。更进一步地，本发明采用构造辅助函数的对CT重建模型进行求解方式，使全变差项融进了α散度重建迭代公式中，使得低频量和高频量的收敛速度进一步加快。在此基础上，再利用全变差项降低图像噪声并保护了图像的结构边缘，所得到的重建的模糊程度则显著降低。

附图说明：

图1为本发明所述重建方法的流程图；

图2为一幅采用本发明方法在稀疏角投影数据下重建出的CT断层图像；

图3为一幅采用现有技术在稀疏角投影数据下重建出的CT断层图像。

具体实施方式：

本实施例以肺癌患者胸部稀疏角度CT图像为例详细描述本发明所述重建方法的具体实施过程。

参见图1，本实施例的实施过程如下所述。

1、启动GE lightspeed 16排CT机，使CT机的球管旋转一周，依次均匀采集72张肺癌患者胸部稀疏角图片的投影数据，然后所采集的72张图片的投影数据按采集的时间顺序分成6组，每一组为12张图片的投影数据；同时记录下CT机的系统参数。

2、令发明内容中所述的式（Ⅰ）和式（Ⅱ）中的α=-1，β=0.001，并以式（Ⅰ）为重建模型采用辅助函数法求解得到的式（Ⅱ）所示的迭代运算方法按下述步骤进行重建：

2.1）将所述的一个周期内的稀疏角度投影数据进行滤波反投影重建，得到初始化待重建图像数据x_j ⁰；其中所述的滤波反投影重建为现有技术，本发明人具体实施的步骤如下所述：

首先对步骤1所述的72张图片的投影数据进行加权，然后对加权后的投影数据进行水平方向滤波，最后对水平方向滤波后的投影数据进行反投影累加。

2.2）将步骤1）所得到的第一组稀疏角度投影数据和步骤2.1）所得到的初始化待重建图像数据x_j ⁰代入下式（Ⅱ）进行迭代运算，得重建图像数据x_j ¹，再将x_j ¹和第二组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，得重建图像数据x_j ²，然后以同样的方法依次进行下一次迭代运算，直至完成6组稀疏角度投影数据的迭代运算后便进行第二个循环的迭代运算；在第二个循环的开始，则将第一个循环最后一次迭代运算结果x_j ⁶和第一组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，直至完成6组稀疏角度投影数据的迭代运算后便进行第三个循环的迭代运算；依次不断地循环迭代运算，当循环50次数时，得到的重建图像数据x_j ³⁰⁰所所形成的肺癌患者胸部断层图片如图2所示；由图2可见，重建的图像的对比度明显，且结构清晰，故停止迭代运算，并将图2作为重建图像。

为了验证本发明所述重建方法的效果，本发明人将上述72张肺癌患者胸部稀疏角图片的投影数据采用公开号为CN 102314698A的专利申请所述的重建方法进行了重建，经历迭代→全变分300个循环后的迭代结果如图3所示。将图3与图2进行比较可见，在两图的相同位置，即矩形框区域内，图2的组织结构明显清晰，而且整幅图的对比度显著提升，图像质量大为改善。

Claims (2)

1.一种稀疏角度CT图像的重建方法，该方法由以下步骤组成：

1）获取CT机的系统参数和一周期内的稀疏角度投影数据，并将所获得的稀疏角度投影数据按采集的时间顺序等分成若干组；

2）以下式（Ⅰ）为重建模型采用辅助函数法求解得到的下式（Ⅱ）所示的迭代运算方法按下述步骤进行重建：

D_α(y,Wx)+β||x||_TV    （Ⅰ）

式（Ⅰ）中，

$D_{\mathrm{\α}}(y,\mathrm{Wx})=\frac{1}{\mathrm{\α}(1-\mathrm{\α})}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}[{\mathrm{\αy}}_i+{(1-\mathrm{\α})\left(\mathrm{Wx}\right)}_i-{y_i}^{\mathrm{\α}}{{\left(\mathrm{Wx}\right)}_i}^{1-\mathrm{\α}}],$ 式中，y表示步骤1）所述的一周期内的稀疏角度投影数据，x表示待重建图像，W表示CT系统矩阵，i表示CT系统中探测元的编号，j表示待重建图像的像素点编号，m表示每一组投影数据所对应的探测元的数量；α表示调整y和Wx之间差别的参数，α∈(-∞,+∞)； ${\left|\right|x\left|\right|}_{\mathrm{TV}}=\underset{n=1,u=1}{\overset{n_{\max },u_{\max }}{\mathrm{\Σ}}}{\▿x}_{n,u}=\underset{n=1,u=1}{\overset{n_{\max },u_{\max }}{\mathrm{\Σ}}}\sqrt{{(x_{n,u}-x_{n-1,u})}^2+{(x_{n,u}-x_{n,u-1})}^2+\mathrm{\ϵ}},$ 式中，x表示待重建图像，n、u分别表示待重建图像像素点所在的行数和列数，ε为大于0小于10^-5的常数；β为正则化系数；

$x_j^{t+1}=x_j^t*{\left\{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{W_{\mathrm{ij}}}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{W}_{\mathrm{ij}}+\mathrm{\α\β}\▿{\left|\right|x_j\left|\right|}_{\mathrm{TV}|x_j=x_j^t}}{\left[\frac{y_i}{{\left(\mathrm{Wx}\right)}_i}\right]}^{\mathrm{\α}}\right\}}^{1/\mathrm{\α}}---\left(\mathrm{II}\right)$

式（Ⅱ）中，x、i、j、m、y、W、α和β的定义与式（Ⅰ）相同，t为自然数，表示迭代运算的次数，N表示构成待重建图像x的像素点个数，

表示全变差项在xj处的导数；

2.1）将所述的一个周期内的稀疏角度投影数据进行滤波反投影重建，得到初始化待重建图像数据；

2.2）将步骤1）所得到的第一组稀疏角度投影数据和步骤2.1）所得到的初始化待重建图像数据代入下式（Ⅱ）进行迭代运算，再将所得到的迭代运算结果和第二组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，然后以同样的方法依次进行下一次迭代运算，直至完成步骤1）所得到的每一组稀疏角度投影数据的迭代运算后便进行第二个循环的迭代运算；在第二个循环的开始，则将第一个循环最后一次迭代运算结果和第一组稀疏角度投影数据代入下式（Ⅱ）进行迭代运算，直至完成步骤1）所得到的每一组稀疏角度投影数据的迭代运算后便进行第三个循环的迭代运算；依次不断地循环迭代运算，当循环次数达到预设的次数时即停止迭代运算，并将当前所得到的迭代运算结果即作为重建图像。

2.根据权利要求1所述的一种稀疏角度CT图像的重建方法，其特征在于，所述的调整y和Wx之间差别的参数α=-1，所述的正则化系数β=0.001。

Images