Examples with RampedScore uk.ac.sussex.gdsc.core.utils.RampedScore used on opensource projects

Example 1 with RampedScore

use of uk.ac.sussex.gdsc.core.utils.RampedScore in project GDSC-SMLM by aherbert.

the class SpotFinderPreview method run.

private void run(ImageProcessor ip, MaximaSpotFilter filter) {
    if (refreshing) {
        return;
    }
    currentSlice = imp.getCurrentSlice();
    final Rectangle bounds = ip.getRoi();
    // Crop to the ROI
    FloatProcessor fp = ip.crop().toFloat(0, null);
    float[] data = (float[]) fp.getPixels();
    final int width = fp.getWidth();
    final int height = fp.getHeight();
    // Store the mean bias and gain of the region data.
    // This is used to correctly overlay the filtered data on the original image.
    double bias = 0;
    double gain = 1;
    boolean adjust = false;
    // Set weights
    final CameraModel cameraModel = fitConfig.getCameraModel();
    if (!(cameraModel instanceof FakePerPixelCameraModel)) {
        // This should be done on the normalised data
        final float[] w = cameraModel.getNormalisedWeights(bounds);
        filter.setWeights(w, width, height);
        data = data.clone();
        if (data.length < ip.getPixelCount()) {
            adjust = true;
            bias = MathUtils.sum(cameraModel.getBias(bounds)) / data.length;
            gain = MathUtils.sum(cameraModel.getGain(bounds)) / data.length;
        }
        cameraModel.removeBiasAndGain(bounds, data);
    }
    final Spot[] spots = filter.rank(data, width, height);
    data = filter.getPreprocessedData();
    final int size = spots.length;
    if (topNScrollBar != null) {
        topNScrollBar.setMaximum(size);
        selectScrollBar.setMaximum(size);
    }
    fp = new FloatProcessor(width, height, data);
    final FloatProcessor out = new FloatProcessor(ip.getWidth(), ip.getHeight());
    out.copyBits(ip, 0, 0, Blitter.COPY);
    if (adjust) {
        fp.multiply(gain);
        fp.add(bias);
    }
    out.insert(fp, bounds.x, bounds.y);
    final double min = fp.getMin();
    final double max = fp.getMax();
    out.setMinAndMax(min, max);
    final Overlay o = new Overlay();
    o.add(new ImageRoi(0, 0, out));
    if (label != null) {
        // Get results for frame
        final Coordinate[] actual = ResultsMatchCalculator.getCoordinates(actualCoordinates, imp.getCurrentSlice());
        final Coordinate[] predicted = new Coordinate[size];
        for (int i = 0; i < size; i++) {
            predicted[i] = new BasePoint(spots[i].x + bounds.x, spots[i].y + bounds.y);
        }
        // Compute assignments
        final LocalList<FractionalAssignment> fractionalAssignments = new LocalList<>(3 * predicted.length);
        final double matchDistance = settings.distance * fitConfig.getInitialPeakStdDev();
        final RampedScore score = RampedScore.of(matchDistance, matchDistance * settings.lowerDistance / 100, false);
        final double dmin = matchDistance * matchDistance;
        final int nActual = actual.length;
        final int nPredicted = predicted.length;
        for (int j = 0; j < nPredicted; j++) {
            // Centre in the middle of the pixel
            final float x = predicted[j].getX() + 0.5f;
            final float y = predicted[j].getY() + 0.5f;
            // Any spots that match
            for (int i = 0; i < nActual; i++) {
                final double dx = (x - actual[i].getX());
                final double dy = (y - actual[i].getY());
                final double d2 = dx * dx + dy * dy;
                if (d2 <= dmin) {
                    final double d = Math.sqrt(d2);
                    final double s = score.score(d);
                    if (s == 0) {
                        continue;
                    }
                    double distance = 1 - s;
                    if (distance == 0) {
                        // In the case of a match below the distance thresholds
                        // the distance will be 0. To distinguish between candidates all below
                        // the thresholds just take the closest.
                        // We know d2 is below dmin so we subtract the delta.
                        distance -= (dmin - d2);
                    }
                    // Store the match
                    fractionalAssignments.add(new ImmutableFractionalAssignment(i, j, distance, s));
                }
            }
        }
        final FractionalAssignment[] assignments = fractionalAssignments.toArray(new FractionalAssignment[0]);
        // Compute matches
        final RankedScoreCalculator calc = RankedScoreCalculator.create(assignments, nActual - 1, nPredicted - 1);
        final boolean save = settings.showTP || settings.showFP;
        final double[] calcScore = calc.score(nPredicted, settings.multipleMatches, save);
        final ClassificationResult result = RankedScoreCalculator.toClassificationResult(calcScore, nActual);
        // Compute AUC and max jaccard (and plot)
        final double[][] curve = RankedScoreCalculator.getPrecisionRecallCurve(assignments, nActual, nPredicted);
        final double[] precision = curve[0];
        final double[] recall = curve[1];
        final double[] jaccard = curve[2];
        final double auc = AucCalculator.auc(precision, recall);
        // Show scores
        final String scoreLabel = String.format("Slice=%d, AUC=%s, R=%s, Max J=%s", imp.getCurrentSlice(), MathUtils.rounded(auc), MathUtils.rounded(result.getRecall()), MathUtils.rounded(MathUtils.maxDefault(0, jaccard)));
        setLabel(scoreLabel);
        // Plot
        String title = TITLE + " Performance";
        Plot plot = new Plot(title, "Spot Rank", "");
        final double[] rank = SimpleArrayUtils.newArray(precision.length, 0, 1.0);
        plot.setLimits(0, nPredicted, 0, 1.05);
        plot.setColor(Color.blue);
        plot.addPoints(rank, precision, Plot.LINE);
        plot.setColor(Color.red);
        plot.addPoints(rank, recall, Plot.LINE);
        plot.setColor(Color.black);
        plot.addPoints(rank, jaccard, Plot.LINE);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, scoreLabel);
        final WindowOrganiser windowOrganiser = new WindowOrganiser();
        ImageJUtils.display(title, plot, 0, windowOrganiser);
        title = TITLE + " Precision-Recall";
        plot = new Plot(title, "Recall", "Precision");
        plot.setLimits(0, 1, 0, 1.05);
        plot.setColor(Color.red);
        plot.addPoints(recall, precision, Plot.LINE);
        plot.drawLine(recall[recall.length - 1], precision[recall.length - 1], recall[recall.length - 1], 0);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, scoreLabel);
        ImageJUtils.display(title, plot, 0, windowOrganiser);
        windowOrganiser.tile();
        // Create Rois for TP and FP
        if (save) {
            final double[] matchScore = RankedScoreCalculator.getMatchScore(calc.getScoredAssignments(), nPredicted);
            int matches = 0;
            for (int i = 0; i < matchScore.length; i++) {
                if (matchScore[i] != 0) {
                    matches++;
                }
            }
            if (settings.showTP) {
                final float[] x = new float[matches];
                final float[] y = new float[x.length];
                int count = 0;
                for (int i = 0; i < matchScore.length; i++) {
                    if (matchScore[i] != 0) {
                        final BasePoint p = (BasePoint) predicted[i];
                        x[count] = p.getX() + 0.5f;
                        y[count] = p.getY() + 0.5f;
                        count++;
                    }
                }
                addRoi(0, o, x, y, count, Color.green);
            }
            if (settings.showFP) {
                final float[] x = new float[nPredicted - matches];
                final float[] y = new float[x.length];
                int count = 0;
                for (int i = 0; i < matchScore.length; i++) {
                    if (matchScore[i] == 0) {
                        final BasePoint p = (BasePoint) predicted[i];
                        x[count] = p.getX() + 0.5f;
                        y[count] = p.getY() + 0.5f;
                        count++;
                    }
                }
                addRoi(0, o, x, y, count, Color.red);
            }
        }
    } else {
        final WindowOrganiser wo = new WindowOrganiser();
        // Option to show the number of neighbours within a set pixel box radius
        final int[] count = spotFilterHelper.countNeighbours(spots, width, height, settings.neighbourRadius);
        // Show as histogram the totals...
        new HistogramPlotBuilder(TITLE, StoredData.create(count), "Neighbours").setIntegerBins(true).setPlotLabel("Radius = " + settings.neighbourRadius).show(wo);
        // TODO - Draw n=0, n=1 on the image overlay
        final LUT lut = LutHelper.createLut(LutColour.FIRE_LIGHT);
        // These are copied by the ROI
        final float[] x = new float[1];
        final float[] y = new float[1];
        // Plot the intensity
        final double[] intensity = new double[size];
        final double[] rank = SimpleArrayUtils.newArray(size, 1, 1.0);
        final int top = (settings.topN > 0) ? settings.topN : size;
        final int size_1 = size - 1;
        for (int i = 0; i < size; i++) {
            intensity[i] = spots[i].intensity;
            if (i < top) {
                x[0] = spots[i].x + bounds.x + 0.5f;
                y[0] = spots[i].y + bounds.y + 0.5f;
                final Color c = LutHelper.getColour(lut, size_1 - i, size);
                addRoi(0, o, x, y, 1, c, 2, 1);
            }
        }
        final String title = TITLE + " Intensity";
        final Plot plot = new Plot(title, "Rank", "Intensity");
        plot.setColor(Color.blue);
        plot.addPoints(rank, intensity, Plot.LINE);
        if (settings.topN > 0 && settings.topN < size) {
            plot.setColor(Color.magenta);
            plot.drawLine(settings.topN, 0, settings.topN, intensity[settings.topN - 1]);
        }
        if (settings.select > 0 && settings.select < size) {
            plot.setColor(Color.yellow);
            final int index = settings.select - 1;
            final double in = intensity[index];
            plot.drawLine(settings.select, 0, settings.select, in);
            x[0] = spots[index].x + bounds.x + 0.5f;
            y[0] = spots[index].y + bounds.y + 0.5f;
            final Color c = LutHelper.getColour(lut, size_1 - settings.select, size);
            addRoi(0, o, x, y, 1, c, 3, 3);
            plot.setColor(Color.black);
            plot.addLabel(0, 0, "Selected spot intensity = " + MathUtils.rounded(in));
        }
        ImageJUtils.display(title, plot, 0, wo);
        wo.tile();
    }
    imp.setOverlay(o);
}

Example 2 with RampedScore

use of uk.ac.sussex.gdsc.core.utils.RampedScore in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method readResults.

private MultiPathFitResults[] readResults() {
    // Extract all the results in memory into a list per frame. This can be cached
    boolean update = false;
    Pair<Integer, TIntObjectHashMap<UniqueIdPeakResult[]>> coords = coordinateCache.get();
    if (coords.getKey() != simulationParameters.id) {
        coords = Pair.of(simulationParameters.id, getCoordinates(results));
        coordinateCache.set(coords);
        update = true;
    }
    actualCoordinates = coords.getValue();
    spotFitResults = BenchmarkSpotFit.getBenchmarkSpotFitResults();
    FitResultData localFitResultData = fitResultDataCache.get();
    final SettingsList scoreSettings = new SettingsList(settings.partialMatchDistance, settings.upperMatchDistance, settings.partialSignalFactor, settings.upperSignalFactor);
    final boolean equalScoreSettings = scoreSettings.equals(localFitResultData.scoreSettings);
    if (update || localFitResultData.fittingId != spotFitResults.id || !equalScoreSettings || localFitResultData.differentSettings(settings)) {
        IJ.showStatus("Reading results ...");
        if (localFitResultData.fittingId < 0) {
            // Copy the settings from the fitter if this is the first run.
            // This just starts the plugin with sensible settings.
            // Q. Should this be per new simulation or fitting result instead?
            final FitEngineConfiguration config = BenchmarkSpotFit.getFitEngineConfiguration();
            settings.failCount = config.getFailuresLimit();
            settings.duplicateDistance = config.getDuplicateDistance();
            settings.duplicateDistanceAbsolute = config.getDuplicateDistanceAbsolute();
            settings.residualsThreshold = (BenchmarkSpotFit.getComputeDoublets()) ? BenchmarkSpotFit.getMultiFilter().residualsThreshold : 1;
        }
        // This functionality is for choosing the optimum filter for the given scoring metric.
        if (!equalScoreSettings) {
            filterAnalysisResult.scores.clear();
        }
        localFitResultData = new FitResultData(spotFitResults.id, scoreSettings, settings);
        // @formatter:off
        // -=-=-=-
        // The scoring is designed to find the best fitter+filter combination for the given spot
        // candidates. The ideal combination would correctly fit+pick all the candidate positions
        // that are close to a localisation.
        // 
        // Use the following scoring scheme for all candidates:
        // 
        // Candidates
        // +----------------------------------------+
        // |   Actual matches                       |
        // |  +-----------+                TN       |
        // |  |  FN       |                         |
        // |  |      +----------                    |
        // |  |      | TP |    | Fitted             |
        // |  +-----------+    | spots              |
        // |         |     FP  |                    |
        // |         +---------+                    |
        // +----------------------------------------+
        // 
        // Candidates     = All the spot candidates
        // Actual matches = Any spot candidate or fitted spot candidate that matches a localisation
        // Fitted spots   = Any spot candidate that was successfully fitted
        // 
        // TP = A spot candidate that was fitted and matches a localisation and is accepted
        // FP = A spot candidate that was fitted but does not match a localisation and is accepted
        // FN = A spot candidate that failed to be fitted but matches a localisation
        // = A spot candidate that was fitted and matches a localisation and is rejected
        // TN = A spot candidate that failed to be fitted and does not match a localisation
        // = A spot candidate that was fitted and does not match a localisation and is rejected
        // 
        // When fitting only produces one result it is possible to compute the TN score.
        // Since unfitted candidates can only be TN or FN we could accumulate these scores and cache
        // them. This was the old method of benchmarking single spot fitting and allowed more scores
        // to be computed.
        // 
        // When fitting produces multiple results then we have to score each fit result against all
        // possible actual results and keep a record of the scores. These can then be assessed when
        // the specific results have been chosen by result filtering.
        // 
        // Using a distance ramped scoring function the degree of match can be varied from 0 to 1.
        // Using a signal-factor ramped scoring function the degree of fitted can be varied from 0
        // to 1. When using ramped scoring functions the fractional allocation of scores using the
        // above scheme is performed, i.e. candidates are treated as if they both match and unmatch.
        // This results in an equivalent to multiple analysis using different thresholds and averaging
        // of the scores.
        // 
        // The totals TP+FP+TN+FN must equal the number of spot candidates. This allows different
        // fitting methods to be compared since the total number of candidates is the same.
        // 
        // Precision = TP / (TP+FP)    : This is always valid as a minimum criteria score
        // Recall    = TP / (TP+FN)    : This is valid between different fitting methods since a
        // method that fits more spots will have a potentially lower FN
        // Jaccard   = TP / (TP+FN+FP) : This is valid between fitting methods
        // 
        // -=-=-=-
        // As an alternative scoring system, different fitting methods can be compared using the same
        // TP value but calculating FN = localisations - TP and FP as Positives - TP. This creates a
        // score against the original number of simulated molecules using everything that was passed
        // through the filter (Positives). This score is comparable when a different spot candidate
        // filter has been used and the total number of candidates is different, e.g. Mean filtering
        // vs. Gaussian filtering
        // -=-=-=-
        // @formatter:on
        final RampedScore distanceScore = RampedScore.of(spotFitResults.distanceInPixels * settings.upperMatchDistance / 100.0, spotFitResults.distanceInPixels * settings.partialMatchDistance / 100.0, false);
        localFitResultData.lowerDistanceInPixels = distanceScore.edge1;
        localFitResultData.distanceInPixels = distanceScore.edge0;
        final double matchDistance = MathUtils.pow2(localFitResultData.distanceInPixels);
        localFitResultData.resultsPrefix3 = "\t" + MathUtils.rounded(distanceScore.edge1 * simulationParameters.pixelPitch) + "\t" + MathUtils.rounded(distanceScore.edge0 * simulationParameters.pixelPitch);
        localFitResultData.limitRange = ", d=" + MathUtils.rounded(distanceScore.edge1 * simulationParameters.pixelPitch) + "-" + MathUtils.rounded(distanceScore.edge0 * simulationParameters.pixelPitch);
        // Signal factor must be greater than 1
        final RampedScore signalScore;
        final double spotSignalFactor = BenchmarkSpotFit.getSignalFactor();
        if (spotSignalFactor > 0 && settings.upperSignalFactor > 0) {
            signalScore = RampedScore.of(spotSignalFactor * settings.upperSignalFactor / 100.0, spotSignalFactor * settings.partialSignalFactor / 100.0, false);
            localFitResultData.lowerSignalFactor = signalScore.edge1;
            localFitResultData.signalFactor = signalScore.edge0;
            localFitResultData.resultsPrefix3 += "\t" + MathUtils.rounded(signalScore.edge1) + "\t" + MathUtils.rounded(signalScore.edge0);
            localFitResultData.limitRange += ", s=" + MathUtils.rounded(signalScore.edge1) + "-" + MathUtils.rounded(signalScore.edge0);
        } else {
            signalScore = null;
            localFitResultData.resultsPrefix3 += "\t0\t0";
            localFitResultData.lowerSignalFactor = localFitResultData.signalFactor = 0;
        }
        // Store all the results
        final ArrayList<MultiPathFitResults> multiPathFitResults = new ArrayList<>(spotFitResults.fitResults.size());
        final List<MultiPathFitResults> syncResults = Collections.synchronizedList(multiPathFitResults);
        // This could be multi-threaded ...
        final int nThreads = getThreads(spotFitResults.fitResults.size());
        final BlockingQueue<Job> jobs = new ArrayBlockingQueue<>(nThreads * 2);
        final List<FitResultsWorker> workers = new LinkedList<>();
        final List<Thread> threads = new LinkedList<>();
        final AtomicInteger uniqueId = new AtomicInteger();
        final CoordinateStore localCoordinateStore = createCoordinateStore();
        final Ticker ticker = ImageJUtils.createTicker(spotFitResults.fitResults.size(), nThreads, null);
        for (int i = 0; i < nThreads; i++) {
            final FitResultsWorker worker = new FitResultsWorker(jobs, syncResults, matchDistance, distanceScore, signalScore, uniqueId, localCoordinateStore.newInstance(), ticker, actualCoordinates);
            final Thread t = new Thread(worker);
            workers.add(worker);
            threads.add(t);
            t.start();
        }
        spotFitResults.fitResults.forEachEntry((frame, candidates) -> {
            put(jobs, new Job(frame, candidates));
            return true;
        });
        // Finish all the worker threads by passing in a null job
        for (int i = 0; i < threads.size(); i++) {
            put(jobs, new Job(0, null));
        }
        // Wait for all to finish
        for (int i = 0; i < threads.size(); i++) {
            try {
                threads.get(i).join();
                final FitResultsWorker worker = workers.get(i);
                localFitResultData.matches += worker.matches;
                localFitResultData.fittedResults += worker.included;
                localFitResultData.totalResults += worker.total;
                localFitResultData.notDuplicateCount += worker.notDuplicateCount;
                localFitResultData.newResultCount += worker.newResultCount;
                localFitResultData.countActual += worker.includedActual;
                if (i == 0) {
                    localFitResultData.depthStats = worker.depthStats;
                    localFitResultData.depthFitStats = worker.depthFitStats;
                    localFitResultData.signalFactorStats = worker.signalFactorStats;
                    localFitResultData.distanceStats = worker.distanceStats;
                } else {
                    localFitResultData.depthStats.add(worker.depthStats);
                    localFitResultData.depthFitStats.add(worker.depthFitStats);
                    localFitResultData.signalFactorStats.add(worker.signalFactorStats);
                    localFitResultData.distanceStats.add(worker.distanceStats);
                }
            } catch (final InterruptedException ex) {
                Thread.currentThread().interrupt();
                throw new ConcurrentRuntimeException("Unexpected interrupt", ex);
            }
        }
        threads.clear();
        ImageJUtils.finished();
        localFitResultData.maxUniqueId = uniqueId.get();
        localFitResultData.resultsList = multiPathFitResults.toArray(new MultiPathFitResults[0]);
        Arrays.sort(localFitResultData.resultsList, (o1, o2) -> Integer.compare(o1.getFrame(), o2.getFrame()));
        MultiPathFilter.resetValidationFlag(localFitResultData.resultsList);
        fitResultDataCache.set(localFitResultData);
    }
    fitResultData = localFitResultData;
    return localFitResultData.resultsList;
}