Examples with StoredDataStatistics gdsc.core.utils.StoredDataStatistics used on opensource projects

Example 11 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class BaseFunctionSolverTest method canFitSingleGaussianBetter.

void canFitSingleGaussianBetter(FunctionSolver solver, boolean applyBounds, FunctionSolver solver2, boolean applyBounds2, String name, String name2, NoiseModel noiseModel) {
    double[] noise = getNoise(noiseModel);
    if (solver.isWeighted())
        solver.setWeights(getWeights(noiseModel));
    int LOOPS = 5;
    randomGenerator.setSeed(seed);
    StoredDataStatistics[] stats = new StoredDataStatistics[6];
    String[] statName = { "Signal", "X", "Y" };
    int[] betterPrecision = new int[3];
    int[] totalPrecision = new int[3];
    int[] betterAccuracy = new int[3];
    int[] totalAccuracy = new int[3];
    int i1 = 0, i2 = 0;
    for (double s : signal) {
        double[] expected = createParams(1, s, 0, 0, 1);
        double[] lower = null, upper = null;
        if (applyBounds || applyBounds2) {
            lower = createParams(0, s * 0.5, -0.2, -0.2, 0.8);
            upper = createParams(3, s * 2, 0.2, 0.2, 1.2);
        }
        if (applyBounds)
            solver.setBounds(lower, upper);
        if (applyBounds2)
            solver2.setBounds(lower, upper);
        for (int loop = LOOPS; loop-- > 0; ) {
            double[] data = drawGaussian(expected, noise, noiseModel);
            for (int i = 0; i < stats.length; i++) stats[i] = new StoredDataStatistics();
            for (double db : base) for (double dx : shift) for (double dy : shift) for (double dsx : factor) {
                double[] p = createParams(db, s, dx, dy, dsx);
                double[] fp = fitGaussian(solver, data, p, expected);
                i1 += solver.getEvaluations();
                double[] fp2 = fitGaussian(solver2, data, p, expected);
                i2 += solver2.getEvaluations();
                // Get the mean and sd (the fit precision)
                compare(fp, expected, fp2, expected, Gaussian2DFunction.SIGNAL, stats[0], stats[1]);
                compare(fp, expected, fp2, expected, Gaussian2DFunction.X_POSITION, stats[2], stats[3]);
                compare(fp, expected, fp2, expected, Gaussian2DFunction.Y_POSITION, stats[4], stats[5]);
            // Use the distance
            //stats[2].add(distance(fp, expected));
            //stats[3].add(distance(fp2, expected2));
            }
            // two sided
            double alpha = 0.05;
            for (int i = 0; i < stats.length; i += 2) {
                double u1 = stats[i].getMean();
                double u2 = stats[i + 1].getMean();
                double sd1 = stats[i].getStandardDeviation();
                double sd2 = stats[i + 1].getStandardDeviation();
                TTest tt = new TTest();
                boolean diff = tt.tTest(stats[i].getValues(), stats[i + 1].getValues(), alpha);
                int index = i / 2;
                String msg = String.format("%s vs %s : %.1f (%s) %s %f +/- %f vs %f +/- %f  (N=%d) %b", name2, name, s, noiseModel, statName[index], u2, sd2, u1, sd1, stats[i].getN(), diff);
                if (diff) {
                    // Different means. Check they are roughly the same
                    if (DoubleEquality.almostEqualRelativeOrAbsolute(u1, u2, 0.1, 0)) {
                        // Basically the same. Check which is more precise
                        if (!DoubleEquality.almostEqualRelativeOrAbsolute(sd1, sd2, 0.05, 0)) {
                            if (sd2 < sd1) {
                                betterPrecision[index]++;
                                println(msg + " P*");
                            } else
                                println(msg + " P");
                            totalPrecision[index]++;
                        }
                    } else {
                        // Check which is more accurate (closer to zero)
                        u1 = Math.abs(u1);
                        u2 = Math.abs(u2);
                        if (u2 < u1) {
                            betterAccuracy[index]++;
                            println(msg + " A*");
                        } else
                            println(msg + " A");
                        totalAccuracy[index]++;
                    }
                } else {
                    // The same means. Check that it is more precise
                    if (!DoubleEquality.almostEqualRelativeOrAbsolute(sd1, sd2, 0.05, 0)) {
                        if (sd2 < sd1) {
                            betterPrecision[index]++;
                            println(msg + " P*");
                        } else
                            println(msg + " P");
                        totalPrecision[index]++;
                    }
                }
            }
        }
    }
    int better = 0, total = 0;
    for (int index = 0; index < statName.length; index++) {
        better += betterPrecision[index] + betterAccuracy[index];
        total += totalPrecision[index] + totalAccuracy[index];
        test(name2, name, statName[index] + " P", betterPrecision[index], totalPrecision[index], printBetterDetails);
        test(name2, name, statName[index] + " A", betterAccuracy[index], totalAccuracy[index], printBetterDetails);
    }
    test(name2, name, String.format("All (eval [%d] [%d]) : ", i2, i1), better, total, true);
}

Example 12 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class PoissonFunctionTest method cumulativeProbabilityIsOneWithInteger.

private int[] cumulativeProbabilityIsOneWithInteger(final double gain, final double mu) {
    final double o = mu * gain;
    PoissonFunction f = new PoissonFunction(1.0 / gain, false);
    double p = 0;
    int x = 0;
    StoredDataStatistics stats = new StoredDataStatistics();
    double maxp = 0;
    int maxc = 0;
    // At large mu it is approximately normal so use 3 sqrt(mu) for the range added to the mean
    if (mu > 0) {
        int max = (int) Math.ceil(o + 3 * Math.sqrt(o));
        for (; x <= max; x++) {
            final double pp = f.likelihood(x, o);
            //System.out.printf("x=%d, p=%f\n", x, pp);
            p += pp;
            stats.add(p);
            if (maxp < pp) {
                maxp = pp;
                maxc = x;
            }
        }
        if (p > 1.01)
            Assert.fail("P > 1: " + p);
    }
    // We have most of the probability density. 
    // Now keep evaluating up and down until no difference
    final double changeTolerance = 1e-6;
    for (; ; x++) {
        final double pp = f.likelihood(x, o);
        //System.out.printf("x=%d, p=%f\n", x, pp);
        p += pp;
        stats.add(p);
        if (maxp < pp) {
            maxp = pp;
            maxc = x;
        }
        if (pp / p < changeTolerance)
            break;
    }
    // Find the range for 99.5% of the sum
    double[] h = stats.getValues();
    int minx = 0, maxx = h.length - 1;
    while (h[minx + 1] < 0.0025) minx++;
    while (h[maxx - 1] > 0.9975) maxx--;
    System.out.printf("g=%f, mu=%f, o=%f, p=%f, min=%d, %f @ %d, max=%d\n", gain, mu, o, p, minx, maxp, maxc, maxx);
    Assert.assertEquals(String.format("g=%f, mu=%f", gain, mu), 1, p, 0.02);
    return new int[] { minx, maxx };
}

Example 13 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics in project GDSC-SMLM by aherbert.

the class PoissonGammaGaussianFunctionTest method fasterThan.

private void fasterThan(PoissonGammaGaussianFunction f1, PoissonGammaGaussianFunction f2) {
    // Generate realistic data from the probability mass function
    double[][] samples = new double[photons.length][];
    for (int j = 0; j < photons.length; j++) {
        int start = (int) (4 * -s);
        int u = start;
        StoredDataStatistics stats = new StoredDataStatistics();
        while (stats.getSum() < 0.995) {
            stats.add(f1.likelihood(u, photons[j]));
            u++;
        }
        // Generate cumulative probability
        double[] data = stats.getValues();
        for (int i = 1; i < data.length; i++) data[i] += data[i - 1];
        // Sample
        RandomGenerator rand = new Well19937c();
        double[] sample = new double[1000];
        for (int i = 0; i < sample.length; i++) {
            final double p = rand.nextDouble();
            int x = 0;
            while (x < data.length && data[x] < p) x++;
            sample[i] = x;
        }
        samples[j] = sample;
    }
    // Warm-up
    run(f1, samples, photons);
    run(f2, samples, photons);
    long t1 = 0;
    for (int i = 0; i < 5; i++) t1 += run(f1, samples, photons);
    long t2 = 0;
    for (int i = 0; i < 5; i++) t2 += run(f2, samples, photons);
    System.out.printf("%s  %d -> %s  %d = %f x\n", getName(f1), t1, getName(f2), t2, (double) t1 / t2);
}

Example 14 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics 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("");
}

Example 15 with StoredDataStatistics

use of gdsc.core.utils.StoredDataStatistics 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));
    }
}