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 index^th row of the sinogram.

Here is the JavaScript implementation

+++ JavaScript IJ snippet +++

	// Direct Fourier Reconstruction
	// Jean-Christophe Taveau
	// http://crazybiocomputing.blogspot.com
	// 1- Get sinogram
	var imp = IJ.getImage();
	var ip=imp.getProcessor();
	var w=imp.getWidth();
	var nbProj= imp.getHeight();
	// Global variable
	var DEG = 180.0 / Math.PI;
	// Create ouput
	var FT_sino = IJ.createImage("1D-DHT", "32-bit Black", w, nbProj, 1);
	var FT_sino_p = FT_sino.getProcessor();
	//
	// M A I N
	//
	IJ.log("Sinogram -> 1DFHT(row)");
	for (var y=0;y<nbProj;y++)
	{
	IJ.log("Slice #"+y);
	// 1- Extract a 1D-projection (one row of sinogram)
	var row = getRow(ip,y);
	// 2- 1D-DHT
	row = DHT(row);
	// 3- Swap data
	for (var j=0;j<w;j++)
	{
	var r = (j + w - w/2) % w;
	if (r%2==0)
	FT_sino_p.putPixelValue(r,y,row[j]);
	else
	FT_sino_p.putPixelValue(r,y,-row[j]);
	}
	}
	FT_sino.show();
	// F U N C T I O N S
	function getRow(ip,iRow) {
	var data= [];
	for (var x=0;x<ip.getWidth();x++)
	data[x]= ip.getPixelValue(x,iRow);
	return data;
	}
	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;
	}

view raw fourierRecSinogram.js hosted with ❤ by GitHub

+++ 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]