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