一种基于伪极坐标tv最小化直线轨迹ct图像重建方法
【专利摘要】本发明公开了一种基于伪极坐标TV最小化直线轨迹CT图像重建方法，克服了现有技术中，直线轨迹计算机断层成像（linearcomputedtomography,LCT）技术的有限角度图像重建的问题。该发明包含以下步骤——步骤1：建立TV最小化重建模型；步骤2：利用ADM最小化TV模型；步骤3：利用PPFFT实现图像空-频域变换；步骤4：实现并运行算法，获得重建图像。该LCT重建技术基于交替方向法设计了TV最小化模型的求解算法，具有稳定的收敛性；并且，由于采用了伪极快速傅里叶变换，该算法具有优异的重建精度和计算效率。基于伪极坐标TV最小化LCT图像重建技术，在LCT技术投入实用化中具有重要意义。
【专利说明】一种基于伪极坐标TV最小化直线轨迹CT图像重建方法
【技术领域】
[0001] 该发明涉及一种CT图像重建方法，特别是涉及一种基于伪极坐标TV最小化直线 轨迹CT图像重建方法。
【背景技术】
[0002] 传统CT技术一般采用旋转轨迹（如圆轨迹，螺旋轨迹等）获取不同角度下的投 影数据，需要对物体或者扫描系统做旋转操作。在一些应用场合，如工业在线检测、海关 集装箱检查、机场、车站的行李安全检查等，这些场合对扫描速率要求较高，基于旋转轨 迹的CT适用性并不佳。目前最新发展起来的一种采用大张角射线源、大面积阵列探测 器和直线扫描轨迹技术的新型X射线断层成像系统-直线轨迹CT(linear computed tomography，LCT)受到了广泛关注。在LCT扫描过程中，探测器和射线源保持相对静止，待 检测物体相对于探测器和射线源作直线运动，通过采集不同位置上的投影来重建物体的断 层图像。与传统CT成像方式相比，直线轨迹CT能够实现快速检测，在对扫描速度要求较高 的应用场合具有重要的实用价值。
[0003] 然而，受射线源张角和探测器阵列尺寸限制，LCT投影数据采集通常限定在一个有 限的范围内（小于180° )，是一个典型的有限角度问题。有限角度扫描无法提供满足经典 解析型重建算法所要求的投影数据集，因此解析重建算法难以获得高质量的重建图像。图 像重建质量是整个LCT系统投入实际应用的关键指标，因此LCT图像重建技术成为了近年 来CT成像领域研究的热点问题。基于压缩感知理论和正则化方法发展起来的重建算法，能 够对较低采样或有限角度下投影数据进行有效重建。基于压缩感知理论发展起来的TV最 小化重建算法基于图像内在的梯度稀疏特性，并利用该特性提升了有限角度重建的精度。 另一方面，LCT的投影频域采样具有伪极坐标分布的特性；该特性表明，基于伪极坐标的快 速傅里叶变换（pseudo polar fast Fourier transform, PPFFT)能够有效避免传统频域重 建中的傅里叶域的插值问题，可以获得更高的频域重建精度。因此，怎样结合TV模型和LCT 投影采样频域分布特性设计高效高精度、具有稳定收敛性的图像重建算法对于LCT技术的 发展关键问题，并且具有重要意义。
[0004] 针对直线轨迹计算机断层成像（linear computed tomography, LCT)技术的有 限角度图像重建问题，根据LCT投影频域采样分布特性，提出基于伪极坐标总变分（total variation，TV)最小化图像重建技术。该技术分为投影数据预处理，建立重建模型，设计求 解算法三部分：首先对采集到的投影数据进行频域变换得到频域观测数据，然后基于图像 梯度稀疏特性建立TV最小化重建优化模型，最后采用交替方向法和伪极坐标快速傅里叶 变换（PPFFT)设计求解算法。
[0005] 对于采样投影数据不足的重建问题，通常采用迭代类算法进行处理。传统的迭 代算法主要有基于数据域的变换类算法、代数重建技术以及统计迭代类算法。虽然传统 的迭代算法较解析类算法有较好重建质量，但是对于有限角度的重建质量通常难以满 足实用。近年来，基于压缩感知理论的重建算法分析了图像内在的稀疏特性，并利用该 特性提升了有限角度重建的精度。2007年，高河伟等人针对LCT的有限角度问题结合 GPEL Gerchberg-Papoul is-type extrapolation using Li no gram)和 TV 正则化技术 提出了 GPEL-TV算法，在抑制噪声的同时具有一定的边缘保持作用。在已有的凸集投影 (projection onto convex sets, P0CS)算法基础上，潘晓川等人以待重建图像TV最小化作 为目标提出了 ASD-POCS (adaptive-steepest-descent-projection onto convex sets)算 法，针对不完全数据CT图像重建问题取得了优于传统迭代类算法的效果。
[0006] 采用交替方向法（alternating direction method, ADM)的优化框架的重建技术 在TV最小化模型的求解中表现出了优异性能。2011年，Vandeghinste等人针对TV最小化 重建模型，将spIit-Bregman优化算法应用于稀疏角度采样的CT图像重建中，取得了较高 精度的重建。2013年，Zhang等人针对LCT有限角度问题，采用各向异性TV优化模型重建 提出了交替方向 TV最小化（alternating direction TV minimization, ADTVM)重建算法。 该算法采用交替方向法（alternating direction method, ADM)，将原优化问题的增广拉格 朗日函数拆分为两个具有解析解的子问题，求解算法较为高效稳定，对有限角度问题的解 决起到了一定的推动作用。
[0007] LCT投影采样频域分布具有特殊性的性质：所有的采样点分布于以原点为中心的 同心矩阵上，并且沿中心成对称辐射状，每条过中心的片段所在直线呈等斜率间隔分布。 LCT采样的频域分布类似于极坐标分布但又相区别，故称其为伪极坐标分布。该种频域分布 表明，采用基于伪极坐标的PPFFT技术能够实现图像空间和采样频域空间的无插值精确快 速变化。受频域重建启发，基于TV最小化重建模型，采用ADM并结合PPFFT技术是设计高 效高精度LCT重建算法的有效方法。
【发明内容】
[0008] 本发明克服了现有技术中直线轨迹计算机断层成像 (linearcomputedtomography, LCT)技术的有限角度图像重建的问题，提供一种使用效果较 好的基于伪极坐标TV最小化直线轨迹CT图像重建方法。本发明的技术解决方案是，提供 一种基于伪极坐标TV最小化直线轨迹CT图像重建方法，包括以下步骤：
[0009] 步骤1 :建立TV最小化重建模型；
[0010] 步骤2 :利用ADM最小化TV模型；
[0011] 步骤3 :利用PPFFT实现图像空-频域变换；
[0012] 步骤4 :实现并运行算法，获得重建图像。
[0013] 所述建立TV最小化重建模型包括关于变量1的傅里叶变换：即获得函数f(x，y) 的二维傅里叶变换；每个探元之间的长度间隔为At，且都能接收到光源的辐射；设光源有
【权利要求】
1. 一种基于伪极坐标TV最小化直线轨迹CT图像重建方法，其特征在于：所述方法包 括以下步骤： 步骤1:建立TV最小化重建模型； 步骤2:利用ADM最小化TV模型； 步骤3 :利用PPFFT实现图像空-频域变换； 步骤4:实现并运行算法，获得重建图像。
2. 根据权利要求1所述的基于伪极坐标TV最小化直线轨迹CT图像重建方法，其特征 在于：所述建立TV最小化重建模型包括关于变量1的傅里叶变换：即获得函数f(x，y)的二 维傅里叶变换；每个探元之间的长度间隔为At，且都能接收到光源的辐射；设光源有2N个 采样位置，其索引为n，采样步进为Λ1 ;从而〇分别离散化表达为于是 AiJyU 图像重建可以表达为寻找如下线性方程组的解： F1J ，其中，为观测到的傅里叶采样数据，/eR4A"为向量化的离散物体函数；组 合系数矩阵&eC4mvx4a"表示对f作离散伪极傅里叶变换。建立LCT的TV最小化重建模 型，也就是求解如下带约束条件的优化问题：f* =argmin| |f| |τν ； s.t.Fpf=F(p), 其中，II/LqiML，AeM4.?^为沿图像i方向的差分算子，f(p)表示对投影数 据做一维傅里叶变换。
3. 根据权利要求1所述的基于伪极坐标TV最小化直线轨迹CT图像重建方法，其特征 在于：所述的利用ADM最小化TV模型包括： (1) 建立无约束重建模型：采用ADM框架给出其求解算法；令9/ = 5，F(P) = /，则该 TV重建模型可以转化为： O.L·. J.
其中/IeR为保真项因子，用来调整目标函数中观测数据的一致性程度；式为带约束优 化模型，为了将其转化为无约束优化问题，利用增广拉格朗日函数进行如下转化：
其中标量AeR二次项惩罚系数，eC4"2为乘子； (2) 利用ADM求解模型：基于ADM框架采用分离变量的方法将转化为两个单目标优化 问题，也即
采用ADM分别对子问题么（/)和h(A)求取其最小值点，于是，基于ADM的求解迭代公 式为：
4.根据权利要求1所述的基于伪极坐标TV最小化直线轨迹CT图像重建方法，其特征 在于：所述的利用PPFFT实现图像空-频域变换包括： (1) 伪极傅里叶变换定义=LCT重建中所应用的Fnf及定义为：
现对Fpf(矽/类似）分析并推导快速计算方法，F/实际上可简写为：
其中α 7，为方便计算且不失一般性，已令Δt= 1，ΔI= 1 ; (2) 对阵列f(I^k2)的每列做2N点一维快速傅里叶变换：展开二重求和式，可以得到：
其中，C= 为与求和过程无关的常量；分析内层关于变量Ic1的求和 271 y2NJ 式：
其计算过程实际为针对阵列fGc1,k2)的每列做2N点一维DFT，其计算可使用FFT技术 实现快速计算； (3)对阵列./；[?，&]的每行做2N点线性调频Z变换：外层对变量k2的求和：
此式中若ηα= 1，则该式实际为对阵列的每行做2N点一维DFT，此时快速计算 仍可以利用FFT实现；当ηα尹1时，式实际上是对序列又"[*2]做线性调频Z变换（chirp-Z transform,CZT);令ηα=η，将式简写为：
注意到2/"Α':="Γ + A:j - (/? - Α:2),将该关系式代入式可得：
/ · \ 若定义4w] =cxp| ，则上式可以进一步改写为： \2NJ
式表明又M可由三步操作算得：首先用序列s[k2]乘以又，从]得到又"[々2].啦2]， 然后利用又》[^2]·5!^]与s[k2]计算卷积，最后再乘以s[m]即可得到又[m];在/"|>]的 计算过程中，卷积运算占用了大部分计算时间，利用一维FFT可以替代此过程：分别先对 又,,[&]·啦:]与s[k2]计算一维FFT，对应相乘后做一维IFFT即可（为了避免FFT带来的周 期循环效应，在做一维FFT之前应先对各序列补零至二倍长度）；利用FFT可以使原卷积运 算〇(N2)复杂度降为NlogN。
5.根据权利要求1所述的基于伪极坐标TV最小化直线轨迹CT图像重建方法，其特征 在于：所述的实现并运行算法，获得重建图像包括： 输入数据P，初始化λ，pi> 〇. ?,=#)，A且Zw = /；_.设定最大迭代次数 N，且令k= 0. m Sr据雨々卜神.
执行如下迭代过程： (2) 审新f，公式如下：
(3) 更新Zi,公式如下：
(4) 更新Ui,公式如下：
若"未达到迭代设定的最大次数，"则返回第2步，否则停止；其中，Fp和< 由PPFFT 实现图像空-频域变换实现快速计算；迭代达到最大迭代轮数退出时，f(N)即为输出重建图 像。
【文档编号】G06T11/00GK104240272SQ201410338497
【公开日】2014年12月24日 申请日期:2014年7月16日 优先权日:2014年7月16日
【发明者】闫镔, 蔡爱龙, 王林元, 张瀚铭, 李磊, 陈健, 陈建林, 曾磊 申请人:中国人民解放军信息工程大学