Learning Tomography: WBP 3D reconstruction

In this series Learning Tomography, I always work with 2D slice ... Here is an example of a 3D reconstruction calculated with the weighted (or filtered) back-projection algorithm. The input data correspond to the stack of sinograms (previously calculated here [Link]).

1- Modifying and improving the script

Compared to the previous scripts I published, I just add another loop to scan all the slices of the sinogram stack. Here is a synopsis of the script in pseudo-code:

  3D_rec = new stack
  for each sinogram
    2D_rec = new slice
    for each row of this sinogram
      bp = back-project(row)
      bp.rotate
      2D_rec = 2D_rec + bp
    end for
    3D_rec.add_slice(2D_rec)
  end for
   
Moreover, I tried to optimize the JavaScript code to speed up the process of 2D reconstruction. This new version − detailed in this post [Link] − is located in the function calcRec2D(one_sinogram, proj_num, angles_array).

Here is the script.

+++ IJ JavaScript snippet: WBP_3Drec.js +++

	// 3D WBP: 3D Weighted BackProjection
	// Jean-Christophe Taveau
	// 22 oct 2012
	// http://crazybiocomputing.blogspot.com
	// Get sinograms stack and information
	var imp = IJ.getImage();
	var w=imp.getWidth();
	var nbProj=imp.getHeight();
	var step=180/nbProj;
	var nz=imp.getNSlices();
	// Output volume
	var stack = IJ.createImage("Rec3D", "32-bit Black", w, w, 2);
	// Global variables
	var tmp = IJ.createImage("tmp", "32-bit Black", w, w, 1);
	var sinogram = IJ.createImage("sino", "32-bit Black", w, w, 1);
	var row = new FloatProcessor(w,1);
	var ic = new ImageCalculator();
	// Precompute angles
	var angles=[];
	for (var i=0;i<nbProj;i++)
	angles[i]=-i*step;
	//
	// M A I N
	//
	for (var j=0;j<nz;j++) {
	// Sinogram -> 2D Slice in the final volume
	var a_sinogram = imp.getStack().getProcessor(j+1);
	if (a_sinogram.getStatistics().mean != 0.0 ) {
	var slice = calcRec2D(a_sinogram,nbProj,step);
	stack.getStack().addSlice(slice.getProcessor() );
	IJ.log("Slice # " + (j+1) );
	}
	else
	IJ.log("Slice # " + j + " no computation");
	}
	stack.show();
	// F U N C T I O N S
	function calcRec2D(sino_ip,num,angle) {
	var w=sino_ip.getWidth();
	var bp = null;
	var rec = IJ.createImage("Rec2D", "32-bit Black", w, w, 1);
	for (var i in angles) {
	// 1- get ith row
	for (var k=0;k<w;k++)
	row.putPixelValue(k,0,sino_ip.getPixelValue(k,i));
	// 2- back-project
	bp = row.resize(w,w);
	// 3- rotate accordingly
	bp.rotate(angles[i]);
	// 4- add to 2D rec slice
	tmp.setProcessor(bp);
	ic.run("Add",rec,tmp);
	}
	return rec;
	}

view raw WBP_3Drec.js hosted with ❤ by GitHub

+++ End of IJ JavaScript snippet: WBP_3Drec.js +++

2- Computing the 3D reconstruction

1. Copy and paste the script in a JavaScript window (Plugins > New > JavaScript), click on the sinogram stack (available as a 32-bit TIFF ramp-filtered stack [Link] or web page [Link]) and run the script (in JavaScript window, Macros > Run Macro). 
2. Once executed, choose the middle slice of the Rec3D window, draw a small circle in the center, and normalize all the slices of the volume at 0% of saturated pixels (Process > Enhance Contrast... > Normalize ).

Note: The stack was previously filtered with a truncated ramp filter with the following formula [if d < 100.0) v=d;]. You can try with more sophisticated filters.

Here is a surface representation calculated with Volume Viewer plugin of the final volume 3DRec.

<< Previous: Computing the sinograms

Other crazybiocomputing posts

Further readings are available in ...

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