Search in sources :

Example 21 with Plot2

use of ij.gui.Plot2 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 22 with Plot2

use of ij.gui.Plot2 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 23 with Plot2

use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.

the class PCPALMClusters method fitBinomial.

/**
	 * Fit a zero-truncated Binomial to the cumulative histogram
	 * 
	 * @param histogramData
	 * @return
	 */
private double[] fitBinomial(HistogramData histogramData) {
    // Get the mean and sum of the input histogram
    double mean;
    double sum = 0;
    count = 0;
    for (int i = 0; i < histogramData.histogram[1].length; i++) {
        count += histogramData.histogram[1][i];
        sum += histogramData.histogram[1][i] * i;
    }
    mean = sum / count;
    String name = "Zero-truncated Binomial distribution";
    Utils.log("Mean cluster size = %s", Utils.rounded(mean));
    Utils.log("Fitting cumulative " + name);
    // Convert to a normalised double array for the binomial fitter
    double[] histogram = new double[histogramData.histogram[1].length];
    for (int i = 0; i < histogramData.histogram[1].length; i++) histogram[i] = histogramData.histogram[1][i] / count;
    // Plot the cumulative histogram
    String title = TITLE + " Cumulative Distribution";
    Plot2 plot = null;
    if (showCumulativeHistogram) {
        // Create a cumulative histogram for fitting
        double[] cumulativeHistogram = new double[histogram.length];
        sum = 0;
        for (int i = 0; i < histogram.length; i++) {
            sum += histogram[i];
            cumulativeHistogram[i] = sum;
        }
        double[] values = Utils.newArray(histogram.length, 0.0, 1.0);
        plot = new Plot2(title, "N", "Cumulative Probability", values, cumulativeHistogram);
        plot.setLimits(0, histogram.length - 1, 0, 1.05);
        plot.addPoints(values, cumulativeHistogram, Plot2.CIRCLE);
        Utils.display(title, plot);
    }
    // Do fitting for different N
    double bestSS = Double.POSITIVE_INFINITY;
    double[] parameters = null;
    int worse = 0;
    int N = histogram.length - 1;
    int min = minN;
    final boolean customRange = (minN > 1) || (maxN > 0);
    if (min > N)
        min = N;
    if (maxN > 0 && N > maxN)
        N = maxN;
    Utils.log("Fitting N from %d to %d%s", min, N, (customRange) ? " (custom-range)" : "");
    // Since varying the N should be done in integer steps do this
    // for n=1,2,3,... until the SS peaks then falls off (is worse then the best 
    // score several times in succession)
    BinomialFitter bf = new BinomialFitter(new IJLogger());
    bf.setMaximumLikelihood(maximumLikelihood);
    for (int n = min; n <= N; n++) {
        PointValuePair solution = bf.fitBinomial(histogram, mean, n, true);
        if (solution == null)
            continue;
        double p = solution.getPointRef()[0];
        Utils.log("Fitted %s : N=%d, p=%s. SS=%g", name, n, Utils.rounded(p), solution.getValue());
        if (bestSS > solution.getValue()) {
            bestSS = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSS < Double.POSITIVE_INFINITY) {
            if (++worse >= 3)
                break;
        }
        if (showCumulativeHistogram)
            addToPlot(n, p, title, plot, new Color((float) n / N, 0, 1f - (float) n / N));
    }
    // Add best it in magenta
    if (showCumulativeHistogram && parameters != null)
        addToPlot((int) parameters[0], parameters[1], title, plot, Color.magenta);
    return parameters;
}
Also used : Color(java.awt.Color) BinomialFitter(gdsc.smlm.fitting.BinomialFitter) Plot2(ij.gui.Plot2) ClusterPoint(gdsc.core.clustering.ClusterPoint) IJLogger(gdsc.core.ij.IJLogger) PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 24 with Plot2

use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.

the class BlinkEstimator method plot.

private void plot(String xAxisTitle, String yAxisTitle, double[] x, double[] y) {
    String title = TITLE + " " + yAxisTitle;
    Plot2 plot = new Plot2(title, xAxisTitle, yAxisTitle, x, y);
    Utils.display(title, plot);
}
Also used : Plot2(ij.gui.Plot2)

Example 25 with Plot2

use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.

the class BlinkEstimator method computeBlinkingRate.

double computeBlinkingRate(MemoryPeakResults results, boolean verbose) {
    parameters = null;
    increaseNFittedPoints = false;
    // Calculate the counts verses dark time curve
    double[] Ntd = calculateCounts(results, maxDarkTime, searchDistance, relativeDistance, verbose);
    double[] td = calculateTd(Ntd);
    if (verbose)
        Utils.log("  Estimate %.0f molecules at td = %.0f ms", Ntd[0], td[0]);
    Ntd = shift(Ntd);
    td = shift(td);
    // Fit curve
    parameters = fit(td, Ntd, nFittedPoints, verbose);
    if (parameters == null)
        return 0;
    // Display
    if (showPlots) {
        String title = TITLE + " Molecule Counts";
        Plot2 plot = new Plot2(title, "td (ms)", "Count", td, Ntd);
        Utils.display(title, plot);
        plot.setColor(Color.red);
        plot.addPoints(blinkingModel.getX(), blinkingModel.value(parameters), Plot2.CIRCLE);
        // Add the rest that is not fitted
        double[] xOther = new double[td.length - blinkingModel.size()];
        double[] yOther = new double[xOther.length];
        for (int i = 0, t = blinkingModel.size(); i < xOther.length; i++, t++) {
            xOther[i] = td[t];
            yOther[i] = blinkingModel.evaluate(td[t], parameters);
        }
        plot.setColor(Color.blue);
        plot.addPoints(xOther, yOther, Plot2.CROSS);
        Utils.display(title, plot);
    }
    // Check if the fitted curve asymptotes above the real curve
    if (blinkingModel.evaluate(td[Ntd.length - 1], parameters) < Ntd[Ntd.length - 1]) {
        if (verbose) {
            Utils.log("  *** Warning ***");
            Utils.log("  Fitted curve does not asymptote above real curve. Increase the number of fitted points to sample more of the overcounting regime");
            Utils.log("  ***************");
        }
        increaseNFittedPoints = true;
    }
    // Blinking rate is 1 + nBlinks
    double blinkingRate = 1 + parameters[1];
    if (verbose)
        Utils.log("  Blinking rate = %s", Utils.rounded(blinkingRate, 4));
    return blinkingRate;
}
Also used : Plot2(ij.gui.Plot2)

Aggregations

Plot2 (ij.gui.Plot2)42 PlotWindow (ij.gui.PlotWindow)17 Point (java.awt.Point)9 BasePoint (gdsc.core.match.BasePoint)8 Statistics (gdsc.core.utils.Statistics)6 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)6 MemoryPeakResults (gdsc.smlm.results.MemoryPeakResults)5 WeightedObservedPoint (org.apache.commons.math3.fitting.WeightedObservedPoint)5 ClusterPoint (gdsc.core.clustering.ClusterPoint)4 PeakResultPoint (gdsc.smlm.ij.plugins.ResultsMatchCalculator.PeakResultPoint)4 PeakResult (gdsc.smlm.results.PeakResult)4 StoredData (gdsc.core.utils.StoredData)3 WindowOrganiser (ij.plugin.WindowOrganiser)3 Rectangle (java.awt.Rectangle)3 ArrayList (java.util.ArrayList)3 LoessInterpolator (org.apache.commons.math3.analysis.interpolation.LoessInterpolator)3 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)3 DescriptiveStatistics (org.apache.commons.math3.stat.descriptive.DescriptiveStatistics)3 Cluster (gdsc.core.clustering.Cluster)2 ClusteringEngine (gdsc.core.clustering.ClusteringEngine)2