Examples with ImageExtractor gdsc.core.utils.ImageExtractor used on opensource projects

Example 1 with ImageExtractor

use of gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PSFCreator method run.

/*
	 * (non-Javadoc)
	 * 
	 * @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor)
	 */
public void run(ImageProcessor ip) {
    loadConfiguration();
    BasePoint[] spots = getSpots();
    if (spots.length == 0) {
        IJ.error(TITLE, "No spots without neighbours within " + (boxRadius * 2) + "px");
        return;
    }
    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);
    fitConfig.setDuplicateDistance(0);
    ArrayList<double[]> centres = new ArrayList<double[]>(spots.length);
    int iterations = 1;
    LoessInterpolator loess = new LoessInterpolator(smoothing, 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
    Utils.log(TITLE + ": " + imp.getTitle());
    Utils.log("Finding spot locations...");
    Utils.log("  %d spot%s without neighbours within %dpx", spots.length, ((spots.length == 1) ? "" : "s"), (boxRadius * 2));
    StoredDataStatistics averageSd = new StoredDataStatistics();
    StoredDataStatistics averageA = new StoredDataStatistics();
    Statistics averageRange = new Statistics();
    MemoryPeakResults allResults = new MemoryPeakResults();
    allResults.setName(TITLE);
    allResults.setBounds(new Rectangle(0, 0, width, height));
    MemoryPeakResults.addResults(allResults);
    for (int n = 1; n <= spots.length; n++) {
        BasePoint spot = spots[n - 1];
        final int x = (int) spot.getX();
        final int y = (int) spot.getY();
        MemoryPeakResults results = fitSpot(stack, width, height, x, y);
        allResults.addAllf(results.getResults());
        if (results.size() < 5) {
            Utils.log("  Spot %d: Not enough fit results %d", n, results.size());
            continue;
        }
        // Get the results for the spot centre and width
        double[] z = new double[results.size()];
        double[] xCoord = new double[z.length];
        double[] yCoord = new double[z.length];
        double[] sd = new double[z.length];
        double[] a = new double[z.length];
        int i = 0;
        for (PeakResult peak : results.getResults()) {
            z[i] = peak.getFrame();
            xCoord[i] = peak.getXPosition() - x;
            yCoord[i] = peak.getYPosition() - y;
            sd[i] = FastMath.max(peak.getXSD(), peak.getYSD());
            a[i] = peak.getAmplitude();
            i++;
        }
        // Smooth the amplitude plot
        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, stop = smoothA.length - 1;
        double limit = smoothA[maximumIndex] * amplitudeFraction;
        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];
        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
        double cx = smoothX[maximumIndex] + x;
        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)) {
            Utils.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];
        Utils.log("  Spot %d => x=%.2f, y=%.2f, z=%d, sd=%.2f, A=%.2f\n", n, cx, cy, cz, csd, ca);
        centres.add(new double[] { cx, cy, cz, csd, n });
    }
    if (interactiveMode) {
        imp.setSlice(currentSlice);
        imp.setOverlay(null);
        // Hide the amplitude and spot plots
        Utils.hide(TITLE_AMPLITUDE);
        Utils.hide(TITLE_PSF_PARAMETERS);
    }
    if (centres.isEmpty()) {
        String msg = "No suitable spots could be identified centres";
        Utils.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 (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.
    ImageStack psf = createStack(stack, minz, maxz, magnification);
    // For each spot
    Statistics stats = new Statistics();
    boolean ok = true;
    for (int i = 0; ok && i < centres.size(); i++) {
        double progress = (double) i / centres.size();
        final double increment = 1.0 / (stack.getSize() * centres.size());
        IJ.showProgress(progress);
        double[] centre = centres.get(i);
        // Extract the spot
        float[][] spot = new float[stack.getSize()][];
        Rectangle regionBounds = null;
        for (int slice = 1; slice <= stack.getSize(); slice++) {
            ImageExtractor ie = new ImageExtractor((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);
        }
        int n = (int) centre[4];
        final float b = getBackground(n, spot);
        if (!subtractBackgroundAndWindow(spot, b, regionBounds.width, regionBounds.height, centre, loess)) {
            Utils.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, magnification, psf, centre, spot, regionBounds, progress, increment, centreEachSlice);
    }
    if (interactiveMode) {
        Utils.hide(TITLE_INTENSITY);
    }
    IJ.showProgress(1);
    if (threadPool != null) {
        threadPool.shutdownNow();
        threadPool = null;
    }
    if (!ok || stats.getN() == 0)
        return;
    final double avSd = getAverage(averageSd, averageA, 2);
    Utils.log("  Average background = %.2f, Av. SD = %s px", stats.getMean(), Utils.rounded(avSd, 4));
    normalise(psf, maxz, avSd * magnification, false);
    IJ.showProgress(1);
    psfImp = Utils.display("PSF", psf);
    psfImp.setSlice(maxz);
    psfImp.resetDisplayRange();
    psfImp.updateAndDraw();
    double[][] fitCom = new double[2][psf.getSize()];
    Arrays.fill(fitCom[0], Double.NaN);
    Arrays.fill(fitCom[1], Double.NaN);
    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
    double[][] com = calculateCentreOfMass(psf, fitCom, nmPerPixel / magnification);
    double[] slice = Utils.newArray(psf.getSize(), 1, 1.0);
    String title = TITLE + " CoM Drift";
    Plot2 plot = new Plot2(title, "Slice", "Drift (nm)");
    plot.addLabel(0, 0, "Red = X; Blue = Y");
    //double[] limitsX = Maths.limits(com[0]);
    //double[] limitsY = Maths.limits(com[1]);
    double[] limitsX = getLimits(com[0]);
    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);
    Utils.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", XmlUtils.toXML(new PSFSettings(maxz, nmPerPixel / magnification, nmPerSlice, stats.getN(), fwhm, createNote())));
    Utils.log("%s : z-centre = %d, nm/Pixel = %s, nm/Slice = %s, %d images, PSF SD = %s nm, FWHM = %s px\n", psfImp.getTitle(), maxz, Utils.rounded(nmPerPixel / magnification, 3), Utils.rounded(nmPerSlice, 3), stats.getN(), Utils.rounded(fittedSd * nmPerPixel, 4), Utils.rounded(fwhm));
    createInteractivePlots(psf, maxz, nmPerPixel / magnification, fittedSd * nmPerPixel);
    IJ.showStatus("");
}

Example 2 with ImageExtractor

use of gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class PSFCreator method fitSpot.

private MemoryPeakResults fitSpot(ImageStack stack, final int width, final int height, final int x, final int y) {
    Rectangle regionBounds = null;
    // Create a fit engine
    MemoryPeakResults results = new MemoryPeakResults();
    results.setSortAfterEnd(true);
    results.begin();
    FitEngine engine = new FitEngine(config, results, Prefs.getThreads(), FitQueue.BLOCKING);
    List<ParameterisedFitJob> jobItems = new ArrayList<ParameterisedFitJob>(stack.getSize());
    for (int slice = 1; slice <= stack.getSize(); slice++) {
        // Extract the region from each frame
        ImageExtractor ie = new ImageExtractor((float[]) stack.getPixels(slice), width, height);
        if (regionBounds == null)
            regionBounds = ie.getBoxRegionBounds(x, y, boxRadius);
        float[] region = ie.crop(regionBounds);
        // Fit only a spot in the centre
        FitParameters params = new FitParameters();
        params.maxIndices = new int[] { boxRadius * regionBounds.width + boxRadius };
        ParameterisedFitJob job = new ParameterisedFitJob(slice, params, slice, region, regionBounds);
        jobItems.add(job);
        engine.run(job);
    }
    engine.end(false);
    results.end();
    return results;
}

Example 3 with ImageExtractor

use of gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class GaussianFit method runFinal.

/**
	 * Perform fitting using the chosen maxima. Update the overlay if successful.
	 * 
	 * @param ip
	 *            The input image
	 */
private void runFinal(ImageProcessor ip) {
    ip.reset();
    Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    float[] data = ImageConverter.getData(ip);
    int width = bounds.width;
    int height = bounds.height;
    // Sort the maxima
    float[] smoothData = data;
    if (getSmooth() > 0) {
        // Smoothing destructively modifies the data so create a copy
        smoothData = Arrays.copyOf(data, width * height);
        AverageFilter filter = new AverageFilter();
        //filter.blockAverage(smoothData, width, height, smooth);
        if (smooth <= border)
            filter.stripedBlockAverageInternal(smoothData, width, height, (float) smooth);
        else
            filter.stripedBlockAverage(smoothData, width, height, (float) smooth);
    }
    Sort.sort(maxIndices, smoothData);
    // Show the candidate peaks
    if (maxIndices.length > 0) {
        String message = String.format("Identified %d peaks", maxIndices.length);
        if (isLogProgress()) {
            IJ.log(message);
            for (int index : maxIndices) {
                IJ.log(String.format("  %.2f @ [%d,%d]", data[index], bounds.x + index % width, bounds.y + index / width));
            }
        }
        // Check whether to run if the number of peaks is large
        if (maxIndices.length > 10) {
            GenericDialog gd = new GenericDialog("Warning");
            gd.addMessage(message + "\nDo you want to fit?");
            gd.showDialog();
            if (gd.wasCanceled())
                return;
        }
    } else {
        IJ.log("No maxima identified");
        return;
    }
    results = new IJTablePeakResults(showDeviations, imp.getTitle() + " [" + imp.getCurrentSlice() + "]");
    results.begin();
    // Perform the Gaussian fit
    long ellapsed = 0;
    if (!singleFit) {
        if (isLogProgress())
            IJ.log("Combined fit");
        // Estimate height from smoothed data
        double[] estimatedHeights = new double[maxIndices.length];
        for (int i = 0; i < estimatedHeights.length; i++) estimatedHeights[i] = smoothData[maxIndices[i]];
        FitConfiguration config = new FitConfiguration();
        setupPeakFiltering(config);
        long time = System.nanoTime();
        double[] params = fitMultiple(data, width, height, maxIndices, estimatedHeights);
        ellapsed = System.nanoTime() - time;
        if (params != null) {
            // Copy all the valid parameters into a new array
            double[] validParams = new double[params.length];
            int c = 0;
            int validPeaks = 0;
            validParams[c++] = params[0];
            double[] initialParams = convertParameters(fitResult.getInitialParameters());
            double[] paramsDev = convertParameters(fitResult.getParameterStdDev());
            Rectangle regionBounds = new Rectangle();
            int[] xpoints = new int[maxIndices.length];
            int[] ypoints = new int[maxIndices.length];
            int nMaxima = 0;
            for (int i = 1, n = 0; i < params.length; i += 6, n++) {
                int y = maxIndices[n] / width;
                int x = maxIndices[n] % width;
                // Check the peak is a good fit
                if (filterResults && config.validatePeak(n, initialParams, params) != FitStatus.OK)
                    continue;
                if (showFit) {
                    // Copy the valid parameters
                    validPeaks++;
                    for (int ii = i, j = 0; j < 6; ii++, j++) validParams[c++] = params[ii];
                }
                double[] peakParams = extractParams(params, i);
                double[] peakParamsDev = extractParams(paramsDev, i);
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, nMaxima, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            }
            setOverlay(nMaxima, xpoints, ypoints);
            // Draw the fit
            if (showFit && validPeaks != 0) {
                double[] pixels = new double[data.length];
                EllipticalGaussian2DFunction f = new EllipticalGaussian2DFunction(validPeaks, width, height);
                invertParameters(validParams);
                f.initialise(validParams);
                for (int x = 0; x < pixels.length; x++) pixels[x] = f.eval(x);
                FloatProcessor fp = new FloatProcessor(width, height, pixels);
                // Insert into a full size image
                FloatProcessor fp2 = new FloatProcessor(ip.getWidth(), ip.getHeight());
                fp2.insert(fp, bounds.x, bounds.y);
                Utils.display(TITLE, fp2);
            }
        } else {
            if (isLogProgress()) {
                IJ.log("Failed to fit " + Utils.pleural(maxIndices.length, "peak") + getReason(fitResult));
            }
            imp.setOverlay(null);
        }
    } else {
        if (isLogProgress())
            IJ.log("Individual fit");
        int nMaxima = 0;
        int[] xpoints = new int[maxIndices.length];
        int[] ypoints = new int[maxIndices.length];
        // Extract each peak and fit individually
        ImageExtractor ie = new ImageExtractor(data, width, height);
        float[] region = null;
        Gaussian2DFitter gf = createGaussianFitter(filterResults);
        for (int n = 0; n < maxIndices.length; n++) {
            int y = maxIndices[n] / width;
            int x = maxIndices[n] % width;
            long time = System.nanoTime();
            Rectangle regionBounds = ie.getBoxRegionBounds(x, y, singleRegionSize);
            region = ie.crop(regionBounds, region);
            int newIndex = (y - regionBounds.y) * regionBounds.width + x - regionBounds.x;
            if (isLogProgress()) {
                IJ.log("Fitting peak " + (n + 1));
            }
            double[] peakParams = fitSingle(gf, region, regionBounds.width, regionBounds.height, newIndex, smoothData[maxIndices[n]]);
            ellapsed += System.nanoTime() - time;
            // Output fit result
            if (peakParams != null) {
                double[] peakParamsDev = null;
                if (showDeviations) {
                    peakParamsDev = convertParameters(fitResult.getParameterStdDev());
                }
                addResult(bounds, regionBounds, data, peakParams, peakParamsDev, n, x, y, data[maxIndices[n]]);
                // Add fit result to the overlay - Coords are updated with the region offsets in addResult
                double xf = peakParams[3];
                double yf = peakParams[4];
                xpoints[nMaxima] = (int) (xf + 0.5);
                ypoints[nMaxima] = (int) (yf + 0.5);
                nMaxima++;
            } else {
                if (isLogProgress()) {
                    IJ.log("Failed to fit peak " + (n + 1) + getReason(fitResult));
                }
            }
        }
        // Update the overlay
        if (nMaxima > 0)
            setOverlay(nMaxima, xpoints, ypoints);
        else
            imp.setOverlay(null);
    }
    results.end();
    if (isLogProgress())
        IJ.log("Time = " + (ellapsed / 1000000.0) + "ms");
}

Example 4 with ImageExtractor

use of gdsc.core.utils.ImageExtractor in project GDSC-SMLM by aherbert.

the class FitWorker method run.

/**
	 * Locate all the peaks in the image specified by the fit job
	 * <p>
	 * WARNING: The FitWorker fits a sub-region of the data for each maxima. It then updates the FitResult parameters
	 * with an offset reflecting the position. The initialParameters are not updated with this offset unless configured.
	 * 
	 * @param job
	 *            The fit job
	 */
public void run(FitJob job) {
    final long start = System.nanoTime();
    job.start();
    this.job = job;
    benchmarking = false;
    this.slice = job.slice;
    // Used for debugging
    //if (logger == null) logger = new gdsc.fitting.logging.ConsoleLogger();
    // Crop to the ROI
    cc = new CoordinateConverter(job.bounds);
    final int width = cc.dataBounds.width;
    final int height = cc.dataBounds.height;
    borderLimitX = width - border;
    borderLimitY = height - border;
    data = job.data;
    FitParameters params = job.getFitParameters();
    this.endT = (params != null) ? params.endT : -1;
    candidates = indentifySpots(job, width, height, params);
    if (candidates.size == 0) {
        finish(job, start);
        return;
    }
    fittedBackground = 0;
    // Always get the noise and store it with the results.
    if (params != null && !Float.isNaN(params.noise)) {
        noise = params.noise;
        fitConfig.setNoise(noise);
    } else if (calculateNoise) {
        noise = estimateNoise(width, height);
        fitConfig.setNoise(noise);
    }
    //System.out.printf("Slice %d : Noise = %g\n", slice, noise);
    if (logger != null)
        logger.info("Slice %d: Noise = %f", slice, noise);
    final ImageExtractor ie = new ImageExtractor(data, width, height);
    double[] region = null;
    if (params != null && params.fitTask == FitTask.MAXIMA_IDENITIFICATION) {
        final float sd0 = (float) fitConfig.getInitialPeakStdDev0();
        final float sd1 = (float) fitConfig.getInitialPeakStdDev1();
        for (int n = 0; n < candidates.getSize(); n++) {
            // Find the background using the perimeter of the data.
            // TODO - Perhaps the Gaussian Fitter should be used to produce the initial estimates but no actual fit done.
            // This would produce coords using the centre-of-mass.
            final Candidate candidate = candidates.get(n);
            final int x = candidate.x;
            final int y = candidate.y;
            final Rectangle regionBounds = ie.getBoxRegionBounds(x, y, fitting);
            region = ie.crop(regionBounds, region);
            final float b = (float) Gaussian2DFitter.getBackground(region, regionBounds.width, regionBounds.height, 1);
            // Offset the coords to the centre of the pixel. Note the bounds will be added later.
            // Subtract the background to get the amplitude estimate then convert to signal.
            final float amplitude = candidate.intensity - ((relativeIntensity) ? 0 : b);
            final float signal = (float) (amplitude * 2.0 * Math.PI * sd0 * sd1);
            final float[] peakParams = new float[] { b, signal, 0, x + 0.5f, y + 0.5f, sd0, sd1 };
            final int index = y * width + x;
            sliceResults.add(createResult(cc.fromDataToGlobalX(x), cc.fromDataToGlobalY(y), data[index], 0, noise, peakParams, null, n));
        }
    } else {
        initialiseFitting();
        // Smooth the data to provide initial background estimates
        final BlockAverageDataProcessor processor = new BlockAverageDataProcessor(1, 1);
        final float[] smoothedData = processor.process(data, width, height);
        final ImageExtractor ie2 = new ImageExtractor(smoothedData, width, height);
        // Allow the results to be filtered for certain peaks
        if (params != null && params.filter != null) {
            resultFilter = new DistanceResultFilter(params.filter, params.distanceThreshold, candidates.getlength());
        //filter = new OptimumDistanceResultFilter(params.filter, params.distanceThreshold, maxIndices.length);
        }
        // The SpotFitter is used to create a dynamic MultiPathFitResult object.
        // This is then passed to a multi-path filter. Thus the same fitting decision process 
        // is used when benchmarking and when running on actual data.
        // Note: The SpotFitter labels each PreprocessedFitResult using the offset in the FitResult object.
        // The initial params and deviations can then be extracted for the results that pass the filter.
        MultiPathFilter filter;
        IMultiPathFitResults multiPathResults = this;
        SelectedResultStore store = this;
        coordinateStore = coordinateStore.resize(width, height);
        if (params != null && params.fitTask == FitTask.BENCHMARKING) {
            // Run filtering as normal. However in the event that a candidate is missed or some
            // results are not generated we must generate them. This is done in the complete(int)
            // method if we set the benchmarking flag.
            benchmarking = true;
            // Filter using the benchmark filter
            filter = params.benchmarkFilter;
            if (filter == null) {
                // Create a default filter using the standard FitConfiguration to ensure sensible fits
                // are stored as the current slice results.
                // Note the current fit configuration for benchmarking may have minimal filtering settings
                // so we do not use that object.
                final FitConfiguration tmp = new FitConfiguration();
                final double residualsThreshold = 0.4;
                filter = new MultiPathFilter(tmp, createMinimalFilter(), residualsThreshold);
            }
        } else {
            // Filter using the configuration
            filter = new MultiPathFilter(fitConfig, createMinimalFilter(), config.getResidualsThreshold());
        }
        // If we are benchmarking then do not generate results dynamically since we will store all 
        // results in the fit job
        dynamicMultiPathFitResult = new DynamicMultiPathFitResult(ie, ie2, !benchmarking);
        // Debug where the fit config may be different between benchmarking and fitting
        if (slice == -1) {
            SettingsManager.saveFitEngineConfiguration(config, String.format("/tmp/config.%b.xml", benchmarking));
            Utils.write(String.format("/tmp/filter.%b.xml", benchmarking), filter.toXML());
            //filter.setDebugFile(String.format("/tmp/fitWorker.%b.txt", benchmarking));
            StringBuilder sb = new StringBuilder();
            sb.append((benchmarking) ? ((gdsc.smlm.results.filter.Filter) filter.getFilter()).toXML() : fitConfig.getSmartFilter().toXML()).append("\n");
            sb.append(((gdsc.smlm.results.filter.Filter) filter.getMinimalFilter()).toXML()).append("\n");
            sb.append(filter.residualsThreshold).append("\n");
            sb.append(config.getFailuresLimit()).append("\n");
            sb.append(fitConfig.getDuplicateDistance()).append("\n");
            if (spotFilter != null)
                sb.append(spotFilter.getDescription()).append("\n");
            sb.append("MaxCandidate = ").append(candidates.getSize()).append("\n");
            for (int i = 0; i < candidates.list.length; i++) {
                sb.append(String.format("Fit %d [%d,%d = %.1f]\n", i, candidates.get(i).x, candidates.get(i).y, candidates.get(i).intensity));
            }
            Utils.write(String.format("/tmp/candidates.%b.xml", benchmarking), sb.toString());
        }
        filter.select(multiPathResults, config.getFailuresLimit(), true, store, coordinateStore);
        if (logger != null)
            logger.info("Slice %d: %d / %d", slice, success, candidates.getSize());
        // Result filter post-processing
        if (resultFilter != null) {
            resultFilter.finalise();
            job.setResults(resultFilter.getResults());
            job.setIndices(resultFilter.getMaxIndices());
            for (int i = 0; i < resultFilter.getFilteredCount(); i++) {
                job.setFitResult(i, resultFilter.getFitResults()[i]);
            }
            sliceResults.clear();
            sliceResults.addAll(resultFilter.getResults());
        }
    }
    // Add the ROI bounds to the fitted peaks
    final float offsetx = cc.dataBounds.x;
    final float offsety = cc.dataBounds.y;
    for (int i = 0; i < sliceResults.size(); i++) {
        final PeakResult result = sliceResults.get(i);
        result.params[Gaussian2DFunction.X_POSITION] += offsetx;
        result.params[Gaussian2DFunction.Y_POSITION] += offsety;
    }
    this.results.addAll(sliceResults);
    finish(job, start);
}