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

Example 16 with Gaussian2DFunction

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

the class PoissonLikelihoodWrapperTest method functionComputesGradientPerDatum.

private void functionComputesGradientPerDatum(int flags) {
    Gaussian2DFunction f1 = GaussianFunctionFactory.create2D(1, maxx, maxx, flags, null);
    // Setup
    testbackground = testbackground_;
    testamplitude1 = testamplitude1_;
    testangle1 = testangle1_;
    testcx1 = testcx1_;
    testcy1 = testcy1_;
    testw1 = testw1_;
    if (!f1.evaluatesBackground()) {
        testbackground = new double[] { testbackground[0] };
    }
    if (!f1.evaluatesSignal()) {
        testamplitude1 = new double[] { testamplitude1[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 17 with Gaussian2DFunction

use of 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");
    PeakFit fitter = createFitter();
    // Use the fit configuration to generate a Gaussian function to test what is being evaluated
    Gaussian2DFunction gf = config.getFitConfiguration().createGaussianFunction(1, 1, 1, new double[] { 0, 10, initialPeakAngle, 0, 0, initialPeakStdDev0, initialPeakStdDev1 });
    createResultsWindow();
    int iteration = 0;
    ignore[ANGLE] = !gf.evaluatesShape();
    ignore[X] = !gf.evaluatesSD0();
    ignore[Y] = !gf.evaluatesSD1();
    double[] params = new double[] { gf.evaluatesShape() ? initialPeakAngle : 0, gf.evaluatesSD0() ? initialPeakStdDev0 : 0, gf.evaluatesSD1() ? initialPeakStdDev1 : 0, 0, 0 };
    double[] params_dev = new double[3];
    boolean[] identical = new boolean[4];
    double[] p = new double[] { Double.NaN, Double.NaN, Double.NaN, Double.NaN };
    addToResultTable(iteration++, 0, params, params_dev, p);
    if (!calculateStatistics(fitter, params, params_dev))
        return (Utils.isInterrupted()) ? ABORTED : INSUFFICIENT_PEAKS;
    if (!addToResultTable(iteration++, size(), params, params_dev, p))
        return BAD_ESTIMATE;
    boolean tryAgain = false;
    do {
        if (!calculateStatistics(fitter, params, params_dev))
            return (Utils.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(iteration++, size(), params, params_dev, 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.updatePreferences) {
                    config.getFitConfiguration().setInitialPeakStdDev0((float) params[X]);
                    config.getFitConfiguration().setInitialPeakStdDev1((float) params[Y]);
                    config.getFitConfiguration().setInitialAngle((float) params[ANGLE]);
                }
                break;
            }
            if (IJ.escapePressed()) {
                IJ.beep();
                IJ.showStatus("Aborted");
                return ABORTED;
            }
        } catch (Exception e) {
            e.printStackTrace();
            return EXCEPTION;
        }
    } while (true);
    return (tryAgain) ? TRY_AGAIN : COMPLETE;
}

Example 18 with Gaussian2DFunction

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

the class FitConfiguration method createFunctionSolver.

private BaseFunctionSolver createFunctionSolver() {
    if (gaussianFunction == null) {
        // Other code may want to call getFunctionSolver() to see if exceptions are thrown
        // so create a dummy function so we can return a function solver.
        gaussianFunction = createGaussianFunction(1, 1, 1, null);
    }
    // Remove noise model
    gaussianFunction.setNoiseModel(null);
    NonLinearFit nlinfit;
    switch(fitSolver) {
        case MLE:
            // Only the Poisson likelihood function supports gradients
            if (searchMethod.usesGradients() && modelCamera) {
                throw new IllegalArgumentException(String.format("The derivative based search method '%s' can only be used with the " + "'%s' likelihood function, i.e. no model camera noise", searchMethod, MaximumLikelihoodFitter.LikelihoodFunction.POISSON));
            }
            MaximumLikelihoodFitter fitter = new MaximumLikelihoodFitter(gaussianFunction);
            fitter.setRelativeThreshold(relativeThreshold);
            fitter.setAbsoluteThreshold(absoluteThreshold);
            fitter.setMaxEvaluations(maxFunctionEvaluations);
            fitter.setMaxIterations(maxIterations);
            fitter.setSearchMethod(searchMethod);
            fitter.setGradientLineMinimisation(gradientLineMinimisation);
            // Specify the likelihood function to use
            if (modelCamera) {
                // Set the camera read noise
                fitter.setSigma(readNoise);
                if (emCCD) {
                    // EMCCD = Poisson+Gamma+Gaussian
                    fitter.setLikelihoodFunction(MaximumLikelihoodFitter.LikelihoodFunction.POISSON_GAMMA_GAUSSIAN);
                } else {
                    // CCD = Poisson+Gaussian
                    fitter.setLikelihoodFunction(MaximumLikelihoodFitter.LikelihoodFunction.POISSON_GAUSSIAN);
                }
            } else {
                fitter.setLikelihoodFunction(MaximumLikelihoodFitter.LikelihoodFunction.POISSON);
            }
            // All models use the amplification gain (i.e. how many ADUs/electron)
            if (amplification <= 0) {
                throw new IllegalArgumentException("The amplification is required for the " + fitSolver.getName());
            }
            fitter.setAlpha(1.0 / amplification);
            return fitter;
        case BOUNDED_LVM_WEIGHTED:
            gaussianFunction.setNoiseModel(getNoiseModel());
        case BOUNDED_LVM:
        case LVM_MLE:
            if (fitSolver == FitSolver.LVM_MLE && gain <= 0) {
                throw new IllegalArgumentException("The gain is required for the " + fitSolver.getName());
            }
            if (useClamping) {
                bounds = new ParameterBounds(gaussianFunction);
                setClampValues(bounds);
            }
            BoundedNonLinearFit bnlinfit = new BoundedNonLinearFit(gaussianFunction, getStoppingCriteria(), bounds);
            nlinfit = bnlinfit;
            break;
        case LVM_QUASI_NEWTON:
            // This only works with a Gaussian2DFunction
            if (gaussianFunction instanceof Gaussian2DFunction) {
                ApacheLVMFitter apacheNLinFit = new ApacheLVMFitter((Gaussian2DFunction) gaussianFunction);
                apacheNLinFit.setMaxEvaluations(maxIterations);
                // TODO - Configure stopping criteria ...
                return apacheNLinFit;
            }
        case LVM_WEIGHTED:
        case LVM:
        default:
            // Only set the weighting function if necessary
            if (fitSolver == FitSolver.LVM_WEIGHTED)
                gaussianFunction.setNoiseModel(getNoiseModel());
            nlinfit = new NonLinearFit(gaussianFunction, getStoppingCriteria());
    }
    nlinfit.setInitialLambda(getLambda());
    if (fitSolver == FitSolver.LVM_MLE) {
        nlinfit.setMLE(true);
    }
    return nlinfit;
}

Example 19 with Gaussian2DFunction

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

the class BaseFunctionSolverTest method drawGaussian.

/**
	 * Draw a Gaussian with Poisson shot noise and Gaussian read noise.
	 *
	 * @param params
	 *            The Gaussian parameters
	 * @param noise
	 *            The read noise
	 * @param noiseModel
	 *            the noise model
	 * @return The data
	 */
double[] drawGaussian(double[] params, double[] noise, NoiseModel noiseModel) {
    double[] data = new double[size * size];
    int n = params.length / 6;
    Gaussian2DFunction f = GaussianFunctionFactory.create2D(n, size, size, flags, null);
    f.initialise(params);
    // Poisson noise
    for (int i = 0; i < data.length; i++) {
        CustomPoissonDistribution dist = new CustomPoissonDistribution(randomGenerator, 1);
        double e = f.eval(i);
        if (e > 0) {
            dist.setMeanUnsafe(e);
            data[i] = dist.sample();
        }
    }
    // Simulate EM-gain
    if (noiseModel == NoiseModel.EMCCD) {
        // Use a gamma distribution
        // Since the call random.nextGamma(...) creates a Gamma distribution 
        // which pre-calculates factors only using the scale parameter we 
        // create a custom gamma distribution where the shape can be set as a property.
        CustomGammaDistribution dist = new CustomGammaDistribution(randomGenerator, 1, emGain);
        for (int i = 0; i < data.length; i++) {
            if (data[i] > 0) {
                dist.setShapeUnsafe(data[i]);
                // The sample will amplify the signal so we remap to the original scale
                data[i] = dist.sample() / emGain;
            }
        }
    }
    // Read-noise
    if (noise != null) {
        for (int i = 0; i < data.length; i++) {
            data[i] += randomGenerator.nextGaussian() * noise[i];
        }
    }
    //gdsc.core.ij.Utils.display("Spot", data, size, size);
    return data;
}

Example 20 with Gaussian2DFunction

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

the class ErfGaussian2DFunctionTest method functionIsFasterThanEquivalentGaussian2DFunction.

// Speed test verses equivalent Gaussian2DFunction
@Test
public void functionIsFasterThanEquivalentGaussian2DFunction() {
    int flags = this.flags & ~GaussianFunctionFactory.FIT_ERF;
    final Gaussian2DFunction gf = GaussianFunctionFactory.create2D(1, maxx, maxy, flags, zModel);
    boolean zDepth = (flags & GaussianFunctionFactory.FIT_Z) != 0;
    final TurboList<double[]> params = new TurboList<double[]>();
    final TurboList<double[]> params2 = new TurboList<double[]>();
    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) {
        double[] a = createParameters(background, amplitude1, shape1, cx1, cy1, w1[0], w1[1]);
        params.add(a);
        if (zDepth) {
            // Change to a standard free circular function
            a = a.clone();
            a[Gaussian2DFunction.X_SD] *= zModel.getSx(a[Gaussian2DFunction.SHAPE]);
            a[Gaussian2DFunction.Y_SD] *= zModel.getSy(a[Gaussian2DFunction.SHAPE]);
            a[Gaussian2DFunction.SHAPE] = 0;
            params2.add(a);
        }
    }
    double[][] x = params.toArray(new double[params.size()][]);
    double[][] x2 = (zDepth) ? params2.toArray(new double[params2.size()][]) : x;
    int runs = 10000 / x.length;
    TimingService ts = new TimingService(runs);
    ts.execute(new FunctionTimingTask(gf, x2, 1));
    ts.execute(new FunctionTimingTask(gf, x2, 0));
    ts.execute(new FunctionTimingTask(f1, x, 2));
    ts.execute(new FunctionTimingTask(f1, x, 1));
    ts.execute(new FunctionTimingTask(f1, x, 0));
    int size = ts.getSize();
    ts.repeat(size);
    ts.report();
    // Sometimes this fails, probably due to JVM optimisations, so skip for now
    if (!TestSettings.ASSERT_SPEED_TESTS)
        return;
    int n = ts.getSize() - 1;
    Assert.assertTrue("0 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
    n--;
    Assert.assertTrue("1 order", ts.get(n).getMean() < ts.get(n - 3).getMean());
}