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

Example 1 with PoissonGammaGaussianFisherInformation

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

the class CameraModelFisherInformationAnalysis method createPoissonGammaGaussianFisherInformation.

private static PoissonGammaGaussianFisherInformation createPoissonGammaGaussianFisherInformation(double gain, double sd) {
    if (sd < 0) {
        IJ.error(TITLE, "EM CCD noise must be positive");
        return null;
    }
    if (gain <= 0) {
        IJ.error(TITLE, "EM CCD gain must be positive");
        return null;
    }
    final int sampling = PoissonGammaGaussianFisherInformation.DEFAULT_SAMPLING;
    // sampling = 2;
    final PoissonGammaGaussianFisherInformation fi = new PoissonGammaGaussianFisherInformation(gain, sd, sampling);
    // fi.setRelativeProbabilityThreshold(1e-6);
    // fi.setMinRange(10);
    // fi.setMaxIterations(3);
    // fi.setMaxRange(10);
    // fi.setUse38(false);
    fi.setMeanThreshold(Double.MAX_VALUE);
    return fi;
}

Example 2 with PoissonGammaGaussianFisherInformation

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

the class CameraModelFisherInformationAnalysis method analyse.

private void analyse() {
    final CameraType type1 = CameraType.forNumber(settings.getCamera1Type());
    final CameraType type2 = CameraType.forNumber(settings.getCamera2Type());
    final FiKey key1 = new FiKey(type1, settings.getCamera1Gain(), settings.getCamera1Noise());
    final FiKey key2 = new FiKey(type2, settings.getCamera2Gain(), settings.getCamera2Noise());
    final BasePoissonFisherInformation f1 = createPoissonFisherInformation(key1);
    if (f1 == null) {
        return;
    }
    final BasePoissonFisherInformation f2 = createPoissonFisherInformation(key2);
    if (f2 == null) {
        return;
    }
    final double[] exp = createExponents();
    if (exp == null) {
        return;
    }
    final double[] photons = new double[exp.length];
    for (int i = 0; i < photons.length; i++) {
        photons[i] = Math.pow(10, exp[i]);
    }
    // Get the alpha.
    // This may be from the cache or computed shutdown the executor if it was created.
    double[] alpha1;
    double[] alpha2;
    try {
        alpha1 = getAlpha(photons, exp, f1, key1);
        alpha2 = getAlpha(photons, exp, f2, key2);
    } finally {
        if (es != null) {
            es.shutdownNow();
        }
    }
    // Compute the Poisson Fisher information
    final double[] fi1 = getFisherInformation(alpha1, photons);
    final double[] fi2 = getFisherInformation(alpha2, photons);
    // ==============
    if (debug && f2 instanceof PoissonGammaGaussianFisherInformation) {
        final PoissonGammaGaussianFisherInformation pgg = (PoissonGammaGaussianFisherInformation) f2;
        final double t = 200;
        final double fi = pgg.getFisherInformation(t);
        final double alpha = fi * t;
        final double[][] data1 = pgg.getFisherInformationFunction(false);
        final double[][] data2 = pgg.getFisherInformationFunction(true);
        final double[] fif = data1[1];
        int max = 0;
        for (int j = 1; j < fif.length; j++) {
            if (fif[max] < fif[j]) {
                max = j;
            }
        }
        ImageJUtils.log("PGG(p=%g) max=%g", t, data1[0][max]);
        final String title = TITLE + " photons=" + MathUtils.rounded(t) + " alpha=" + MathUtils.rounded(alpha);
        final Plot plot = new Plot(title, "Count", "FI function");
        double yMax = MathUtils.max(data1[1]);
        yMax = MathUtils.maxDefault(yMax, data2[1]);
        plot.setLimits(data2[0][0], data2[0][data2[0].length - 1], 0, yMax);
        plot.setColor(Color.red);
        plot.addPoints(data1[0], data1[1], Plot.LINE);
        plot.setColor(Color.blue);
        plot.addPoints(data2[0], data2[1], Plot.LINE);
        ImageJUtils.display(title, plot);
    }
    // ==============
    final Color color1 = Color.BLUE;
    final Color color2 = Color.RED;
    final WindowOrganiser wo = new WindowOrganiser();
    // Test if we can use ImageJ support for a X log scale
    final boolean logScaleX = ((float) photons[0] != 0);
    final double[] x = (logScaleX) ? photons : exp;
    final String xTitle = (logScaleX) ? "photons" : "log10(photons)";
    // Get interpolation for alpha. Convert to base e.
    final double[] logU = exp.clone();
    final double scale = Math.log(10);
    for (int i = 0; i < logU.length; i++) {
        logU[i] *= scale;
    }
    final BasePoissonFisherInformation if1 = getInterpolatedPoissonFisherInformation(type1, logU, alpha1, f1);
    final BasePoissonFisherInformation if2 = getInterpolatedPoissonFisherInformation(type2, logU, alpha2, f2);
    // Interpolate with 5 points per sample for smooth curve
    final int n = 5 * exp.length;
    final double[] iexp = new double[n + 1];
    final double[] iphotons = new double[iexp.length];
    final double h = (exp[exp.length - 1] - exp[0]) / n;
    for (int i = 0; i <= n; i++) {
        iexp[i] = exp[0] + i * h;
        iphotons[i] = Math.pow(10, iexp[i]);
    }
    final double[] ix = (logScaleX) ? iphotons : iexp;
    final double[] ialpha1 = getAlpha(if1, iphotons);
    final double[] ialpha2 = getAlpha(if2, iphotons);
    final int pointShape = getPointShape(settings.getPointOption());
    final String name1 = getName(key1);
    final String name2 = getName(key2);
    final String legend = name1 + "\n" + name2;
    String title = "Relative Fisher Information";
    Plot plot = new Plot(title, xTitle, "Noise coefficient (alpha)");
    plot.setLimits(x[0], x[x.length - 1], -0.05, 1.05);
    if (logScaleX) {
        plot.setLogScaleX();
    }
    plot.setColor(color1);
    plot.addPoints(ix, ialpha1, Plot.LINE);
    plot.setColor(color2);
    plot.addPoints(ix, ialpha2, Plot.LINE);
    plot.setColor(Color.BLACK);
    plot.addLegend(legend);
    // Option to show nodes
    if (pointShape != -1) {
        plot.setColor(color1);
        plot.addPoints(x, alpha1, pointShape);
        plot.setColor(color2);
        plot.addPoints(x, alpha2, pointShape);
        plot.setColor(Color.BLACK);
    }
    ImageJUtils.display(title, plot, 0, wo);
    // The approximation should not produce an infinite computation
    double[] limits = new double[] { fi2[fi2.length - 1], fi2[fi2.length - 1] };
    limits = limits(limits, fi1);
    limits = limits(limits, fi2);
    // Check if we can use ImageJ support for a Y log scale
    final boolean logScaleY = ((float) limits[1] <= Float.MAX_VALUE);
    if (!logScaleY) {
        for (int i = 0; i < fi1.length; i++) {
            fi1[i] = Math.log10(fi1[i]);
            fi2[i] = Math.log10(fi2[i]);
        }
        limits[0] = Math.log10(limits[0]);
        limits[1] = Math.log10(limits[1]);
    }
    final String yTitle = (logScaleY) ? "Fisher Information" : "log10(Fisher Information)";
    title = "Fisher Information";
    plot = new Plot(title, xTitle, yTitle);
    plot.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
    if (logScaleX) {
        plot.setLogScaleX();
    }
    if (logScaleY) {
        plot.setLogScaleY();
    }
    plot.setColor(color1);
    plot.addPoints(x, fi1, Plot.LINE);
    plot.setColor(color2);
    plot.addPoints(x, fi2, Plot.LINE);
    plot.setColor(Color.BLACK);
    plot.addLegend(legend);
    // // Option to show nodes
    // This gets messy as the lines are straight
    // if (pointShape != -1)
    // {
    // plot.setColor(color1);
    // plot.addPoints(x, pgFI, pointShape);
    // plot.setColor(color3);
    // plot.addPoints(x, pggFI, pointShape);
    // plot.setColor(Color.BLACK);
    // }
    ImageJUtils.display(title, plot, 0, wo);
    wo.tile();
}