use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.
the class EMGainAnalysis method plotPMF.
private void plotPMF() {
if (!showPMFDialog())
return;
final int gaussWidth = 5;
int dummyBias = (int) Math.max(500, gaussWidth * _noise + 1);
double[] pmf = pdf(0, _photons, _gain, _noise, dummyBias);
double[] x = Utils.newArray(pmf.length, 0, 1.0);
double yMax = Maths.max(pmf);
// Truncate x
int max = 0;
double sum = 0;
double p = 1 - tail;
while (sum < p && max < pmf.length) {
sum += pmf[max];
if (sum > 0.5 && pmf[max] == 0)
break;
max++;
}
int min = pmf.length;
sum = 0;
p = 1 - head;
while (sum < p && min > 0) {
min--;
sum += pmf[min];
if (sum > 0.5 && pmf[min] == 0)
break;
}
//int min = (int) (dummyBias - gaussWidth * _noise);
pmf = Arrays.copyOfRange(pmf, min, max);
x = Arrays.copyOfRange(x, min, max);
// Get the approximation
double[] f = new double[x.length];
LikelihoodFunction fun;
double myNoise = _noise;
switch(approximation) {
case 3:
fun = new PoissonFunction(1.0 / _gain, true);
break;
case 2:
// The mean does not matter so just use zero
fun = PoissonGaussianFunction.createWithStandardDeviation(1.0 / _gain, 0, _noise);
break;
case 1:
myNoise = 0;
case 0:
default:
PoissonGammaGaussianFunction myFun = new PoissonGammaGaussianFunction(1.0 / _gain, myNoise);
myFun.setMinimumProbability(0);
fun = myFun;
}
double expected = _photons;
if (offset != 0)
expected += offset * expected / 100.0;
expected *= _gain;
//double sum2 = 0;
for (int i = 0; i < f.length; i++) {
// Adjust the x-values to remove the dummy bias
x[i] -= dummyBias;
f[i] = fun.likelihood(x[i], expected);
//sum += pmf[i];
//sum2 += f[i];
}
//System.out.printf("Approximation sum = %f : %f\n", sum ,sum2);
if (showApproximation)
yMax = Maths.maxDefault(yMax, f);
String label = String.format("Gain=%s, noise=%s, photons=%s", Utils.rounded(_gain), Utils.rounded(_noise), Utils.rounded(_photons));
Plot2 plot = new Plot2("PMF", "ADUs", "p");
plot.setLimits(x[0], x[x.length - 1], 0, yMax);
plot.setColor(Color.red);
plot.addPoints(x, pmf, Plot2.LINE);
if (showApproximation) {
plot.setColor(Color.blue);
plot.addPoints(x, f, Plot2.LINE);
}
plot.setColor(Color.magenta);
plot.drawLine(_photons * _gain, 0, _photons * _gain, yMax);
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
PlotWindow win1 = Utils.display("PMF", plot);
// Plot the difference between the actual and approximation
double[] delta = new double[f.length];
for (int i = 0; i < f.length; i++) {
if (pmf[i] == 0 && f[i] == 0)
continue;
if (relativeDelta)
delta[i] = DoubleEquality.relativeError(f[i], pmf[i]) * Math.signum(f[i] - pmf[i]);
else
delta[i] = f[i] - pmf[i];
}
Plot2 plot2 = new Plot2("PMF delta", "ADUs", (relativeDelta) ? "Relative delta" : "delta");
double[] limits = Maths.limits(delta);
plot2.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
plot2.setColor(Color.red);
plot2.addPoints(x, delta, Plot2.LINE);
plot2.setColor(Color.magenta);
plot2.drawLine(_photons * _gain, limits[0], _photons * _gain, limits[1]);
plot2.setColor(Color.black);
plot2.addLabel(0, 0, label + ((offset == 0) ? "" : ", expected = " + Utils.rounded(expected / _gain)));
PlotWindow win2 = Utils.display("PMF delta", plot2);
if (Utils.isNewWindow()) {
Point p2 = win2.getLocation();
p2.y += win1.getHeight();
win2.setLocation(p2);
}
}
use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.
the class DoubletAnalysis method showHistogram.
/**
* Show a cumulative histogram of the data.
*
* @param values
* the values
* @param xTitle
* The name of plotted statistic
*/
public void showHistogram(double[] values, String xTitle) {
double[][] h = Maths.cumulativeHistogram(values, false);
String title = TITLE + " " + xTitle + " Cumulative";
Plot2 plot = new Plot2(title, xTitle, "Frequency");
double xMax = h[0][h[0].length - 1];
double yMax = h[1][h[1].length - 1];
plot.setLimits(0, xMax, 0, 1.05 * yMax);
plot.setColor(Color.blue);
plot.addPoints(h[0], h[1], Plot2.BAR);
display(title, plot);
}
use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.
the class DiffusionRateTest method showExample.
private void showExample(int totalSteps, double diffusionSigma, RandomGenerator[] random) {
MoleculeModel m = new MoleculeModel(0, new double[3]);
float[] xValues = new float[totalSteps];
float[] x = new float[totalSteps];
float[] y = new float[totalSteps];
final double[] axis = (settings.getDiffusionType() == DiffusionType.LINEAR_WALK) ? nextVector() : null;
for (int j = 0; j < totalSteps; j++) {
if (settings.getDiffusionType() == DiffusionType.GRID_WALK)
m.walk(diffusionSigma, random);
else if (settings.getDiffusionType() == DiffusionType.LINEAR_WALK)
m.slide(diffusionSigma, axis, random[0]);
else
m.move(diffusionSigma, random);
x[j] = (float) (m.getX());
y[j] = (float) (m.getY());
xValues[j] = (float) ((j + 1) / settings.stepsPerSecond);
}
// Plot x and y coords on a timeline
String title = TITLE + " example coordinates";
Plot2 plot = new Plot2(title, "Time (seconds)", "Distance (um)");
float[] xUm = convertToUm(x);
float[] yUm = convertToUm(y);
float[] limits = Maths.limits(xUm);
limits = Maths.limits(limits, yUm);
plot.setLimits(0, totalSteps / settings.stepsPerSecond, limits[0], limits[1]);
plot.setColor(Color.red);
plot.addPoints(xValues, xUm, Plot2.LINE);
plot.setColor(Color.blue);
plot.addPoints(xValues, yUm, Plot2.LINE);
Utils.display(title, plot);
// Scale up and draw 2D position
for (int j = 0; j < totalSteps; j++) {
x[j] *= magnification;
y[j] *= magnification;
}
float[] limitsx = getLimits(x);
float[] limitsy = getLimits(y);
int width = (int) (limitsx[1] - limitsx[0]);
int height = (int) (limitsy[1] - limitsy[0]);
// Ensure we draw something, even it is a simple dot at the centre for no diffusion
if (width == 0) {
width = (int) (32 * magnification);
limitsx[0] = -width / 2;
}
if (height == 0) {
height = (int) (32 * magnification);
limitsy[0] = -height / 2;
}
ImageProcessor ip = new ByteProcessor(width, height);
// Adjust x and y using the minimum to centre
x[0] -= limitsx[0];
y[0] -= limitsy[0];
for (int j = 1; j < totalSteps; j++) {
// Adjust x and y using the minimum to centre
x[j] -= limitsx[0];
y[j] -= limitsy[0];
// Draw a line
ip.setColor(32 + (223 * j) / (totalSteps - 1));
ip.drawLine(round(x[j - 1]), round(y[j - 1]), round(x[j]), round(y[j]));
}
// Draw the final position
ip.putPixel((int) round(x[totalSteps - 1]), (int) round(y[totalSteps - 1]), 255);
ImagePlus imp = Utils.display(TITLE + " example", ip);
// Apply the fire lookup table
WindowManager.setTempCurrentImage(imp);
LutLoader lut = new LutLoader();
lut.run("fire");
WindowManager.setTempCurrentImage(null);
}
use of ij.gui.Plot2 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
* @param steps
* the steps
*/
private void plotJumpDistances(String title, DoubleData jumpDistances, int dimensions) {
// Cumulative histogram
// --------------------
double[] values = jumpDistances.values();
String title2 = title + " Cumulative Jump Distance " + dimensions + "D";
double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(values);
Plot2 jdPlot = new Plot2(title2, "Distance (um^2)", "Cumulative Probability", jdHistogram[0], jdHistogram[1]);
PlotWindow pw2 = Utils.display(title2, jdPlot);
if (Utils.isNewWindow())
idList[idCount++] = pw2.getImagePlus().getID();
// 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)
double estimatedD = simpleD * simpleSteps;
double max = Maths.max(values);
double[] x = Utils.newArray(1000, 0, max / 1000);
double k = dimensions / 2.0;
double mean = 4 * estimatedD;
GammaDistribution dist = new GammaDistribution(k, mean);
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);
Utils.display(title2, jdPlot);
// Histogram
// ---------
title2 = title + " Jump " + dimensions + "D";
int plotId = Utils.showHistogram(title2, jumpDistances, "Distance (um^2)", 0, 0, Math.max(20, values.length / 1000));
if (Utils.isNewWindow())
idList[idCount++] = plotId;
// Recompute the expected function
for (int i = 0; i < x.length; i++) y[i] = dist.density(x[i]);
// Scale to have the same area
if (Utils.xValues.length > 1) {
final double area1 = jumpDistances.size() * (Utils.xValues[1] - Utils.xValues[0]);
final double area2 = dist.cumulativeProbability(x[x.length - 1]);
final double scaleFactor = area1 / area2;
for (int i = 0; i < y.length; i++) y[i] *= scaleFactor;
}
jdPlot = Utils.plot;
jdPlot.setColor(Color.red);
jdPlot.addPoints(x, y, Plot.LINE);
Utils.display(WindowManager.getImage(plotId).getTitle(), jdPlot);
}
use of ij.gui.Plot2 in project GDSC-SMLM by aherbert.
the class BenchmarkSpotFit method showDoubleHistogram.
private double[] showDoubleHistogram(StoredDataStatistics[][] stats, final int i, WindowOrganiser wo, double[][] matchScores, double nPredicted) {
String xLabel = filterCriteria[i].name;
LowerLimit lower = filterCriteria[i].lower;
UpperLimit upper = filterCriteria[i].upper;
double[] j = null;
double[] metric = null;
double maxJ = 0;
if (i <= FILTER_PRECISION && (showFilterScoreHistograms || upper.requiresJaccard || lower.requiresJaccard)) {
// Jaccard score verses the range of the metric
Arrays.sort(matchScores, new Comparator<double[]>() {
public int compare(double[] o1, double[] o2) {
if (o1[i] < o2[i])
return -1;
if (o1[i] > o2[i])
return 1;
return 0;
}
});
final int scoreIndex = FILTER_PRECISION + 1;
int n = results.size();
double tp = 0;
double fp = 0;
j = new double[matchScores.length + 1];
metric = new double[j.length];
for (int k = 0; k < matchScores.length; k++) {
final double score = matchScores[k][scoreIndex];
tp += score;
fp += (1 - score);
j[k + 1] = tp / (fp + n);
metric[k + 1] = matchScores[k][i];
}
metric[0] = metric[1];
maxJ = Maths.max(j);
if (showFilterScoreHistograms) {
String title = TITLE + " Jaccard " + xLabel;
Plot plot = new Plot(title, xLabel, "Jaccard", metric, j);
// Remove outliers
double[] limitsx = Maths.limits(metric);
Percentile p = new Percentile();
double l = p.evaluate(metric, 25);
double u = p.evaluate(metric, 75);
double iqr = 1.5 * (u - l);
limitsx[1] = Math.min(limitsx[1], u + iqr);
plot.setLimits(limitsx[0], limitsx[1], 0, Maths.max(j));
PlotWindow pw = Utils.display(title, plot);
if (Utils.isNewWindow())
wo.add(pw);
}
}
// [0] is all
// [1] is matches
// [2] is no match
StoredDataStatistics s1 = stats[0][i];
StoredDataStatistics s2 = stats[1][i];
StoredDataStatistics s3 = stats[2][i];
if (s1.getN() == 0)
return new double[4];
DescriptiveStatistics d = s1.getStatistics();
double median = 0;
Plot2 plot = null;
String title = null;
if (showFilterScoreHistograms) {
median = d.getPercentile(50);
String label = String.format("n = %d. Median = %s nm", s1.getN(), Utils.rounded(median));
int id = Utils.showHistogram(TITLE, s1, xLabel, filterCriteria[i].minBinWidth, (filterCriteria[i].restrictRange) ? 1 : 0, 0, label);
if (id == 0) {
IJ.log("Failed to show the histogram: " + xLabel);
return new double[4];
}
if (Utils.isNewWindow())
wo.add(id);
title = WindowManager.getImage(id).getTitle();
// Reverse engineer the histogram settings
plot = Utils.plot;
double[] xValues = Utils.xValues;
int bins = xValues.length;
double yMin = xValues[0];
double binSize = xValues[1] - xValues[0];
double yMax = xValues[0] + (bins - 1) * binSize;
if (s2.getN() > 0) {
double[] values = s2.getValues();
double[][] hist = Utils.calcHistogram(values, yMin, yMax, bins);
if (hist[0].length > 0) {
plot.setColor(Color.red);
plot.addPoints(hist[0], hist[1], Plot2.BAR);
Utils.display(title, plot);
}
}
if (s3.getN() > 0) {
double[] values = s3.getValues();
double[][] hist = Utils.calcHistogram(values, yMin, yMax, bins);
if (hist[0].length > 0) {
plot.setColor(Color.blue);
plot.addPoints(hist[0], hist[1], Plot2.BAR);
Utils.display(title, plot);
}
}
}
// Do cumulative histogram
double[][] h1 = Maths.cumulativeHistogram(s1.getValues(), true);
double[][] h2 = Maths.cumulativeHistogram(s2.getValues(), true);
double[][] h3 = Maths.cumulativeHistogram(s3.getValues(), true);
if (showFilterScoreHistograms) {
title = TITLE + " Cumul " + xLabel;
plot = new Plot2(title, xLabel, "Frequency");
// Find limits
double[] xlimit = Maths.limits(h1[0]);
xlimit = Maths.limits(xlimit, h2[0]);
xlimit = Maths.limits(xlimit, h3[0]);
// Restrict using the inter-quartile range
if (filterCriteria[i].restrictRange) {
double q1 = d.getPercentile(25);
double q2 = d.getPercentile(75);
double iqr = (q2 - q1) * 2.5;
xlimit[0] = Maths.max(xlimit[0], median - iqr);
xlimit[1] = Maths.min(xlimit[1], median + iqr);
}
plot.setLimits(xlimit[0], xlimit[1], 0, 1.05);
plot.addPoints(h1[0], h1[1], Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(h2[0], h2[1], Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(h3[0], h3[1], Plot.LINE);
}
// Determine the maximum difference between the TP and FP
double maxx1 = 0;
double maxx2 = 0;
double max1 = 0;
double max2 = 0;
// We cannot compute the delta histogram, or use percentiles
if (s2.getN() == 0) {
upper = UpperLimit.ZERO;
lower = LowerLimit.ZERO;
}
final boolean requireLabel = (showFilterScoreHistograms && filterCriteria[i].requireLabel);
if (requireLabel || upper.requiresDeltaHistogram() || lower.requiresDeltaHistogram()) {
if (s2.getN() != 0 && s3.getN() != 0) {
LinearInterpolator li = new LinearInterpolator();
PolynomialSplineFunction f1 = li.interpolate(h2[0], h2[1]);
PolynomialSplineFunction f2 = li.interpolate(h3[0], h3[1]);
for (double x : h1[0]) {
if (x < h2[0][0] || x < h3[0][0])
continue;
try {
double v1 = f1.value(x);
double v2 = f2.value(x);
double diff = v2 - v1;
if (diff > 0) {
if (max1 < diff) {
max1 = diff;
maxx1 = x;
}
} else {
if (max2 > diff) {
max2 = diff;
maxx2 = x;
}
}
} catch (OutOfRangeException e) {
// Because we reached the end
break;
}
}
} else {
// Switch to percentiles if we have no delta histogram
if (upper.requiresDeltaHistogram())
upper = UpperLimit.NINETY_NINE_PERCENT;
if (lower.requiresDeltaHistogram())
lower = LowerLimit.ONE_PERCENT;
}
// System.out.printf("Bounds %s : %s, pos %s, neg %s, %s\n", xLabel, Utils.rounded(getPercentile(h2, 0.01)),
// Utils.rounded(maxx1), Utils.rounded(maxx2), Utils.rounded(getPercentile(h1, 0.99)));
}
if (showFilterScoreHistograms) {
// We use bins=1 on charts where we do not need a label
if (requireLabel) {
String label = String.format("Max+ %s @ %s, Max- %s @ %s", Utils.rounded(max1), Utils.rounded(maxx1), Utils.rounded(max2), Utils.rounded(maxx2));
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
}
PlotWindow pw = Utils.display(title, plot);
if (Utils.isNewWindow())
wo.add(pw.getImagePlus().getID());
}
// Now compute the bounds using the desired limit
double l, u;
switch(lower) {
case ONE_PERCENT:
l = getPercentile(h2, 0.01);
break;
case MAX_NEGATIVE_CUMUL_DELTA:
l = maxx2;
break;
case ZERO:
l = 0;
break;
case HALF_MAX_JACCARD_VALUE:
l = getValue(metric, j, maxJ * 0.5);
break;
default:
throw new RuntimeException("Missing lower limit method");
}
switch(upper) {
case MAX_POSITIVE_CUMUL_DELTA:
u = maxx1;
break;
case NINETY_NINE_PERCENT:
u = getPercentile(h2, 0.99);
break;
case NINETY_NINE_NINE_PERCENT:
u = getPercentile(h2, 0.999);
break;
case ZERO:
u = 0;
break;
case MAX_JACCARD2:
u = getValue(metric, j, maxJ) * 2;
//System.out.printf("MaxJ = %.4f @ %.3f\n", maxJ, u / 2);
break;
default:
throw new RuntimeException("Missing upper limit method");
}
double min = getPercentile(h1, 0);
double max = getPercentile(h1, 1);
return new double[] { l, u, min, max };
}
Aggregations