Transcript Diapositiva 1

Image Reconstruction
from Projections
Antti Tuomas Jalava
Jaime Garrido Ceca
Overview

Reconstruction methods







Fourier slice theorem & Fourier
method
Backprojection
Filtered backprojection
Algebraic reconstruction
Diffractive tomography
Display of CT images
Tissue characterization with CT
Projection Geometry

Problem: Reconstructing 2D Image.


Given parallel-ray projections.
1D projection


 
 f x, y ds    f x, y  x cos  y sin   t dxdy
1
AB

(Radon transform).
Density distribution f x, y 
x cos  y sin   t1
Ray AB
p t1  

p t 
  
Integral evaluated for different values of the ray offset t1.
1D projection or Radon transform.
The Fourier Slice Theorem

1D Fourier Transform of 1D projection of 2D image
is equal to the radial section (slice or profile) of the 2D Fourier
Transform of the 2D image at the angle of the projection.

P w   p t exp j 2wt dt
F u, v 

P w    F u, v  F w, 
 
  f x, yexp j 2 ux  vydxdy

u  w cos
v  w sin 
The Fourier Slice Theorem

How to obtain f(x,y) applying Fourier Slice Theorem:

Assumption: we have projections available at all angles from 0º to 180º.
1.
2.
3.
From projections, we take their 1D Fourier transform.
Fill the 2D Fourier Space with the corresponding radial sections.
Take an inverse 2D Fourier transform to obtain f x, y 

Problem: finite number of projections available

Solution: Interpolation is needed in 2D Fourier space.
Backprojection


Simplest reconstruction procedure
Assumptions:



Rays: Ideal straight lines.
Image: dimensionless points.
Procedure

Estimate of the density at a point by simply summing (integrating )
all the rays that pass through it at various angles.
f  x, y  


p t d
0
where
t  x cos  y sin 
Problem:
•Finite number of rays per projection
•Finite number of projections
Interpolation is required.
Backprojection




BP produces a spoke-line pattern
blurring details.
Finite number of projections produces streaking artifacts.
Reconstructed image modeled by convolution between PSF
(impulse response) and the original image.
Solution:

Applying deconvolution filters to the reconstructed image.

Filtered BP technique.
Point density function
Filtered Backprojection

After some manipulations, we get:
f  x, y  

 q t d
0
where
Filter is
represented by
this function:
Ramp filter

q t    P w w exp j 2wt dw


In practice, smoothing window should be applied to reduce
the amplification of high-frequency noise.
Discrete Filtered Backprojection

Projection
p t 
in frequency domain is manipulated:
 2W  1 N / 2
 m   2 
P  k 

p

 exp  j mk


N
 N  2W m N / 2  2W  

Frequency
axis
discretized
Finite
number of
samples
Samples at the sampling
rate 2W
Discrete Filtered Backprojection

The filtered projection

q t 
q t    P w w exp j 2wt dw 

 m  2W
q 

 2W  N
2W
N
may then be obtained as:
N /2
2W 
 2W  2W

exp j 2k
t
k
N  N
N


 P  k
N / 2
 2W  2W
 2

P  k
exp j
m k
k

 N  N
 N

N / 2
N /2
 m  2W N / 2  2W  2W  2W   2 
q   
P  k
G k
k
 exp j mk

2
W
N
N
N
N
 



  N 
N / 2
•
Problem: control noise
enhancement
•
Solution we apply
hamming window:
 2W 
 2W  
G k   0,54  046cos k

 N 
 N W
Discrete Filtered Backprojection

Finally, we get this expression :
~
 L
f x, y    ql x cosl  y sin l 
L l 1

Algorithmic:
Measure projection.
2. Compute filtered projection.
3. Backproject the filtered.
4. Repeat 1-3 all projection angles
1.
Original
Back
projection
1°
Filtered
back
projection
10°
Filtered
back
projection
1°
Algebraic Reconstruction Techniques



Projections seen as set of simultaneous equations.
Kaczmarz method
 Iterative method.
 Implemented easily.
 Assumptions:
 Discrete pixels.
 Image density is constant within each cell.
Equations
w11 f1  w12 f 2    w1N f N  p1
w21 f1  w22 f 2    w2 N f N  p 2

wM 1 f1  wM 2 f 2    wMN f M  p M
Contribution factor of
the nth image element
to the mth ray sum.
Algebraic Reconstruction Techniques


Karzmarz method take the approach of successively and iteratively
projecting an initial guess and its successors from one hyperplane to
the next.
In general, the mth estimate is obtained from the (m-1)th estimate
as:
f

m 
f
m1
 f m1  wm  pm 
  wm
 
wm  wm


Because the image is updated by altering the pixels along each
individual ray sum, the index of the updated estimate or of the
iteration is equal to the index of the latest ray sum used.
Algebraic Reconstruction Techniques

Characteristics worth:
 ART
proceed ray by ray and it is iterative
 Small angles between hyperplanes


Large number or iterations
It should be reduced by using optimized ray-access schemes.
 M>N
noisy measurements
oscillate in the
neighborhood of the intersections of the hyperplanes.
 M<N under-determined.
 Any a priori information about image is easily
introduced into the iterative procedure.
Approximations to the Kaczmarz method

We could rewrite reconstruction step at the nth pixel level as:
True ray sum

p q 
f nm   max0, f nm1  m m 
Nm 

Computed ray sum
Number of pixels
crossed by the
mth ray.

Corrections could also be multiplicative:
f nm   f nm 1
pm
qm
Approximations to the Kaczmarz method

Generic ART procedure:
1.
2.
3.
4.
5.
0 
Prepare an initial estimate f
Compute ray sum qm
Obtain difference between true ray sum and the computed ray
sum and apply the correction.
Perform Steps 2 and 3 for all rays available.
Repeat Steps 2-4 as many times as required.
Original
1.
178 angles
dt =
1 voxel width
2.
3.
4.
5.
Original again
6.
Imaging with Diffraction Sources

Non ionizing radiation
 Ultrasonic
 Electromagnetic
(optical or thermal)

Refraction and diffraction

Fourier diffraction theorem
Imaging with Diffraction Sources
When an object,
f(x,y), is illuminated
with a plane wave the
Fourier transform of
the forward scattered
fields measured on
line TT’ gives the
values of the 2-D
transform, F(w1,w2),
of the object along a
circular arc in the
frequency domain, as
shown in the right
half of the figure.
Display of CT Images
 µ


HU  K 
 1
 µ




µ = measured attenuation coefficient.
µ = attenuation coefficient of water
When K = 1000 units are called Hounsfield Units

Air: -1000 HU
 Water: 0 HU
 Bone 1000 HU

Study





86 healthy infants aged 0-5 years
White matter: 15 HU to 22 HU
Gray matter: 23 HU to 30 HU
Difference between grey and white matter exactly 8 HU (In all measurements)
Boris P, Bundgaard F, Olsen A. Childs Nerv Syst. 1987;3(3):175-7
Microtomography

µ-scale CT
Volume: few cm3

Nanotomography already introduced.

Biomedical use:

 Both dead and alive (in-vivo) rat and mouse scanning.
 Human skin samples, small tumors, mice bone for
osteoporosis research.
Estimation of Tissue Components with CT

Manual segmentation of tumor by radiologist

Parametric model for the tissue composition
 Gaussian

mixture model
Method to estimate the parameters of the model
 EM
algorithm
Gaussian Mixture Model (i)

Fit M gaussian
kernels to intensity
histogram
Gaussian Mixture Model (ii)

Intensity value for voxel is a Gaussian random variable.

pi  x |  i  

  x  µi 2 
1
exp

2
2

2  i
i


M number of different tissues in tumor

i = the fraction of belonging to ith tissue (probability).



Parameters for ith tissue: i  µi ,  i
 Probability that voxel belonging to that tissue gets value x

M
i 1
i  1
Tumor as whole: PDF is a mixture of M Gaussians
  i ,..., M ,1 ,..., M 
px |    i 1 i pi x |  i 
M
Gaussian Mixture Model (iii)

Tumor as whole: PDF is a mixture of M Gaussians
  i ,..., M ,1 ,..., M 
px |    i 1 i pi x |  i 
M

Probability of parameter set 


p px | 
p | x  

p x 
 If nothing is known about 


p | x   c px | 
Find  that maximizes likelihood
N


L | x  ˆ px |    px j | 
j 1
Gaussian Mixture Model (iv)

Probability that jth voxel with value x j belongs
to the ith tissue type
pi |  px j | i,   i px j | i 
pi | x j ,  

px j , 
px j , 

EM algorithm (iterative, chapter 8) ->
inew , inew , inew
Ending Remarks





Some image manipulation tasks can be performed in 1D
in radon domain (edge detection etc.).
Reconstruction heavily dependent on reconstruction
algorithm (method).
MRI images are usually reconstructed with Fourier
method (according to book).
CT allows fast 3D imaging
So does MRI. MRI has better sensitivity especially with
soft tissues.