Learning Tomography: ART (Part II)

 After introducing the principle of ART according to Kaczmarz method in the last post [Link], it's time to implement it in ImageJ with JavaScript language...


1- Implementation

In ImageJ, the ART implementation is really simple as the projection computation is done by the Plot Profile function. The algorithm is schematically:

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

There is something new compared to the previous post: the relaxation factor λ (lambda) whose value is comprised between 0.0 and 1.0. This parameter allows to modify the convergence speed. A small value yields a more precise solution with a slow convergence. However, a value close to 1.0 allows a higher speed at the expense of quality

+++ IJ JavaScript snippet: ART_simple.js +++

	// ART: Algebraic Reconstruction Technique
	// 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();
	// 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);
	}
	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.js hosted with ❤ by GitHub

+++ End of JavaScript snippet +++

2-Results

To try the previous script, I use the sinogram of Fig. 1 that you can download [Link] as a TIFF 32-bit image without any compression.

Fig.1: Sinogram used to try the ART script.

Open the sinogram and the ART_simple.js script and run it with Macros > Run Macros ... The image of Fig.2 is then created.

Fig.2: Resulting 2D ART reconstruction

Note: If you want to see this script in action uncomment line #27 // rec.show();

It's also possible to see the result of the reconstruction after each iteration  since they are stored in the stack 'Rec2D_Cycles' as shown in Fig.3.

Fig.3: Results of the reconstruction along the five cycles.

In this kind of iterative techniques, the low frequencies appear in the first cycles and then the high frequencies (details and noise) are represented. In the gallery of Fig. 3, Lena appears blurred (mainly composed of low frequencies) and then, the details (high frequencies) are appearing.

Playing with ART parameters

By default, the ART is calculated for 5 cycles and a relaxation factor  of 0.33. However, you can easily modify the iterations number (line #6: var nbCycles=5;) and the relaxation factor  (line #7: var relax=0.33;)

Other crazybiocomputing posts

Further readings are available in ...

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