10/19/2006

ENGG 167

MEDICAL IMAGING
Lecture 12: Friday, Oct. 19 Image Reconstruction from Projections
Reference: Chapters 12 & 13, The Essential Physics of Medical Imaging, Bushberg Computed Tomography, Kalender, Verlag, 2000. Chapter 12, Intermediate Physics for Medicine and Biology, 3rd Ed., Hobbie. Chapter 3, Principles of Computerized Tomographic Imaging, Kak and Slaney, IEEE Press Medical Physics and Biomedical Engineering, Brown, et al, IoP Publishing.

Image Reconstruction
1) 2) 2) 3) 2-D Fourier transform review Filtering in space and frequency domain Backprojection Filtered Backprojection

Ref: Hobbie

10/19/2006

1) Basics Fourier Transform of an image f(x,y) Transformed image is described in k-space, where kx= 2/x, ky=2/y are spatial frequencies

Fourier amplitude

Real space
Ref: Hobbie
3

1) What does the DFT image represent ?

Highest Frequencies kx,max = (Nx/Sx) = 1/[2(pixel width)] ky,max = (Ny/Sy) lowest frequencies kx,min = 1/Sx ky,min = 1/Sy Maximal negative frequency kx,max = - (Nx/Sx) = -1/[2(pixel width)] ky,max = - (Ny/Sy)

N number of pixels S size of image (mm)

Ref: Hobbie
4

10/19/2006

2) Fourier Slice Theorem example

X=abs(fft2(F));imagesc((X(:,:,1)));colorbar;colormap('bone');

X=abs(fft2(F));imagesc(log(X(:,:,1)));colorbar;colormap('bone');

F=imread('abdomen visible human gray.jpg'); imagesc((F));colorbar; colormap('bone');

Note DFT has same number of pixels, but two images magnitude and phase

X=abs(fft2(F));imagesc(fftshift(log(X(:,:,1))));colorbar;colormap('bone');

P=angle(fft2(F));imagesc(fftshift(P(:,:,1)));colorbar 5

3) Fourier Transform of an image f(x,y) in MATLAB

Note DFT image has both amplitude and phase information

F=imread('kpaulsen.jpg');imagesc((F));colorbar

X=abs(fft2(F));imagesc(log(X(:,:,1)));colorbar P=angle(fft2(F));imagesc(P(:,:,1));colorbar

K-space image with K=0 at center

X=abs(fft2(F));imagesc(fftshift(log(X(:,:,1))));colorbar P=angle(fft2(F));imagesc(fftshift(P(:,:,1)));colorbar

10/19/2006

1) Review: Fourier Transform of an image f(x,y) in MATLAB

A=fft2(F);F1=uint8(abs(fft2(A./20500))); F2=imrotate(F1,180);imagesc(F2)

2) Filtering space domain convolution

Convolution integral g(x) = filtered data g(x) = initial data h(x-x) = filter function

g '( x) =

g ( x ')h( x x ')dx '

g '( x) = g ( x) h( x)

[1 2 1]

smoothing filter

g=[1 3 4 2 3 5 6 7 8 9 5 6 7 8 9 7 5 4 3 2 3 2 4 2 1 1 1 ]; plot(g); h=[1 2 1]; gp=conv(g,h); Figure; plot(gp)

[0.25 0.5 0.25] =

h=[.25 .5 .25]; gp=conv(g,h); plot(gp) 8

10/19/2006

2) Filtering space domain convolution

Convolution integral g(x) = filtered data g(x) = initial data h(x-x) = filter function

g '( x) =

g ( x ')h( x x ')dx '

g=[1 3 4 2 3 5 6 7 8 9 5678975432 3 2 4 2 1 1 1 ]; plot(g); h=[1 2 1]; gp=conv(g,h); Figure; plot(gp)

g '( x) = g ( x) h( x)

[-1 2 -1] =

[-1 2 -1] =

[-1 4 -1] =

g=[1 2 2 2 2 7 7 7 77777444 4 4 2 2 2 2 2 2 2 2 1 ]; plot(g) gp=conv(g,h); plot(gp) h=[-1 4 -1]; gp=conv(g,h); plot(gp)

2) Filtering of Images space domain

Convolution integral g(x,y) = image f(x,y) = test object h(x-x,y-y) = filter function

g '( x, y ) =

g ( x ', y ')h( x x ', y y ')dx ' dy '

g '( x, y ) = g ( x, y ) h( x, y )
X=imread('happyface.jpg'); imagesc(X); B=[1 2 1; 2 4 2; 1 2 1]; Y=conv2(X,B); Imagesc(Y);colormap(bone);

121 242 121

-1 -2 -1 -2 8 -2 = -1 -2 -1 1 1 1 0 0 0 = -1 -1 -1

imagesc(X); B=[-1 -2 -1; -2 8 -2; -1 -2 -1]; Y=conv2(X,B); Imagesc(Y);colormap(bone);

imagesc(X); B=[1 2 1; 0 0 0; -1 -1 -1]; Y=conv2(X,B); Imagesc(Y);colormap(bone); Prewitt Filter 10 h = fspecial(type,parameters) - predefined filters of any shape

10/19/2006

2) Filtering of Images frequency domain

The Fourier transform of a convolution is the multiplication of the transformed functions

This is also true in 2-D (images)

g '( x, y ) = g ( x, y ) h( x, y )
g '( k x , k y ) = g ( k x , k y ) h(k x , k y )

g '( x) =

g '( x) = g ( x) h( x) g '(k x ) = g (k x )h(k x )

121 242 121
FFT2

g ( x ')h( x x ') dx '

FFT2

=
FFT2

X=imread('happyface.jpg'); imagesc(X); B=[1 2 1; 2 4 2; 1 2 1]; Y=conv2(X,B); Imagesc(Y);colormap(bone);

=
11

3) Image Reconstruction : projection data

Ref: Kalendar 12

10/19/2006

3) Backprojection - Linear single backprojection Detector Source

inverse = 1 back projection Detector Source

13

3) Backprojection two linear projections Source Detector 1

Detector 2

Detector 2

+
Source

Detector 1

14

10/19/2006

3) Backprojection Multiple linear projections 3 projections 4 projections many projections

Original object
15

3) Image Reconstruction concept of backprojection

Ref: Bushberg 16

10/19/2006

3) Image Reconstruction Algebraic Reconstruction Technique (ART)

1100 1010 1001 0110 0101 0011 7 6 5 9 8 7

A B C D

Over-determined non-square matrix K Multiply by KT first Invert square matrix Larger problems must be solved iteratively using standard methods for solving large matrix operation problems.

K x =b KTK x = KTb x = [KTK]-1 KTb

Ref: Bushberg 17

3) Relating to x-y coordinate back-projection

Ref: Hobbie

18

10/19/2006

3) Filtered Backprojection - parallel and cone beams

Optional Read through equations from Chapter 3 in book by Kak and Slaney. Available on the web at: http://rvl4.ecn.purdue.edu/~malcolm/pct/pct-toc.html

Ref: Kak and Slaney

19

3) The Radon Transform Kak and Slaney

Ref: Kak and Slaney

20

10

10/19/2006

3) Fourier Slice Theorem Kak and Slaney

Ref: Kak and Slaney

21

3) Fourier Slice Theorem Kak and Slaney

Ref: Kak and Slaney

22

11

10/19/2006

2) Fourier Slice Theorem what does the DFT image represent ?

Projection profile DFT of projection profile

1-D DFT

Creates line Central slice of Entire DFT image

Integrate intensities along x-direction

Ref: Hobbie

23

2) Fourier Slice Theorem what does the DFT image represent ?

Projection profiles

1-D DFTs

Create lines in central slice of entire DFT image

Integrate intensities along x-direction

Ref: Hobbie

The more angles used, the better the Fourier space image is filled

24

12

10/19/2006

2) Fourier Slice Theorem what does the DFT image represent ?

DFT image represents integration of original projections DFT transformed and summed together.

1-D DFTs at each projection

This is the fast way to create the DFT image from projection data. The more projections taken, the more complete the sampling.
Ref: Hobbie
25

3) Radon Transform sinogram data backprojection

Ref: Brown et al.

26

13

10/19/2006

3) Backprojection of Radon Transform

Ref: Brown et al.

27

3) Backprojection filtering required to bring out high

Ref: Brown et al.

28

14

10/19/2006

3) Image Reconstruction - concept of backprojection and filtered backprojection

Ref: Bushberg 29

3) Reconstruction

Ref: Kalendar

30

15

10/19/2006

3) Parallel Beams and Fan Beams

Ref: Kak and Slaney

Chapter 3 in book by Kak and Slaney. http://rvl4.ecn.purdue.edu/~malcolm/pct/pct-toc.html

31

Filtered Backprojection - Review

32

16

10/19/2006

33

MATLAB implementation of a test phantom

P = phantom('Modified Shepp-Logan',200); imshow(P) Shepp-Logan phantom is a standard head CT simulation. Described in detail in Kak & Slaney, pg 102.

Ref: Matlab Help Guide

34

17

10/19/2006

MATLAB implementation of the Radon transform (i.e. projection data)

R = radon(I,theta) If you omit theta, it defaults to 0:179.
Imaging a square iptsetpref('ImshowAxesVisible','on') I = zeros(100,100); I(25:75,25:75) = 1; theta = 0:180; [RI,xp] = radon(I,theta); imshow(theta,xp,RI,[],'notruesize'), colormap(bone), colorbar

Imaging the Modified Shepp-Logan phantom

P = phantom('Modified Shepp-Logan',200); imshow(P)

[RP,xp] = radon(P,theta); imshow(theta,xp,RP,[],'notruesize'), colormap(bone), colorbar Ref: Matlab Help Guide

35

MATLAB implementation of the Inverse Radon transform

(i.e. image reconstruction from projection data)

I = iradon(P,theta) I = iradon(P,theta,interp,filter,d,n)
interp specifies the type of interpolation to use in the backprojection. The available options are listed in order of increasing accuracy and computational complexity: 'nearest' - nearest neighbor interpolation 'linear' - linear interpolation (default) 'spline' - spline interpolation filter specifies the filter to use for frequency domain filtering. filter is a string that specifies any of the following standard filters: 'Ram-Lak' - The cropped Ram-Lak or ramp filter (default). The frequency response of this filter is | f |. Because this filter is sensitive to noise in the projections, one of the filters listed below may be preferable. These filters multiply the Ram-Lak filter by a window that de-emphasizes high frequencies. 'Shepp-Logan' - The Shepp-Logan filter multiplies the Ram-Lak filter by a sinc function. 'Cosine' - The cosine filter multiplies the Ram-Lak filter by a cosine function. 'Hamming' - The Hamming filter multiplies the Ram-Lak filter by a Hamming window. 'Hann' - The Hann filter multiplies the Ram-Lak filter by a Hann window. d is a scalar in the range (0,1] that modifies the filter by rescaling its frequency axis. The default is 1. If d is less than 1, the filter is compressed to fit into the frequency range [0,d], in normalized frequencies; all frequencies above d are set to 0. n is a scalar that specifies the number of rows and columns in the reconstructed image. If n is not specified, the size is determined from the length of the projections. n = 2*floor(size(P,1)/(2*sqrt(2))) 36

Ref: Matlab Help Guide

18

10/19/2006

MATLAB implementation of the Inverse Radon transform

(i.e. image reconstruction from projection data)

I = iradon(P,theta) I = iradon(P,theta,interp,filter,d,n)

KI = iradon(PI,theta) imshow(theta,xp,KI,[],'notruesize'), colormap(bone), colorbar

KP = iradon(P,theta) imshow(theta,xp,KP,[],'notruesize'), colormap(bone), colorbar

37

Ref: Matlab Help Guide

Effect of projection number

theta = 0:40:170

theta = 0:20:170

theta = 0:10:170

theta = 0:1:170

P = phantom('Modified Shepp-Logan',200);[RP,xp] = radon(P,theta);KP = iradon(RP,theta); imshow(theta,xp,KP,[],'notruesize'), colormap(bone), colorbar

38

Ref: Matlab Help Guide

19

10/19/2006

Effect of noise in the projection data

P = phantom('Modified Shepp-Logan',200); [RP,xp] = radon(P,theta); RN1=RP.*(1+sd*randn(size(RP)));

imshow(theta,xp,RN1,[],'notruesize'), colormap(jet), colorbar; KP1 = iradon(RN1,theta); imshow(KP1), colormap(jet), colorbar

sd=0

sd=0.05

sd=0.10

sd=0.20

sd=0.50
39

Ref: Matlab Help Guide

20