Learning Tomography: Multiplicative ART

There are several variants of the Algebraic Reconstruction Techniques (ART): One of them is Multiplicative ART (MART) based on the same principle but the error is calculated as a ratio.

1- Algorithm

In MART, the algorithm becomes.

for each iteration
  for each row of sinogram
    p ← row;
    a ← calc_projection(rec, angle);
    err ← (p / a)^λ;
    rec  ← rec * err;
  end for
end for

2- Implementation

Nothing special about the implementation, this is exactly the same script as ART_simple.js but step #4 replaced by ...

    // 4- Compute the error as a ratio
  ic.run("Divide",proj,a); 
  IJ.run(proj, "Macro...", "code=[v=pow(v,"+relax+");]");

And the reconstruction update is replaced by a multiplication:

  ic.run("Multiply",rec,proj);

Moreover, to avoid a division by zero at the first iteration, the 2D reconstruction image is initialized with the mean value of the sinogram.

  var stats = imp.getStatistics();
  var avg=stats.mean;
  rec.getProcessor().setValue(avg);
  rec.getProcessor().fill();

Of course, you can add the constraints as described in a previous post [Link].

+++ IJ JavaScript snippet: MART_simple.js +++

	// MART: Multiplicative ART
	// Jean-Christophe Taveau
	// 5 nov 2012
	// http://crazybiocomputing.blogspot.com
	// ART Parameters
	var nbCycles=5;
	var relax=0.33;
	// 1- Get sinogram information
	var imp = IJ.getImage();
	var w=imp.getWidth();
	var nbProj=imp.getHeight();
	var step=180/nbProj;
	// Create output and temporary images/stack
	var rec = IJ.createImage("Rec2D", "32-bit Black", w, w, 1);
	var recStack = IJ.createImage("Rec2D_Cycles", "32-bit Black", w, w, nbCycles);
	var a = IJ.createImage("tmp", "32-bit Black", w, 1, 1);
	var ic = new ImageCalculator();
	var stats = imp.getStatistics();
	var avg=stats.mean;
	rec.getProcessor().setValue(avg);
	rec.getProcessor().fill();
	// M a i n l o o p
	// Comment out the following lines to see the ART in action
	// rec.show();
	// recStack.show();
	for (var i=0;i<nbCycles;i++)
	{
	for (var j=0;j<nbProj;j++)
	{
	// 2- Compute projection from temporary 2D-rec
	calcProj(rec,step*j);
	// 3- Extract a 1D-projection (one row of sinogram)
	imp.setRoi(0,j,imp.getWidth(), 1);
	var proj = imp.duplicate();
	// 4- Compute the error between sinogram and calculated projections
	ic.run("Divide",proj,a);
	IJ.run(proj, "Macro...", "code=[v=pow(v,"+relax+");]");
	// 5- Create the backprojected ray from this err pixel and add it to rec;
	IJ.run(proj, "Size...", "width="+w+" height="+w+" average interpolation=None");
	IJ.run(proj, "Rotate... ", "angle="+(-j * step)+" grid=1 interpolation=Bilinear fill");
	ic.run("Multiply",rec,proj);
	}
	recStack.setSlice(i+1);
	recStack.setProcessor(rec.getProcessor().duplicate());
	}
	rec.show();
	recStack.show();
	// F U N C T I O N S
	function calcProj(img, rotangle)
	{
	var tmp = img.duplicate();
	tmp.getProcessor().setBackgroundValue(0.0);
	tmp.getProcessor().rotate(rotangle);
	tmp.setRoi(0,0,tmp.getWidth(), tmp.getHeight());
	var plot = new ProfilePlot(tmp);
	var data = plot.getProfile();
	for (var x in data)
	a.getProcessor().putPixelValue(x,0,data[x]);
	tmp.close();
	}

view raw MART_simple.js hosted with ❤ by GitHub

+++ End of IJ snippet +++

3- Some results

From this sinogram [Link], the MART script yields the 2D reconstruction of Fig. 1.

Fig. 1: Five cycles of MART with a relaxation factor of 0.33.

Other crazybiocomputing posts

Further readings are available in ...

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