io.github.mianalysis.mia.process.analysis.CurvatureCalculator Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of mia-algorithms Show documentation
Show all versions of mia-algorithms Show documentation
ModularImageAnalysis (MIA) is an ImageJ plugin which provides a modular framework for assembling image and object analysis workflows. Detected objects can be transformed, filtered, measured and related. Analysis workflows are batch-enabled by default, allowing easy processing of high-content datasets.
package io.github.mianalysis.mia.process.analysis;
import java.awt.Color;
import java.util.ArrayList;
import java.util.Iterator;
import java.util.TreeMap;
import org.apache.commons.math3.analysis.interpolation.LoessInterpolator;
import org.apache.commons.math3.analysis.interpolation.SplineInterpolator;
import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
import ij.ImagePlus;
import ij.gui.Line;
import ij.gui.Overlay;
import io.github.mianalysis.mia.object.coordinates.Point;
/**
* Created by sc13967 on 26/01/2018.
*/
public class CurvatureCalculator {
private ArrayList> path;
private PolynomialSplineFunction[] splines = null;
private TreeMap curvature = null;
private FittingMethod fittingMethod = FittingMethod.STANDARD;
private boolean isLoop = false;
public enum FittingMethod {
STANDARD, LOESS
}
private int NNeighbours = 5;
private double loessBandwidth = 0.04;
private int loessIterations = 0;
private double loessAccuracy = 100;
public CurvatureCalculator(ArrayList> path, boolean isLoop) {
this.path = path;
this.isLoop = isLoop;
}
public void calculateCurvature() {
// Checking there are enough points for fitting
int nKnots = path.size();
if (nKnots < (NNeighbours + 1))
return;
if (isLoop)
nKnots = nKnots + 2 * NNeighbours;
// Preparing line coordinates for spline fitting
double[] t = new double[nKnots];
double[] x = new double[nKnots];
double[] y = new double[nKnots];
int count = 0;
Point prevVertex = null;
// If a loop, starting with the final few points from the opposite end
if (isLoop) {
int len = path.size();
for (int i = NNeighbours; i > 0; i--) {
Point vertex = path.get((len - i));
t[count] = count == 0 ? 0 : t[count - 1] + vertex.calculateDistanceToPoint(prevVertex);
x[count] = vertex.getX();
y[count++] = vertex.getY();
prevVertex = vertex;
}
}
for (Point vertex : path) {
t[count] = count == 0 ? 0 : t[count - 1] + vertex.calculateDistanceToPoint(prevVertex);
x[count] = vertex.getX();
y[count++] = vertex.getY();
prevVertex = vertex;
}
// If a loop, ending with the first few points from the opposite end
if (isLoop) {
for (int i = 0; i < NNeighbours; i++) {
Point vertex = path.get(i);
t[count] = count == 0 ? 0 : t[count - 1] + vertex.calculateDistanceToPoint(prevVertex);
x[count] = vertex.getX();
y[count++] = vertex.getY();
prevVertex = vertex;
}
}
// Storing splines
splines = new PolynomialSplineFunction[6];
// Fitting the spline pair
switch (fittingMethod) {
case LOESS:
LoessInterpolator loessInterpolator = new LoessInterpolator(loessBandwidth, loessIterations,
loessAccuracy);
splines[0] = loessInterpolator.interpolate(t, x);
splines[1] = loessInterpolator.interpolate(t, y);
break;
case STANDARD:
SplineInterpolator splineInterpolator = new SplineInterpolator();
splines[0] = splineInterpolator.interpolate(t, x);
splines[1] = splineInterpolator.interpolate(t, y);
break;
}
// Calculating curvature
splines[2] = splines[0].polynomialSplineDerivative(); // dx
splines[3] = splines[2].polynomialSplineDerivative(); // ddx
splines[4] = splines[1].polynomialSplineDerivative(); // dy
splines[5] = splines[4].polynomialSplineDerivative(); // ddy
// Extracting the gradients as a function of position along the curve
double[] knots = splines[0].getKnots();
curvature = new TreeMap<>();
double w = (double) NNeighbours / 2d;
int startIdx = isLoop ? NNeighbours : 0;
int endIdx = isLoop ? knots.length - NNeighbours : knots.length;
for (int i = startIdx; i < endIdx; i++) {
double pos = knots[i];
double minPos = Math.max(knots[0], pos - w);
double maxPos = Math.min(knots[knots.length - 1], pos + w);
double width = maxPos - minPos;
double dx = (splines[0].value(maxPos) - splines[0].value(minPos)) / width;
double dy = (splines[1].value(maxPos) - splines[1].value(minPos)) / width;
double ddx = (splines[2].value(maxPos) - splines[2].value(minPos)) / (0.5 * width);
double ddy = (splines[3].value(maxPos) - splines[3].value(minPos)) / (0.5 * width);
double k = (dx * ddy - dy * ddx) / Math.pow((dx * dx + dy * dy), 3d / 2d);
curvature.put(knots[i], k);
}
}
public TreeMap getCurvature() {
if (curvature == null)
calculateCurvature();
return curvature;
}
public void showOverlay(ImagePlus ipl, int[] position, int lineWidth) {
if (curvature == null)
calculateCurvature();
// If fitting failed, curvature will still be null
if (curvature == null)
return;
// Calculating maximum curvature
double maxCurvature = Double.MIN_VALUE;
for (double currentCurvature : curvature.values()) {
maxCurvature = Math.max(Math.abs(currentCurvature), maxCurvature);
}
showOverlay(ipl, maxCurvature, position, lineWidth);
}
public void showOverlay(ImagePlus ipl, double maxCurvature, int[] position, double lineWidth) {
if (curvature == null)
calculateCurvature();
Overlay ovl = ipl.getOverlay();
if (ovl == null)
ovl = new Overlay();
if (curvature.size() < 2) {
ipl.setOverlay(ovl);
return;
}
Iterator iterator = curvature.keySet().iterator();
double p2 = iterator.next();
while (iterator.hasNext()) {
double p1 = p2;
p2 = iterator.next();
double x1 = splines[0].value(p1);
double y1 = splines[1].value(p1);
double x2 = splines[0].value(p2);
double y2 = splines[1].value(p2);
double k1 = Math.abs(curvature.get(p1));
double k2 = Math.abs(curvature.get(p2));
double k = (k1 + k2) / 2;
double b = k / (maxCurvature * 1.5);
Color color = Color.getHSBColor((float) b, 1f, 1f);
Line line = new Line(x1, y1, x2, y2);
line.setStrokeWidth(lineWidth);
line.setStrokeColor(color);
line.setPosition(position[2]);
ovl.addElement(line);
}
ipl.setOverlay(ovl);
}
public int getNNeighbours() {
return NNeighbours;
}
public void setNNeighbours(int NNeighbours) {
this.NNeighbours = NNeighbours;
// Calculating the bandwidth
loessBandwidth = (double) NNeighbours / (double) path.size();
// Keeping the bandwidth within limits
if (loessBandwidth < 0)
loessBandwidth = 0;
if (loessBandwidth > 1)
loessBandwidth = 1;
}
public double getLoessBandwidth() {
return loessBandwidth;
}
public int getLoessIterations() {
return loessIterations;
}
public void setLoessIterations(int loessIterations) {
this.loessIterations = loessIterations;
}
public double getLoessAccuracy() {
return loessAccuracy;
}
public void setLoessAccuracy(double loessAccuracy) {
this.loessAccuracy = loessAccuracy;
}
public FittingMethod getFittingMethod() {
return fittingMethod;
}
public void setFittingMethod(FittingMethod fittingMethod) {
this.fittingMethod = fittingMethod;
}
public ArrayList> getSpline() {
if (splines == null)
calculateCurvature();
ArrayList> spline = new ArrayList<>();
double[] knots = splines[0].getKnots();
int startIdx = isLoop ? NNeighbours : 0;
int endIdx = isLoop ? knots.length - NNeighbours : knots.length;
for (int i = startIdx; i < endIdx; i++) {
double t = knots[i];
int x = (int) Math.round(splines[0].value(t));
int y = (int) Math.round(splines[1].value(t));
spline.add(new Point(x, y, 0));
}
return spline;
}
}