Examples with HeightResultProcedure uk.ac.sussex.gdsc.smlm.results.procedures.HeightResultProcedure used on opensource projects

Example 1 with HeightResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.HeightResultProcedure in project GDSC-SMLM by aherbert.

the class PsfCreator method fitPsf.

/**
 * Fit the new PSF image and show a graph of the amplitude/width.
 *
 * @param psfStack the psf stack
 * @param loess the loess
 * @param cz the cz
 * @param averageRange the average range
 * @param fitCom the fit com
 * @return The width of the PSF in the z-centre
 */
private double fitPsf(ImageStack psfStack, LoessInterpolator loess, int cz, double averageRange, final double[][] fitCom) {
    IJ.showStatus("Fitting final PSF");
    // is not appropriate for a normalised PSF.
    if (fitConfig.getFitSolver() != FitSolver.LVM_LSE) {
        ImageJUtils.log("  " + FitProtosHelper.getName(fitConfig.getFitSolver()) + " is not appropriate for final PSF fitting.");
        ImageJUtils.log("  Switching to Least Square Estimation");
        fitConfig.setFitSolver(FitSolver.LVM_LSE);
        if (settings.getInteractiveMode()) {
            // This assumes the LVM does not need the calibration
            PeakFit.configureFitSolver(config, null, null, 0);
        }
    }
    // Update the box radius since this is used in the fitSpot method.
    boxRadius = psfStack.getWidth() / 2;
    final int x = boxRadius;
    final int y = boxRadius;
    final double shift = fitConfig.getCoordinateShiftFactor();
    // Scale the PSF
    final PSF.Builder localPsf = fitConfig.getPsf().toBuilder();
    for (int i = 0; i < localPsf.getParametersCount(); i++) {
        final PSFParameter param = localPsf.getParameters(i);
        if (param.getUnit() == PSFParameterUnit.DISTANCE) {
            final PSFParameter.Builder b = param.toBuilder();
            b.setValue(b.getValue() * settings.getMagnification());
            localPsf.setParameters(i, b);
        }
    }
    fitConfig.setPsf(localPsf.build());
    // Need to be updated after the widths have been set
    fitConfig.setCoordinateShiftFactor(shift);
    fitConfig.setBackgroundFitting(false);
    // Since the PSF will be normalised remove the camera calibration
    fitConfig.setCameraType(CameraType.CAMERA_TYPE_NA);
    fitConfig.setMinPhotons(0);
    fitConfig.setBias(0);
    fitConfig.setGain(1);
    // No complex filtering so we get a fit. It should be easy to fit anyway.
    fitConfig.setPrecisionThreshold(0);
    fitConfig.setDirectFilter(null);
    // fitConfig.setDisableSimpleFilter(true);
    // Use this for debugging the fit
    // fitConfig.setLog(uk.ac.sussex.gdsc.core.ij.ImageJPluginLoggerHelper.getDefaultLogger());
    final MemoryPeakResults results = fitSpot(psfStack, psfStack.getWidth(), psfStack.getHeight(), x, y);
    if (results.size() < 5) {
        ImageJUtils.log("  Final PSF: Not enough fit results %d", results.size());
        return 0;
    }
    // Get the results for the spot centre and width
    final double[] z = new double[results.size()];
    final double[] xCoord = new double[z.length];
    final double[] yCoord = new double[z.length];
    final double[] sd = new double[z.length];
    final double[] a = new double[z.length];
    // Set limits for the fit
    final float maxWidth = (float) (Math.max(fitConfig.getInitialXSd(), fitConfig.getInitialYSd()) * settings.getMagnification() * 4);
    // PSF is normalised to 1
    final float maxSignal = 2;
    final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
    wp.getWxWy();
    final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
    hp.getH();
    final Counter counter = new Counter();
    final Counter counterOk = new Counter();
    // We have fit the results so they will be in the preferred units
    results.forEach((PeakResultProcedure) peak -> {
        int index = counter.getAndIncrement();
        final float w = Math.max(wp.wx[index], wp.wy[index]);
        if (peak.getIntensity() > maxSignal || w > maxWidth) {
            return;
        }
        index = counterOk.getAndIncrement();
        z[index] = peak.getFrame();
        fitCom[0][peak.getFrame() - 1] = xCoord[index] = peak.getXPosition() - x;
        fitCom[1][peak.getFrame() - 1] = yCoord[index] = peak.getYPosition() - y;
        sd[index] = w;
        a[index] = hp.heights[index];
    });
    // Truncate
    final double[] z2 = Arrays.copyOf(z, counter.getCount());
    final double[] xCoord2 = Arrays.copyOf(xCoord, z2.length);
    final double[] yCoord2 = Arrays.copyOf(yCoord, z2.length);
    final double[] sd2 = Arrays.copyOf(sd, z2.length);
    final double[] a2 = Arrays.copyOf(a, z2.length);
    // Extract the average smoothed range from the individual fits
    final int r = (int) Math.ceil(averageRange / 2);
    int start = 0;
    int stop = z2.length - 1;
    for (int j = 0; j < z2.length; j++) {
        if (z2[j] > cz - r) {
            start = j;
            break;
        }
    }
    for (int j = z2.length; j-- > 0; ) {
        if (z2[j] < cz + r) {
            stop = j;
            break;
        }
    }
    // Extract xy centre coords and smooth
    double[] smoothX = new double[stop - start + 1];
    double[] smoothY = new double[smoothX.length];
    double[] smoothSd = new double[smoothX.length];
    double[] smoothA = new double[smoothX.length];
    final double[] newZ = new double[smoothX.length];
    int smoothCzIndex = 0;
    for (int j = start, k = 0; j <= stop; j++, k++) {
        smoothX[k] = xCoord2[j];
        smoothY[k] = yCoord2[j];
        smoothSd[k] = sd2[j];
        smoothA[k] = a2[j];
        newZ[k] = z2[j];
        if (newZ[k] == cz) {
            smoothCzIndex = k;
        }
    }
    smoothX = loess.smooth(newZ, smoothX);
    smoothY = loess.smooth(newZ, smoothY);
    smoothSd = loess.smooth(newZ, smoothSd);
    smoothA = loess.smooth(newZ, smoothA);
    // Update the widths and positions using the magnification
    final double scale = 1.0 / settings.getMagnification();
    for (int j = 0; j < xCoord2.length; j++) {
        xCoord2[j] *= scale;
        yCoord2[j] *= scale;
        sd2[j] *= scale;
    }
    for (int j = 0; j < smoothX.length; j++) {
        smoothX[j] *= scale;
        smoothY[j] *= scale;
        smoothSd[j] *= scale;
    }
    showPlots(z2, a2, newZ, smoothA, xCoord2, yCoord2, sd2, newZ, smoothX, smoothY, smoothSd, cz);
    // Store the data for replotting
    this.z = z2;
    this.amplitude = a2;
    this.smoothAz = newZ;
    this.smoothA = smoothA;
    this.xCoord = xCoord2;
    this.yCoord = yCoord2;
    this.sd = sd2;
    this.newZ = newZ;
    this.smoothX = smoothX;
    this.smoothY = smoothY;
    this.smoothSd = smoothSd;
    // maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
    return smoothSd[smoothCzIndex];
}

Example 2 with HeightResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.HeightResultProcedure in project GDSC-SMLM by aherbert.

the class PsfCreator method runUsingFitting.

private void runUsingFitting() {
    if (!showFittingDialog()) {
        return;
    }
    if (!loadConfiguration()) {
        return;
    }
    final BasePoint[] spots = getSpots(0, true);
    if (spots.length == 0) {
        IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
        return;
    }
    final ImageStack stack = getImageStack();
    final int width = imp.getWidth();
    final int height = imp.getHeight();
    final int currentSlice = imp.getSlice();
    // Adjust settings for a single maxima
    config.setIncludeNeighbours(false);
    final ArrayList<double[]> centres = new ArrayList<>(spots.length);
    final int iterations = 1;
    final LoessInterpolator loess = new LoessInterpolator(settings.getSmoothing(), iterations);
    // TODO - The fitting routine may not produce many points. In this instance the LOESS
    // interpolator
    // fails to smooth the data very well. A higher bandwidth helps this but perhaps
    // try a different smoothing method.
    // For each spot
    ImageJUtils.log(TITLE + ": " + imp.getTitle());
    ImageJUtils.log("Finding spot locations...");
    ImageJUtils.log("  %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
    final StoredDataStatistics averageSd = new StoredDataStatistics();
    final StoredDataStatistics averageA = new StoredDataStatistics();
    final Statistics averageRange = new Statistics();
    final MemoryPeakResults allResults = new MemoryPeakResults();
    allResults.setCalibration(fitConfig.getCalibration());
    allResults.setPsf(fitConfig.getPsf());
    allResults.setName(TITLE);
    allResults.setBounds(new Rectangle(0, 0, width, height));
    MemoryPeakResults.addResults(allResults);
    for (int n = 1; n <= spots.length; n++) {
        final BasePoint spot = spots[n - 1];
        final int x = (int) spot.getX();
        final int y = (int) spot.getY();
        final MemoryPeakResults results = fitSpot(stack, width, height, x, y);
        allResults.add(results);
        if (results.size() < 5) {
            ImageJUtils.log("  Spot %d: Not enough fit results %d", n, results.size());
            continue;
        }
        // Get the results for the spot centre and width
        final double[] z = new double[results.size()];
        final double[] xCoord = new double[z.length];
        final double[] yCoord = new double[z.length];
        final double[] sd;
        final double[] a;
        final Counter counter = new Counter();
        // We have fit the results so they will be in the preferred units
        results.forEach(new PeakResultProcedure() {

            @Override
            public void execute(PeakResult peak) {
                final int i = counter.getAndIncrement();
                z[i] = peak.getFrame();
                xCoord[i] = peak.getXPosition() - x;
                yCoord[i] = peak.getYPosition() - y;
            }
        });
        final WidthResultProcedure wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
        wp.getW();
        sd = SimpleArrayUtils.toDouble(wp.wx);
        final HeightResultProcedure hp = new HeightResultProcedure(results, IntensityUnit.COUNT);
        hp.getH();
        a = SimpleArrayUtils.toDouble(hp.heights);
        // Smooth the amplitude plot
        final double[] smoothA = loess.smooth(z, a);
        // Find the maximum amplitude
        int maximumIndex = findMaximumIndex(smoothA);
        // Find the range at a fraction of the max. This is smoothed to find the X/Y centre
        int start = 0;
        int stop = smoothA.length - 1;
        final double limit = smoothA[maximumIndex] * settings.getAmplitudeFraction();
        for (int j = 0; j < smoothA.length; j++) {
            if (smoothA[j] > limit) {
                start = j;
                break;
            }
        }
        for (int j = smoothA.length; j-- > 0; ) {
            if (smoothA[j] > limit) {
                stop = j;
                break;
            }
        }
        averageRange.add(stop - start + 1);
        // Extract xy centre coords and smooth
        double[] smoothX = new double[stop - start + 1];
        double[] smoothY = new double[smoothX.length];
        double[] smoothSd = new double[smoothX.length];
        final double[] newZ = new double[smoothX.length];
        for (int j = start, k = 0; j <= stop; j++, k++) {
            smoothX[k] = xCoord[j];
            smoothY[k] = yCoord[j];
            smoothSd[k] = sd[j];
            newZ[k] = z[j];
        }
        smoothX = loess.smooth(newZ, smoothX);
        smoothY = loess.smooth(newZ, smoothY);
        smoothSd = loess.smooth(newZ, smoothSd);
        // Since the amplitude is not very consistent move from this peak to the
        // lowest width which is the in-focus spot.
        maximumIndex = findMinimumIndex(smoothSd, maximumIndex - start);
        // Find the centre at the amplitude peak
        final double cx = smoothX[maximumIndex] + x;
        final double cy = smoothY[maximumIndex] + y;
        int cz = (int) newZ[maximumIndex];
        double csd = smoothSd[maximumIndex];
        double ca = smoothA[maximumIndex + start];
        // The average should weight the SD using the signal for each spot
        averageSd.add(smoothSd[maximumIndex]);
        averageA.add(ca);
        if (ignoreSpot(n, z, a, smoothA, xCoord, yCoord, sd, newZ, smoothX, smoothY, smoothSd, cx, cy, cz, csd)) {
            ImageJUtils.log("  Spot %d was ignored", n);
            continue;
        }
        // Store result - it may have been moved interactively
        maximumIndex += this.slice - cz;
        cz = (int) newZ[maximumIndex];
        csd = smoothSd[maximumIndex];
        ca = smoothA[maximumIndex + start];
        ImageJUtils.log("  Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f", n, cx, cy, cz, csd, ca);
        centres.add(new double[] { cx, cy, cz, csd, n });
    }
    if (settings.getInteractiveMode()) {
        imp.setSlice(currentSlice);
        imp.setOverlay(null);
        // Hide the amplitude and spot plots
        ImageJUtils.hide(TITLE_AMPLITUDE);
        ImageJUtils.hide(TITLE_PSF_PARAMETERS);
    }
    if (centres.isEmpty()) {
        final String msg = "No suitable spots could be identified";
        ImageJUtils.log(msg);
        IJ.error(TITLE, msg);
        return;
    }
    // Find the limits of the z-centre
    int minz = (int) centres.get(0)[2];
    int maxz = minz;
    for (final double[] centre : centres) {
        if (minz > centre[2]) {
            minz = (int) centre[2];
        } else if (maxz < centre[2]) {
            maxz = (int) centre[2];
        }
    }
    IJ.showStatus("Creating PSF image");
    // Create a stack that can hold all the data.
    final ImageStack psf = createStack(stack, minz, maxz, settings.getMagnification());
    // For each spot
    final Statistics stats = new Statistics();
    boolean ok = true;
    for (int i = 0; ok && i < centres.size(); i++) {
        final double increment = 1.0 / (stack.getSize() * centres.size());
        setProgress((double) i / centres.size());
        final double[] centre = centres.get(i);
        // Extract the spot
        final float[][] spot = new float[stack.getSize()][];
        Rectangle regionBounds = null;
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            final ImageExtractor ie = ImageExtractor.wrap((float[]) stack.getPixels(slice), width, height);
            if (regionBounds == null) {
                regionBounds = ie.getBoxRegionBounds((int) centre[0], (int) centre[1], boxRadius);
            }
            spot[slice - 1] = ie.crop(regionBounds);
        }
        if (regionBounds == null) {
            // Empty stack
            continue;
        }
        final int n = (int) centre[4];
        final float b = getBackground(n, spot);
        if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
            ImageJUtils.log("  Spot %d was ignored", n);
            continue;
        }
        stats.add(b);
        // Adjust the centre using the crop
        centre[0] -= regionBounds.x;
        centre[1] -= regionBounds.y;
        // This takes a long time so this should track progress
        ok = addToPsf(maxz, settings.getMagnification(), psf, centre, spot, regionBounds, increment, settings.getCentreEachSlice());
    }
    if (settings.getInteractiveMode()) {
        ImageJUtils.hide(TITLE_INTENSITY);
    }
    IJ.showProgress(1);
    if (!ok || stats.getN() == 0) {
        return;
    }
    final double avSd = getAverage(averageSd, averageA, 2);
    ImageJUtils.log("  Average background = %.2f, Av. SD = %s px", stats.getMean(), MathUtils.rounded(avSd, 4));
    normalise(psf, maxz, avSd * settings.getMagnification(), false);
    IJ.showProgress(1);
    psfImp = ImageJUtils.display(TITLE_PSF, psf);
    psfImp.setSlice(maxz);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    final double[][] fitCom = new double[2][psf.getSize()];
    Arrays.fill(fitCom[0], Double.NaN);
    Arrays.fill(fitCom[1], Double.NaN);
    final double fittedSd = fitPsf(psf, loess, maxz, averageRange.getMean(), fitCom);
    // Compute the drift in the PSF:
    // - Use fitted centre if available; otherwise find CoM for each frame
    // - express relative to the average centre
    final double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / settings.getMagnification());
    final double[] slice = SimpleArrayUtils.newArray(psf.getSize(), 1, 1.0);
    final String title = TITLE + " CoM Drift";
    final Plot plot = new Plot(title, "Slice", "Drift (nm)");
    plot.addLabel(0, 0, "Red = X; Blue = Y");
    // double[] limitsX = Maths.limits(com[0]);
    // double[] limitsY = Maths.limits(com[1]);
    final double[] limitsX = getLimits(com[0]);
    final double[] limitsY = getLimits(com[1]);
    plot.setLimits(1, psf.getSize(), Math.min(limitsX[0], limitsY[0]), Math.max(limitsX[1], limitsY[1]));
    plot.setColor(Color.red);
    plot.addPoints(slice, com[0], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[0]), Plot.LINE);
    plot.setColor(Color.blue);
    plot.addPoints(slice, com[1], Plot.DOT);
    plot.addPoints(slice, loess.smooth(slice, com[1]), Plot.LINE);
    ImageJUtils.display(title, plot);
    // TODO - Redraw the PSF with drift correction applied.
    // This means that the final image should have no drift.
    // This is relevant when combining PSF images. It doesn't matter too much for simulations
    // unless the drift is large.
    // Add Image properties containing the PSF details
    final double fwhm = getFwhm(psf, maxz);
    psfImp.setProperty("Info", ImagePsfHelper.toString(ImagePsfHelper.create(maxz, nmPerPixel / settings.getMagnification(), settings.getNmPerSlice(), stats.getN(), fwhm, createNote())));
    ImageJUtils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, " + "PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, MathUtils.rounded(nmPerPixel / settings.getMagnification(), 3), MathUtils.rounded(settings.getNmPerSlice(), 3), stats.getN(), MathUtils.rounded(fittedSd * nmPerPixel, 4), MathUtils.rounded(fwhm));
    createInteractivePlots(psf, maxz, nmPerPixel / settings.getMagnification(), fittedSd * nmPerPixel);
    IJ.showStatus("");
}

Example 3 with HeightResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.HeightResultProcedure in project GDSC-SMLM by aherbert.

the class SpotInspector method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        IJ.showStatus("");
        return;
    }
    // Check if the original image is open
    ImageSource source = results.getSource();
    if (source == null) {
        IJ.error(TITLE, "Unknown original source image");
        return;
    }
    source = source.getOriginal();
    source.setReadHint(ReadHint.NONSEQUENTIAL);
    if (!source.open()) {
        IJ.error(TITLE, "Cannot open original source image: " + source.toString());
        return;
    }
    final float stdDevMax = getStandardDeviation(results);
    if (stdDevMax < 0) {
        // TODO - Add dialog to get the initial peak width
        IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
        return;
    }
    // Rank spots
    rankedResults = new ArrayList<>(results.size());
    // Data for the sorting
    final PrecisionResultProcedure pp;
    if (settings.sortOrderIndex == 1) {
        pp = new PrecisionResultProcedure(results);
        pp.getPrecision();
    } else {
        pp = null;
    }
    // Build procedures to get:
    // Shift = position in pixels - originXY
    final StandardResultProcedure sp;
    if (settings.sortOrderIndex == 9) {
        sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
        sp.getXyr();
    } else {
        sp = null;
    }
    // SD = gaussian widths only for Gaussian PSFs
    final WidthResultProcedure wp;
    if (settings.sortOrderIndex >= 6 && settings.sortOrderIndex <= 8) {
        wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
        wp.getWxWy();
    } else {
        wp = null;
    }
    // Amplitude for Gaussian PSFs
    final HeightResultProcedure hp;
    if (settings.sortOrderIndex == 2) {
        hp = new HeightResultProcedure(results, IntensityUnit.PHOTON);
        hp.getH();
    } else {
        hp = null;
    }
    final Counter c = new Counter();
    results.forEach((PeakResultProcedure) result -> {
        final float[] score = getScore(result, c.getAndIncrement(), pp, sp, wp, hp, stdDevMax);
        rankedResults.add(new PeakResultRank(result, score[0], score[1]));
    });
    Collections.sort(rankedResults, PeakResultRank::compare);
    // Prepare results table
    final ImageJTablePeakResults table = new ImageJTablePeakResults(false, results.getName(), true);
    table.copySettings(results);
    table.setTableTitle(TITLE);
    table.setAddCounter(true);
    table.setShowZ(results.is3D());
    // TODO - Add to settings
    table.setShowFittingData(true);
    table.setShowNoiseData(true);
    if (settings.showCalibratedValues) {
        table.setDistanceUnit(DistanceUnit.NM);
        table.setIntensityUnit(IntensityUnit.PHOTON);
    } else {
        table.setDistanceUnit(DistanceUnit.PIXEL);
        table.setIntensityUnit(IntensityUnit.COUNT);
    }
    table.begin();
    // Add a mouse listener to jump to the frame for the clicked line
    textPanel = table.getResultsWindow().getTextPanel();
    // We must ignore old instances of this class from the mouse listeners
    id = currentId.incrementAndGet();
    textPanel.addMouseListener(new MouseAdapter() {

        @Override
        public void mouseClicked(MouseEvent event) {
            SpotInspector.this.mouseClicked(event);
        }
    });
    // Add results to the table
    int count = 0;
    for (final PeakResultRank rank : rankedResults) {
        rank.rank = count++;
        table.add(rank.peakResult);
    }
    table.end();
    if (settings.plotScore || settings.plotHistogram) {
        // Get values for the plots
        float[] xValues = null;
        float[] yValues = null;
        double yMin;
        double yMax;
        int spotNumber = 0;
        xValues = new float[rankedResults.size()];
        yValues = new float[xValues.length];
        for (final PeakResultRank rank : rankedResults) {
            xValues[spotNumber] = spotNumber + 1;
            yValues[spotNumber++] = recoverScore(rank.score);
        }
        // Set the min and max y-values using 1.5 x IQR
        final DescriptiveStatistics stats = new DescriptiveStatistics();
        for (final float v : yValues) {
            stats.addValue(v);
        }
        if (settings.removeOutliers) {
            final double lower = stats.getPercentile(25);
            final double upper = stats.getPercentile(75);
            final double iqr = upper - lower;
            yMin = Math.max(lower - iqr, stats.getMin());
            yMax = Math.min(upper + iqr, stats.getMax());
            IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
        } else {
            yMin = stats.getMin();
            yMax = stats.getMax();
            IJ.log(String.format("Data range: %f - %f", yMin, yMax));
        }
        plotScore(xValues, yValues, yMin, yMax);
        plotHistogram(yValues);
    }
    // Extract spots into a stack
    final int w = source.getWidth();
    final int h = source.getHeight();
    final int size = 2 * settings.radius + 1;
    final ImageStack spots = new ImageStack(size, size, rankedResults.size());
    // To assist the extraction of data from the image source, process them in time order to allow
    // frame caching. Then set the appropriate slice in the result stack
    Collections.sort(rankedResults, (o1, o2) -> Integer.compare(o1.peakResult.getFrame(), o2.peakResult.getFrame()));
    for (final PeakResultRank rank : rankedResults) {
        final PeakResult r = rank.peakResult;
        // Extract image
        // Note that the coordinates are relative to the middle of the pixel (0.5 offset)
        // so do not round but simply convert to int
        final int x = (int) (r.getXPosition());
        final int y = (int) (r.getYPosition());
        // Extract a region but crop to the image bounds
        int minX = x - settings.radius;
        int minY = y - settings.radius;
        final int maxX = Math.min(x + settings.radius + 1, w);
        final int maxY = Math.min(y + settings.radius + 1, h);
        int padX = 0;
        int padY = 0;
        if (minX < 0) {
            padX = -minX;
            minX = 0;
        }
        if (minY < 0) {
            padY = -minY;
            minY = 0;
        }
        final int sizeX = maxX - minX;
        final int sizeY = maxY - minY;
        float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
        // Prevent errors with missing data
        if (data == null) {
            data = new float[sizeX * sizeY];
        }
        ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
        // Pad if necessary, i.e. the crop is too small for the stack
        if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
            final ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
            spotIp2.insert(spotIp, padX, padY);
            spotIp = spotIp2;
        }
        final int slice = rank.rank + 1;
        spots.setPixels(spotIp.getPixels(), slice);
        spots.setSliceLabel(MathUtils.rounded(rank.originalScore), slice);
    }
    source.close();
    // Reset to the rank order
    Collections.sort(rankedResults, PeakResultRank::compare);
    final ImagePlus imp = ImageJUtils.display(TITLE, spots);
    imp.setRoi((PointRoi) null);
    // Make bigger
    for (int i = 10; i-- > 0; ) {
        imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
    }
}