一种不完全投影图像重建的快速迭代算法

收稿日期!!""# 9" !"&!!!!浙江大学学报"工学版#网址!’’’&()*+,-./&0(*&12*&3,!1,4 基金项目!国家"85##高技术研究发展计划资助项目$!""!::9#7"7"%& 作者简介!高欣$9%67& %’男’吉林人’博士’从事成像理论与脑神经网络研究&<$=->.(/>,44-*!H-D))&3)=&3, 第#8卷第%期 !"";年%月 浙!江!大!学!学!报"工学版# B)*+,-.)?CD1(>-,4E,>F1+/>GH"<,4>,11+>,4I3>1,31# J).K)&% I1L&!""; 一种不完全投影图像重建的快速迭代算法 高!欣9!夏顺仁!!汪元美!!罗戎蕾! $9&浙江大学V:XdVS国家重点实验室’浙江 杭州#9""!6) !&浙江大学 生物医学工程教育部重点实验室’浙江 杭州#9""!6% 摘!要!为了保证不完全投影数据的重建图像质量’同时提高重建速度’提出了一种基于共轭梯度法的快速迭代算 法&通过实时获取投影矩阵分量’以固定步长替代共轭梯度法中一维搜索最优迭代步长’在确保质量的同时缩短重 建时间&利用模拟投影数据和实际导弹断层扫描数据进行图像重建’结果表明’与卷积反投影和代数重建法相比’此 算法特别适用于扇形扫描的不完全投影数据的图像重建’在保证重建图像拟合度的同时’大大提高了重建速度& 关键词!快速迭代)图像重建)共轭梯度法)不完全投影数据)扇形扫描 中图分类号!PK%99&6#!!!!!文献标识码!:!!!!!文章编号!9""8$%6#R$!"";%"%$99"8$"; K’#$($"&’$(4"’*+%&($,-.%&(-’+"&"/%)#$&B/$(%) .&%-()/%-E*"$"E&%L"/$(%)# S:YR>,9’RU:ID*,$+1,!’b:KST*-,$=1>!’\EYN),4$.1>! $9e-&$&"L"9M$N)%$&)%9)*+O! d+?’34"5/$(26(/7"%8/&9’:$(2;4),#9""!6’+4/($) !eL"9M$N)%$&)%9)*@/)’"0/.$G1(2/(""%/(2)*P/(/8&%9)*10,.$&/)(’:$(2;4),#9""!6’+4/($% 50#$&’/$(:,)F1.?-/G>G1+-G>F1-.4)+>GD=?)+>=-4>,4+13),/G+*3G>),@-/12),3),(*4-G14+-2>1,GG13D,>[*1 ’-/L+1/1,G&PD1+1-.$G>=1)@G->,=1,G)?GD13)=L),1,G)?L+)(13G>),=-G+>A-,2GD1/*@/G>G*G>),)?3),$ /G-,G/G1L?)+)LG>=-./G1L>,GD1),1$2>=1,/>),/1-+3D>,3),(*4-G14+-2>1,G-+1-2)LG12G)1,/*+1GD11$ [*-.>GH)?+13),/G+*3G12>=-41?+)=>,3)=L.1G1L+)(13G>),/-,2/D)+G1,GD1+13),/G+*3G>),G>=1/>=*.G-,1$ )*/.H&PD11AL1+>=1,G-.+1/*.G/)?+13),/G+*3G>),?+)=3)=L*G1+$41,1+-G12L+)(13G>),/-,2+1-.=>//>.1 G)=)4+-LDHL+)(13G>),/D-F1/*331//?*..H/D)’,GD-GGD1?-/G-.4)+>GD=4*-+-,G11/GD1?>G,1//G))+>4>,-. >=-41’/L112/GD1>G1+-G>),G>=1’-,23-,@1/L13>-..H*/12?)++13),/G+*3G>),?)+>,3)=L.1G1L+)(13G>),/>, O:K/3-,,>,4& 6"78%&1#(?-/G>G1+-G>F1-.4)+>GD=)>=-41+13),/G+*3G>),)3),(*4-G14+-2>1,GG13D,>[*1)>,3)=L.1G1L+)$ (13G>),/)?-,/3-,,>,4 !!图像重建的方法主要有两种(一类是对横截剖 面的直接数学反计算’称为解析重建)另一类是一系 列的区域迭代’称为迭代重建&比如滤波反投影$?>.$ G1+@-3ZL+)(13G>),’OcQ%和卷积反投影$3),F).*$ G>),@-3ZL+)(13G>),’VcQ%属于第一类’其理论基 础是O)*+>1+切片定理*9+&第二类算法采用的是一 系列迭代技术’常被称为代数重建法$-.41@+->3+1$ 3),/G+*3G>),G13D,>[*1’:NP%&在实际应用中经常 采用不完全投影数据’这时解析重建法往往得不到 较理想的图像质量)而迭代算法可以弥补这一缺陷’ 同时可以根据不同 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 选取不同的目标函数进行迭 代重建’从而实现目标优化&本文首先介绍以最小模 估计为目标函数的 :NP算法’接着介绍最速下降 重建法和共轭梯度重建法’在此基础上引出快速迭 代算法&采用广泛应用的扇形扫描投影数据’对所有 算法进行程序模拟’并对结果进行分析比较& 9!迭代算法 9&9!成像方程 断层成像的数学模型可以表述为!设(B "Q9# Q!#$#QK%P 为被重建的向量#) B &$/5’PCK 为 P CK 维投影矩阵#其中元素$/5 为加权因子#表示 第5个像素对第/根射线线积分的贡献#亦即第/根 射线与第5个像素相交的面积&则投影重建图像问 题可用如下方程描述! ’ K 5B9 $/5Q5 B9/(!/B9#!#$#P< "9% 令*B "99#9!#$#9P%P 为已知P 个投影数据&问 题是如何由已知投影数据*求出图像向量(&注意 到大多数$/5为"#只有少数像素单元对任意给定的 射线和有贡献#因而式"9%的求解为高度病态求解& 9&:!代数重建法 代数重建法是著名的迭代重建算法之一)!*#现 在通常意义上的:NP是指a-30=-+0算法)9*&事实 上式"9%还可以写为 )+(B*+ "!% 为了解决这一问题#通常将重建问题转化为最 小化问题!=>,()+(A*(!<:NP则采用线性代 数中的逐次超松弛迭代法"/*331//>F1)F1+$+1.-A-$ G>),#IfN%来求解重建问题&此过程可以由如下方 程来描述! (/ B(/A9A( /A9+)/A9/ )/+)/ +)/< "#% 式中!)/B"$/9#$/!#$#$/K%#而)/+)/为)/与其自 身的内积& :NP迭代过程是从向量初始近似("#!( 开 始<第/]9次迭代(/]9修正时仅考虑一根射线及变 化与该射线相交的像素值#求出原投影数据与(/]9 沿该射线的射线和之差#并以该射线中加权沿此射 线的各像素再分布#使沿某条射线的像素值在不改 变其他像素的情况下与该射线对应的测量投影数据 一致#以此产生新的迭代估计(/& :NP算法的迭代形式非常简洁#计算也较快# 但是#每一次迭代的改进有限#因为:NP算法的一 次迭代只计算了一个投影角度&在实验中发现迭代次 数与质量有很大关系#迭代次数至少要在PCK,!以 上#才能重建出可以辨认的图像#并要求重建图像的 大小最好为PCK<为了便于算法之间的比较#本文 将PCK个投影方程迭代一遍称为一次迭代& 只有当迭代次数达到较高数值时代数重建法才 可以重建出质量较好的图像#否则总是存在区域灰 度不均的情况)#*#此外还存在较多的噪声#这将大大 影响图像的判断& 9&1,G#VS%来求解&共轭梯度法是最重要的共轭 方向法之一#其特点是利用上一次的搜索方向和本 次出发点的负梯度的线性组合来生成共轭方向#计 算量和存储量都较少#与最速下降法差不多);*& 当目标函数是具有正定 W1//1,矩阵的二次函 数且每次迭代所依据的步长&J 都是精确的最优步 长时#按共轭方向构造方法构造搜索方向,JD9 进行 搜索#最多可在(步之内达到极小点&但是如果(JD9 不是按精确的最优步长产生#就有可能导致,JD9 不 是下降方向&另外#若目标函数不是具有正定 W1/$ /1,矩阵的二次函数#在仅仅(次迭代后一般都不 能达到精确的极小点&这时需要令当前所得向量(( "或者((D9%为初始点#重新进行新一轮迭代&本文 采用的是O.1G3D1+$N11F1/共轭梯度法);#7*& !!快速共轭梯度法 一个重建算法的好坏应该综合评定#其中包括! 控制噪声的能力-收敛速度和能否容易地把既定的 数据包括进去&对于 :NP 算法#虽然重建时间较 短#但对于不完全投影数据重建质量不理想(而在大 体相当的时间里#采用最速下降法和共轭梯度法却 %"99第%期 高欣!等"一种不完全投影图像重建的快速迭代算法 可以在质量上有所提高&在消除图像噪声方面!共轭 梯度法能够得到最好的结果和最快的收敛&本文提 出了一种快速共轭梯度法"?-/G3),(*4-G14+-2>1,G! OVS#!它能够在保证图像质量"抑制噪声#的同时! 进一步减少重建时间& 对于迭代算法而言!最大的困难是计算$存储以 及快速检索$/5!即投影矩阵)的获取&考虑如下情 况%由98"g9!8投影数据重建一个9!8g9!8的图 像!投影矩阵的大小为PCK!其中P B98"C9!8 为投影射线的条数!KB9!8C9!8为重建像素点的 个数&系数$/5的总数量将达到9"8!这必将导致存储 和检索上的耗时!要想完整地存储这么大的矩阵是 不可能的&对于扇形射束投影模式!通常的做法是将 其重排成平行射束!利用平行射束的对称性减少投 影矩阵)的存储量&7!5’&本文提出另一种方法%每次 迭代实时计算投影矩阵)/分量!即第/条射线同重 建图像像素相交的信息"坐标和相交长度(区域面 积#&这样不仅可以避免开辟大量的内存空间!同时 也加快了迭代速度!当然速度加快的程度还受制于 获取)/ 的速度&该方法的主要思想是%对某一条射 线进行操作时!只开辟足够存放该射线与图像像素 相交信息的内存空间&对于重建图像大小为9!8g 9!8!$/5以 2)*@.1型存储的情况为例!只需开辟 9!8g9!8g8个字节存放相应信息!用完后!此段内 存再存放下一条射线与图像像素的相交信息!直到 所有的操作结束&如果再结合扇形投影的对称性!内 存节省率将达到86h&由于每根射线同像素相交的 区域是非常有限的!)/ 可以根据图形学&6’的相关知 识进行快速获取& 共轭梯度法在每次迭代中的收敛速度快于最速 下降法!意味着它的寻优能力较最速下降法好!但共 轭梯度法比最速下降法耗时长&就O.1G3D1+$N11F1/ 法而言!主要耗时包括两部分%求目标函数的梯度 +*"(J#)求解最优步长&J 进行的一维搜索& 实际应用中常常根据共轭梯度法最优步长的区 间性!用一个固定值来代替一维搜索得到的最优值! 虽然有偏差但对于收敛性没有太大影响&实验表明! 这种策略并不影响算法的最终收敛性!只是在开始 的两三步中与通过正常一维搜索得到最优步长的收 敛性略有差异"见图9!(为迭代次数!(*(为目标 函数梯度的范数#&随着迭代次数的增加!采用固定 步长得到的目标函数梯度的范数同采用一维搜索得 到的目标函数梯度的范数趋于一致!因此可以大胆 采用大多实际应用中固定步长的做法!这样可以节 省一半的时间& 图9!DJ!AM!KAM目标函数梯度的范数曲线 O>4&9!K)+=3*+F1)?)@(13G>),?*,3G>),>, IX$VS$OVS ! #!结果分析 对于算法IX$VS和OVS!本文进行了计算机 程序模拟&采用自定义ID11L$\)4-,头模型!较 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 模型右下角多了一个扁长椭圆!图像大小为9!8g 9!8!如图!所示& 图:!自定义D,""EC;%+’)头模型 O>4&!!E/1+$21?>,12ID11L$\)4-,D1-2=)21. ! 利用 98"g9!8 扇形扫描投影数据!进行 :NP$IX$VS$OVS的重建!算法结果如图#所示! 图中为迭代!"次的结果& 表9给出了重建所需的迭代次数($时间&!以 及和自定义ID11L$\)4-,头模型之间的三个误差测 度"$0$*作为重建质量好坏的三种标准<"$0$*定 义如下% !!"B’ 5 "Q5AQ5i#!(’ 5 "Q5#!! "7# 表9!重建误差及时间比较 P-@&9!V)=L-+>/),)?+13),/G+*3G>),1++)+/-,2G>=1/ 算法 ( " 0 * &(/ :NP !" "&55" "&%7! ] ;%;&!5! IX !" " ! "&57# ;&%;% ;%#&667 VS !" "&;;# "&689 "&5!8 79"&8#6 OVS !" "b "&559 #&"99 ##%&!5! "999 浙!江!大!学!学!报!工学版" 第#8卷 图4&#!N13),/G+*3G>),+1/*.G/-?G1+!"G>=1/>G1+-G>), @H2>??1+1,G-.4)+>GD=/ ! !! 0B’ 5 !Q5AQ5i"!#’ 5 !Q5AQ"!$ !5" *B ()%(]*(!!& !6" 式中&Q5为标准元素$Q5i为重建元素$Q为标准元素 的均值& 图;给出了图#中图像与自定义的ID11L$\)4$ -,图像在9BA"<9处的剖面灰度值比较&由图; 看出IX与标准图像在9_]"<9处的灰度值较VS 吻合得更好$表9中误差"的值也验证了这一点&但 表9中*是目标函数的值$对此项VS对应的值确 !!!’’原始图像!!!!!!!!! 重建图像 图=!不同方法迭代!"次后在9BA"<9处同 原图的灰度值比较 O>4&;!V)=L-+>/),)?+13),/G+*3G>),+1/*.G/-,2)+>4>,-. >=-41-G9_]"&9-?G1+!">G1+-G>),G>=1/ ! 实比IX的值小$从而验证了VS法比IX法的收敛 速度快&这主要是投影过程中存在的误差所致&从图 像上来看$VS重建结果的清晰度比IX要好&而本 文提出的OVS法$在拟合度和收敛速度两方面兼顾 了以上两种方法的优点& 为进一步验证OVS的有效性$本文采用实际导 弹断面的扫描数据$此数据采用扇束扫描方式$参数 为9!8g9!8$分别用N-=$\-Z卷积反投影和OVS 算法重建$结果如图7所示&显然OVS的重建质量 要优于N-=$\-Z卷积反投影& 图>!实际导弹断面的投影重建 O>4&7!N13),/G+*3G>),+1/*.G/?+)=+1-.=>//>.1 G)=)4+-LDHL+)(13G>),/ ! ;!结!语 从以上分析可以看出共轭梯度法求解基于最小 模估计的重建质量明显优于 :NP和最速下降法$ 而本文提出的快速共轭梯度法在保证重建图像质量 的同时大大加快了迭代速度$为迭代法求解少数据 投影重建问题提供了有效方法& 参考文献!!"."&")/"#"# (9)a:a:V$I\:K,4$b:KSS1&V),F1+41,31)?GD1/>=*.G-,1)*/ -.41@+->3+13),/G+*3G>),G13D,>[*1 !I:NP"(B)&@HHH I&’)#’/$(%)%)2"1(/’*@-’+()+$!""#$9!!8"&%76 %59& (#)W3+13),/G+*3G>), G13D,>[*1/3-,@1=-213)=L*G-G>),-..H1??>3>1,G(B)&@HHH I&’)#’/$(%)#%)2"1(/’*@-’+()+$9%%#$9!!#"&5"" 5"%& (;)\E/), \1/.1H$9%8%& (7)邓先礼&最优化技术(M)&重庆&重庆大学出版社$9%%8& (5)WEW*>&M*.G> /.>31D1.>3-.VP&I3-,-,2+13),/G+*3$ G>),(B)&2"1(/’*F,7#(/#$9%%%$!5!9"&% 9"& (6)孙家广&计算机图形学(M)&北京&清华大学出版社$ 9%%8&957 96"& 9999第%期 高欣!等"一种不完全投影图像重建的快速迭代算法