Learning Tomography: Adding Constraints

The computation time for iterative reconstruction techniques (IRT) is typically orders of magnitude longer than for Weighted (Filtered) Back-Projection algorithm [Link], however, IRT are powerful  because during the iteration process, you are able to add constraints to enhance convergence.


1- Non-negativity and box constraints

In most cases, having negative pixel values are non-sense during ART process, thus, it is easy to constrain the ART algorithm to clamp all the negative pixel values to zero. An extension of the latter is to restrain values between a minimum and a maximum.

In our script ART_simple.js, we only add a sixth step within the inner loop.

     //6- Add constraints
    IJ.run(rec, "Macro...", codeMin);
    IJ.run(rec, "Macro...", codeMax);


where codeMin and codeMax are respectively defined as

  var codeMin="code=[v=maxOf(v,"+min+")];";
  var codeMax="code=[v=minOf(v,"+max+")];";

Other parameters

Depending of the prior knowledge of the sample, you can add what you want to speed up the convergence or yield a better solution if you have a small number of projections, for example.
Here, we know that the sample is surrounded by a mask of radius=120.0. This can be injected as a new constraint.

   IJ.run(rec, "Macro...", "code=[if (d > 120.0) v=0.0;]");

2- Script and example

Here is the updated script ART_simple_constraints.js

+++ IJ JavaScript snippet: ART_simple_constraints.js +++

	// ART: Algebraic Reconstruction Technique
	// Jean-Christophe Taveau
	// 5 nov 2012
	// http://crazybiocomputing.blogspot.com
	// ART Parameters
	var nbCycles=5;
	var relax=0.33;
	var min=0.0;
	var max=2.0;
	var codeMin="code=[v=maxOf(v,"+min+")];";
	var codeMax="code=[v=minOf(v,"+max+")];";
	// 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();
	// 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("Subtract",proj,a);
	IJ.run(proj, "Multiply...", "value="+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("Add",rec,proj);
	//6- Add constraints
	IJ.run(rec, "Macro...", codeMin);
	IJ.run(rec, "Macro...", codeMax);
	}
	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 ART_simple_constraints.js hosted with ❤ by GitHub

+++ End of IJ snippet +++

From our reference sinogram [Link], this new script calculates Lena's portrait as shown in Fig.1.

Fig.1: Reconstruction calculated with 3 cycles, a relaxation factor of 0.33 and constraints (min=0.0, max=2.0, and radius=120.0).

Other crazybiocomputing posts

Further readings are available in ...

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