Curvature radius simple calculation for ImageJ/fiji #Fiji #Beanshell #Curvature #ImageJ

curvature_radius.bsh

/*
 * Simple script that takes a selection and computes the curvature radius for the whole perimeter of the shape.
 * 
 * By Olivier Burri, 
 * BioImaging and Optics Platform, BIOP
 * Ecole Polytechnique Fédérale de Lausanne (EPFL)
 * Last update: May 2017
 * 
 * 
 */
import ij.ImagePlus;
import ij.IJ;
import ij.process.FloatPolygon;
import ij.gui.Roi;
import ij.gui.Plot;
import ij.measure.ResultsTable;
import ij.gui.OvalRoi;
import ij.plugin.frame.RoiManager;
// Discrete curvature

//Get ROI and then get all coordinates
ImagePlus imp 	= IJ.getImage();	//	get the active image
IJ.run(imp, "Interpolate", "interval=1 smooth");
rm = RoiManager.getInstance();
if (rm==null) {
	rm = new RoiManager();
}
	rm.reset();

FloatPolygon r = imp.getRoi().getFloatPolygon();
rm.addRoi(imp.getRoi());

double[] k = new double[r.npoints];
double[] x = new double[r.npoints];
rt = new ResultsTable();
int step = 30;
double tol = 25000;
	double px = 0;
	double py = 0;
	double pr = 0;
	double lx = 0;
	double ly = 0;
	double lr = 0;
	int n =0;
	int ln=0;
	
for (int i=0; i<r.npoints; i++) {
	int p0 = (i<step) ? r.npoints-step+i : i-step;
	int p1 = i;
	int p2 = (i+step)%r.npoints;
	double [] centrad = new double [3];

	double[] pa = {r.xpoints[p0], r.ypoints[p0]};
	double[] pb = {r.xpoints[p1], r.ypoints[p1]};
	double[] pc = {r.xpoints[p2], r.ypoints[p2]};


	// returns 3 double values: the centre (x,y) coordinates & radius
	// of the circle passing through 3 points pa, pb and pc
		double a, b, c, d, e, f, g;

		if ((pa[0]==pb[0] && pb[0]==pc[0]) || (pa[1]==pb[1] && pb[1]==pc[1])){ //colinear coordinates
			centrad[0]=0; //x
			centrad[1]=0; //y
			centrad[2]=-1; //radius
			return;
		}

		a = pb[0] - pa[0];
		b = pb[1] - pa[1];
		c = pc[0] - pa[0];
		d = pc[1] - pa[1];

		e = a*(pa[0] + pb[0]) + b*(pa[1] + pb[1]);
		f = c*(pa[0] + pc[0]) + d*(pa[1] + pc[1]);

		g = 2.0*(a*(pc[1] - pb[1])-b*(pc[0] - pb[0]));
		//  If g is 0 then the three points are colinear and no finite-radius
		//  circle through them exists. Return -1 for the radius. Somehow this does not
		// work as it should (representation of double number?), so it is trapped earlier..

		if (g==0.0){
			centrad[0]=0; //x
			centrad[1]=0; //y
			centrad[2]=-1; //radius
		}
		else { //return centre and radius of the circle
			centrad[0] = (d * e - b * f) / g;
			centrad[1] = (a * f - c * e) / g;
			centrad[2] = Math.sqrt(Math.pow((pa[0] - centrad[0]),2) + Math.pow((pa[1] - centrad[1]),2));
		}


		// Get similarity to previous measurement
		double sim = Math.sqrt(Math.pow(px - centrad[0],2) + Math.pow(py - centrad[1],2) + Math.pow(pr + centrad[2], 4));
		
		if(sim <= tol && centrad[2]!=-1) {
			px = (px + centrad[0]) / 2;
			py = (py + centrad[1]) / 2;
			pr = (pr + centrad[2]) / 2;
			n++;
		} else {
			n=0; 
			px = centrad[0];
			py = centrad[1];
			pr = centrad[2];
		}
		// Keep largest
		if(n > ln) {
			lx = px;
			ly = py;
			lr = pr;
			ln=n;
			IJ.log("Largest r = "+lr+" x: "+lx+" y: "+ly+" n="+ln);
		}

		// Get the sign of the radius, it just depends if we are above or below the circle
		double sign = Math.signum (pb[1] - centrad[1]);
	//k[i] = 2*Math.abs((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)) / Math.sqrt(((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1))*((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1))*((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)));
	//x[i] = i;
	rt.incrementCounter();
	rt.addValue("CenterX", centrad[0]);
	rt.addValue("CenterY", centrad[1]);
	rt.addValue("R", centrad[2] );
	rt.addValue("R Signed", centrad[2] * sign );
	if (centrad[2] < 1e5 && centrad[2] != -1 ) {
		Roi r = new OvalRoi(centrad[0] - centrad[2], centrad[1] - centrad[2], 2*centrad[2], 2*centrad[2]);
		r.setName("Circle at Index "+IJ.pad(rt.getCounter(),5));
		rm.addRoi(r);
	}
	
}
Roi l = new OvalRoi(lx - lr, ly - lr, 2*lr, 2*lr);
l.setName("Most Constant Circle");
rm.addRoi(l);
l.setImage(imp);
rt.show("SS");
rm.select(imp, rm.getCount()-1);

a bit like this

I apologize, as I am a new to the script using. Anyone who could guide me, how to use the above script for imageJ. Should I make a txt file and then install plugin or some other way. I have tried but it shows errors. Please help me out.

@teodej

Hi, Thanks for this great script. Do you know a way to modify it to get the sign of the curvature knowing that it turns anticlockwise and using a normal line to the radius?

Sorry for the very long response time. Gist comments do not show up on my warnings. I think I managed to implement this, simply by checking the sign of the vector defined by the circle radius and the point in question. Let me know if it works.

Dear @metabbas
The simplest is to use the ImageJ forum when you have this sort of question and to familiarize yourself with the ImageJ documentation
It is not a plugin but a script. Just download the file or copy it to a text file and change its extension to .bsh (Beanshell). You can then just drag and drop this file into ImageJ. It should open the script editor. Then you can simply hit the Run button.

Please note that the plugin expects a Region of Interest (ROI) to be drawn on the image. Freehand or area makes the most sense for curvature measurements.

Dear @lacan
I have a question about the script provided. This one works well I have each of the curves which are represented on my image with the data in ROI manager. My question is what are the results to use to calculate a total or average curvature attributed to a form or how to represent it in fiji? or in your script this may already be provided and I did not pay attention
thank you in advance