use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class BackgroundEstimator method plot.
private void plot(WindowOrganiser wo, double[] xValues, double[] data1, double[] data2, double[] data3, String title, String title1, String title2, String title3) {
// Get limits
double[] a = Maths.limits(xValues);
double[] b = Maths.limits(data1);
b = Maths.limits(b, data2);
if (data3 != null)
b = Maths.limits(b, data3);
title = imp.getTitle() + " " + title;
Plot2 plot = new Plot2(title, "Slice", title);
double range = b[1] - b[0];
if (range == 0)
range = 1;
plot.setLimits(a[0], a[1], b[0] - 0.05 * range, b[1] + 0.05 * range);
plot.setColor(Color.blue);
plot.addPoints(xValues, data1, Plot2.LINE);
plot.draw();
Statistics stats = new Statistics(data1);
String label = String.format("%s (Blue) = %s +/- %s", title1, Utils.rounded(stats.getMean()), Utils.rounded(stats.getStandardDeviation()));
plot.setColor(Color.red);
plot.addPoints(xValues, data2, Plot2.LINE);
stats = new Statistics(data2);
label += String.format(", %s (Red) = %s +/- %s", title2, Utils.rounded(stats.getMean()), Utils.rounded(stats.getStandardDeviation()));
if (data3 != null) {
plot.setColor(Color.green);
plot.addPoints(xValues, data3, Plot2.LINE);
stats = new Statistics(data3);
label += String.format(", %s (Green) = %s +/- %s", title3, Utils.rounded(stats.getMean()), Utils.rounded(stats.getStandardDeviation()));
}
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
PlotWindow pw = Utils.display(title, plot);
if (Utils.isNewWindow())
wo.add(pw);
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class LSQLVMGradientProcedureTest method gradientProcedureComputesSameOutputWithBias.
@Test
public void gradientProcedureComputesSameOutputWithBias() {
ErfGaussian2DFunction func = new SingleFreeCircularErfGaussian2DFunction(blockWidth, blockWidth);
int nparams = func.getNumberOfGradients();
int iter = 100;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
ArrayList<double[]> alphaList = new ArrayList<double[]>(iter);
ArrayList<double[]> betaList = new ArrayList<double[]>(iter);
ArrayList<double[]> xList = new ArrayList<double[]>(iter);
// Manipulate the background
double defaultBackground = Background;
try {
Background = 1e-2;
createData(1, iter, paramsList, yList, true);
EJMLLinearSolver solver = new EJMLLinearSolver(1e-5, 1e-6);
for (int i = 0; i < paramsList.size(); i++) {
double[] y = yList.get(i);
double[] a = paramsList.get(i);
BaseLSQLVMGradientProcedure p = LSQLVMGradientProcedureFactory.create(y, func);
p.gradient(a);
double[] beta = p.beta;
alphaList.add(p.getAlphaLinear());
betaList.add(beta.clone());
for (int j = 0; j < nparams; j++) {
if (Math.abs(beta[j]) < 1e-6)
System.out.printf("[%d] Tiny beta %s %g\n", i, func.getName(j), beta[j]);
}
// Solve
if (!solver.solve(p.getAlphaMatrix(), beta))
throw new AssertionError();
xList.add(beta);
//System.out.println(Arrays.toString(beta));
}
//for (int b = 1; b < 1000; b *= 2)
for (double b : new double[] { -500, -100, -10, -1, -0.1, 0, 0.1, 1, 10, 100, 500 }) {
Statistics[] rel = new Statistics[nparams];
Statistics[] abs = new Statistics[nparams];
for (int i = 0; i < nparams; i++) {
rel[i] = new Statistics();
abs[i] = new Statistics();
}
for (int i = 0; i < paramsList.size(); i++) {
double[] y = add(yList.get(i), b);
double[] a = paramsList.get(i).clone();
a[0] += b;
BaseLSQLVMGradientProcedure p = LSQLVMGradientProcedureFactory.create(y, func);
p.gradient(a);
double[] beta = p.beta;
double[] alpha2 = alphaList.get(i);
double[] beta2 = betaList.get(i);
double[] x2 = xList.get(i);
Assert.assertArrayEquals("Beta", beta2, beta, 1e-10);
Assert.assertArrayEquals("Alpha", alpha2, p.getAlphaLinear(), 1e-10);
// Solve
solver.solve(p.getAlphaMatrix(), beta);
Assert.assertArrayEquals("X", x2, beta, 1e-10);
for (int j = 0; j < nparams; j++) {
rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
abs[j].add(Math.abs(x2[j] - beta[j]));
}
}
for (int i = 0; i < nparams; i++) System.out.printf("Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g\n", b, func.getName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation());
}
} finally {
Background = defaultBackground;
}
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class GradientCalculatorSpeedTest method gradientCalculatorComputesSameOutputWithBias.
@Test
public void gradientCalculatorComputesSameOutputWithBias() {
Gaussian2DFunction func = new SingleEllipticalGaussian2DFunction(blockWidth, blockWidth);
int nparams = func.getNumberOfGradients();
GradientCalculator calc = new GradientCalculator(nparams);
int n = func.size();
int iter = 100;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
ArrayList<double[]> yList = new ArrayList<double[]>(iter);
ArrayList<double[][]> alphaList = new ArrayList<double[][]>(iter);
ArrayList<double[]> betaList = new ArrayList<double[]>(iter);
ArrayList<double[]> xList = new ArrayList<double[]>(iter);
// Manipulate the background
double defaultBackground = Background;
try {
Background = 1e-2;
createData(1, iter, paramsList, yList, true);
EJMLLinearSolver solver = new EJMLLinearSolver(1e-5, 1e-6);
for (int i = 0; i < paramsList.size(); i++) {
double[] y = yList.get(i);
double[] a = paramsList.get(i);
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
calc.findLinearised(n, y, a, alpha, beta, func);
alphaList.add(alpha);
betaList.add(beta.clone());
for (int j = 0; j < nparams; j++) {
if (Math.abs(beta[j]) < 1e-6)
System.out.printf("[%d] Tiny beta %s %g\n", i, func.getName(j), beta[j]);
}
// Solve
if (!solver.solve(alpha, beta))
throw new AssertionError();
xList.add(beta);
//System.out.println(Arrays.toString(beta));
}
double[][] alpha = new double[nparams][nparams];
double[] beta = new double[nparams];
//for (int b = 1; b < 1000; b *= 2)
for (double b : new double[] { -500, -100, -10, -1, -0.1, 0, 0.1, 1, 10, 100, 500 }) {
Statistics[] rel = new Statistics[nparams];
Statistics[] abs = new Statistics[nparams];
for (int i = 0; i < nparams; i++) {
rel[i] = new Statistics();
abs[i] = new Statistics();
}
for (int i = 0; i < paramsList.size(); i++) {
double[] y = add(yList.get(i), b);
double[] a = paramsList.get(i).clone();
a[0] += b;
calc.findLinearised(n, y, a, alpha, beta, func);
double[][] alpha2 = alphaList.get(i);
double[] beta2 = betaList.get(i);
double[] x2 = xList.get(i);
Assert.assertArrayEquals("Beta", beta2, beta, 1e-10);
for (int j = 0; j < nparams; j++) {
Assert.assertArrayEquals("Alpha", alpha2[j], alpha[j], 1e-10);
}
// Solve
solver.solve(alpha, beta);
Assert.assertArrayEquals("X", x2, beta, 1e-10);
for (int j = 0; j < nparams; j++) {
rel[j].add(DoubleEquality.relativeError(x2[j], beta[j]));
abs[j].add(Math.abs(x2[j] - beta[j]));
}
}
for (int i = 0; i < nparams; i++) System.out.printf("Bias = %.2f : %s : Rel %g +/- %g: Abs %g +/- %g\n", b, func.getName(i), rel[i].getMean(), rel[i].getStandardDeviation(), abs[i].getMean(), abs[i].getStandardDeviation());
}
} finally {
Background = defaultBackground;
}
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class PSFCreator method run.
/*
* (non-Javadoc)
*
* @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
*/
public void run(ImageProcessor ip) {
loadConfiguration();
BasePoint[] spots = getSpots();
if (spots.length == 0) {
IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
return;
}
ImageStack stack = getImageStack();
final int width = imp.getWidth();
final int height = imp.getHeight();
final int currentSlice = imp.getSlice();
// Adjust settings for a single maxima
config.setIncludeNeighbours(false);
fitConfig.setDuplicateDistance(0);
ArrayList<double[]> centres = new ArrayList<double[]>(spots.length);
int iterations = 1;
LoessInterpolator loess = new LoessInterpolator(smoothing, iterations);
// TODO - The fitting routine may not produce many points. In this instance the LOESS interpolator
// fails to smooth the data very well. A higher bandwidth helps this but perhaps
// try a different smoothing method.
// For each spot
Utils.log(TITLE + ": " + imp.getTitle());
Utils.log("Finding spot locations...");
Utils.log(" %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
StoredDataStatistics averageSd = new StoredDataStatistics();
StoredDataStatistics averageA = new StoredDataStatistics();
Statistics averageRange = new Statistics();
MemoryPeakResults allResults = new MemoryPeakResults();
allResults.setName(TITLE);
allResults.setBounds(new Rectangle(0, 0, width, height));
MemoryPeakResults.addResults(allResults);
for (int n = 1; n <= spots.length; n++) {
BasePoint spot = spots[n - 1];
final int x = (int) spot.getX();
final int y = (int) spot.getY();
MemoryPeakResults results = fitSpot(stack, width, height, x, y);
allResults.addAllf(results.getResults());
if (results.size() < 5) {
Utils.log(" Spot %d: Not enough fit results %d", n, results.size());
continue;
}
// Get the results for the spot centre and width
double[] z = new double[results.size()];
double[] xCoord = new double[z.length];
double[] yCoord = new double[z.length];
double[] sd = new double[z.length];
double[] a = new double[z.length];
int i = 0;
for (PeakResult peak : results.getResults()) {
z[i] = peak.getFrame();
xCoord[i] = peak.getXPosition() - x;
yCoord[i] = peak.getYPosition() - y;
sd[i] = FastMath.max(peak.getXSD(), peak.getYSD());
a[i] = peak.getAmplitude();
i++;
}
// Smooth the amplitude plot
double[] smoothA = loess.smooth(z, a);
// Find the maximum amplitude
int maximumIndex = findMaximumIndex(smoothA);
// Find the range at a fraction of the max. This is smoothed to find the X/Y centre
int start = 0, stop = smoothA.length - 1;
double limit = smoothA[maximumIndex] * amplitudeFraction;
for (int j = 0; j < smoothA.length; j++) {
if (smoothA[j] > limit) {
start = j;
break;
}
}
for (int j = smoothA.length; j-- > 0; ) {
if (smoothA[j] > limit) {
stop = j;
break;
}
}
averageRange.add(stop - start + 1);
// Extract xy centre coords and smooth
double[] smoothX = new double[stop - start + 1];
double[] smoothY = new double[smoothX.length];
double[] smoothSd = new double[smoothX.length];
double[] newZ = new double[smoothX.length];
for (int j = start, k = 0; j <= stop; j++, k++) {
smoothX[k] = xCoord[j];
smoothY[k] = yCoord[j];
smoothSd[k] = sd[j];
newZ[k] = z[j];
}
smoothX = loess.smooth(newZ, smoothX);
smoothY = loess.smooth(newZ, smoothY);
smoothSd = loess.smooth(newZ, smoothSd);
// Since the amplitude is not very consistent move from this peak to the
// lowest width which is the in-focus spot.
maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
// Find the centre at the amplitude peak
double cx = smoothX[maximumIndex] + x;
double cy = smoothY[maximumIndex] + y;
int cz = (int) newZ[maximumIndex];
double csd = smoothSd[maximumIndex];
double ca = smoothA[maximumIndex + start];
// The average should weight the SD using the signal for each spot
averageSd.add(smoothSd[maximumIndex]);
averageA.add(ca);
if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
Utils.log(" Spot %d was ignored", n);
continue;
}
// Store result - it may have been moved interactively
maximumIndex += this.slice - cz;
cz = (int) newZ[maximumIndex];
csd = smoothSd[maximumIndex];
ca = smoothA[maximumIndex + start];
Utils.log(" Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f\n", n, cx, cy, cz, csd, ca);
centres.add(new double[] { cx, cy, cz, csd, n });
}
if (interactiveMode) {
imp.setSlice(currentSlice);
imp.setOverlay(null);
// Hide the amplitude and spot plots
Utils.hide(TITLE_AMPLITUDE);
Utils.hide(TITLE_PSF_PARAMETERS);
}
if (centres.isEmpty()) {
String msg = "No suitable spots could be identified centres";
Utils.log(msg);
IJ.error(TITLE, msg);
return;
}
// Find the limits of the z-centre
int minz = (int) centres.get(0)[2];
int maxz = minz;
for (double[] centre : centres) {
if (minz > centre[2])
minz = (int) centre[2];
else if (maxz < centre[2])
maxz = (int) centre[2];
}
IJ.showStatus("Creating PSF image");
// Create a stack that can hold all the data.
ImageStack psf = createStack(stack, minz, maxz, magnification);
// For each spot
Statistics stats = new Statistics();
boolean ok = true;
for (int i = 0; ok && i < centres.size(); i++) {
double progress = (double) i / centres.size();
final double increment = 1.0 / (stack.getSize() * centres.size());
IJ.showProgress(progress);
double[] centre = centres.get(i);
// Extract the spot
float[][] spot = new float[stack.getSize()][];
Rectangle regionBounds = null;
for (int slice = 1; slice <= stack.getSize(); slice++) {
ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height);
if (regionBounds == null)
regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
spot[slice - 1] = ie.crop(regionBounds);
}
int n = (int) centre[4];
final float b = getBackground(n, spot);
if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
Utils.log(" Spot %d was ignored", n);
continue;
}
stats.add(b);
// Adjust the centre using the crop
centre[0] -= regionBounds.x;
centre[1] -= regionBounds.y;
// This takes a long time so this should track progress
ok = addToPSF(maxz, magnification, psf, centre, spot, regionBounds, progress, increment, centreEachSlice);
}
if (interactiveMode) {
Utils.hide(TITLE_INTENSITY);
}
IJ.showProgress(1);
if (threadPool != null) {
threadPool.shutdownNow();
threadPool = null;
}
if (!ok || stats.getN() == 0)
return;
final double avSd = getAverage(averageSd, averageA, 2);
Utils.log(" Average background = %.2f, Av. SD = %s px", stats.getMean(), Utils.rounded(avSd, 4));
normalise(psf, maxz, avSd * magnification, false);
IJ.showProgress(1);
psfImp = Utils.display("PSF", psf);
psfImp.setSlice(maxz);
psfImp.resetDisplayRange();
psfImp.updateAndDraw();
double[][] fitCom = new double[2][psf.getSize()];
Arrays.fill(fitCom[0], Double.NaN);
Arrays.fill(fitCom[1], Double.NaN);
double fittedSd = fitPSF(psf, loess, maxz, averageRange.getMean(), fitCom);
// Compute the drift in the PSF:
// - Use fitted centre if available; otherwise find CoM for each frame
// - express relative to the average centre
double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / magnification);
double[] slice = Utils.newArray(psf.getSize(), 1, 1.0);
String title = TITLE + " CoM Drift";
Plot2 plot = new Plot2(title, "Slice", "Drift (nm)");
plot.addLabel(0, 0, "Red = X; Blue = Y");
//double[] limitsX = Maths.limits(com[0]);
//double[] limitsY = Maths.limits(com[1]);
double[] limitsX = getLimits(com[0]);
double[] limitsY = getLimits(com[1]);
plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
plot.setColor(Color.red);
plot.addPoints(slice, com[0], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(slice, com[1], Plot.DOT);
plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
Utils.display(title, plot);
// TODO - Redraw the PSF with drift correction applied.
// This means that the final image should have no drift.
// This is relevant when combining PSF images. It doesn't matter too much for simulations
// unless the drift is large.
// Add Image properties containing the PSF details
final double fwhm = getFWHM(psf, maxz);
psfImp.setProperty("Info", XmlUtils.toXML(new PSFSettings(maxz, nmPerPixel / magnification, nmPerSlice, stats.getN(), fwhm, createNote())));
Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, Utils.rounded(nmPerPixel / magnification, 3), Utils.rounded(nmPerSlice, 3), stats.getN(), Utils.rounded(fittedSd * nmPerPixel, 4), Utils.rounded(fwhm));
createInteractivePlots(psf, maxz, nmPerPixel / magnification, fittedSd * nmPerPixel);
IJ.showStatus("");
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class Noise method drawPlot.
/**
* Build a plot of the noise estimate from the current frame.
* Limit the preview to 100 frames.
*/
private void drawPlot() {
NoiseEstimator.Method method1 = NoiseEstimator.Method.values()[algorithm];
NoiseEstimator.Method method2 = NoiseEstimator.Method.values()[algorithm2];
IJ.showStatus("Estimating noise ...");
boolean twoMethods = method1 != method2;
boolean preserveResiduals = method1.name().contains("Residuals") && method2.name().contains("Residuals") && twoMethods;
int start = imp.getCurrentSlice();
int end = FastMath.min(imp.getStackSize(), start + 100);
int size = end - start + 1;
double[] xValues = new double[size];
double[] yValues1 = new double[size];
double[] yValues2 = (twoMethods) ? new double[size] : null;
ImageStack stack = imp.getImageStack();
Rectangle bounds = imp.getProcessor().getRoi();
float[] buffer = null;
for (int slice = start, i = 0; slice <= end; slice++, i++) {
IJ.showProgress(i, size);
final ImageProcessor ip = stack.getProcessor(slice);
buffer = ImageConverter.getData(ip.getPixels(), ip.getWidth(), ip.getHeight(), bounds, buffer);
final NoiseEstimator ne = new NoiseEstimator(buffer, bounds.width, bounds.height);
ne.preserveResiduals = preserveResiduals;
ne.setRange(lowestPixelsRange);
xValues[i] = slice;
yValues1[i] = ne.getNoise(method1);
if (twoMethods)
yValues2[i] = ne.getNoise(method2);
}
IJ.showProgress(1);
IJ.showStatus("Plotting noise ...");
// Get limits
double[] a = Tools.getMinMax(xValues);
double[] b1 = Tools.getMinMax(yValues1);
if (twoMethods) {
double[] b2 = Tools.getMinMax(yValues2);
b1[0] = FastMath.min(b1[0], b2[0]);
b1[1] = FastMath.max(b1[1], b2[1]);
}
String title = imp.getTitle() + " Noise";
Plot2 plot = new Plot2(title, "Slice", "Noise", xValues, yValues1);
double range = b1[1] - b1[0];
if (range == 0)
range = 1;
plot.setLimits(a[0], a[1], b1[0] - 0.05 * range, b1[1] + 0.05 * range);
plot.setColor(Color.blue);
plot.draw();
String label = String.format("Blue = %s", Utils.rounded(new Statistics(yValues1).getMean()));
if (twoMethods) {
plot.setColor(Color.red);
plot.addPoints(xValues, yValues2, Plot2.LINE);
label += String.format(", Red = %s", Utils.rounded(new Statistics(yValues2).getMean()));
}
plot.addLabel(0, 0, label);
Utils.display(title, plot);
IJ.showStatus("");
}
Aggregations