Examples with Plot ij.gui.Plot used on opensource projects

Example 56 with Plot

use of ij.gui.Plot in project GDSC-SMLM by aherbert.

the class TraceLengthAnalysis method draw.

private void draw(WindowOrganiser wo) {
    lastMsdThreshold = settings.msdThreshold;
    lastNormalise = settings.normalise;
    // Find the index in the MSD array
    int index = Arrays.binarySearch(msds, lastMsdThreshold);
    if (index < 0) {
        index = -(index + 1);
    }
    lastIndex = index;
    // Histogram the distributions
    computeHistogram(0, index, lengths, h1);
    computeHistogram(index, lengths.length, lengths, h2);
    final int sum1 = (int) MathUtils.sum(h1);
    final int sum2 = (int) MathUtils.sum(h2);
    final float max1 = createHistogramValues(h1, (lastNormalise) ? sum1 : 1, y1);
    final float max2 = createHistogramValues(h2, (lastNormalise) ? sum2 : 1, y2);
    final String title = "Trace length distribution";
    final Plot plot = new Plot(title, "Length", "Frequency");
    plot.setColor(Color.red);
    plot.addPoints(x1, y1, Plot.BAR);
    plot.setColor(Color.blue);
    plot.addPoints(x1, y2, Plot.BAR);
    plot.setColor(Color.black);
    final double p = 100.0 * sum1 / (sum1 + sum2);
    plot.addLabel(0, 0, String.format("Fixed (red) = %d (%s%%), Moving (blue) = %d (%s%%)", sum1, MathUtils.rounded(p), sum2, MathUtils.rounded(100 - p)));
    plot.setLimits(Double.NaN, Double.NaN, 0, Math.max(max1, max2) * 1.05);
    ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
}

Example 57 with Plot

use of ij.gui.Plot in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    settings = Settings.load();
    // Saved by reference so just save now
    settings.save();
    // Read in multiple traced datasets
    // All datasets must have the same pixel pitch and exposure time
    // Get parameters
    // Convert datasets to tracks
    // For each track compute the 4 local track features using the configured window
    // 
    // Optional:
    // Fit a multi-variate Gaussian mixture model to the data
    // (using the configured number of components/populations)
    // Assign each point in the track using the model.
    // Smooth the assignments.
    // 
    // The alternative is to use the localisation category to assign populations.
    // 
    // Plot histograms of each track parameter, coloured by component
    final List<MemoryPeakResults> combinedResults = new LocalList<>();
    if (!showInputDialog(combinedResults)) {
        return;
    }
    final boolean hasCategory = showHasCategoryDialog(combinedResults);
    if (!showDialog(hasCategory)) {
        return;
    }
    ImageJUtils.log(TITLE + "...");
    final List<Trace> tracks = getTracks(combinedResults, settings.window, settings.minTrackLength);
    if (tracks.isEmpty()) {
        IJ.error(TITLE, "No tracks. Please check the input data and min track length setting.");
        return;
    }
    final Calibration cal = combinedResults.get(0).getCalibration();
    final CalibrationReader cr = new CalibrationReader(cal);
    // Use micrometer / second
    final TypeConverter<DistanceUnit> distanceConverter = cr.getDistanceConverter(DistanceUnit.UM);
    final double exposureTime = cr.getExposureTime() / 1000.0;
    final Pair<int[], double[][]> trackData = extractTrackData(tracks, distanceConverter, exposureTime, hasCategory);
    final double[][] data = trackData.getValue();
    // Histogram the raw data.
    final Array2DRowRealMatrix raw = new Array2DRowRealMatrix(data, false);
    final WindowOrganiser wo = new WindowOrganiser();
    // Store the histogram data for plotting the components
    final double[][] columns = new double[FEATURE_NAMES.length][];
    final double[][] limits = new double[FEATURE_NAMES.length][];
    // Get column data
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        columns[i] = raw.getColumn(i);
        if (i == FEATURE_D) {
            // Plot using a logarithmic scale
            SimpleArrayUtils.apply(columns[i], Math::log10);
        }
        limits[i] = MathUtils.limits(columns[i]);
    }
    // Compute histogram bins
    final int[] bins = new int[FEATURE_NAMES.length];
    if (settings.histogramBins > 0) {
        Arrays.fill(bins, settings.histogramBins);
    } else {
        for (int i = 0; i < FEATURE_NAMES.length; i++) {
            bins[i] = HistogramPlot.getBins(StoredData.create(columns[i]), BinMethod.FD);
        }
        // Use the maximum so all histograms look the same
        Arrays.fill(bins, MathUtils.max(bins));
    }
    // Compute plots
    final Plot[] plots = new Plot[FEATURE_NAMES.length];
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        final double[][] hist = HistogramPlot.calcHistogram(columns[i], limits[i][0], limits[i][1], bins[i]);
        plots[i] = new Plot(TITLE + " " + FEATURE_NAMES[i], getFeatureLabel(i, i == FEATURE_D), "Frequency");
        plots[i].addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(plots[i].getTitle(), plots[i], ImageJUtils.NO_TO_FRONT, wo);
    }
    wo.tile();
    // The component for each data point
    int[] component;
    // The number of components
    int numComponents;
    // Data used to fit the Gaussian mixture model
    double[][] fitData;
    // The fitted model
    MixtureMultivariateGaussianDistribution model;
    if (hasCategory) {
        // Use the category as the component.
        // No fit data and no output model
        fitData = null;
        model = null;
        // The component is stored at the end of the raw track data.
        final int end = data[0].length - 1;
        component = Arrays.stream(data).mapToInt(d -> (int) d[end]).toArray();
        numComponents = MathUtils.max(component) + 1;
        // In the EM algorithm the probability of each data point is computed and normalised to
        // sum to 1. The normalised probabilities are averaged to create the weights.
        // Note the probability of each data point uses the previous weight and the algorithm
        // iterates.
        // This is not a fitted model but the input model so use
        // zero weights to indicate no fitting was performed.
        final double[] weights = new double[numComponents];
        // Remove the trailing component to show the 'model' in a table.
        createModelTable(Arrays.stream(data).map(d -> Arrays.copyOf(d, end)).toArray(double[][]::new), weights, component);
    } else {
        // Multivariate Gaussian mixture EM
        // Provide option to not use the anomalous exponent in the population mix.
        int sortDimension = SORT_DIMENSION;
        if (settings.ignoreAlpha) {
            // Remove index 0. This shifts the sort dimension.
            sortDimension--;
            fitData = Arrays.stream(data).map(d -> Arrays.copyOfRange(d, 1, d.length)).toArray(double[][]::new);
        } else {
            fitData = SimpleArrayUtils.deepCopy(data);
        }
        final MultivariateGaussianMixtureExpectationMaximization mixed = fitGaussianMixture(fitData, sortDimension);
        if (mixed == null) {
            IJ.error(TITLE, "Failed to fit a mixture model");
            return;
        }
        model = sortComponents(mixed.getFittedModel(), sortDimension);
        // For the best model, assign to the most likely population.
        component = assignData(fitData, model);
        // Table of the final model using the original data (i.e. not normalised)
        final double[] weights = model.getWeights();
        numComponents = weights.length;
        createModelTable(data, weights, component);
    }
    // Output coloured histograms of the populations.
    final LUT lut = LutHelper.createLut(settings.lutIndex);
    IntFunction<Color> colourMap;
    if (LutHelper.getColour(lut, 0).equals(Color.BLACK)) {
        colourMap = i -> LutHelper.getNonZeroColour(lut, i, 0, numComponents - 1);
    } else {
        colourMap = i -> LutHelper.getColour(lut, i, 0, numComponents - 1);
    }
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        // Extract the data for each component
        final double[] col = columns[i];
        final Plot plot = plots[i];
        for (int n = 0; n < numComponents; n++) {
            final StoredData feature = new StoredData();
            for (int j = 0; j < component.length; j++) {
                if (component[j] == n) {
                    feature.add(col[j]);
                }
            }
            if (feature.size() == 0) {
                continue;
            }
            final double[][] hist = HistogramPlot.calcHistogram(feature.values(), limits[i][0], limits[i][1], bins[i]);
            // Colour the points
            plot.setColor(colourMap.apply(n));
            plot.addPoints(hist[0], hist[1], Plot.BAR);
        }
        plot.updateImage();
    }
    createTrackDataTable(tracks, trackData, fitData, model, component, cal, colourMap);
// Analysis.
// Assign the original localisations to their track component.
// Q. What about the start/end not covered by the window?
// Save tracks as a dataset labelled with the sub-track ID?
// Output for the bound component and free components track parameters.
// Compute dwell times.
// Other ...
// Track analysis plugin:
// Extract all continuous segments of the same component.
// Produce MSD plot with error bars.
// Fit using FBM model.
}

Example 58 with Plot

use of ij.gui.Plot 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();
}

Example 59 with Plot

use of ij.gui.Plot in project GDSC-SMLM by aherbert.

the class BlinkEstimator method computeBlinkingRate.

/**
 * Compute blinking rate.
 *
 * @param results the results
 * @param verbose the verbose
 * @return the blinking rate
 */
double computeBlinkingRate(MemoryPeakResults results, boolean verbose) {
    parameters = null;
    increaseNFittedPoints = false;
    // Calculate the counts verses dark time curve
    double[] ntd = calculateCounts(results, settings.maxDarkTime, settings.searchDistance, settings.relativeDistance, verbose);
    double[] td = calculateTd(ntd);
    if (verbose) {
        ImageJUtils.log("  Estimate %.0f molecules at td = %.0f ms", ntd[0], td[0]);
    }
    ntd = shift(ntd);
    td = shift(td);
    // Fit curve
    parameters = fit(td, ntd, settings.numberOfFittedPoints, verbose);
    if (parameters == null) {
        return 0;
    }
    // Display
    if (showPlots) {
        final String title = TITLE + " Molecule Counts";
        final Plot plot = new Plot(title, "td (ms)", "Count");
        plot.addPoints(td, ntd, Plot.LINE);
        ImageJUtils.display(title, plot);
        plot.setColor(Color.red);
        plot.addPoints(blinkingModel.getX(), blinkingModel.value(parameters), Plot.CIRCLE);
        // Add the rest that is not fitted
        final double[] xOther = new double[td.length - blinkingModel.size()];
        final double[] yOther = new double[xOther.length];
        for (int i = 0, t = blinkingModel.size(); i < xOther.length; i++, t++) {
            xOther[i] = td[t];
            yOther[i] = blinkingModel.evaluate(td[t], parameters);
        }
        plot.setColor(Color.blue);
        plot.addPoints(xOther, yOther, Plot.CROSS);
        ImageJUtils.display(title, plot);
    }
    // Check if the fitted curve asymptotes above the real curve
    if (blinkingModel.evaluate(td[ntd.length - 1], parameters) < ntd[ntd.length - 1]) {
        if (verbose) {
            ImageJUtils.log("  *** Warning ***");
            ImageJUtils.log("  Fitted curve does not asymptote above real curve. Increase the number" + " of fitted points to sample more of the overcounting regime");
            ImageJUtils.log("  ***************");
        }
        increaseNFittedPoints = true;
    }
    // Blinking rate is 1 + nBlinks
    final double blinkingRate = 1 + parameters[1];
    if (verbose) {
        ImageJUtils.log("  Blinking rate = %s", MathUtils.rounded(blinkingRate, 4));
    }
    return blinkingRate;
}

Example 60 with Plot

use of ij.gui.Plot in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method convolveHistogram.

/**
 * Convolve the histogram. The output is a discrete probability distribution.
 *
 * @param settings the settings
 * @return The histogram
 */
private static double[][] convolveHistogram(CameraModelAnalysisSettings settings) {
    // Find the range of the Poisson
    final PoissonDistribution poisson = new PoissonDistribution(settings.getPhotons());
    final int maxn = poisson.inverseCumulativeProbability(UPPER);
    final double gain = getGain(settings);
    final double noise = getReadNoise(settings);
    final boolean debug = false;
    // Build the Probability Mass/Density Function (PDF) of the distribution:
    // either a Poisson (PMF) or Poisson-Gamma (PDF). The PDF is 0 at all
    // values apart from the step interval.
    // Note: The Poisson-Gamma is computed without the Dirac delta contribution
    // at c=0. This allows correct convolution with the Gaussian of the dirac delta
    // and the rest of the Poisson-Gamma (so matching the simulation).
    final TDoubleArrayList list = new TDoubleArrayList();
    double step;
    String name;
    int upsample = 100;
    // Store the Dirac delta value at c=0. This must be convolved separately.
    double dirac = 0;
    // EM-CCD
    if (settings.getMode() == MODE_EM_CCD) {
        name = "Poisson-Gamma";
        final double m = gain;
        final double p = settings.getPhotons();
        dirac = PoissonGammaFunction.dirac(p);
        // Chose whether to compute a discrete PMF or a PDF using the approximation.
        // Note: The delta function at c=0 is from the PMF of the Poisson. So it is
        // a discrete contribution. This is omitted from the PDF and handled in
        // a separate convolution.
        // noise != 0;
        final boolean discrete = false;
        if (discrete) {
            // Note: This is obsolete as the Poisson-Gamma function is continuous.
            // Sampling it at integer intervals is not valid, especially for low gain.
            // The Poisson-Gamma PDF should be integrated to form a discrete PMF.
            step = 1.0;
            double upper;
            if (settings.getPhotons() < 20) {
                upper = maxn;
            } else {
                // Approximate reasonable range of Poisson as a Gaussian
                upper = settings.getPhotons() + 3 * Math.sqrt(settings.getPhotons());
            }
            final GammaDistribution gamma = new GammaDistribution(null, upper, gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
            final int maxc = (int) gamma.inverseCumulativeProbability(0.999);
            final int minn = Math.max(1, poisson.inverseCumulativeProbability(LOWER));
            // See Ulbrich & Isacoff (2007). Nature Methods 4, 319-321, SI equation 3.
            // Note this is not a convolution of a single Gamma distribution since the shape
            // is modified not the count. So it is a convolution of a distribution made with
            // a gamma of fixed count and variable shape.
            // The count=0 is a special case.
            list.add(PoissonGammaFunction.poissonGammaN(0, p, m));
            final long total = (maxn - minn) * (long) maxc;
            if (total < 1000000) {
                // Full computation
                // G(c) = sum n { (1 / n!) p^n e^-p (1 / ((n-1!)m^n)) c^n-1 e^-c/m }
                // Compute as a log
                // - log(n!) + n*log(p)-p -log((n-1)!) - n * log(m) + (n-1) * log(c) -c/m
                // Note: Both methods work
                LogFactorialCache lfc = LOG_FACTORIAL_CACHE.get().get();
                if (lfc == null) {
                    lfc = new LogFactorialCache();
                    LOG_FACTORIAL_CACHE.set(new SoftReference<>(lfc));
                }
                lfc.ensureRange(minn - 1, maxn);
                final double[] f = new double[maxn + 1];
                final double logm = Math.log(m);
                final double logp = Math.log(p);
                for (int n = minn; n <= maxn; n++) {
                    f[n] = -lfc.getLogFactorialUnsafe(n) + n * logp - p - lfc.getLogFactorialUnsafe(n - 1) - n * logm;
                }
                // double total2 = total;
                for (int c = 1; c <= maxc; c++) {
                    double sum = 0;
                    final double c_m = c / m;
                    final double logc = Math.log(c);
                    for (int n = minn; n <= maxn; n++) {
                        sum += StdMath.exp(f[n] + (n - 1) * logc - c_m);
                    }
                    // sum2 += pd[n] * gd[n].density(c);
                    list.add(sum);
                // total += sum;
                // This should match the approximation
                // double approx = PoissonGammaFunction.poissonGamma(c, p, m);
                // total2 += approx;
                // System.out.printf("c=%d sum=%g approx=%g error=%g\n", c, sum2, approx,
                // uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(sum2, approx));
                }
            // System.out.printf("sum=%g approx=%g error=%g\n", total, total2,
            // uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(total, total2));
            } else {
                // Approximate
                for (int c = 1; c <= maxc; c++) {
                    list.add(PoissonGammaFunction.poissonGammaN(c, p, m));
                }
            }
        } else {
            // This integrates the PDF using the approximation and up-samples together.
            // Compute the sampling interval.
            step = 1.0 / upsample;
            // Reset
            upsample = 1;
            // Compute the integral of [-step/2:step/2] for each point.
            // Use trapezoid integration.
            final double step_2 = step / 2;
            double prev = PoissonGammaFunction.poissonGammaN(0, p, m);
            double next = PoissonGammaFunction.poissonGammaN(step_2, p, m);
            list.add((prev + next) * 0.25);
            double max = 0;
            for (int i = 1; ; i++) {
                // Each remaining point is modelling a PMF for the range [-step/2:step/2]
                prev = next;
                next = PoissonGammaFunction.poissonGammaN(i * step + step_2, p, m);
                final double pp = (prev + next) * 0.5;
                if (max < pp) {
                    max = pp;
                }
                if (pp / max < 1e-5) {
                    // Use this if non-zero since it has been calculated
                    if (pp != 0) {
                        list.add(pp);
                    }
                    break;
                }
                list.add(pp);
            }
        }
        // Ensure the combined sum of PDF and Dirac is 1
        final double expected = 1 - dirac;
        // Require an odd number to get an even number (n) of sub-intervals:
        if (list.size() % 2 == 0) {
            list.add(0);
        }
        final double[] g = list.toArray();
        // Number of sub intervals
        final int n = g.length - 1;
        // h = (a-b) / n = sub-interval width
        final double h = 1;
        double sum2 = 0;
        double sum4 = 0;
        for (int j = 1; j <= n / 2 - 1; j++) {
            sum2 += g[2 * j];
        }
        for (int j = 1; j <= n / 2; j++) {
            sum4 += g[2 * j - 1];
        }
        final double sum = (h / 3) * (g[0] + 2 * sum2 + 4 * sum4 + g[n]);
        // Check
        // System.out.printf("Sum=%g Expected=%g\n", sum * step, expected);
        SimpleArrayUtils.multiply(g, expected / sum);
        list.resetQuick();
        list.add(g);
    } else {
        name = "Poisson";
        // Apply fixed gain. Just change the step interval of the PMF.
        step = gain;
        for (int n = 0; n <= maxn; n++) {
            list.add(poisson.probability(n));
        }
        final double p = poisson.probability(list.size());
        if (p != 0) {
            list.add(p);
        }
    }
    // Debug
    if (debug) {
        final String title = name;
        final Plot plot = new Plot(title, "x", "y");
        plot.addPoints(SimpleArrayUtils.newArray(list.size(), 0, step), list.toArray(), Plot.LINE);
        ImageJUtils.display(title, plot);
    }
    double zero = 0;
    double[] pg = list.toArray();
    // Sample Gaussian
    if (noise > 0) {
        step /= upsample;
        pg = list.toArray();
        // Convolve with Gaussian kernel
        final double[] kernel = GaussianKernel.makeGaussianKernel(Math.abs(noise) / step, 6, true);
        if (upsample != 1) {
            // Use scaled convolution. This is faster that zero filling distribution g.
            pg = Convolution.convolve(kernel, pg, upsample);
        } else if (dirac > 0.01) {
            // The Poisson-Gamma may be stepped at low mean causing wrap artifacts in the FFT.
            // This is a problem if most of the probability is in the Dirac.
            // Otherwise it can be ignored and the FFT version is OK.
            pg = Convolution.convolve(kernel, pg);
        } else {
            pg = Convolution.convolveFast(kernel, pg);
        }
        // The convolution will have created a larger array so we must adjust the offset for this
        final int radius = kernel.length / 2;
        zero -= radius * step;
        // Add convolution of the dirac delta function.
        if (dirac != 0) {
            // We only need to convolve the Gaussian at c=0
            for (int i = 0; i < kernel.length; i++) {
                pg[i] += kernel[i] * dirac;
            }
        }
        // Debug
        if (debug) {
            String title = "Gaussian";
            Plot plot = new Plot(title, "x", "y");
            plot.addPoints(SimpleArrayUtils.newArray(kernel.length, -radius * step, step), kernel, Plot.LINE);
            ImageJUtils.display(title, plot);
            title = name + "-Gaussian";
            plot = new Plot(title, "x", "y");
            plot.addPoints(SimpleArrayUtils.newArray(pg.length, zero, step), pg, Plot.LINE);
            ImageJUtils.display(title, plot);
        }
        zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1.0);
        pg = list.toArray();
        zero = (int) Math.floor(zero);
        step = 1.0;
    // No convolution means we have the Poisson PMF/Poisson-Gamma PDF already
    } else if (step != 1) {
        // Sample to 1.0 pixel step interval.
        if (settings.getMode() == MODE_EM_CCD) {
            // Poisson-Gamma PDF
            zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1 - dirac);
            pg = list.toArray();
            zero = (int) Math.floor(zero);
            // Add the dirac delta function.
            if (dirac != 0) {
                // Note: zero is the start of the x-axis. This value should be -1.
                assert (int) zero == -1;
                // Use as an offset to find the actual zero.
                pg[-(int) zero] += dirac;
            }
        } else {
            // Poisson PMF
            // Simple non-interpolated expansion.
            // This should be used when there is no Gaussian convolution.
            final double[] pd = pg;
            list.resetQuick();
            // Account for rounding.
            final Round round = getRound(settings);
            final int maxc = round.round(pd.length * step + 1);
            pg = new double[maxc];
            for (int n = pd.length; n-- > 0; ) {
                pg[round.round(n * step)] += pd[n];
            }
            if (pg[0] != 0) {
                list.add(0);
                list.add(pg);
                pg = list.toArray();
                zero--;
            }
        }
        step = 1.0;
    } else {
        // Add the dirac delta function.
        list.setQuick(0, list.getQuick(0) + dirac);
    }
    return new double[][] { SimpleArrayUtils.newArray(pg.length, zero, step), pg };
}