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);