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

Example 1 with PoissonGammaGaussianFunction

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

the class CameraModelAnalysis method getLikelihoodFunction.

private static LikelihoodFunction getLikelihoodFunction(CameraModelAnalysisSettings settings) {
    final double alpha = 1.0 / getGain(settings);
    final double noise = getReadNoise(settings);
    final Model model = Model.forNumber(settings.getModel());
    switch(model) {
        case POISSON_PMF:
            return new PoissonFunction(alpha);
        case POISSON_DISRECTE:
            return new InterpolatedPoissonFunction(alpha, false);
        case POISSON_CONTINUOUS:
            return new InterpolatedPoissonFunction(alpha, true);
        case POISSON_GAUSSIAN_PDF:
        case POISSON_GAUSSIAN_PMF:
            final PoissonGaussianConvolutionFunction f1 = PoissonGaussianConvolutionFunction.createWithStandardDeviation(alpha, noise);
            f1.setComputePmf(model == Model.POISSON_GAUSSIAN_PMF);
            return f1;
        case POISSON_GAUSSIAN_APPROX:
            return PoissonGaussianFunction2.createWithStandardDeviation(alpha, noise);
        case POISSON_POISSON:
            return PoissonPoissonFunction.createWithStandardDeviation(alpha, noise);
        case POISSON_GAMMA_GAUSSIAN_PDF_CONVOLUTION:
            return PoissonGammaGaussianConvolutionFunction.createWithStandardDeviation(alpha, noise);
        case POISSON_GAMMA_PMF:
            return PoissonGammaFunction.createWithAlpha(alpha);
        case POISSON_GAMMA_GAUSSIAN_APPROX:
        case POISSON_GAMMA_GAUSSIAN_PDF_INTEGRATION:
        case POISSON_GAMMA_GAUSSIAN_PMF_INTEGRATION:
        case POISSON_GAMMA_GAUSSIAN_SIMPSON_INTEGRATION:
        case POISSON_GAMMA_GAUSSIAN_LEGENDRE_GAUSS_INTEGRATION:
            final PoissonGammaGaussianFunction f2 = new PoissonGammaGaussianFunction(alpha, noise);
            f2.setMinimumProbability(0);
            f2.setConvolutionMode(getConvolutionMode(model));
            // The function should return a PMF/PDF depending on how it is used
            f2.setPmfMode(!settings.getSimpsonIntegration());
            return f2;
        default:
            throw new IllegalStateException();
    }
}

Example 2 with PoissonGammaGaussianFunction

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

the class EmGainAnalysis method createPoissonGammaGaussianFunction.

private LikelihoodFunction createPoissonGammaGaussianFunction(double noise) {
    final PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / settings.settingGain, noise);
    fun.setMinimumProbability(0);
    return fun;
}

Example 3 with PoissonGammaGaussianFunction

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

the class EmGainAnalysis method fit.

/**
 * Fit the EM-gain distribution (Gaussian * Gamma).
 *
 * @param histogram The distribution
 */
private void fit(int[] histogram) {
    final int[] limits = limits(histogram);
    final double[] x = getX(limits);
    final double[] y = getY(histogram, limits);
    Plot plot = new Plot(TITLE, "ADU", "Frequency");
    double yMax = MathUtils.max(y);
    plot.setLimits(limits[0], limits[1], 0, yMax);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    ImageJUtils.display(TITLE, plot);
    // Estimate remaining parameters.
    // Assuming a gamma_distribution(shape,scale) then mean = shape * scale
    // scale = gain
    // shape = Photons = mean / gain
    double mean = getMean(histogram) - settings.bias;
    // Note: if the bias is too high then the mean will be negative. Just move the bias.
    while (mean < 0) {
        settings.bias -= 1;
        mean += 1;
    }
    double photons = mean / settings.gain;
    if (settings.settingSimulate) {
        ImageJUtils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.settingBias, MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    }
    ImageJUtils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    final int max = (int) x[x.length - 1];
    double[] pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot.setColor(Color.blue);
    plot.addPoints(x, pg, Plot.LINE);
    ImageJUtils.display(TITLE, plot);
    // Perform a fit
    final CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16);
    final double[] startPoint = new double[] { photons, settings.gain, settings.noise, settings.bias };
    int maxEval = 3000;
    final String[] paramNames = { "Photons", "Gain", "Noise", "Bias" };
    // Set bounds
    final double[] lower = new double[] { 0, 0.5 * settings.gain, 0, settings.bias - settings.noise };
    final double[] upper = new double[] { 2 * photons, 2 * settings.gain, settings.gain, settings.bias + settings.noise };
    // Restart until converged.
    // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data.
    PointValuePair solution = null;
    for (int iter = 0; iter < 3; iter++) {
        IJ.showStatus("Fitting histogram ... Iteration " + iter);
        try {
            // Basic Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef()));
            if (solution == null || optimum.getValue() < solution.getValue()) {
                final double[] point = optimum.getPointRef();
                // Check the bounds
                for (int i = 0; i < point.length; i++) {
                    if (point[i] < lower[i] || point[i] > upper[i]) {
                        throw new ComputationException(String.format("Fit out of of estimated range: %s %f", paramNames[i], point[i]));
                    }
                }
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
        try {
            // Bounded Powell optimiser
            final MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            final MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef())));
            final double[] point = adapter.unboundedToBounded(optimum.getPointRef());
            optimum = new PointValuePair(point, optimum.getValue());
            if (solution == null || optimum.getValue() < solution.getValue()) {
                solution = optimum;
            }
        } catch (final Exception ex) {
            IJ.log("Bounded Powell error: " + ex.getMessage());
            if (ex instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
    }
    ImageJUtils.finished();
    if (solution == null) {
        ImageJUtils.log("Failed to fit the distribution");
        return;
    }
    final double[] point = solution.getPointRef();
    photons = point[0];
    settings.gain = point[1];
    settings.noise = point[2];
    settings.bias = (int) Math.round(point[3]);
    final String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) settings.bias, MathUtils.rounded(settings.gain), MathUtils.rounded(settings.noise), MathUtils.rounded(photons));
    ImageJUtils.log(label);
    if (settings.settingSimulate) {
        ImageJUtils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", MathUtils.rounded(relativeError(settings.bias, settings.settingBias)), MathUtils.rounded(relativeError(settings.gain, settings.settingGain)), MathUtils.rounded(relativeError(settings.noise, settings.settingNoise)), MathUtils.rounded(relativeError(photons, settings.settingPhotons)));
    }
    // Show the PoissonGammaGaussian approximation
    double[] approxValues = null;
    if (settings.showApproximation) {
        approxValues = new double[x.length];
        final PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / settings.gain, settings.noise);
        final double expected = photons * settings.gain;
        for (int i = 0; i < approxValues.length; i++) {
            approxValues[i] = fun.likelihood(x[i] - settings.bias, expected);
        }
        yMax = MathUtils.maxDefault(yMax, approxValues);
    }
    // Replot
    pg = pdf(max, photons, settings.gain, settings.noise, (int) settings.bias);
    plot = new Plot(TITLE, "ADU", "Frequency");
    plot.setLimits(limits[0], limits[1], 0, yMax * 1.05);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot.DOT);
    plot.setColor(Color.red);
    plot.addPoints(x, pg, Plot.LINE);
    plot.addLabel(0, 0, label);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, approxValues, Plot.LINE);
    }
    ImageJUtils.display(TITLE, plot);
}