Learning Tomography: Advanced DFR Implementation

The Direct Fourier Reconstruction (DFR) implementation previously published (Part II:  [Link]) works correctly. However, the polar to cartesian conversion is not very efficient because there are missing pixels in the 2D Fourier space. Modifying the source code can easily fixed this problem...


The basic implementation is not perfect because when the 2D Fourier space is filled by the 1D-FT of each sinogram row, some pixels are not filled as shown in Fig. 1 displaying the coverage of the 2D Fourier space.

Fig.1: Fill of the 2D Fourier space. As shown in the red zoomed area (B), missing pixels (displayed in black) are never filled by the values of the 1D-FTs during the polar to cartesian conversion.

Thus, to be sure that all the pixels of the Fourier space are correctly filled, we must modify the filling algorithm. Rather than scanning each value of sinogram and computing the X,Y-coordinates of the 2D-Fourier space, it is better (and more flexible) to invert the process by calculating the polar coordinates of the sinogram from the cartesian coordinates of the 2D-Fourier space (Fig.2).

 << Figure >>

The synopsis written in pseudo-code now becomes...

for x=0 to 2DFT_space.height
  for y=0 to 2DFT_space.width
    {radius,theta} = compute_polar_coordinates_from(x,y)
    value = getInterpolatedPixel(radius,theta,FT_sinogram)
    2DFT_space.setPixel(x,y,value)
  end for
end for

The conversion of the cartesian coordinate (x,y) to polar (ρ,θ) is defined by the following formula:

ρ = √(x²+y²) with ρ ∈ [0; +maxRadius]
θ = atan2(y,x) with θ ∈ ]-π; π]

Finally, we have to check that the polar coordinates (ρ,θ) are comprised in the range [0;width[ and [0,height[ of the sinogram, respectively.

The corresponding lines of code are ...

   var theta = -Math.atan2(yc,xc) * DEG;
   var rad = Math.sqrt(xc * xc + yc * yc);
   if (theta < 0 ) {
     rad = -rad;
     theta=Math.min(theta+180.0,180.0);
   }



Finally, the new JavaScript implementation is ...

+++ IJ JavaScript snippet: FourierRec.js +++

	// Direct Fourier Reconstruction
	// Jean-Christophe Taveau
	// http://crazybiocomputing.blogspot.com
	// 1- Get FT(sinogram) and create the final stack
	var FT_sino = IJ.getImage();
	var FT_sino_p = FT_sino.getProcessor();
	var w=FT_sino.getWidth();
	var nbProj= FT_sino.getHeight();
	var step=180/nbProj;
	var orgx = w/2.0;
	var orgy = orgx;
	// Global variable
	var DEG = 180.0 / Math.PI;
	// Create output 2D reconstruction
	var rec = IJ.createImage("Rec2D DirectFourier", "32-bit Black", w, w, 1);
	var rec_p = rec.getProcessor();
	//
	// M A I N
	//
	IJ.log("Fill 2D Fourier Space with 1D-FT central slice");
	// 4- Cartesian -> Polar: nearest scheme interpolation
	for (var y=0;y<rec_p.getHeight();y++) {
	var yc = y-orgy;
	for (var x=0;x<rec_p.getWidth();x++) {
	var xc = x-orgx;
	// Get polar coordinates (theta and radius)
	var 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;
	if (rad < w/2.0) {
	pix = interpolate(sino_x,sino_y,FT_sino_p);
	rec_p.putPixelValue(x,y,pix);
	}
	}
	}
	IJ.log("Inverse 2D Fourier Transform");
	IJ.run(rec, "Swap Quadrants", "");
	var fht2d = new FHT(rec.getProcessor(),true );
	fht2d.inverseTransform();
	rec.setProcessor(fht2d);
	IJ.run(rec, "Swap Quadrants", "");
	rec.show();
	// F U N C T I O N S
	function interpolate(x,y,iproc) {
	var ix = Math.round(x);
	var iy = Math.round(y);
	return iproc.getPixelValue(ix,iy);
	}

view raw FourierRec.js hosted with ❤ by GitHub

+++ End of IJ JavaScript snippet +++

Thanks to this new implementation, it is now possible to improve the type of interpolation scheme.

To be continued...


<< Previous: Zero padding

Next:Better interpolation scheme >>

Other crazybiocomputing posts

Further readings are available in ...

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