CN104574458A - 基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法 - Google Patents

Abstract

本发明公开了一种基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法，克服了现有技术中，图像重建方法仍有缺陷的问题。该发明的步骤如下：对平行束CT采集到的投影数据进行一维FFT变换，得到对应的极坐标下的投影频域数据；在基于极坐标的投影频域基础上，利用NUFFT技术实现图像空频域变换，以避免频域插值造成的精度损失，并建立TV最小化图像重建模型；在从投影频域恢复待重建图像的过程中，设计了基于TV最小化的频域优化模型，利用交替方向法对TV最小化模型通过增广Lagrangian函数法和交替方向乘子法进行迭代求解。本发明结合NUFFT技术和优化策略中先进的交替方向思想，能避免频域插值、计算和存储资源需求小，收敛性能好。

Description

基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法

技术领域

本发明涉及一种CT图像重建方法，特别是涉及一种基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法。

背景技术

计算机断层成像扫描(ComputedTomography,CT)是一种先进的透视成像技术，凭借其无损、精确、三维可视化等优点，CT已在工业无损检测和医疗诊断等各个领域获得了广泛的应用。随着CT设备性能的不断提高，其在临床应用上产生的X射线辐射问题也引起了人们的广泛关注，有效降低CT扫描的辐射剂量是目前CT领域亟待解决的热点和难点。

近年来利用稀疏角度投影进行CT重建，通过在360角度范围内按一定间隔等扫描角度采样的CT稀疏角度图像重建是一种有效的降低辐射剂量的方法，它可以明显缩短X射线扫描的照射时间。然而研究表明，由于投影数据不满足Tuy-Smith数据完备性条件，解析型重建算法无法得到良好的重建质量，重建图像往往伪影严重。目前针对稀疏角度的CT图像重建方法基本都是基于优化思路的迭代方法，然而迭代型重建算法所需的计算资源、存储资源都非常大，随着重建规模的增加，重建的存储需求和计算量迅速增加，在很多情况下，很难满足实际应用的需求。以迭代算法中通用的反投影计算为例：设待重建三维图像的各维尺寸都为N，则相应反投影的计算复杂度将高达O(N⁴)，重建一个分辨率为5123的三维图像，计算将循环687.2亿次，在普通PC机上完成如此大的计算量是非常耗时的，难以满足实际应用的要求。

中心切片定理揭示了投影空间与Fourier空间的内在联系，为直接Fourier方法(DirectFourierMethod,DFM)提供了理论基础。从中心切片定理提出开始，DFM方法就一直受到研究者关注。传统的Fourier重建方法在坐标系转化时，利用插值的方法得到笛卡尔坐标系下的频域采样点，这种方法得到的频域采样不可避免的引入了插值误差，由于时频变换的全局特性，经过Fourier反变换后误差会扩散到整个原始图像，从而造成较大的重建误差。

非标准快速Fourier变换(NonUniformFastFourierTransform，NUFFT)是在非标准空间进行傅立叶变换的一系列快速算法。该类算法首先使用窗函数对原始数据进行加权，然后计算过采样FFT，最后再使用窗函数与过采样FFT的结果进行卷积实现对非均匀坐标空间下的傅里叶变换。由于窗函数的平滑作用减小了插值误差，从而使得NUFFT算法具有较高的精确性。

经过现有技术的文献检索发现，利用非标准快速Fourier变换实现平行束CT重建的文献采取的实现方式存在对投影数据集的完整性要求高、迭代收敛慢等不足。2004年，Matej S，Fessler J和Kazantsev G在IEEE Transactions on Medical Imaging期刊上发表论文《Iterativetomographic image reconstruction using Fourier-based forward and back-projectors》，在平行束投影的迭代重建中，采用NUFFT方法来计算投影和反投影，从而加快了重建速度，该方法的不足在于迭代框架采用的是传统的CG方法，该方法收敛速度慢，且对投影数据的角度数目要求较高，针对稀疏角度的重建情形难以得到高质量的重建图像。

TV最小化方法：

近年来基于压缩感知理论和正则化方法发展起来的重建算法能够对较低采样下的投影数据进行有效重建。根据待重建目标通常具有的近似分片常数特性，可利用全变分(TotalVariation，TV)最小化的方法进行重建。

2008年，Sidky E和Pan X.在Physics in Medical and Biology期刊上发表论文《Imagereconstruction in circular cone-beam computed tomography by constrained,total-variationminimization》，提出了的基于TV最速下降和凸集投影约束相结合的ASD-POCS算法，较传统的迭代型算法，在重建质量上有较大的改善。

2011年，Vandeghinste B，Goossens B和Beenhouwer J等人在International Meeting on FullyThree-Dimensional Reconstruction in Radiology and Nuclear Medicine上发表论文《Split-Bregman-based sparse-view CT reconstruction》，用Split-Bregman算法应用于CT图像重建中，并取得了优于ASD-POCS算法的重建效果。

2013年，Zhang Han-Ming,Wang Lin-Yuan,Yan Bin等人在Chinese Physics B期刊上发表论文《Image reconstruction based on alternating direction total variation minimization in linear scancomputed tomography》，提出了一种基于TV最小化和交替方向法的ADTVM(alternatingdirection total variation minimization)迭代重建算法，在直线CT重建问题中取得了较好的重建结果。

上述基于TV最小化的重建方法均是空域迭代型重建算法，虽然迭代框架结构优于传统的迭代算法，但依旧存在计算资源和存储资源消耗大的问题。

发明内容：

本发明要解决的技术问题是：克服现有技术的缺陷，结合NUFFT技术和优化策略中先进的交替方向思想，设计一种既能避免频域插值、对计算和存储资源需求小，又具有收敛性能好、重建质量高特点的基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法。

本发明的技术方案通过如下步骤实现：

步骤1：对平行束CT采集到的投影数据进行一维FFT变换，得到对应的极坐标下的投影频域数据；

步骤2：在基于极坐标的投影频域基础上，利用NUFFT技术实现图像空频域变换，以避免频域插值造成的精度损失，并建立TV最小化图像重建模型；

步骤3：在从投影频域恢复待重建图像的过程中，设计了基于TV最小化的频域优化模型，利用交替方向法对TV最小化模型通过增广Lagrangian函数法和交替方向乘子法进行迭代求解。

步骤2的具体流程包括：在平行束CT系统中，设待重建图像为f，令θ为物体旋转的角度，s为探测器上的探元索引，探测器采集的投影为p_θ(s)。根据中心切片定理，投影数据p_θ(s)一维Fourier变换后的投影频域为P，它和极坐标分布下图像频域具有等价对应关系。令代表非标准快速Fourier变换算子，用非标准快速Fourier变换实现待重建图像f和其极坐标下图像频域之间的空频域变换，可建立如下重建模型：

结合TV正则化手段，对该问题进行求解。以TV最小化为目标函数，建立带约束的优化模型如下：

其中，||f||_TV为图像f的梯度的L1范数。

步骤3的求解流程包括：令向量w代表f的离散梯度，算子D为梯度求解运算，则变量代换后等价的约束优化模型如下：

构建相应的增广Lagrangian乘子函数：

其中，ν和λ为Lagrangian乘子，β和μ为惩罚项平衡系数。原TV最小化问题则等价于增广Lagrange函数的最小化求解问题。采用交替方向优化策略，通过分离变量将该问题分解为2个单目标优化的子问题进行交替求解，步骤如下：

(a)初始化：f⁰＝0,ν＝0,λ＝0,β＝β₀,μ＝μ₀,并令k＝0。

(b)求解w子问题：

$\underset{w}{\min }({\left|\right|w\left|\right|}_1-v^T({\mathrm{Df}}^k-w)+\frac{\mathrm{\β}}{2}\·{\left|\right|{\mathrm{Df}}^k-w\left|\right|}_2^2).$

该子问题可以利用shrinkage算子实现快速求解：

$w^{k+1}=\mathrm{shrinkage}({\mathrm{Df}}^k,\mathrm{\β},v)=\max \{{\left|\right|{\mathrm{Df}}^k-\frac{v}{\mathrm{\β}}\left|\right|}_2-\frac{1}{\mathrm{\β}},0\}\frac{{\mathrm{Df}}^k-v/\mathrm{\β}}{{\left|\right|{\mathrm{Df}}^k-v/\mathrm{\β}\left|\right|}_2}$

(c)求解f子问题：

对D_k(f)关于f求偏导并令其等于零可得f的求解表达如下：

其中M⁺表示M的Moore-Penrose伪逆，表示NUFFT算子的伴随变换算子。

(d)更新乘子ν和λ。

(e)收敛条件：

设定ε＞0，当||f^k+1-f^k||²＜ε时，则停止迭代，并输出重建图像f:＝f^k+1；否则，令k＝k+1并转回步骤(b)。

与现有技术相比，本发明基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法具有以下优点：克服现有技术的缺陷，结合NUFFT技术和优化策略中先进的交替方向思想，能避免频域插值、对计算和存储资源需求小，是一种具有收敛性能好、重建质量高的平行束CT稀疏角度采样的图像重建方法。

附图说明

图1是基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法的整体框架图；

图2是基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法中的FBP重建结果；

图3是基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法中的TV-SART重建结果；

图4是基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法的重建结果。

具体实施方式

下面结合附图和具体实施方式对本发明基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法作进一步说明：本发明的总体架构如图所示，具体步骤如下：

步骤1：对平行束CT采集到的投影数据进行一维FFT变换，得到对应的极坐标下的投影频域数据；

步骤2：在基于极坐标的投影频域基础上，利用NUFFT技术实现图像空频域变换，并建立TV最小化图像重建模型；

步骤3：利用交替方向法对TV最小化模型求解。

设f(x,y)为待重建的二维图像目标，为f对应的经过二维Fourier变换得到的图像频域，它们具有如下的关系：

$\hat{f}(u,v)=\underset{-\∞}{\overset{+\∞}{\∫}}\underset{-\∞}{\overset{+\∞}{\∫}}f(x,y)e^{-j2\mathrm{\π}(\mathrm{xu}+\mathrm{yv})}\mathrm{dxdy}.---\left(1\right)$

在平行束CT系统中，投影即为对f(x,y)的Radon变换，令θ为物体旋转的角度，s为探测器上的探元索引，则投影的表达式如下：

$p_{\mathrm{\θ}}\left(s\right)=\underset{-\∞}{\overset{+\∞}{\∫}}f(s\cos \mathrm{\θ}-t\sin \mathrm{\θ},s\sin \mathrm{\θ}+t\cos \mathrm{\θ})\mathrm{dt}.---\left(2\right)$

根据中心切片定理，投影p_θ(s)的一维Fourier变换P_θ(ρ)等于f(x,y)二维Fourier变换沿与探测器平行方向过原点的片段，因此，它们具有如下关系：

$\begin{array}{c}{P}_{\mathrm{\θ}}\left(\mathrm{\ρ}\right)=\underset{-\∞}{\overset{+\∞}{\∫}}{p}_{\mathrm{\θ}}\left(s\right)e^{-j2\mathrm{\π\ρs}}\mathrm{ds}\\ =\hat{f}{(u,v)|}_{\underset{v=\mathrm{\ρ}\sin \mathrm{\θ}}{u=\mathrm{\ρ}\cos \mathrm{\θ}}}\end{array}.---\left(3\right)$

从上式中可看到，投影数据一维Fourier变换后的投影频域，和极坐标分布下图像频域具有等价对应关系。传统的频域迭代重建算法利用标准快速Fourier变换得到待重建对象对应的笛卡尔坐标系下的图像频域数据，再通过频域插值映射到邻近的投影频域空间，以实现正投影过程；类似的，反投影过程中，由极坐标下的投影频域需先对频域空间进行笛卡尔网格插值，再变换到空域图像。这些步骤中的频域插值会引入误差，频域空间的误差会影响到空域图像的全局，因而重建的精度较低，重建图像中易出现伪影。

针对该问题，本方法采用非标准快速Fourier变换实现传统算法中的正、反投影过程，NUFFT及其伴随变换可以实现无插值的空频域间快速Fourier计算，直接得到基于极坐标系下的图像频域，由于和投影数据得到的投影频域的坐标点一一对应，避免了频域插值过程，从而大幅提高重建精度。令代表非标准快速Fourier变换算子，则有如下关系式：

根据上述关系，令P为投影频域，f为待重建图像，可建立如下重建模型：

结合目前本领域性能较佳的TV正则化手段，对该问题进行求解。以TV最小化为目标函数，建立优化模型并如下：

令变量w代表f的离散梯度，算子D为梯度求解运算，则变量代换后等价的约束优化模型如下：

构建相应的增广Lagrangian乘子函数：

其中，ν和λ为Lagrangian乘子，β和μ为惩罚项平衡系数。原TV最小化问题则等价于增广Lagrange函数的最小化求解问题。采用交替方向优化策略，通过分离变量将该问题分解为2个单目标优化的子问题进行交替求解，迭代步骤如下：

(a)初始化：f⁰＝0,ν＝0,λ＝0,β＝β₀,μ＝μ₀,并令k＝0。

(b)求解w子问题：

$\underset{w}{\min }({\left|\right|w\left|\right|}_1-v^T({\mathrm{Df}}^k-w)+\frac{\mathrm{\β}}{2}\·{\left|\right|{\mathrm{Df}}^k-w\left|\right|}_2^2).---\left(9\right)$

该子问题可以利用shrinkage算子实现快速求解：

$w^{k+1}=\mathrm{shrinkage}({\mathrm{Df}}^k,\mathrm{\β},v)=\max \{{\left|\right|{\mathrm{Df}}^k-\frac{v}{\mathrm{\β}}\left|\right|}_2-\frac{1}{\mathrm{\β}},0\}\frac{{\mathrm{Df}}^k-v/\mathrm{\β}}{{\left|\right|{\mathrm{Df}}^k-v/\mathrm{\β}\left|\right|}_2}---\left(10\right)$

(c)求解f子问题：

对D_k(f)关于f求偏导并令其等于零可得f的求解表达如下：

其中M⁺表示M的Moore-Penrose伪逆，表示NUFFT算子的伴随变换算子。

(d)更新乘子ν和λ。

(e)收敛条件：

设定ε＞0，当||f^k+1-f^k||²＜ε时，则停止迭代，并输出重建图像f:＝f^k+1；否则，令k＝k+1并转回步骤(b)。

本发明方法利用交替方向乘子法的思想，通过分离变量，将基于NUFFT的成像模型对应的TV最小化求解转换为两个子问题交替求解，由于每个子问题都具有解析解，因此本方法具有快速高效的特点。

以下用仿真实验验证本方法的的有效性，对比算法为主流的空域迭代算法TV-SART算法。实验中算法的编程和运行软件为Matlab 2012a，运行PC配置为Intel i7-3770CPU 3.40GHz和24GB内存。

用标准Shepp-Logan数字体模进行平行束CT的实验仿真，重建大小为256×256，探测器大小为512，每个维度的重建分辨率和探测器探元分辨率都为1mm，采样条件为180度范围内等角度间隔稀疏采样18个投影。实验中各算法的参数设置如下：TV-SART算法中SART的松弛因子为1.0，一轮总体迭代中，TV最速下降的迭代轮数为20；本发明提出方法的平衡参数μ为1024，β为32，NUFFT变换后每个角度下的频域采样点数为512。

FBP解析算法、TV-SART算法和本发明方法的重建结果如图2所示，其中迭代型算法的迭代次数都为400。从图中可以看出，本发明提出的方法在稀疏角度条件下能有效重建目标，保持较好的图像细节。TV-SART算法和本方法在重建精度和重建耗时的对比如表1所示。从表格的对比数据中可以看出，一方面，由于采用交替方向的优化求解思路，本方法在迭代收敛性能上优于传统的基于TV最速下降和SART相结合的TV-SART方法；另一方面，由于使用频域NUFFT算子替代耗时的空域正反投影过程，本发明提出方法在计算效率上远优于空域迭代型算法。综上，本发明提出的方法具有重建质量高、收敛速度快及计算效率高的特点，在平行束CT稀疏角度重建问题中具有很好实用性。

表1 不同方法的精度和时间对比

Claims (3)

1.一种基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法，其特征是：通过如下步骤实现：

步骤1：对平行束CT采集到的投影数据进行一维FFT变换，得到对应的极坐标下的投影频域数据；

步骤2：在基于极坐标的投影频域基础上，利用NUFFT技术实现图像空频域变换，以避免频域插值造成的精度损失，并建立TV最小化图像重建模型；

步骤3：在从投影频域恢复待重建图像的过程中，设计了基于TV最小化的频域优化模型，利用交替方向法对TV最小化模型通过增广Lagrangian函数法和交替方向乘子法进行迭代求解。

2.根据权利要求1所述的基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法，其特征是：步骤2的具体流程包括：在平行束CT系统中，设待重建图像为f，令θ为物体旋转的角度，s为探测器上的探元索引，探测器采集的投影为p_θ(s)。根据中心切片定理，投影数据p_θ(s)一维Fourier变换后的投影频域为P，它和极坐标分布下图像频域具有等价对应关系。令代表非标准快速Fourier变换算子，用非标准快速Fourier变换实现待重建图像f和其极坐标下图像频域之间的空频域变换，可建立如下重建模型：

结合TV正则化手段，对该问题进行求解。以TV最小化为目标函数，建立带约束的优化模型如下：

f^*＝arg min||f|_TV,

其中，||f||_TV为图像f的梯度的L1范数。

3.根据权利要求1所述的基于非标准快速Fourier变换和交替方向法的平行束CT稀疏角度重建方法，其特征是：步骤3的求解流程包括：令向量w代表f的离散梯度，算子D为梯度求解运算，则变量代换后等价的约束优化模型如下：

min||w||₁，Df＝w and f≥0.

构建相应的增广Lagrangian乘子函数：

其中，ν和λ为Lagrangian乘子，β和μ为惩罚项平衡系数。原TV最小化问题则等价于增广Lagrange函数的最小化求解问题。采用交替方向优化策略，通过分离变量将该问题分解为2个单目标优化的子问题进行交替求解，步骤如下：

(a)初始化：f⁰＝0,ν＝0,λ＝0,β＝β₀,μ＝μ₀,并令k＝0。

(b)求解w子问题：

$\underset{w}{\min }({\left|\right|w\left|\right|}_1-v^T({\mathrm{Df}}^k-w)+\frac{\mathrm{\β}}{2}\·{\left|\right|{\mathrm{Df}}^k-w\left|\right|}_2^2).$

该子问题可以利用shrinkage算子实现快速求解：

$w^{k+1}=\mathrm{shrinkage}({\mathrm{Df}}^k,\mathrm{\β},v)=\max \{{\left|\right|{\mathrm{Df}}^k-\frac{v}{\mathrm{\β}}\left|\right|}_2-\frac{1}{\mathrm{\β}},0\}\frac{{\mathrm{Df}}^k-v/\mathrm{\β}}{{\left|\right|{\mathrm{Df}}^k-v/\mathrm{\β}\left|\right|}_2}$

(c)求解f子问题：

对D_k(f)关于f求偏导并令其等于零可得f的求解表达如下：

其中M⁺表示M的Moore-Penrose伪逆，F_N ^T表示NUFFT算子的伴随变换算子。

(d)更新乘子ν和λ。

ν'＝ν-β(Df^k+1-w^k+1)

(e)收敛条件：

设定ε＞0，当||f^k+1-f^k||²＜ε时，则停止迭代，并输出重建图像f:＝f^k+1；否则，令k＝k+1并转回步骤(b)。