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);
}
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;
}
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");
}
}
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;
}
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);
}
Aggregations