Examples with NormalizedGaussianSampler org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler used on opensource projects

Example 11 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler in project GDSC-SMLM by aherbert.

the class DynamicMultipleTargetTracingTest method testTraceMolecules.

/**
 * Test trace molecules using 2 molecules. One is fixed and the other moves across it. The tracing
 * should assign the fixed molecule correctly as it has a low local diffusion rate and different
 * intensity.
 */
@Test
void testTraceMolecules() {
    // The test is not very robust and fails 10% of the time. A fixed seed corrects this.
    final UniformRandomProvider rng = RngUtils.create(0x12345L);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    // localisation precision (in pixels)
    final double s = 0.1;
    final SharedStateContinuousSampler intensity1 = SamplerUtils.createGaussianSampler(rng, 1000, 100);
    final SharedStateContinuousSampler intensity2 = SamplerUtils.createGaussianSampler(rng, 500, 50);
    final MemoryPeakResults results = new MemoryPeakResults(100);
    final CalibrationWriter writer = results.getCalibrationWriterSafe();
    // 0.1 um pixels, 1 second exposure time
    writer.setDistanceUnit(DistanceUnit.PIXEL);
    writer.setNmPerPixel(100);
    writer.setExposureTime(1000);
    results.setCalibration(writer.getCalibration());
    // First molecule diffuses roughly across the field from top-left to bottom-right.
    // 5 frames is the default for local stats, 15 frames for trajectory removal.
    // Use 20 so we build local stats and can expire a trajectory.
    final int size = 20;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (i + gauss.sample() * s), (float) (i + gauss.sample() * s), (float) (intensity1.sample())));
    }
    // Second molecule is fixed in the centre with a lower intensity (allow
    // correct matching when tracks overlap)
    final int x = size / 2;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (x + gauss.sample() * s), (float) (x + gauss.sample() * s), (float) (intensity2.sample())));
    }
    // Add a single molecule that will not connect to anything in the second frame.
    // This should create a trajectory that will expire.
    results.add(new PeakResult(1, size, size, (float) (intensity1.sample())));
    // 1 diffuses top-left to bottom-right.
    // 2 is fixed in the centre.
    // 3 is in the bottom-right for 1 frame.
    // 
    // 1
    // 1
    // 1
    // 12
    // 1
    // 1
    // 13
    // 
    // Molecule 3 can sometimes connect to the long lifetime molecules once they have been alive
    // long enough to create a local probability model. The default lifetime is 5 frames.
    // Setting this to 10 frames allows a better local model to be created.
    // Move centre to centre each jump => sqrt(2 * 0.1^2) = 0.141 um or 0.02 um^2
    // MSD = 4D => D = 0.02 / 4 = 0.005
    final DmttConfiguration config = DmttConfiguration.newBuilder(0.005).setTemporalWindow(10).build();
    final List<Trace> traces = new DynamicMultipleTargetTracing(results).traceMolecules(config);
    // Should have 3 traces
    Assertions.assertEquals(3, traces.size());
    // Assert ids start from 1
    for (int i = 0; i < traces.size(); i++) {
        Assertions.assertEquals(i + 1, traces.get(i).getId());
    }
    // Traces should be 2 full length and 1 single peak
    Assertions.assertEquals(size, traces.get(0).size());
    Assertions.assertEquals(size, traces.get(1).size());
    Assertions.assertEquals(1, traces.get(2).size());
    // Do an analysis on the actual tracks.
    // One should be based in the centre and the other should have parts close to position (i,i)
    // for each frame i.
    final PeakResult[] peaks = results.toArray();
    // Assume traces are initially created using the input order of the results.
    final Trace t1 = traces.get(0);
    final Trace t2 = traces.get(1);
    for (int i = 0; i < size; i++) {
        Assertions.assertSame(peaks[i], t1.get(i));
        Assertions.assertSame(peaks[i + size], t2.get(i));
    }
}

Example 12 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    pluginSettings = Settings.load();
    pluginSettings.save();
    if (IJ.controlKeyDown()) {
        simpleTest();
        return;
    }
    extraOptions = ImageJUtils.isExtraOptions();
    if (!showDialog()) {
        return;
    }
    lastSimulation.set(null);
    final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond());
    conversionFactor = 1000000.0 / (settings.getPixelPitch() * settings.getPixelPitch());
    // Diffusion rate is um^2/sec. Convert to pixels per simulation frame.
    final double diffusionRateInPixelsPerSecond = settings.getDiffusionRate() * conversionFactor;
    final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.getStepsPerSecond();
    final double precisionInPixels = myPrecision / settings.getPixelPitch();
    final boolean addError = myPrecision != 0;
    ImageJUtils.log(TITLE + " : D = %s um^2/sec, Precision = %s nm", MathUtils.rounded(settings.getDiffusionRate(), 4), MathUtils.rounded(myPrecision, 4));
    ImageJUtils.log("Mean-displacement per dimension = %s nm/sec", MathUtils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.getDiffusionRate()), 4));
    if (extraOptions) {
        ImageJUtils.log("Step size = %s, precision = %s", MathUtils.rounded(ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep)), MathUtils.rounded(precisionInPixels));
    }
    // Convert diffusion co-efficient into the standard deviation for the random walk
    final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
    final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep);
    ImageJUtils.log("Simulation step-size = %s nm", MathUtils.rounded(settings.getPixelPitch() * diffusionSigma, 4));
    // Move the molecules and get the diffusion rate
    IJ.showStatus("Simulating ...");
    final long start = System.nanoTime();
    final UniformRandomProvider random = UniformRandomProviders.create();
    final Statistics[] stats2D = new Statistics[totalSteps];
    final Statistics[] stats3D = new Statistics[totalSteps];
    final StoredDataStatistics jumpDistances2D = new StoredDataStatistics(totalSteps);
    final StoredDataStatistics jumpDistances3D = new StoredDataStatistics(totalSteps);
    for (int j = 0; j < totalSteps; j++) {
        stats2D[j] = new Statistics();
        stats3D[j] = new Statistics();
    }
    final SphericalDistribution dist = new SphericalDistribution(settings.getConfinementRadius() / settings.getPixelPitch());
    final Statistics asymptote = new Statistics();
    // Save results to memory
    final MemoryPeakResults results = new MemoryPeakResults(totalSteps);
    results.setCalibration(CalibrationHelper.create(settings.getPixelPitch(), 1, 1000.0 / settings.getStepsPerSecond()));
    results.setName(TITLE);
    results.setPsf(PsfHelper.create(PSFType.CUSTOM));
    int peak = 0;
    // Store raw coordinates
    final ArrayList<Point> points = new ArrayList<>(totalSteps);
    final StoredData totalJumpDistances1D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances2D = new StoredData(settings.getParticles());
    final StoredData totalJumpDistances3D = new StoredData(settings.getParticles());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    for (int i = 0; i < settings.getParticles(); i++) {
        if (i % 16 == 0) {
            IJ.showProgress(i, settings.getParticles());
            if (ImageJUtils.isInterrupted()) {
                return;
            }
        }
        // Increment the frame so that tracing analysis can distinguish traces
        peak++;
        double[] origin = new double[3];
        final int id = i + 1;
        final MoleculeModel m = new MoleculeModel(id, origin.clone());
        if (addError) {
            origin = addError(origin, precisionInPixels, gauss);
        }
        if (pluginSettings.useConfinement) {
            // Note: When using confinement the average displacement should asymptote
            // at the average distance of a point from the centre of a ball. This is 3r/4.
            // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM
            // The equivalent in 2D is 2r/3. However although we are plotting 2D distance
            // this is a projection of the 3D position onto the plane and so the particles
            // will not be evenly spread (there will be clustering at centre caused by the
            // poles)
            final double[] axis = (diffusionType == DiffusionType.LINEAR_WALK) ? nextVector(gauss) : null;
            for (int j = 0; j < totalSteps; j++) {
                double[] xyz = m.getCoordinates();
                final double[] originalXyz = xyz.clone();
                for (int n = pluginSettings.confinementAttempts; n-- > 0; ) {
                    if (diffusionType == DiffusionType.GRID_WALK) {
                        m.walk(diffusionSigma, random);
                    } else if (diffusionType == DiffusionType.LINEAR_WALK) {
                        m.slide(diffusionSigma, axis, random);
                    } else {
                        m.move(diffusionSigma, random);
                    }
                    if (!dist.isWithin(m.getCoordinates())) {
                        // Reset position
                        for (int k = 0; k < 3; k++) {
                            xyz[k] = originalXyz[k];
                        }
                    } else {
                        // The move was allowed
                        break;
                    }
                }
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
            asymptote.add(distance(m.getCoordinates()));
        } else if (diffusionType == DiffusionType.GRID_WALK) {
            for (int j = 0; j < totalSteps; j++) {
                m.walk(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else if (diffusionType == DiffusionType.LINEAR_WALK) {
            final double[] axis = nextVector(gauss);
            for (int j = 0; j < totalSteps; j++) {
                m.slide(diffusionSigma, axis, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        } else {
            for (int j = 0; j < totalSteps; j++) {
                m.move(diffusionSigma, random);
                double[] xyz = m.getCoordinates();
                points.add(new Point(id, xyz));
                if (addError) {
                    xyz = addError(xyz, precisionInPixels, gauss);
                }
                peak = record(xyz, id, peak, stats2D[j], stats3D[j], jumpDistances2D, jumpDistances3D, origin, results);
            }
        }
        // Debug: record all the particles so they can be analysed
        // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ());
        final double[] xyz = m.getCoordinates();
        double d2 = 0;
        totalJumpDistances1D.add(d2 = xyz[0] * xyz[0]);
        totalJumpDistances2D.add(d2 += xyz[1] * xyz[1]);
        totalJumpDistances3D.add(d2 += xyz[2] * xyz[2]);
    }
    final long nanoseconds = System.nanoTime() - start;
    IJ.showProgress(1);
    MemoryPeakResults.addResults(results);
    simulation = new SimulationData(results.getName(), myPrecision);
    // Convert pixels^2/step to um^2/sec
    final double msd2D = (jumpDistances2D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    final double msd3D = (jumpDistances3D.getMean() / conversionFactor) / (results.getCalibrationReader().getExposureTime() / 1000);
    ImageJUtils.log("Raw data D=%s um^2/s, Precision = %s nm, N=%d, step=%s s, mean2D=%s um^2, " + "MSD 2D = %s um^2/s, mean3D=%s um^2, MSD 3D = %s um^2/s", MathUtils.rounded(settings.getDiffusionRate()), MathUtils.rounded(myPrecision), jumpDistances2D.getN(), MathUtils.rounded(results.getCalibrationReader().getExposureTime() / 1000), MathUtils.rounded(jumpDistances2D.getMean() / conversionFactor), MathUtils.rounded(msd2D), MathUtils.rounded(jumpDistances3D.getMean() / conversionFactor), MathUtils.rounded(msd3D));
    aggregateIntoFrames(points, addError, precisionInPixels, gauss);
    IJ.showStatus("Analysing results ...");
    if (pluginSettings.showDiffusionExample) {
        showExample(totalSteps, diffusionSigma, random);
    }
    // Plot a graph of mean squared distance
    final double[] xValues = new double[stats2D.length];
    final double[] yValues2D = new double[stats2D.length];
    final double[] yValues3D = new double[stats3D.length];
    final double[] upper2D = new double[stats2D.length];
    final double[] lower2D = new double[stats2D.length];
    final double[] upper3D = new double[stats3D.length];
    final double[] lower3D = new double[stats3D.length];
    final SimpleRegression r2D = new SimpleRegression(false);
    final SimpleRegression r3D = new SimpleRegression(false);
    final int firstN = (pluginSettings.useConfinement) ? pluginSettings.fitN : totalSteps;
    for (int j = 0; j < totalSteps; j++) {
        // Convert steps to seconds
        xValues[j] = (j + 1) / settings.getStepsPerSecond();
        // Convert values in pixels^2 to um^2
        final double mean2D = stats2D[j].getMean() / conversionFactor;
        final double mean3D = stats3D[j].getMean() / conversionFactor;
        final double sd2D = stats2D[j].getStandardDeviation() / conversionFactor;
        final double sd3D = stats3D[j].getStandardDeviation() / conversionFactor;
        yValues2D[j] = mean2D;
        yValues3D[j] = mean3D;
        upper2D[j] = mean2D + sd2D;
        lower2D[j] = mean2D - sd2D;
        upper3D[j] = mean3D + sd3D;
        lower3D[j] = mean3D - sd3D;
        if (j < firstN) {
            r2D.addData(xValues[j], yValues2D[j]);
            r3D.addData(xValues[j], yValues3D[j]);
        }
    }
    // TODO - Fit using the equation for 2D confined diffusion:
    // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2)
    // s = localisation precision
    // R = confinement radius
    // D = 2D diffusion coefficient
    // t = time
    final PolynomialFunction fitted2D;
    final PolynomialFunction fitted3D;
    if (r2D.getN() > 0) {
        // Do linear regression to get diffusion rate
        final double[] best2D = new double[] { r2D.getIntercept(), r2D.getSlope() };
        fitted2D = new PolynomialFunction(best2D);
        final double[] best3D = new double[] { r3D.getIntercept(), r3D.getSlope() };
        fitted3D = new PolynomialFunction(best3D);
        // For 2D diffusion: d^2 = 4D
        // where: d^2 = mean-square displacement
        double diffCoeff = best2D[1] / 4.0;
        final String msg = "2D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")";
        IJ.showStatus(msg);
        ImageJUtils.log(msg);
        diffCoeff = best3D[1] / 6.0;
        ImageJUtils.log("3D Diffusion rate = " + MathUtils.rounded(diffCoeff, 4) + " um^2 / sec (" + TextUtils.nanosToString(nanoseconds) + ")");
    } else {
        fitted2D = fitted3D = null;
    }
    // Create plots
    plotMsd(totalSteps, xValues, yValues2D, lower2D, upper2D, fitted2D, 2);
    plotMsd(totalSteps, xValues, yValues3D, lower3D, upper3D, fitted3D, 3);
    plotJumpDistances(TITLE, jumpDistances2D, 2, 1);
    plotJumpDistances(TITLE, jumpDistances3D, 3, 1);
    // Show the total jump length for debugging
    // plotJumpDistances(TITLE + " total", totalJumpDistances1D, 1, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances2D, 2, totalSteps);
    // plotJumpDistances(TITLE + " total", totalJumpDistances3D, 3, totalSteps);
    windowOrganiser.tile();
    if (pluginSettings.useConfinement) {
        ImageJUtils.log("3D asymptote distance = %s nm (expected %.2f)", MathUtils.rounded(asymptote.getMean() * settings.getPixelPitch(), 4), 3 * settings.getConfinementRadius() / 4);
    }
}

Example 13 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler in project GDSC-SMLM by aherbert.

the class DiffusionRateTest method msdAnalysis.

/**
 * Tabulate the observed MSD for different jump distances.
 *
 * @param points the points
 */
private void msdAnalysis(ArrayList<Point> points) {
    if (myMsdAnalysisSteps == 0) {
        return;
    }
    IJ.showStatus("MSD analysis ...");
    IJ.showProgress(1, myMsdAnalysisSteps);
    // This will only be fast if the list is an array
    final Point[] list = points.toArray(new Point[0]);
    // Compute the base MSD
    final Point origin = new Point(0, 0, 0);
    double sum = origin.distance2(list[0]);
    int count = 1;
    for (int i = 1; i < list.length; i++) {
        final Point last = list[i - 1];
        final Point current = list[i];
        if (last.id == current.id) {
            sum += last.distance2(current);
        } else {
            sum += origin.distance2(current);
        }
        count++;
    }
    // Create a new set of points that have coordinates that
    // are the rolling average over the number of aggregate steps
    final DoubleRollingArray x = new DoubleRollingArray(pluginSettings.aggregateSteps);
    final DoubleRollingArray y = new DoubleRollingArray(pluginSettings.aggregateSteps);
    int id = 0;
    int length = 0;
    for (final Point p : points) {
        if (p.id != id) {
            x.clear();
            y.clear();
        }
        id = p.id;
        x.add(p.x);
        y.add(p.y);
        // Only create a point if the full aggregation size is reached
        if (x.isFull()) {
            list[length++] = new Point(id, x.getAverage(), y.getAverage());
        }
    }
    // Q - is this useful?
    final double p = myPrecision / settings.getPixelPitch();
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    final int totalSteps = (int) Math.ceil(settings.getSeconds() * settings.getStepsPerSecond() - pluginSettings.aggregateSteps);
    final int limit = Math.min(totalSteps, myMsdAnalysisSteps);
    final Ticker ticker = ImageJUtils.createTicker(limit, 1);
    final TextWindow msdTable = createMsdTable((sum / count) * settings.getStepsPerSecond() / conversionFactor);
    try (BufferedTextWindow bw = new BufferedTextWindow(msdTable)) {
        bw.setIncrement(0);
        for (int step = 1; step <= myMsdAnalysisSteps; step++) {
            sum = 0;
            count = 0;
            for (int i = step; i < length; i++) {
                final Point last = list[i - step];
                final Point current = list[i];
                if (last.id == current.id) {
                    if (p == 0) {
                        sum += last.distance2(current);
                        count++;
                    } else {
                        // is the same if enough samples are present
                        for (int ii = 1; ii-- > 0; ) {
                            sum += last.distance2(current, p, gauss);
                            count++;
                        }
                    }
                }
            }
            if (count == 0) {
                break;
            }
            bw.append(addResult(step, sum, count));
            ticker.tick();
        }
    }
    IJ.showProgress(1);
}

Example 14 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method simulateFromPoissonGammaGaussian.

/**
 * Randomly generate a histogram from poisson-gamma-gaussian samples.
 *
 * @return The histogram
 */
private int[] simulateFromPoissonGammaGaussian() {
    // Randomly sample
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final PoissonSampler poisson = new PoissonSampler(rng, settings.settingPhotons);
    final MarsagliaTsangGammaSampler gamma = new MarsagliaTsangGammaSampler(rng, settings.settingPhotons, settings.settingGain);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    final int steps = settings.simulationSize;
    final int[] samples = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0) {
            IJ.showProgress(n, steps);
        }
        // Poisson
        double sample = poisson.sample();
        // Gamma
        if (sample > 0) {
            gamma.setAlpha(sample);
            sample = gamma.sample();
        }
        // Gaussian
        sample += settings.settingNoise * gauss.sample();
        // Convert the sample to a count
        samples[n] = (int) Math.round(sample + settings.settingBias);
    }
    final int max = MathUtils.max(samples);
    final int[] histogram = new int[max + 1];
    for (final int s : samples) {
        histogram[s]++;
    }
    return histogram;
}

Example 15 with NormalizedGaussianSampler

use of org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler 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 (settings.blinkingDistribution == 3) {
        log("  - Clusters = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Molecules/cluster = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
        log("  - p-Value = %s", MathUtils.rounded(settings.pvalue, 4));
    } else {
        log("  - Molecules = %d", settings.numberOfMolecules);
        log("  - Simulation size = %s um", MathUtils.rounded(settings.simulationSize, 4));
        log("  - Blinking rate = %s", MathUtils.rounded(settings.blinkingRate, 4));
        log("  - Blinking distribution = %s", Settings.BLINKING_DISTRIBUTION[settings.blinkingDistribution]);
    }
    log("  - Average precision = %s nm", MathUtils.rounded(settings.sigmaS, 4));
    log("  - Clusters simulation = " + Settings.CLUSTER_SIMULATION[settings.clusterSimulation]);
    if (settings.clusterSimulation > 0) {
        log("  - Cluster number = %s +/- %s", MathUtils.rounded(settings.clusterNumber, 4), MathUtils.rounded(settings.clusterNumberStdDev, 4));
        log("  - Cluster radius = %s nm", MathUtils.rounded(settings.clusterRadius, 4));
    }
    final double nmPerPixel = 100;
    final double width = settings.simulationSize * 1000.0;
    final UniformRandomProvider rng = UniformRandomProviders.create();
    final UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, rng.nextInt());
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    settings.molecules = new ArrayList<>(settings.numberOfMolecules);
    // Create some dummy results since the calibration is required for later analysis
    settings.results = new MemoryPeakResults(PsfHelper.create(PSFType.CUSTOM));
    settings.results.setCalibration(CalibrationHelper.create(nmPerPixel, 1, 100));
    settings.results.setSource(new NullSource("Molecule Simulation"));
    settings.results.begin();
    int count = 0;
    // Generate a sequence of coordinates
    final ArrayList<double[]> xyz = new ArrayList<>((int) (settings.numberOfMolecules * 1.1));
    final Statistics statsRadius = new Statistics();
    final Statistics statsSize = new Statistics();
    final String maskTitle = TITLE + " Cluster Mask";
    ByteProcessor bp = null;
    double maskScale = 0;
    if (settings.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 = settings.lowResolutionImageSize;
        int[] mask = null;
        // scale is in nm/pixel
        maskScale = width / maskSize;
        final ArrayList<double[]> clusterCentres = new ArrayList<>();
        int totalSteps = 1 + (int) Math.ceil(settings.numberOfMolecules / settings.clusterNumber);
        if (settings.clusterSimulation == 2 || settings.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, rng);
            double[] centre;
            IJ.showStatus("Computing clusters mask");
            final int roiRadius = (int) Math.round((settings.clusterRadius * 2) / maskScale);
            if (settings.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 * MathUtils.pow2(settings.clusterRadius / maskScale)));
            }
            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
                final double cx = centre[0] / maskScale;
                final double cy = centre[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0);
                try {
                    maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                } catch (final IllegalArgumentException ex) {
                    // 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;
                }
            }
            ImageJUtils.finished();
        } else {
            // Pick centres randomly from the distribution
            while (clusterCentres.size() < totalSteps) {
                clusterCentres.add(dist.next());
            }
        }
        final double scaledRadius = settings.clusterRadius / maskScale;
        if (settings.showClusterMask || settings.clusterSimulation == 3) {
            // Show the mask for the clusters
            if (mask == null) {
                mask = new int[maskSize * maskSize];
            } else {
                Arrays.fill(mask, 0);
            }
            final int roiRadius = (int) Math.round(scaledRadius);
            for (final double[] c : clusterCentres) {
                final double cx = c[0] / maskScale;
                final double cy = c[1] / maskScale;
                fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1);
            }
            if (settings.clusterSimulation == 3) {
                // We have the mask. Now pick points at random from the mask.
                final MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, rng);
                // Allocate each molecule position to a parent circle so defining clusters.
                final int[][] clusters = new int[clusterCentres.size()][];
                final int[] clusterSize = new int[clusters.length];
                for (int i = 0; i < settings.numberOfMolecules; i++) {
                    final 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++) {
                        final 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
                    final double[] com = new double[2];
                    for (int n = 0; n < size; n++) {
                        final 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++) {
                        final double dx = xyz.get(clusters[j][n])[0] - com[0];
                        final double dy = xyz.get(clusters[j][n])[1] - com[1];
                        statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                    }
                }
            }
            if (settings.showClusterMask) {
                bp = new ByteProcessor(maskSize, maskSize);
                for (int i = 0; i < mask.length; i++) {
                    if (mask[i] != 0) {
                        bp.set(i, 128);
                    }
                }
                ImageJUtils.display(maskTitle, bp);
            }
        }
        // Use the simulated cluster centres to create clusters of the desired size
        if (settings.clusterSimulation == 1 || settings.clusterSimulation == 2) {
            for (final double[] clusterCentre : clusterCentres) {
                final int clusterN = (int) Math.round((settings.clusterNumberStdDev > 0) ? settings.clusterNumber + gauss.sample() * settings.clusterNumberStdDev : settings.clusterNumber);
                if (clusterN < 1) {
                    continue;
                }
                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.
                    final double[] com = new double[3];
                    int size = 0;
                    while (size < clusterN) {
                        // Generate a random point within a circle uniformly
                        // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly
                        final double t = 2.0 * Math.PI * rng.nextDouble();
                        final double u = rng.nextDouble() + rng.nextDouble();
                        final double r = settings.clusterRadius * ((u > 1) ? 2 - u : u);
                        final double x = r * Math.cos(t);
                        final double y = r * Math.sin(t);
                        final double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y };
                        xyz.add(xy);
                        for (int k = 0; k < 2; k++) {
                            com[k] += xy[k];
                        }
                        size++;
                    }
                    // 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(size);
                    if (size == 1) {
                        statsRadius.add(0);
                    } else {
                        for (int k = 0; k < 2; k++) {
                            com[k] /= size;
                        }
                        while (size > 0) {
                            final double dx = xyz.get(xyz.size() - size)[0] - com[0];
                            final double dy = xyz.get(xyz.size() - size)[1] - com[1];
                            statsRadius.add(Math.sqrt(dx * dx + dy * dy));
                            size--;
                        }
                    }
                }
            }
        }
    } else {
        // Random distribution
        for (int i = 0; i < settings.numberOfMolecules; 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 = settings.sigmaS / Math.sqrt(2);
    // Show optional histograms
    StoredDataStatistics intraDistances = null;
    StoredData blinks = null;
    if (settings.showHistograms) {
        final int capacity = (int) (xyz.size() * settings.blinkingRate);
        intraDistances = new StoredDataStatistics(capacity);
        blinks = new StoredData(capacity);
    }
    final Statistics statsSigma = new Statistics();
    for (int i = 0; i < xyz.size(); i++) {
        int occurrences = getBlinks(rng, settings.blinkingRate);
        if (blinks != null) {
            blinks.add(occurrences);
        }
        final int size = settings.molecules.size();
        // Get coordinates in nm
        final double[] moleculeXyz = xyz.get(i);
        if (bp != null && occurrences > 0) {
            bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255);
        }
        while (occurrences-- > 0) {
            final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2);
            // Add random precision
            if (sigma1D > 0) {
                final double dx = gauss.sample() * sigma1D;
                final double dy = gauss.sample() * 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];
            settings.molecules.add(new Molecule(x, y, i, 1));
            // Store in pixels
            final float xx = (float) (x / nmPerPixel);
            final float yy = (float) (y / nmPerPixel);
            final float[] params = PeakResult.createParams(0, 0, xx, yy, 0);
            settings.results.add(i + 1, (int) xx, (int) yy, 0, 0, 0, 0, params, null);
        }
        if (settings.molecules.size() > size) {
            count++;
            if (intraDistances != null) {
                final int newCount = settings.molecules.size() - size;
                if (newCount == 1) {
                    // No intra-molecule distances
                    continue;
                }
                // Get the distance matrix between these molecules
                final double[][] matrix = new double[newCount][newCount];
                for (int ii = size, x = 0; ii < settings.molecules.size(); ii++, x++) {
                    for (int jj = size + 1, y = 1; jj < settings.molecules.size(); jj++, y++) {
                        final double d2 = settings.molecules.get(ii).distance2(settings.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));
            }
        }
    }
    settings.results.end();
    if (bp != null) {
        final ImagePlus imp = ImageJUtils.display(maskTitle, bp);
        final Calibration cal = imp.getCalibration();
        cal.setUnit("nm");
        cal.pixelWidth = cal.pixelHeight = maskScale;
    }
    log("Simulation results");
    log("  * Molecules = %d (%d activated)", xyz.size(), count);
    log("  * Blinking rate = %s", MathUtils.rounded((double) settings.molecules.size() / xyz.size(), 4));
    log("  * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? MathUtils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0");
    if (intraDistances != null) {
        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);
            final double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false);
            // Determine 95th and 99th percentile
            // Will not be null as we requested a non-integer histogram.
            int p99 = intraHist[0].length - 1;
            final double limit1 = 0.99 * intraHist[1][p99];
            final 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)", MathUtils.rounded(intraDistances.getMean(), 4), MathUtils.rounded(intraHist[0][p95], 4), MathUtils.rounded(intraHist[0][p99], 4), MathUtils.rounded(intraHist[0][intraHist[0].length - 1], 4));
            if (settings.distanceAnalysis) {
                performDistanceAnalysis(intraHist, p99);
            }
        }
    }
    if (settings.clusterSimulation > 0) {
        log("  * Cluster number = %s +/- %s", MathUtils.rounded(statsSize.getMean(), 4), MathUtils.rounded(statsSize.getStandardDeviation(), 4));
        log("  * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", MathUtils.rounded(statsRadius.getMean(), 4), MathUtils.rounded(statsRadius.getStandardDeviation(), 4));
    }
}