Examples with Gaussian2DFunction gdsc.smlm.function.gaussian.Gaussian2DFunction used on opensource projects

Example 1 with Gaussian2DFunction

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

the class Gaussian2DFunctionTest method functionComputesGaussian.

@Test
public void functionComputesGaussian() {
    double background = 0;
    int maxx = 30;
    Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, zModel);
    Gaussian2DFunction f2;
    if ((flags & GaussianFunctionFactory.FIT_ERF) == 0)
        f2 = GaussianFunctionFactory.create2D(1, maxx, maxx, GaussianFunctionFactory.FIT_ELLIPTICAL, zModel);
    else
        f2 = GaussianFunctionFactory.create2D(1, maxx, maxx, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, zModel);
    boolean zDepth = (flags & GaussianFunctionFactory.FIT_Z) != 0;
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : new double[] { maxx / 2 + 0.373f }) for (double cy1 : new double[] { maxx / 2 + 0.876f }) for (double[] w1 : testw1) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        f.initialise(a);
        if (zDepth) {
            // Change to a standard free circular function
            a[Gaussian2DFunction.X_SD] *= zModel.getSx(a[Gaussian2DFunction.SHAPE]);
            a[Gaussian2DFunction.Y_SD] *= zModel.getSy(a[Gaussian2DFunction.SHAPE]);
            a[Gaussian2DFunction.SHAPE] = 0;
        }
        f2.initialise(a);
        double sum = 0;
        for (int index = maxx * maxx; index-- > 0; ) {
            double r1 = f.eval(index);
            double r2 = f2.eval(index);
            //System.out.printf("%d,%d r1=%f\n", index%maxx, index/maxx, r1);
            sum += r1;
            final boolean ok = eq2.almostEqualRelativeOrAbsolute(r1, r2);
            if (!ok)
                Assert.assertTrue(String.format("%g != %g @ [%d,%d]", r1, r2, index / maxx, index % maxx), ok);
        }
        Assert.assertTrue(sum + " != " + amplitude1, eq3.almostEqualRelativeOrAbsolute((double) sum, amplitude1));
    }
}

Example 2 with Gaussian2DFunction

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

the class Gaussian2DFunctionTest method functionComputesTargetGradient.

private void functionComputesTargetGradient(int targetParameter) {
    int gradientIndex = findGradientIndex(f1, targetParameter);
    double[] dyda = new double[f1.gradientIndices().length];
    double[] dyda2 = new double[dyda.length];
    double[] a;
    Gaussian2DFunction f1a = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    Gaussian2DFunction f1b = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    Statistics s = new Statistics();
    for (double background : testbackground) // Peak 1
    for (double amplitude1 : testamplitude1) for (double shape1 : testshape1) for (double cx1 : testcx1) for (double cy1 : testcy1) for (double[] w1 : testw1) {
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        f1.initialise(a);
        // Numerically solve gradient. 
        // Calculate the step size h to be an exact numerical representation
        final double xx = a[targetParameter];
        // Get h to minimise roundoff error
        double h = Precision.representableDelta(xx, h_);
        // Evaluate at (x+h) and (x-h)
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        a[targetParameter] = xx + h;
        f1a.initialise(a);
        a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        a[targetParameter] = xx - h;
        f1b.initialise(a);
        for (int x : testx) for (int y : testy) {
            int i = y * maxx + x;
            f1.eval(i, dyda);
            double value2 = f1a.eval(i, dyda2);
            double value3 = f1b.eval(i, dyda2);
            double gradient = (value2 - value3) / (2 * h);
            double error = DoubleEquality.relativeError(gradient, dyda2[gradientIndex]);
            s.add(error);
            Assert.assertTrue(gradient + " sign != " + dyda2[gradientIndex], (gradient * dyda2[gradientIndex]) >= 0);
            //System.out.printf("[%d,%d] %f == [%d] %f? (%g)\n", x, y, gradient,
            //		gradientIndex, dyda2[gradientIndex], error);
            //System.out.printf("[%d,%d] %f == [%d] %f?\n", x, y, gradient, gradientIndex, dyda[gradientIndex]);
            Assert.assertTrue(gradient + " != " + dyda[gradientIndex], eq.almostEqualRelativeOrAbsolute(gradient, dyda[gradientIndex]));
        }
    }
    System.out.printf("functionComputesTargetGradient %s %s (error %s +/- %s)\n", f1.getClass().getSimpleName(), f1.getName(targetParameter), Utils.rounded(s.getMean()), Utils.rounded(s.getStandardDeviation()));
}

Example 3 with Gaussian2DFunction

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

the class SCMOSLikelihoodWrapperTest method functionComputesGradient.

private void functionComputesGradient(int flags) {
    Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, null);
    // Setup
    testbackground = testbackground_;
    testsignal1 = testsignal1_;
    testangle1 = testangle1_;
    testcx1 = testcx1_;
    testcy1 = testcy1_;
    testw1 = testw1_;
    if (!f1.evaluatesBackground()) {
        testbackground = new double[] { testbackground[0] };
    }
    if (!f1.evaluatesSignal()) {
        testsignal1 = new double[] { testsignal1[0] };
    }
    if (!f1.evaluatesShape()) {
        testangle1 = new double[] { 0 };
    }
    // Position is always evaluated
    boolean noSecondWidth = false;
    if (!f1.evaluatesSD0()) {
        // Just use 1 width
        testw1 = new double[][] { testw1[0] };
        // If no width 0 then assume we have no width 1 as well
        noSecondWidth = true;
    } else if (!f1.evaluatesSD1()) {
        // No evaluation of second width needs only variation in width 0 so truncate 
        testw1 = Arrays.copyOf(testw1, 2);
        noSecondWidth = true;
    }
    if (noSecondWidth) {
        for (int i = 0; i < testw1.length; i++) {
            testw1[i][1] = testw1[i][0];
        }
    }
    double fraction = 90;
    if (f1.evaluatesBackground())
        functionComputesTargetGradient(f1, Gaussian2DFunction.BACKGROUND, fraction);
    if (f1.evaluatesSignal())
        functionComputesTargetGradient(f1, Gaussian2DFunction.SIGNAL, fraction);
    if (f1.evaluatesShape())
        functionComputesTargetGradient(f1, Gaussian2DFunction.SHAPE, fraction);
    functionComputesTargetGradient(f1, Gaussian2DFunction.X_POSITION, fraction);
    functionComputesTargetGradient(f1, Gaussian2DFunction.Y_POSITION, fraction);
    if (f1.evaluatesSD0())
        functionComputesTargetGradient(f1, Gaussian2DFunction.X_SD, fraction);
    if (f1.evaluatesSD1())
        functionComputesTargetGradient(f1, Gaussian2DFunction.Y_SD, fraction);
}

Example 4 with Gaussian2DFunction

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

the class SCMOSLikelihoodWrapperTest method functionComputesGradientPerDatum.

private void functionComputesGradientPerDatum(int flags) {
    Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, null);
    // Setup
    testbackground = testbackground_;
    testsignal1 = testsignal1_;
    testangle1 = testangle1_;
    testcx1 = testcx1_;
    testcy1 = testcy1_;
    testw1 = testw1_;
    if (!f1.evaluatesBackground()) {
        testbackground = new double[] { testbackground[0] };
    }
    if (!f1.evaluatesSignal()) {
        testsignal1 = new double[] { testsignal1[0] };
    }
    if (!f1.evaluatesShape()) {
        testangle1 = new double[] { 0 };
    }
    // Position is always evaluated
    boolean noSecondWidth = false;
    if (!f1.evaluatesSD0()) {
        // Just use 1 width
        testw1 = new double[][] { testw1[0] };
        // If no width 0 then assume we have no width 1 as well
        noSecondWidth = true;
    } else if (!f1.evaluatesSD1()) {
        // No evaluation of second width needs only variation in width 0 so truncate 
        testw1 = Arrays.copyOf(testw1, 2);
        noSecondWidth = true;
    }
    if (noSecondWidth) {
        for (int i = 0; i < testw1.length; i++) {
            testw1[i][1] = testw1[i][0];
        }
    }
    if (f1.evaluatesBackground())
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.BACKGROUND);
    if (f1.evaluatesSignal())
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.SIGNAL);
    if (f1.evaluatesShape())
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.SHAPE);
    functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.X_POSITION);
    functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.Y_POSITION);
    if (f1.evaluatesSD0())
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.X_SD);
    if (f1.evaluatesSD1())
        functionComputesTargetGradientPerDatum(f1, Gaussian2DFunction.Y_SD);
}

Example 5 with Gaussian2DFunction

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

the class SpeedTest method f1ComputesSameAsf2.

void f1ComputesSameAsf2(int npeaks, int flags1, int flags2) {
    DoubleEquality eq = new DoubleEquality(1e-2, 1e-10);
    int iter = 2000;
    ArrayList<double[]> paramsList2 = (npeaks == 1) ? copyList(paramsListSinglePeak, iter) : copyList(paramsListDoublePeak, iter);
    Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags1, null);
    Gaussian2DFunction f2 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, flags2, null);
    double[] dyda1 = new double[1 + npeaks * 6];
    double[] dyda2 = new double[1 + npeaks * 6];
    int[] gradientIndices = f1.gradientIndices();
    int[] g1 = new int[gradientIndices.length];
    int[] g2 = new int[gradientIndices.length];
    int nparams = 0;
    for (int i = 0; i < gradientIndices.length; i++) {
        int index1 = f1.findGradientIndex(g1[i]);
        int index2 = f2.findGradientIndex(g2[i]);
        if (index1 >= 0 && index2 >= 0) {
            g1[nparams] = index1;
            g2[nparams] = index2;
            nparams++;
        }
    }
    for (int i = 0; i < paramsList2.size(); i++) {
        f1.initialise(paramsList2.get(i));
        f2.initialise(paramsList2.get(i));
        for (int j = 0; j < x.length; j++) {
            double y1 = f1.eval(x[j], dyda1);
            double y2 = f2.eval(x[j], dyda2);
            Assert.assertTrue("Not same y[" + j + "] @ " + i + " " + y1 + " != " + y2, eq.almostEqualRelativeOrAbsolute(y1, y2));
            for (int ii = 0; ii < nparams; ii++) Assert.assertTrue("Not same dyda[" + j + "] @ " + gradientIndices[g1[ii]] + ": " + dyda1[g1[ii]] + " != " + dyda2[g2[ii]], eq.almostEqualRelativeOrAbsolute(dyda1[g1[ii]], dyda2[g2[ii]]));
        }
    }
}