一种基于GRENGEAT公式的锥束CT重建算法研究(1)

第13卷第9期2008年9月
中国图象图形学报
Jour na l of I m ag e and G raphi cs
V01．13，N o．9
S ep．，2008
一种基于G r engeat公式的锥束C T重建算法研究
张东平张定华张丰收
(西北工业大学现代设计与集成制造技术教育部重点实验室，西安710072)
摘要基于单圆轨迹的Fel dkam p(FD K)重建算法只有在小锥角条件下才能取得较好重建效果。

随着锥角增大，图像质量迅速下降。

因此在实际应用中受到一定限制。

在分析了物体R adon数据与单圆扫描获得的R adon数据间的关系后，结合空间可变滤波反投影重建算法(SV—FB P)，提出了一种基于G r engeat公式的叠加算法。

该重建算．法由两部分重建结果叠加而成，第1部分结果首先采用FD K算法通过对单圆扫描的投影数据进行重建来获得；然后采用外推方法获得缺失的数据。

并利用SV—FB P进行重建得到第2部分结果；最后将两部分结果进行叠加。

实验结果表明，该算法不仅有效地抑制了FD K算法重建的伪影，而且使锥角的使用范围比FD K算法提高了3—4倍。

这种新的叠加重建算法在大长物体的重建中，具有重要的理论和应用价值。

关键词单圆轨迹FD K算法空间可变滤波反投影重建算法叠加算法
中图法分类号：T P391．41文献标识码：A文章编号：1006-8961(2008)09—1649-06
A C one—beam C T R e cons t r uct i on A l gor i t hm
B ased on G r angea t’S For m ul a
Z H A N G D on g—p i ng，Z H A N G D i ng—hua，ZH A N G Feng shou
(T h e K ey Labor at ory of C ont em por ar y D e si gn and如t egr at ed M an uf act u r i ng，M i ni s tr y of Ed uc at i on，
N o ahw est e m Po l yt e c hn i c al U n i ver s i t y，X i’A n710072)
A b st r ac t FD K al gor i t h m ba sed on s i ng l e ci rcular orbi t ca n onl y be use d w hen t he con e angl e i s sm al l．T he i m age s
r econs t r uct ed by Fel dk a m p al gor i t hm bec om e bl ur r ed a nd di s tor t ed w hen t he co ne angl e i n cr ea ses．A s a r es ult t her e i s a l i m i t at ion f or pr act i cal appl i cat i on．A f t er t he r el at i onshi p be t w een t he R a d on t r ans f or m of t he object a nd t he R adon dat a acqui r ed by t he sca nni ng of t he s ingle ci r cu l ar orbi t i s anal yzed i n t he pa per a nd t he sup er po si ng al gor i t h m i s pr e sent ed ba sed on G r e ngea t’8f or m ul a acc or di ng t o t he s pace-var i a nt f il te r ed bac kpr oj ec t i on(SV—F B P)．The fi rst s ub-r es ul t i s cal cul at ed by appl yi ng t he FD K al gor i t hm t o t he pr oj ect i on dat a acqu i r ed by a s i ng l e cir cul ar orbi t s ca n．The se con d s ub—r es ult i s a cor r ect i on r ecei ved by appl y i ng SV—F B P al gor i t h m t o t he m i s si ng R adon dat a acqu i r ed by ext r a pol at i on．T he f inal r es ult i s t he s am of t he t w o sub-r esul t s．T he r es ult dem onst r a t es t he al gor i t h m i s a val i d m et hod t o r ed uce a rt if a ct s of t he
i m a ges r eco nst r uct ed by FD K al go r i t hm，and t he appl i cabl e con e angl e r ange c an be3t o4t i m e s l ar g er t han t hat of t he FD K
al g or i t hm．T hi s new sup er po si ng m e t ho d w i l l posses s i m p or t ant
t heo r et i cal a nd appl i cat i on val ue i n t he f ie l d s of bi g a nd l ong obj ect r econst r ucti on．
K eyw or ds s i n gl e cir cul ar or bi t，F D K al gor i t hm，s pace-va r i ant f il t er ed backpr oj ect i on al g or i t hm，sup er posi ng al gor i t h m
1引言
在计算机层析成像(C T)技术13趋成熟的背景下，关于3维物体重建方法的研究取得了突破性的进展，3维图像重建的理论体系也逐步确立起来。

其中，FD K算法…是一种基于圆轨迹扫描的近似重建算法，其是由Fel dkam p等人于1984年提
基金项目：国家自然科学基金项目(50375126)；航空科学基金项目(04153069)
收稿日期：2006—07一18；改回日期：2007—03一t5
第一作者简介：张东平(1980一)。

男。

现为西北工业大学硕士研究生。

主要进行工业CT无损检测与锥束CT重建算法的研究，已发表论文两篇。

E—m ai l：z dp2801395@si n a．eom
中国图象图形学报第13卷
出的。

FD K算法实际上是2维扇束滤波反投影算
法的3维扩展。

该方法具有简单、有效、快速等特
点，其在小锥角情况下(锥角小于±4。

)可取得较
好的重建效果。

许多实用的C T产品采用的就是
这种方法。

然而，由于该方法操作在不完全的投
影数据上，它本身又是一种近似的重建方法，因此
只有对源轨迹所在的平面进行的重建才是完全正
确的。

随着锥角的增大。

由于图像的质量急剧下
降，说明F D K算法的重建效果不佳，因此为了提高
图像的重建质量和减少重建误差，许多改进的方法，如G．FD K，T—FD K∞J，C C．FD K"1等相继被提出，经验证，这些方法可提高重建图像的质量，取得了较好的效果。

本文提出了一种基于G r engeat公式H o的叠加重建算法，以提高FD K算法在中大锥角扫描重建时的图像质量。

该算法先利用外推的方法获得R adon空间阴影区(缺失)数据；然后通过对阴影区的数据进行滤波反投影重建来获得阴影区数据的重建切片；最后将FD K重建结果与阴影区数据重建结果进行叠加，即得到最终的切片图像。

仿真实验表明，该算法能较好地抑制FD K算法在中大锥角扫描重建时所产生的伪影，从而提高了重建图像的质量。

图1锥束C T的3维R a d on变换
F i g．1The3D R a d on T r ansf or m of C o ne-bea m
G T
据得到物体R adon空间的各点，这样就可以利用R adon逆变换公式来精确重建物体。

但是由于锥束投影无法直接得到3维R adon空间的点，因此无法直接用R a don逆变换来重建物体。

1987年，G r angeat发现锥束投影转换为R adon 变换的导数是可行的，进而给出了锥束投影和3维R adon变换的导数的关系，并提出了基于R adon变换导数的精确重建方法。

后来，K udo提出的空间可变滤波反投影(SV—FBP)重建算法”1和A xe l s son和D a ni e l ss on提出的直接傅里叶重建算法(D FM)¨’川都是基于此公式展开研究的。

2R a don变换与锥束投影的关系3SV_FB P重建算法
R adon变换及逆变换是奥地利数学家R adon于1917年提出。

2维R adon变换就是直线积分，每条直线的积分对应于Radon空间的一个点。

3维R adon变换(见图1)是2维R adon变换的推广，其每个平面的积分对应于R adon空间的一个点。

R a don变换具体表述如下：
R f(p，／i t)=J J八z)dx(1)
其中以工)是物体各点的密度函数；P(P，露)是由P，n所确定的平面，P是原点到平面的距离，／'j是平面的法向量，R为Rodon变换。

R adon逆变换为
f(x)一专岸…，州n(2)其中，积分区域是以原点到j的距离为直径的球面。

这样，物体的重建就是先通过非平面源轨迹来获取物体各位置的锥束投影，然后根投影
当点源轨迹满足图像重建的完全条件(即如果与物体相交的每个平面都与源轨迹相交，则可精确重建该物体)时，SV—FB P重建算法则可以精确重建该物体。

单圆源轨迹不满足精确重建的完全条件，而只有非平面源轨迹，如直线+圆、垂直双圆等扫描轨迹才满足完全条件¨。

SV—FB P算法是采用与FD K算法相似的滤波反投影算法，但该算法在滤波部分与FD K算法不同，即FD K算法是采用1维斜坡滤波，而SV—F B P 算法则采用一系列加权、2维R adon变换、导数滤波、2维反投影等的2维可变滤波处理，最后通过3维加权反投影重建出断层图像。

具体步骤如下(如图2所示)：
(1)获得投影数据后，再对每幅图像乘以与FD K算法相同的加权因子；
(2)采用前向投影算法对每幅2维投影图像进行R a don 变换；
第9期
张东平等：一种基于G r engeat 公式的锥束C T 重建算法研究
(3)对R adon 空间的数据与空域中的导数滤波器沿径向方向进行卷积。

以完成导数滤波，对投影的2维R adon 变换的数据进行求导也就是对物体的3维R a don 变换进行求导(G r a ngea t 公式)；
(4)对数据冗余进行处理，乘以加权函数肘(f)，而膨(f )是由源轨迹确定的，每个R adon 空间的点f 代表一个与源轨迹相交，且与物体相交的
平面的面积分。

若用m (f )表示积分平面与源轨迹相交点的个数，则M (f )表示对1／m (f )进行平滑处理；
(5)通过2维反投影来得到滤波后的探测器数据；
(6)在轨迹的切线方向完成第2次导数滤波；(7)通过3维加权反投影得到待重建的物体。

图2空间可变滤波反投影重建算法流程图
Fi g ．2
T he s pace-var i a nt fi l t er e d
backpr o j ect i on al gor i t h m
4基于G r engeat 公式叠加重建算法
4．1
重建算法
单圆源轨迹扫描获得的投影数据不能精确重建
物体，这是由于部分与物体相交的平面无法与源轨迹相交，因此不满足完全条件。

对于这些缺失的数
数据
w
据可先采用插值与外推方法来补全，然后再进行重
建。

本文所采用的算法是将FD K 算法重建的结果与SV —F B P 算法重建的结果进行叠加，以提高重建图像的质量。

阴影区外的数据由单圆源轨迹扫描获得，可采用FD K 算法进行重建；而阴影区内的缺失数据则采用空间可变滤波反投影算法重建。

重建流程如图3所示，即首先将锥束投影数据近似假设为
一④一
图3
基于G r en geat 公式的重建算法流程图(R 代表R adon 变换，D 代表导数滤波，w 代表加权处理，
z 代表找出阴影区，2D B P 代表2维平行束反投影，3D B P 代表3维平行束反投影)
Fig ．3
T he r econs t ruct i on al gor i t h m ba sed
on
G r en geat ’8
f or m ul a (R ：R adon T r an sf or m ，D ：D e r i vat i ve fi l t eri n
g ，
W ：W ei g ht i n g ，Z ：Z e r o out s i de s hadow z one ',．2D B
P ：2D par al l el backp r oj ect i on 。

3D B P ：3D par al l el baekpr oj ect i on)甲⑨
，一
二
醪图
2一
、二
：“：一习叫i ；
i 至=
墨Ⅲl|n_．II ll ．一
；--．．：=f--．==：Ii ；-_I
．一l-．．．—E ㈠比F 拦止。

一
，一卜，一1、、
、一
一一
，_
一
一一
一一～
一、
需盟目⑨下图
中国图象图形学报第13卷
平行束投影数据，同时对平行束投影数据采用2维R adon变换来得到R adon空间的数据点，此时，R adon空间的数据是完全的；然后通过进行导数滤波、数据冗余加权处理及导数滤波运算来得到经过两次微分的R adon空间点，由于假设的平行束投影数据仅用来外推阴影区内的缺失数据，因此将阴影区外的数据全部置零后，再进行2维平行束反投影即可得到滤波后的平行束投影；接着通过对滤波后的投影数据进行3维平行束反投影计算来得到空间可变滤波反投影算法重建的结果；最后将两部分重建结果进行叠加，即得到最终的重建结果。

4．2叠加算法关键技术的研究
4．2．1数据冗余的处理
满足完全条件的非平面源轨迹的冗余度可近似地由所有与重建物体相交的平面与源轨迹相交点个数的平均数来加权处理。

平均数越大，冗余度越大。

在叠加重建算法中，由于可将锥束投影假设为平行束投影，而锥束投影则是在一个完整的圆周上扫描获得的，因此平行束投影的点源中(A)和与其对称的点源中(A)的投影完全相同，故冗余函数M。

(n)毫2。

FD K算法在冗余数据加权处理时也取M。

(n)；2(权因子为0．5)，然而由于在中心平面这种处理才是完全正确的，离中心平面越远，误差就会越大，因此也会使重建的伪影加剧。

4．2．2阴影区域的确定
如果物体以x)在一个以原点为中心的球形支撑域中，则由物体的R adon变换形成的R adon空间也是以原点为中心的同样大小的球。

采用单圆轨迹扫描方式不能够获得R adon空间的全部数据，只能获得圆环面内所包含的数据。

由于圆环面无法完全包含物体的R adon空间，因此不满足精确重建的完全条件。

这些与物体相交，且与单圆不相交的平面就构成了物体R adon空间的阴影区(图4)。

从图4中可以看出，当扫描半径较大，锥角较小时，由于阴影区域很小，因此可以将这些缺失数据忽略，直接进行重建，也可得到较好的近似效果。

当扫描半径减小，锥角增大时，由于R adon 空间缺失的数据增加，因此图像重建质量就迅速下降。

为了精确确定阴影区的大小，假设平面P用原点到平面的法矢n(0，t p)和原点到平面的距离P来表示。

以(0，t p)={si nocos‘p，s i nosi nt p，cos0}(3)
在xoy平面内单圆源轨迹的半径为r。

向过Z 轴和向量n所确定的平面进行投影，即得到如图5所示的投影面。

此时，单圆就投影成中。

与咖：之间的一段直线，而平面P则投影成一条直线。

从图5可以看出，与物体相交，且与单圆相交的平面集合为
S={(n(0，妒)，P)l0<0<1T，0≤妒<叮r，
l PI≤r。

s in0}(4)与物体相交，但与单圆不相交的平面集合为
s={(n(0，妒)，P)10<0<订，0≤妒<订，
h s i n0|<P≤r2}(5)
由于源轨迹及物体的R adon空间沿z轴具有旋转对称性，因此，当9从0到1r变化时，集合S是一个圆环体。

从物体的R adon空间中减去圆环体，就得到了阴影区S(如图4所示)。

图4R adon空间阴影区示意图
Fi g．4T he s hadow zone of t he R adon s pa ce
图5单圆源轨迹下的剖面图
Fi g．5The s ect i on of t he s i ng l e ci r cu l ar orbi
t
第9期张东平等：一种基于G r engeat公式的锥束C T重建算法研究1653
5实验与分析
为了验证本文提出的叠加重建算法的效果，通过设计3维头模型与圆盘模型对算法进行了重建试验。

3维头模型由不同材质的9个椭球体构成(如图6所示)。

扫描源轨迹的半径为171．67m m，锥角为4-200。

圆盘模型由7个相同材质的椭圆盘构成。

扫描源轨迹半径为244．58m m，锥角为4-15。

以上模型1的投影图像为256幅，投影图像大小为256×256pi xel s，像素大小为0．512m m。

实验时，采用仿真方法产生投影数据后，先采用FD K算法重建，再采用S V—FB P算法对重建图像进行校正。

从图7、图8可以看出，FD K算法在中心平面附近重建图像的质量较好，而离中心平面越远，图像伪影越大，图像质量越差。

这主要是由于锥角较大、致使阴影区缺失的数据较多而造成的。

由于叠加算法能够较好地补偿R adon空间不足的数据，因此能有效地抑制FD K算法在中大锥角扫描重建时产生的伪影。

叠加算法采用了一系列的2维可变滤波，特别是对于投影图像的2维R adon变换与2维反投影变换，由于引入了较大的计算量，从而使重建算法的效率降低(见表1)。

但随着硬件设备的快速发展，可采用单指令多数据、分布式计算等并行加速方法，相信会提高重建算法的效率。

图63维头模型
Fig．63D m odel of t he pha nt om
图7不同算法重建图像的质量比较Fig．7T he co m pa r i so n of t he i m ag
es
1654中国图象图形学报第13卷
表1 Tab．1
算法的时间和误差
T i m e cos t s of a l gor i t hm
辍
1龋
擎
督
像素
图8圆盘模型及不同重建算法得到的吸收系数曲线
Fi g．8I nt ens i t y prof il es of ph a nt om a nd s li c es of di s k m odel 图9是采用N PU—C B V C T扫描装置对航空发动机空心涡轮叶片进行扫描(扫描锥角为±15。

)，然后分别采用FD K算法和叠加算法重建结果的对比。

通过实际工程对比，由于本算法有效抑制了图像的伪影，从而提高了图像重建的质量。

图9FD K算法与叠加算法重建图像的对比Fig．9C o m p ar i so n be t w een r econs t ruct i on s li ce by FD K al gor i t h m a nd r econs t ruct i on s li ce by sup er po si ng al gor i t h m 6结论
本文根据锥束投影与物体R adon变换的导数关系，在分析了单圆扫描获得R a don空间的数据与物体R adon空间的数据关系后，提出了一种叠加重建算法。

该算法将F D K算法和空间可变滤波反投影重建算法进行了综合，并将两者的重建结果进行叠加，以抑制FD K算法在中大锥角扫描重建时所产生的伪影。

仿真实验与实际工程用例表明，该算法可有效地抑制重建伪影和提高重建图像的质量。

由于该算法无需增加其他扫描轨迹，因此扫描机构简单。

由于其在获得一幅投影图像后就可以进行滤波重建，因此扫描与重建可同时进行，这样就可有效提高扫描重建的效率。

●
参考文献(Ref er e nce s)
1Fel dk am p L A，D a vi s L C，K res s J W．Pr act i cal cone-beam a l gori t h m [J]．J our n al of t he O p t i c al So c i e t y o f A m er i ca A≈O pt i cs，I m age Sci ence，and V is i on，1984，1，612—619．
2G r a s s M．3-D cone-beam C T r ec o nst r u c t i o n f or ci rcu l a r t r aj e ct o ri e s [J]．Physi cal M edical B i o l ogy，2000，45(1)：329—347．
3Z eng K ai，C hen Zhi-q i ang，Z hang L i，e t a1．C one be am r eco n st r u ct i o n
a l gori t h m bas ed o n t w o c o nce nt r i c ci r cl es[J]．Jour nal T s i n【ghua
U ni v er si t y(Sci enc e＆T ec hnol ogy)，2004，44(6)：725—727．[曾凯，陈志强，张丽等．基于同心圆轨道的锥形束CT重建算法[J]．清华大学学报(自然科学版)，2004，44(6)：725—727．]
4G r an geat P．M a t he m a t i c a l f r am e w or k o f c o ne—be am3-D r ec o nst r u c t i o n vi a t he f i r st d eri v at i v e o f t h e r a d on t r ans f orm[A]．I n：G．T．H er m an，A．K．L oui s，F．N a t te r e r，Ed s．M at he m a t i ca l M et hods i n T om ogr aphy
[C]，New Y o r k：Spr inger-V erl ag，1991：66～97．
5K ud o H．Sai t o T．D er i vati on a n d i m p l e m e n t a t ion o f B con e beam r ec o nst r u c t i o n a l go r i t h m fo r nonplanar or bi t s[J]．I EEE Tra nsa c ti on s M edical I m ag i ng，1994，13(1)：196—211．
6A xel ss on C，D aniel s s on P E．Thr e e-di m e ns i on a l r ec o nst r u c t i o n f rom
c on e—be a m
d a t a i n0(N’l agN)t i m e[J]．Physi cal M edical B i ol ogy，
1994，39(1)：477—491．
7A x e ls s o n—J ac ob s o n C。

D e fr i se M，D aniel s s on P et a1．3D re con st ru ct i o n usi n g c o ne—be am backpr uj ect i on，t he R ad on
t ransf orm and l i n ogra m t echni ques[J]．Physi cal M edical Biol ogy，1995，40(3)：1321—1325．
8Sm i t h B D．I m a ge r ec o nst r u c t i o n f r o m c o ne—be am pr oj e ct i on：nec e ss a ry and su ffi c i en t c o nd i t i on s and r ec onst r uc ti on m et hods[J]．
I EEE Transacti ons o n M edi cal I m agi m g，1985，4(1)：14—25
．。