Skip to content

Instantly share code, notes, and snippets.

@mutterer
Forked from lacan/curvature_radius.bsh
Created September 10, 2020 12:13
Show Gist options
  • Save mutterer/efee4738eb6eb92ccbf68b05ecad91d4 to your computer and use it in GitHub Desktop.
Save mutterer/efee4738eb6eb92ccbf68b05ecad91d4 to your computer and use it in GitHub Desktop.
Curvature radius simple calculation for ImageJ/fiji #Fiji #Beanshell #Curvature #ImageJ
/*
* 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);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment