Examples with PreprocessedPeakResult uk.ac.sussex.gdsc.smlm.results.filter.PreprocessedPeakResult used on opensource projects

Example 1 with PreprocessedPeakResult

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

the class FitWorker method add.

@Override
public void add(SelectedResult selectedResult) {
    // TODO - Print the current state of the dynamicMultiPathFitResult to file.
    // This will allow debugging what is different between the benchmark fit and the PeakFit.
    // Output:
    // slice
    // candidate Id
    // Initial and final params for each fit result.
    // Details of the selected result.
    // Add to the slice results.
    final PreprocessedPeakResult[] results = selectedResult.results;
    if (results == null) {
        if (logger != null) {
            final int candidateId = dynamicMultiPathFitResult.getCandidateId();
            // @formatter:off
            LoggerUtils.log(logger, Level.INFO, "Not fit %d (%d,%d) %s", candidateId, cc.fromDataToGlobalX(candidates.get(candidateId).x), cc.fromDataToGlobalY(candidates.get(candidateId).y), (selectedResult.fitResult != null && selectedResult.fitResult.data != null) ? ((FitResult) selectedResult.fitResult.data).getStatus().toString() : "");
        // @formatter:on
        }
        // Reporting
        if (this.counter != null) {
            final FitType fitType = dynamicMultiPathFitResult.fitType;
            addFitType(fitType);
        }
        return;
    }
    final int candidateId = dynamicMultiPathFitResult.getCandidateId();
    final FitResult fitResult = (FitResult) selectedResult.fitResult.data;
    // The background for each result was the local background. We want the fitted global background
    final float background = (float) fitResult.getParameters()[0];
    final double[] dev = fitResult.getParameterDeviations();
    if (queueSize != 0) {
        throw new RuntimeException("There are results queued already!");
    }
    for (int i = 0; i < results.length; i++) {
        final PreprocessedPeakResult peak = results[i];
        if (peak.isExistingResult()) {
            continue;
        }
        if (peak.isNewResult()) {
            final double[] p = peak.toGaussian2DParameters();
            // Store slice results relative to the data frame (not the global bounds)
            // Convert back so that 0,0 is the top left of the data bounds
            p[Gaussian2DFunction.X_POSITION] -= cc.dataBounds.x;
            p[Gaussian2DFunction.Y_POSITION] -= cc.dataBounds.y;
            final float[] params = new float[p.length];
            params[Gaussian2DFunction.BACKGROUND] = background;
            for (int j = 1; j < p.length; j++) {
                params[j] = (float) p[j];
            }
            final float[] paramDevs;
            if (dev == null) {
                paramDevs = null;
            } else {
                paramDevs = new float[p.length];
                paramDevs[Gaussian2DFunction.BACKGROUND] = (float) dev[Gaussian2DFunction.BACKGROUND];
                final int offset = peak.getId() * Gaussian2DFunction.PARAMETERS_PER_PEAK;
                for (int j = 1; j < p.length; j++) {
                    paramDevs[j] = (float) dev[offset + j];
                }
            }
            addSingleResult(peak.getCandidateId(), params, paramDevs, fitResult.getError(), peak.getNoise(), fitConfig.getLocationVariance(peak));
            if (logger != null) {
                // Show the shift, signal and width spread
                LoggerUtils.log(logger, Level.INFO, "Fit OK %d (%.1f,%.1f) [%d]: Shift = %.3f,%.3f : SNR = %.2f : Width = %.2f,%.2f", peak.getCandidateId(), peak.getX(), peak.getY(), peak.getId(), Math.sqrt(peak.getXRelativeShift2()), Math.sqrt(peak.getYRelativeShift2()), peak.getSnr(), peak.getXSdFactor(), peak.getYSdFactor());
            }
        } else {
        // This is a candidate that passed validation. Store the estimate as passing the primary
        // filter.
        // We now do this is the pass() method.
        // storeEstimate(results[i].getCandidateId(), results[i], FILTER_RANK_PRIMARY);
        }
    }
    job.setFitResult(candidateId, fitResult);
    // Reporting
    if (this.counter != null) {
        final FitType fitType = dynamicMultiPathFitResult.fitType;
        if (selectedResult.fitResult.getStatus() == 0) {
            fitType.setOk(true);
            if (dynamicMultiPathFitResult.getSuperMultiFitResult() == selectedResult.fitResult) {
                fitType.setMultiOk(true);
            } else if (dynamicMultiPathFitResult.getSuperMultiDoubletFitResult() == selectedResult.fitResult) {
                fitType.setMultiDoubletOk(true);
            } else if (dynamicMultiPathFitResult.getSuperDoubletFitResult() == selectedResult.fitResult) {
                fitType.setDoubletOk(true);
            }
        }
        addFitType(fitType);
    }
    if (logger != null) {
        switch(fitResult.getStatus()) {
            case OK:
                // We log good results in the loop above.
                break;
            case BAD_PARAMETERS:
            case FAILED_TO_ESTIMATE_WIDTH:
                logger.log(Level.INFO, () -> "Bad parameters: " + Arrays.toString(fitResult.getInitialParameters()));
                break;
            default:
                logger.log(Level.INFO, "{0}", fitResult.getStatus());
                break;
        }
    }
    // Debugging
    if (debugLogger != null) {
        double[] peakParams = fitResult.getParameters();
        if (peakParams != null) {
            // Parameters are the raw values from fitting the region. Convert for logging.
            peakParams = Arrays.copyOf(peakParams, peakParams.length);
            final int npeaks = peakParams.length / Gaussian2DFunction.PARAMETERS_PER_PEAK;
            for (int i = 0; i < npeaks; i++) {
                peakParams[i * Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.X_POSITION] += cc.fromFitRegionToGlobalX();
                peakParams[i * Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.Y_POSITION] += cc.fromFitRegionToGlobalY();
                if (fitConfig.isAngleFitting()) {
                    peakParams[i * Gaussian2DFunction.PARAMETERS_PER_PEAK + Gaussian2DFunction.ANGLE] *= 180.0 / Math.PI;
                }
            }
        }
        final int x = candidates.get(candidateId).x;
        final int y = candidates.get(candidateId).y;
        LoggerUtils.log(debugLogger, Level.INFO, "%d:%d [%d,%d] %s (%s) = %s", slice, candidateId, cc.fromDataToGlobalX(x), cc.fromDataToGlobalY(y), fitResult.getStatus(), fitResult.getStatusData(), Arrays.toString(peakParams));
    }
}

Example 2 with PreprocessedPeakResult

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

the class FitConfiguration method validatePeak.

/**
 * Check peak to see if the fit was sensible.
 *
 * @param n The peak number
 * @param initialParams The initial peak parameters
 * @param params The fitted peak parameters
 * @param paramDevs the fitted peak parameter variances (can be null)
 * @return True if the fit fails the criteria
 */
public FitStatus validatePeak(int n, double[] initialParams, double[] params, double[] paramDevs) {
    if (isDirectFilter()) {
        // Use the global noise.
        if (peakResultValidationData != null) {
            peakResultValidationData.setResult(n, initialParams, params, paramDevs);
        }
        final PreprocessedPeakResult peak = createPreprocessedPeakResult(0, n, initialParams, params, paramDevs, peakResultValidationData, ResultType.NEW, 0, 0, false);
        if (directFilter.accept(peak)) {
            return setValidationResult(FitStatus.OK, null);
        }
        if (log != null) {
            log.info(() -> String.format("Bad peak %d: %s", peak.getId(), DirectFilter.getStatusMessage(peak, directFilter.getResult())));
        }
        if (DirectFilter.anySet(directFilter.getResult(), FilterValidationFlag.X_SD_FACTOR | FilterValidationFlag.Y_SD_FACTOR)) {
            return setValidationResult(FitStatus.WIDTH_DIVERGED, null);
        }
        // At the moment we do not get any other validation data
        return setValidationResult(FitStatus.FAILED_SMART_FILTER, null);
    }
    // where additional peaks will be neighbours. In the future we may want to control this better.
    if (isRegionValidation()) {
        final int offset = n * Gaussian2DFunction.PARAMETERS_PER_PEAK;
        final double x = params[Gaussian2DFunction.X_POSITION + offset] + coordinateOffset;
        final double y = params[Gaussian2DFunction.Y_POSITION + offset] + coordinateOffset;
        if (x <= 0 || x >= fitRegionWidth || y <= 0 || y >= fitRegionHeight) {
            if (log != null) {
                log.info(() -> String.format("Bad peak %d: Coordinates outside fit region (x=%g,y=%g) <> %d,%d", n, x, y, fitRegionWidth, fitRegionHeight));
            }
            return setValidationResult(FitStatus.OUTSIDE_FIT_REGION, new double[] { x, y, fitRegionWidth, fitRegionHeight });
        }
    }
    if (isDisableSimpleFilter()) {
        return setValidationResult(FitStatus.OK, null);
    }
    final int offset = n * Gaussian2DFunction.PARAMETERS_PER_PEAK;
    // Check spot movement
    final double xShift = params[Gaussian2DFunction.X_POSITION + offset] - initialParams[Gaussian2DFunction.X_POSITION + offset];
    final double yShift = params[Gaussian2DFunction.Y_POSITION + offset] - initialParams[Gaussian2DFunction.Y_POSITION + offset];
    final double maxShift = coordinateShift;
    if (Math.abs(xShift) > maxShift || Math.abs(yShift) > maxShift) {
        if (log != null) {
            log.info(() -> String.format("Bad peak %d: Fitted coordinates moved (x=%g,y=%g) > %g", n, xShift, yShift, maxShift));
        }
        return setValidationResult(FitStatus.COORDINATES_MOVED, new double[] { xShift, yShift });
    }
    if (zEnabled) {
        final double z = params[Gaussian2DFunction.Z_POSITION + offset];
        if (z < getMinZ() || z > getMaxZ()) {
            return setValidationResult(FitStatus.Z_MOVED, z);
        }
    }
    // Check signal threshold.
    // The threshold should be set in the same units as those used during fitting.
    final double signal = params[Gaussian2DFunction.SIGNAL + offset];
    // Compare the signal to the desired signal strength
    if (signal < minSignal) {
        if (log != null) {
            log.info(() -> String.format("Bad peak %d: Insufficient signal %g", n, signal));
        }
        // System.out.printf("Bad peak %d: Insufficient signal (%g)\n", n, signal);
        return setValidationResult(FitStatus.INSUFFICIENT_SIGNAL, signal);
    }
    double xsd;
    double ysd;
    // Map the width parameters using the z-model if present
    if (getAstigmatismZModel() != null) {
        final double z = params[Gaussian2DFunction.Z_POSITION + offset];
        xsd = astigmatismZModel.getSx(z);
        ysd = astigmatismZModel.getSy(z);
        // Check widths. This may be the only filter used even if z-fitting
        // (i.e. a z-depth filter is not used)
        final double xFactor = xsd / initialParams[Gaussian2DFunction.X_SD + offset];
        final double yFactor = ysd / initialParams[Gaussian2DFunction.Y_SD + offset];
        final double s2 = xFactor * yFactor;
        if (s2 > widthFactor || s2 < minWidthFactor) {
            if (log != null) {
                log.info(() -> String.format("Bad peak %d: Fitted width diverged (x=%gx,y=%gx)", n, xFactor, yFactor));
            }
            return setValidationResult(FitStatus.WIDTH_DIVERGED, new double[] { xFactor, yFactor });
        }
    } else {
        xsd = params[Gaussian2DFunction.X_SD + offset];
        ysd = params[Gaussian2DFunction.Y_SD + offset];
    }
    double noise = this.noise;
    if (peakResultValidationData != null) {
        peakResultValidationData.setResult(n, initialParams, params, paramDevs);
        noise = peakResultValidationData.getNoise();
    }
    // Check SNR threshold. Assume the noise has been set in the same units as the fit result.
    final double snr = Gaussian2DPeakResultHelper.getMeanSignalUsingP05(signal, xsd, ysd) / noise;
    // Compare the signal to the desired signal strength
    if (snr < getSignalStrength()) {
        if (log != null) {
            log.info(() -> String.format("Bad peak %d: Insufficient SNR %g", n, snr));
        }
        // System.out.printf("Bad peak %d: Insufficient SNR (%g)\n", n, snr);
        return setValidationResult(FitStatus.INSUFFICIENT_SNR, snr);
    }
    // Check widths
    if (isXSdFitting()) {
        boolean badWidth = false;
        double xfactor = 0;
        double yfactor = 0;
        xfactor = xsd / initialParams[Gaussian2DFunction.X_SD + offset];
        if (isTwoAxisGaussian2D) {
            yfactor = ysd / initialParams[Gaussian2DFunction.Y_SD + offset];
            final double s2 = xfactor * yfactor;
            badWidth = (s2 > widthFactor || s2 < minWidthFactor);
        } else {
            badWidth = (xfactor > widthFactor || xfactor < minWidthFactor);
            yfactor = xfactor;
        }
        if (badWidth) {
            if (log != null) {
                final double localXFactor = xfactor;
                final double localYFactor = yfactor;
                log.info(() -> String.format("Bad peak %d: Fitted width diverged (x=%gx,y=%gx)", n, localXFactor, localYFactor));
            }
            return setValidationResult(FitStatus.WIDTH_DIVERGED, new double[] { xfactor, yfactor });
        }
    }
    // Check precision. This is above zero if a threshold is present.
    if (precisionThreshold > 0) {
        final double variance;
        switch(getPrecisionMethodValue()) {
            case PrecisionMethod.MORTENSEN_VALUE:
            case PrecisionMethod.MORTENSEN_LOCAL_BACKGROUND_VALUE:
                final double sd = (isTwoAxisGaussian2D) ? Gaussian2DPeakResultHelper.getStandardDeviation(xsd, ysd) : xsd;
                final double localBackground = isPrecisionUsingBackground() && peakResultValidationData != null ? peakResultValidationData.getLocalBackground() : params[Gaussian2DFunction.BACKGROUND];
                variance = getVariance(localBackground, params[Gaussian2DFunction.SIGNAL + offset] * signalToPhotons, sd, isPrecisionUsingBackground());
                break;
            case PrecisionMethod.POISSON_CRLB_VALUE:
                variance = getVariance(paramDevs, n);
                break;
            default:
                // This should not happen
                throw new IllegalStateException("Unknown precision method: " + getPrecisionMethod());
        }
        if (variance > precisionThreshold) {
            if (log != null) {
                final double precision = Math.sqrt(variance);
                log.info(() -> String.format("Bad peak %d: Insufficient precision (%gx)", n, precision));
            }
            return setValidationResult(FitStatus.INSUFFICIENT_PRECISION, variance);
        }
    }
    return setValidationResult(FitStatus.OK, null);
}

Example 3 with PreprocessedPeakResult

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

the class BenchmarkFilterAnalysis method createResults.

/**
 * Create peak results.
 *
 * @param filterResults The results from running the filter (or null)
 * @param filter the filter
 */
private MemoryPeakResults createResults(PreprocessedPeakResult[] filterResults, DirectFilter filter, boolean withBorder) {
    if (filterResults == null) {
        final MultiPathFilter multiPathFilter = createMpf(filter, defaultMinimalFilter);
        filterResults = filterResults(multiPathFilter);
    }
    final MemoryPeakResults newResults = new MemoryPeakResults();
    newResults.copySettings(this.results);
    newResults.setName(TITLE);
    if (withBorder) {
        // To produce the same results as the PeakFit plugin we must implement the border
        // functionality used in the FitWorker. This respects the border of the spot filter.
        final FitEngineConfiguration config = new FitEngineConfiguration();
        updateAllConfiguration(config);
        final MaximaSpotFilter spotFilter = config.createSpotFilter();
        final int border = spotFilter.getBorder();
        final Rectangle bounds = getBounds();
        final int borderLimitX = bounds.x + bounds.width - border;
        final int borderLimitY = bounds.y + bounds.height - border;
        for (final PreprocessedPeakResult spot : filterResults) {
            if (spot.getX() > border && spot.getX() < borderLimitX && spot.getY() > border && spot.getY() < borderLimitY) {
                final double[] p = spot.toGaussian2DParameters();
                final float[] params = new float[p.length];
                for (int j = 0; j < p.length; j++) {
                    params[j] = (float) p[j];
                }
                final int frame = spot.getFrame();
                final int origX = (int) p[Gaussian2DFunction.X_POSITION];
                final int origY = (int) p[Gaussian2DFunction.Y_POSITION];
                newResults.add(frame, origX, origY, 0, 0, spot.getNoise(), spot.getMeanSignal(), params, null);
            }
        }
    } else {
        for (final PreprocessedPeakResult spot : filterResults) {
            final double[] p = spot.toGaussian2DParameters();
            final float[] params = new float[p.length];
            for (int j = 0; j < p.length; j++) {
                params[j] = (float) p[j];
            }
            final int frame = spot.getFrame();
            final int origX = (int) p[Gaussian2DFunction.X_POSITION];
            final int origY = (int) p[Gaussian2DFunction.Y_POSITION];
            newResults.add(frame, origX, origY, 0, 0, spot.getNoise(), spot.getMeanSignal(), params, null);
        }
    }
    return newResults;
}

Example 4 with PreprocessedPeakResult

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

the class BenchmarkFilterAnalysis method showOverlay.

/**
 * Show overlay.
 *
 * <ul>
 *
 * <li>Green = TP
 *
 * <li>Red = FP
 *
 * <li>Magenta = FP (Ignored from analysis)
 *
 * <li>Yellow = FN
 *
 * <li>Orange = FN (Outside border)
 *
 * </ul>
 *
 * @param allAssignments The assignments generated from running the filter (or null)
 * @param filter the filter
 * @return The results from running the filter (or null)
 */
@Nullable
@SuppressWarnings("null")
private PreprocessedPeakResult[] showOverlay(ArrayList<FractionalAssignment[]> allAssignments, DirectFilter filter) {
    final ImagePlus imp = CreateData.getImage();
    if (imp == null) {
        return null;
    }
    // Run the filter manually to get the results that pass.
    if (allAssignments == null) {
        allAssignments = getAssignments(filter);
    }
    final Overlay o = new Overlay();
    // Do TP
    final TIntHashSet actual = new TIntHashSet();
    final TIntHashSet predicted = new TIntHashSet();
    for (final FractionalAssignment[] assignments : allAssignments) {
        if (assignments == null || assignments.length == 0) {
            continue;
        }
        float[] tx = null;
        float[] ty = null;
        int count = 0;
        if (settings.showTP) {
            tx = new float[assignments.length];
            ty = new float[assignments.length];
        }
        int frame = 0;
        for (int i = 0; i < assignments.length; i++) {
            final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
            final UniqueIdPeakResult peak = (UniqueIdPeakResult) c.peak;
            final BasePreprocessedPeakResult spot = (BasePreprocessedPeakResult) c.peakResult;
            actual.add(peak.uniqueId);
            predicted.add(spot.getUniqueId());
            frame = spot.getFrame();
            if (settings.showTP) {
                tx[count] = spot.getX();
                ty[count++] = spot.getY();
            }
        }
        if (settings.showTP) {
            SpotFinderPreview.addRoi(frame, o, tx, ty, count, Color.green);
        }
    }
    float[] x = new float[10];
    float[] y = new float[x.length];
    float[] x2 = new float[10];
    float[] y2 = new float[x2.length];
    // Do FP (all remaining results that are not a TP)
    PreprocessedPeakResult[] filterResults = null;
    if (settings.showFP) {
        final MultiPathFilter multiPathFilter = createMpf(filter, defaultMinimalFilter);
        filterResults = filterResults(multiPathFilter);
        int frame = 0;
        int c1 = 0;
        int c2 = 0;
        for (int i = 0; i < filterResults.length; i++) {
            if (frame != filterResults[i].getFrame()) {
                if (c1 != 0) {
                    SpotFinderPreview.addRoi(frame, o, x, y, c1, Color.red);
                }
                if (c2 != 0) {
                    SpotFinderPreview.addRoi(frame, o, x2, y2, c2, Color.magenta);
                }
                c1 = c2 = 0;
            }
            frame = filterResults[i].getFrame();
            if (predicted.contains(filterResults[i].getUniqueId())) {
                continue;
            }
            if (filterResults[i].ignore()) {
                if (x2.length == c2) {
                    x2 = Arrays.copyOf(x2, c2 * 2);
                    y2 = Arrays.copyOf(y2, c2 * 2);
                }
                x2[c2] = filterResults[i].getX();
                y2[c2++] = filterResults[i].getY();
            } else {
                if (x.length == c1) {
                    x = Arrays.copyOf(x, c1 * 2);
                    y = Arrays.copyOf(y, c1 * 2);
                }
                x[c1] = filterResults[i].getX();
                y[c1++] = filterResults[i].getY();
            }
        }
        if (c1 != 0) {
            SpotFinderPreview.addRoi(frame, o, x, y, c1, Color.red);
        }
        if (c2 != 0) {
            SpotFinderPreview.addRoi(frame, o, x2, y2, c2, Color.magenta);
        }
    }
    // Do FN (all remaining peaks that have not been matched)
    if (settings.showFN) {
        final boolean checkBorder = (filterResult.analysisBorder != null && filterResult.analysisBorder.x != 0);
        final float border;
        final float xlimit;
        final float ylimit;
        if (checkBorder) {
            final Rectangle lastAnalysisBorder = filterResult.analysisBorder;
            border = lastAnalysisBorder.x;
            xlimit = lastAnalysisBorder.x + lastAnalysisBorder.width;
            ylimit = lastAnalysisBorder.y + lastAnalysisBorder.height;
        } else {
            border = xlimit = ylimit = 0;
        }
        // Add the results to the lists
        actualCoordinates.forEachEntry(new CustomTIntObjectProcedure(x, y, x2, y2) {

            @Override
            public boolean execute(int frame, UniqueIdPeakResult[] results) {
                int c1 = 0;
                int c2 = 0;
                if (x.length <= results.length) {
                    x = new float[results.length];
                    y = new float[results.length];
                }
                if (x2.length <= results.length) {
                    x2 = new float[results.length];
                    y2 = new float[results.length];
                }
                for (int i = 0; i < results.length; i++) {
                    // Ignore those that were matched by TP
                    if (actual.contains(results[i].uniqueId)) {
                        continue;
                    }
                    if (checkBorder && outsideBorder(results[i], border, xlimit, ylimit)) {
                        x2[c2] = results[i].getXPosition();
                        y2[c2++] = results[i].getYPosition();
                    } else {
                        x[c1] = results[i].getXPosition();
                        y[c1++] = results[i].getYPosition();
                    }
                }
                if (c1 != 0) {
                    SpotFinderPreview.addRoi(frame, o, x, y, c1, Color.yellow);
                }
                if (c2 != 0) {
                    SpotFinderPreview.addRoi(frame, o, x2, y2, c2, Color.orange);
                }
                return true;
            }
        });
    }
    imp.setOverlay(o);
    return filterResults;
}

Example 5 with PreprocessedPeakResult

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

the class BenchmarkSpotFit method runFitting.

private BenchmarkSpotFitResult runFitting() {
    // Extract all the results in memory into a list per frame. This can be cached
    boolean refresh = false;
    Pair<Integer, TIntObjectHashMap<List<Coordinate>>> coords = coordinateCache.get();
    if (coords.getKey() != simulationParameters.id) {
        // Do not get integer coordinates
        // The Coordinate objects will be PeakResultPoint objects that store the original PeakResult
        // from the MemoryPeakResults
        coords = Pair.of(simulationParameters.id, ResultsMatchCalculator.getCoordinates(results, false));
        coordinateCache.set(coords);
        refresh = true;
    }
    final TIntObjectHashMap<List<Coordinate>> actualCoordinates = coords.getValue();
    // Extract all the candidates into a list per frame. This can be cached if the settings have not
    // changed
    final int width = (config.isIncludeNeighbours()) ? config.getFittingWidth() : 0;
    CandidateData candidateData = candidateDataCache.get();
    if (refresh || candidateData == null || candidateData.differentSettings(filterResult.id, settings, width)) {
        candidateData = subsetFilterResults(filterResult.filterResults, width);
        candidateDataCache.set(candidateData);
    }
    final StopWatch stopWatch = StopWatch.createStarted();
    final ImageStack stack = imp.getImageStack();
    clearFitResults();
    // Save results to memory
    final MemoryPeakResults peakResults = new MemoryPeakResults();
    peakResults.copySettings(this.results);
    peakResults.setName(TITLE);
    config.configureOutputUnits();
    final FitConfiguration fitConfig = config.getFitConfiguration();
    peakResults.setCalibration(fitConfig.getCalibration());
    MemoryPeakResults.addResults(peakResults);
    // Create a pool of workers
    final int nThreads = Prefs.getThreads();
    final BlockingQueue<Integer> jobs = new ArrayBlockingQueue<>(nThreads * 2);
    final List<Worker> workers = new LinkedList<>();
    final List<Thread> threads = new LinkedList<>();
    final Ticker ticker = ImageJUtils.createTicker(stack.getSize(), nThreads, "Fitting frames ...");
    final PeakResults syncResults = SynchronizedPeakResults.create(peakResults, nThreads);
    for (int i = 0; i < nThreads; i++) {
        final Worker worker = new Worker(jobs, stack, actualCoordinates, candidateData.filterCandidates, syncResults, ticker);
        final Thread t = new Thread(worker);
        workers.add(worker);
        threads.add(t);
        t.start();
    }
    // Fit the frames
    final long startTime = System.nanoTime();
    for (int i = 1; i <= stack.getSize(); i++) {
        put(jobs, i);
    }
    // Finish all the worker threads by passing in a null job
    for (int i = 0; i < threads.size(); i++) {
        put(jobs, -1);
    }
    // Wait for all to finish
    for (int i = 0; i < threads.size(); i++) {
        try {
            threads.get(i).join();
        } catch (final InterruptedException ex) {
            Thread.currentThread().interrupt();
            throw new ConcurrentRuntimeException(ex);
        }
    }
    final long runTime = System.nanoTime() - startTime;
    threads.clear();
    ImageJUtils.finished();
    if (ImageJUtils.isInterrupted()) {
        return null;
    }
    stopWatch.stop();
    final String timeString = stopWatch.toString();
    IJ.log("Spot fit time : " + timeString);
    IJ.showStatus("Collecting results ...");
    if (fitConfig.isFitCameraCounts()) {
        // Convert to photons for consistency
        results.convertToPreferredUnits();
    }
    final TIntObjectHashMap<FilterCandidates> fitResults = new TIntObjectHashMap<>();
    for (final Worker w : workers) {
        fitResults.putAll(w.results);
    }
    // Assign a unique ID to each result
    int count = 0;
    // Materialise into an array since we use it twice
    final FilterCandidates[] candidates = fitResults.values(new FilterCandidates[fitResults.size()]);
    for (final FilterCandidates result : candidates) {
        for (int i = 0; i < result.fitResult.length; i++) {
            final MultiPathFitResult fitResult = result.fitResult[i];
            count += count(fitResult.getSingleFitResult());
            count += count(fitResult.getMultiFitResult());
            count += count(fitResult.getDoubletFitResult());
            count += count(fitResult.getMultiDoubletFitResult());
        }
    }
    final PreprocessedPeakResult[] preprocessedPeakResults = new PreprocessedPeakResult[count];
    count = 0;
    for (final FilterCandidates result : candidates) {
        for (int i = 0; i < result.fitResult.length; i++) {
            final MultiPathFitResult fitResult = result.fitResult[i];
            count = store(fitResult.getSingleFitResult(), count, preprocessedPeakResults);
            count = store(fitResult.getMultiFitResult(), count, preprocessedPeakResults);
            count = store(fitResult.getDoubletFitResult(), count, preprocessedPeakResults);
            count = store(fitResult.getMultiDoubletFitResult(), count, preprocessedPeakResults);
        }
    }
    final BenchmarkSpotFitResult newSpotFitResults = new BenchmarkSpotFitResult(simulationParameters.id, fitResults);
    newSpotFitResults.distanceInPixels = distanceInPixels;
    newSpotFitResults.lowerDistanceInPixels = lowerDistanceInPixels;
    newSpotFitResults.stopWatch = stopWatch;
    summariseResults(newSpotFitResults, runTime, preprocessedPeakResults, count, candidateData, actualCoordinates);
    IJ.showStatus("");
    spotFitResults.set(newSpotFitResults);
    return newSpotFitResults;
}