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

Example 31 with Gaussian2DFunction

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

the class FisherInformationMatrixTest method createFisherInformationMatrix.

private static FisherInformationMatrix createFisherInformationMatrix(UniformRandomProvider rg, int columns, int zeroColumns) {
    final int maxx = 10;
    final int size = maxx * maxx;
    // Use a real Gaussian function here to compute the Fisher information.
    // The matrix may be sensitive to the type of equation used.
    int npeazeroColumnss = 1;
    Gaussian2DFunction fun = createFunction(maxx, npeazeroColumnss);
    while (fun.getNumberOfGradients() < columns) {
        npeazeroColumnss++;
        fun = createFunction(maxx, npeazeroColumnss);
    }
    final double[] a = new double[1 + npeazeroColumnss * Gaussian2DFunction.PARAMETERS_PER_PEAK];
    a[Gaussian2DFunction.BACKGROUND] = nextUniform(rg, 1, 5);
    for (int i = 0, j = 0; i < npeazeroColumnss; i++, j += Gaussian2DFunction.PARAMETERS_PER_PEAK) {
        a[j + Gaussian2DFunction.SIGNAL] = nextUniform(rg, 100, 300);
        // Non-overlapping peazeroColumnss otherwise the Crlb are poor
        a[j + Gaussian2DFunction.X_POSITION] = nextUniform(rg, 2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.Y_POSITION] = nextUniform(rg, 2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.X_SD] = nextUniform(rg, 1.5, 2);
        a[j + Gaussian2DFunction.Y_SD] = nextUniform(rg, 1.5, 2);
    }
    fun.initialise(a);
    final GradientCalculator calc = GradientCalculatorUtils.newCalculator(fun.getNumberOfGradients());
    double[][] matrixI = calc.fisherInformationMatrix(size, a, fun);
    // Reduce to the desired size
    matrixI = Arrays.copyOf(matrixI, columns);
    for (int i = 0; i < columns; i++) {
        matrixI[i] = Arrays.copyOf(matrixI[i], columns);
    }
    // Zero selected columns
    if (zeroColumns > 0) {
        final int[] zero = RandomUtils.sample(zeroColumns, columns, rg);
        for (final int i : zero) {
            for (int j = 0; j < columns; j++) {
                matrixI[i][j] = matrixI[j][i] = 0;
            }
        }
    }
    // Create matrix
    return new FisherInformationMatrix(matrixI, 1e-3);
}

Example 32 with Gaussian2DFunction

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

the class UnivariateLikelihoodFisherInformationCalculatorTest method computePoissonFisherInformation.

private static void computePoissonFisherInformation(UniformRandomProvider rng, Model model) {
    // Create function
    final Gaussian2DFunction func = GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_CIRCLE, null);
    final double[] params = new double[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    params[Gaussian2DFunction.BACKGROUND] = nextUniform(rng, 0.1, 0.3);
    params[Gaussian2DFunction.SIGNAL] = nextUniform(rng, 100, 300);
    params[Gaussian2DFunction.X_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.Y_POSITION] = nextUniform(rng, 4, 6);
    params[Gaussian2DFunction.X_SD] = nextUniform(rng, 1, 1.3);
    Gradient1Function f1 = func;
    FisherInformation fi;
    switch(model) {
        // Get a variance
        case POISSON_GAUSSIAN:
            final double var = 0.9 + 0.2 * rng.nextDouble();
            fi = new PoissonGaussianApproximationFisherInformation(Math.sqrt(var));
            f1 = (Gradient1Function) OffsetFunctionFactory.wrapFunction(func, SimpleArrayUtils.newDoubleArray(func.size(), var));
            break;
        case POISSON:
            fi = new PoissonFisherInformation();
            break;
        case HALF_POISSON:
            fi = new HalfPoissonFisherInformation();
            break;
        default:
            throw new IllegalStateException();
    }
    // This introduces a dependency on a different package, and relies on that
    // computing the correct answer. However that code predates this and so the
    // test ensures that the FisherInformationCalculator functions correctly.
    final PoissonGradientProcedure p1 = PoissonGradientProcedureUtils.create(f1);
    p1.computeFisherInformation(params);
    final double[] e = p1.getLinear();
    final FisherInformationCalculator calc = new UnivariateLikelihoodFisherInformationCalculator(func, fi);
    final FisherInformationMatrix I = calc.compute(params);
    final double[] o = I.getMatrix().data;
    final boolean emCcd = model == Model.HALF_POISSON;
    if (emCcd) {
        // Assumes half the poisson fisher information
        SimpleArrayUtils.multiply(e, 0.5);
    }
    Assertions.assertArrayEquals(e, o, 1e-6);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(5e-2, 0);
    if (model == Model.POISSON || model == Model.HALF_POISSON) {
        // Get the Mortensen approximation for fitting Poisson data with a Gaussian.
        // Set a to 100 for the square pixel adjustment.
        final double a = 100;
        final double s = params[Gaussian2DFunction.X_SD] * a;
        final double N = params[Gaussian2DFunction.SIGNAL];
        final double b2 = params[Gaussian2DFunction.BACKGROUND];
        double var = Gaussian2DPeakResultHelper.getMLVarianceX(a, s, N, b2, emCcd);
        // Convert expected variance to pixels
        var /= (a * a);
        // Get the limits by inverting the Fisher information
        final double[] crlb = I.crlb();
        TestAssertions.assertTest(var, crlb[2], predicate);
        TestAssertions.assertTest(var, crlb[3], predicate);
    }
}

Example 33 with Gaussian2DFunction

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

the class Gaussian2DPeakResultHelperTest method canComputePixelAmplitude.

@Test
void canComputePixelAmplitude() {
    final float[] x = new float[] { 0f, 0.1f, 0.3f, 0.5f, 0.7f, 1f };
    final float[] s = new float[] { 0.8f, 1f, 1.5f, 2.2f };
    final float[] paramsf = new float[1 + Gaussian2DFunction.PARAMETERS_PER_PEAK];
    paramsf[Gaussian2DFunction.BACKGROUND] = 0;
    paramsf[Gaussian2DFunction.SIGNAL] = 105;
    final Gaussian2DFunction f = GaussianFunctionFactory.create2D(1, 1, 1, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    final SimpleRegression r = new SimpleRegression(false);
    for (final float tx : x) {
        for (final float ty : x) {
            for (final float sx : s) {
                for (final float sy : s) {
                    paramsf[Gaussian2DFunction.X_POSITION] = tx;
                    paramsf[Gaussian2DFunction.Y_POSITION] = ty;
                    paramsf[Gaussian2DFunction.X_SD] = sx;
                    paramsf[Gaussian2DFunction.Y_SD] = sy;
                    // Get the answer using a single pixel image
                    // Note the Gaussian2D functions set the centre of the pixel as 0,0 so offset
                    final double[] params = SimpleArrayUtils.toDouble(paramsf);
                    params[Gaussian2DFunction.X_POSITION] -= 0.5;
                    params[Gaussian2DFunction.Y_POSITION] -= 0.5;
                    f.initialise0(params);
                    final double e = f.eval(0);
                    final PSF psf = PsfHelper.create(PSFType.TWO_AXIS_GAUSSIAN_2D);
                    final CalibrationWriter calibration = new CalibrationWriter();
                    calibration.setCountPerPhoton(1);
                    calibration.setIntensityUnit(IntensityUnit.PHOTON);
                    calibration.setNmPerPixel(1);
                    calibration.setDistanceUnit(DistanceUnit.PIXEL);
                    final Gaussian2DPeakResultCalculator calc = Gaussian2DPeakResultHelper.create(psf, calibration, Gaussian2DPeakResultHelper.AMPLITUDE | Gaussian2DPeakResultHelper.PIXEL_AMPLITUDE);
                    final double o1 = calc.getAmplitude(paramsf);
                    final double o2 = calc.getPixelAmplitude(paramsf);
                    // logger.fine(FunctionUtils.getSupplier("e=%f, o1=%f, o2=%f", e, o1, o2));
                    Assertions.assertEquals(e, o2, 1e-3);
                    r.addData(e, o1);
                }
            }
        }
    }
    // logger.fine(FunctionUtils.getSupplier("Regression: pixel amplitude vs amplitude = %f,
    // slope=%f, n=%d", r.getR(), r.getSlope(),
    // r.getN()));
    // The simple amplitude over estimates the actual pixel amplitude
    Assertions.assertTrue(r.getSlope() > 1);
}