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

Example 86 with Plot

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

the class PsfCreator method runUsingAlignment.

private void runUsingAlignment() {
    if (!showAlignmentDialog()) {
        return;
    }
    boxRadius = (int) Math.ceil(settings.getRadius());
    final CalibrationReader calibration = new CalibrationReader(settings.getCalibration());
    // Limit this
    final int halfBoxRadius = boxRadius / 2;
    settings.setAnalysisWindow(Math.min(settings.getAnalysisWindow(), halfBoxRadius));
    // Find the selected PSF spots x,y,z centre
    // We offset the centre to the middle of pixel.
    BasePoint[] centres = getSpots(0.5f, false);
    if (centres.length == 0) {
        IJ.error(TITLE, "No PSFs");
        return;
    }
    CameraModel cameraModel = null;
    if (calibration.isScmos()) {
        cameraModel = CameraModelManager.load(calibration.getCameraModelName());
        if (cameraModel == null) {
            IJ.error(TITLE, "No camera model");
            return;
        }
        cameraModel = PeakFit.cropCameraModel(cameraModel, IJImageSource.getBounds(imp), null, true);
    } else {
        cameraModel = new CcdCameraModel(calibration.getBias(), 1);
    }
    // Extract the image data for processing as float
    final float[][] image = CreateData.extractImageStack(imp, 0, imp.getStackSize() - 1);
    for (final float[] data : image) {
        cameraModel.removeBiasAndGain(data);
    }
    zSelector = new PsfCentreSelector();
    // Relocate the initial centres
    ImageJUtils.showStatus("Relocating initial centres");
    centres = relocateCentres(image, centres);
    if (centres == null) {
        return;
    }
    zRadius = (int) Math.ceil(settings.getAlignmentZRadius());
    // Check the region overlap in 3D and exclude overlapping PSFs
    boolean[] bad = findSpotOverlap(centres, null);
    centres = getNonBadSpots(centres, bad);
    if (centres.length == 0) {
        IJ.error(TITLE, "No PSFs without neighbours within the box region");
        return;
    }
    // Multi-thread for speed
    if (threadPool == null) {
        threadPool = Executors.newFixedThreadPool(Prefs.getThreads());
    }
    // Extract each PSF into a scaled PSF
    ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", 0));
    ExtractedPsf[] psfs = extractPsfs(image, centres);
    Point location = null;
    // Iterate until centres have converged
    boolean converged = false;
    for (int iter = 0; !converged && iter < settings.getMaxIterations(); iter++) {
        if (ImageJUtils.isInterrupted()) {
            return;
        }
        // Combine all PSFs
        ImageJUtils.showStatus(String.format("[%d] Aligning PSFs", iter + 1));
        final ExtractedPsf combined = combine(psfs);
        combined.createProjections();
        // Get the current combined z-centre.
        // This is used to get the centre of mass for repositioning.
        // It also effects the alignment so do it for the first iteration.
        zSelector.setPsf(combined);
        if (iter == 0) {
            // TODO - check if the z-centre should be guessed here.
            // We assume that the combined PSF may be easier to guess if the initial
            // guess for each individual PSF was OK. It may not be necessary since all
            // the PSFs are combined around their z-centres. Once alignment has
            // started we skip this step.
            zSelector.analyse();
            zSelector.guessZCentre();
        }
        if (settings.getInteractiveMode()) {
            if (iter != 0) {
                zSelector.analyse();
            }
            // zSelector.guessZCentre();
            final double dz = zSelector.run("Update combined PSF z-centre", true, false, false, null);
            if (dz < 0) {
                return;
            }
        }
        // Align each to the combined PSF
        final float[][] translation = align(combined, psfs);
        if (ImageJUtils.isInterrupted()) {
            return;
        }
        // Find the new centre using the old centre plus the alignment shift
        for (int j = 0; j < psfs.length; j++) {
            centres[j] = psfs[j].updateCentre(translation[j]);
            // Update to get the correct scale
            translation[j][0] = centres[j].getX() - psfs[j].centre.getX();
            translation[j][1] = centres[j].getY() - psfs[j].centre.getY();
            translation[j][2] = centres[j].getZ() - psfs[j].centre.getZ();
            ImageJUtils.log("[%d] Centre %d : Shift X = %s : Shift Y = %s : Shift Z = %s", iter, j + 1, rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
        }
        final boolean[] excluded = new boolean[psfs.length];
        if (checkAlignments) {
            combined.show(TITLE_PSF);
            // Ask about each centre in turn.
            // Update Point ROI using float coordinates and set image slice to
            // correct z-centre.
            // imp.saveRoi();
            imp.killRoi();
            final ImageCanvas ic = imp.getCanvas();
            // ic.setMagnification(16);
            int reject = 0;
            final float box = boxRadius + 0.5f;
            final int n = imp.getStackSize();
            for (int j = 0; j < centres.length; j++) {
                psfs[j].show(TITLE_SPOT_PSF);
                final Overlay o = new Overlay();
                o.add(createRoi(psfs[j].centre.getX(), psfs[j].centre.getY(), Color.RED));
                final float cx = centres[j].getX();
                final float cy = centres[j].getY();
                o.add(createRoi(cx, cy, Color.GREEN));
                final Roi roi = new Roi(cx - box, cy - box, 2 * box, 2 * box);
                o.add(roi);
                // The centre is absolute within the original stack
                imp.setSlice(MathUtils.clip(1, n, Math.round(centres[j].getZ())));
                final Rectangle r = ic.getSrcRect();
                final int x = centres[j].getXint();
                final int y = centres[j].getYint();
                if (!r.contains(x, y)) {
                    r.x = x - r.width / 2;
                    r.y = y - r.height / 2;
                    ic.setSourceRect(r);
                }
                imp.setOverlay(o);
                imp.updateAndDraw();
                final NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
                ImageJUtils.addMessage(gd, "Shift X = %s\nShift Y = %s\nShift Z = %s", rounder.toString(translation[j][0]), rounder.toString(translation[j][1]), rounder.toString(translation[j][2]));
                final int spotIndex = j;
                gd.addAndGetButton("Exclude spot", event -> {
                    if (excluded[spotIndex]) {
                        ImageJUtils.log("Included spot %d", spotIndex + 1);
                        excluded[spotIndex] = false;
                    } else {
                        ImageJUtils.log("Excluded spot %d", spotIndex + 1);
                        excluded[spotIndex] = true;
                    }
                });
                gd.enableYesNoCancel("Accept", "Reject");
                if (location != null) {
                    gd.setLocation(location.x, location.y);
                }
                gd.showDialog();
                if (gd.wasCanceled()) {
                    resetImp();
                    return;
                }
                final boolean failed = excluded[spotIndex] || !gd.wasOKed();
                if (failed) {
                    reject++;
                    centres[j] = psfs[j].centre;
                    // For RMSD computation
                    Arrays.fill(translation[j], 0f);
                }
                location = gd.getLocation();
            }
            resetImp();
            if (reject == psfs.length) {
                IJ.error(TITLE, "No PSF translations were accepted");
                return;
            }
        }
        bad = findSpotOverlap(centres, excluded);
        final int badCount = count(bad);
        final int excludedCount = count(excluded);
        int ok = bad.length - badCount - excludedCount;
        if (ok < bad.length) {
            if (badCount != 0 && settings.getInteractiveMode()) {
                final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
                gd.addMessage("Warning: Regions now overlap!");
                gd.addMessage("OK = " + TextUtils.pleural(ok, "PSF"));
                gd.addMessage("Overlapping = " + TextUtils.pleural(badCount, "PSF"));
                // gd.addMessage("Excluded = " + TextUtils.pleural(excludedCount, "PSF"));
                gd.enableYesNoCancel("Exclude", "Include");
                if (location != null) {
                    gd.setLocation(location.x, location.y);
                }
                gd.showDialog();
                if (gd.wasCanceled()) {
                    resetImp();
                    return;
                }
                if (!gd.wasOKed()) {
                    // allow bad spots
                    Arrays.fill(bad, false);
                    ok = bad.length;
                }
                location = gd.getLocation();
            }
            if (ok == 0) {
                IJ.error(TITLE, "No PSFs remaining");
                resetImp();
                return;
            }
        }
        // Merge bad and excluded to get new centres
        for (int i = 0; i < bad.length; i++) {
            if (excluded[i]) {
                bad[i] = true;
            }
        }
        ok = bad.length - count(bad);
        final BasePoint[] newCentres = getNonBadSpots(centres, bad);
        // Find the change in centres
        final double[] rmsd = new double[2];
        for (int j = 0; j < psfs.length; j++) {
            if (bad[j]) {
                continue;
            }
            rmsd[0] += MathUtils.pow2(translation[j][0]) + MathUtils.pow2(translation[j][1]);
            rmsd[1] += MathUtils.pow2(translation[j][2]);
        }
        for (int j = 0; j < 2; j++) {
            rmsd[j] = Math.sqrt(rmsd[j] / ok);
        }
        ImageJUtils.showStatus(String.format("[%d] Checking combined PSF", iter + 1));
        // Compute CoM shift using the current z-centre and z-window
        final double[] shift = combined.getCentreOfMassXyShift(zSelector.getCentreSlice());
        final double shiftd = Math.sqrt(shift[0] * shift[0] + shift[1] * shift[1]);
        ImageJUtils.log("[%d] RMSD XY = %s : RMSD Z = %s : Combined CoM shift = %s,%s (%s)", iter, rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
        if (settings.getInteractiveMode()) {
            // Ask if OK to continue?
            final GenericDialog gd = new GenericDialog(TITLE);
            ImageJUtils.addMessage(gd, "RMSD XY = %s\nRMSD Z = %s\nCombined CoM shift = %s,%s (%s)", rounder.toString(rmsd[0]), rounder.toString(rmsd[1]), rounder.toString(shift[0]), rounder.toString(shift[1]), rounder.toString(shiftd));
            // Check if we can do more iterations
            if (iter + 1 < settings.getMaxIterations()) {
                gd.enableYesNoCancel("Continue", "Converged");
            } else {
                gd.setOKLabel("Converged");
            }
            gd.showDialog();
            if (gd.wasCanceled()) {
                return;
            }
            converged = !gd.wasOKed();
        } else {
            // Check convergence thresholds
            converged = rmsd[0] < settings.getRmsdXyThreshold() && rmsd[1] < settings.getRmsdZThreshold() && shiftd < settings.getComShiftThreshold();
        }
        // For the next round we move to the non-overlapping spots
        centres = newCentres;
        // Update the centres using the centre-of-mass of the combined PSF
        centres = updateUsingCentreOfMassXyShift(shift, shiftd, combined, centres);
        // Extract each PSF into a scaled PSF
        ImageJUtils.showStatus(String.format("[%d] Extracting PSFs", iter + 1));
        psfs = extractPsfs(image, centres);
    }
    // Combine all
    ExtractedPsf combined = combine(psfs);
    // Show an interactive dialog for cropping the PSF and choosing the
    // final output
    final PsfOutputSelector cropSelector = new PsfOutputSelector(combined);
    combined = cropSelector.run();
    if (combined == null) {
        return;
    }
    if (settings.getUpdateRoi()) {
        final float[] ox = new float[centres.length];
        final float[] oy = new float[centres.length];
        for (int i = 0; i < centres.length; i++) {
            ox[i] = centres[i].getX();
            oy[i] = centres[i].getY();
        }
        imp.setRoi(new OffsetPointRoi(ox, oy));
    }
    // For an image PSF we can just enlarge the PSF and window.
    // For a CSpline then we already have the 3D cubic spline function.
    // However we want to post-process the function to allow windowing and
    // normalisation. So we enlarge by 3 in each dimension.
    // The CSpline can be created by solving the coefficients for the
    // 4x4x4 (64) sampled points on each node.
    int magnification;
    if (settings.getOutputType() == OUTPUT_TYPE_IMAGE_PSF) {
        magnification = settings.getPsfMagnification();
    } else {
        magnification = 3;
    }
    // Enlarge the combined PSF for final processing
    ExtractedPsf finalPsf = combined.enlarge(magnification, threadPool);
    // Show a dialog to collect final z-centre interactively
    ImageJUtils.showStatus("Analysing PSF");
    zSelector.setPsf(finalPsf);
    zSelector.analyse();
    // zSelector.guessZCentre(); // No need to guess the centre
    final double dz = zSelector.run("Finalise PSF", true, true, true, null);
    if (dz < 0) {
        return;
    }
    zCentre = zSelector.getCentreSlice();
    if (settings.getCropToZCentre()) {
        finalPsf = finalPsf.cropToZCentre(zCentre);
        // Back to 1-based index
        zCentre = finalPsf.stackZCentre + 1;
    }
    // When click ok the background is subtracted from the PSF
    // All pixels below the background are set to zero
    // Apply a Tukey window to roll-off to zero at the outer pixels
    ImageJUtils.showStatus("Windowing PSF");
    final double[] wx = ImageWindow.tukeyEdge(finalPsf.maxx, settings.getWindow());
    final double[] wz = ImageWindow.tukeyEdge(finalPsf.psf.length, settings.getWindow());
    // Normalisation so the max intensity frame is one
    final float[][] psf = finalPsf.psf;
    final int maxz = psf.length;
    final double[] sum = new double[maxz];
    for (int z = 0; z < maxz; z++) {
        sum[z] = applyWindow(psf[z], z, wx, wz, zSelector.background);
    }
    // Smooth the intensity
    ImageJUtils.showStatus("Normalising PSF");
    final Smoother smoother = zSelector.ssmoother;
    final double[] ssum = smoother.smooth(sum).getDSmooth();
    // Compute normalisation and apply.
    SimpleArrayUtils.multiply(ssum, 1.0 / MathUtils.max(ssum));
    for (int z = 0; z < psf.length; z++) {
        if (sum[z] != 0) {
            SimpleArrayUtils.multiply(psf[z], ssum[z] / sum[z]);
        }
        sum[z] = MathUtils.sum(psf[z]);
    }
    // Show the final intensity profile
    final double[] slice = SimpleArrayUtils.newArray(maxz, 1, 1.0);
    final Plot plot = new Plot(TITLE_SIGNAL, "Slice", "Signal");
    final double[] range = MathUtils.limits(sum);
    plot.setLimits(1, maxz, range[0], range[1]);
    plot.setColor(Color.black);
    plot.addPoints(slice, sum, Plot.LINE);
    ImageJUtils.display(TITLE_SIGNAL, plot);
    // Create a new extracted PSF and show
    ImageJUtils.showStatus("Displaying PSF");
    magnification = finalPsf.magnification;
    finalPsf = new ExtractedPsf(psf, finalPsf.maxx, finalPsf.centre, magnification);
    finalPsf.createProjections();
    psfOut = finalPsf.show(TITLE_PSF, zCentre);
    psfImp = psfOut[0];
    // Add image info
    final int imageCount = centres.length;
    final ImagePSF.Builder imagePsf = ImagePsfHelper.create(zCentre, nmPerPixel / magnification, settings.getNmPerSlice() / magnification, imageCount, 0, createNote()).toBuilder();
    // Add the CoM
    // Find the XY centre around the z centre
    final double[] com = getCentreOfMassXy(finalPsf.psf, finalPsf.maxx, finalPsf.maxy, zCentre - 1, settings.getComWindow(), getComXyBorder(finalPsf.maxx, finalPsf.maxy));
    imagePsf.setXCentre(com[0]);
    imagePsf.setYCentre(com[1]);
    imagePsf.setZCentre(zCentre - 1);
    psfImp.setProperty("Info", ImagePsfHelper.toString(imagePsf));
    psfImp.setRoi(new OffsetPointRoi(com[0], com[1]));
    psfImp.setSlice(zCentre);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    ImageJUtils.log("Final Centre-of-mass = %s,%s\n", rounder.toString(com[0]), rounder.toString(com[1]));
    ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images\n", psfImp.getTitle(), zCentre, MathUtils.rounded(imagePsf.getPixelSize(), 3), MathUtils.rounded(imagePsf.getPixelDepth(), 3), imageCount);
    if (settings.getOutputType() == OUTPUT_TYPE_CSPLINE) {
        // Ask this again as it is important
        // if (TextUtils.isNullOrEmpty(settings.getSplineFilename()))
        // {
        final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
        gd.addFilenameField("Spline_filename", settings.getSplineFilename());
        gd.showDialog(true);
        if (gd.wasCanceled()) {
            return;
        }
        settings.setSplineFilename(gd.getNextString());
        // }
        if (!TextUtils.isNullOrEmpty(settings.getSplineFilename())) {
            // Save the result ...
            IJ.showStatus("Creating cubic spline");
            final CubicSplinePsf cubicSplinePsf = CubicSplineManager.createCubicSpline(imagePsf, psfImp.getImageStack(), settings.getSinglePrecision());
            IJ.showStatus("Saving cubic spline");
            CubicSplineManager.save(cubicSplinePsf, settings.getSplineFilename());
            final String msg = "Spline saved to " + settings.getSplineFilename();
            IJ.showStatus(msg);
            IJ.log(msg);
            // To leave the status message
            return;
        }
    }
    IJ.showStatus("");
}

Example 87 with Plot

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

the class PsfCreator method plotSignalAtSpecifiedSd.

/**
 * Show a plot of the amount of signal within N x SD for each z position. This indicates how much
 * the PSF has spread from the original Gaussian shape.
 *
 * @param psf The PSF
 * @param fittedSd The width of the PSF (in pixels)
 * @param factor The factor to use
 * @param slice The slice used to create the label
 */
private void plotSignalAtSpecifiedSd(ImageStack psf, double fittedSd, double factor, int slice) {
    if (signalZ == null) {
        // Get the bounds
        final int radius = (int) Math.round(fittedSd * factor);
        final int min = Math.max(0, psf.getWidth() / 2 - radius);
        final int max = Math.min(psf.getWidth() - 1, psf.getWidth() / 2 + radius);
        // Create a circle mask of the PSF projection
        final ByteProcessor circle = new ByteProcessor(max - min + 1, max - min + 1);
        circle.setColor(255);
        circle.fillOval(0, 0, circle.getWidth(), circle.getHeight());
        final byte[] mask = (byte[]) circle.getPixels();
        // Sum the pixels within the mask for each slice
        signalZ = new double[psf.getSize()];
        signal = new double[psf.getSize()];
        for (int i = 0; i < psf.getSize(); i++) {
            double sum = 0;
            final float[] data = (float[]) psf.getProcessor(i + 1).getPixels();
            for (int y = min, ii = 0; y <= max; y++) {
                int index = y * psf.getWidth() + min;
                for (int x = min; x <= max; x++, ii++, index++) {
                    if (mask[ii] != 0 && data[index] > 0) {
                        sum += data[index];
                    }
                }
            }
            double total = 0;
            for (final float f : data) {
                if (f > 0) {
                    total += f;
                }
            }
            signalZ[i] = i + 1;
            signal[i] = 100 * sum / total;
        }
        signalTitle = String.format("%% PSF signal at %s x SD", MathUtils.rounded(factor, 3));
        signalLimits = MathUtils.limits(signal);
    }
    // Plot the sum
    final boolean alignWindows = (WindowManager.getFrame(signalTitle) == null);
    final double total = signal[slice - 1];
    final Plot plot = new Plot(signalTitle, "z", "Signal");
    plot.addPoints(signalZ, signal, Plot.LINE);
    plot.addLabel(0, 0, String.format("Total = %s. z = %s nm", MathUtils.rounded(total), MathUtils.rounded((slice - zCentre) * settings.getNmPerSlice())));
    plot.setColor(Color.green);
    plot.drawLine(slice, signalLimits[0], slice, signalLimits[1]);
    plot.setColor(Color.blue);
    final PlotWindow plotWindow = ImageJUtils.display(signalTitle, plot);
    if (alignWindows && plotWindow != null) {
        final PlotWindow otherWindow = getPlot(TITLE_AMPLITUDE);
        if (otherWindow != null) {
            // Put the two plots tiled together so both are visible
            final Point l = plotWindow.getLocation();
            l.x = otherWindow.getLocation().x + otherWindow.getWidth();
            l.y = otherWindow.getLocation().y;
            plotWindow.setLocation(l);
        }
    }
}

Example 88 with Plot

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

the class Noise method drawPlot.

/**
 * Build a plot of the noise estimate from the current frame. Limit the preview to 100 frames.
 */
private void drawPlot() {
    final NoiseEstimatorMethod[] values = SettingsManager.getNoiseEstimatorMethodValues();
    final NoiseEstimator.Method method1 = FitProtosHelper.convertNoiseEstimatorMethod(values[settings.algorithm]);
    final NoiseEstimator.Method method2 = FitProtosHelper.convertNoiseEstimatorMethod(values[settings.algorithm2]);
    IJ.showStatus("Estimating noise ...");
    final boolean twoMethods = method1 != method2;
    final boolean preserveResiduals = method1.name().contains("Residuals") && method2.name().contains("Residuals") && twoMethods;
    final int current = imp.getCurrentSlice();
    final int stackSize = imp.getStackSize();
    final int preview = 100;
    int start = current;
    int end = current + preview;
    if (end > stackSize) {
        final int shift = end - stackSize;
        start -= shift;
        end = stackSize;
        start = Math.max(1, start);
    }
    final int size = end - start + 1;
    final double[] xValues = new double[size];
    final double[] yValues1 = new double[size];
    final double[] yValues2 = (twoMethods) ? new double[size] : null;
    final ImageStack stack = imp.getImageStack();
    final Rectangle bounds = imp.getProcessor().getRoi();
    float[] buffer = null;
    for (int slice = start, i = 0; slice <= end; slice++, i++) {
        IJ.showProgress(i, size);
        final ImageProcessor ip = stack.getProcessor(slice);
        buffer = ImageJImageConverter.getData(ip.getPixels(), ip.getWidth(), ip.getHeight(), bounds, buffer);
        cameraModel.removeBiasAndGain(bounds, buffer);
        final NoiseEstimator ne = NoiseEstimator.wrap(buffer, bounds.width, bounds.height);
        ne.setPreserveResiduals(preserveResiduals);
        ne.setRange(settings.lowestPixelsRange);
        xValues[i] = slice;
        yValues1[i] = ne.getNoise(method1);
        if (yValues2 != null) {
            yValues2[i] = ne.getNoise(method2);
        }
    }
    IJ.showProgress(1);
    IJ.showStatus("Plotting noise ...");
    // Get limits
    final double[] a = Tools.getMinMax(xValues);
    final double[] b1 = Tools.getMinMax(yValues1);
    if (twoMethods) {
        final double[] b2 = Tools.getMinMax(yValues2);
        b1[0] = Math.min(b1[0], b2[0]);
        b1[1] = Math.max(b1[1], b2[1]);
    }
    final String title = imp.getTitle() + " Noise";
    final Plot plot = new Plot(title, "Slice", yAxisTitle);
    double range = b1[1] - b1[0];
    if (range == 0) {
        range = 1;
    }
    plot.setLimits(a[0], a[1], b1[0] - 0.05 * range, b1[1] + 0.05 * range);
    plot.setColor(Color.blue);
    plot.addPoints(xValues, yValues1, Plot.LINE);
    String label = String.format("%s (Blue) = %s", trim(method1.getName()), MathUtils.rounded(Statistics.create(yValues1).getMean()));
    if (twoMethods) {
        plot.setColor(Color.red);
        plot.addPoints(xValues, yValues2, Plot.LINE);
        label += String.format(", %s (Red) = %s", trim(method2.getName()), MathUtils.rounded(Statistics.create(yValues2).getMean()));
    }
    plot.setColor(Color.black);
    plot.addLabel(0, 0, label);
    ImageJUtils.display(title, plot);
    IJ.showStatus("");
}

Example 89 with Plot

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

the class MeanVarianceTest method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    settings = Settings.load();
    settings.save();
    String helpKey = "mean-variance-test";
    if (ImageJUtils.isExtraOptions()) {
        final ImagePlus imp = WindowManager.getCurrentImage();
        if (imp.getStackSize() > 1) {
            final GenericDialog gd = new GenericDialog(TITLE);
            gd.addMessage("Perform single image analysis on the current image?");
            gd.addNumericField("Bias", settings.bias, 0);
            gd.addHelp(HelpUrls.getUrl(helpKey));
            gd.showDialog();
            if (gd.wasCanceled()) {
                return;
            }
            singleImage = true;
            settings.bias = Math.abs(gd.getNextNumber());
        } else {
            IJ.error(TITLE, "Single-image mode requires a stack");
            return;
        }
    }
    List<ImageSample> images;
    String inputDirectory = "";
    if (singleImage) {
        IJ.showStatus("Loading images...");
        images = getImages();
        if (images.size() == 0) {
            IJ.error(TITLE, "Not enough images for analysis");
            return;
        }
    } else {
        inputDirectory = IJ.getDirectory("Select image series ...");
        if (inputDirectory == null) {
            return;
        }
        final SeriesOpener series = new SeriesOpener(inputDirectory);
        series.setVariableSize(true);
        if (series.getNumberOfImages() < 3) {
            IJ.error(TITLE, "Not enough images in the selected directory");
            return;
        }
        if (!IJ.showMessageWithCancel(TITLE, String.format("Analyse %d images, first image:\n%s", series.getNumberOfImages(), series.getImageList()[0]))) {
            return;
        }
        IJ.showStatus("Loading images");
        images = getImages(series);
        if (images.size() < 3) {
            IJ.error(TITLE, "Not enough images for analysis");
            return;
        }
        if (images.get(0).exposure != 0) {
            IJ.error(TITLE, "First image in series must have exposure 0 (Bias image)");
            return;
        }
    }
    final boolean emMode = (arg != null && arg.contains("em"));
    GenericDialog gd = new GenericDialog(TITLE);
    gd.addMessage("Set the output options:");
    gd.addCheckbox("Show_table", settings.showTable);
    gd.addCheckbox("Show_charts", settings.showCharts);
    if (emMode) {
        // Ask the user for the camera gain ...
        gd.addMessage("Estimating the EM-gain requires the camera gain without EM readout enabled");
        gd.addNumericField("Camera_gain (Count/e-)", settings.cameraGain, 4);
    }
    if (emMode) {
        helpKey += "-em-ccd";
    }
    gd.addHelp(HelpUrls.getUrl(helpKey));
    gd.showDialog();
    if (gd.wasCanceled()) {
        return;
    }
    settings.showTable = gd.getNextBoolean();
    settings.showCharts = gd.getNextBoolean();
    if (emMode) {
        settings.cameraGain = gd.getNextNumber();
    }
    IJ.showStatus("Computing mean & variance");
    final double nImages = images.size();
    for (int i = 0; i < images.size(); i++) {
        IJ.showStatus(String.format("Computing mean & variance %d/%d", i + 1, images.size()));
        images.get(i).compute(singleImage, i / nImages, (i + 1) / nImages);
    }
    IJ.showProgress(1);
    IJ.showStatus("Computing results");
    // Allow user to input multiple bias images
    int start = 0;
    final Statistics biasStats = new Statistics();
    final Statistics noiseStats = new Statistics();
    final double bias;
    if (singleImage) {
        bias = settings.bias;
    } else {
        while (start < images.size()) {
            final ImageSample sample = images.get(start);
            if (sample.exposure == 0) {
                biasStats.add(sample.means);
                for (final PairSample pair : sample.samples) {
                    noiseStats.add(pair.variance);
                }
                start++;
            } else {
                break;
            }
        }
        bias = biasStats.getMean();
    }
    // Get the mean-variance data
    int total = 0;
    for (int i = start; i < images.size(); i++) {
        total += images.get(i).samples.size();
    }
    if (settings.showTable && total > 2000) {
        gd = new GenericDialog(TITLE);
        gd.addMessage("Table output requires " + total + " entries.\n \nYou may want to disable the table.");
        gd.addCheckbox("Show_table", settings.showTable);
        gd.showDialog();
        if (gd.wasCanceled()) {
            return;
        }
        settings.showTable = gd.getNextBoolean();
    }
    final TextWindow results = (settings.showTable) ? createResultsWindow() : null;
    double[] mean = new double[total];
    double[] variance = new double[mean.length];
    final Statistics gainStats = (singleImage) ? new StoredDataStatistics(total) : new Statistics();
    final WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = (singleImage) ? 0 : start, j = 0; i < images.size(); i++) {
        final StringBuilder sb = (settings.showTable) ? new StringBuilder() : null;
        final ImageSample sample = images.get(i);
        for (final PairSample pair : sample.samples) {
            if (j % 16 == 0) {
                IJ.showProgress(j, total);
            }
            mean[j] = pair.getMean();
            variance[j] = pair.variance;
            // Gain is in Count / e
            double gain = variance[j] / (mean[j] - bias);
            gainStats.add(gain);
            obs.add(mean[j], variance[j]);
            if (emMode) {
                gain /= (2 * settings.cameraGain);
            }
            if (sb != null) {
                sb.append(sample.title).append('\t');
                sb.append(sample.exposure).append('\t');
                sb.append(pair.slice1).append('\t');
                sb.append(pair.slice2).append('\t');
                sb.append(IJ.d2s(pair.mean1, 2)).append('\t');
                sb.append(IJ.d2s(pair.mean2, 2)).append('\t');
                sb.append(IJ.d2s(mean[j], 2)).append('\t');
                sb.append(IJ.d2s(variance[j], 2)).append('\t');
                sb.append(MathUtils.rounded(gain, 4)).append("\n");
            }
            j++;
        }
        if (results != null && sb != null) {
            results.append(sb.toString());
        }
    }
    IJ.showProgress(1);
    if (singleImage) {
        StoredDataStatistics stats = (StoredDataStatistics) gainStats;
        ImageJUtils.log(TITLE);
        if (emMode) {
            final double[] values = stats.getValues();
            MathArrays.scaleInPlace(0.5, values);
            stats = StoredDataStatistics.create(values);
        }
        if (settings.showCharts) {
            // Plot the gain over time
            final String title = TITLE + " Gain vs Frame";
            final Plot plot = new Plot(title, "Slice", "Gain");
            plot.addPoints(SimpleArrayUtils.newArray(gainStats.getN(), 1, 1.0), stats.getValues(), Plot.LINE);
            final PlotWindow pw = ImageJUtils.display(title, plot);
            // Show a histogram
            final String label = String.format("Mean = %s, Median = %s", MathUtils.rounded(stats.getMean()), MathUtils.rounded(stats.getMedian()));
            final WindowOrganiser wo = new WindowOrganiser();
            final PlotWindow pw2 = new HistogramPlotBuilder(TITLE, stats, "Gain").setRemoveOutliersOption(1).setPlotLabel(label).show(wo);
            if (wo.isNotEmpty()) {
                final Point point = pw.getLocation();
                point.y += pw.getHeight();
                pw2.setLocation(point);
            }
        }
        ImageJUtils.log("Single-image mode: %s camera", (emMode) ? "EM-CCD" : "Standard");
        final double gain = stats.getMedian();
        if (emMode) {
            final double totalGain = gain;
            final double emGain = totalGain / settings.cameraGain;
            ImageJUtils.log("  Gain = 1 / %s (Count/e-)", MathUtils.rounded(settings.cameraGain, 4));
            ImageJUtils.log("  EM-Gain = %s", MathUtils.rounded(emGain, 4));
            ImageJUtils.log("  Total Gain = %s (Count/e-)", MathUtils.rounded(totalGain, 4));
        } else {
            settings.cameraGain = gain;
            ImageJUtils.log("  Gain = 1 / %s (Count/e-)", MathUtils.rounded(settings.cameraGain, 4));
        }
    } else {
        IJ.showStatus("Computing fit");
        // Sort
        final int[] indices = rank(mean);
        mean = reorder(mean, indices);
        variance = reorder(variance, indices);
        // Compute optimal coefficients.
        // a - b x
        final double[] init = { 0, 1 / gainStats.getMean() };
        final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(2).withStartPoint(init);
        final double[] best = fitter.fit(obs.toList());
        // Construct the polynomial that best fits the data.
        final PolynomialFunction fitted = new PolynomialFunction(best);
        if (settings.showCharts) {
            // Plot mean verses variance. Gradient is gain in Count/e.
            final String title = TITLE + " results";
            final Plot plot = new Plot(title, "Mean", "Variance");
            final double[] xlimits = MathUtils.limits(mean);
            final double[] ylimits = MathUtils.limits(variance);
            double xrange = (xlimits[1] - xlimits[0]) * 0.05;
            if (xrange == 0) {
                xrange = 0.05;
            }
            double yrange = (ylimits[1] - ylimits[0]) * 0.05;
            if (yrange == 0) {
                yrange = 0.05;
            }
            plot.setLimits(xlimits[0] - xrange, xlimits[1] + xrange, ylimits[0] - yrange, ylimits[1] + yrange);
            plot.setColor(Color.blue);
            plot.addPoints(mean, variance, Plot.CROSS);
            plot.setColor(Color.red);
            plot.addPoints(new double[] { mean[0], mean[mean.length - 1] }, new double[] { fitted.value(mean[0]), fitted.value(mean[mean.length - 1]) }, Plot.LINE);
            ImageJUtils.display(title, plot);
        }
        final double avBiasNoise = Math.sqrt(noiseStats.getMean());
        ImageJUtils.log(TITLE);
        ImageJUtils.log("  Directory = %s", inputDirectory);
        ImageJUtils.log("  Bias = %s +/- %s (Count)", MathUtils.rounded(bias, 4), MathUtils.rounded(avBiasNoise, 4));
        ImageJUtils.log("  Variance = %s + %s * mean", MathUtils.rounded(best[0], 4), MathUtils.rounded(best[1], 4));
        if (emMode) {
            // The gradient is the observed gain of the noise.
            // In an EM-CCD there is a noise factor of 2.
            // Q. Is this true for a correct noise factor calibration:
            // double noiseFactor = (Read Noise EM-CCD) / (Read Noise CCD)
            // Em-gain is the observed gain divided by the noise factor multiplied by camera gain
            final double emGain = best[1] / (2 * settings.cameraGain);
            // Compute total gain
            final double totalGain = emGain * settings.cameraGain;
            final double readNoise = avBiasNoise / settings.cameraGain;
            // Effective noise is standard deviation of the bias image divided by the total gain (in
            // Count/e-)
            final double readNoiseE = avBiasNoise / totalGain;
            ImageJUtils.log("  Read Noise = %s (e-) [%s (Count)]", MathUtils.rounded(readNoise, 4), MathUtils.rounded(avBiasNoise, 4));
            ImageJUtils.log("  Gain = 1 / %s (Count/e-)", MathUtils.rounded(1 / settings.cameraGain, 4));
            ImageJUtils.log("  EM-Gain = %s", MathUtils.rounded(emGain, 4));
            ImageJUtils.log("  Total Gain = %s (Count/e-)", MathUtils.rounded(totalGain, 4));
            ImageJUtils.log("  Effective Read Noise = %s (e-) (Read Noise/Total Gain)", MathUtils.rounded(readNoiseE, 4));
        } else {
            // The gradient is the observed gain of the noise.
            settings.cameraGain = best[1];
            // Noise is standard deviation of the bias image divided by the gain (in Count/e-)
            final double readNoise = avBiasNoise / settings.cameraGain;
            ImageJUtils.log("  Read Noise = %s (e-) [%s (Count)]", MathUtils.rounded(readNoise, 4), MathUtils.rounded(avBiasNoise, 4));
            ImageJUtils.log("  Gain = 1 / %s (Count/e-)", MathUtils.rounded(1 / settings.cameraGain, 4));
        }
    }
    IJ.showStatus("");
}