Search in sources :

Example 76 with Plot

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

the class DiffusionRateTest method plotJumpDistances.

/**
 * Plot a cumulative histogram and standard histogram of the jump distances.
 *
 * @param title the title
 * @param jumpDistances the jump distances
 * @param dimensions the number of dimensions for the jumps
 */
private void plotJumpDistances(String title, DoubleData jumpDistances, int dimensions) {
    // Cumulative histogram
    // --------------------
    final double[] values = jumpDistances.values();
    String title2 = title + " Cumulative Jump Distance " + dimensions + "D";
    final double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(values);
    Plot jdPlot = new Plot(title2, "Distance (um^2)", "Cumulative Probability");
    jdPlot.addPoints(jdHistogram[0], jdHistogram[1], Plot.LINE);
    ImageJUtils.display(title2, jdPlot, windowOrganiser);
    // Plot the expected function
    // This is the Chi-squared distribution: The sum of the squares of k independent
    // standard normal random variables with k = dimensions. It is a special case of
    // the gamma distribution. If the normals have non-unit variance the distribution
    // is scaled.
    // Chi ~ Gamma(k/2, 2) // using the scale parameterisation of the gamma
    // s^2 * Chi ~ Gamma(k/2, 2*s^2)
    // So if s^2 = 2D:
    // 2D * Chi ~ Gamma(k/2, 4D)
    final double estimatedD = pluginSettings.simpleD * pluginSettings.simpleSteps;
    final double max = MathUtils.max(values);
    final double[] x = SimpleArrayUtils.newArray(1000, 0, max / 1000);
    final double k = dimensions / 2.0;
    final double mean = 4 * estimatedD;
    final GammaDistribution dist = new GammaDistribution(null, k, mean);
    final double[] y = new double[x.length];
    for (int i = 0; i < x.length; i++) {
        y[i] = dist.cumulativeProbability(x[i]);
    }
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(title2, jdPlot);
    // Histogram
    // ---------
    title2 = title + " Jump " + dimensions + "D";
    final HistogramPlot histogramPlot = new HistogramPlotBuilder(title2, jumpDistances, "Distance (um^2)").build();
    // Assume the plot works
    histogramPlot.show(windowOrganiser);
    // Recompute the expected function
    for (int i = 0; i < x.length; i++) {
        y[i] = dist.density(x[i]);
    }
    // Scale to have the same area
    final double[] xvalues = histogramPlot.getPlotXValues();
    if (xvalues.length > 1) {
        final double area1 = jumpDistances.size() * (xvalues[1] - xvalues[0]);
        final double area2 = dist.cumulativeProbability(x[x.length - 1]);
        final double scale = area1 / area2;
        for (int i = 0; i < y.length; i++) {
            y[i] *= scale;
        }
    }
    jdPlot = histogramPlot.getPlot();
    jdPlot.setColor(Color.red);
    jdPlot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(histogramPlot.getPlotTitle(), jdPlot);
}
Also used : HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) GammaDistribution(org.apache.commons.math3.distribution.GammaDistribution)

Example 77 with Plot

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

the class EmGainAnalysis method simulateFromPdf.

/**
 * Random sample from the cumulative probability distribution function that is used during
 * fitting.
 *
 * @return The histogram
 */
private int[] simulateFromPdf() {
    final double step = getStepSize(settings.settingPhotons, settings.settingGain, settings.settingNoise);
    final Pdf pdf = pdf(0, step, settings.settingPhotons, settings.settingGain, settings.settingNoise);
    // Debug this
    final double[] g = pdf.probability;
    final double[] x = pdf.x;
    final Plot plot = new Plot(TITLE + " PDF", "ADU", "P");
    plot.addPoints(x, Arrays.copyOf(g, g.length), Plot.LINE);
    ImageJUtils.display(TITLE + " PDF", plot);
    // Get cumulative probability
    double sum = 0;
    for (int i = 0; i < g.length; i++) {
        final double p = g[i];
        g[i] += sum;
        sum += p;
    }
    for (int i = 0; i < g.length; i++) {
        g[i] /= sum;
    }
    // Ensure value of 1 at the end
    g[g.length - 1] = 1;
    // Randomly sample
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final int bias = (int) settings.settingBias;
    final int[] bins = new int[x.length];
    for (int i = 0; i < x.length; i++) {
        bins[i] = bias + (int) x[i];
    }
    final int[] h = new int[bins[bins.length - 1] + 1];
    final int steps = settings.simulationSize;
    final Ticker ticker = ImageJUtils.createTicker(steps, 1);
    for (int n = 0; n < steps; n++) {
        int index = binarySearch(g, rng.nextDouble());
        if (index < 0) {
            index = -(index + 1);
        }
        h[bins[index]]++;
        ticker.tick();
    }
    return h;
}
Also used : Plot(ij.gui.Plot) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) Point(java.awt.Point)

Example 78 with Plot

use of ij.gui.Plot 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
    final ArrayList<ClusterPoint> points = new ArrayList<>(settings.molecules.size());
    for (final Molecule m : settings.molecules) {
        // Precision was used to store the molecule ID
        points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons));
    }
    final ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, SimpleImageJTrackProgress.getInstance());
    IJ.showStatus("Clustering to check inter-molecule distances");
    engine.setTrackJoins(true);
    final List<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]);
    IJ.showStatus("");
    if (clusters != null) {
        final double[] intraIdDistances = engine.getIntraIdDistances();
        final double[] interIdDistances = engine.getInterIdDistances();
        final int all = interIdDistances.length + intraIdDistances.length;
        log("  * Fraction of inter-molecule particle linkage @ %s nm = %s %%", MathUtils.rounded(intraHist[0][p99], 4), (all > 0) ? MathUtils.rounded(100.0 * interIdDistances.length / all, 4) : "0");
        // Show a double cumulative histogram plot
        final double[][] intraIdHist = MathUtils.cumulativeHistogram(intraIdDistances, false);
        final double[][] interIdHist = MathUtils.cumulativeHistogram(interIdDistances, false);
        // Plot
        final String title = TITLE + " molecule linkage distance";
        final Plot plot = new Plot(title, "Distance", "Frequency");
        plot.addPoints(intraIdHist[0], intraIdHist[1], Plot.LINE);
        double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0;
        if (interIdHist[1].length > 0) {
            max = Math.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], Plot.LINE);
        plot.setColor(Color.black);
        plot.addLegend("Intra-molecule\nInter-molecule");
        ImageJUtils.display(title, plot);
    } else {
        log("Aborted clustering to check inter-molecule distances");
    }
}
Also used : Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) ArrayList(java.util.ArrayList) Cluster(uk.ac.sussex.gdsc.core.clustering.Cluster) WeightedObservedPoint(org.apache.commons.math3.fitting.WeightedObservedPoint) ClusterPoint(uk.ac.sussex.gdsc.core.clustering.ClusterPoint) ClusteringEngine(uk.ac.sussex.gdsc.core.clustering.ClusteringEngine) ClusterPoint(uk.ac.sussex.gdsc.core.clustering.ClusterPoint)

Example 79 with Plot

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

the class PcPalmClusters method fitBinomial.

/**
 * Fit a zero-truncated Binomial to the cumulative histogram.
 *
 * @param histogramData the histogram data
 * @return the double[]
 */
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;
    final String name = "Zero-truncated Binomial distribution";
    ImageJUtils.log("Mean cluster size = %s", MathUtils.rounded(mean));
    ImageJUtils.log("Fitting cumulative " + name);
    // Convert to a normalised double array for the binomial fitter
    final 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
    final String title = TITLE + " Cumulative Distribution";
    Plot plot = null;
    if (settings.showCumulativeHistogram) {
        // Create a cumulative histogram for fitting
        final double[] cumulativeHistogram = new double[histogram.length];
        sum = 0;
        for (int i = 0; i < histogram.length; i++) {
            sum += histogram[i];
            cumulativeHistogram[i] = sum;
        }
        final double[] values = SimpleArrayUtils.newArray(histogram.length, 0.0, 1.0);
        plot = new Plot(title, "N", "Cumulative Probability");
        plot.addPoints(values, cumulativeHistogram, Plot.LINE);
        plot.addPoints(values, cumulativeHistogram, Plot.CIRCLE);
        // plot.addPoints(values, cumulativeHistogram, Plot.BAR);
        plot.setLimits(0, histogram.length - 1, 0, 1.05);
        ImageJUtils.display(title, plot);
    }
    // Do fitting for different N
    double bestSs = Double.POSITIVE_INFINITY;
    double[] parameters = null;
    int worse = 0;
    int countN = histogram.length - 1;
    int min = settings.minN;
    final boolean customRange = (settings.minN > 1) || (settings.maxN > 0);
    if (min > countN) {
        min = countN;
    }
    if (settings.maxN > 0 && countN > settings.maxN) {
        countN = settings.maxN;
    }
    ImageJUtils.log("Fitting N from %d to %d%s", min, countN, (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)
    final BinomialFitter bf = new BinomialFitter(ImageJPluginLoggerHelper.getLogger(getClass()));
    bf.setMaximumLikelihood(settings.maximumLikelihood);
    for (int n = min; n <= countN; n++) {
        final PointValuePair solution = bf.fitBinomial(histogram, mean, n, true);
        if (solution == null) {
            continue;
        }
        final double p = solution.getPointRef()[0];
        ImageJUtils.log("Fitted %s : N=%d, p=%s. SS=%g", name, n, MathUtils.rounded(p), solution.getValue());
        if (bestSs > solution.getValue()) {
            bestSs = solution.getValue();
            parameters = new double[] { n, p };
            worse = 0;
        } else if (bestSs < Double.POSITIVE_INFINITY && ++worse >= 3) {
            break;
        }
        if (settings.showCumulativeHistogram) {
            addToPlot(n, p, title, plot, new Color((float) n / countN, 0, 1f - (float) n / countN));
        }
    }
    // Add best it in magenta
    if (settings.showCumulativeHistogram && parameters != null) {
        addToPlot((int) parameters[0], parameters[1], title, plot, Color.magenta);
    }
    return parameters;
}
Also used : Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) Color(java.awt.Color) BinomialFitter(uk.ac.sussex.gdsc.smlm.fitting.BinomialFitter) ClusterPoint(uk.ac.sussex.gdsc.core.clustering.ClusterPoint) PointValuePair(org.apache.commons.math3.optim.PointValuePair)

Example 80 with Plot

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

the class PcPalmClusters method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (!showDialog()) {
        return;
    }
    PcPalmMolecules.logSpacer();
    ImageJUtils.log(TITLE);
    PcPalmMolecules.logSpacer();
    final long start = System.currentTimeMillis();
    HistogramData histogramData;
    if (fileInput) {
        histogramData = loadHistogram(settings.histogramFile);
    } else {
        histogramData = doClustering();
    }
    if (histogramData == null) {
        return;
    }
    final float[][] hist = histogramData.histogram;
    // Create a histogram of the cluster sizes
    String title = TITLE + " Molecules/cluster";
    final String xTitle = "Molecules/cluster";
    final String yTitle = "Frequency";
    // Plot the histogram
    Plot plot = new Plot(title, xTitle, yTitle);
    plot.addPoints(hist[0], hist[1], Plot.BAR);
    ImageJUtils.display(title, plot);
    final HistogramData noiseData = loadNoiseHistogram(histogramData);
    if (noiseData != null && subtractNoise(histogramData, noiseData)) {
        // Update the histogram
        title += " (noise subtracted)";
        plot = new Plot(title, xTitle, yTitle);
        plot.addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(title, plot);
        // Automatically save
        if (autoSave) {
            final String newFilename = FileUtils.replaceExtension(histogramData.filename, ".noise.tsv");
            if (saveHistogram(histogramData, newFilename)) {
                ImageJUtils.log("Saved noise-subtracted histogram to " + newFilename);
            }
        }
    }
    // Fit the histogram
    final double[] fitParameters = fitBinomial(histogramData);
    if (fitParameters != null) {
        // Add the binomial to the histogram
        final int n = (int) fitParameters[0];
        final double p = fitParameters[1];
        ImageJUtils.log("Optimal fit : N=%d, p=%s", n, MathUtils.rounded(p));
        final BinomialDistribution dist = new BinomialDistribution(n, p);
        // A zero-truncated binomial was fitted.
        // pi is the adjustment factor for the probability density.
        final 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 = (numberOfMolecules / p) / n;
            ImageJUtils.log("Estimated number of clusters : (%d / %s) / %d = %s", numberOfMolecules, MathUtils.rounded(p), n, MathUtils.rounded(actual));
        }
        final double[] x = new double[n + 2];
        final 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;
            y[i] = dist.probability(i) * normalisingFactor;
        }
        x[n + 1] = n + 1;
        y[n + 1] = 0;
        // Redraw the plot since the limits may have changed
        plot.setColor(Color.magenta);
        plot.addPoints(x, y, Plot.LINE);
        plot.addPoints(x, y, Plot.CIRCLE);
        plot.setColor(Color.black);
        plot.updateImage();
        ImageJUtils.display(title, plot);
    }
    final double seconds = (System.currentTimeMillis() - start) / 1000.0;
    final String msg = TITLE + " complete : " + seconds + "s";
    IJ.showStatus(msg);
    ImageJUtils.log(msg);
}
Also used : Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) BinomialDistribution(org.apache.commons.math3.distribution.BinomialDistribution) ClusterPoint(uk.ac.sussex.gdsc.core.clustering.ClusterPoint)

Aggregations

Plot (ij.gui.Plot)89 HistogramPlot (uk.ac.sussex.gdsc.core.ij.HistogramPlot)20 Point (java.awt.Point)19 PlotWindow (ij.gui.PlotWindow)17 Color (java.awt.Color)13 WindowOrganiser (uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser)13 HistogramPlotBuilder (uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder)12 BasePoint (uk.ac.sussex.gdsc.core.match.BasePoint)12 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)11 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)11 Rectangle (java.awt.Rectangle)9 ArrayList (java.util.ArrayList)9 GenericDialog (ij.gui.GenericDialog)8 NonBlockingExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog)7 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)7 Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)7 StoredData (uk.ac.sussex.gdsc.core.utils.StoredData)7 StoredDataStatistics (uk.ac.sussex.gdsc.core.utils.StoredDataStatistics)7 ImagePlus (ij.ImagePlus)6 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)5