Learning Tomography: DFR with Fourier Transform

In this series dedicated to Direct Fourier Reconstruction (DFR), here is a new version based on Discrete Fourier Transform rather than the Hartley Transform...


Compared to the Hartley Transform, the Fourier transform (FT) uses complex numbers [Wikipedia].

1 - Complex numbers

1-1- JavaScript and complex numbers

Unfortunately in JavaScript, there is no class Complex, but you can easily create one with the following constructor:

function Complex(real,imaginary) {
  this.re = real; // Real part
  this.im = imaginary; // Imaginary part
}

but, as we don't need many functionalities, a simple associative array - in this case - is enough and a complex number can be defined as follows:

var complex = {"re":real,"im":imaginary};

1-2- ImageJ and FFT

In ImageJ, a FT can be stored in a 2-slices 32-bit stack: a first slice entitled 'Real' and a second slice 'Imaginary'. Moreover, if this stack contains a title beginning with 'Complex of ', it is automatically recognized as a frequency stack by ImageJ and an inverse FFT can be applied.

var FT = IJ.createImage("Complex of A", "32-bit Black", w, w, 2);
FT.getStack().setSliceLabel("Real",1);
FT.getStack().setSliceLabel("Imaginary",2);

Note: Functions getVoxel(...) and setVoxel(...) can be used to access the pixel values within the stack.

2- Script: Computing FT(sinogram)

The algorithm is exactly the same than previously [Link], the only differences are:

1. The management of a 2-slices stack (real and imaginary).
2. The function DHT is replaced by DFT (Discrete Fourier Transform)

+++ JavaScript IJ snippet: FourierRec_FFTsino.js +++

	// Direct Fourier Reconstruction
	// Jean-Christophe Taveau
	// http://crazybiocomputing.blogspot.com
	// 1- Get sinogram and create the final stack
	var imp = IJ.getImage();
	var ip=imp.getProcessor();
	var w=imp.getWidth();
	var nbProj=imp.getHeight();
	var step=180/nbProj;
	var orgx = w/2;
	var orgy = w/2;
	// Output complex image of sinogram
	var FTsino = IJ.createImage("Complex of sinogram", "32-bit Black", w, nbProj, 2);
	FTsino.getStack().setSliceLabel("Real",1);
	FTsino.getStack().setSliceLabel("Imaginary",2);
	var FTsino_p = FTsino.getStack();
	//
	// M A I N
	//
	IJ.log("Part I: 1D-FFT(sinogram) i.e. stack of 1D-FFT(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- Init complex numbers of this row
	var tmp = new Array(row.length);
	for (var j in row)
	tmp[j] = {"re":row[j], "im": 0.0};
	// 3- 1D-DFT
	row = DFT(tmp,false);
	// 4- Swap data
	for (var j=0;j<w;j++)
	{
	var r = (j + w - w/2) % w;
	if (r % 2 == 0) {
	FTsino_p.setVoxel(r,y,0,row[j].re);
	FTsino_p.setVoxel(r,y,1,row[j].im);
	}
	else {
	FTsino_p.setVoxel(r,y,0,-row[j].re);
	FTsino_p.setVoxel(r,y,1,-row[j].im);
	}
	}
	}
	FTsino.show();
	// F U N C T I O N S
	// Slow Fourier Transform
	// From http://arachnoid.com/signal_processing/dft.html
	function DFT(input,inverse)
	{
	var N = input.length;
	var dft = new Array(N);
	var pi2 = (inverse)? 2.0 * Math.PI:-2.0 * Math.PI;
	var a=0.0;
	var invs = 1.0 / N;
	for(var y = 0;y < N;y++) {
	dft[y] = {"re":0.0, "im":0.0};
	for(var x = 0;x < N;x++) {
	a = pi2 * y * x * invs;
	dft[y].re += input[x].re * Math.cos(a) - input[x].im * Math.sin(a);
	dft[y].im += input[x].re * Math.sin(a) + input[x].im * Math.cos(a);
	}
	if(!inverse) {
	dft[y].re *= invs;
	dft[y].im *= invs;
	}
	}
	return dft;
	}

view raw FourierRec_FFTsino.js hosted with ❤ by GitHub

+++ End of JavaScript IJ snippet +++

>> Next: Part II - Fourier reconstruction

Other crazybiocomputing posts

Further readings are available in ...

• Learning Tomography Series  [Link]
• Image Processing TOC [Link]