投影视野受限的图像重建及其误差分析

1、投影视野受限的图像重建及其误差分析赵双任1,扬新铁21多伦多DTI 成像行2西北工业大学摘 要:CT(Computer Tomography设备的探头宽度有限可引起探头视野受限 (Limited Field of View, LFOV,投影数据会因此被截断。由被截断的投影数据进行图像重建会产生截断膺像(Truncation Artifacts。外插算法(Extrapolation可以减小截断膺像。但此种算法常常会矫枉过正或矫正不足。传统的重建重投影(Reconstruction Reprojection迭代算法适合视角受限(Limited Angle of View的情况,但对投影视野受限情况

2、并没有明显的优势。最近出现了一种新的重建重投影迭代算法,该算法对投影视野受限情况很有效。本文为此新迭代算法建立理论基础,并讨论该算法的误差及合理性。关键词:膺像,截断, 图像重建,投影,视野,断层扫描,CT1.引言一些常用CT探头的宽度为40cm。探头视野大约为26cm取决于CT系统的摄像的成像比例。常用CT系统采用扇形或锥形光束。探头视野是CT系统旋转中心处探头对应的宽度。成像比例指光源到CT系统旋转中心的距离与光源到探头的距离的比值。病人的身体宽度一般大于26cm,所以探头视野小于病人的身体的宽度。此种情况称为探头视野受限。探头视野受限会引起由探头测量所得到的投影数据值在探头边缘处突然降为

3、零,此投影被称为截断投影。滤波投影(Filtered Backprojection, FPB算法1常常用于对扇形和锥形光束CT投影数据的图像重建。若将FPB算法直接用于截断的投影数据,重建的图像将包含截断膺像。截断膺像是在感兴趣区域边缘处内的亮环。此处感兴趣区域即要进行图像重建的区域,此区域一般为圆形,正好落在所有CT投影的视野之内。感兴趣区域的直径应小于CT视野的宽度,即感兴趣区域加上一个环状边缘区域所合成区域的直径正好与CT视野宽度相同。文献2采用的边缘宽度为10个像素点。但截断膺像主要集中在感兴趣区域的边缘处,增加环状区域的宽度只是选择更容易的问题去求解,而没有解决问题本身。因此本文假定

4、上述环状边缘区域的宽度为零,即感兴趣区域的直径与CT视野的宽度正好相同。CT感兴趣区域以外处重建的图像的质量很差,一般不予以考虑。消除截断膺像的常用算法为外插算法3,4。外插算法可估计投影视野外丢失了的投影数据。该算法的结果依赖于用于外插的函数。例如,常数外插,线性函数外插,多项式函数外插,COS函数外插,指数函数外插等等。外插算法重建的图形常常不是矫枉过正就是矫正不足。如果矫枉过正,感兴趣区域边缘处的亮环变成暗环。如果矫正不足此亮环有所减弱,但仍不消失。小波(Wavelet算法也常被用来减轻截断膺像。但是小波算法常常接合外插算法一起使用2,3。文献3承认在所谓小波算法中抑制截断膺像的实际上还

5、是外插算法。已知病人某些先决条件的方法也被用来解决探头视野受限的问题5。如果病人的身体落在感兴趣区域外只是一小部分,比如两个胳膊,这时Hibert变换算法可给出解析解6。局部断层扫描(LocalTomography算法7,8可完全消除截断膺像。但该算法在消除截断膺像的同时使图像凹陷而成杯状。图像重建重投影的迭代算法10,11给出了对视角受限(Limited angle ofview问题的计算结果并提到该算法对探头视野受限问题也有效,但没有给出这方面的计算结果。至今也没有看到利用此算法得到比外插算法更好的计算结果。文献11将该算法用于解决偏移探头(Offset Detector的投影数据截断问题

6、。此种情况下,对应一个扫描角度的投影,探头有两个位置。探头的视野由两个位置叠加来定,因此其视野一般并不受限。此种情况问题是要消除两个位置投影数据的重叠,而非视野受限。文献12推出了一个新的重建重投影迭代算法,该算法对解决投影视野受限的问题已取得比外插方法更好的结果。但是因为没有一个较为完善的理论基础,此算法很容易湮没在众多解决投影视野受限的算法之中,而不被人们重视。文献13试图为此算法建立一个理论基础。本文修改文献13的一些笔误,并增加了对该算法的误差和合理性的详细讨论。进一步加深对此迭代算法的理解。2.要解决的问题定义图像X的投影为。图像p (x X X =,其所在的区域为,即。为一个像素点

7、。图像x x X 称为原像。投影,(u p p =。是投影角,为一个投影点。投影所在的区域为u ,u 。记投影算子为P ,图像重建算子或解算子为,所以,R r X 为X所重建的图像,他们之间有关系为Eq 1XP p =Eq 2p R X r=以平行光束为例,其投影算子P 为Radon 变换。其对应的解算子可解析得给出,即FBP 算法。R 定义,图像所在的区域被分成两个区域,一个是感兴趣的区域,记为,其上图像记为;另一个是感兴趣区域之外的区域,记为,iiX o,其上的图像记为。因此有o X o i +=,。 =o i X X X定义,是投影视野之内的区域。t h是投影视野之外的区域。上的投影为,

8、可称为截断的投影。区域ttp tp h 对应的投影为,可称为空心投影。投影可写为,。hp hp p =h t p p p 式(Eq 1, Eq 2关系可改写为,Eq 3=o i ho hi to ti h t X X P P P P p p 和Eq 4=h t oh ot ih it r o r i p p R R R R X X 投影数据视野受限,投影数据被截断的图像重建问题是可归结为:假设我们已知投影矩阵,解矩阵和截断投影,但空心投影未知。在此情况下,求解。ho hito tiP P P P oh otih it R R R R tp hp riX3.断层扫描系统的特性(1 X在区域内,对

9、于平行光束和扇形光束断层扫描系统的图像重建rX应该与它的原像X 相同(至少在理论上如此,即Eq 5X X P R =上面公式在任意区域内成立,因此对区域i 也有效,即 若图像在区域为内,有 i X i Eq 6ii ti it ri X X P R X =“”为等号,考虑的任意性上式可表示为 i X 引理 1 Eq 7iti it I P R =式中是在区域上的单位算子。上式的一个特殊情况为,该式在上成立。iI iIP R = (2 因为X 光射线是直线。对区域iX h上没有贡献。这可表示为,Eq 8=i hi X P考虑到的任意性有iX 引理 2 Eq 9=hi P上式中为零算子。0(3 另

10、一方面为内区iti X P i上图像的投影,由此投影进行图像重建得到的图像不应对外区有贡献,oEq 100=i ti ot X P R或得 引理 3 Eq 11=ti ot P R上式中为零算子。0 (4 如果已知,在投影区域o X oX t上的贡献为。从截断投影数据中排除这部分贡献可得o to X P t p o to t X P p 。此投影数据只含内区域的贡献。因此是非截断投影。由它可对内区图像进行完全的重建。即,i o to t X P p 引理 4Eq 12(o to t it i ri X P p R X X =(5假定是图像'X 是对图像X作了一些改变的图像。但是该改变只

11、发生在区域上,所以有oEq 13 i i X X ='Eq 14+=X X X o o ' 此处是对原图像的改变量。记+X Eq 15 +X P p to t Eq 16+X P p ho h式(Eq 15,Eq 16是在视野内区域+X t 和视野之外区域对投影的贡献。它们都不为零。h 'X 的投影可表示为Eq 17+=t t t p p p ' Eq 18+=h h h p p p '由此也可进行图像重建Eq 19=''''h t oh otih it r o r i p p R R R R X X 在此情况下应有 Eq

12、 20ii r i X X X =''上式中第一个等号是因为在区域i 内重建的图像应与原图像相同(Eq 6。第二个等式是因为假定原图像在区域i内没有发生改变(Eq 13。(Eq 20可总结为引理 5假定是图像'X 对图像X 作了改变的图像。但是改变只发生在区域o上。由(Eq 19, Eq 20有,Eq 21=''h t ih iti p p R R X 其中由式(Eq 17,Eq 18给出。 ''h t p p4. 投影视野受限情况下的图像重建 引理 5告诉我们如果图像的改变只发生在外区域+X o 上,由式(Eq 21有Eq 22'

13、;'h ih t it i p R p R X +=引理 5只要求图像在区域+X o 上,不对作其它限制。因此可取满足+X +X Eq 230'=+X X X o o式(Eq 23是本文的关键。应注意的是并非应有0'=o X ,而是可以通过对的选择实现。式(Eq 3可写成,+X 0'=o X Eq 24'''o ho i hi h X P X P p +=式(Eq 9告诉我们;式(Eq 23告诉我们=hi P 0'=o X ,因此有Eq 25'=h p蒋上式代入式(Eq 22得, Eq 26't it i p R

14、X =但是怎样才可以使?式(Eq 19 上半部分可写为, 0'=o XEq 27'''h oh t ot o p R p R X += 考虑到式(Eq 23和式(Eq 250'=o X 0'=h p 有Eq 28'=t ot p R考虑到式(Eq 17得, Eq 290(=+t t ot p p R或Eq 30 tot t ot p R p R =+上面方程的解一般不是唯一的。其中一个解可写为,Eq 31t ot ot t p R R p +=上标 “+”代表广义逆,是对的近似值。考虑式(Eq 26和式(Eq17+t p+t p Eq 3

15、2(+=t t it i p p R X将式(Eq 31 代入上式中的得+t p+t p Eq 33(t ot ot t it i p R R p R X +=i X 是的计算值,它有可能不是精确的。i X5. 上述算法的实施为了实施上述算法,还得做3个方面的改进。 (1外插方法外插可减少区域边缘的不连续性。设E 为外插算子。最简单的外插方法是常数外插方法。考虑投影数据为(u ,p p =,u ,分别为投影角和投影坐标值。常数外插定义为Eq 34(|(|a >u if a ,p au if u ,p =u ,p E指数形式外插定义为Eq 35(|(|a >u if a u a ,p

16、 a u if u ,p =u ,p E (exp(2 为参数,它可根据不同情况选择。复杂的外插方法还有多项式形式,余玄函数形式等等。此处就不赘述。在考虑外插后,应对式(Eq 33中其它一些变量业的作相应调整,使其定义域也包括感兴趣区域之外的区域。其中可改写为,可改写为,可改写为it R R ot R o R i X X,因此有Eq 36(too t p E R R p RE X +=式中,上式中只要求图像,oh ot o R R R =X在区域接近即可,而对区域的图像不做要求。因为在i i X o o 区域上无论怎样努力也不可能得到好的重建图像。(2对广义逆进行适当近似广义逆计算很费时。引理

17、 1告诉我们,因此有,是算子+oR IP R =1=P R 1P P 的普通逆。普通逆是广义逆的特殊情况。因此也有。显然。此处但是可考虑有 P R =+o o P R +=ho to o P P P Eq 37o o P R +或 Eq 38R H P R P R R o o o o =+其中分块算子的定义为HEq 39oi x ifX x if =X H 0式(Eq 36可改写成 Eq 40(t t p E R H P p E R X =(3考虑迭代算法因为上面采用了近似表达式(Eq 37,式(Eq 40很可能不够精确。迭代算法可补充上式为近似解的不足。以式(Eq 40为基础,进一步构成如下

18、迭代算法Eq 41t t pp =0(Eq 420(0(t p E R X =Eq 43(1(1(=n n tn X H P p pEq 44(n n pE R X =以上式(Eq 41, Eq 42给出迭代的起始值。 式(Eq 43, Eq 44迭代的主体。两部分合起来便是在文献12提出的迭代算法。其中为外插算子,为分块算子,EH(n X 给出X 的近似值。在区域i上(n X 将接近i X6.算法的其它形式同式(Eq 33类似,可得到其它形式的解。对式(Eq 30两边乘以可以推出oP Eq 45tot o t ot o p R P p R P =+其中,由此解出, =ho to o P P

19、P Eq 46tot o ot o t p R P R P p +=Eq 47t ot o ot o it i p R P R P I R X (+=同样道理,在式(Eq 45 用置换得 ot P o PEq 48t ot ot ot ot it i p R P R P I R X (+=交换上面矩阵的次序即 Eq 49otot ot ot ot ot ot ot R P R P R P R P +符号“”并不代表推导,而是某种猜测,它的正确性有待后面证明。因此有Eq 50t ot to ot to it i p R P R P I R X (+=考虑Eq 51ot ot ot ot ot o

20、t ot ot R R P P R P R P +=作下面变换 Eq 52+otot ot ot ot ot P P R R P P 得Eq 53t to to it i p P P I R X (+=式(Eq 33, Eq 47, Eq 48, Eq 50,Eq 53是不同形式但很类似的解,它们彼此不一定相等。其中式(Eq 33, Eq 53比较简单,因此是下面讨论的重点。由式(Eq 33到 式(Eq 53的“推导”是不严格的。其实由一个类似得到式(Eq 33的过程也可以得到式(Eq 53。不过那样过于麻烦,因此不再赘述。7.误差分析首先讨论式(Eq 33解的误差。式(Eq 3可写成 Eq

21、54oto i ti t X P X P p +=所以有 Eq 55o to ot i ti ot t ot X P R X P R p R +=考虑(Eq 11Eq 56oto ot t ot X P R p R =将上式代入式(Eq 33得Eq 57(o to ot ot t it i X P R R p R X +=定义解的误差为Eq 58 i i id X X X =重写式(Eq 12Eq 59(o to t it i X P p R X = 由定义(Eq 58并考虑式(Eq 59与式(Eq 57的差得Eq 60o to to ot ot it id X P P R R R X (=+

22、上式可用来估计式(Eq 33误差。从上式可看出Eq 61to ot ot to id PR R P X +其次考虑式(Eq 53解的误差。考虑同式(Eq 56相类似的原因考虑有Eq 62o to to t to XP P p P +上式中使用“”符号,因为上式紧紧是猜测,并没有给出一个完整地证明。同式(Eq 60类似,对应式(Eq 53的误差应为Eq 63o to to to to it id XP P P P R X (=+对应式(Eq 61有Eq 64to to to to id P P P P X +考虑SVD(Singular Value Decomposition方法矩阵可以被正交分

23、解 Eq 65UDVP to =其中为正交矩阵,即U Vm I UU =,n I VV =。,为单大小为,的单位矩阵。为对角矩阵,且有 m I n I n m D Eq 66=000d Dd 也为分块对角矩阵。,为SVD特异值。为分块零矩阵。相似的有=n d d d d .211d 2d n d 0Eq 67U VD P to +=此处Eq 68=+0001d D 因此有Eq 69U I U P P toto =+0001式中或分块单位矩阵。因此有1IEq 70to to to to P P P P =+由式(Eq 64有Eq 710=id X 这意味着在理论上该方法的误差为零。因此对应式(E

24、q 53的解有Eq 72i i X X =上式的成立基于(Eq 62的近似公式,因此受到一定限制。但至少上式应近似的成立。其实如果与有相同的SVD 特异值或有更多的特异值,类似上面关于(Eq 70的讨论,就有ot R toP Eq 73to to ot ot P P R R =+因此式(Eq 61中误差 为零。另一方面在此情况下可以由展开即id X to P ot R Eq 74kotk to R A P =应用式(Eq 56,式(Eq 62中近似号“”可用等号“=”代替。因此(Eq 72 的成立也就为精确而非近似的了。因此如果与有相同SVD 特异值,本文算法得到的解(Eq 33,Eq 53的

25、误差为零,即它们不是近似解而是精确解。ot R to P7.算法的合理性的讨论前面由式(Eq 33得到了迭代算法。但推导过程似乎有些武断。本节由式(Eq 53出发,看可推出什么算法。(Eq 53可改写成Eq 75t to o p P X+Eq 76(oto t it i X P p R X =对于方程(Eq 75中的可由下列线性方程组 o X Eq 77t o to pX P = 求解。而这个方程的解又可由迭代实现。外插对这类方程求解也会有帮助,因此也予以考虑。起始方程为Eq 78t o o pE R X=0(其中考虑了外插后应调整为,这是因为此算子作用的区域由改变。是投影空间到图像区域ot

26、R o R t o R o 的图像重建算子。迭代的投影为Eq 791(=n oto t n tX P p p o X 得迭代解可表示为Eq 80(n to n o p E R X =o X 可表示为下列的和Eq 81.(1(0(+=n o o o o X X X X o X 是式(Eq 77的解。对式(Eq 76考虑外插有Eq 82(o to t i i X P p E R X =同上的原理,考虑外插后调整为。是由投影空间到图像区域的算子。仔细分析发现式(Eq 78, Eq 79, Eq 80, Eq 81, Eq 82的迭代算法与式(Eq 42,Eq 43,Eq 44完全相同。在得到算法(E

27、q 42,Eq 43, Eq 44时我们作了两个近似,(A(Eq 37;(B改变式(Eq 40为式(Eq 42,Eq 43和Eq 44的迭代算法。这两个近似似乎很武断。算法(Eq 78, Eq 79, Eq 80, Eq 81, Eq 82得到的很自然,由此可更进一步加深对本文迭代算法合理性的理解。it R i R i R i8.结论本文总结了投影受限情况下图像重建的理论基础。讨论了该算法的误差。证明了如果与有相同的SVD 特异值,本文算法的误差为零,因此得到是精确解。本文最后详细由广意逆推导出迭代算法。由此说明此种情况下图像重建重投影的迭代算法合理性。至于与是否有相同SVD 特异值的证明将另

28、文讨论。ot R toP ot R toP参考文献1Kak A C and Slaney M 1988 Principles of computerized tomographic imaging IEEE Press2Rashid-Farrokhi F, Liu K J R, Berenstein C A and Walnut D 1997 Wavelet-based multiresolution local tomography IEEE Transactions on Image Processing 6 1412-303Nilsson M 2003 Local Tomography

29、at glance Thesis, Mathematics, Centre for Mathematical Sciences, LundUniversity, ISSN 1404-028X, ISBN 91-628-5741-X, LUTFMA-20074Seger M M 2002 Ramp filter implementation on truncated projection data, application to 3D linear tomography for logs Proceedings SSAB Symposium on Image Analysis, Lund, Sweden March 7-8, Editor Astrom5Ruchala K J, Olivera G H, Kapatoes J M, Reckwerdt P J, Mack T R 2002 Methods for improving limited field-of-view radiotherapy reconstructions using