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 Gaussian2DFunctionTest method functionComputesTargetGradient.
private void functionComputesTargetGradient(int targetParameter) {
int gradientIndex = findGradientIndex(f1, targetParameter);
double[] dyda = new double[f1.gradientIndices().length];
double[] dyda2 = new double[dyda.length];
double[] a;
Gaussian2DFunction f1a = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
Gaussian2DFunction f1b = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
Statistics s = new Statistics();
for (double background : testbackground) // Peak 1
for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
f1.initialise(a);
// Numerically solve gradient.
// Calculate the step size h to be an exact numerical representation
final double xx = a[targetParameter];
// Get h to minimise roundoff error
double h = Precision.representableDelta(xx, h_);
// Evaluate at (x+h) and (x-h)
a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
a[targetParameter] = xx + h;
f1a.initialise(a);
a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
a[targetParameter] = xx - h;
f1b.initialise(a);
for (int x : testx) for (int y : testy) {
int i = y * maxx + x;
f1.eval(i, dyda);
double value2 = f1a.eval(i, dyda2);
double value3 = f1b.eval(i, dyda2);
double gradient = (value2 - value3) / (2 * h);
double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
s.add(error);
Assert.assertTrue(gradient + " sign != " + dyda2[gradientIndex], (gradient * dyda2[gradientIndex]) >= 0);
//System.out.printf("[%d,%d] %f == [%d] %f? (%g)\n", x, y, gradient,
// gradientIndex, dyda2[gradientIndex], error);
//System.out.printf("[%d,%d] %f == [%d] %f?\n", x, y, gradient, gradientIndex, dyda[gradientIndex]);
Assert.assertTrue(gradient + " != " + dyda[gradientIndex], eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]));
}
}
System.out.printf("functionComputesTargetGradient %s %s (error %s +/- %s)\n", f1.getClass().getSimpleName(), f1.getName(targetParameter), Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()));
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class BenchmarkFit method run.
private void run() {
// Initialise the answer. Convert to units of the image (ADUs and pixels)
answer[Gaussian2DFunction.BACKGROUND] = benchmarkParameters.getBackground() * benchmarkParameters.gain;
answer[Gaussian2DFunction.SIGNAL] = benchmarkParameters.getSignal() * benchmarkParameters.gain;
answer[Gaussian2DFunction.X_POSITION] = benchmarkParameters.x;
answer[Gaussian2DFunction.Y_POSITION] = benchmarkParameters.y;
answer[Gaussian2DFunction.X_SD] = benchmarkParameters.s / benchmarkParameters.a;
answer[Gaussian2DFunction.Y_SD] = benchmarkParameters.s / benchmarkParameters.a;
// Set up the fit region. Always round down since 0.5 is the centre of the pixel.
int x = (int) benchmarkParameters.x;
int y = (int) benchmarkParameters.y;
region = new Rectangle(x - regionSize, y - regionSize, 2 * regionSize + 1, 2 * regionSize + 1);
if (!new Rectangle(0, 0, imp.getWidth(), imp.getHeight()).contains(region)) {
// Check if it is incorrect by only 1 pixel
if (region.width <= imp.getWidth() + 1 && region.height <= imp.getHeight() + 1) {
Utils.log("Adjusting region %s to fit within image bounds (%dx%d)", region.toString(), imp.getWidth(), imp.getHeight());
region = new Rectangle(0, 0, imp.getWidth(), imp.getHeight());
} else {
IJ.error(TITLE, "Fit region does not fit within the image");
return;
}
}
// Adjust the centre & account for 0.5 pixel offset during fitting
x -= region.x;
y -= region.y;
answer[Gaussian2DFunction.X_POSITION] -= (region.x + 0.5);
answer[Gaussian2DFunction.Y_POSITION] -= (region.y + 0.5);
// Configure for fitting
fitConfig.setBackgroundFitting(backgroundFitting);
fitConfig.setNotSignalFitting(!signalFitting);
fitConfig.setComputeDeviations(false);
final ImageStack stack = imp.getImageStack();
// Create a pool of workers
int nThreads = Prefs.getThreads();
BlockingQueue<Integer> jobs = new ArrayBlockingQueue<Integer>(nThreads * 2);
List<Worker> workers = new LinkedList<Worker>();
List<Thread> threads = new LinkedList<Thread>();
for (int i = 0; i < nThreads; i++) {
Worker worker = new Worker(jobs, stack, region, fitConfig);
Thread t = new Thread(worker);
workers.add(worker);
threads.add(t);
t.start();
}
final int totalFrames = benchmarkParameters.frames;
// Store all the fitting results
results = new double[totalFrames * getNumberOfStartPoints()][];
resultsTime = new long[results.length];
// Fit the frames
totalProgress = totalFrames;
stepProgress = Utils.getProgressInterval(totalProgress);
progress = 0;
for (int i = 0; i < totalFrames; i++) {
// Only fit if there were simulated photons
if (benchmarkParameters.p[i] > 0) {
put(jobs, i);
}
}
// Finish all the worker threads by passing in a null job
for (int i = 0; i < threads.size(); i++) {
put(jobs, -1);
}
// Wait for all to finish
for (int i = 0; i < threads.size(); i++) {
try {
threads.get(i).join();
} catch (InterruptedException e) {
e.printStackTrace();
}
}
threads.clear();
if (comFitting)
Utils.log(TITLE + ": CoM within start offset = %d / %d (%s%%)", comValid.intValue(), totalFrames, Utils.rounded((100.0 * comValid.intValue()) / totalFrames));
IJ.showProgress(1);
IJ.showStatus("Collecting results ...");
// Collect the results
Statistics[] stats = new Statistics[NAMES.length];
for (int i = 0; i < workers.size(); i++) {
Statistics[] next = workers.get(i).stats;
for (int j = 0; j < next.length; j++) {
if (stats[j] == null)
stats[j] = next[j];
else
stats[j].add(next[j]);
}
}
workers.clear();
// Show a table of the results
summariseResults(stats);
// Optionally show histograms
if (showHistograms) {
IJ.showStatus("Calculating histograms ...");
int[] idList = new int[NAMES.length];
int count = 0;
double[] convert = getConversionFactors();
boolean requireRetile = false;
for (int i = 0; i < NAMES.length; i++) {
if (displayHistograms[i] && convert[i] != 0) {
// We will have to convert the values...
double[] tmp = ((StoredDataStatistics) stats[i]).getValues();
for (int j = 0; j < tmp.length; j++) tmp[j] *= convert[i];
StoredDataStatistics tmpStats = new StoredDataStatistics(tmp);
idList[count++] = Utils.showHistogram(TITLE, tmpStats, NAMES[i], 0, 0, histogramBins, String.format("%s +/- %s", Utils.rounded(tmpStats.getMean()), Utils.rounded(tmpStats.getStandardDeviation())));
requireRetile = requireRetile || Utils.isNewWindow();
}
}
if (count > 0 && requireRetile) {
idList = Arrays.copyOf(idList, count);
new WindowOrganiser().tileWindows(idList);
}
}
if (saveRawData) {
String dir = Utils.getDirectory("Data_directory", rawDataDirectory);
if (dir != null)
saveData(stats, dir);
}
IJ.showStatus("");
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class PCPALMMolecules method runSimulation.
private void runSimulation(boolean resultsAvailable) {
if (resultsAvailable && !showSimulationDialog())
return;
startLog();
log("Simulation parameters");
if (blinkingDistribution == 3) {
log(" - Clusters = %d", nMolecules);
log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4));
log(" - Molecules/cluster = %s", Utils.rounded(blinkingRate, 4));
log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]);
log(" - p-Value = %s", Utils.rounded(p, 4));
} else {
log(" - Molecules = %d", nMolecules);
log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4));
log(" - Blinking rate = %s", Utils.rounded(blinkingRate, 4));
log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]);
}
log(" - Average precision = %s nm", Utils.rounded(sigmaS, 4));
log(" - Clusters simulation = " + CLUSTER_SIMULATION[clusterSimulation]);
if (clusterSimulation > 0) {
log(" - Cluster number = %s +/- %s", Utils.rounded(clusterNumber, 4), Utils.rounded(clusterNumberSD, 4));
log(" - Cluster radius = %s nm", Utils.rounded(clusterRadius, 4));
}
final double nmPerPixel = 100;
double width = simulationSize * 1000.0;
// Allow a border of 3 x sigma for +/- precision
//if (blinkingRate > 1)
width -= 3 * sigmaS;
RandomGenerator randomGenerator = new Well19937c(System.currentTimeMillis() + System.identityHashCode(this));
RandomDataGenerator dataGenerator = new RandomDataGenerator(randomGenerator);
UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, randomGenerator.nextInt());
molecules = new ArrayList<Molecule>(nMolecules);
// Create some dummy results since the calibration is required for later analysis
results = new MemoryPeakResults();
results.setCalibration(new gdsc.smlm.results.Calibration(nmPerPixel, 1, 100));
results.setSource(new NullSource("Molecule Simulation"));
results.begin();
int count = 0;
// Generate a sequence of coordinates
ArrayList<double[]> xyz = new ArrayList<double[]>((int) (nMolecules * 1.1));
Statistics statsRadius = new Statistics();
Statistics statsSize = new Statistics();
String maskTitle = TITLE + " Cluster Mask";
ByteProcessor bp = null;
double maskScale = 0;
if (clusterSimulation > 0) {
// Simulate clusters.
// Note: In the Veatch et al. paper (Plos 1, e31457) correlation functions are built using circles
// with small radii of 4-8 Arbitrary Units (AU) or large radii of 10-30 AU. A fluctuations model is
// created at T = 1.075 Tc. It is not clear exactly how the particles are distributed.
// It may be that a mask is created first using the model. The particles are placed on the mask using
// a specified density. This simulation produces a figure to show either a damped cosine function
// (circles) or an exponential (fluctuations). The number of particles in each circle may be randomly
// determined just by density. The figure does not discuss the derivation of the cluster size
// statistic.
//
// If this plugin simulation is run with a uniform distribution and blinking rate of 1 then the damped
// cosine function is reproduced. The curve crosses g(r)=1 at a value equivalent to the average
// distance to the centre-of-mass of each drawn cluster, not the input cluster radius parameter (which
// is a hard upper limit on the distance to centre).
final int maskSize = lowResolutionImageSize;
int[] mask = null;
// scale is in nm/pixel
maskScale = width / maskSize;
ArrayList<double[]> clusterCentres = new ArrayList<double[]>();
int totalSteps = 1 + (int) Math.ceil(nMolecules / clusterNumber);
if (clusterSimulation == 2 || clusterSimulation == 3) {
// Clusters are non-overlapping circles
// Ensure the circles do not overlap by using an exclusion mask that accumulates
// out-of-bounds pixels by drawing the last cluster (plus some border) on an image. When no
// more pixels are available then stop generating molecules.
// This is done by cumulatively filling a mask and using the MaskDistribution to select
// a new point. This may be slow but it works.
// TODO - Allow clusters of different sizes...
mask = new int[maskSize * maskSize];
Arrays.fill(mask, 255);
MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator);
double[] centre;
IJ.showStatus("Computing clusters mask");
int roiRadius = (int) Math.round((clusterRadius * 2) / maskScale);
if (clusterSimulation == 3) {
// Generate a mask of circles then sample from that.
// If we want to fill the mask completely then adjust the total steps to be the number of
// circles that can fit inside the mask.
totalSteps = (int) (maskSize * maskSize / (Math.PI * Math.pow(clusterRadius / maskScale, 2)));
}
while ((centre = maskDistribution.next()) != null && clusterCentres.size() < totalSteps) {
IJ.showProgress(clusterCentres.size(), totalSteps);
// The mask returns the coordinates with the centre of the image at 0,0
centre[0] += width / 2;
centre[1] += width / 2;
clusterCentres.add(centre);
// Fill in the mask around the centre to exclude any more circles that could overlap
double cx = centre[0] / maskScale;
double cy = centre[1] / maskScale;
fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0);
//Utils.display("Mask", new ColorProcessor(maskSize, maskSize, mask));
try {
maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator);
} catch (IllegalArgumentException e) {
// This can happen when there are no more non-zero pixels
log("WARNING: No more room for clusters on the mask area (created %d of estimated %d)", clusterCentres.size(), totalSteps);
break;
}
}
IJ.showProgress(1);
IJ.showStatus("");
} else {
// Pick centres randomly from the distribution
while (clusterCentres.size() < totalSteps) clusterCentres.add(dist.next());
}
if (showClusterMask || clusterSimulation == 3) {
// Show the mask for the clusters
if (mask == null)
mask = new int[maskSize * maskSize];
else
Arrays.fill(mask, 0);
int roiRadius = (int) Math.round((clusterRadius) / maskScale);
for (double[] c : clusterCentres) {
double cx = c[0] / maskScale;
double cy = c[1] / maskScale;
fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1);
}
if (clusterSimulation == 3) {
// We have the mask. Now pick points at random from the mask.
MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator);
// Allocate each molecule position to a parent circle so defining clusters.
int[][] clusters = new int[clusterCentres.size()][];
int[] clusterSize = new int[clusters.length];
for (int i = 0; i < nMolecules; i++) {
double[] centre = maskDistribution.next();
// The mask returns the coordinates with the centre of the image at 0,0
centre[0] += width / 2;
centre[1] += width / 2;
xyz.add(centre);
// Output statistics on cluster size and number.
// TODO - Finding the closest cluster could be done better than an all-vs-all comparison
double max = distance2(centre, clusterCentres.get(0));
int cluster = 0;
for (int j = 1; j < clusterCentres.size(); j++) {
double d2 = distance2(centre, clusterCentres.get(j));
if (d2 < max) {
max = d2;
cluster = j;
}
}
// Assign point i to cluster
centre[2] = cluster;
if (clusterSize[cluster] == 0) {
clusters[cluster] = new int[10];
}
if (clusters[cluster].length <= clusterSize[cluster]) {
clusters[cluster] = Arrays.copyOf(clusters[cluster], (int) (clusters[cluster].length * 1.5));
}
clusters[cluster][clusterSize[cluster]++] = i;
}
// Generate real cluster size statistics
for (int j = 0; j < clusterSize.length; j++) {
final int size = clusterSize[j];
if (size == 0)
continue;
statsSize.add(size);
if (size == 1) {
statsRadius.add(0);
continue;
}
// Find centre of cluster and add the distance to each point
double[] com = new double[2];
for (int n = 0; n < size; n++) {
double[] xy = xyz.get(clusters[j][n]);
for (int k = 0; k < 2; k++) com[k] += xy[k];
}
for (int k = 0; k < 2; k++) com[k] /= size;
for (int n = 0; n < size; n++) {
double dx = xyz.get(clusters[j][n])[0] - com[0];
double dy = xyz.get(clusters[j][n])[1] - com[1];
statsRadius.add(Math.sqrt(dx * dx + dy * dy));
}
}
}
if (showClusterMask) {
bp = new ByteProcessor(maskSize, maskSize);
for (int i = 0; i < mask.length; i++) if (mask[i] != 0)
bp.set(i, 128);
Utils.display(maskTitle, bp);
}
}
// Use the simulated cluster centres to create clusters of the desired size
if (clusterSimulation == 1 || clusterSimulation == 2) {
for (double[] clusterCentre : clusterCentres) {
int clusterN = (int) Math.round((clusterNumberSD > 0) ? dataGenerator.nextGaussian(clusterNumber, clusterNumberSD) : clusterNumber);
if (clusterN < 1)
continue;
//double[] clusterCentre = dist.next();
if (clusterN == 1) {
// No need for a cluster around a point
xyz.add(clusterCentre);
statsRadius.add(0);
statsSize.add(1);
} else {
// Generate N random points within a circle of the chosen cluster radius.
// Locate the centre-of-mass and the average distance to the centre.
double[] com = new double[3];
int j = 0;
while (j < clusterN) {
// Generate a random point within a circle uniformly
// http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
double t = 2.0 * Math.PI * randomGenerator.nextDouble();
double u = randomGenerator.nextDouble() + randomGenerator.nextDouble();
double r = clusterRadius * ((u > 1) ? 2 - u : u);
double x = r * Math.cos(t);
double y = r * Math.sin(t);
double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y };
xyz.add(xy);
for (int k = 0; k < 2; k++) com[k] += xy[k];
j++;
}
// Add the distance of the points from the centre of the cluster.
// Note this does not account for the movement due to precision.
statsSize.add(j);
if (j == 1) {
statsRadius.add(0);
} else {
for (int k = 0; k < 2; k++) com[k] /= j;
while (j > 0) {
double dx = xyz.get(xyz.size() - j)[0] - com[0];
double dy = xyz.get(xyz.size() - j)[1] - com[1];
statsRadius.add(Math.sqrt(dx * dx + dy * dy));
j--;
}
}
}
}
}
} else {
// Random distribution
for (int i = 0; i < nMolecules; i++) xyz.add(dist.next());
}
// The Gaussian sigma should be applied so the overall distance from the centre
// ( sqrt(x^2+y^2) ) has a standard deviation of sigmaS?
final double sigma1D = sigmaS / Math.sqrt(2);
// Show optional histograms
StoredDataStatistics intraDistances = null;
StoredData blinks = null;
if (showHistograms) {
int capacity = (int) (xyz.size() * blinkingRate);
intraDistances = new StoredDataStatistics(capacity);
blinks = new StoredData(capacity);
}
Statistics statsSigma = new Statistics();
for (int i = 0; i < xyz.size(); i++) {
int nOccurrences = getBlinks(dataGenerator, blinkingRate);
if (showHistograms)
blinks.add(nOccurrences);
final int size = molecules.size();
// Get coordinates in nm
final double[] moleculeXyz = xyz.get(i);
if (bp != null && nOccurrences > 0) {
bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255);
}
while (nOccurrences-- > 0) {
final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2);
// Add random precision
if (sigma1D > 0) {
final double dx = dataGenerator.nextGaussian(0, sigma1D);
final double dy = dataGenerator.nextGaussian(0, sigma1D);
localisationXy[0] += dx;
localisationXy[1] += dy;
if (!dist.isWithinXY(localisationXy))
continue;
// Calculate mean-squared displacement
statsSigma.add(dx * dx + dy * dy);
}
final double x = localisationXy[0];
final double y = localisationXy[1];
molecules.add(new Molecule(x, y, i, 1));
// Store in pixels
float[] params = new float[7];
params[Gaussian2DFunction.X_POSITION] = (float) (x / nmPerPixel);
params[Gaussian2DFunction.Y_POSITION] = (float) (y / nmPerPixel);
results.addf(i + 1, (int) x, (int) y, 0, 0, 0, params, null);
}
if (molecules.size() > size) {
count++;
if (showHistograms) {
int newCount = molecules.size() - size;
if (newCount == 1) {
//intraDistances.add(0);
continue;
}
// Get the distance matrix between these molecules
double[][] matrix = new double[newCount][newCount];
for (int ii = size, x = 0; ii < molecules.size(); ii++, x++) {
for (int jj = size + 1, y = 1; jj < molecules.size(); jj++, y++) {
final double d2 = molecules.get(ii).distance2(molecules.get(jj));
matrix[x][y] = matrix[y][x] = d2;
}
}
// Get the maximum distance for particle linkage clustering of this molecule
double max = 0;
for (int x = 0; x < newCount; x++) {
// Compare to all-other molecules and get the minimum distance
// needed to join at least one
double linkDistance = Double.POSITIVE_INFINITY;
for (int y = 0; y < newCount; y++) {
if (x == y)
continue;
if (matrix[x][y] < linkDistance)
linkDistance = matrix[x][y];
}
// Check if this is larger
if (max < linkDistance)
max = linkDistance;
}
intraDistances.add(Math.sqrt(max));
}
}
}
results.end();
if (bp != null)
Utils.display(maskTitle, bp);
// Used for debugging
//System.out.printf(" * Molecules = %d (%d activated)\n", xyz.size(), count);
//if (clusterSimulation > 0)
// System.out.printf(" * Cluster number = %s +/- %s. Radius = %s +/- %s\n",
// Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4),
// Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4));
log("Simulation results");
log(" * Molecules = %d (%d activated)", xyz.size(), count);
log(" * Blinking rate = %s", Utils.rounded((double) molecules.size() / xyz.size(), 4));
log(" * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? Utils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0");
if (showHistograms) {
if (intraDistances.getN() == 0) {
log(" * Mean Intra-Molecule particle linkage distance = 0 nm");
log(" * Fraction of inter-molecule particle linkage @ 0 nm = 0 %%");
} else {
plot(blinks, "Blinks/Molecule", true);
double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false);
// Determine 95th and 99th percentile
int p99 = intraHist[0].length - 1;
double limit1 = 0.99 * intraHist[1][p99];
double limit2 = 0.95 * intraHist[1][p99];
while (intraHist[1][p99] > limit1 && p99 > 0) p99--;
int p95 = p99;
while (intraHist[1][p95] > limit2 && p95 > 0) p95--;
log(" * Mean Intra-Molecule particle linkage distance = %s nm (95%% = %s, 99%% = %s, 100%% = %s)", Utils.rounded(intraDistances.getMean(), 4), Utils.rounded(intraHist[0][p95], 4), Utils.rounded(intraHist[0][p99], 4), Utils.rounded(intraHist[0][intraHist[0].length - 1], 4));
if (distanceAnalysis) {
performDistanceAnalysis(intraHist, p99);
}
}
}
if (clusterSimulation > 0) {
log(" * Cluster number = %s +/- %s", Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4));
log(" * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4));
}
}
use of gdsc.core.utils.Statistics in project GDSC-SMLM by aherbert.
the class CreateData method showSummary.
private double showSummary(List<? extends FluorophoreSequenceModel> fluorophores, List<LocalisationModel> localisations) {
IJ.showStatus("Calculating statistics ...");
createSummaryTable();
Statistics[] stats = new Statistics[NAMES.length];
for (int i = 0; i < stats.length; i++) {
stats[i] = (settings.showHistograms || alwaysRemoveOutliers[i]) ? new StoredDataStatistics() : new Statistics();
}
// Find the largest timepoint
ImagePlus outputImp = WindowManager.getImage(benchmarkImageId);
int nFrames;
if (outputImp == null) {
sortLocalisationsByTime(localisations);
nFrames = localisations.get(localisations.size() - 1).getTime();
} else {
nFrames = outputImp.getStackSize();
}
int[] countHistogram = new int[nFrames + 1];
// Use the localisations that were drawn to create the sampled on/off times
rebuildNeighbours(localisations);
// Assume that there is at least one localisation
LocalisationModel first = localisations.get(0);
// The current localisation
int currentId = first.getId();
// The last time this localisation was on
int lastT = first.getTime();
// Number of blinks
int blinks = 0;
// On-time of current pulse
int currentT = 0;
double signal = 0;
final double centreOffset = settings.size * 0.5;
// Used to convert the sampled times in frames into seconds
final double framesPerSecond = 1000.0 / settings.exposureTime;
final double gain = (settings.getTotalGain() > 0) ? settings.getTotalGain() : 1;
for (LocalisationModel l : localisations) {
if (l.getData() == null)
System.out.println("No localisation data. This should not happen!");
final double noise = (l.getData() != null) ? l.getData()[1] : 1;
final double intensity = (l.getData() != null) ? l.getData()[4] : l.getIntensity();
final double intensityInPhotons = intensity / gain;
// Q. What if the noise is zero, i.e. no background photon / read noise?
// Just ignore it at current.
final double snr = intensity / noise;
stats[SIGNAL].add(intensityInPhotons);
stats[NOISE].add(noise / gain);
if (noise != 0)
stats[SNR].add(snr);
//if (l.isContinuous())
if (l.getNext() != null && l.getPrevious() != null) {
stats[SIGNAL_CONTINUOUS].add(intensityInPhotons);
if (noise != 0)
stats[SNR_CONTINUOUS].add(snr);
}
int id = l.getId();
// Check if this a new fluorophore
if (currentId != id) {
// Add previous fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Reset
blinks = 0;
currentT = 1;
currentId = id;
signal = intensityInPhotons;
} else {
signal += intensityInPhotons;
// Check if the current fluorophore pulse is broken (i.e. a blink)
if (l.getTime() - 1 > lastT) {
blinks++;
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
currentT = 1;
stats[SAMPLED_T_OFF].add(((l.getTime() - 1) - lastT) / framesPerSecond);
} else {
// Continuous on-time
currentT++;
}
}
lastT = l.getTime();
countHistogram[lastT]++;
stats[X].add((l.getX() - centreOffset) * settings.pixelPitch);
stats[Y].add((l.getY() - centreOffset) * settings.pixelPitch);
stats[Z].add(l.getZ() * settings.pixelPitch);
}
// Final fluorophore
stats[SAMPLED_BLINKS].add(blinks);
stats[SAMPLED_T_ON].add(currentT / framesPerSecond);
stats[TOTAL_SIGNAL].add(signal);
// Samples per frame
for (int t = 1; t < countHistogram.length; t++) stats[SAMPLES].add(countHistogram[t]);
if (fluorophores != null) {
for (FluorophoreSequenceModel f : fluorophores) {
stats[BLINKS].add(f.getNumberOfBlinks());
// On-time
for (double t : f.getOnTimes()) stats[T_ON].add(t);
// Off-time
for (double t : f.getOffTimes()) stats[T_OFF].add(t);
}
} else {
// show no blinks
stats[BLINKS].add(0);
stats[T_ON].add(1);
//stats[T_OFF].add(0);
}
if (results != null) {
final boolean emCCD = (settings.getEmGain() > 1);
// Convert depth-of-field to pixels
final double depth = settings.depthOfField / settings.pixelPitch;
for (PeakResult r : results.getResults()) {
final double precision = r.getPrecision(settings.pixelPitch, gain, emCCD);
stats[PRECISION].add(precision);
// The error stores the z-depth in pixels
if (Math.abs(r.error) < depth)
stats[PRECISION_IN_FOCUS].add(precision);
stats[WIDTH].add(r.getSD());
}
// Compute density per frame. Multithread for speed
if (settings.densityRadius > 0) {
IJ.showStatus("Calculating density ...");
ExecutorService threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
List<Future<?>> futures = new LinkedList<Future<?>>();
final ArrayList<float[]> coords = new ArrayList<float[]>();
int t = results.getHead().getFrame();
final Statistics densityStats = stats[DENSITY];
final float radius = (float) (settings.densityRadius * getHWHM());
final Rectangle bounds = results.getBounds();
currentIndex = 0;
finalIndex = results.getTail().getFrame();
// Store the density for each result.
int[] allDensity = new int[results.size()];
int allIndex = 0;
for (PeakResult r : results.getResults()) {
if (t != r.getFrame()) {
allIndex += runDensityCalculation(threadPool, futures, coords, densityStats, radius, bounds, allDensity, allIndex);
}
coords.add(new float[] { r.getXPosition(), r.getYPosition() });
t = r.getFrame();
}
runDensityCalculation(threadPool, futures, coords, densityStats, radius, bounds, allDensity, allIndex);
Utils.waitForCompletion(futures);
threadPool.shutdownNow();
threadPool = null;
IJ.showProgress(1);
// Split results into singles (density = 0) and clustered (density > 0)
MemoryPeakResults singles = copyMemoryPeakResults("No Density");
MemoryPeakResults clustered = copyMemoryPeakResults("Density");
int i = 0;
for (PeakResult r : results.getResults()) {
// Store density in the original value field
r.origValue = allDensity[i];
if (allDensity[i++] == 0)
singles.add(r);
else
clustered.add(r);
}
}
}
StringBuilder sb = new StringBuilder();
sb.append(datasetNumber).append("\t");
sb.append((fluorophores == null) ? localisations.size() : fluorophores.size()).append("\t");
sb.append(stats[SAMPLED_BLINKS].getN() + (int) stats[SAMPLED_BLINKS].getSum()).append("\t");
sb.append(localisations.size()).append("\t");
sb.append(nFrames).append("\t");
sb.append(Utils.rounded(areaInUm)).append("\t");
sb.append(Utils.rounded(localisations.size() / (areaInUm * nFrames), 4)).append("\t");
sb.append(Utils.rounded(getHWHM(), 4)).append("\t");
double s = getPsfSD();
sb.append(Utils.rounded(s, 4)).append("\t");
s *= settings.pixelPitch;
final double sa = PSFCalculator.squarePixelAdjustment(s, settings.pixelPitch) / settings.pixelPitch;
sb.append(Utils.rounded(sa, 4)).append("\t");
// Width not valid for the Image PSF
int nStats = (imagePSF) ? stats.length - 1 : stats.length;
for (int i = 0; i < nStats; i++) {
double centre = (alwaysRemoveOutliers[i]) ? ((StoredDataStatistics) stats[i]).getStatistics().getPercentile(50) : stats[i].getMean();
sb.append(Utils.rounded(centre, 4)).append("\t");
}
if (java.awt.GraphicsEnvironment.isHeadless()) {
IJ.log(sb.toString());
return stats[SIGNAL].getMean();
} else {
summaryTable.append(sb.toString());
}
// Show histograms
if (settings.showHistograms) {
IJ.showStatus("Calculating histograms ...");
boolean[] chosenHistograms = getChoosenHistograms();
WindowOrganiser wo = new WindowOrganiser();
boolean requireRetile = false;
for (int i = 0; i < NAMES.length; i++) {
if (chosenHistograms[i]) {
wo.add(Utils.showHistogram(TITLE, (StoredDataStatistics) stats[i], NAMES[i], (integerDisplay[i]) ? 1 : 0, (settings.removeOutliers || alwaysRemoveOutliers[i]) ? 2 : 0, settings.histogramBins * ((integerDisplay[i]) ? 100 : 1)));
requireRetile = requireRetile || Utils.isNewWindow();
}
}
wo.tile();
}
IJ.showStatus("");
return stats[SIGNAL].getMean();
}
Aggregations