Thursday, January 3, 2013

Learning Tomography: Fourier Rec. - Part I



First part of this series dedicated to the implementation of Direct Fourier Reconstruction technique: the computation of the Fourier transform of the sinogram.


In this first part, we are going to transform the sinogram (a stack of 1D projections) into a stack of 1D Fourier transforms as shown in Fig. 1.
Fig.1: Main steps required to compute the Fourier Transform of the sinogram.

1- 1D Fourier Transform

First, we need a 1D Fourier Transform and I choose a (slow) basic implementation of the 1D Discrete Hartley transform (a Fourier-related transform) [Wikipedia] used by ImageJ.

function DHT(input) {
  var N = input.length;
  var dht = new Array(N);
  var w=0.0;
  for (var k = 0; k < N; k++) {
    dht[k] = 0;
    for (var n = 0; n < N; n++) {
      w = 2 * Math.PI * n / N;
      dht[k] += input[n] * (Math.cos(k * w) + Math.sin(k * w));
    }
  }
  return dht;
}

Note: The main advantage of the Hartley transform is its simplicity and the fact that there is no complex number.

2- Implementation

The script calculates for each row of the sinogram, the corresponding 1D Discrete Hartley Transform (DHT) and then swap the two halves of the DHT data to put the origin at the middle of the row. The pseudo-code is as follows:

for index = 0 to sinogram.height
   data = getRow(index)
   dht = compute 1D-DHT(data)
   output[index] = swap_two_halves(dht)
end for


In the script, a function getRow(processor, index) is added to conveniently extract the indexth row of the sinogram.

Here is the JavaScript implementation

+++ JavaScript IJ snippet +++ +++ End of IJ snippet +++

3- Results

From our 32-bit test sinogram of Lena [Link], running the script yields the following image (Fig. 2).

Fig. 2: Hartley transform of sinogram. For display, a logarithm was applied (Process > Math > Log) and a normalize at 0% (Process > Enhance Contrast)


<< Previous: Introduction

>> Next: Fill the 2D-Fourier space

Other crazybiocomputing posts

Further readings are available in ...
  • Learning Tomography Series  [Link]
  • Image Processing TOC [Link]

No comments:

Post a Comment