"Hard Field" Methods

Hands-on Practical Reconstruction using Matlab

Prof. Bill Lionheart, Dr. Marta Betcke, David Szotten

In the following we will use Matlab with Imaging Toolbox. Whenever you want to learn more about a syntax of a command type

help XXX

or type helpdesk and search for the command in the graphical interface.

You will need some additional M-Files. Download them and save in the working directory.

The following short program generates an image containing 3 points.

N = 101;
I = zeros(N);
I(10,60) = 1;
I(80,90) = 1;
I(ceil(N/2),ceil(N/2)) = 1;

Plot the image and its sinogram (Radon transform).

%Radon transform for angles: 0:179 degree
RI = radon(I,0:179);
figure(1),subplot(2,2,1); imagesc(I); colormap gray; title('Image')
figure(1),subplot(2,2,2); imagesc(RI); colormap gray; title('Sinogram')

Notice, the appearance of the sinogram of the point in the center of the image.

Now, we want to see sinogram of a more complicated object. The following code generates the 2D Shepp-Logan phantom and takes its Radon transform.

I = phantom;
%Radon transform for angles: 0:179 degree
RI = radon(I,0:179);
figure(1),subplot(2,3,1); imagesc(I); title('Shepp Logan phantom')
figure(1),subplot(2,3,2); imagesc(RI); title('Sinogram')


Topic 2: Backprojection, Filtered Backprojection

The following code takes the backprojection of the sinogram. It uses linear interpolation and no filtering

Now, we take a Fourier transform along the direction "s". We plot the real and imaginary part of the complex Fourier transform separately.

%Length of the Fourier transform
NFFT = size(RI,1);
%Fourier transform along the first variable
fftRI = fft(RI,NFFT,1);
figure(2),subplot(2,2,1); imagesc(real(fftshift(fftRI,1))); title('Real part of Fourier(Radon(I))')
figure(2),subplot(2,2,2); imagesc(imag(fftshift(fftRI,1))); title('Imaginary part of Fourier(Radon(I))')

Here for convenience, we used a Matlab command fftshift, which shifts the 0 frequency component to the center of each column.

Next, we compute a discrete Ram-Lak filter (Ramp filter with low pass filter) and plot a sinogram filtered along the "n" direction (in frequency domain)

%Ram-Lak Filter in frequency domain
filt = abs(-floor(NFFT/2):floor(NFFT/2));
%Filtered sinogram in frequency domain
fftfiltRI = fftshift(fftRI,1).*(filt'*ones(1,180));
figure(2),subplot(2,2,3); imagesc(real(fftfiltRI)); title('Real part of Filter(Fourier(Radon(I)))')
figure(2),subplot(2,2,4); imagesc(imag(fftfiltRI)); title('Imaginary part of Filter(Fourier(Radon(I)))')

Now, we compute the inverse Fourier transform and plot the filtered sinogram

%Filtered sinogram in space domain
filtRI = ifft(ifftshift(fftfiltRI,1),NFFT,1);
figure(1),subplot(2,3,4); imagesc(filtRI); title('Filtered sinogram')

The filtered sinogram can now be backprojected without a filter. The same result can be obtained by calling iradon with Ram-Lak filter

%Backprojection of the filtered sinogram
FIR = iradon(filtRI,0:179,'linear','none');
figure(1),subplot(2,3,5); imagesc(FIR); title('Backprojected filtered sinogram')

%Filtered backprojection with iradon with Ram-Lak filter.
%You can use also other windows (see help iradon).
FIR2 = iradon(RI,0:179,'linear','Ram-Lak');
figure(1),subplot(2,3,6); imagesc(FIR2); title('Filtered backprojection (iradon)')

As a further exercise generate the backprojection and Ram-Lak filtered backprojection for the 3-point phantom using iradon.

Topic 3: Algebraic reconstruction (ART)

Construct a small version of a Shepp-Logan phantom

I = phantom(30);

We take forward projection with Siddon algorithm. You will need an additional file ray_tracer_2d.m

[rays, sino] = siddon(I);

Remove all the rays, that did not hit the image and prepare the right hand side b.

ind = find(sum(rays,2));
A = rays(ind,:);
b = sino(:); b = b(ind);

Reconstruct with the Landweber algorithm

[x,rs] = landweber(A,b,500);

Reconstruct with the Kaczmarz algorithm

[x,rs] = kaczmarz(A,b,25,1:size(A,1));

See what happens if Kaczmarz is applied with a random order of projections

[x,rs] = kaczmarz(A,b,15,randperm(size(A,1)));

Reconstruct with block ART, with 30 ordered sets (OS), each set contains one projection per angle.

[x,rs] = sart(A,b, 50, 30, 1:30:size(A,1)-30)

Topic 4: Noise

You can carry out the reconstruction in points 2 and 3 in the presence of noise. The following code generates a sinogram with 5% normally distributed noise.

randn('state',0);
sinon = sino + mean(sino(:))*0.05*var(sino(:))*randn(size(sino));
imagesc(sinon);
bn = sinon(:);
bn = bn(ind);

Topic 5: Limited data

We now investigate how the reconstruction is affected by data truncation. The examples below deal with truncation resulting from scanning an object larger than the detector. Other common truncation modes include limited angle, where for certain angles no data is available. If you have time you can try this yourself.

We begin by generating a phantom again, along with a complete sinogram:

%define a phantom
I = phantom;

%calculate its radon transform
RI = radon(I,0:179);

%reconstruct image
RI_fbp = iradon(RI, 0:179,'linear','Ram-Lak');

We now simulate a smaller detector size by setting data outside the detector to zero. We attempt to reconstruct using the same algorithm and compare the results:

%truncate sinogram
RI_limited = RI; RI_limited( [1:109, 259:end], :) = 0;

%reconstruct using truncated sinogram
RI_limited_fbp = iradon(RI_limited, 0:179,'linear','Ram-Lak');

figure(2), clf
subplot(2,2,1), imagesc(RI), colormap gray, title('Original sinogram')
subplot(2,2,2), imagesc(RI_limited), title('Truncated sinogram')
subplot(2,2,3), imagesc(RI_fbp), axis image, title('Reconstructed from original sinogram')
subplot(2,2,4), imagesc(RI_limited_fbp), axis image, title('Reconstructed from truncated sinogram')

To get a feel for how this works you can try changing some of the parameters, such as the level of truncation. You can download a Matlab applet, which lets you do this interactively. You need the two files listed below:

limited_data_app.m
limited_data_app.fig

To run, simply type
limited_data_app

This uses the standard phantom. To use a custom image, give it as an argument:
square = zeros(100); square(10:20,30:40) = 1; limited_data_app(square);