Transcript Lecture 8 - University of Delaware

ELEG 479
Lecture #8
Mark Mirotznik, Ph.D.
Associate Professor
The University of Delaware
Summary of Last Lecture
X-ray Radiography
 Overview of different systems for projection radiography
 Instrumentation
 Overall system layout
 X-ray sources
 grids and filters
 detectors
 Imaging Equations
 Basic equations
 Geometrical distortions
 More complicated imaging equations
Hounsfield’s Experimental CT
Lets look at how CT works!
Example
y
= xray attenuation of 0
x

= xray attenuation of 2.5
= xray attenuation of 5

Our First Projection
g (  0, l )


l
Our First Projection
g (  0, l )


l

l

Rotate and Take Another Projection
g (  45 , l )
o
l
  45

o
l


Rotate and Take Another Projection
g (  90 , l )
o
l
  90
l
o



This is called a sinogram
l


This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l


Sinogram
This is called a sinogram
l

Sinogram

This is called a sinogram
l

Sinogram

This is called a sinogram
l

Sinogram

This is called a sinogram
l


The sinogram is what is measured by a CT machine.
The real trick is how do we reconstruct the
unknown image from the sinogram data?
Radon Transform
Id
) and f ( x, y )   ( x, y )
Io
x( s)  l  cos( )  s  sin(  )
y ( s)  l  sin(  )  s  cos( )
Given g (l , )   ln(
g (l , ) 

 f ( x(s), y(s))ds

In CT we measure g (l ,  )   ln(
Id
)
Io
and need to find f ( x, y )   ( x, y )
using
g (l , ) 

 f ( x(s), y(s))ds

x( s )  l  cos( )  s  sin(  )
y ( s )  l  sin(  )  s  cos( )
Radon Transform
In CT we measure
g (l ,  )
and need to find
f ( x, y )
We use
g (l ,  ) 
 
  f ( x, y)   ( x cos   y sin   l )dxdy
  
Reconstruction
The Problem
In imaging we measure g(l,) and need to determine f(x,y)
g (l ,  ) 
 
  f ( x, y)   ( x cos   y sin   l )dxdy
  

y
p
??
g(,l)
0
x
l
f(x,y)
Back Projection Method
A little trick that almost works!
Object
Back Projection
Back Projection Method
A little trick that almost works!
Back Projection
Object
We do this for every angle and then add
together all the back projected images
Back Projection Method
Step #1: Generate a complete an image for each projection (e.g. for
each angle )
b ( x, y )  g ( x cos( )  y sin( ), )
These are called back projected images
Step #2: Add all the back projected images together
p
fb ( x, y)   b ( x, y)d
0
Back Projection Method
Original object
Reconstructed object
Reconstructed Image
Kind of worked but we need to do better than
this. Need to come up with a better
reconstruction algorithm.
Projection-Slice Theorem
This is a very important theorem in CT imaging
First take the 1D Fourier transform a projection g(l,)
G (  , )  1D g (l ,  

 j 2pl
g
(
l
,

)
e
dl


g (l ,  )
Projection-Slice Theorem
This is a very important theorem in CT imaging
First take the 1D Fourier transform a projection g(l,)
G (  , )  1D g (l ,  

 j 2pl
g
(
l
,

)
e
dl


Next we substitute the Radon transform for g(l,)
g (l ,  ) 
 
  f ( x, y)   ( x cos   y sin   l )dxdy
  
G(  , ) 
  

   
f ( x, y )   ( x cos   y sin   l )e  j 2pl dxdydl
Projection-Slice Theorem
This is a very important theorem in CT imaging
First take the 1D Fourier transform a projection g(l,)
G (  , )  1D g (l ,  

 j 2pl
g
(
l
,

)
e
dl


Next we substitute the Radon transform for g(l,)
g (l ,  ) 
 
  f ( x, y)   ( x cos   y sin   l )dxdy
  
G(  , ) 
  

f ( x, y )   ( x cos   y sin   l )e  j 2pl dxdydl
   
Next we do a little rearranging


 j 2pl
G (  ,  )    f ( x, y )   ( x cos   y sin   l )e
dl dxdy
  
 

 
Projection-Slice Theorem
This is a very important theorem in CT imaging
Next we do a little rearranging


 j 2pl
G (  ,  )    f ( x, y )   ( x cos   y sin   l )e
dl dxdy
  
 

 
Applying the properties of the delta function
G (  , ) 
 





f ( x, y ) e  j 2p  x cos  y sin   dxdy
  
G (  , ) 
 

f ( x, y ) e  j 2p  x cos  y sin   dxdy
  
What does this look like?
Projection-Slice Theorem
This is a very important theorem in CT imaging
 

G (  , ) 




f ( x, y ) e  j 2p  x cos  y sin   dxdy
  
 

G (  , ) 
f ( x, y ) e  j 2p  x cos  y sin   dxdy
  
What does this look like?
 
This looks a lot like F (u, v) 

f ( x, y )e  j 2p  xu  yv dxdy
  
with
u   cos( ), v   sin(  )
Projection-Slice Theorem
This is a very important theorem in CT imaging
G(  , ) 
 

f ( x, y )e  j 2p  x cos  y sin dxdy
  
G (  , )  F (  cos( ),  sin(  ))
So what does this mean?
Projection-Slice Theorem
This is a very important theorem in CT imaging
G(  , ) 
 

f ( x, y )e  j 2p  x cos  y sin dxdy
  
G (  , )  F (  cos( ),  sin(  ))
Question: So what does this mean?
Answer: If I take the 1D FT of a projection at an angle the
result is the same as a slice of the 2D FT of the original
object f(x,y)
Projection-Slice Theorem
This is a very important theorem in CT imaging
G(  , ) 
 
  f ( x, y ) e
 j 2p  x cos  y sin 
dxdy
  
G (  , )  F (  cos( ),  sin(  ))
So what does this mean?
If I take the 1D FT of a projection at an angle the result is the same as a
slice of the 2D FT of the original object f(x,y)
Projection-Slice Theorem
If I take the 1D FT of a projection at an angle the result is the same as a
slice of the 2D FT of the original object f(x,y)
u   cos( ), v   sin(  )
2D FT
o
f(x,y)
g (l , o )
F(u,v)
o
G (  ,  o )  1D g (l ,  o ) 

 j 2pl
g
(
l
,

)
e
dl
o


u   cos( ), v   sin(  )
v
u
  u 2  v 2 ,   tan 1  
The Fourier Reconstruction Method
2D IFT

f(x,y)
g (l , o )
F(u,v)
o
u   cos( ), v   sin(  )
G (  ,  o )  1D g (l ,  o ) 

 j 2pl
g
(
l
,

)
e
dl
o


Take projections at all angles .
Take 1D FT of each projection to build F(u,v) one slice at a time.
Take the 2D inverse FT to reconstruct the original object based on F(u,v)
f ( x, y)  21D G(  , )
Image Reconstruction Using Filtered
Backprojection
g ( , l )
Filter
t

Backprojection
Filtered Back Projection
The Fourier method is not widely used in CT because of the
computational issues with creating the 2D FT from projections.
However, the method does lead to a popular technique called
filtered back projection.
 
f ( x, y ) 

j 2p  xu  yv 
F
(
u
,
v
)
e
dudv

  
In polar coordinates the inverse Fourier transform can be written as
2p 
f ( x, y) 
  F (  cos ,  sin  )e
0 
with u   cos( ),
v   sin(  )
j 2p  x cos  y sin 
dd
Filtered Back Projection
The Fourier method is not widely used in CT because of the
computational issues with creating the 2D FT from projections.
However, the method does lead to a popular technique called
filtered back projection.
In polar coordinates the inverse Fourier transform can be written as
2p 
f ( x, y) 

j 2p  x cos  y sin 
F
(

cos

,

sin

)
e
dd

0 
with u   cos( ),
v   sin(  )
From the projection theorem G (  ,  )  F (  cos( ),  sin(  ))
2p 
We can write this as f ( x, y) 

j 2p  x cos  y sin 
G
(

,

)
e
dd

0 
Filtered Back Projection
The Fourier method is not widely used in CT because of the computational
issues with creating the 2D FT from projections. However, the method does
lead to a popular technique called filtered back projection.
2p 
We can write this as f ( x, y ) 

j 2p  x cos  y sin  
G
(

,

)
e
dd

0 
Since g (l , )  g (l ,  p ) you can show
p 
f ( x, y)    G(  , )e j 2p x cos  y sin   dd
0 
which can be rewritten as
p


j 2p l
f ( x, y )     G (  ,  )e
 d 
d
0 
 l  x cos( )  y sin( )
Filtered Back Projection
verses Back Projection
A. Back Projection
b ( x, y )  g ( x cos( )  y sin( ), )
p
f b ( x, y )   b ( x, y )d
0
B. Filtered Back Projection
G (  , )  1D g (l ,  

 j 2pl
g
(
l
,

)
e
dl


p


j 2p l
f ( x, y)     G(  , )e
 d 
d
0 
 l  x cos( ) y sin( )
Filtered Back Projection Method
Filtered Back Projection
This always works!
Object
Digital Filter
1) take 1D FFT of projection
2) multiply by ramp filter
3) take 1D inverse FFT
4) make a back projection
Filtered Back Projection Method
Filtered Back Projection
Always works!
Object
Digital Filter
1) take 1D FFT of projection
2) multiply by ramp filter
3) take 1D inverse FFT
4) make a back projection
Filtered Back Projection Method
Filtered Back Projection
Always works!
Object
We do this for every angle
and then add together all
the filtered back projected
images
Digital Filter
1) take 1D FFT of projection
2) multiply by ramp filter
3) take 1D inverse FFT
4) make a back projection
Filtered Back Projection
verses Back Projection
A. Back Projection
Matlab Demo
b ( x, y )  g ( x cos( )  y sin( ), )
Your Assignment
p
(b)
Write a matlab function that reconstructs an image
f b ( x, y )  b ( x, y )d
using the filtered back projection method

0
B. Filtered Back Projection
G (  , )  1D g (l ,  

 j 2pl
g
(
l
,

)
e
dl


p


j 2p l
f ( x, y)     G(  , )e
 d 
d
0 
 l  x cos( ) y sin( )
Convolution Back Projection
From the filtered back projection algorithm we get
p


j 2p l
f ( x, y )     G (  ,  )e
 d 
d
0 
 l  x cos( )  y sin( )
It may be easier computationally to compute the inner 1D IFT
using a convolution
recall
1D1 F1 ( )  F2 ( )  f1 ( x)  f 2 ( x)
p

f ( x, y)   g (l , )  1D1 (  )
0

l  x cos( )  y sin( )
d
Convolution Back Projection
p

f ( x, y)   g (l , )  1D1 (  )

l  x cos( )  y sin( )
0
Let
c(l )  1D1 (  )
p
f ( x, y )   g (l ,  )  c(l ) l  x cos( )  y sin( ) d
0
d
Convolution Back Projection
c(l )  1D1 (  )
p
f ( x, y )   g (l ,  )  c(l ) l  x cos( )  y sin( ) d
0
The problem is c(l )  1D1 (  ) 



 e j 2pl d does not exist
Convolution Back Projection
c(l )  1D1 (  )
p
f ( x, y )   g (l ,  )  c(l ) l  x cos( )  y sin( ) d
0
The problem is c(l )  1D1 (  ) 


 e j 2pl d does not exist

The solution c~ (l )  1D1 (  W (  )) 


 W (  )e j 2pl d

where W (  ) is called a weighting function
p
f ( x, y )   g (l ,  )  c~ (l )d
0
c~ (l )  1D1 (   W (  )) 



 W (  )e j 2pl d
Convolution Back Projection
p
f ( x, y )   g (l ,  )  c~ (l )d
0
c~ (l )  1D1 (   W (  )) 


 W (  )e j 2pl d

Common window functions
 Hamming window
 Lanczos window (sinc function)
 Simple rectangular window
 Ram-Lak window
 Kaiser window
 Shepp-Logan window
•
•
•
•
Incorporated
linear array of
30 detectors
More data
acquired to
improve
image quality
(600 rays x
540 views)
Shortest scan
time was 18
seconds/slice
Narrow fan
beam allows
more
scattered
radiation to
be detected
•
•
•
•
Number of
detectors
increased
substantially (to
more than 800
detectors)
Angle of fan beam
increased to cover
entire patient
– Eliminated
need for
translational
motion
Mechanically
joined x-ray tube
and detector
array rotate
together
Newer systems
have scan times
of ½ second
2G
3G
Ring artifacts
• The rotate/rotate geometry of 3rd generation
scanners leads to a situation in which each
detector is responsible for the data
corresponding to a ring in the image
• Drift in the signal levels of the detectors over
time affects the t values that are
backprojected to produce the CT image,
causing ring artifacts
Ring artifacts
•
•
Designed to overcome
the problem of ring
artifacts
Stationary ring of
about 4,800 detectors
•
•
Designed to overcome
the problem of ring
artifacts
Stationary ring of
about 4,800 detectors
•
•
•
•
Developed
specifically for
cardiac
tomographic
imaging
No conventional
x-ray tube; large
arc of tungsten
encircles patient
and lies directly
opposite to the
detector ring
Electron beam
steered around
the patient to
strike the annular
tungsten target
Capable of 50msec scan times;
can produce fastframe-rate CT
movies of the
beating heart
•
•
•
•
Helical CT scanners acquire data while the table is moving
By avoiding the time required to translate the patient table, the total scan time
required to image the patient can be much shorter
Allows the use of less contrast agent and increases patient throughput
In some instances the entire scan be done within a single breath-hold of the
patient
Computer Assignment
1. Write a MATLAB program that reconstructs an image from its projections using
the back projection method. Your program should allow the user to input a
phantom object and a set (e.g. vector) of projection angle. Your program should
then: (a) compute the sinogram of the object (you can use Matlab’s radon.m
command to do this), (b) compute the reconstructed image from the sinogram
and vector of projection angles, (c) try your program out for several different
objects and several different ranges of projection angles
2. Do the same as #1 using the filter back projection method.
3. (grad students only) Do the same with the convolution back projection method