圖像重建的分塊迭代算法

1、圖像重建的分塊迭代算法趙雙任1,揚新鐵21多倫多DTI 成像行2西北工業(yè)大學摘 要:圖像重建和重新投影(Reconstruction and Re-projection的迭代算法已在投影視野受限(Limited Field of View, LFOV的情下使用。該情況下投影數(shù)據(jù)被截斷。此處投影采用的平行光束,扇形光束或錐形光束。該迭代算法對消除由于投影數(shù)據(jù)某處被截斷所產(chǎn)生的截斷膺像(Truncation Artifacts非常有效。最近發(fā)現(xiàn)該算法不僅對消除截斷膺像有效而且對消除條紋狀膺像也有一些效果。條紋狀膺像是在投影數(shù)據(jù)未被截斷情況下由濾波反投影(Filterred Backprojecti

2、on, FBP 算法所產(chǎn)生的膺像。條紋狀膺像與一個完整掃描所含的投影數(shù)及一個投影內(nèi)的點數(shù)有關(guān),該數(shù)目越小膺像就越嚴重。此發(fā)現(xiàn)引起了如下猜想,即該迭代算法或許也適用于投影數(shù)據(jù)未被截斷的情況。本文對該迭代算法也做了進一步改進,使它更有利于減小條紋狀膺像。主要的改進是把圖像劃分成較小的區(qū)域,以及剔除原迭代算法中用于數(shù)據(jù)被截斷時的外插(Extrapolation技術(shù)。本文投影數(shù)據(jù)采用平行光束。迭代算法包含兩個或者更多個圖像重建過程。第一個過程采用FBP算法。第二次圖像重建過程中第一次重建的圖像被切成小塊;每一小塊圖像都被重新投影并再次重建;最終合并成一個完整圖像。分塊迭代算法的計算結(jié)果與FBP 算法的

3、以及傳統(tǒng)的迭代法的結(jié)果做了比較。關(guān)鍵詞:膺像,斷層掃描,迭代,圖像重建,投影,視野,截斷,平行光束,CT一.簡介對平行光束扇形光束和錐形光束的CT(Computer Tomography系統(tǒng)已經(jīng)有許多圖像重建算法1-7。如果探頭寬度足夠大,物體全部落在CT 探頭的視野之內(nèi),測量所得到的投影數(shù)據(jù)是完整的。相反如果物體的一部分落在CT 探頭的視野之外,測量所得到的投影數(shù)據(jù)在探頭的邊界處被截斷。由截斷了的投影數(shù)據(jù)所得的重建圖像包含截斷膺像(Truncation Artifacts。將局部斷層掃描(Local Tomography方法8-12,圖像重建再投影(Reconstruction and Re

4、projection的方法13-16和外差(Extrapolation方法17-19相結(jié)合,文獻20研究了探頭的視野受限,投影數(shù)據(jù)被截斷的情況下的圖像重建問題,并建立了一個新穎的迭代算法。該算法將物體分成兩部份。一部分為感興趣的區(qū)域,另一部分是感興趣之外的區(qū)域。感興趣的區(qū)域剛好全部落在探頭視野之內(nèi),該區(qū)域的圖像重建是探頭視野受限情況下所關(guān)注的。該迭代算法由至少兩次圖像重建構(gòu)成。第一次采用FBP方法包括對在探頭之外未知投影數(shù)據(jù)進行外差。用FBP算法得到感興趣區(qū)域之外部份的圖像。將該部分圖像再次投影。將原投影數(shù)據(jù)與再次投影數(shù)據(jù)相結(jié)合。使用由此所得的投影數(shù)據(jù)對感興趣的區(qū)域進行第二次圖像重建。該算法對

5、減少截斷膺像非常有效。最近發(fā)現(xiàn),其實該算法也對投影數(shù)據(jù)未被截斷的情況起作用。跟據(jù)這一發(fā)現(xiàn),本文將該迭代算法進一步改造,使之更適應(yīng)投影數(shù)據(jù)未被截斷的情況。原算法中適用投影數(shù)據(jù)被截斷部份比如外差方法被剔除。原算法中感興趣的區(qū)域為圓形,改造后成成為方形。原算法中只有一塊感興趣的區(qū)域。改造后圖像所在區(qū)域被分成許多子區(qū)域。這些子區(qū)域都是感興趣的區(qū)域。迭代算法被用于每一個子區(qū)域。子區(qū)域內(nèi)的圖像最終被拼接成整塊圖像。由此構(gòu)成所謂分塊迭代算法。與FBP方法比較,分塊迭代算法的抑制膺像作用很顯著。但因為其他迭代算法也對膺像有一定抑制作用。本文還將分塊迭代算法與一種傳統(tǒng)迭代算法進行了比較。此處提到傳統(tǒng)迭代方法是一

6、種常見的線性方程迭代解法。盡管傳統(tǒng)迭代方法對抑制膺像也有貢獻,但分塊迭代算法的抑制膺像作用要更上一層。二.傳統(tǒng)的圖像重建迭代算法設(shè)X 為二維物體,或稱原像。p 為測量的或仿真的投影數(shù)據(jù),此處假定CT探頭寬度足夠大因此投影數(shù)據(jù)未被截斷。P 為平行光束投影算子或Radon變換1。投影數(shù)據(jù)與投影算子之間關(guān)系為X P =p (1 如果X 已知,上式可對投影數(shù)據(jù)進行仿真。若X 未知,一般X 可由濾波反投影(FBP方法1求得(p R =X 0 (2 此處(0X 是對原像的重建。為與后面迭代算法保持一致,此處采用上標0表示第零次迭代。R 為FBP 圖像重建算子,它先對投影數(shù)據(jù)進行濾波,再將濾波后的數(shù)據(jù)反投影

7、。FBP 算法求得的解已廣泛在臨床使用。對此本文不再贅述。值得一提的是FBP 算法所得的解(0X 仍含有一些膺像。該膺像呈條文狀且與一個完整CT 掃描所包含的投影數(shù)目及每一個投影上的點數(shù)及有關(guān)。投影數(shù)目和投影上的點數(shù)越少,膺像越嚴重。減少膺像的辦法之一是迭代算法。為了同本文后面要講述的分塊迭代算法有所區(qū)別,稱其為傳統(tǒng)迭代算法。下面對傳統(tǒng)迭代算法作簡單敘述。該算法把上述結(jié)果式(2重新投影(1p (01X P =p (3再從投影數(shù)據(jù)減去重新投影得到兩者之差稱為投影誤差d p(1p p =p d (4 再對投影誤差進行圖像重建,得到圖像的改進量d Xd d p R =X(5 將第一次重建的圖像與圖像

8、改進量相加便得傳統(tǒng)迭代解 (d X +X =X 01 (6將上面式(3-6步驟合并得(001PX p R +X =X (7 考慮到迭代可能不至一次,該式可改寫為 (11n n n X P p R +X =X (8式(2與式(8結(jié)合給出傳統(tǒng)的圖像迭代算法。該算法同解線性方程組時在已知近似解算子進一步求精確解的迭代算法一樣。這種算法并沒有什么奇特之處。此處給出它的目的是為了提供對下面分塊迭代算法的一個參考。三. 圖像重建的分塊迭代算法由于CT探頭視野寬度不足可引起投影數(shù)據(jù)被截斷。對截斷的投影作FBP圖像重建可產(chǎn)生嚴重截斷膺像。文獻21提出了抑制截斷膺像的迭代算法。對上述迭代算法稍加改造,可得到分塊

9、迭代算法。分塊迭代算法適用于投影數(shù)據(jù)未被截斷的情況。下面詳細講述分塊迭代算法。為了討論方便先定義分塊算子i T 。i T 可從整塊圖像中切出一小塊感興趣的區(qū)域。iT 的定義為(9 =i i i ROI x if ROI x if x T 01(其中x 是圖像中的一個像素點。ROI(Region of Interest 表示感興趣的區(qū)域;i ROI 為第i 塊感興趣的區(qū)域,本文假定它是方形或矩形。例如原圖像可以被分成3X3=9 塊。設(shè)i X 是第i 塊感興趣區(qū)域上的圖像,i X 可表示為X T =X i i(10 感興趣區(qū)域i ROI 以外處的圖像可表是為 X T (I =Y i i (11 此

10、處I 為恒等算子,即滿足X I =X 。i Y 是從圖像X 中挖去i X 后的圖像。分塊迭代的關(guān)鍵是將感興趣區(qū)域i ROI 以外的區(qū)域上的圖像重新投影 (i PY =p 1 (12再從測量投影數(shù)據(jù)中減去重新投影得到兩者的差d p(1p p =p d (13 再將此進行圖像重建得到對應(yīng)圖像的差d Xd d p R =X (14 值得一提,此差含義與式(5不同。此處圖像的差d X 在i ROI 以外的區(qū)域與式(5相同。但在i ROI 區(qū)域上同式(6或式(10有些接近,是對X 的圖像的某種重建?？紤]上述原因,不妨利用分塊算子i T 切去i ROI 以外的圖像得到在第i 塊感興趣區(qū)域內(nèi)圖像(d i i

11、 X T =X 1(15再將所有分塊區(qū)域上的圖像拼接起來可得到一個總的圖像 (N =i i X =X111(16 以上給出了分塊迭代算法的步驟。合并以上步驟,并用0(X 取代式(10中的X 得(011X T (I P p R T =Xi N =i i (17 式中0(X 由式(2定義。式(2和式(17構(gòu)成分塊迭代算法。由此算法重建的圖像在塊與塊之間的邊界處有一接縫。為了消除接縫,對上述公式作適當改進(X T P(I R(p T =X+i N =i i 011 (18此處+i T 稱為加邊分塊算子,其定義為MARGIN+(ROI x if MARGIN+(ROI x if =(xT i i +i

12、 01(19MARGIN 是邊緣的意思,上述算子對感興趣區(qū)域增加一圈邊緣。邊緣的寬度可考慮為2到10個像素。由式(18可消除塊與塊之間的接縫?？紤]到迭代可能不至一次,式(18可進一步改寫為 (11n +i N =i i n X T P(I p R T =X(20將此式展開 (11111n +i N =i i n N =i i N =i i n X RPT T +RPIXT Rp T =X(21 考慮一個整塊圖像切塊后在拼接起來應(yīng)與原圖像完全相同。即如果X 為整塊圖像則有 XT =X N =i i 1 (22 式(21中前兩項Rp 和(1n RPIX 都是一整塊圖像,對此可應(yīng)用式(22簡化。最后

13、一項中(1n +i X RPT 是N 個不同的圖像因此不能簡化。由此式(21化為(111n n +i N =i i n PX p R +X RPT T =X(23定義一個分塊投影重建算子+i N =i i RPT T =K 1 (24 該算子對圖像作如下處理,(1對圖像進行分塊,并考慮適當邊緣,即作+i T 算子運算;(2作投影運算,即作P 算子運算;(3圖像重建比如FBP 算法,即作R 算子運算;(4用分塊運算i T 切除感興趣區(qū)域以外的圖像;(5合并分塊,即N =i 1的運算。利用式(24,式(23可改寫成(11n n n X p R +KX =X (25式(2和式(25構(gòu)成分塊迭代算法。

14、式(25同式(8非常相似。式(25中僅在第一項中多了分塊投影重建算子K ,它是該算法的關(guān)鍵。四.結(jié)論本文提出了用于CT 圖像重建的分塊迭代算法。并將該算法與FBP 算法和傳統(tǒng)迭代算法作了比較?？梢钥隙?分塊迭代算法的膺像總體上比其它兩種算法小,因此是一種非常有前途的圖像重建算法。參考文獻1 A. C. Kak and Malcolm Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, 19882 L. A. Feldkamp, L. C. Davis, and j. W. kress, Practical co

15、ne-beam algorithm. J. Opt.Soc. Am. A/Vol. 1, No. 6/June 1984, pp 612-619.3 B. D. Smith, Image Reconstruction from cone-beam projections: necessary and sufficient conditions. IEEE transactions onMedical Image, Vol. MI-4, No. 1 Mach 1985, pp 14-25.4 P. Grangeat, Mathematical framework of cone beam 3D

16、reconstruction via the first derivative of the Radon transform, G. T.Herman, A.K. Louis, and F. Natterer (eds., in Mathematical methods in Tomography, Lecture notes in Mathematics No. 1497, Springer Verlag, 66-97.5 Defrise, M. Clack, R. (1994 A cone-beam Reconstruction algorithm using shift Variant

17、Filtering and cone-beam backprojection,IEEE trans. On Medical Imaging. 13. 186-195.6 Shuangren Zhao, H. Halling, Image reconstruction for Fan beam X-ray Tomography Using a new Integral Transform pair.7 Shuangren Zhao, H. Halling, Reconstruction of cone beam projection with free source path by a Gene

18、ralized Fourier method,Proceedings of the 1995 International meeting on fully three-Dimensional Image Reconstruction in radiology and Nuclear medicine, p323.8 Adel Faridani, erk L. Ritman, and kennan T. Smith: “Local Tomography” SIAM. J. APPL. MATH. Vol 52, No 2 pp 459-484,April 1992.9 Alexander Kat

19、sevich, Cone beam Local Tomography, SIAM J. APPL. MATH. Vol. 59, No. 6, pp. 2224-2246. 1999.10 Alexander Katsevich, Theoretically exact Filtered backprojection-type Inversion algorithm for spiral CT. SIAM J. App. MATH.Vol. 62, No. 6, pp. 2012-202611 F. Rashid-Farrokhi, K. J. R. Liu, C. A. Berenstein

20、, and D. Walnut, “Wavelet-based Multiresolution Local Tomography”, IEEETransactions on Image Processing, 6 (1997 pp. 1412-1430.12 khi, K. J. R. Liu, C. A. Berenstein, and D. Walnut, “Wavelet-based Multiresolution Local Tomography”, IEEE Transactions onImage Processing, 6 (1997 pp. 1412-1430.13 Nassi

21、, M.; Brody, W.R.; Medoff, B.P.; Macovski, A Iterative reconstruction-reprojection: an algorithm for limited data cardiac-computed tomography. IEEE trans. Biomed. Engineering BME-295:333-340; 1982.14 Paul S. Cho, Anthony D. Rudd and Roger H. Johnson, Cone-beam CT from width truncated projections, Computerized MedicalImaging and Graphics, Vol 20, No. 1, pp. 49-57, 1996. p49-57.15Mattias Nilsson, Thesis, Department: Mathematics, Lund Institute of Technology, Lund University, Centre for Mathematical Sciences.16J. H. Kim, K. Y. KWAK, S. B. Park, and Z. H. Cho, Projection Space Iteration Reconstr