Examples with ErfGaussian2DFunction uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction used on opensource projects

Example 11 with ErfGaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class LsqLvmGradientProcedureTest method doubleCreateGaussianData.

/**
 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude,
 * angle, xpos, ypos, xwidth, ywidth }
 *
 * @param rng the random
 * @param npeaks the npeaks
 * @param params set on output
 * @param randomiseParams Set to true to randomise the params
 * @return the double[]
 */
private double[] doubleCreateGaussianData(UniformRandomProvider rng, int npeaks, double[] params, boolean randomiseParams) {
    final int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    final ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    params[0] = random(rng, background);
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = random(rng, signal);
        params[j + 2] = random(rng, xpos);
        params[j + 3] = random(rng, ypos);
        params[j + 4] = random(rng, xwidth);
        params[j + 5] = random(rng, ywidth);
    }
    final double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        final double u = func.eval(i);
        y[i] = GdscSmlmTestUtils.createPoissonSampler(rng, u).sample();
    }
    if (randomiseParams) {
        params[0] = random(rng, params[0]);
        for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
            params[j] = random(rng, params[j]);
            params[j + 2] = random(rng, params[j + 2]);
            params[j + 3] = random(rng, params[j + 3]);
            params[j + 4] = random(rng, params[j + 4]);
            params[j + 5] = random(rng, params[j + 5]);
        }
    }
    return y;
}

Example 12 with ErfGaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class LvmGradientProcedureTest method doubleCreateGaussianData.

/**
 * Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
 * Only the chosen parameters are randomised and returned for a maximum of (background, amplitude,
 * angle, xpos, ypos, xwidth, ywidth }
 *
 * @param npeaks the npeaks
 * @param params set on output
 * @param randomiseParams Set to true to randomise the params
 * @return the double[]
 */
private double[] doubleCreateGaussianData(UniformRandomProvider rng, int npeaks, double[] params, boolean randomiseParams) {
    final int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    final ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    params[0] = random(rng, background);
    for (int i = 0, j = 0; i < npeaks; i++, j += Gaussian2DFunction.PARAMETERS_PER_PEAK) {
        params[j + Gaussian2DFunction.SIGNAL] = random(rng, signal);
        params[j + Gaussian2DFunction.X_POSITION] = random(rng, xpos);
        params[j + Gaussian2DFunction.Y_POSITION] = random(rng, ypos);
        params[j + Gaussian2DFunction.X_SD] = random(rng, xwidth);
        params[j + Gaussian2DFunction.Y_SD] = random(rng, ywidth);
    }
    if (npeaks > 1) {
        // Move the peaks around so they do not overlap
        final double[] shift = SimpleArrayUtils.newArray(npeaks, -2, 4.0 / (npeaks - 1));
        RandomUtils.shuffle(shift, rng);
        for (int i = 0, j = 0; i < npeaks; i++, j += Gaussian2DFunction.PARAMETERS_PER_PEAK) {
            params[j + Gaussian2DFunction.X_POSITION] += shift[i];
        }
        RandomUtils.shuffle(shift, rng);
        for (int i = 0, j = 0; i < npeaks; i++, j += Gaussian2DFunction.PARAMETERS_PER_PEAK) {
            params[j + Gaussian2DFunction.Y_POSITION] += shift[i];
        }
    }
    final double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        final double u = func.eval(i);
        y[i] = GdscSmlmTestUtils.createPoissonSampler(rng, u).sample();
    }
    if (randomiseParams) {
        params[0] = random(rng, params[0]);
        for (int i = 0, j = 0; i < npeaks; i++, j += Gaussian2DFunction.PARAMETERS_PER_PEAK) {
            params[j + Gaussian2DFunction.SIGNAL] = random(rng, params[j + Gaussian2DFunction.SIGNAL]);
            params[j + Gaussian2DFunction.X_POSITION] = random(rng, params[j + Gaussian2DFunction.X_POSITION]);
            params[j + Gaussian2DFunction.Y_POSITION] = random(rng, params[j + Gaussian2DFunction.Y_POSITION]);
            params[j + Gaussian2DFunction.X_SD] = random(rng, params[j + Gaussian2DFunction.X_SD]);
            params[j + Gaussian2DFunction.Y_SD] = random(rng, params[j + Gaussian2DFunction.Y_SD]);
        }
    }
    return y;
}

Example 13 with ErfGaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class LsqVarianceGradientProcedureTest method crlbIsHigherWithPrecomputed.

@SeededTest
void crlbIsHigherWithPrecomputed(RandomSeed seed) {
    final int iter = 10;
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    final ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final double[] a = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    final int n = func.getNumberOfGradients();
    // Get a background
    final double[] b = new double[func.size()];
    for (int i = 0; i < b.length; i++) {
        b[i] = nextUniform(rng, 1, 2);
    }
    for (int i = 0; i < iter; i++) {
        a[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
        a[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
        a[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
        a[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
        a[Gaussian2DFunction.Y_SD] = nextUniform(rng, 1, 1.3);
        final LsqVarianceGradientProcedure p1 = LsqVarianceGradientProcedureUtils.create(func);
        p1.variance(a);
        final LsqVarianceGradientProcedure p2 = LsqVarianceGradientProcedureUtils.create(OffsetGradient1Function.wrapGradient1Function(func, b));
        p2.variance(a);
        final double[] crlb1 = p1.variance;
        final double[] crlb2 = p2.variance;
        Assertions.assertNotNull(crlb1);
        Assertions.assertNotNull(crlb2);
        // Arrays.toString(crlb2));
        for (int j = 0; j < n; j++) {
            Assertions.assertTrue(crlb1[j] < crlb2[j]);
        }
    }
}

Example 14 with ErfGaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class BaseSteppingFunctionSolverTest method getSolver.

SteppingFunctionSolver getSolver(SteppingFunctionSolverClamp clamp, SteppingFunctionSolverType type, ToleranceChecker tc) {
    final ErfGaussian2DFunction f = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, size, size, flags, null);
    final ParameterBounds bounds = new ParameterBounds(f);
    switch(clamp) {
        case DYNAMIC_CLAMP:
            bounds.setDynamicClamp(true);
            bounds.setClampValues(defaultClampValues);
            break;
        case CLAMP:
            bounds.setClampValues(defaultClampValues);
            break;
        case NO_CLAMP:
        default:
            break;
    }
    SteppingFunctionSolver solver;
    switch(type) {
        case LSELVM:
            solver = new LseLvmSteppingFunctionSolver(f, tc, bounds);
            break;
        case MLELVM:
        case FastLogMLELVM:
            final MleLvmSteppingFunctionSolver mleSolver = new MleLvmSteppingFunctionSolver(f, tc, bounds);
            solver = mleSolver;
            // MLE requires a positive function value so use a lower bound
            solver.setBounds(getLb(), null);
            // For testing the fast log version
            if (type == FastLogMLELVM) {
                mleSolver.setFastLog(FastLogFactory.getFastLog());
            }
            break;
        case WLSELVM:
            solver = new WLseLvmSteppingFunctionSolver(f, tc, bounds);
            break;
        case FastMLE:
            solver = new FastMleSteppingFunctionSolver(f, tc, bounds);
            // MLE requires a positive function value so use a lower bound
            solver.setBounds(getLb(), null);
            break;
        default:
            throw new NotImplementedException();
    }
    if (solver instanceof LvmSteppingFunctionSolver) {
        ((LvmSteppingFunctionSolver) solver).setInitialLambda(1);
    }
    return solver;
}

Example 15 with ErfGaussian2DFunction

use of uk.ac.sussex.gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.

the class PsfModelGradient1FunctionTest method canComputeValueAndGradient.

@Test
void canComputeValueAndGradient() {
    // Use a reasonable z-depth function from the Smith, et al (2010) paper (page 377)
    final double sx = 1.08;
    final double sy = 1.01;
    final double gamma = 0.389;
    final double d = 0.531;
    final double Ax = -0.0708;
    final double Bx = -0.073;
    final double Ay = 0.164;
    final double By = 0.0417;
    final AstigmatismZModel zModel = HoltzerAstigmatismZModel.create(sx, sy, gamma, d, Ax, Bx, Ay, By);
    // Small size ensure the PSF model covers the entire image
    final int maxx = 11;
    final int maxy = 11;
    final double[] ve = new double[maxx * maxy];
    final double[] vo = new double[maxx * maxy];
    final double[][] ge = new double[maxx * maxy][];
    final double[][] go = new double[maxx * maxy][];
    final PsfModelGradient1Function psf = new PsfModelGradient1Function(new GaussianPsfModel(zModel), maxx, maxy);
    final ErfGaussian2DFunction f = new SingleAstigmatismErfGaussian2DFunction(maxx, maxy, zModel);
    f.setErfFunction(ErfFunction.COMMONS_MATH);
    final double[] a2 = new double[Gaussian2DFunction.PARAMETERS_PER_PEAK + 1];
    final DoubleDoubleBiPredicate equality = TestHelper.doublesAreClose(1e-8, 0);
    final double c = maxx * 0.5;
    for (int i = -1; i <= 1; i++) {
        final double x0 = c + i * 0.33;
        for (int j = -1; j <= 1; j++) {
            final double x1 = c + j * 0.33;
            for (int k = -1; k <= 1; k++) {
                final double x2 = k * 0.33;
                for (final double in : new double[] { 23.2, 405.67 }) {
                    // Background is constant for gradients so just use 1 value
                    final double[] a = new double[] { 2.2, in, x0, x1, x2 };
                    psf.initialise1(a);
                    psf.forEach(new Gradient1Procedure() {

                        int index = 0;

                        @Override
                        public void execute(double value, double[] dyDa) {
                            vo[index] = value;
                            go[index] = dyDa.clone();
                            index++;
                        }
                    });
                    a2[Gaussian2DFunction.BACKGROUND] = a[0];
                    a2[Gaussian2DFunction.SIGNAL] = a[1];
                    a2[Gaussian2DFunction.X_POSITION] = a[2] - 0.5;
                    a2[Gaussian2DFunction.Y_POSITION] = a[3] - 0.5;
                    a2[Gaussian2DFunction.Z_POSITION] = a[4];
                    f.initialise1(a2);
                    f.forEach(new Gradient1Procedure() {

                        int index = 0;

                        @Override
                        public void execute(double value, double[] dyDa) {
                            ve[index] = value;
                            ge[index] = dyDa.clone();
                            index++;
                        }
                    });
                    for (int ii = 0; ii < ve.length; ii++) {
                        TestAssertions.assertTest(ve[ii], vo[ii], equality);
                        TestAssertions.assertArrayTest(ge[ii], go[ii], equality);
                    }
                }
            }
        }
    }
}