Search in sources :

Example 6 with ErfGaussian2DFunction

use of 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 npeaks
	 *            the npeaks
	 * @param params
	 *            set on output
	 * @param randomiseParams
	 *            Set to true to randomise the params
	 * @return the double[]
	 */
private double[] doubleCreateGaussianData(int npeaks, double[] params, boolean randomiseParams) {
    int n = blockWidth * blockWidth;
    // Generate a 2D Gaussian
    ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    params[0] = random(Background);
    for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
        params[j] = random(Signal);
        params[j + 2] = random(Xpos);
        params[j + 3] = random(Ypos);
        params[j + 4] = random(Xwidth);
        params[j + 5] = random(Ywidth);
    }
    double[] y = new double[n];
    func.initialise(params);
    for (int i = 0; i < y.length; i++) {
        // Add random Poisson noise
        y[i] = rdg.nextPoisson(func.eval(i));
    }
    if (randomiseParams) {
        params[0] = random(params[0]);
        for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
            params[j] = random(params[j]);
            params[j + 2] = random(params[j + 2]);
            params[j + 3] = random(params[j + 3]);
            params[j + 4] = random(params[j + 4]);
            params[j + 5] = random(params[j + 5]);
        }
    }
    return y;
}
Also used : ErfGaussian2DFunction(gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) SingleFreeCircularErfGaussian2DFunction(gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction)

Example 7 with ErfGaussian2DFunction

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

the class FastMLEGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.

@Test
public void gradientProcedureComputesSameWithPrecomputed() {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    double[] a1 = new double[7];
    double[] a2 = new double[13];
    final double[] x = new double[f1.size()];
    final double[] b = new double[f1.size()];
    for (int i = 0; i < iter; i++) {
        a2[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(0.1, 0.3);
        a2[Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a2[Gaussian2DFunction.X_POSITION] = rdg.nextUniform(3, 5);
        a2[Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(3, 5);
        a2[Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
        a2[Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
        a2[6 + Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a2[6 + Gaussian2DFunction.X_POSITION] = rdg.nextUniform(5, 7);
        a2[6 + Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(5, 7);
        a2[6 + Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
        a2[6 + Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
        // Simulate Poisson data
        f2.initialise0(a2);
        f1.forEach(new ValueProcedure() {

            int k = 0;

            public void execute(double value) {
                x[k++] = (value > 0) ? rdg.nextPoisson(value) : 0;
            }
        });
        // Precompute peak 2 (no background)
        a1[Gaussian2DFunction.BACKGROUND] = 0;
        for (int j = 1; j < 7; j++) a1[j] = a2[6 + j];
        f1.initialise0(a1);
        f1.forEach(new ValueProcedure() {

            int k = 0;

            public void execute(double value) {
                b[k++] = value;
            }
        });
        // Reset to peak 1
        for (int j = 0; j < 7; j++) a1[j] = a2[j];
        // Compute peak 1+2
        FastMLEGradient2Procedure p12 = FastMLEGradient2ProcedureFactory.create(x, f2);
        p12.computeSecondDerivative(a2);
        double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
        double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
        // Compute peak 1+(precomputed 2)
        FastMLEGradient2Procedure p1b2 = FastMLEGradient2ProcedureFactory.create(x, PrecomputedGradient2Function.wrapGradient2Function(f1, b));
        p1b2.computeSecondDerivative(a1);
        double[] d12 = p1b2.d1;
        double[] d22 = p1b2.d2;
        Assert.assertArrayEquals(" Result: Not same @ " + i, p12.u, p1b2.u, 1e-10);
        Assert.assertArrayEquals(" D1: Not same @ " + i, d11, d12, 1e-10);
        Assert.assertArrayEquals(" D2: Not same @ " + i, d21, d22, 1e-10);
        double[] v1 = p12.computeValue(a2);
        double[] v2 = p1b2.computeValue(a1);
        Assert.assertArrayEquals(" Value: Not same @ " + i, v1, v2, 1e-10);
        double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
        double[] d2 = p1b2.computeFirstDerivative(a1);
        Assert.assertArrayEquals(" 1st derivative: Not same @ " + i, d1, d2, 1e-10);
    }
}
Also used : ValueProcedure(gdsc.smlm.function.ValueProcedure) ErfGaussian2DFunction(gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction) SingleFreeCircularErfGaussian2DFunction(gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction) SingleAstigmatismErfGaussian2DFunction(gdsc.smlm.function.gaussian.erf.SingleAstigmatismErfGaussian2DFunction) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) Well19937c(org.apache.commons.math3.random.Well19937c) Test(org.junit.Test)

Aggregations

ErfGaussian2DFunction (gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction)7 SingleFreeCircularErfGaussian2DFunction (gdsc.smlm.function.gaussian.erf.SingleFreeCircularErfGaussian2DFunction)5 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)3 Well19937c (org.apache.commons.math3.random.Well19937c)3 Test (org.junit.Test)3 SingleAstigmatismErfGaussian2DFunction (gdsc.smlm.function.gaussian.erf.SingleAstigmatismErfGaussian2DFunction)2 NotImplementedException (gdsc.core.utils.NotImplementedException)1 Statistics (gdsc.core.utils.Statistics)1 FisherInformationMatrix (gdsc.smlm.fitting.FisherInformationMatrix)1 EJMLLinearSolver (gdsc.smlm.fitting.linear.EJMLLinearSolver)1 ValueProcedure (gdsc.smlm.function.ValueProcedure)1 ArrayList (java.util.ArrayList)1