VIEWS: 7 PAGES: 87 CATEGORY: Company Agreements POSTED ON: 8/4/2011 Public Domain
Computer Tomography: Computational theory and methods University of Bergen 15 October 2004, Bergen. Borislav V. Minchev Borko.Minchev@ii.uib.no http://www.ii.uib.no/ borko Department of computer science University of Bergen, Norway Computer Tomography: Computational theory and methods – p.1/28 Outline Introduction ¡ Fields of applications ¡ Principles of X-ray CT ¡ Reconstruction techniques in 2D ¡ - Transform based methods (FBP) - Algebraic reconstruction techniques Con-beam reconstructions ¡ - Circular - Helical Conclusions ¡ Computer Tomography: Computational theory and methods – p.2/28 Deﬁnition Tomography: technique for imaging cross-sections or slices from a set of external measurements of a spatially varying function. Computer Tomography: Computational theory and methods – p.3/28 Fields of applications X-ray CT ¡ - medical imaging - industrial nondestructive evaluation - airport screening - microtomography Computer Tomography: Computational theory and methods – p.4/28 Fields of applications X-ray CT ¡ - medical imaging - industrial nondestructive evaluation - airport screening - microtomography Magnetic Resonance Imaging (MRI) ¡ - electromagnetic waves Computer Tomography: Computational theory and methods – p.4/28 Fields of applications X-ray CT ¡ - medical imaging - industrial nondestructive evaluation - airport screening - microtomography Magnetic Resonance Imaging (MRI) ¡ - electromagnetic waves Transmission electron microscopy ¡ - an electron beem Computer Tomography: Computational theory and methods – p.4/28 Fields of applications X-ray CT ¡ - medical imaging - industrial nondestructive evaluation - airport screening - microtomography Magnetic Resonance Imaging (MRI) ¡ - electromagnetic waves Transmission electron microscopy ¡ - an electron beem CT principles are also applied in: ¡ - oceanography - astronomy - geophysics - optics Computer Tomography: Computational theory and methods – p.4/28 Some hisrory Radon 1917 ¡ - reconstruction of a function from its projections Computer Tomography: Computational theory and methods – p.5/28 Some hisrory Radon 1917 ¡ - reconstruction of a function from its projections Hounsﬁeld 1972 ¡ - patented the ﬁrst CT scanner Computer Tomography: Computational theory and methods – p.5/28 Some hisrory Radon 1917 ¡ - reconstruction of a function from its projections Hounsﬁeld 1972 ¡ - patented the ﬁrst CT scanner Cormack 1963 - 1964 ¡ - contributionto mathematics of X-ray CT Computer Tomography: Computational theory and methods – p.5/28 Some hisrory Radon 1917 ¡ - reconstruction of a function from its projections Hounsﬁeld 1972 ¡ - patented the ﬁrst CT scanner Cormack 1963 - 1964 ¡ - contributionto mathematics of X-ray CT Hounsﬁeld and Cormack shared the 1979 Nobel Prize for medicine ¡ Computer Tomography: Computational theory and methods – p.5/28 Some hisrory Radon 1917 ¡ - reconstruction of a function from its projections Hounsﬁeld 1972 ¡ - patented the ﬁrst CT scanner Cormack 1963 - 1964 ¡ - contributionto mathematics of X-ray CT Hounsﬁeld and Cormack shared the 1979 Nobel Prize for medicine ¡ since then the ﬁeld is steadily evolving and givis rise to never ceasing ﬂow of new ¡ algorithms Computer Tomography: Computational theory and methods – p.5/28 Principles of X-ray CT Computer Tomography: Computational theory and methods – p.6/28 Principles of X-ray CT X-ray traverses an ¡ object along stight line attenuated signals from ¡ various directions are recorded reconstruction algorithms ¡ are used to convert the projections Computer Tomography: Computational theory and methods – p.6/28 Projections Let represents the density of a two-dimensional object ¥£ ¨ ¢ ¤ § ¦ applied intensity form the X-ray source © the attenuated intensity of a ray as it propagates throught the object © along the line £ ¨ ¦ Computer Tomography: Computational theory and methods – p.7/28 Projections Let represents the density of a two-dimensional object ¥£ ¨ ¢ ¤ § ¦ applied intensity form the X-ray source © the attenuated intensity of a ray as it propagates throught the object © along the line £ ¨ ¦ then © ¤¥£ "¨ £ ¨ ¢ '% 0( § # & ) $ ¦ © ! where - projection angle - detector position Computer Tomography: Computational theory and methods – p.7/28 Projections £ ¨ ¤ £ "¨ ¢ ) ( § # ¦ ! Computer Tomography: Computational theory and methods – p.7/28 Projections £ ¨ ¤ £ "¨ ¢ ) ( § # ¦ ! The line integral can be written as Radon Transform 1 1 £ ¨ ¤¥£ 3¨ ¥£ ¤ £ ¨ "¨ ¢ " ) ( § 4 ¦ § ¤ ¦ § ¦ ¦ 21 21 where is the Dirac delta function and 3 ¤ £ ¨ &9 8 567 7 4 ¦ § § $ ¦ (the equation of the X-ray) Computer Tomography: Computational theory and methods – p.7/28 Fourier Slice Theorem Consider the two-dimensional FT of the object function 1 1 G " S 2 DE QI ¥£ ¨ ¥£ ¥¨ H P R F ¢ " @ A B ¤ § C ¤ § ¦ ¦ 21 21 and the one-dimensional FT of the projection £ ¨ ) ( 1 2 DE W ¥£ ¨ £ ¥¨ V F ) T " ) ( U C 21 Computer Tomography: Computational theory and methods – p.8/28 Fourier Slice Theorem Consider the two-dimensional FT of the object function 1 1 G " S 2 DE QI ¥£ ¨ ¥£ ¥¨ H P R F ¢ " @ A B ¤ § C ¤ § ¦ ¦ 21 21 and the one-dimensional FT of the projection £ ¨ ) ( 1 2 DE W ¥£ ¨ £ ¥¨ V F ) T " ) ( U C Simplest case, when ( ) 21 X X B Y ` 1 1 1 1 2 DE 2 DE ¥£ ¨ ¤ £ ¥¨ ¤ £ "¨ H H F F ¢ ¢ " " " @ X A § C ¤ § § § C ¤ ¦ ¦ ¦ 21 21 21 21 Y ` 1 1 2 DE 2 DE ¥£ "¨ ¤ £ ¥¨ A £ ¨ H H F F ¢ ) T " " ) ( ¤ § # C ¤ C ¤ b b a a ¦ ! 21 21 Computer Tomography: Computational theory and methods – p.8/28 Fourier Slice Theorem In general for any anglr ¥£ ¨ ¥£ ¨ U £ d¨ ) T ¦ @ @ &9 567 7 U U ¦ U c Computer Tomography: Computational theory and methods – p.8/28 Fourier reconstruction algorithm Take projection of an object function at angles ¡ ¦ ¦ e E f ¦c c c Fourier transform each projection to obtain the values of ¡ A £ ¨ @ B ¦ Recover the object function by using the inverse FT ¡ ¤ £ ¨ ¢ § ¦ 1 1 G " S DE PI ¤ £ ¨ A £ ¥¨ H R F ¢ " @ § B C A B ¦ ¦ 21 21 Computer Tomography: Computational theory and methods – p.9/28 Fourier reconstruction algorithm Take projection of an object function at angles ¡ ¦ ¦ e E f ¦c c c Fourier transform each projection to obtain the values of ¡ A £ ¨ @ B ¦ Recover the object function by using the inverse FT ¡ ¤ £ ¨ ¢ § ¦ 1 1 G " S DE PI ¤ £ ¨ A £ ¥¨ H R F ¢ " @ § B C A B ¦ ¦ 21 21 Note:inﬁnite number of projections on an interval with lenght allows exact g reconstuction. Computer Tomography: Computational theory and methods – p.9/28 Fourier reconstruction algorithm Take projection of an object function at angles ¡ ¦ ¦ e E f ¦c c c Fourier transform each projection to obtain the values of ¡ A £ ¨ @ B ¦ Recover the object function by using the inverse FT ¡ ¤ £ ¨ ¢ § ¦ 1 1 G " S DE PI ¤ £ ¨ A £ ¥¨ H R F ¢ " @ § B C A B ¦ ¦ 21 21 Discretising on the square h h h h ip s ip ip ip q q q q ¤ § r $ $ Ew Ew v v u y G G w HS G w S RS DE I ¤ £ ¨ x F ¢ @ t § C ¦ ¦ ¦ hE h h Ew Ew 2 v 2 v x a a where is even integer. The values of F(u,v) on a square grid are obtained from the available values along the radial lines by interpolation. Computer Tomography: Computational theory and methods – p.9/28 Reconstruction for Parallel Projection 1 1 G " S DE QI ¥£ ¨ A £ C ¨ H P R F ¢ " @ ¤ § B A B ¦ ¦ 21 21 ¥£ ¨ E 1 F &9 8 DE 567 7 ¤ § U £ ¥¨ V F ¦ " " @ C U U b b Computer Tomography: Computational theory and methods – p.10/28 ¢ ¤ ¥£ ¦ § ¨ b b E F F 21 1 b 21 1 1 21 1 ) T @ U £ U £ @ A £ ¨ ¦ ¦ U G ¥¨ C B S C C ¨ DE DE DE F F F V G V ¤ ¥£ H " U P QI 567 " R " S A " 8 7 § B &9 U ¨ " U " Reconstruction for Parallel Projection Computer Tomography: Computational theory and methods – p.10/28 ¢ ¤ ¥£ ¦ § ¨ b b E F F 21 1 b 21 1 1 21 1 ) T @ U £ U £ @ A £ ¨ ¦ ¦ U G ¥¨ C B S C C ¨ DE DE DE F F F V G V ¤ ¥£ H " U P QI 567 " R " S A " 8 7 § B b &9 F U ¨ 21 1 " U " 21 1 ) ( 0£ C ¥¨ 2 DE F V " U C DE F V " U " Reconstruction for Parallel Projection Computer Tomography: Computational theory and methods – p.10/28 ¢ ¤ ¥£ ¦ § ¨ b b b E F F F 21 1 e d b 21 1 21 1 1 21 1 ) T @ ) f( U £ U £ @ A £ 0£ ¦ ¨ ¦ ¨ U G C ¥¨ B S C C ¨ DE DE DE 21 1 F F F U V G V ¤ ¥£ H C " g U P QI G DE 567 " R F " S 2 S V A G " 8 7 § B 2" S b &9 U F " U ¨ 21 1 " j i i i ei h " U " 21 1 ) ( b 0£ F Y ¥¨ C 21 1 2 DE F ) ( V " £ k¨ U £ C DE $ F "¨ V " U ` " " Reconstruction for Parallel Projection Computer Tomography: Computational theory and methods – p.10/28 ¢ ¤ ¥£ ¦ § ¨ b b b b E F F F £ F 21 1 e d ) f( b lk 21 1 21 1 1 21 1 ¨ ) T @ ¥£ ¤ ) f( U £ U £ @ A £ 0£ 567 ¨ ¦ ¦ ¨ U G C ¥¨ B C ¨ S C DE DE 8 DE 21 1 F 7 F F § U V G V ¤ ¥£ H &9 C " g U QI P G DE 567 "¨ " R F " S 2 S V A G " 8 7 § B b 2" S b &9 F U F " ) m U ¨ ¤¥£ 21 1 " j i i i ei h " U 567 " 21 1 8 ) ( 7 b § 0£ F &9 Y C ¥¨ "¨ 21 1 2 DE F ) ( V " £ k¨ U £ C DE $ F "¨ V " U ` " " Reconstruction for Parallel Projection Computer Tomography: Computational theory and methods – p.10/28 ¢ ¤ ¥£ ¦ § ¨ b b b b E F F F F 21 1 £ d e ) f( b 21 1 lk 1 1 21 1 21 ) T Filtering ¨ @ ¥£ ¤ f( ) U £ U £ @ A £ 0£ ¦ 567 ¨ ¦ ¨ U G C ¥¨ B S C C ¨ DE DE DE 8 21 1 F F Backprojection 7 U F V G § V ¤ ¥£ H &9 C " g U P QI G DE 567 " "¨ R F " S 2 S V A G " 8 7 § B b 2" S b &9 F U F " ) m U ¨ ¥£ ¤ 21 1 " j i i i ei h " U 567 " 21 1 8 ) ( 7 b 0£ § F &9 Y ¥¨ C "¨ 2 DE 21 1 F ) ( V " £ k¨ U £ C DE $ F "¨ V " U ` " " Reconstruction for Parallel Projection Computer Tomography: Computational theory and methods – p.10/28 The FBP algorithm Step 1: Filtering 1 DE " £ ¨ ¥£ ¨ C V F ) m ) T U U U 21 Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering 1 DE " £ ¨ ¥£ ¨ C V F ) m ) T U U U 21 Use badlimited FT to compute i.e. U £ ¨ ) T Ew v 2 e n p qn p o p u G w S f v 2 DE ¥£ ¨ x F ) T ) T ( t U C ) $ ¦ ¦c c c ¦ ¦ o o p p p p Ew f 2 v a where is a frequency band limit and is some (large) number corresponding to the o number of the sampling points. (FFT) Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering 1 DE " £ ¨ ¥£ ¨ C V F ) m ) T U U U 21 Use badlimited FT to compute i.e. U £ ¨ ) T Ew v 2 e n p qn p o p u G w S f v 2 DE ¥£ ¨ x F ) T ) T ( t U C ) $ ¦ ¦c c c ¦ ¦ o o p p p p Ew f 2 v a where is a frequency band limit and is some (large) number corresponding to the o number of the sampling points. (FFT) Therefore Ew v ssr ssr qn p pn p n p o o o p p G w S f v DE r r r r rC r x F ) m ) T q t $ ¦ ¦c c c ¦ c o p p p Ew 2 v x a Inverse DFT of the product and . £ i ¨ ¨ i ¨ ) T o o £p £p Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering Superior results are obtained by deemphasizing the high frequencies with the Hamming window function . t Ew v ssr ssr qn p pn p n p n p o o o o p p p G w S f v DE r r r r r r x F ) m ) T t t C o p Ew 2 v x a Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering Superior results are obtained by deemphasizing the high frequencies with the Hamming window function . t qn p pn p qn p qn p o ) m u ( t l ) ¦ o o o p p p where is the inverse DFT of . £ vi ¨ i ¨ £ i ¨ ¨ qu o o t o £p £p Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering Superior results are obtained by deemphasizing the high frequencies with the Hamming window function . t qn p pn p qn p qn p o ) m u ( t l ) ¦ o o o p p p where is the inverse DFT of . £ vi ¨ i ¨ £ i ¨ ¨ qu o o t o £p £p Step 2: Backprojection Reconstruc by the formula ¤ £ ¨ ¢ § ¦ x g ¤¥£ ¨ ¤ £ ¨ ¢ z) m &9 8 t 567 7 § § ¦ y w e a where is the number of the projection angles . w Note The values of in Step 2 can be approxiamted from the values obtained in Step 1 ) m by interpolation (linear). Computer Tomography: Computational theory and methods – p.11/28 The FBP algorithm Step 1: Filtering Superior results are obtained by deemphasizing the high frequencies with the Hamming window function . t qn p pn p qn p qn p o ) m u ( t l ) ¦ o o o p p p where is the inverse DFT of . £ vi ¨ i ¨ £ i ¨ ¨ qu o o t o £p £p Step 2: Backprojection Reconstruc by the formula ¤ £ ¨ ¢ § ¦ x g ¤¥£ ¨ ¤ £ ¨ ¢ z) m &9 8 t 567 7 § § ¦ y w e a where is the number of the projection angles . w Note The values of in Step 2 can be approxiamted from the values obtained in Step 1 ) m by interpolation (linear). Computational cost: Backprojection summation —> | £ ¨ { Fast backprojections (Nilsson’96) —> E £ ¨ { 6 % } Computer Tomography: Computational theory and methods – p.11/28 Reconstruction for Fan Projections Equiangular Rays Let be a fan projection ¥£ ¨ f~ is the angle between the source and a reference axis. gives the location of a ray |SO| source to origin dist. - source to point dist. ¤ £ ¨ § ¦ In polar coordinates 4 £ ¨ u ¦ , ¨E ¨E 47 £ 4 £ &9 8 8 567 where . u $ Computer Tomography: Computational theory and methods – p.12/28 Reconstruction for Fan Projections Equiangular Rays 1 F ¥£ ¨ £ k¨ ¥£ "¨ ¢ " f( &9 8 567 7 ¤ § ¤ § ) $ ¦ b 21 E 1 F u 4 £ ¨ £ k¨ 4 £ £ ¨ "¨ ¢ u u " ( 567 ) $ $ ¦ p b 21 Computer Tomography: Computational theory and methods – p.12/28 where ¢ ¢ 4 £ ¤ ¥£ ¦ ¦ u § ¨ ¨ c c p c u b E b F ~ m F b E £ £ £ E u F ¨ ¨ ¨ 21 1 m 21 1 ) f( p £ £ u ~ ~ ( ) k¨ 7 "¨ n ¥£ ¥£ £ ¤ ¥£ &9 ¨ ¨ k¨ ¦ l 4 £ 567 pE ¥£ 567 k 567 ¨ £ 8 £ 7 Equiangular Rays § ¨ $ &9 u ¨ $ $ "¨ "¨ " " Reconstruction for Fan Projections Computer Tomography: Computational theory and methods – p.12/28 Weighted Backprojection Algorithm Equiangular Rays is the sampling interval of the projection Let are the angles at which projections are taken Step 1 For each , generate the corresponding modiﬁed projection £ ¨ y f~ ¦ y ~ £ ¨ £ ¨ ~ 567 Computer Tomography: Computational theory and methods – p.13/28 Weighted Backprojection Algorithm Equiangular Rays is the sampling interval of the projection Let are the angles at which projections are taken Step 1 For each , generate the corresponding modiﬁed projection £ ¨ y f~ ¦ y ~ £ ¨ £ ¨ ~ 567 E Step 2 Convolve with to generate the e y ~ £ ¨ £ ¨ £ ¨ k E corresponding ﬁltered projection (FFT) y ~ £ ¨ £ ¨ £ ¨ y m l Computer Tomography: Computational theory and methods – p.13/28 Weighted Backprojection Algorithm Equiangular Rays is the sampling interval of the projection Let are the angles at which projections are taken Step 1 For each , generate the corresponding modiﬁed projection £ ¨ y f~ ¦ y ~ £ ¨ £ ¨ ~ 567 E Step 2 Convolve with to generate the e y ~ £ ¨ £ ¨ £ ¨ k E corresponding ﬁltered projection (FFT) y ~ £ ¨ £ ¨ £ ¨ £ ¨ y m q l l ¦ where is a smoothing ﬁlter. £ ¨ q Computer Tomography: Computational theory and methods – p.13/28 Weighted Backprojection Algorithm Equiangular Rays is the sampling interval of the projection Let are the angles at which projections are taken Step 1 For each , generate the corresponding modiﬁed projection £ ¨ y f~ ¦ y ~ £ ¨ £ ¨ ~ 567 E Step 2 Convolve with to generate the e y ~ £ ¨ £ ¨ £ ¨ k E corresponding ﬁltered projection (FFT) y ~ £ ¨ £ ¨ £ ¨ £ ¨ y m q l l ¦ where is a smoothing ﬁlter. £ ¨ q Step 3 Perform a weighted backprojection u ¤ £ ¨ ¥£ ¨ ¢ y m w t § ¦ ¦ E ¤ £ ¨ ¦ § e ¦ e a where is the angle of a ray that passes through the point and . ¤ £ ¨ i w p § g ¦ Computer Tomography: Computational theory and methods – p.13/28 Reconstruction for Fan Projections Equally Spaced Detectors Let be a fan projection £ ¨ f~ # is the angle between the source and a reference axis. is the distance along the detector # |SO| source to origin dist. - the dist. source to projection of a pixel on the central ray In polar coordinates 4 £ ¨ u ¦ £ ¨ u &9 47 8 $ Computer Tomography: Computational theory and methods – p.14/28 where ¢ ¢ 4 £ ¤ ¥£ ¦ ¦ u § ¨ ¨ c c p c u b E b ~ m F F b E £ £ £ # # # E u F ¨ ¨ ¨ 21 1 m p 21 1 0( ) k u ~ ~ £ £ # f( ) £ £ £ k¨ "¨ # # # £ ¨ ¨ ¨ ¤ ¥£ k¨ l ¦ 4 £ 567 £ # 567 E ¨ £ 8 8 7 # E § $ &9 u ¨ Equally Spaced Detectors $ $ "¨ "¨ " " Reconstruction for Fan Projections Computer Tomography: Computational theory and methods – p.14/28 Algebraic reconstruction techniques Advantages - easy to handle different rays geometry - more adaptable to missing data - better image quality - metal artifacts are reduced - less radiation dose Computer Tomography: Computational theory and methods – p.15/28 Algebraic reconstruction techniques Advantages - easy to handle different rays geometry - more adaptable to missing data - better image quality - metal artifacts are reduced - less radiation dose Disadvantage Higher computational cost! Computer Tomography: Computational theory and methods – p.15/28 Image and Projection Representation The idea: Assume that the cross section consists of array of unknowns Computer Tomography: Computational theory and methods – p.16/28 Image and Projection Representation The idea: Assume that the cross section consists of array of unknowns - total number of cells is the value of ¤¥£¡ § ¢ ¦ in the -th cell ¨ - the -th line integral ray-sum ¬ «© ª Computer Tomography: Computational theory and methods – p.16/28 Image and Projection Representation The idea: Assume that the cross section consists of array of unknowns - total number of cells is the value of ¤¥£¡ § ¢ ¦ in the -th cell ¨ - the -th line integral ray-sum ¬ «© ª µ ¬ ² ³ ª± ª© ¤ ¤ ¤ ´ ´¤ ´ ¯ ° ® where is the total number of rays µ represents the contribution ª± -th -th of the cell to the ray integral ¨ ¬ Computer Tomography: Computational theory and methods – p.16/28 Projection methods Therefore we have to solve the following linear system of equations e ¢ E ¢ ¢ U8 8 U8 U ¶ e e eE e v v e e ¢ E ¢ ¢ U8 8 U8 U E¶ E e E E v v E . . . e ¢ E ¢ ¢ U8 8 U8 U ¶ · E· · v v · e for large values of and . ¸ Computer Tomography: Computational theory and methods – p.17/28 Projection methods Therefore we have to solve the following linear system of equations e ¢ E ¢ ¢ U8 8 U8 U ¶ e e eE e v v e e ¢ E ¢ ¢ U8 8 U8 U E¶ E e E E v v E . . . e ¢ E ¢ ¢ U8 8 U8 U ¶ · E· · v v · e for large values of and . ¸ Use iterative methods - projestions (Kaczmarz’37, Tanabe’71) y G S G S G S b b b Let is the initial guess G bS ¢ ¢ E ¢ ¢ e v ¦ ¦c c c ¹ Â ¼¯y » À¾ ½ º yV 2Á y ¿ G S G 2 eS ¢ ¢ Ã ¦ up U $ ¦c c c ¦ yV yV ¿ where U £ ¨ U ¦ U U e ÄÅ E v ¦c c c Computer Tomography: Computational theory and methods – p.17/28 Projection methods Theorem (Tanabe’71) ÇG S If there exists a unique solution to the system of equations, then ¢ Æ G S ÇG S f · Æ ¢ ¢ s9 % È Éf 1 Computer Tomography: Computational theory and methods – p.18/28 Projection methods Theorem (Tanabe’71) ÇG S If there exists a unique solution to the system of equations, then ¢ Æ G S ÇG S f · Æ ¢ ¢ s9 % È Éf 1 The rate of convergence depends of the angles between the hyperplanes - Full orthogonalization is computationally not feasible - Pairwise orthogonalization scheme (Ramakrishnan’79) - Carefully choose the order of the hyperplanes (Hounsﬁeld’72) Computer Tomography: Computational theory and methods – p.18/28 Projection methods Theorem (Tanabe’71) ÇG S If there exists a unique solution to the system of equations, then ¢ Æ G S ÇG S f · Æ ¢ ¢ s9 % È Éf 1 The rate of convergence depends of the angles between the hyperplanes - Full orthogonalization is computationally not feasible - Pairwise orthogonalization scheme (Ramakrishnan’79) - Carefully choose the order of the hyperplanes (Hounsﬁeld’72) When the process oscillate in the neighborhood of the intersections of the ¸Ë Ê hyperplanes When the process converges to a solution , such that ÆÇG S ¢ ¸ q Ê is minimized. G bS ÆÇG S ¢ ¢ $ Computer Tomography: Computational theory and methods – p.18/28 where where Therefore D ¢ G Consider again S D ¢ G ¢ G It can also be written as S 2 eS Ì 8 ¢ G Í ¢ f ¶ 2 eS G v D ¢ G a eU $ $ D ¢ G S e 2 S E Ì º ¹ S U f ¼ ¯y » U D ¢ G y D ¿V ½ Ð ¿ À¾ yÁ Ô WÒÓÑ 2 eS yV yV Projection methods ½ Í Î yÏ f yV 2 8 v y Õ a 2Á Ò U e ¦ up U Â ¢ D G ¢ f G ¦c cD S c ¦ c e Ã U 2 S ¦ cf ¦ ¦ up Ã ¦c c c ¦ ¦ up ¦c c c ¦ Computer Tomography: Computational theory and methods – p.19/28 Projection methods Consider again ¹ Â ¯y » À¾ ½ º ¼ yV 2Á y ¿ G S G 2 eS ¢ ¢ Ã ¦ up U $ ¦c c c ¦ ¿V yV y It can also be written as ¶ Ì D ¢ G S D ¢ G 2 eS $ ¦ Î ¦ up Ã ¦ up 8 U D ¦c c c ¦ ¦c c c ¦ v E Í eU f f a where G U 2 S e v G 2 S e Í ¢ ¢ Ì U cf f e f a Therefore D ¢ G S D ¢ G 2 eS D G S ¢ 8 ¦ where D ¢ G S yÁ yÏ 2 U cD WÒÓÑ Õ Ð yV Ò ½ Ô The main computational diffuculty is in the calculation, storage and retrieval of the weight coefﬁcients . U D Computer Tomography: Computational theory and methods – p.19/28 Algebraic Reconstruction Tech.(ART) Replace by 0 or 1, depending wether the center of the -th cell is within the -th ray. Î Ã Ê U D Thus we can use ¶ Ì D G S $ ¢ for all cells whos centers are within the -th ray. Ã -th v is the number of cells whose centers are within the ray. E Í Ã eU f f a Computer Tomography: Computational theory and methods – p.20/28 Algebraic Reconstruction Tech.(ART) Replace by 0 or 1, depending wether the center of the -th cell is within the -th ray. Î Ã Ê U D Thus we can use ¶ Ì D G S $ ¢ for all cells whos centers are within the -th ray. Ã -th v is the number of cells whose centers are within the ray. E Í Ã eU f f a Superior approxiamtion is Ê ¶ Ì D G S ¢ $ ¦ where is the lenght of the -th ray trought the reconstruction region. Ã Computer Tomography: Computational theory and methods – p.20/28 Algebraic Reconstruction Tech.(ART) Replace by 0 or 1, depending wether the center of the -th cell is within the -th ray. Î Ã Ê U D Thus we can use ¶ Ì D G S $ ¢ for all cells whos centers are within the -th ray. Ã -th v is the number of cells whose centers are within the ray. E Í Ã eU f f a Superior approxiamtion is Ê ¶ Ì D G S ¢ $ ¦ where is the lenght of the -th ray trought the reconstruction region. Ã - ART suffer from salt and pepper noise, caused by the approximations used for . U D Computer Tomography: Computational theory and methods – p.20/28 Algebraic Reconstruction Tech.(ART) Replace by 0 or 1, depending wether the center of the -th cell is within the -th ray. Î Ã Ê U D Thus we can use ¶ Ì D G S $ ¢ for all cells whos centers are within the -th ray. Ã -th v is the number of cells whose centers are within the ray. E Í Ã eU f f a Superior approxiamtion is Ê ¶ Ì D G S ¢ $ ¦ where is the lenght of the -th ray trought the reconstruction region. Ã - ART suffer from salt and pepper noise, caused by the approximations used for . U D - It is possible to reduce the noce by D ¢ G S relaxation (i.e. update a pixel by ) Ö Ö u q ¦ SIRT - change the cell values after going through all of the equations Computer Tomography: Computational theory and methods – p.20/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê Insted of solving v D ¢ c ¸ Ã ¦ up ¶ U D ¦ ¦c c c ¦ D e a We use again the Radon Transform 1 1 ¤ £ ¨ ¤ £ 3¨ ¥£ ¤ £ ¨ "¨ ¢ ¢ " ~ ¶ § § 4 § ¤ ¦ § ¦ ¦ ¦ 21 21 where is the equation of the -th ray and is the projection operator. ¤ £ ¨ ~ X Ã 4 § ¦ Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê Insted of solving v D ¢ c ¸ Ã ¦ up ¶ U D ¦ ¦c c c ¦ D e a We use again the Radon Transform 1 1 ¤ £ ¨ ¤ £ 3¨ ¥£ ¤ £ ¨ "¨ ¢ ¢ " ~ ¶ § § 4 § ¤ ¦ § ¦ ¦ ¦ 21 21 where is the equation of the -th ray and is the projection operator. ¤ £ ¨ ~ X Ã 4 § ¦ v Let be basis functions and let Ø× ¥£ ¨Ù ¤ § D D e ¦ a ¤ £Ú v ¤ £ ¨ ¨ Í ¤ £ ¨ ¢ ¢ D Ø t § § e § D D ¦ ¦ ¦ a then v ¤ £Ú ¥£ ¨ ¨ ¢ ¢ ~ ~ t ¶ ¤ § § Û D D ¦ ¦ D e a wher . ¥£ ¨ ~Ø Û ¤ § D D ¦ Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê Insted of solving v D ¢ c ¸ Ã ¦ up ¶ U D ¦ ¦c c c ¦ D e a We solve v c ¸ Ã ¦ up ¶ Û D D ¦ ¦c c c ¦ D e a - The two systems are the same for the pixel basis ÞÜ inside the -th pixel u Î Ý ¥£ ¨ Ø ¤ § D ¦ everywhere else X ß - Superior reconstructions are obtained by using bilinear elements pyramid shaped basis Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê bilinear elements insted of the traditional pixel basis Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê bilinear elements insted of the traditional pixel basis Correction terms are simultaneously applied to all the rays Ê as for the SIRT methods Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê bilinear elements insted of the traditional pixel basis Correction terms are simultaneously applied to all the rays Ê as for the SIRT methods Hamming window Ê emphasizes the corrections applied near the middle of a ray relative to those applied near the ends. Computer Tomography: Computational theory and methods – p.21/28 SART (Simultaneous ART) First reported by Anderson and Kak ’84 Reduce the error in the approximation of the ray integrals Ê bilinear elements insted of the traditional pixel basis Correction terms are simultaneously applied to all the rays Ê as for the SIRT methods Hamming window Ê emphasizes the corrections applied near the middle of a ray relative to those applied near the ends. The overall SART iterative scheme is · G S f Ö ¶ Û G S G S f f f e I $ ¦ Î ¦ up 8 Û D D D ¦ ¦c c c ¦ · v Í Í Û Û D D e D e e a a a where is the -th row of and are relaxation coefﬁcients. h Ö Ã X fË Û Computer Tomography: Computational theory and methods – p.21/28 More iterative algorithms Censor and Elfving’02 The idea: Use generalized projections and diagonal weighting matrices to reﬂect the sparsity of the matrix. h A family of diagonal matrices with diagonal elements · Deﬁnition: à× Ù â X DË á e a is called Sparsity Pattern Oriented (SPO) with respect to matix if: h á¸ · 1) 2) iff for every . Í e à © ¸ X X Ã ¦ up Û D D ¦c c c ¦ a Computer Tomography: Computational theory and methods – p.22/28 More iterative algorithms Censor and Elfving’02 The idea: Use generalized projections and diagonal weighting matrices to reﬂect the sparsity of the matrix. h A family of diagonal matrices with diagonal elements · Deﬁnition: à× Ù â X DË á e a is called Sparsity Pattern Oriented (SPO) with respect to matix if: h á¸ · 1) 2) iff for every . Í e à © ¸ X X Ã ¦ up Û D D ¦c c c ¦ a Diagonal Weighting (DWE) · G S f ¶ Û G S G S f f e I $ Ö c Î ¦ up 8 Û f D D D ¦ ¦c c c ¦ è Õ yê èév Í y ãä½ Ô ½ Ô èë y æ ç Ô yå ãæ ç Ô è y Computer Tomography: Computational theory and methods – p.22/28 More iterative algorithms Censor and Elfving’02 The idea: Use generalized projections and diagonal weighting matrices to reﬂect the sparsity of the matrix. h A family of diagonal matrices with diagonal elements · Deﬁnition: à× Ù â X DË á e a is called Sparsity Pattern Oriented (SPO) with respect to matix if: h á¸ · 1) 2) iff for every . Í e à © ¸ X X Ã ¦ up Û D D ¦c c c ¦ a Special case ÞÜ if e ì X Û D Ý åÆ D if X X Û ß D where is the number of the nonzero elements in column of . h Î # D Component Averaging (CAV) · G S f ¶ Û G S G S f f e I $ Ö c Î ¦ up 8 Û D f D D ¦ ¦c c c ¦ a íîv íÛ E Í # e í e a Computer Tomography: Computational theory and methods – p.22/28 More iterative algorithms Censor and Elfving’02 The idea: Use generalized projections and diagonal weighting matrices to reﬂect the sparsity of the matrix. h A family of diagonal matrices with diagonal elements · Deﬁnition: à× Ù â X DË á e a is called Sparsity Pattern Oriented (SPO) with respect to matix if: h á¸ · 1) 2) iff for every . Í e à © ¸ X X Ã ¦ up Û D D ¦c c c ¦ a Special case ÞÜ if e ì X Û D Ý åÆ D if X X Û ß D where is the number of the nonzero elements in column of . h Î # D Component Averaging (CAV) · G S f ¶ Û G S G S f f e I $ Ö c Î ¦ up 8 Û D f D D ¦ ¦c c c ¦ a íîv íÛ E Í # e í e a Jiang and Wang’03 ART, SART, CAV, DWE are special cases of the Landweber method Ê G S G S G S f f f e 2 e I ð £ ¨ h h Ö ï o 8 ¶ f $ Computer Tomography: Computational theory and methods – p.22/28 3D Reconstructions Slice by Slice ¡ Disadvantages - low speed - discrete nature - higher radiation dose Computer Tomography: Computational theory and methods – p.23/28 3D Reconstructions Slice by Slice ¡ Disadvantages - low speed - discrete nature - higher radiation dose Cone-Beam Reconstructions ¡ - Circular - Helical Computer Tomography: Computational theory and methods – p.23/28 Circular Cone-Beam Reconstructions Computer Tomography: Computational theory and methods – p.24/28 Circular Cone-Beam Reconstructions Apply fan-beam FBP to The Idea: each plane in the cone Feldkamp, Davis, Kress ’84 (FDK) Computer Tomography: Computational theory and methods – p.24/28 Circular Cone-Beam Reconstructions Apply fan-beam FBP to The Idea: each plane in the cone Feldkamp, Davis, Kress ’84 (FDK) A ray from the cone is described by &9 8 567 7 ¤ § £ 7 ¨ &9 &9 8 $ ¤7 567 4 § $ ò8 567 ñ where locate a ray in a fan £ ¨ ¦ specify the location ot the fan 4 £ ¨ ¦ - projection angle - fan angle —> (rotation by degrees around the -axes) ¤ £ ¨ £ ¨ ¦ §ñ ¦ #ñ ñ ¦ ¦ Computer Tomography: Computational theory and methods – p.24/28 The FDK Algorithm Based on the FBP Algorithm for equispatial rays E E n p 1 F u £ ¨ £ ¨ ¢ ó ¦ k " ¶ " ~ ¦ #ñ ¦ ¶ ¶ $ ¦ ô ¨E £ p E óE ¶ E 8 8 # # b $ $ 21 where and . ö Tõ ó 2 ñ ö Æ Computer Tomography: Computational theory and methods – p.25/28 The FDK Algorithm Based on the FBP Algorithm for equispatial rays E E n p 1 F u £ ¨ £ ¨ ¢ ó ¦ k " ¶ " ~ ¦ #ñ ¦ ¶ ¶ $ ¦ ô ¨E £ p E óE ¶ E 8 8 # # b $ $ 21 where and . ö Tõ ó 2 ñ ö Æ The Algorithm Step 1 Compute the modiﬁed projection ~ £ ¨ £ ¨ ó ó ~ ¶ ¦ ¶ ¦ ô E óE ¶ E 8 8 Step 2 Convolve with ~ £ ¨ £ ¨ ó k ip ¶ ¶ ¦ u ~ £ ¨ £ ¨ £ ¨ m ó ó k l ¦ ¶ ¶ ¶ ¦ p Step 3 Backproject over the 3D reconstruction grid E E ñ F £ ¨ £ "¨ ¢ m ¦ #ñ ¦ ¦ ¨E £ # # # b $ $ $ Computer Tomography: Computational theory and methods – p.25/28 The Tuy–Smith sufﬁciency condition Exact reconstruction is possible if all planes intersecting the object also intersect the source trajectory at least once. Computer Tomography: Computational theory and methods – p.26/28 The Tuy–Smith sufﬁciency condition Exact reconstruction is possible if all planes intersecting the object also intersect the source trajectory at least once. Circular trajectory does not satisfay this condition since a plane parallel to the trajectory may also intersect the object. Computer Tomography: Computational theory and methods – p.26/28 The Tuy–Smith sufﬁciency condition Exact reconstruction is possible if all planes intersecting the object also intersect the source trajectory at least once. Circular trajectory does not satisfay this condition since a plane parallel to the trajectory may also intersect the object. Exact reconstruction form helical cone-beam data is possible. Tam’95 Kudo, Noo, Defrise ’98 Computer Tomography: Computational theory and methods – p.26/28 Conclusions The challenges: Three dimentional images ¡ Better quality ¡ Faster algorithms ¡ Restricted in time ÷ Computer Tomography: Computational theory and methods – p.27/28 References K. Tanabe, Projection Method for Solving a Singular Systems of Linear Equations ¡ and its Applications, Numer. Math., 17, 203-214 (1971). H. Tuy, An Inversion Formula for Cone-Beam Reconstruction, SIAM J. Appl. Math. ¡ 42 546–552 (1983). L. Feldkamp, L. Davis and J. Kress, Practical Cone-Beam Algorithm, Journal of the ¡ Optical Society of America 1, 612-619 (1984). K. Tam, Three-Dimensional Computerized Tomography Scanning Method and ¡ System for Large Objects with Smaller Area Detectors, US Patent 5,390,112, Feb 14.(1995). S. Nilsson, Fast Backprojection, Technical report LiTH-ISY-R-1865, Linköping ¡ University (1996). A. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press. ¡ (1998). H.Kudo, F. Noo and M. Defris, Cone-Beam Filtered Backprojection Algorithm for ¡ Truncated Helical Data, Physics in Medicine and Biology 43, 2885–2909 (1998). Computer Tomography: Computational theory and methods – p.28/28 References G.Wang and M. Vannier, Computerized Tomography, (1998) available at: ¡ http://dolphin.radiology.uiowa.edu/ge/Teaching/ct/ct.html Y. Saad, H. van der Vorst, Iterative Solution of Linear Systems in the 20th Century, ¡ J. Comput. Appl. Math. 123, 1–33 (2000). F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction, ¡ SIAM, (2001). H. Turbell, Cone-Beam Reconstruction Using Filtered Backprojectio, Dissertation ¡ No.672, Linköping Sudies in Science and Technology (2001). Y. Censor and T. Elfving, Block-Iterative Algorithms with Diagonally Scaled Oblique ¡ Projections for the Linear Feasibility Problems, Siam J. Matrix Anal. Appl., 24, 40-58 (2002). M. Jiang and G. Wang, Convergence Studies on Iterative Algorithms for Image ¡ Reconstruction, IEEE Transactions on Medical Imaging , 22, 569-579 (2003). Computer Tomography: Computational theory and methods – p.28/28