Skip to content

Instantly share code, notes, and snippets.

@lacan
Last active October 13, 2022 15:18
Show Gist options
  • Save lacan/42f4abe856f697e664d1062c200fd21f to your computer and use it in GitHub Desktop.
Save lacan/42f4abe856f697e664d1062c200fd21f 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);
@lacan
Copy link
Author

lacan commented Feb 13, 2018

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.

@cyril-gdj
Copy link

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment