Examples with CalibrationReader uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader used on opensource projects

Example 16 with CalibrationReader

use of uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader in project GDSC-SMLM by aherbert.

the class ResultsImageSampler method getSample.

/**
 * Gets the sample image. The image is a stack of the samples with an overlay of the localisation
 * positions. The info property is set with details of the localisations and the image is
 * calibrated.
 *
 * @param countNo the number of samples with no localisations
 * @param countLow the number of samples with low localisations
 * @param countHigh the number of samples with high localisations
 * @return the sample image (could be null if no samples were made)
 */
public ImagePlus getSample(int countNo, int countLow, int countHigh) {
    final ImageStack out = new ImageStack(size, size);
    if (!isValid()) {
        return null;
    }
    sampleList.clear();
    // empty
    for (final int i : RandomUtils.sample(countNo, no.length, rng)) {
        sampleList.add(ResultsSample.createEmpty(no[i]));
    }
    // low
    for (final int i : RandomUtils.sample(countLow, lower, rng)) {
        sampleList.add(data[i]);
    }
    // high
    for (final int i : RandomUtils.sample(countHigh, upper, rng)) {
        sampleList.add(data[i + lower]);
    }
    if (sampleList.isEmpty()) {
        return null;
    }
    double nmPerPixel = 1;
    if (results.hasCalibration()) {
        final CalibrationReader calibration = results.getCalibrationReader();
        if (calibration.hasNmPerPixel()) {
            nmPerPixel = calibration.getNmPerPixel();
        }
    }
    // Sort descending by number in the frame
    final ResultsSample[] samples = sampleList.toArray(new ResultsSample[0]);
    Arrays.sort(samples, ReverseCountComparator.INSTANCE);
    final int[] xyz = new int[3];
    final Rectangle stackBounds = new Rectangle(stack.getWidth(), stack.getHeight());
    final Overlay overlay = new Overlay();
    final StringBuilder sb = new StringBuilder();
    if (nmPerPixel == 1) {
        sb.append("Sample X Y Z Signal\n");
    } else {
        sb.append("Sample X(nm) Y(nm) Z(nm) Signal\n");
    }
    for (final ResultsSample sample : samples) {
        getXyz(sample.index, xyz);
        // Construct the region to extract
        Rectangle target = new Rectangle(xyz[0], xyz[1], size, size);
        target = target.intersection(stackBounds);
        if (target.width == 0 || target.height == 0) {
            continue;
        }
        // Extract the frame
        final int slice = xyz[2];
        final ImageProcessor ip = stack.getProcessor(slice);
        // Cut out the desired pixels (some may be blank if the block overruns the source image)
        final ImageProcessor ip2 = ip.createProcessor(size, size);
        for (int y = 0; y < target.height; y++) {
            for (int x = 0, i = y * size, index = (y + target.y) * ip.getWidth() + target.x; x < target.width; x++, i++, index++) {
                ip2.setf(i, ip.getf(index));
            }
        }
        final int sampleSize = sample.size();
        if (sampleSize > 0) {
            final float[] ox = new float[sampleSize];
            final float[] oy = new float[sampleSize];
            final int position = out.getSize() + 1;
            // Create an ROI with the localisations
            for (int i = 0; i < sampleSize; i++) {
                final PeakResult p = sample.list.get(i);
                ox[i] = p.getXPosition() - xyz[0];
                oy[i] = p.getYPosition() - xyz[1];
                sb.append(position).append(' ');
                sb.append(MathUtils.rounded(ox[i] * nmPerPixel)).append(' ');
                sb.append(MathUtils.rounded(oy[i] * nmPerPixel)).append(' ');
                sb.append(MathUtils.rounded(p.getZPosition() * nmPerPixel)).append(' ');
                sb.append(MathUtils.rounded(p.getIntensity())).append('\n');
            }
            final PointRoi roi = new OffsetPointRoi(ox, oy, sampleSize);
            roi.setPosition(position);
            overlay.add(roi);
        }
        out.addSlice(String.format("Frame=%d @ %d,%d px (n=%d)", slice, xyz[0], xyz[1], sampleSize), ip2.getPixels());
    }
    if (out.getSize() == 0) {
        return null;
    }
    final ImagePlus imp = new ImagePlus("Sample", out);
    imp.setOverlay(overlay);
    // Note: Only the info property can be saved to a TIFF file
    imp.setProperty("Info", sb.toString());
    if (nmPerPixel != 1) {
        final ij.measure.Calibration cal = new ij.measure.Calibration();
        cal.setUnit("nm");
        cal.pixelHeight = cal.pixelWidth = nmPerPixel;
        imp.setCalibration(cal);
    }
    return imp;
}

Example 17 with CalibrationReader

use of uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader in project GDSC-SMLM by aherbert.

the class PeakFit method addCameraOptions.

/**
 * Adds the camera options.
 *
 * @param gd the dialog
 * @param options the options
 * @param calibrationProvider the calibration provider
 */
public static void addCameraOptions(final ExtendedGenericDialog gd, final int options, final CalibrationProvider calibrationProvider) {
    final CalibrationReader calibration = new CalibrationReader(calibrationProvider.getCalibration());
    gd.addChoice("Camera_type", SettingsManager.getCameraTypeNames(), CalibrationProtosHelper.getName(calibration.getCameraType()), new OptionListener<Integer>() {

        @Override
        public boolean collectOptions(Integer field) {
            final CalibrationWriter calibration = new CalibrationWriter(calibrationProvider.getCalibration());
            final CameraType t = SettingsManager.getCameraTypeValues()[field];
            if (calibration.getCameraType() != t) {
                calibration.setCameraType(t);
                calibrationProvider.saveCalibration(calibration.getCalibration());
            }
            return collectOptions(false);
        }

        @Override
        public boolean collectOptions() {
            return collectOptions(true);
        }

        private boolean collectOptions(boolean silent) {
            final CalibrationWriter calibration = new CalibrationWriter(calibrationProvider.getCalibration());
            final ExtendedGenericDialog egd = new ExtendedGenericDialog("Camera type options", null);
            if (calibration.isCcdCamera()) {
                egd.addNumericField("Camera_bias", calibration.getBias(), 2, 6, "Count");
                if (BitFlagUtils.anyNotSet(options, FLAG_NO_GAIN)) {
                    egd.addNumericField("Gain", calibration.getCountPerPhoton(), 4, 6, "Count/photon");
                }
                if (BitFlagUtils.anyNotSet(options, FLAG_NO_READ_NOISE)) {
                    egd.addNumericField("Read_noise", calibration.getReadNoise(), 4, 6, "Count");
                }
                if (BitFlagUtils.areSet(options, FLAG_QUANTUM_EFFICIENCY)) {
                    egd.addNumericField("Quantum_efficiency", calibration.getQuantumEfficiency(), 4, 6, "electron/photon");
                }
            } else if (calibration.isScmos()) {
                final String[] models = CameraModelManager.listCameraModels(true);
                egd.addChoice("Camera_model_name", models, calibration.getCameraModelName());
                if (BitFlagUtils.areSet(options, FLAG_QUANTUM_EFFICIENCY)) {
                    egd.addNumericField("Quantum_efficiency", calibration.getQuantumEfficiency(), 4, 6, "electron/photon");
                }
            } else {
                IJ.error("Unsupported camera type " + CalibrationProtosHelper.getName(calibration.getCameraType()));
                return false;
            }
            egd.setSilent(silent);
            egd.showDialog(true, gd);
            if (egd.wasCanceled()) {
                return false;
            }
            final Calibration old = calibration.getCalibration();
            if (calibration.isCcdCamera()) {
                calibration.setBias(Math.abs(egd.getNextNumber()));
                if (BitFlagUtils.anyNotSet(options, FLAG_NO_GAIN)) {
                    calibration.setCountPerPhoton(Math.abs(egd.getNextNumber()));
                }
                if (BitFlagUtils.anyNotSet(options, FLAG_NO_READ_NOISE)) {
                    calibration.setReadNoise(Math.abs(egd.getNextNumber()));
                }
                if (BitFlagUtils.areSet(options, FLAG_QUANTUM_EFFICIENCY)) {
                    calibration.setQuantumEfficiency(Math.abs(egd.getNextNumber()));
                }
            } else if (calibration.isScmos()) {
                // Note: Since this does not go through the FitConfiguration object the
                // camera model is not invalidated. However any code using this function
                // should later call configureFitSolver(...) which will set the camera model
                // using the camera model name.
                calibration.setCameraModelName(egd.getNextChoice());
                if (BitFlagUtils.areSet(options, FLAG_QUANTUM_EFFICIENCY)) {
                    calibration.setQuantumEfficiency(Math.abs(egd.getNextNumber()));
                }
            }
            final Calibration current = calibration.getCalibration();
            final boolean changed = !old.equals(current);
            if (changed) {
                calibrationProvider.saveCalibration(current);
            }
            return changed;
        }
    });
}

Example 18 with CalibrationReader

use of uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader in project GDSC-SMLM by aherbert.

the class PeakFit method showDialog.

private int showDialog(ImagePlus imp) {
    // Executing as an ImageJ plugin.
    // Load the settings
    resultsSettings = SettingsManager.readResultsSettings(0).toBuilder();
    // Settings are within the FitEngineSettings
    config = SettingsManager.readFitEngineConfiguration(0);
    fitConfig = config.getFitConfiguration();
    settings = Settings.load();
    if (simpleFit) {
        return showSimpleDialog();
    }
    // Note: The bounds are not set when running in the fit maxima option (since all candidates
    // have been identified already the crop is not required).
    final boolean isCrop = (bounds != null && imp != null && (bounds.width < imp.getWidth() || bounds.height < imp.getHeight()));
    // Some options are not always needed
    if (extraOptions || isCrop) {
        extraSettings = ExtraSettings.load();
        ignoreBoundsForNoise = extraSettings.optionIgnoreBoundsForNoise;
    }
    if (!extraOptions) {
        resultsSettings.getResultsImageSettingsBuilder().setRollingWindowSize(0);
        fitConfig.setBackgroundFitting(true);
        fitConfig.setNoise(0);
        config.setNoiseMethod(NoiseEstimatorMethod.QUICK_RESIDUALS_LEAST_MEAN_OF_SQUARES);
    }
    final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
    String helpKey;
    if (maximaIdentification) {
        helpKey = "spot-finder";
        gd.addMessage("Identify candidate maxima");
    } else {
        helpKey = "peak-fit";
        gd.addMessage("Fit 2D Gaussian to identified maxima");
    }
    // Note: Currently is is not useful to append "-series" when running for a series image
    // source since the params in this dialog do not concern the image input.
    gd.addHelp(HelpUrls.getUrl(helpKey));
    final String[] templates = ConfigurationTemplate.getTemplateNames(true);
    gd.addChoice("Template", templates, templates[0]);
    final CalibrationReader calibration = fitConfig.getCalibrationReader();
    addCameraOptions(gd, 0, fitConfig);
    gd.addNumericField("Calibration", calibration.getNmPerPixel(), 2, 6, "nm/px");
    gd.addNumericField("Exposure_time", calibration.getExposureTime(), 2, 6, "ms");
    if (isCrop) {
        gd.addCheckbox("Ignore_bounds_for_noise", ignoreBoundsForNoise);
    }
    final FitConfigurationProvider fitConfigurationProvider = () -> fitConfig;
    final FitEngineConfigurationProvider fitEngineConfigurationProvider = () -> config;
    addPsfOptions(gd, fitConfigurationProvider);
    addDataFilterOptions(gd, fitEngineConfigurationProvider);
    addSearchOptions(gd, fitEngineConfigurationProvider);
    addBorderOptions(gd, fitEngineConfigurationProvider);
    addFittingOptions(gd, fitEngineConfigurationProvider);
    if (extraOptions && !fitMaxima) {
        gd.addCheckbox("Interlaced_data", extraSettings.interlacedData);
        gd.addSlider("Integrate_frames", 1, 5, extraSettings.integrateFrames);
    }
    // Special case top get the slider since the GenericDialog does not provide access to this.
    Scrollbar sliderCoordinateShiftFactor = null;
    final boolean isShowGenericDialog = ImageJUtils.isShowGenericDialog();
    if (!maximaIdentification) {
        gd.addMessage("--- Gaussian fitting ---");
        gd.addChoice("Fit_solver", SettingsManager.getFitSolverNames(), FitProtosHelper.getName(fitConfig.getFitSolver()));
        if (extraOptions) {
            gd.addCheckbox("Fit_background", fitConfig.isBackgroundFitting());
        }
        // Parameters specific to each Fit solver are collected in a second dialog
        gd.addNumericField("Fail_limit", config.getFailuresLimit(), 0);
        gd.addNumericField("Pass_rate", config.getPassRate(), 2);
        gd.addCheckbox("Include_neighbours", config.isIncludeNeighbours());
        gd.addSlider("Neighbour_height", 0.01, 1, config.getNeighbourHeightThreshold());
        gd.addSlider("Residuals_threshold", 0.01, 1, config.getResidualsThreshold());
        addDuplicateDistanceOptions(gd, fitEngineConfigurationProvider);
        gd.addMessage("--- Peak filtering ---\nDiscard fits that shift; are too low; or expand/contract");
        gd.addCheckbox("Smart_filter", fitConfig.isSmartFilter());
        gd.addCheckbox("Disable_simple_filter", fitConfig.isDisableSimpleFilter());
        gd.addSlider("Shift_factor", 0.0, 2.5, fitConfig.getCoordinateShiftFactor());
        if (isShowGenericDialog) {
            sliderCoordinateShiftFactor = gd.getLastScrollbar();
        }
        gd.addNumericField("Signal_strength", fitConfig.getSignalStrength(), 2);
        gd.addNumericField("Min_photons", fitConfig.getMinPhotons(), 0);
        if (extraOptions) {
            gd.addNumericField("Noise", fitConfig.getNoise(), 2);
            gd.addChoice("Noise_method", SettingsManager.getNoiseEstimatorMethodNames(), config.getNoiseMethod().ordinal());
        }
        gd.addSlider("Min_width_factor", 0, 0.99, fitConfig.getMinWidthFactor());
        gd.addSlider("Width_factor", 1, 4.5, fitConfig.getMaxWidthFactor());
        addPrecisionOptions(gd, fitConfigurationProvider);
    // Q. Add dynamically displayed options for z-filtering here?
    }
    gd.addMessage("--- Results ---");
    gd.addCheckbox("Log_progress", resultsSettings.getLogProgress());
    if (!maximaIdentification) {
        gd.addCheckbox("Show_deviations", resultsSettings.getShowDeviations());
    }
    ResultsManager.addTableResultsOptions(gd, resultsSettings);
    ResultsManager.addImageResultsOptions(gd, resultsSettings, (extraOptions) ? ResultsManager.FLAG_EXTRA_OPTIONS : 0);
    if (extraOptions) {
        gd.addCheckbox("Show_processed_frames", extraSettings.showProcessedFrames);
    }
    ResultsManager.addFileResultsOptions(gd, resultsSettings, ResultsManager.FLAG_RESULTS_DIRECTORY);
    ResultsManager.addInMemoryResultsOptions(gd, resultsSettings);
    if (extraOptions) {
        gd.addMessage("--- Misc ---");
        gd.addSlider("Fraction_of_threads", 0.1, 1, settings.fractionOfThreads);
    }
    // Add a mouse listener to the config file field
    if (isShowGenericDialog) {
        new ItemDialogListener(sliderCoordinateShiftFactor).attach(gd, isCrop);
    }
    gd.showDialog();
    if (gd.wasCanceled() || !readDialog(gd, isCrop)) {
        return DONE;
    }
    if (imp != null) {
        // Store whether the user selected to process all the images.
        final int flags = IJ.setupDialog(imp, pluginFlags);
        // Check if cancelled
        if ((flags & DONE) != 0) {
            return DONE;
        }
        if ((flags & DOES_STACKS) == 0) {
            // Save the slice number for the overlay
            singleFrame = imp.getCurrentSlice();
            // Account for interlaced data
            if (extraSettings.interlacedData) {
                int start = singleFrame;
                // Calculate the first frame that is not skipped
                while (ignoreFrame(start) && start > extraSettings.dataStart) {
                    start--;
                }
                if (start < extraSettings.dataStart) {
                    log("The current frame (%d) is before the start of the interlaced data", singleFrame);
                    return DONE;
                }
                if (start != singleFrame) {
                    log("Updated the current frame (%d) to a valid interlaced data frame (%d)", singleFrame, start);
                }
                singleFrame = start;
            }
            // Account for integrated frames
            int endFrame = singleFrame;
            if (extraSettings.integrateFrames > 1) {
                int totalFrames = 1;
                while (totalFrames < extraSettings.integrateFrames) {
                    endFrame++;
                    if (!ignoreFrame(endFrame)) {
                        totalFrames++;
                    }
                }
                log("Updated the image end frame (%d) to %d allow %d integrated frames", singleFrame, endFrame, extraSettings.integrateFrames);
            }
            // Create a new image source with the correct frames
            setSource(new IJImageSource(imp, singleFrame, endFrame - singleFrame));
            // Store the image so the results can be added as an overlay
            this.imp = imp;
            this.imp.setOverlay(null);
        }
    }
    // Allow interlaced data by wrapping the image source
    if (extraSettings.interlacedData) {
        setSource(new InterlacedImageSource(this.source, extraSettings.dataStart, extraSettings.dataBlock, extraSettings.dataSkip));
    }
    // Allow frame aggregation by wrapping the image source
    if (extraSettings.integrateFrames > 1) {
        setSource(new AggregatedImageSource(this.source, extraSettings.integrateFrames));
    }
    // Ask if the user wants to log progress on multiple frame images
    if (resultsSettings.getLogProgress() && source.getFrames() > 1) {
        final ExtendedGenericDialog egd = new ExtendedGenericDialog(TITLE);
        egd.addMessage("Warning: Log progress on multiple-frame image will be slow");
        egd.addCheckbox("Log_progress", true);
        egd.showDialog();
        if (egd.wasCanceled()) {
            return DONE;
        }
        if (!egd.getNextBoolean()) {
            resultsSettings.setLogProgress(false);
            SettingsManager.writeSettings(resultsSettings.build());
        }
    }
    // single call to be made.
    return pluginFlags;
}

Example 19 with CalibrationReader

use of uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader in project GDSC-SMLM by aherbert.

the class PeakFit method configureDataFilter.

/**
 * Show a dialog to configure the data filter. The data filter type and the first data filter must
 * ALREADY be set in the configuration. The subsequent filters are then configured, e.g. for
 * difference and jury filters.
 *
 * <p>The updated settings are saved to the settings file. An error message is shown if the dialog
 * is cancelled or the configuration is invalid.
 *
 * <p>If the configuration is for a per-pixel camera type (e.g. sCMOS) then the camera model will
 * be loaded using the configured camera model name. This will be used to validate the filter to
 * check the filter supports the per-pixel camera type.
 *
 * @param config the config
 * @param flags the flags
 * @return True if the configuration succeeded
 */
public static boolean configureDataFilter(final FitEngineConfiguration config, int flags) {
    int numberOfFilters = 1;
    int filterCount;
    switch(config.getDataFilterType()) {
        case JURY:
            filterCount = Integer.MAX_VALUE;
            break;
        case DIFFERENCE:
            filterCount = 2;
            break;
        case SINGLE:
        default:
            filterCount = 1;
    }
    final String[] filterNames = SettingsManager.getDataFilterMethodNames();
    final DataFilterMethod[] filterValues = SettingsManager.getDataFilterMethodValues();
    // Check we have at least the first filter.
    if (config.getDataFiltersCount() == 0) {
        throw new IllegalStateException("No primary filter is configured");
    }
    final FitEngineConfigurationProvider fitEngineConfigurationProvider = () -> config;
    for (int i = 1; i < filterCount; i++) {
        final int filter = i + 1;
        final int ii = i;
        final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
        if (filter == filterCount) {
            // This is maximum filter count so no continue option
            ImageJUtils.addMessage(gd, "Configure the %s filter.", FitProtosHelper.getName(config.getDataFilterType()));
        } else {
            gd.enableYesNoCancel("Add", "Continue");
            ImageJUtils.addMessage(gd, "Configure the %s filter.\nClick continue to proceed with the current set of %d.", FitProtosHelper.getName(config.getDataFilterType()), i);
        }
        final String fieldName = "Spot_filter" + filter;
        if (IJ.isMacro()) {
            // Use blank default value so bad macro parameters return nothing
            gd.addStringField(fieldName, "");
        } else {
            gd.addChoice(fieldName, filterNames, filterNames[config.getDataFilterMethod(ii, config.getDataFilterMethod(ii - 1)).ordinal()]);
        }
        addRelativeParameterOptions(gd, new RelativeParameterProvider(0, 4.5, "Smoothing" + filter, fitEngineConfigurationProvider, true) {

            @Override
            void setAbsolute(boolean absolute) {
                // Get the current settings
                final FitEngineConfiguration c = fitEngineConfigurationProvider.getFitEngineConfiguration();
                final DataFilterMethod m = c.getDataFilterMethod(ii);
                final double smooth = c.getDataFilterParameter(ii).getValue();
                // Reset with the new absolute value
                c.setDataFilter(m, smooth, absolute, ii);
            }

            @Override
            boolean isAbsolute() {
                final FitEngineConfiguration c = fitEngineConfigurationProvider.getFitEngineConfiguration();
                return c.getDataFilterParameterAbsolute(ii, c.getDataFilterParameterAbsolute(ii - 1));
            }

            @Override
            double getValue() {
                final FitEngineConfiguration c = fitEngineConfigurationProvider.getFitEngineConfiguration();
                return c.getDataFilterParameterValue(ii, c.getDataFilterParameterValue(ii - 1));
            }
        });
        gd.showDialog();
        if (gd.wasCanceled()) {
            return false;
        }
        if (gd.wasOKed()) {
            int filterIndex = -1;
            if (IJ.isMacro()) {
                final String filterName = gd.getNextString();
                for (int j = 0; j < filterNames.length; j++) {
                    if (filterNames[j].equals(filterName)) {
                        filterIndex = j;
                        break;
                    }
                }
                if (filterIndex < 0) {
                    break;
                }
            } else {
                filterIndex = gd.getNextChoiceIndex();
            }
            // Note: The absolute flag is set in extra options
            config.setDataFilter(filterValues[filterIndex], Math.abs(gd.getNextNumber()), i);
            gd.collectOptions();
            numberOfFilters++;
        } else {
            break;
        }
    }
    config.setNumberOfFilters(numberOfFilters);
    if (BitFlagUtils.anyNotSet(flags, FLAG_NO_SAVE)) {
        saveFitEngineSettings(config);
    }
    final FitConfiguration fitConfig = config.getFitConfiguration();
    final CalibrationReader calibration = fitConfig.getCalibrationReader();
    if (calibration.isScmos()) {
        fitConfig.setCameraModel(CameraModelManager.load(fitConfig.getCameraModelName()));
    }
    try {
        config.createSpotFilter();
    } catch (final IllegalStateException ex) {
        IJ.error(TITLE, ex.getMessage());
        return false;
    }
    return true;
}

Example 20 with CalibrationReader

use of uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader 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("");
}