Learning Tomography: DFR with Fourier - Part II

Second (and final) part of the implementation of the Direct Fourier Reconstruction with Fourier  − instead of Hartley − Transform...

1- Script:  Computing 2D reconstruction image

Fig. 1: Scheme of the various steps for a Direct Fourier Reconstruction. Part I: From the sinogram, 1D Fourier Transform (FT) is applied to each row. Part II: Each 1D FT is inserted as a central line into a 2D Fourier space according to its rotation angle. Finally, this 3D Fourier space is converted into an image thanks to an inverse FT.

As summarized in Fig. 1, the script takes as input the stack of 1D FTs  − calculated previously from the sinogram [Link] − and then fills the 2D Fourier Space depending of interpolation scheme − Nearest or Bilinear whose implementation is described in this post [Link]− .
For sake of convenience, two additional methods getComplex(...) and setComplex(...) are added to read/write the real and imaginary parts in the 2D Fourier stack.
The output image is the 2D Fourier stack.

+++ JavaScript IJ snippet: FourierRec_FFT.js +++

	// Direct Fourier Reconstruction
	// Use of Fourier Transform
	// Jean-Christophe Taveau
	// http://crazybiocomputing.blogspot.com
	// 1- Get sinogram and create the final stack
	var FTsino = IJ.getImage();
	var FTsino_p = FTsino.getStack()
	var w = FTsino.getWidth();
	var nbProj = FTsino.getHeight();
	var step = 180/nbProj;
	var orgx = w/2;
	var orgy = w/2;
	// Global variable
	var DEG = 180.0 / Math.PI;
	var NEAREST = 1;
	var BILINEAR = 2;
	var mode = NEAREST;
	// Output 2D reconstruction
	var rec = IJ.createImage("Complex of Rec2Direct", "32-bit Black", w, w, 2);
	rec.getStack().setSliceLabel("Real",1);
	rec.getStack().setSliceLabel("Imaginary",2);
	var rec_p = rec.getStack();
	//
	// M A I N
	//
	IJ.log("Part II: Fill 2D Fourier Space with 1D-DFT central slice");
	var theta = 0.0;
	var pix={};
	for (var y=0;y<rec_p.getHeight();y++) {
	if (y%10==0)
	IJ.log("Row # " + y);
	var yc = y - orgy;
	for (var x=0;x<rec_p.getWidth();x++) {
	var xc = x - orgx;
	theta = Math.atan2(yc,xc) * DEG;
	var rad = Math.sqrt(xc * xc + yc * yc);
	if (theta <0.0) {
	rad = -rad;
	theta=Math.min(theta+180.0,180.0)
	}
	sino_x = rad + orgx;
	sino_y = theta;
	pix = interpolate(sino_x,sino_y,FTsino_p);
	setComplex(rec_p,x,y,pix);
	}
	}
	rec.show();
	// F U N C T I O N S
	function interpolate(x,y,iproc) {
	switch (mode) {
	case NEAREST : return nearest(x,y,iproc);break;
	case BILINEAR: return bilinear(x,y,iproc);break;
	default : return nearest(x,y,iproc);break;
	}
	}
	function nearest(x,y,iproc) {
	var ix = Math.round(x);
	var iy = Math.round(y);
	return getComplex(iproc,ix,iy);
	}
	function bilinear(x,y,iproc) {
	var ix = Math.floor(x);
	var iy = Math.floor(y);
	var dx = x - ix;
	var dy = y - iy;
	v00 = getComplex(iproc,ix,iy);
	v10 = getComplex(iproc,ix+1,iy );
	v01 = getComplex(iproc,ix,iy+1);
	v11 = getComplex(iproc,ix+1,iy+1);
	var value={};
	value.re = v00.re*(1-dx)*(1-dy) + v10.re*dx*(1-dy) + v01.re*(1-dx)*dy + v11.re*dx*dy;
	value.im = v00.im*(1-dx)*(1-dy) + v10.im*dx*(1-dy) + v01.im*(1-dx)*dy + v11.im*dx*dy;
	return value;
	}
	function getComplex(ip,x,y) {
	var re = ip.getVoxel(x,y,0);
	var im = ip.getVoxel(x,y,1);
	return {"re":re,"im":im};
	}
	function setComplex(ip,x,y,cmplx) {
	ip.setVoxel(x,y,0,cmplx.re);
	ip.setVoxel(x,y,1,cmplx.im);
	}

view raw FourierRec_FFT.js hosted with ❤ by GitHub

+++ End of JavaScript IJ snippet +++

2- Results

To get the final image, click on the stack entitled "Complex of RecDirect" and process an inverse FFT (Process > FFT > Inverse FFT), to get the image of Fig. 2. Finally, swap the quadrants (Process > FFT > Swap Quadrants).

Fig.2: Image obtained after computing an Inverse Fourier Transform of the image created by the script FourierRec_FFT.j. In this case, the sinogram was previously zero-padded.

<<: Part I - Fourier Transform of sinogram

Other crazybiocomputing posts

Further readings are available in ...

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