Search in sources :

Example 6 with NormalizedGaussianSampler

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

the class CameraModelAnalysis method simulatePoissonGaussian.

/**
 * Randomly generate a histogram from poisson-gaussian samples.
 *
 * @param settings the settings
 * @param random the random
 * @return The histogram
 */
private static IntHistogram simulatePoissonGaussian(CameraModelAnalysisSettings settings, UniformRandomProvider random) {
    final int[] poissonSample = simulatePoisson(settings, random);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
    final double gain = getGain(settings);
    final double noise = getReadNoise(settings);
    final int samples = settings.getSamples();
    final int noiseSamples = (noise > 0) ? settings.getNoiseSamples() : 1;
    final int[] sample = new int[samples * noiseSamples];
    int count = 0;
    final Round round = getRound(settings);
    for (int n = poissonSample.length; n-- > 0; ) {
        // Fixed camera gain
        final double d = poissonSample[n] * gain;
        // Over-sample the Gaussian
        for (int j = noiseSamples; j-- > 0; ) {
            // Convert the sample to a count
            sample[count++] = round.round(d + noise * gauss.sample());
        }
    }
    return createHistogram(sample, count);
}
Also used : NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)

Example 7 with NormalizedGaussianSampler

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

the class DynamicMultipleTargetTracingTest method testTraceMoleculesDisableIntensityModel.

/**
 * Test trace molecules using 2 molecules. One is fixed and the other moves past it. The tracing
 * should assign the fixed molecule correctly as it has a low local diffusion rate.
 */
@Test
void testTraceMoleculesDisableIntensityModel() {
    final UniformRandomProvider rng = RngUtils.create(125631236L);
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    // localisation precision (in pixels)
    final double s = 0.1;
    final SharedStateContinuousSampler intensity = SamplerUtils.createGaussianSampler(rng, 1000, 100);
    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;
    final float x1 = size / 2 + 0.5f;
    for (int i = 0; i < size; i++) {
        results.add(new PeakResult(i, (float) (x1 + gauss.sample() * s), (float) (i + gauss.sample() * s), (float) (intensity.sample())));
    }
    // Second molecule is fixed in the centre with a same intensity
    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) (intensity.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, x1, size, (float) (intensity.sample())));
    // Move centre to centre each jump => 0.1 um or 0.01 um^2
    // MSD = 4D => D = 0.01 / 4 = 0.0025
    final DmttConfiguration config = DmttConfiguration.newBuilder(0.0025).setDisableIntensityModel(true).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));
    }
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) DmttConfiguration(uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration) CalibrationWriter(uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest) Test(org.junit.jupiter.api.Test)

Example 8 with NormalizedGaussianSampler

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

the class GaussianPsfModel method sample.

/**
 * Sample from a Gaussian distribution.
 *
 * @param n The number of samples
 * @param x0 The Gaussian centre in dimension 0
 * @param x1 The Gaussian centre in dimension 1
 * @param s0 The Gaussian standard deviation dimension 0
 * @param s1 The Gaussian standard deviation dimension 1
 * @param rng The random generator to use for sampling
 * @return The sample x and y values
 */
public double[][] sample(final int n, final double x0, final double x1, final double s0, final double s1, UniformRandomProvider rng) {
    this.s0 = s0;
    this.s1 = s1;
    final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(rng);
    final double[] x = sample(n, x0, s0, gauss);
    final double[] y = sample(n, x1, s1, gauss);
    return new double[][] { x, y };
}
Also used : NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)

Example 9 with NormalizedGaussianSampler

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

the class MoleculeModel method move.

/**
 * Move the molecule using a random Gaussian shift with standard deviation of the given diffusion
 * rate.
 *
 * <p>Note: The array provided by {@link #getCoordinates()} is updated and returned.
 *
 * @param diffusionRate Diffusion rate for each dimension
 * @param random Random generator
 * @return The new coordinates
 */
public double[] move(double diffusionRate, UniformRandomProvider random) {
    final double[] coords = getCoordinates();
    if (diffusionRate > 0) {
        final NormalizedGaussianSampler gauss = SamplerUtils.createNormalizedGaussianSampler(random);
        for (int i = 0; i < 3; i++) {
            final double shift = gauss.sample() * diffusionRate;
            // Clip the movement?
            // shift = MathUtils.clip(-5*diffusionRate,5*diffusionRate, shift)
            coords[i] += shift;
        }
    }
    return coords;
}
Also used : NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)

Example 10 with NormalizedGaussianSampler

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

the class DynamicMultipleTargetTracingTest method testTraceMoleculesDisableLocalDiffusionModel.

/**
 * Test trace molecules using 2 molecules. One is fixed and the other moves past it. The tracing
 * should assign the fixed molecule correctly as it has a low local diffusion rate.
 */
@Test
void testTraceMoleculesDisableLocalDiffusionModel() {
    // The test is not very robust and fails 20% 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())));
    // 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).setDisableLocalDiffusionModel(true).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));
    }
}
Also used : SharedStateContinuousSampler(org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler) DmttConfiguration(uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration) CalibrationWriter(uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) NormalizedGaussianSampler(org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler) SeededTest(uk.ac.sussex.gdsc.test.junit5.SeededTest) Test(org.junit.jupiter.api.Test)

Aggregations

NormalizedGaussianSampler (org.apache.commons.rng.sampling.distribution.NormalizedGaussianSampler)16 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)8 ArrayList (java.util.ArrayList)3 SharedStateContinuousSampler (org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler)3 Test (org.junit.jupiter.api.Test)3 StoredDataStatistics (uk.ac.sussex.gdsc.core.utils.StoredDataStatistics)3 MarsagliaTsangGammaSampler (uk.ac.sussex.gdsc.core.utils.rng.MarsagliaTsangGammaSampler)3 CalibrationWriter (uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter)3 DmttConfiguration (uk.ac.sussex.gdsc.smlm.results.DynamicMultipleTargetTracing.DmttConfiguration)3 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)3 Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)2 StoredData (uk.ac.sussex.gdsc.core.utils.StoredData)2 SeededTest (uk.ac.sussex.gdsc.test.junit5.SeededTest)2 TDoubleArrayList (gnu.trove.list.array.TDoubleArrayList)1 TIntHashSet (gnu.trove.set.hash.TIntHashSet)1 IJ (ij.IJ)1 ImagePlus (ij.ImagePlus)1 Plot (ij.gui.Plot)1 Calibration (ij.measure.Calibration)1 PlugIn (ij.plugin.PlugIn)1