Gabor T. Herman. “Algorithms for Computed Tomography”
2000 CRC Press LLC. <>.
AlgorithmsforComputed
Tomography
GaborT.Herman
UniversityofPennsylvania
26.1Introduction
26.2TheReconstructionProblem
26.3TransformMethods
26.4FilteredBackprojection(FBP)
26.5TheLinogramMethod
26.6SeriesExpansionMethods
26.7AlgebraicReconstructionTechniques(ART)
26.8ExpectationMaximization(EM)
26.9ComparisonofthePerformanceofAlgorithms
References
26.1 Introduction
Computedtomographyistheprocessofreconstructingtheinteriorsofobjectsfromdatacollected
basedontransmittedoremittedradiation.Theproblemoccursinawiderangeofapplicationareas.
Herewediscussthecomputeralgorithmsusedforachievingthereconstructions.
26.2 TheReconstructionProblem
Wewanttosolvethefollowinggeneralproblem.Thereisathree-dimensionalstructurewhose
internalcompositionisunknowntous.Wesubjectthisstructuretosomekindofradiation,eitherby
transmittingtheradiationthroughthestructureorbyintroducingtheemitteroftheradiationintothe
structure.Wemeasuretheradiationtransmittedthroughoremittedfromthestructureatanumber
ofpoints.Computedtomography(CT)istheprocessofobtainingfromthesemeasurementsthe
distributionofthephysicalparameter(s)insidethestructurethathaveaneffectonthemeasurements.
Theproblemoccursinawiderangeofareas,suchasx-rayCT,emissiontomography,photon
migrationimaging,andelectronmicroscopicreconstruction;see,e.g.,[1,2].Alloftheseareinverse
problemsofvarioussorts;see,e.g.,[3].
Whereitisnototherwisestated,wewillbediscussingthespecialreconstructionproblemof

26.3 Transform Methods
The Radon transform has an inverse, R
−1
, defined as follows. For a function p of  and θ,

R
−1
p

(r, φ) =
1
2π
2

π
0

∞
−∞
1
r cos(θ − φ) − 
p
1
(,θ)ddθ ,
(26.1)
where p
1
(, θ ) denotes the partial derivative of p with respect to its first variable . [Note that it is
intrinsically assumed in this definition that p is sufficiently smooth for the existence of the integral
in Eq. (26.1)]. It is known [1] that for any function f which satisfies some physically reasonable

−∞
p(, θ )q



− , θ

d
(26.3)
are carried out, using a convolving function q (of one variable) whose exact choice will have an
important influence on the appearance of the final image. Second, our estimate f
∗
of f is obtained
by backprojection as follows:
f
∗
(r, φ) =

π
0
[
p ∗
Y
q
]
(
r cos(θ − φ),θ
)
dθ.
(26.4)

the mth estimate by a two-step process:
1. For −N ≤ n

≤ N, calculate
p
c

n

d,m

= d
N

n=−N
p
(
nd, m
)
q

n

− n

d

,
(26.6)
using the measured values of p(nd, m) and precalculated values (same for all m)of


,m

.
(26.7)
This is a discretization of Eq. (26.4). To do it, we need to interpolate in the first variable
of p
c
from the values calculated in Eq. (26.6) to obtain the values needed in Eq. (26.7).
Inpractice, once f
m+1
(r
j
,φ
j
) has been calculated, f
m
(r
j
,φ
j
)is no longer neededand the
computercanreusethesamememorylocationfor f
0
(r
j
,φ
j
),...,f
M−1

points whose locations will be precisely specified below); if they were collected otherwise, we need
to interpolate prior to reconstruction. If the function is to be estimated at an array of points with
rectangular coordinates {(id, j d)|−N ≤ i ≤ N,−N ≤ j ≤ N} (this array is assumed to cover the
object to be reconstructed), then the data function p needs to be known at points
(
nd
m
,θ
m
)
, −2N − 1 ≤ n ≤ 2N + 1, −2N − 1 ≤ m ≤ 2N + 1
(26.8)
c

1999 by CRC Press LLC
and at points

nd
m
,
π
2
+ θ
m

, −2N − 1 ≤ n ≤ 2N + 1, −2N − 1 ≤ m ≤ 2N + 1,
(26.9)
where
θ
m

(4N+3)d
tan θ
m

, −2N − 1 ≤ k ≤ 2N + 1 ,
−2N − 1 ≤ m ≤ 2N + 1
(26.11)
and at points (also in a rectangular coordinate system)

k
(4N+3)d
tan

π
2
+ θ
m

,
k
(4N+3)d

, −2N − 1 ≤ k ≤ 2N + 1 ,
−2N − 1 ≤ m ≤ 2N + 1 .
(26.12)
2. Windowing — At this point we may suppress those frequencies which we suspect to be
noise-dominatedbymultiplyingwith awindow function (correspondingtotheconvolving
function in FBP).
3. Separating into two functions — The sampled Fourier transform F of the object to be
reconstructed is written as the sum of two functions, G and H. G has the same values