use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class PcPalmMolecules method calculateAveragePrecision.
/**
* Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision
* distribution.
*
* <p>A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does
* not fit within 3 SDs of the simple mean then the simple mean is returned.
*
* @param molecules the molecules
* @param title the plot title (null if no plot should be displayed)
* @param histogramBins the histogram bins
* @param logFitParameters Record the fit parameters to the ImageJ log
* @param removeOutliers The distribution is created using all values within 1.5x the
* inter-quartile range (IQR) of the data
* @return The average precision
*/
public double calculateAveragePrecision(List<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
// Plot histogram of the precision
final float[] data = new float[molecules.size()];
final DescriptiveStatistics stats = new DescriptiveStatistics();
double yMin = Double.NEGATIVE_INFINITY;
double yMax = 0;
for (int i = 0; i < data.length; i++) {
data[i] = (float) molecules.get(i).precision;
stats.addValue(data[i]);
}
// Set the min and max y-values using 1.5 x IQR
if (removeOutliers) {
final double lower = stats.getPercentile(25);
final double upper = stats.getPercentile(75);
if (Double.isNaN(lower) || Double.isNaN(upper)) {
if (logFitParameters) {
ImageJUtils.log("Error computing IQR: %f - %f", lower, upper);
}
} else {
final double iqr = upper - lower;
yMin = Math.max(lower - iqr, stats.getMin());
yMax = Math.min(upper + iqr, stats.getMax());
if (logFitParameters) {
ImageJUtils.log(" Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
}
}
}
if (yMin == Double.NEGATIVE_INFINITY) {
yMin = stats.getMin();
yMax = stats.getMax();
if (logFitParameters) {
ImageJUtils.log(" Data range: %f - %f", yMin, yMax);
}
}
int bins;
if (histogramBins <= 0) {
bins = (int) Math.ceil((stats.getMax() - stats.getMin()) / HistogramPlot.getBinWidthScottsRule(stats.getStandardDeviation(), (int) stats.getN()));
} else {
bins = histogramBins;
}
final float[][] hist = HistogramPlot.calcHistogram(data, yMin, yMax, bins);
Plot plot = null;
if (title != null) {
plot = new Plot(title, "Precision", "Frequency");
final float[] xValues = hist[0];
final float[] yValues = hist[1];
plot.addPoints(xValues, yValues, Plot.BAR);
ImageJUtils.display(title, plot);
}
// Extract non-zero data
float[] x = Arrays.copyOf(hist[0], hist[0].length);
float[] y = Arrays.copyOf(hist[1], hist[1].length);
int count = 0;
for (int i = 0; i < y.length; i++) {
if (y[i] > 0) {
x[count] = x[i];
y[count] = y[i];
count++;
}
}
x = Arrays.copyOf(x, count);
y = Arrays.copyOf(y, count);
// Sense check to fitted data. Get mean and SD of histogram
final double[] stats2 = HistogramPlot.getHistogramStatistics(x, y);
double mean = stats2[0];
if (logFitParameters) {
log(" Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
}
// Standard Gaussian fit
final double[] parameters = fitGaussian(x, y);
if (parameters == null) {
log(" Failed to fit initial Gaussian");
return mean;
}
double newMean = parameters[1];
double error = Math.abs(stats2[0] - newMean) / stats2[1];
if (error > 3) {
log(" Failed to fit Gaussian: %f standard deviations from histogram mean", error);
return mean;
}
if (newMean < yMin || newMean > yMax) {
log(" Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
return mean;
}
mean = newMean;
if (logFitParameters) {
log(" Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
}
final double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
// Fit to a skewed Gaussian (or appropriate function)
final double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
if (skewParameters == null) {
log(" Failed to fit Skewed Gaussian");
return mean;
}
final SkewNormalFunction sn = new SkewNormalFunction(skewParameters);
if (logFitParameters) {
log(" Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance()));
}
newMean = sn.getMean();
error = Math.abs(stats2[0] - newMean) / stats2[1];
if (error > 3) {
log(" Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error);
return mean;
}
if (newMean < yMin || newMean > yMax) {
log(" Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
return mean;
}
// Use original histogram x-axis to maintain all the bins
if (plot != null) {
plot.setColor(Color.red);
addToPlot(plot, hist[0], skewParameters, Plot.LINE);
plot.setColor(Color.black);
ImageJUtils.display(title, plot);
}
// Return the average precision from the fitted curve
return newMean;
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class DoubletAnalysis method plotJaccard.
private void plotJaccard(ResidualsScore residualsScore, ResidualsScore residualsScoreReference) {
final String title = TITLE + " Score";
final Plot plot = new Plot(title, "Residuals", "Score");
double max = getMax(0.01, residualsScore);
if (residualsScoreReference != null) {
max = getMax(max, residualsScore);
}
plot.setLimits(0, 1, 0, max * 1.05);
addLines(plot, residualsScore, 1);
if (residualsScoreReference != null) {
addLines(plot, residualsScoreReference, 0.5);
}
plot.setColor(Color.black);
plot.addLabel(0, 0, String.format("Residuals %s; Jaccard %s (Black); Precision %s (Blue); Recall %s (Red)", MathUtils.rounded(residualsScore.residuals[residualsScore.maxJaccardIndex]), MathUtils.rounded(residualsScore.jaccard[residualsScore.maxJaccardIndex]), MathUtils.rounded(residualsScore.precision[residualsScore.maxJaccardIndex]), MathUtils.rounded(residualsScore.recall[residualsScore.maxJaccardIndex])));
ImageJUtils.display(title, plot, windowOrganiser);
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class DoubletAnalysis method showCumulativeHistogram.
/**
* Show a cumulative histogram of the data.
*
* @param values the values
* @param xTitle The name of plotted statistic
*/
private void showCumulativeHistogram(double[] values, String xTitle) {
final double[][] h = MathUtils.cumulativeHistogram(values, false);
final String title = TITLE + " " + xTitle + " Cumulative";
final Plot plot = new Plot(title, xTitle, "Frequency");
plot.setColor(Color.blue);
plot.addPoints(h[0], h[1], Plot.BAR);
ImageJUtils.display(title, plot, windowOrganiser);
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class PsfDrift method displayPlot.
private double[][] displayPlot(String title, String yLabel, double[] x, double[] y, double[] se, LoessInterpolator loess, int start, int end) {
// Extract non NaN numbers
double[] newX = new double[x.length];
double[] newY = new double[x.length];
int count = 0;
for (int i = 0; i < x.length; i++) {
if (!Double.isNaN(y[i])) {
newX[count] = x[i];
newY[count] = y[i];
count++;
}
}
newX = Arrays.copyOf(newX, count);
newY = Arrays.copyOf(newY, count);
title = TITLE + " " + title;
final Plot plot = new Plot(title, "z (nm)", yLabel);
final double[] limitsx = MathUtils.limits(x);
double[] limitsy = new double[2];
if (se != null) {
if (count > 0) {
limitsy = new double[] { newY[0] - se[0], newY[0] + se[0] };
for (int i = 1; i < newY.length; i++) {
limitsy[0] = MathUtils.min(limitsy[0], newY[i] - se[i]);
limitsy[1] = MathUtils.max(limitsy[1], newY[i] + se[i]);
}
}
} else if (count > 0) {
limitsy = MathUtils.limits(newY);
}
final double rangex = Math.max(0.05 * (limitsx[1] - limitsx[0]), 0.1);
final double rangey = Math.max(0.05 * (limitsy[1] - limitsy[0]), 0.1);
plot.setLimits(limitsx[0] - rangex, limitsx[1] + rangex, limitsy[0] - rangey, limitsy[1] + rangey);
if (loess == null) {
addPoints(plot, Plot.LINE, newX, newY, x[start], x[end]);
} else {
addPoints(plot, Plot.DOT, newX, newY, x[start], x[end]);
newY = loess.smooth(newX, newY);
addPoints(plot, Plot.LINE, newX, newY, x[start], x[end]);
}
if (se != null) {
plot.setColor(Color.magenta);
for (int i = 0; i < x.length; i++) {
if (!Double.isNaN(y[i])) {
plot.drawLine(x[i], y[i] - se[i], x[i], y[i] + se[i]);
}
}
// Draw the start and end lines for the valid range
plot.setColor(Color.green);
plot.drawLine(x[start], limitsy[0], x[start], limitsy[1]);
plot.drawLine(x[end], limitsy[0], x[end], limitsy[1]);
} else {
// draw a line for the recall limit
plot.setColor(Color.magenta);
plot.drawLine(limitsx[0] - rangex, settings.recallLimit, limitsx[1] + rangex, settings.recallLimit);
}
ImageJUtils.display(title, plot, windowOrganiser);
return new double[][] { newX, newY };
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class PsfDrift method showHwhm.
private void showHwhm() {
// Build a list of suitable images
final List<String> titles = createImageList(false);
if (titles.isEmpty()) {
IJ.error(TITLE, "No suitable PSF images");
return;
}
final GenericDialog gd = new GenericDialog(TITLE);
gd.addMessage("Approximate the volume of the PSF as a Gaussian and\n" + "compute the equivalent Gaussian width.");
settings = Settings.load();
gd.addChoice("PSF", titles.toArray(new String[0]), settings.title);
gd.addCheckbox("Use_offset", settings.useOffset);
gd.addSlider("Smoothing", 0, 0.5, settings.smoothing);
gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
gd.showDialog();
if (gd.wasCanceled()) {
return;
}
settings.title = gd.getNextChoice();
settings.useOffset = gd.getNextBoolean();
settings.smoothing = gd.getNextNumber();
settings.save();
imp = WindowManager.getImage(settings.title);
if (imp == null) {
IJ.error(TITLE, "No PSF image for image: " + settings.title);
return;
}
psfSettings = getPsfSettings(imp);
if (psfSettings == null) {
IJ.error(TITLE, "No PSF settings for image: " + settings.title);
return;
}
final int size = imp.getStackSize();
final ImagePsfModel psf = createImagePsf(1, size, 1);
final double[] w0 = psf.getAllHwhm0();
final double[] w1 = psf.getAllHwhm1();
// Get current centre
final int centre = psfSettings.getCentreImage();
// Extract valid values (some can be NaN)
double[] sw0 = new double[w0.length];
double[] sw1 = new double[w1.length];
final TDoubleArrayList s0 = new TDoubleArrayList(w0.length);
final TDoubleArrayList s1 = new TDoubleArrayList(w0.length);
int c0 = 0;
int c1 = 0;
for (int i = 0; i < w0.length; i++) {
if (Double.isFinite(w0[i])) {
s0.add(i + 1);
sw0[c0++] = w0[i];
}
if (Double.isFinite(w1[i])) {
s1.add(i + 1);
sw1[c1++] = w1[i];
}
}
if (c0 == 0 && c1 == 0) {
IJ.error(TITLE, "No computed HWHM for image: " + settings.title);
return;
}
double[] slice0 = s0.toArray();
sw0 = Arrays.copyOf(sw0, c0);
double[] slice1 = s1.toArray();
sw1 = Arrays.copyOf(sw1, c1);
// Smooth
if (settings.smoothing > 0) {
final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
sw0 = loess.smooth(slice0, sw0);
sw1 = loess.smooth(slice1, sw1);
}
final TDoubleArrayList minWx = new TDoubleArrayList();
final TDoubleArrayList minWy = new TDoubleArrayList();
for (int i = 0; i < w0.length; i++) {
double weight = 0;
if (Double.isFinite(w0[i])) {
if (Double.isFinite(w1[i])) {
weight = w0[i] * w1[i];
} else {
weight = w0[i] * w0[i];
}
} else if (Double.isFinite(w1[i])) {
weight = w1[i] * w1[i];
}
if (weight != 0) {
minWx.add(i + 1);
minWy.add(Math.sqrt(weight));
}
}
// Smooth the combined line
final double[] cx = minWx.toArray();
double[] cy = minWy.toArray();
if (settings.smoothing > 0) {
final LoessInterpolator loess = new LoessInterpolator(settings.smoothing, 1);
cy = loess.smooth(cx, cy);
}
final int newCentre = SimpleArrayUtils.findMinIndex(cy);
// Convert to FWHM
final double fwhm = psfSettings.getFwhm();
// Widths are in pixels
final String title = TITLE + " HWHM";
final Plot plot = new Plot(title, "Slice", "HWHM (px)");
double[] limits = MathUtils.limits(sw0);
limits = MathUtils.limits(limits, sw1);
final double maxY = limits[1] * 1.05;
plot.setLimits(1, size, 0, maxY);
plot.setColor(Color.red);
plot.addPoints(slice0, sw0, Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(slice1, sw1, Plot.LINE);
plot.setColor(Color.magenta);
plot.addPoints(cx, cy, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, "X=red; Y=blue, Combined=Magenta");
final PlotWindow pw = ImageJUtils.display(title, plot);
// Show a non-blocking dialog to allow the centre to be updated ...
// Add a label and dynamically update when the centre is moved.
final NonBlockingExtendedGenericDialog gd2 = new NonBlockingExtendedGenericDialog(TITLE);
final double scale = psfSettings.getPixelSize();
// @formatter:off
ImageJUtils.addMessage(gd2, "Update the PSF information?\n \n" + "Current z-centre = %d, FHWM = %s px (%s nm)\n", centre, MathUtils.rounded(fwhm), MathUtils.rounded(fwhm * scale));
// @formatter:on
gd2.addSlider("z-centre", cx[0], cx[cx.length - 1], newCentre);
final TextField tf = gd2.getLastTextField();
gd2.addMessage("");
gd2.addAndGetButton("Reset", event -> tf.setText(Integer.toString(newCentre)));
final Label label = gd2.getLastLabel();
gd2.addCheckbox("Update_centre", settings.updateCentre);
gd2.addCheckbox("Update_HWHM", settings.updateHwhm);
gd2.enableYesNoCancel();
gd2.hideCancelButton();
final UpdateDialogListener dl = new UpdateDialogListener(cx, cy, maxY, newCentre, scale, pw, label);
gd2.addDialogListener(dl);
gd.addHelp(HelpUrls.getUrl("psf-hwhm"));
gd2.showDialog();
if (gd2.wasOKed() && (settings.updateCentre || settings.updateHwhm)) {
final ImagePSF.Builder b = psfSettings.toBuilder();
if (settings.updateCentre) {
b.setCentreImage(dl.centre);
}
if (settings.updateHwhm) {
b.setFwhm(dl.getFwhm());
}
imp.setProperty("Info", ImagePsfHelper.toString(b));
}
}
Aggregations