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

Example 1 with PoissonFunction

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

the class EmGainAnalysis method plotPmf.

@SuppressWarnings("unused")
private void plotPmf() {
    if (!showPmfDialog()) {
        return;
    }
    final double step = getStepSize(settings.settingPhotons, settings.settingGain, settings.settingNoise);
    final Pdf pdf = pdf(0, step, settings.settingPhotons, settings.settingGain, settings.settingNoise);
    double[] pmf = pdf.probability;
    double yMax = MathUtils.max(pmf);
    // Get the approximation
    LikelihoodFunction fun;
    switch(settings.approximationType) {
        case 3:
            fun = new PoissonFunction(1.0 / settings.settingGain);
            break;
        case 2:
            // Use adaptive normalisation
            fun = PoissonGaussianFunction2.createWithStandardDeviation(1.0 / settings.settingGain, settings.settingNoise);
            break;
        case 1:
            // Create Poisson-Gamma (no Gaussian noise)
            fun = createPoissonGammaGaussianFunction(0);
            break;
        case 0:
        default:
            fun = createPoissonGammaGaussianFunction(settings.settingNoise);
    }
    double expected = settings.settingPhotons;
    if (settings.settingOffset != 0) {
        expected += settings.settingOffset * expected / 100.0;
    }
    // Normalise
    final boolean normalise = false;
    if (normalise) {
        final double sum = MathUtils.sum(pmf);
        for (int i = pmf.length; i-- > 0; ) {
            pmf[i] /= sum;
        }
    }
    // Get CDF
    double sum = 0;
    double sum2 = 0;
    double[] x = pdf.x;
    double[] fvalues = new double[x.length];
    double[] cdf1 = new double[pmf.length];
    double[] cdf2 = new double[pmf.length];
    for (int i = 0; i < cdf1.length; i++) {
        sum += pmf[i] * step;
        cdf1[i] = sum;
        fvalues[i] = fun.likelihood(x[i], expected);
        sum2 += fvalues[i] * step;
        cdf2[i] = sum2;
    }
    // Truncate x for plotting
    int max = 0;
    double plimit = 1 - settings.tail;
    while (sum < plimit && max < pmf.length) {
        sum += pmf[max] * step;
        if (sum > 0.5 && pmf[max] == 0) {
            break;
        }
        max++;
    }
    int min = pmf.length;
    sum = 0;
    plimit = 1 - settings.head;
    while (sum < plimit && min > 0) {
        min--;
        sum += pmf[min] * step;
        if (sum > 0.5 && pmf[min] == 0) {
            break;
        }
    }
    pmf = Arrays.copyOfRange(pmf, min, max);
    x = Arrays.copyOfRange(x, min, max);
    fvalues = Arrays.copyOfRange(fvalues, min, max);
    if (settings.showApproximation) {
        yMax = MathUtils.maxDefault(yMax, fvalues);
    }
    final String label = String.format("Gain=%s, noise=%s, photons=%s", MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    final Plot plot = new Plot("PMF", "ADUs", "p");
    plot.setLimits(x[0], x[x.length - 1], 0, yMax);
    plot.setColor(Color.red);
    plot.addPoints(x, pmf, Plot.LINE);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, fvalues, Plot.LINE);
    }
    plot.setColor(Color.magenta);
    plot.drawLine(settings.settingPhotons * settings.settingGain, 0, settings.settingPhotons * settings.settingGain, yMax);
    plot.setColor(Color.black);
    plot.addLabel(0, 0, label);
    final PlotWindow win1 = ImageJUtils.display("PMF", plot);
    // Plot the difference between the actual and approximation
    final double[] delta = new double[fvalues.length];
    for (int i = 0; i < fvalues.length; i++) {
        if (pmf[i] == 0 && fvalues[i] == 0) {
            continue;
        }
        if (settings.relativeDelta) {
            delta[i] = DoubleEquality.relativeError(fvalues[i], pmf[i]) * Math.signum(fvalues[i] - pmf[i]);
        } else {
            delta[i] = fvalues[i] - pmf[i];
        }
    }
    final Plot plot2 = new Plot("PMF delta", "ADUs", (settings.relativeDelta) ? "Relative delta" : "delta");
    final double[] limits = MathUtils.limits(delta);
    plot2.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
    plot2.setColor(Color.red);
    plot2.addPoints(x, delta, Plot.LINE);
    plot2.setColor(Color.magenta);
    plot2.drawLine(settings.settingPhotons * settings.settingGain, limits[0], settings.settingPhotons * settings.settingGain, limits[1]);
    plot2.setColor(Color.black);
    plot2.addLabel(0, 0, label + ((settings.settingOffset == 0) ? "" : ", expected = " + MathUtils.rounded(expected / settings.settingGain)));
    final WindowOrganiser wo = new WindowOrganiser();
    final PlotWindow win2 = ImageJUtils.display("PMF delta", plot2, wo);
    if (wo.isNotEmpty()) {
        final Point p2 = win1.getLocation();
        p2.y += win1.getHeight();
        win2.setLocation(p2);
    }
    // Plot the CDF of each distribution.
    // Compute the Kolmogorov distance as the supremum (maximum)
    // difference between the two cumulative probability distributions.
    // https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
    double kolmogorovDistance = 0;
    double xd = x[0];
    for (int i = 0; i < cdf1.length; i++) {
        final double dist = Math.abs(cdf1[i] - cdf2[i]);
        if (kolmogorovDistance < dist) {
            kolmogorovDistance = dist;
            xd = pdf.x[i];
        }
    }
    cdf1 = Arrays.copyOfRange(cdf1, min, max);
    cdf2 = Arrays.copyOfRange(cdf2, min, max);
    final Plot plot3 = new Plot("CDF", "ADUs", "p");
    yMax = 1.05;
    plot3.setLimits(x[0], x[x.length - 1], 0, yMax);
    plot3.setColor(Color.red);
    plot3.addPoints(x, cdf1, Plot.LINE);
    plot3.setColor(Color.blue);
    plot3.addPoints(x, cdf2, Plot.LINE);
    plot3.setColor(Color.magenta);
    plot3.drawLine(settings.settingPhotons * settings.settingGain, 0, settings.settingPhotons * settings.settingGain, yMax);
    plot3.drawDottedLine(xd, 0, xd, yMax, 2);
    plot3.setColor(Color.black);
    plot3.addLabel(0, 0, label + ", Kolmogorov distance = " + MathUtils.rounded(kolmogorovDistance) + " @ " + xd);
    plot3.addLegend("CDF\nApprox");
    final int size = wo.size();
    final PlotWindow win3 = ImageJUtils.display("CDF", plot3, wo);
    if (size != wo.size()) {
        final Point p2 = win1.getLocation();
        p2.x += win1.getWidth();
        win3.setLocation(p2);
    }
}

Example 2 with PoissonFunction

use of uk.ac.sussex.gdsc.smlm.function.PoissonFunction 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();
    }
}