1. business and industrial
2. energy
3. oil

Back Projection Reconstruction
for CT, MRI and Nuclear Medicine
F33AB5
CT collects Projections
•
•
•
•
•
•
•
•
Introduction
Coordinate systems
Crude BPR
Iterative reconstruction
Fourier Transforms
Central Section Theorem
Direct Fourier Reconstruction
Filtered Reconstruction
To produce an image the
projections are back projected
Crude back projection
• Add up the effect of spreading each
projection back across the image space.
• This assumes equal probability that the
object contributing to a point on the
projection lay at any point along the ray
producing that point.
• This results in a blurred image.
Crude v filtered BPR
90
360
Crude BPR
Filtered BPR
Sinograms
r
r
q
Stack up projections
Solutions
• Two competitive techniques
– Iterative reconstruction
• better where signal to noise ratio is poor
– Filtered BPR
• faster
• Explained by Brooks and di Chiro in Phys.
Med. Biol. 21(5) 689-732 1976.
• This is repeated for
various angles of .
tra
ns
la
s
Xray
b
ea
m

an
sla
ted
Sample
X- r
ray
tub
e,
tr
– parallel rays, at
position r,
– across projection
at angle .
De
tec
tor
,
• Data collected as
series of
ted
Coordinate system
Attenuation of ray along a projection
• Attenuation occurs exponentially in tissue.

I  I 0e
  ( x ) dx
• (x) is the attenuation coefficient at position
x along the ray path.
Definition of a projection
s
Xray
b
ea
m
an
sla
ted
Sample
X- r
ray
tub
e,
tr
=   (r, s) ds
De
tec
tor
,
I
0
p(r ,  ) = ln
I
=   (x, y ) ds
tra
ns
la
ted
• Attenuation of a ray at position r, on the
projection at angle , is given by a line integral.
• s is distance along the ray, at position
r across

the projection at angle .
Coordinate systems
• (x,y) and (r,s) describe the distribution of
attenuation coefficients in 2 coordinate systems
tube
related by .
y
r (along projection)
r  x cos i  y sin i
S (along ray path)
detector

x
• where i =1..M for M different projection
orientations
• angular increment is  = /M.
Crude back projection
• Simply sum effects of back-projected rays
from each projection, at each point in the
image.
M
 (x, y) =  p(r ,  i) 
*
i =1
M
 (x, y) =  p (xcos  i + ysin  i ,  i) 
*
i =1
Crude back projection
• After crude back projection, the resulting
image, *(x,y), is convolution of the object
((x,y)) with a 1/r function.
Convolution
• Mathematical description of smearing.
• Imagine moving a camera during an
exposure. Every point on the object would
now be represented by a series of points on
the film: the image has been convolved with
a function related to the motion of the
camera

Iterative Technique
• Guess at a simulated object on a PxQ grid
(j, where j=1PxQ),
• Use this to produce simulated projections
• Compare simulated projections to
measured projections
• Systematically vary simulated object until
new simulated projections look like the
measured ones.
• For your scanner calculate
jj(r,i), the path length through
the jth voxel for the ray at (r,i)
• j need only be estimated once
at the start of the reconstruction,
• j is zero for most pixels for a
given ray in a projection
1
2 3
4
5 6 7 8
9 10 11 12
13 14 15 16
j=0
2=0.1
7=1.2
• The simulated projections are given by:
P Q
 (r ,  i ) =   j (r ,  i )  j
j 1
• j is mean simulated attenuation coefficient
in the jth voxel.
6
10
15
14
1 2 3
8 9 4
7 6 5
9
19
10
6
21
18
16 17 2
Object and projections
2 2 2
7 7 7
6 6 6
15 15 15
First ‘guess’
From Physics of Medical Imaging by Webb
15
13
To solve
• Analytically, construct P x Q simultaneous
equations putting (r,i) equal to the measured
projections, p(r,i):
•
p (r,  i ) =   j (r,
j
)

i
j
– this produces a huge number of equations
– image noise means that the solution is not exact and
the problem is 'ill posed’
• Instead iterate: modify j until (r,i) looks like
the real projection p(r,i).
Iterating
• Initially estimate j by projecting data in
projection at  = 0 into rows, or even simply
by making whole image grey.
• Calculate (r,i) for each i in turn.
• For each value of r and , calculate the
difference between (r,) and p(r,).
• Modify i by sharing difference equally
between all pixels contributing to ray.
6
10
15
14
1 2 3
8 9 4
7 6 5
16 17 2
Object
19
10
6
21
18
15
10
15
13
2 2 2
7 7 7
6 6 6
21/322/3 1
71/372/3 6
61/362/3 5
15 15 15
16 17 12
First ‘guess’
Next iteration
122/3
Fourier Transforms
• Imagine a note played by a flute.
• It contains a mixture of many frequency sound
waves (different pitched sounds)
• Record the sound (to get a signal that varies in
time)
• Fourier Transforming this signal will give the
frequencies contained in the sound (spectrum)
Time
Frequency
Fourier transforms of images
• A diffraction pattern is the Fourier
transform of the slit giving rise to it
y
ky
FT
kx
x
Central Section theorem
• The 1D Fourier transform of a projection
through an object is the same as a particular
line through 2DFT of the object.
• This particular line lies along the conjugate
of the r axis of the relevant projection.
y
ky
FT
kx
x
Projection
Direct Fourier Reconstruction
• Fourier Transform of each projection can be
used to fill Fourier space description of
object.
ky
y
Fp(r,1)
Fp(r,2)
kx
Inverse
Fourier
Transform
x
Direct Fourier Reconstruction
• BUT this fills in Fourier space with more
data near the centre.
• Must interpolate data in Fourier space back
to rectangular grid before inverse Fourier
transform, which is slow.
Relationship between object and
crude BPR results
• Crude back projection from above:

 (x , y) =  p(xcos  + ysin  ,  ) d  ,
*
o
• Defining inverse transform of projection as:
p(r ,  ) =
• then


Fp (k,  ) e
2 ikr
dk,
-
 
 (x, y) =   Fp (k ,  ) e
*
o -
2 ik(xcos + ysin )
[k]
dk d 
|k|
 
 (x, y) =   Fp (k ,  ) e
*
o -
2 ik(xcos + ysin )
[k]
dk d 
|k|
• The right hand side has been multiplied and
divided by k so that it has the form of a
2DFT in polar coordinates
– k conjugate to r
– k conjugate to r
– the integrating factor is kdrd  dxdy
Fp (k,  )
F (k x , k y) =
k
*
• Crude back projected image is same as the
true image, except Fourier amplitudes have
been multiplied by (magnitude of spatial
frequency)-1.
– Physically because of spherical sampling.
– Mathematically because of changes
coordimates.
in
Filtered BPR
• Multiplying 2 functions together is
equivalent to convolving the Fourier
Transforms of the functions.
• Fourier transform of (1/k) is (1/r)
• Multiplying FT of image with 1/k is same as
convolving real image with 1/r
• ie BPR has effect we supposed.
Filtered BPR
• Therefore there are two possible
approaches to deblurring the
crude BPR images:
• Deconvolve multiplying by f
(1/f x f = 1) in Fourier domain.
• Convolve with Radon filter in the
image domain, to overcome effect
of being filtered with 1/r by crude
BPR.