Search in sources :

Example 6 with Line

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");
    }
}
Also used : TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) ArrayList(java.util.ArrayList) IJTrackProgress(gdsc.core.ij.IJTrackProgress) Cluster(gdsc.core.clustering.Cluster) Plot2(ij.gui.Plot2) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(gdsc.core.clustering.ClusterPoint) ClusteringEngine(gdsc.core.clustering.ClusteringEngine) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Example 7 with Line

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;
}
Also used : Plot2(ij.gui.Plot2) BinomialDistribution(org.apache.commons.math3.distribution.BinomialDistribution) ClusterPoint(gdsc.core.clustering.ClusterPoint)

Example 8 with Line

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]);
        }
    }
}
Also used : Point2D(java.awt.geom.Point2D) Path2D(java.awt.geom.Path2D) PolynomialSplineFunction(org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction) Line2D(java.awt.geom.Line2D)

Example 9 with Line

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();
}
Also used : Line(org.twak.utils.Line) Gaussian(org.apache.commons.math3.analysis.function.Gaussian)

Example 10 with Line

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());
}
Also used : Arrays(java.util.Arrays) Location(org.apache.commons.math3.geometry.partitioning.Region.Location) ConsecutivePairs(org.twak.utils.collections.ConsecutivePairs) Segment(org.apache.commons.math3.geometry.euclidean.twod.Segment) ConvexHull2D(org.apache.commons.math3.geometry.euclidean.twod.hull.ConvexHull2D) MultiMapSet(org.twak.utils.collections.MultiMapSet) HashMap(java.util.HashMap) Graph2D(org.twak.utils.geom.Graph2D) Pair(org.twak.utils.Pair) ConsecutiveItPairs(org.twak.utils.collections.ConsecutiveItPairs) Vector2D(org.apache.commons.math3.geometry.euclidean.twod.Vector2D) ArrayList(java.util.ArrayList) HashSet(java.util.HashSet) Loop(org.twak.utils.collections.Loop) Map(java.util.Map) Mathz(org.twak.utils.Mathz) PaintThing(org.twak.utils.PaintThing) Euclidean2D(org.apache.commons.math3.geometry.euclidean.twod.Euclidean2D) MathIllegalArgumentException(org.apache.commons.math3.exception.MathIllegalArgumentException) LoopL(org.twak.utils.collections.LoopL) UnionWalker(org.twak.utils.geom.UnionWalker) Iterator(java.util.Iterator) Collection(java.util.Collection) Line(org.twak.utils.Line) Set(java.util.Set) MonotoneChain(org.apache.commons.math3.geometry.euclidean.twod.hull.MonotoneChain) Region(org.apache.commons.math3.geometry.partitioning.Region) Vector2d(javax.vecmath.Vector2d) LinearForm(org.twak.utils.geom.LinearForm) Collectors(java.util.stream.Collectors) Point2d(javax.vecmath.Point2d) List(java.util.List) Paint(java.awt.Paint) ItComb(org.twak.utils.collections.ItComb) Loopable(org.twak.utils.collections.Loopable) InsufficientDataException(org.apache.commons.math3.exception.InsufficientDataException) Anglez(org.twak.utils.geom.Anglez) Collections(java.util.Collections) ConsecutiveItPairs(org.twak.utils.collections.ConsecutiveItPairs) ArrayList(java.util.ArrayList) LinearForm(org.twak.utils.geom.LinearForm) Paint(java.awt.Paint) Line(org.twak.utils.Line) Point2d(javax.vecmath.Point2d) HashSet(java.util.HashSet) MultiMapSet(org.twak.utils.collections.MultiMapSet)

Aggregations

ArrayList (java.util.ArrayList)9 Line (org.twak.utils.Line)8 Plot2 (ij.gui.Plot2)7 List (java.util.List)7 Map (java.util.Map)7 Set (java.util.Set)7 Point2d (javax.vecmath.Point2d)7 HashSet (java.util.HashSet)6 Collectors (java.util.stream.Collectors)6 Pair (org.twak.utils.Pair)6 Iterator (java.util.Iterator)5 Point3d (javax.vecmath.Point3d)5 Vector2d (javax.vecmath.Vector2d)5 PointValuePair (org.apache.commons.math3.optim.PointValuePair)5 Tweed (org.twak.tweed.Tweed)5 TweedSettings (org.twak.tweed.TweedSettings)5 Material (com.jme3.material.Material)4 ColorRGBA (com.jme3.math.ColorRGBA)4 Geometry (com.jme3.scene.Geometry)4 Mesh (com.jme3.scene.Mesh)4