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

Example 21 with Gaussian2DFunction

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

the class EjmlLinearSolverTest method runSolverSpeedTest.

private void runSolverSpeedTest(RandomSeed seed, int flags) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    final int n = f0.size();
    final double[] y = new double[n];
    final LocalList<DenseMatrix64F> aList = new LocalList<>();
    final LocalList<DenseMatrix64F> bList = new LocalList<>();
    final double[] testbackground = new double[] { 0.2, 0.7 };
    final double[] testsignal1 = new double[] { 30, 100, 300 };
    final double[] testcx1 = new double[] { 4.9, 5.3 };
    final double[] testcy1 = new double[] { 4.8, 5.2 };
    final double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    final int np = f0.getNumberOfGradients();
    final GradientCalculator calc = GradientCalculatorUtils.newCalculator(np);
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    // double lambda = 10;
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double w1 : testw1) {
                        final double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
                        f0.initialise(p);
                        f0.forEach(new ValueProcedure() {

                            int index = 0;

                            @Override
                            public void execute(double value) {
                                // Poisson data
                                y[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
                            }
                        });
                        final double[][] alpha = new double[np][np];
                        final double[] beta = new double[np];
                        // double ss =
                        calc.findLinearised(n, y, p, alpha, beta, f0);
                        // TestLog.fine(logger,"SS = %f", ss);
                        // As per the LVM algorithm
                        // for (int i = 0; i < np; i++)
                        // alpha[i][i] *= lambda;
                        aList.add(EjmlLinearSolver.toA(alpha));
                        bList.add(EjmlLinearSolver.toB(beta));
                    }
                }
            }
        }
    }
    final DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[0]);
    final DenseMatrix64F[] b = bList.toArray(new DenseMatrix64F[0]);
    final int runs = 100000 / a.length;
    final TimingService ts = new TimingService(runs);
    final LocalList<SolverTimingTask> tasks = new LocalList<>();
    // Added in descending speed order
    tasks.add(new PseudoInverseSolverTimingTask(a, b));
    tasks.add(new LinearSolverTimingTask(a, b));
    tasks.add(new CholeskySolverTimingTask(a, b));
    tasks.add(new CholeskyLdltSolverTimingTask(a, b));
    tasks.add(new DirectInversionSolverTimingTask(a, b));
    for (final SolverTimingTask task : tasks) {
        if (!task.badSolver) {
            ts.execute(task);
        }
    }
    final int size = ts.getSize();
    ts.repeat();
    if (logger.isLoggable(Level.INFO)) {
        logger.info(ts.getReport(size));
    }
    // Just check the PseudoInverse is slowest
    for (int i = 1; i < size; i++) {
        logger.log(TestLogUtils.getTimingRecord(ts.get(-(size)), ts.get(-i)));
    }
    if (np > 2) {
        // The Direct solver may not be faster at size=5
        int i = (np == 5) ? 2 : 1;
        final int size_1 = size - 1;
        for (; i < size_1; i++) {
            logger.log(TestLogUtils.getTimingRecord(ts.get(-(size_1)), ts.get(-i)));
        }
    }
}

Example 22 with Gaussian2DFunction

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

the class EjmlLinearSolverTest method runInversionSpeedTest.

private void runInversionSpeedTest(RandomSeed seed, int flags) {
    Assumptions.assumeTrue(TestSettings.allow(TestComplexity.MEDIUM));
    final Gaussian2DFunction f0 = GaussianFunctionFactory.create2D(1, 10, 10, flags, null);
    final int n = f0.size();
    final double[] y = new double[n];
    final LocalList<DenseMatrix64F> aList = new LocalList<>();
    final double[] testbackground = new double[] { 0.2, 0.7 };
    final double[] testsignal1 = new double[] { 30, 100, 300 };
    final double[] testcx1 = new double[] { 4.9, 5.3 };
    final double[] testcy1 = new double[] { 4.8, 5.2 };
    final double[] testw1 = new double[] { 1.1, 1.2, 1.5 };
    final int np = f0.getNumberOfGradients();
    final GradientCalculator calc = GradientCalculatorUtils.newCalculator(np);
    final UniformRandomProvider rng = RngUtils.create(seed.getSeed());
    // double lambda = 10;
    for (final double background : testbackground) {
        // Peak 1
        for (final double signal1 : testsignal1) {
            for (final double cx1 : testcx1) {
                for (final double cy1 : testcy1) {
                    for (final double w1 : testw1) {
                        final double[] p = new double[] { background, signal1, 0, cx1, cy1, w1, w1 };
                        f0.initialise(p);
                        f0.forEach(new ValueProcedure() {

                            int index = 0;

                            @Override
                            public void execute(double value) {
                                // Poisson data
                                y[index++] = GdscSmlmTestUtils.createPoissonSampler(rng, value).sample();
                            }
                        });
                        final double[][] alpha = new double[np][np];
                        final double[] beta = new double[np];
                        // double ss =
                        calc.findLinearised(n, y, p, alpha, beta, f0);
                        // TestLog.fine(logger,"SS = %f", ss);
                        // As per the LVM algorithm
                        // for (int i = 0; i < np; i++)
                        // alpha[i][i] *= lambda;
                        aList.add(EjmlLinearSolver.toA(alpha));
                    }
                }
            }
        }
    }
    final DenseMatrix64F[] a = aList.toArray(new DenseMatrix64F[0]);
    final boolean[] ignore = new boolean[a.length];
    final double[][] answer = new double[a.length][];
    final int runs = 100000 / a.length;
    final TimingService ts = new TimingService(runs);
    final LocalList<InversionTimingTask> tasks = new LocalList<>();
    // Added in descending speed order
    tasks.add(new PseudoInverseInversionTimingTask(a, ignore, answer));
    tasks.add(new LinearInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyLdltInversionTimingTask(a, ignore, answer));
    tasks.add(new CholeskyInversionTimingTask(a, ignore, answer));
    tasks.add(new DirectInversionInversionTimingTask(a, ignore, answer));
    tasks.add(new DiagonalDirectInversionInversionTimingTask(a, ignore, answer));
    for (final InversionTimingTask task : tasks) {
        if (!task.badSolver) {
            ts.execute(task);
        }
    }
    final int size = ts.getSize();
    ts.repeat();
    if (logger.isLoggable(Level.INFO)) {
        logger.info(ts.getReport(size));
    }
    // When it is present the DiagonalDirect is fastest (n<=5)
    if (np <= 5) {
        for (int i = 2; i <= size; i++) {
            logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-1)));
        }
        if (np < 5) {
            // n < 5 Direct is fastest
            for (int i = 3; i <= size; i++) {
                logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-2)));
            }
        } else {
            // and may not be faster than Direct at n=5 so that comparison is ignored.
            for (int i = 4; i <= size; i++) {
                logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-3)));
            }
        }
    } else {
        // Cholesky should be fastest.
        for (int i = 2; i <= size; i++) {
            logger.log(TestLogUtils.getTimingRecord(ts.get(-i), ts.get(-1)));
        }
    }
}

Example 23 with Gaussian2DFunction

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

the class PoissonGradientProcedureTest method varianceMatchesFormula.

@SeededTest
void varianceMatchesFormula() {
    // Assumptions.assumeTrue(false);
    final double[] n_ = new double[] { 20, 50, 100, 500 };
    final double[] b2_ = new double[] { 0, 1, 2, 4 };
    final double[] s_ = new double[] { 1, 1.2, 1.5 };
    final double[] x_ = new double[] { 4.8, 5, 5.5 };
    final double a = 100;
    final int size = 10;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final PoissonGradientProcedure p = PoissonGradientProcedureUtils.create(f);
    final int ix = f.findGradientIndex(Gaussian2DFunction.X_POSITION);
    final int iy = f.findGradientIndex(Gaussian2DFunction.Y_POSITION);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    for (final double n : n_) {
        params[Gaussian2DFunction.SIGNAL] = n;
        for (final double b2 : b2_) {
            params[Gaussian2DFunction.BACKGROUND] = b2;
            for (final double s : s_) {
                final double ss = s * a;
                params[Gaussian2DFunction.X_SD] = s;
                for (final double x : x_) {
                    params[Gaussian2DFunction.X_POSITION] = x;
                    for (final double y : x_) {
                        params[Gaussian2DFunction.Y_POSITION] = y;
                        p.computeFisherInformation(params);
                        final FisherInformationMatrix m1 = new FisherInformationMatrix(p.getLinear(), p.numberOfGradients);
                        final double[] crlb = m1.crlb();
                        if (crlb == null) {
                            Assertions.fail("No variance");
                        }
                        @SuppressWarnings("null") final double o1 = Math.sqrt(crlb[ix]) * a;
                        final double o2 = Math.sqrt(crlb[iy]) * a;
                        final double e = Gaussian2DPeakResultHelper.getMLPrecisionX(a, ss, n, b2, false);
                        // logger.fine(FunctionUtils.getSupplier("e = %f : o = %f %f", e, o1, o2);
                        Assertions.assertEquals(e, o1, e * 5e-2);
                        Assertions.assertEquals(e, o2, e * 5e-2);
                    }
                }
            }
        }
    }
}

Example 24 with Gaussian2DFunction

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

the class FunctionHelperTest method canGetMeanValueForGaussian.

@Test
void canGetMeanValueForGaussian() {
    final float intensity = 100;
    // Realistic standard deviations.
    // Only test the highest
    final float[] s = { 1, 1.2f, 1.5f, 2, 2.5f };
    final int n_1 = s.length - 1;
    // Flag to indicate that all levels should be run and the difference reported
    final boolean debug = false;
    final Logger logger = (debug) ? Logger.getLogger(FunctionHelperTest.class.getName()) : null;
    for (int i = (debug) ? 0 : n_1; i < s.length; i++) {
        final float sx = s[i];
        for (int j = i; j < s.length; j++) {
            final float sy = s[j];
            final int size = 1 + 2 * (int) Math.ceil(Math.max(sx, sy) * 4);
            final float[] a = Gaussian2DPeakResultHelper.createParams(0.f, intensity, size / 2f, size / 2f, 0.f, sx, sy, 0);
            final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, size, size, GaussianFunctionFactory.FIT_FREE_CIRCLE, null);
            final double[] values = f.computeValues(SimpleArrayUtils.toDouble(a));
            // ImagePlus imp = new ImagePlus("gauss", new FloatProcessor(size, size, values));
            // double cx = size / 2.;
            // Shape shape = new Ellipse2D.Double(cx - sx, cx - sy, 2 * sx, 2 * sy);
            // imp.setRoi(new ShapeRoi(shape));
            // IJ.save(imp, "/Users/ah403/1.tif");
            final double scale = MathUtils.sum(values) / intensity;
            for (int range = 1; range <= 3; range++) {
                final double e = Gaussian2DPeakResultHelper.getMeanSignalUsingR(intensity, sx, sy, range);
                final double o = FunctionHelper.getMeanValue(values.clone(), scale * Gaussian2DPeakResultHelper.cumulative2D(range));
                if (debug) {
                    logger.fine(FunctionUtils.getSupplier("%g,%g   %d  %g %g  %g", sx, sy, range, e, o, DoubleEquality.relativeError(e, o)));
                }
                // Only test the highest
                if (i == n_1) {
                    Assertions.assertEquals(e, o, e * 0.025);
                }
            }
        }
    }
}

Example 25 with Gaussian2DFunction

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

the class PsfEstimator method estimatePsf.

private int estimatePsf() {
    log("Estimating PSF ... Press escape to abort");
    final PeakFit fitter = createFitter();
    // Use the fit configuration to generate a Gaussian function to test what is being evaluated
    final Gaussian2DFunction gf = config.getFitConfiguration().createGaussianFunction(1, 1, 1);
    ignore[ANGLE] = !gf.evaluatesAngle();
    ignore[X] = !gf.evaluatesSD0();
    ignore[Y] = !gf.evaluatesSD1();
    final double[] params = new double[] { gf.evaluatesAngle() ? initialPeakAngle : 0, gf.evaluatesSD0() ? initialPeakStdDev0 : 0, gf.evaluatesSD1() ? initialPeakStdDev1 : 0, 0, 0 };
    final double[] paramsDev = new double[3];
    final boolean[] identical = new boolean[4];
    final double[] p = new double[] { Double.NaN, Double.NaN, Double.NaN, Double.NaN };
    final TextWindow resultsWindow = createResultsWindow();
    addToResultTable(resultsWindow, 0, 0, params, paramsDev, p);
    if (!calculateStatistics(fitter, params, paramsDev)) {
        return (ImageJUtils.isInterrupted()) ? ABORTED : INSUFFICIENT_PEAKS;
    }
    if (!addToResultTable(resultsWindow, 1, size(), params, paramsDev, p)) {
        return BAD_ESTIMATE;
    }
    boolean tryAgain = false;
    int iteration = 2;
    for (; ; ) {
        if (!calculateStatistics(fitter, params, paramsDev)) {
            return (ImageJUtils.isInterrupted()) ? ABORTED : INSUFFICIENT_PEAKS;
        }
        try {
            for (int i = 0; i < 3; i++) {
                getP(i, p, identical);
            }
            if (!ignore[Y]) {
                getPairedP(sampleNew[X], sampleNew[Y], XY, p, identical);
            }
            if (!addToResultTable(resultsWindow, iteration++, size(), params, paramsDev, p)) {
                return BAD_ESTIMATE;
            }
            if ((ignore[ANGLE] || identical[ANGLE] || identical[XY]) && (ignore[X] || identical[X]) && (ignore[Y] || identical[Y])) {
                tryAgain = checkAngleSignificance() || checkXySignificance(identical);
                // Update recommended values. Only use if significant
                params[X] = sampleNew[X].getMean();
                params[Y] = (!ignore[Y] && !identical[XY]) ? sampleNew[Y].getMean() : params[X];
                params[ANGLE] = (!ignore[ANGLE]) ? sampleNew[ANGLE].getMean() : 0;
                // update starting configuration
                initialPeakAngle = (float) params[ANGLE];
                initialPeakStdDev0 = (float) params[X];
                initialPeakStdDev1 = (float) params[Y];
                if (settings.getUpdatePreferences()) {
                    config.getFitConfiguration().setInitialPeakStdDev0((float) params[X]);
                    try {
                        config.getFitConfiguration().setInitialPeakStdDev1((float) params[Y]);
                        config.getFitConfiguration().setInitialAngle((float) params[ANGLE]);
                    } catch (final IllegalStateException | ConfigurationException ex) {
                    // Ignore this as the current PSF is not a 2 axis and theta Gaussian PSF
                    }
                }
                break;
            }
            if (IJ.escapePressed()) {
                IJ.beep();
                IJ.showStatus("Aborted");
                return ABORTED;
            }
        } catch (final Exception ex) {
            Logger.getLogger(getClass().getName()).log(Level.WARNING, "Error while estimating the PSF", ex);
            return EXCEPTION;
        }
    }
    return (tryAgain) ? TRY_AGAIN : COMPLETE;
}