use of org.apache.commons.math3.geometry.euclidean.twod.Line in project GDSC-SMLM by aherbert.
the class PCPALMMolecules method performDistanceAnalysis.
private void performDistanceAnalysis(double[][] intraHist, int p99) {
// We want to know the fraction of distances between molecules at the 99th percentile
// that are intra- rather than inter-molecule.
// Do single linkage clustering of closest pair at this distance and count the number of
// links that are inter and intra.
// Convert molecules for clustering
ArrayList<ClusterPoint> points = new ArrayList<ClusterPoint>(molecules.size());
for (Molecule m : molecules) // Precision was used to store the molecule ID
points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons));
ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, new IJTrackProgress());
IJ.showStatus("Clustering to check inter-molecule distances");
engine.setTrackJoins(true);
ArrayList<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]);
IJ.showStatus("");
if (clusters != null) {
double[] intraIdDistances = engine.getIntraIdDistances();
double[] interIdDistances = engine.getInterIdDistances();
int all = interIdDistances.length + intraIdDistances.length;
log(" * Fraction of inter-molecule particle linkage @ %s nm = %s %%", Utils.rounded(intraHist[0][p99], 4), (all > 0) ? Utils.rounded(100.0 * interIdDistances.length / all, 4) : "0");
// Show a double cumulative histogram plot
double[][] intraIdHist = Maths.cumulativeHistogram(intraIdDistances, false);
double[][] interIdHist = Maths.cumulativeHistogram(interIdDistances, false);
// Plot
String title = TITLE + " molecule linkage distance";
Plot2 plot = new Plot2(title, "Distance", "Frequency", intraIdHist[0], intraIdHist[1]);
double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0;
if (interIdHist[1].length > 0)
max = FastMath.max(max, interIdHist[1][interIdHist[1].length - 1]);
plot.setLimits(0, intraIdHist[0][intraIdHist[0].length - 1], 0, max);
plot.setColor(Color.blue);
plot.addPoints(interIdHist[0], interIdHist[1], Plot2.LINE);
plot.setColor(Color.black);
Utils.display(title, plot);
} else {
log("Aborted clustering to check inter-molecule distances");
}
}
use of org.apache.commons.math3.geometry.euclidean.twod.Line in project GDSC-SMLM by aherbert.
the class PCPALMClusters method run.
/*
* (non-Javadoc)
*
* @see ij.plugin.PlugIn#run(java.lang.String)
*/
public void run(String arg) {
SMLMUsageTracker.recordPlugin(this.getClass(), arg);
if (!showDialog())
return;
PCPALMMolecules.logSpacer();
Utils.log(TITLE);
PCPALMMolecules.logSpacer();
long start = System.currentTimeMillis();
HistogramData histogramData;
if (fileInput) {
histogramData = loadHistogram(histogramFile);
} else {
histogramData = doClustering();
}
if (histogramData == null)
return;
float[][] hist = histogramData.histogram;
// Create a histogram of the cluster sizes
String title = TITLE + " Molecules/cluster";
String xTitle = "Molecules/cluster";
String yTitle = "Frequency";
// Create the data required for fitting and plotting
float[] xValues = Utils.createHistogramAxis(hist[0]);
float[] yValues = Utils.createHistogramValues(hist[1]);
// Plot the histogram
float yMax = Maths.max(yValues);
Plot2 plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
if (xValues.length > 0) {
double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
}
Utils.display(title, plot);
HistogramData noiseData = loadNoiseHistogram(histogramData);
if (noiseData != null) {
if (subtractNoise(histogramData, noiseData)) {
// Update the histogram
title += " (noise subtracted)";
xValues = Utils.createHistogramAxis(hist[0]);
yValues = Utils.createHistogramValues(hist[1]);
yMax = Maths.max(yValues);
plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
if (xValues.length > 0) {
double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, yMax * 1.05);
}
Utils.display(title, plot);
// Automatically save
if (autoSave) {
String newFilename = Utils.replaceExtension(histogramData.filename, ".noise.tsv");
if (saveHistogram(histogramData, newFilename)) {
Utils.log("Saved noise-subtracted histogram to " + newFilename);
}
}
}
}
// Fit the histogram
double[] fitParameters = fitBinomial(histogramData);
if (fitParameters != null) {
// Add the binomial to the histogram
int n = (int) fitParameters[0];
double p = fitParameters[1];
Utils.log("Optimal fit : N=%d, p=%s", n, Utils.rounded(p));
BinomialDistribution dist = new BinomialDistribution(n, p);
// A zero-truncated binomial was fitted.
// pi is the adjustment factor for the probability density.
double pi = 1 / (1 - dist.probability(0));
if (!fileInput) {
// Calculate the estimated number of clusters from the observed molecules:
// Actual = (Observed / p-value) / N
final double actual = (nMolecules / p) / n;
Utils.log("Estimated number of clusters : (%d / %s) / %d = %s", nMolecules, Utils.rounded(p), n, Utils.rounded(actual));
}
double[] x = new double[n + 2];
double[] y = new double[n + 2];
// Scale the values to match those on the histogram
final double normalisingFactor = count * pi;
for (int i = 0; i <= n; i++) {
x[i] = i + 0.5;
y[i] = dist.probability(i) * normalisingFactor;
}
x[n + 1] = n + 1.5;
y[n + 1] = 0;
// Redraw the plot since the limits may have changed
plot = new Plot2(title, xTitle, yTitle, xValues, yValues);
double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.maxDefault(yMax, y) * 1.05);
plot.setColor(Color.magenta);
plot.addPoints(x, y, Plot2.LINE);
plot.addPoints(x, y, Plot2.CIRCLE);
plot.setColor(Color.black);
Utils.display(title, plot);
}
double seconds = (System.currentTimeMillis() - start) / 1000.0;
String msg = TITLE + " complete : " + seconds + "s";
IJ.showStatus(msg);
Utils.log(msg);
return;
}
use of org.apache.commons.math3.geometry.euclidean.twod.Line in project triplea by triplea-game.
the class MapRouteDrawer method drawCurvedPath.
/**
* Draws a smooth curve through the given array of points
*
* <p>
* This algorithm is called Spline-Interpolation
* because the Apache-commons-math library we are using here does not accept
* values but {@code f(x)=y} with x having to increase all the time
* the idea behind this is to use a parameter array - the so called index
* as x array and splitting the points into a x and y coordinates array.
* </p>
*
* <p>
* Finally those 2 interpolated arrays get unified into a single {@linkplain Point2D} array and drawn to the Map
* </p>
*
* @param graphics The {@linkplain Graphics2D} Object to be drawn on
* @param points The Knot Points for the Spline-Interpolator aka the joints
*/
private void drawCurvedPath(final Graphics2D graphics, final Point2D[] points) {
final double[] index = createParameterizedIndex(points);
final PolynomialSplineFunction xcurve = splineInterpolator.interpolate(index, getValues(points, Point2D::getX));
final double[] xcoords = getCoords(xcurve, index);
final PolynomialSplineFunction ycurve = splineInterpolator.interpolate(index, getValues(points, Point2D::getY));
final double[] ycoords = getCoords(ycurve, index);
final List<Path2D> paths = routeCalculator.getAllNormalizedLines(xcoords, ycoords);
for (final Path2D path : paths) {
drawTransformedShape(graphics, path);
}
// draws the Line to the Cursor on every possible screen, so that the line ends at the cursor no matter what...
final List<Point2D[]> finishingPoints = routeCalculator.getAllPoints(new Point2D.Double(xcoords[xcoords.length - 1], ycoords[ycoords.length - 1]), points[points.length - 1]);
final boolean hasArrowEnoughSpace = points[points.length - 2].distance(points[points.length - 1]) > ARROW_LENGTH;
for (final Point2D[] finishingPointArray : finishingPoints) {
drawTransformedShape(graphics, new Line2D.Double(finishingPointArray[0], finishingPointArray[1]));
if (hasArrowEnoughSpace) {
drawArrow(graphics, finishingPointArray[0], finishingPointArray[1]);
}
}
}
use of org.apache.commons.math3.geometry.euclidean.twod.Line in project chordatlas by twak.
the class Slice method findSupportPoints.
public void findSupportPoints() {
supportPoints.clear();
supportMax = 0;
// Map<Line, Double> totalLineSupport = new HashedMap();
double distSigma = 0.2;
Gaussian distanceG = new Gaussian(0, distSigma), angleG = new Gaussian(0, 0.3);
for (Line l : gis.allLines()) {
double length = l.length();
int noPoints = (int) Math.max(2, length / 0.1);
for (int n = 0; n <= noPoints; n++) {
SupportPoint sp = new SupportPoint(l, l.fromPPram(n / (double) noPoints));
// Double.MAX_VALUE;
sp.support = 0;
for (Line nearL : slice.getNear(sp.pt, distSigma * 3)) {
double angle = nearL.absAngle(sp.line);
// Math.min (sp.pt.distance(line.start), sp.pt.distance(line.end)) ;
double dist = nearL.distance(sp.pt, true);
sp.support += distanceG.value(dist) * angleG.value(angle) * nearL.length();
}
supportPoints.add(sp);
}
}
foundLines = null;
repaint();
}
use of org.apache.commons.math3.geometry.euclidean.twod.Line in project chordatlas by twak.
the class Concarnie method apply.
private List<Problem> apply(Problem problem, int count) {
if (depth > 10 || problem.soup.isEmpty()) {
problem.addPortal();
return Collections.emptyList();
}
Set<Line> portals = new HashSet<>();
Set<Line> in = new HashSet<>(problem.soup);
for (Line sl : problem.hull) {
RangeMerge<Line> rm = new RangeMerge<>(P.CON_TOL, P.CON_TOL);
LinearForm lf = new LinearForm(sl);
for (Line l : problem.soup) {
if (onHull(sl, l)) {
double pps = lf.findPParam(l.start), ppe = lf.findPParam(l.end);
rm.add(pps, ppe, l);
in.remove(l);
}
}
List<Double> rmGet = rm.get();
if (rmGet.isEmpty()) {
if (!sl.start.equals(sl.end)) {
if (Double.isNaN(sl.start.x))
System.out.println("help!");
// whole thing is portal
portals.add(sl);
}
} else {
List<Point2d> occupied = new ArrayList();
{
double lf1 = lf.findPParam(sl.start), lf2 = lf.findPParam(sl.end);
if (lf1 > lf2) {
double tmp = lf1;
lf1 = lf2;
lf2 = tmp;
}
for (double d : rmGet) {
d = Mathz.clamp(d, lf1, lf2);
occupied.add(lf.fromPParam(d));
}
}
boolean onHull = false;
{
Point2d snapE = occupied.get(0), snapS = occupied.get(occupied.size() - 1);
if (snapS.distance(sl.start) < P.CON_TOL)
snapS.set(sl.start);
else {
occupied.add(new Point2d(sl.start));
}
if (snapE.distance(sl.end) < P.CON_TOL)
snapE.set(sl.end);
else {
occupied.add(0, new Point2d(sl.end));
onHull = true;
}
}
for (Pair<Point2d, Point2d> pair : new ConsecutiveItPairs<Point2d>(occupied)) {
onHull = !onHull;
if (pair.first().equals(pair.second()))
continue;
Line line = new Line(pair.first(), pair.second());
if (onHull) {
if (depth % 2 == 0)
line = line.reverse();
graph.add(line);
} else {
portals.add(line.reverse());
}
}
}
}
if (in.size() == problem.soup.size()) {
// we didn't do anything! remove something, anything...
List<Line> d = new ArrayList(in);
Collections.sort(d, (a, b) -> Double.compare(a.length(), b.length()));
for (// remove the shortest 1/3 lines
int i = 0; // remove the shortest 1/3 lines
i < Math.max(1, in.size() * 0.33); // remove the shortest 1/3 lines
i++) in.remove(d.get(i));
}
List<Portal> mergedPortals = mergeConsecutive(portals);
// assign each member of in to a portal
MultiMapSet<Portal, Line> subproblems = new MultiMapSet();
if (!mergedPortals.isEmpty()) {
// O(n^3) closest-clique assignment
MultiMapSet<Portal, Line> sub2 = new MultiMapSet();
for (Portal p : mergedPortals) sub2.putAll(p, p.lines);
while (!in.isEmpty()) {
double bestDist = Double.MAX_VALUE;
Portal bestP = null;
Line bestL = null;
for (Line l : in) for (Portal p : sub2.keySet()) for (Line sl : sub2.get(p)) {
double dlsl = l.distance(sl);
if (// ignore lines a long way away
dlsl > Math.max(P.CON_TOL * 3, p.length * 0.5))
continue;
double dist = dlsl + 0.1 * l.distance(p.summary);
if (dist < bestDist) {
bestP = p;
bestDist = dist;
bestL = l;
}
}
if (bestL == null)
break;
in.remove(bestL);
double lenBestL = bestL.length();
if (lenBestL > P.CON_TOL && lenBestL > 0.5 * bestP.summary.length()) {
in.add(new Line(bestL.start, bestL.fromPPram(0.5)));
in.add(new Line(bestL.fromPPram(0.5), bestL.end));
} else {
subproblems.put(bestP, bestL);
sub2.put(bestP, bestL);
}
}
} else {
mergedPortals.add(null);
subproblems.putAll(null, in);
}
return mergedPortals.stream().map(x -> new Problem(x == null ? Collections.emptyList() : x.lines, subproblems.get(x))).collect(Collectors.toList());
}
Aggregations