You are here: MIMS > events > workshops > VCIPT-2007
MIMS

# 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.

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

```IR = iradon(RI,0:179,'linear','none'); figure(1),subplot(2,3,3); imagesc(IR); title('Backprojected sinogram')```

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:

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);`