发明名称 一种基于伪极坐标TV最小化直线轨迹CT图像重建方法
摘要 本发明公开了一种基于伪极坐标TV最小化直线轨迹CT图像重建方法,克服了现有技术中,直线轨迹计算机断层成像(linear computed tomography,LCT)技术的有限角度图像重建的问题。该发明包含以下步骤——步骤1:建立TV最小化重建模型;步骤2:利用ADM最小化TV模型;步骤3:利用PPFFT实现图像空‑频域变换;步骤4:实现并运行算法,获得重建图像。该LCT重建技术基于交替方向法设计了TV最小化模型的求解算法,具有稳定的收敛性;并且,由于采用了伪极快速傅里叶变换,该算法具有优异的重建精度和计算效率。基于伪极坐标TV最小化LCT图像重建技术,在LCT技术投入实用化中具有重要意义。
申请公布号 CN104240272B 申请公布日期 2017.03.15
申请号 CN201410338497.8 申请日期 2014.07.16
申请人 中国人民解放军信息工程大学 发明人 闫镔;蔡爱龙;王林元;张瀚铭;李磊;陈健;陈建林;曾磊
分类号 G06T11/00(2006.01)I 主分类号 G06T11/00(2006.01)I
代理机构 郑州大通专利商标代理有限公司 41111 代理人 陈大通
主权项 一种基于伪极坐标TV最小化直线轨迹CT图像重建方法,其特征在于:所述方法包括以下步骤:步骤1:建立TV最小化重建模型;步骤2:利用ADM最小化TV模型;步骤3:利用PPFFT实现图像空‑频域变换;其包括:(1)伪极傅里叶变换定义:LCT重建中所应用的F<sub>p</sub>f及<img file="FDA0001207067460000011.GIF" wi="104" he="70" />定义为:<img file="FDA0001207067460000012.GIF" wi="998" he="302" /> 及<img file="FDA0001207067460000013.GIF" wi="1678" he="342" />现对F<sub>p</sub>f分析并推导快速计算方法,F<sub>p</sub>f实际上可简写为:<img file="FDA0001207067460000014.GIF" wi="1246" he="145" />其中<img file="FDA0001207067460000015.GIF" wi="166" he="111" />为方便计算且不失一般性,已令Δt=1,Δl=1;(2)对阵列f(k<sub>1</sub>,k<sub>2</sub>)的每列做2N点一维快速傅里叶变换:展开二重求和式,可以得到:<img file="FDA0001207067460000016.GIF" wi="1406" he="294" />其中,<img file="FDA0001207067460000017.GIF" wi="451" he="131" />为与求和过程无关的常量;分析内层关于变量k<sub>1</sub>的求和式:<img file="FDA0001207067460000018.GIF" wi="814" he="146" />其计算过程实际为针对阵列f(k<sub>1</sub>,k<sub>2</sub>)的每列做2N点一维DFT,其计算可使用FFT技术实现快速计算;(3)对阵列<img file="FDA0001207067460000019.GIF" wi="163" he="70" />的每行做2N点线性调频Z变换:外层对变量k<sub>2</sub>的求和:<img file="FDA00012070674600000110.GIF" wi="1605" he="190" />此式中若nα=1,则该式实际为对阵列<img file="FDA00012070674600000111.GIF" wi="158" he="71" />的每行做2N点一维DFT,此时快速计算仍可以利用FFT实现;当nα≠1时,该式实际上是对序列<img file="FDA00012070674600000112.GIF" wi="141" he="75" />做线性调频Z变换(chirp‑Z transform,CZT);令nα=η,将式简写为:<img file="FDA0001207067460000021.GIF" wi="1190" he="221" />注意到<img file="FDA0001207067460000022.GIF" wi="553" he="77" />将该关系式代入上式可得:<img file="dest_path_FDA0000538662490000033.GIF" wi="1528" he="304" />若定义<img file="FDA0001207067460000024.GIF" wi="483" he="130" />则上式可以进一步改写为:<img file="FDA0001207067460000025.GIF" wi="1323" he="225" />上式表明<img file="FDA0001207067460000026.GIF" wi="117" he="69" />可由三步操作算得:首先用序列s[k<sub>2</sub>]乘以<img file="FDA0001207067460000027.GIF" wi="147" he="68" />得到<img file="FDA0001207067460000028.GIF" wi="290" he="69" />然后利用<img file="FDA0001207067460000029.GIF" wi="260" he="70" />与s[k<sub>2</sub>]计算卷积,最后再乘以s[m]即可得到<img file="FDA00012070674600000210.GIF" wi="147" he="69" />在<img file="FDA00012070674600000211.GIF" wi="121" he="71" />的计算过程中,卷积运算占用了大部分计算时间,利用一维FFT可以替代此过程:分别先对<img file="FDA00012070674600000212.GIF" wi="260" he="72" />与s[k<sub>2</sub>]计算一维FFT,对应相乘后做一维IFFT即可;利用FFT可以使原卷积运算O(N<sup>2</sup>)复杂度降为N log N;步骤4:实现并运行算法,获得重建图像,其包括:输入数据p,初始化λ,ρ<sub>i</sub>>0,<img file="FDA00012070674600000213.GIF" wi="370" he="70" />且f<sup>(0)</sup>=f<sub>(FBP)</sub>,设定最大迭代次数N1,且令k=0,(1)数据预处理:<img file="FDA00012070674600000214.GIF" wi="815" he="147" />执行如下迭代过程:(2)更新f,公式如下:<img file="dest_path_FDA00005386624900000316.GIF" wi="1278" he="155" /> (3)更新z<sub>i</sub>,公式如下:<img file="FDA00012070674600000216.GIF" wi="1118" he="166" />(4)更新u<sub>i</sub>,公式如下:<img file="FDA00012070674600000217.GIF" wi="678" he="92" />(5)k←k+1若“未达到迭代设定的最大次数,”则返回步骤4(2),否则停止;其中,F<sub>p</sub>和<img file="FDA00012070674600000218.GIF" wi="70" he="70" />由PPFFT实现图像空‑频域变换实现快速计算;迭代达到最大迭代轮数退出时,f<sup>(N1)</sup>即为输出重建图像,其中,Δt为探元之间的长度间隔,n为光源位置索引,2N为采样位置的数量,Δl为采样步进,<img file="FDA00012070674600000219.GIF" wi="38" he="61" />为观测到的傅里叶采样数据,Di为图像i方向的差分算子,λ为数据保真项比例系数,用来调整目标函数中观测数据的一致性程度,ui为拉格朗日乘子,<img file="FDA00012070674600000220.GIF" wi="175" he="71" />为向量化的离散物体函数,组合系数矩阵<img file="FDA00012070674600000221.GIF" wi="275" he="78" />表示 对f作离散伪极傅里叶变换;点O为物体坐标系O‑xy原点,其在光源运动直线上的投影为O',为分析方便引入一个等价虚拟探测器,该探测器与x轴重合,并且其中心为光源S在x轴上的垂足O",在扫描过程中,光源S的位置用其与O'的距离l标记,探测器探元B的位置用其与O"的距离t标记,ρ<sub>i</sub>为二次项惩罚系数。
地址 450001 河南省郑州市高新区科学大道62号解放军信息工程大学