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

Example 6 with ErfGaussian2DFunction

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

the class FastMleGradient2ProcedureTest 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 = 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);
    }
    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 7 with ErfGaussian2DFunction

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

the class CubicSplineFunctionTest method speedTest.

@SuppressWarnings("null")
private void speedTest(int n, int order) {
    // No assertions, this is just a report
    Assumptions.assumeTrue(logger.isLoggable(Level.INFO));
    // Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final CubicSplineFunction cf = (n == 2) ? f2 : f1;
    Assumptions.assumeTrue(null != cf);
    final CubicSplineFunction cff = (n == 2) ? f2f : f1f;
    final ErfGaussian2DFunction gf = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(n, maxx, maxy, GaussianFunctionFactory.FIT_ASTIGMATISM, zModel);
    final Gaussian2DFunction gf2 = (order < 2) ? GaussianFunctionFactory.create2D(n, maxx, maxy, GaussianFunctionFactory.FIT_SIMPLE_FREE_CIRCLE, zModel) : null;
    final LocalList<double[]> l1 = new LocalList<>();
    final LocalList<double[]> l2 = new LocalList<>();
    final LocalList<double[]> l3 = new LocalList<>();
    final double[] a = new double[1 + n * CubicSplineFunction.PARAMETERS_PER_PEAK];
    final double[] b = new double[1 + n * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    double[] bb = null;
    a[CubicSplineFunction.BACKGROUND] = 0.1;
    b[Gaussian2DFunction.BACKGROUND] = 0.1;
    for (int i = 0; i < n; i++) {
        a[i * CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.SIGNAL] = 10;
        b[i * Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.SIGNAL] = 10;
    }
    if (n == 2) {
        // Fix second peak parameters
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.X_POSITION] = testcx1[0];
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.Y_POSITION] = testcy1[0];
        a[CubicSplineFunction.PARAMETERS_PER_PEAK + CubicSplineFunction.Z_POSITION] = testcz1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] = testcx1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] = testcy1[0];
        b[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Z_POSITION] = testcz1[0];
    }
    if (gf2 != null) {
        bb = b.clone();
        if (n == 2) {
            // Fix second peak parameters
            bb[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_SD] = zModel.getSx(testcz1[0]);
            bb[Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_SD] = zModel.getSy(testcz1[0]);
        }
    }
    for (int x = 0; x <= maxx; x++) {
        a[CubicSplineFunction.X_POSITION] = x;
        b[Gaussian2DFunction.X_POSITION] = x;
        for (int y = 0; y <= maxy; y++) {
            a[CubicSplineFunction.Y_POSITION] = y;
            b[Gaussian2DFunction.Y_POSITION] = y;
            for (int z = -zDepth; z <= zDepth; z++) {
                a[CubicSplineFunction.Z_POSITION] = z;
                b[Gaussian2DFunction.Z_POSITION] = z;
                l1.add(a.clone());
                l2.add(b.clone());
                if (gf2 != null) {
                    bb[Gaussian2DFunction.X_SD] = zModel.getSx(z);
                    bb[Gaussian2DFunction.Y_SD] = zModel.getSy(z);
                    l3.add(bb.clone());
                }
            }
        }
    }
    final double[][] x1 = l1.toArray(new double[0][]);
    final double[][] x2 = l2.toArray(new double[0][]);
    final double[][] x3 = l3.toArray(new double[0][]);
    final TimingService ts = new TimingService(5);
    ts.execute(new FunctionTimingTask(gf, x2, order));
    if (gf2 != null) {
        ts.execute(new FunctionTimingTask(gf2, x3, order));
    }
    ts.execute(new FunctionTimingTask(cf, x1, order));
    ts.execute(new FunctionTimingTask(cff, x1, order, " single-precision"));
    final int size = ts.getSize();
    ts.repeat(size);
    logger.info(ts.getReport(size));
}

Example 8 with ErfGaussian2DFunction

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

the class GaussianPsfModel method createGaussianFunction.

private ErfGaussian2DFunction createGaussianFunction(final int x0range, final int x1range) {
    final ErfGaussian2DFunction f = new SingleAstigmatismErfGaussian2DFunction(x0range, x1range, zModel);
    f.setErfFunction(ErfGaussian2DFunction.ErfFunction.COMMONS_MATH);
    return f;
}

Example 9 with ErfGaussian2DFunction

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

the class GaussianPsfModel method computeValue.

@Override
protected boolean computeValue(int width, int height, double x0, double x1, double x2, double[] value) {
    final double s0 = getS0(x2);
    final double s1 = getS1(x2);
    // Evaluate the Gaussian error function on a pixel grid covering +/- 5 SD
    final int x0min = clip((int) (x0 - range * s0), width);
    final int x1min = clip((int) (x1 - range * s1), height);
    final int x0max = clip((int) Math.ceil(x0 + range * s0), width);
    final int x1max = clip((int) Math.ceil(x1 + range * s1), height);
    final int x0range = x0max - x0min;
    final int x1range = x1max - x1min;
    // min should always be less than max
    ValidationUtils.checkStrictlyPositive(x0range, "Range0");
    ValidationUtils.checkStrictlyPositive(x1range, "Range1");
    // Evaluate using a function.
    // This allows testing the function verses the default gaussian2D() method
    // as they should be nearly identical (the Erf function may be different).
    // Thus the gradient can be evaluated using the same function.
    final ErfGaussian2DFunction f = createGaussianFunction(x0range, x1range);
    final double[] p = new double[Gaussian2DFunction.PARAMETERS_PER_PEAK + 1];
    p[Gaussian2DFunction.SIGNAL] = 1;
    // The function computes the centre of the pixel as 0,0
    p[Gaussian2DFunction.X_POSITION] = x0 - x0min - 0.5;
    p[Gaussian2DFunction.Y_POSITION] = x1 - x1min - 0.5;
    p[Gaussian2DFunction.Z_POSITION] = x2;
    final double[] v = f.computeValues(p);
    for (int y = 0; y < x1range; y++) {
        // Locate the insert location
        int indexTo = (y + x1min) * width + x0min;
        int indexFrom = y * x0range;
        for (int x = 0; x < x0range; x++) {
            value[indexTo++] = v[indexFrom++];
        }
    }
    return true;
}

Example 10 with ErfGaussian2DFunction

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

the class GaussianPsfModel method computeValueAndGradient.

@Override
protected boolean computeValueAndGradient(int width, int height, double x0, double x1, double x2, double[] value, double[][] jacobian) {
    final double s0 = getS0(x2);
    final double s1 = getS1(x2);
    // Evaluate the Gaussian error function on a pixel grid covering +/- 5 SD
    final int x0min = clip((int) (x0 - range * s0), width);
    final int x1min = clip((int) (x1 - range * s1), height);
    final int x0max = clip((int) Math.ceil(x0 + range * s0), width);
    final int x1max = clip((int) Math.ceil(x1 + range * s1), height);
    final int x0range = x0max - x0min;
    final int x1range = x1max - x1min;
    // min should always be less than max
    ValidationUtils.checkStrictlyPositive(x0range, "Range0");
    ValidationUtils.checkStrictlyPositive(x1range, "Range1");
    // Evaluate using a function.
    // This allows testing the function verses the default gaussian2D() method
    // as they should be nearly identical (the Erf function may be different).
    // Thus the gradient can be evaluated using the same function.
    final ErfGaussian2DFunction f = createGaussianFunction(x0range, x1range);
    final double[] p = new double[Gaussian2DFunction.PARAMETERS_PER_PEAK + 1];
    p[Gaussian2DFunction.SIGNAL] = 1;
    // The function computes the centre of the pixel as 0,0.
    // The PSF sets the centre as 0.5,0.5.
    p[Gaussian2DFunction.X_POSITION] = x0 - x0min - 0.5;
    p[Gaussian2DFunction.Y_POSITION] = x1 - x1min - 0.5;
    p[Gaussian2DFunction.Z_POSITION] = x2;
    final int i0 = f.findGradientIndex(Gaussian2DFunction.X_POSITION);
    final int i1 = f.findGradientIndex(Gaussian2DFunction.Y_POSITION);
    final int i2 = f.findGradientIndex(Gaussian2DFunction.Z_POSITION);
    final Pair<double[], double[][]> pair = f.computeValuesAndJacobian(p);
    final double[] v = pair.getKey();
    final double[][] j = pair.getValue();
    for (int y = 0; y < x1range; y++) {
        // Locate the insert location
        int indexTo = (y + x1min) * width + x0min;
        int indexFrom = y * x0range;
        for (int x = 0; x < x0range; x++) {
            value[indexTo] = v[indexFrom];
            jacobian[indexTo][0] = j[indexFrom][i0];
            jacobian[indexTo][1] = j[indexFrom][i1];
            jacobian[indexTo][2] = j[indexFrom][i2];
            indexTo++;
            indexFrom++;
        }
    }
    return true;
}