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

Example 21 with SharedStateContinuousSampler

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

the class WeightedKernelFilterTest method filterPerformsWeightedKernelFiltering.

@SeededTest
void filterPerformsWeightedKernelFiltering(RandomSeed seed) {
    final DataFilter filter = createDataFilter();
    final UniformRandomProvider rg = RngFactory.create(seed.get());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, 2, 0.2);
    final float[] offsets = getOffsets(filter);
    final int[] boxSizes = getBoxSizes(filter);
    final FloatFloatBiPredicate equality = Predicates.floatsAreClose(1e-6, 0);
    for (final int width : primes) {
        for (final int height : new int[] { 29 }) {
            final float[] data = createData(width, height, rg);
            // Uniform weights
            final float[] w1 = new float[width * height];
            Arrays.fill(w1, 0.5f);
            // Weights simulating the variance of sCMOS pixels
            final float[] w2 = new float[width * height];
            for (int i = 0; i < w2.length; i++) {
                w2[i] = (float) (1.0 / Math.max(0.01, gs.sample()));
            }
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        // For each pixel over the range around the pixel (vi).
                        // kernel filter: sum(vi * ki) / sum(ki)
                        // Weighted kernel filter: sum(vi * wi * ki) / sum(ki * wi)
                        // Note: The kernel filter is like a weighted filter
                        // (New kernel = wi * ki)
                        filter.setWeights(null, width, height);
                        // Uniform weights
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w1, boxSize, offset, internal, equality);
                        // Random weights.
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w2, boxSize, offset, internal, equality);
                    }
                }
            }
        }
    }
}

Example 22 with SharedStateContinuousSampler

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

the class FrcTest method canComputeMirrored.

@SeededTest
void canComputeMirrored(RandomSeed seed) {
    // Sample lines through an image to create a structure.
    final int size = 1024;
    final double[][] data = new double[size * 2][];
    final UniformRandomProvider r = RngFactory.create(seed.get());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
    for (int x = 0, y = 0, y2 = size, i = 0; x < size; x++, y++, y2--) {
        data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
        data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
    }
    // Create 2 images
    final Rectangle bounds = new Rectangle(0, 0, size, size);
    ImageJImagePeakResults i1 = createImage(bounds);
    ImageJImagePeakResults i2 = createImage(bounds);
    final int[] indices = SimpleArrayUtils.natural(data.length);
    PermutationSampler.shuffle(r, indices);
    for (final int i : indices) {
        final ImageJImagePeakResults image = i1;
        i1 = i2;
        i2 = image;
        image.add((float) data[i][0], (float) data[i][1], 1);
    }
    i1.end();
    i2.end();
    final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
    final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
    // Test
    final Frc frc = new Frc();
    FloatProcessor[] fft1;
    FloatProcessor[] fft2;
    fft1 = frc.getComplexFft(ip1);
    fft2 = frc.getComplexFft(ip2);
    final float[] dataA1 = (float[]) fft1[0].getPixels();
    final float[] dataB1 = (float[]) fft1[1].getPixels();
    final float[] dataA2 = (float[]) fft2[0].getPixels();
    final float[] dataB2 = (float[]) fft2[1].getPixels();
    final float[] numeratorE = new float[dataA1.length];
    final float[] absFft1E = new float[dataA1.length];
    final float[] absFft2E = new float[dataA1.length];
    Frc.compute(numeratorE, absFft1E, absFft2E, dataA1, dataB1, dataA2, dataB2);
    Assertions.assertTrue(checkSymmetry(numeratorE, size), "numeratorE");
    Assertions.assertTrue(checkSymmetry(absFft1E, size), "absFft1E");
    Assertions.assertTrue(checkSymmetry(absFft2E, size), "absFft2E");
    final float[] numeratorA = new float[dataA1.length];
    final float[] absFft1A = new float[dataA1.length];
    final float[] absFft2A = new float[dataA1.length];
    Frc.computeMirrored(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
    // for (int y=0, i=0; y<size; y++)
    // for (int x=0; x<size; x++, i++)
    // {
    // logger.fine(FormatSupplier.getSupplier("[%d,%d = %d] %f ?= %f", x, y, i, numeratorE[i],
    // numeratorA[i]);
    // }
    Assertions.assertArrayEquals(numeratorE, numeratorA, "numerator");
    Assertions.assertArrayEquals(absFft1E, absFft1A, "absFft1");
    Assertions.assertArrayEquals(absFft2E, absFft2A, "absFft2");
    Frc.computeMirroredFast(size, numeratorA, absFft1A, absFft2A, dataA1, dataB1, dataA2, dataB2);
    // Check this.
    for (int y = 1; y < size; y++) {
        for (int x = 1, i = y * size + 1; x < size; x++, i++) {
            Assertions.assertEquals(numeratorE[i], numeratorA[i], "numerator");
            Assertions.assertEquals(absFft1E[i], absFft1A[i], "absFft1");
            Assertions.assertEquals(absFft2E[i], absFft2A[i], "absFft2");
        }
    }
}

Example 23 with SharedStateContinuousSampler

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

the class FrcTest method computeMirroredIsFaster.

@SeededTest
void computeMirroredIsFaster(RandomSeed seed) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    // Sample lines through an image to create a structure.
    final int N = 2048;
    final double[][] data = new double[N * 2][];
    final UniformRandomProvider r = RngUtils.create(seed.getSeed());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(r, 0, 5);
    for (int x = 0, y = 0, y2 = N, i = 0; x < N; x++, y++, y2--) {
        data[i++] = new double[] { x + gs.sample(), y + gs.sample() };
        data[i++] = new double[] { x + gs.sample(), y2 + gs.sample() };
    }
    // Create 2 images
    final Rectangle bounds = new Rectangle(0, 0, N, N);
    ImageJImagePeakResults i1 = createImage(bounds);
    ImageJImagePeakResults i2 = createImage(bounds);
    final int[] indices = SimpleArrayUtils.natural(data.length);
    PermutationSampler.shuffle(r, indices);
    for (final int i : indices) {
        final ImageJImagePeakResults image = i1;
        i1 = i2;
        i2 = image;
        image.add((float) data[i][0], (float) data[i][1], 1);
    }
    i1.end();
    i2.end();
    final ImageProcessor ip1 = i1.getImagePlus().getProcessor();
    final ImageProcessor ip2 = i2.getImagePlus().getProcessor();
    // Test
    final Frc frc = new Frc();
    FloatProcessor[] fft1;
    FloatProcessor[] fft2;
    fft1 = frc.getComplexFft(ip1);
    fft2 = frc.getComplexFft(ip2);
    final float[] dataA1 = (float[]) fft1[0].getPixels();
    final float[] dataB1 = (float[]) fft1[1].getPixels();
    final float[] dataA2 = (float[]) fft2[0].getPixels();
    final float[] dataB2 = (float[]) fft2[1].getPixels();
    final float[] numerator = new float[dataA1.length];
    final float[] absFft1 = new float[dataA1.length];
    final float[] absFft2 = new float[dataA1.length];
    final TimingService ts = new TimingService(10);
    ts.execute(new MyTimingTask("compute") {

        @Override
        public Object run(Object data) {
            Frc.compute(numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    ts.execute(new MyTimingTask("computeMirrored") {

        @Override
        public Object run(Object data) {
            Frc.computeMirrored(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    ts.execute(new MyTimingTask("computeMirroredFast") {

        @Override
        public Object run(Object data) {
            Frc.computeMirroredFast(N, numerator, absFft1, absFft2, dataA1, dataB1, dataA2, dataB2);
            return null;
        }
    });
    final int size = ts.getSize();
    ts.repeat(size);
    if (logger.isLoggable(Level.INFO)) {
        logger.info(ts.getReport(size));
    }
    // The 'Fast' method may not always be faster so just log results
    final TimingResult slow = ts.get(-3);
    final TimingResult fast = ts.get(-2);
    final TimingResult fastest = ts.get(-1);
    logger.log(TestLogUtils.getTimingRecord(slow, fastest));
    logger.log(TestLogUtils.getTimingRecord(fast, fastest));
    // It should be faster than the non mirrored version
    Assertions.assertTrue(ts.get(-1).getMean() <= ts.get(-3).getMean());
}

Example 24 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project gdsc by aherbert.

the class FindFociTest method combine.

private static short[] combine(UniformRandomProvider rg, float[] data1, float[] data2, float[] data3) {
    // Combine images and add a bias and read noise
    final SharedStateContinuousSampler g = SamplerUtils.createGaussianSampler(rg, BIAS, 5);
    final short[] data = new short[data1.length];
    for (int i = 0; i < data.length; i++) {
        final double mu = data1[i] + data2[i] + data3[i];
        double value = g.sample();
        if (mu != 0) {
            value += new PoissonSampler(rg, mu).sample();
        }
        if (value < 0) {
            value = 0;
        } else if (value > 65535) {
            value = 65535;
        }
        data[i] = (short) value;
    }
    return data;
}

Example 25 with SharedStateContinuousSampler

use of org.apache.commons.rng.sampling.distribution.SharedStateContinuousSampler in project gdsc-smlm by aherbert.

the class WeightedKernelFilterTest method filterPerformsWeightedKernelFiltering.

@SeededTest
void filterPerformsWeightedKernelFiltering(RandomSeed seed) {
    final DataFilter filter = createDataFilter();
    final UniformRandomProvider rg = RngFactory.create(seed.get());
    final SharedStateContinuousSampler gs = SamplerUtils.createGaussianSampler(rg, 2, 0.2);
    final float[] offsets = getOffsets(filter);
    final int[] boxSizes = getBoxSizes(filter);
    final FloatFloatBiPredicate equality = Predicates.floatsAreClose(1e-6, 0);
    for (final int width : primes) {
        for (final int height : new int[] { 29 }) {
            final float[] data = createData(width, height, rg);
            // Uniform weights
            final float[] w1 = new float[width * height];
            Arrays.fill(w1, 0.5f);
            // Weights simulating the variance of sCMOS pixels
            final float[] w2 = new float[width * height];
            for (int i = 0; i < w2.length; i++) {
                w2[i] = (float) (1.0 / Math.max(0.01, gs.sample()));
            }
            for (final int boxSize : boxSizes) {
                for (final float offset : offsets) {
                    for (final boolean internal : checkInternal) {
                        // For each pixel over the range around the pixel (vi).
                        // kernel filter: sum(vi * ki) / sum(ki)
                        // Weighted kernel filter: sum(vi * wi * ki) / sum(ki * wi)
                        // Note: The kernel filter is like a weighted filter
                        // (New kernel = wi * ki)
                        filter.setWeights(null, width, height);
                        // Uniform weights
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w1, boxSize, offset, internal, equality);
                        // Random weights.
                        testfilterPerformsWeightedFiltering(filter, width, height, data, w2, boxSize, offset, internal, equality);
                    }
                }
            }
        }
    }
}