Examples with FractionalAssignment uk.ac.sussex.gdsc.core.match.FractionalAssignment used on opensource projects

Example 6 with FractionalAssignment

use of uk.ac.sussex.gdsc.core.match.FractionalAssignment in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method reportResults.

private ComplexFilterScore reportResults(boolean newResults, List<ComplexFilterScore> filters) {
    if (filters.isEmpty()) {
        IJ.log("Warning: No filters pass the criteria");
        return null;
    }
    getCoordinateStore();
    Collections.sort(filters);
    FractionClassificationResult topFilterClassificationResult = null;
    ArrayList<FractionalAssignment[]> topFilterResults = null;
    String topFilterSummary = null;
    if (settings.showSummaryTable || settings.saveTemplate) {
        final Consumer<String> summaryWindow = createSummaryWindow();
        int count = 0;
        final double range = (settings.summaryDepth / simulationParameters.pixelPitch) * 0.5;
        final int[] np = { 0 };
        fitResultData.depthStats.forEach(depth -> {
            if (Math.abs(depth) < range) {
                np[0]++;
            }
        });
        for (final ComplexFilterScore fs : filters) {
            final ArrayList<FractionalAssignment[]> list = new ArrayList<>(fitResultData.resultsList.length);
            final FractionClassificationResult r = scoreFilter(fs.getFilter(), defaultMinimalFilter, fitResultData.resultsList, list, coordinateStore);
            final StringBuilder sb = createResult(fs.getFilter(), r);
            if (topFilterResults == null) {
                topFilterResults = list;
                topFilterClassificationResult = r;
            }
            // Show the recall at the specified depth. Sum the distance and signal factor of all scored
            // spots.
            int scored = 0;
            double tp = 0;
            double distance = 0;
            double sf = 0;
            double rmsd = 0;
            final SimpleRegression regression = new SimpleRegression(false);
            for (final FractionalAssignment[] assignments : list) {
                if (assignments == null) {
                    continue;
                }
                for (int i = 0; i < assignments.length; i++) {
                    final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
                    if (Math.abs(c.peak.getZPosition()) <= range) {
                        tp += c.getScore();
                    }
                    distance += c.distToTarget;
                    sf += c.getSignalFactor();
                    rmsd += c.distToTarget * c.distToTarget;
                    regression.addData(c.peakResult.getSignal(), c.peak.getIntensity());
                }
                scored += assignments.length;
            }
            final double slope = regression.getSlope();
            sb.append('\t');
            sb.append(MathUtils.rounded(tp / np[0])).append('\t');
            sb.append(MathUtils.rounded(distance / scored)).append('\t');
            sb.append(MathUtils.rounded(sf / scored)).append('\t');
            // RMSD to be the root mean square deviation in a single dimension so divide by 2.
            // (This assumes 2D Euclidean distances.)
            sb.append(MathUtils.rounded(Math.sqrt(MathUtils.div0(rmsd / 2, scored)))).append('\t');
            sb.append(MathUtils.rounded(slope)).append('\t');
            if (fs.atLimit() != null) {
                sb.append(fs.atLimit());
            }
            String text = sb.toString();
            if (topFilterSummary == null) {
                topFilterSummary = text;
                if (!settings.showSummaryTable) {
                    break;
                }
            }
            if (fs.time != 0) {
                sb.append('\t');
                sb.append(fs.algorithm);
                sb.append('\t');
                sb.append(org.apache.commons.lang3.time.DurationFormatUtils.formatDurationHMS(fs.time));
            } else {
                sb.append("\t\t");
            }
            if (fs.paramTime != 0) {
                sb.append('\t');
                sb.append(fs.getParamAlgorithm());
                sb.append('\t');
                sb.append(org.apache.commons.lang3.time.DurationFormatUtils.formatDurationHMS(fs.paramTime));
            } else {
                sb.append("\t\t");
            }
            text = sb.toString();
            summaryWindow.accept(text);
            count++;
            if (settings.summaryTopN > 0 && count >= settings.summaryTopN) {
                break;
            }
        }
        // Add a spacer to the summary table if we have multiple results
        if (count > 1 && settings.showSummaryTable) {
            summaryWindow.accept("");
        }
    }
    final DirectFilter localBestFilter = filters.get(0).getFilter();
    if (settings.saveBestFilter) {
        saveFilter(localBestFilter);
    }
    if (topFilterClassificationResult == null) {
        topFilterResults = new ArrayList<>(fitResultData.resultsList.length);
        scoreFilter(localBestFilter, defaultMinimalFilter, fitResultData.resultsList, topFilterResults, coordinateStore);
    }
    if (newResults || filterAnalysisResult.scores.isEmpty()) {
        filterAnalysisResult.scores.add(new FilterResult(settings.failCount, residualsThreshold, settings.duplicateDistance, settings.duplicateDistanceAbsolute, filters.get(0)));
    }
    if (settings.saveTemplate) {
        saveTemplate(topFilterSummary);
    }
    showPlots();
    calculateSensitivity();
    topFilterResults = depthAnalysis(topFilterResults, localBestFilter);
    topFilterResults = scoreAnalysis(topFilterResults, localBestFilter);
    componentAnalysis(filters.get(0));
    PreprocessedPeakResult[] filterResults = null;
    if (isShowOverlay()) {
        filterResults = showOverlay(topFilterResults, localBestFilter);
    }
    saveResults(filterResults, localBestFilter);
    wo.tile();
    return filters.get(0);
}

Example 7 with FractionalAssignment

use of uk.ac.sussex.gdsc.core.match.FractionalAssignment in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method depthAnalysis.

/**
 * Depth analysis.
 *
 * @param allAssignments The assignments generated from running the filter (or null)
 * @param filter the filter
 * @return the assignments
 */
@Nullable
private ArrayList<FractionalAssignment[]> depthAnalysis(ArrayList<FractionalAssignment[]> allAssignments, DirectFilter filter) {
    if (!settings.depthRecallAnalysis || simulationParameters.fixedDepth) {
        return null;
    }
    // Build a histogram of the number of spots at different depths
    final double[] depths = fitResultData.depthStats.getValues();
    double[] limits = MathUtils.limits(depths);
    final int bins = HistogramPlot.getBinsSqrtRule(depths.length);
    final double[][] h1 = HistogramPlot.calcHistogram(depths, limits[0], limits[1], bins);
    final double[][] h2 = HistogramPlot.calcHistogram(fitResultData.depthFitStats.getValues(), limits[0], limits[1], bins);
    // manually to get the results that pass.
    if (allAssignments == null) {
        allAssignments = getAssignments(filter);
    }
    double[] depths2 = new double[results.size()];
    int count = 0;
    for (final FractionalAssignment[] assignments : allAssignments) {
        if (assignments == null) {
            continue;
        }
        for (int i = 0; i < assignments.length; i++) {
            final CustomFractionalAssignment c = (CustomFractionalAssignment) assignments[i];
            depths2[count++] = c.peak.getZPosition();
        }
    }
    depths2 = Arrays.copyOf(depths2, count);
    // Build a histogram using the same limits
    final double[][] h3 = HistogramPlot.calcHistogram(depths2, limits[0], limits[1], bins);
    // Convert pixel depth to nm
    for (int i = 0; i < h1[0].length; i++) {
        h1[0][i] *= simulationParameters.pixelPitch;
    }
    limits[0] *= simulationParameters.pixelPitch;
    limits[1] *= simulationParameters.pixelPitch;
    // Produce a histogram of the number of spots at each depth
    final String title1 = TITLE + " Depth Histogram";
    final Plot plot1 = new Plot(title1, "Depth (nm)", "Frequency");
    plot1.setLimits(limits[0], limits[1], 0, MathUtils.max(h1[1]));
    plot1.setColor(Color.black);
    plot1.addPoints(h1[0], h1[1], Plot.BAR);
    plot1.addLabel(0, 0, "Black = Spots; Blue = Fitted; Red = Filtered");
    plot1.setColor(Color.blue);
    plot1.addPoints(h1[0], h2[1], Plot.BAR);
    plot1.setColor(Color.red);
    plot1.addPoints(h1[0], h3[1], Plot.BAR);
    plot1.setColor(Color.magenta);
    ImageJUtils.display(title1, plot1, wo);
    // Interpolate
    final double halfBinWidth = (h1[0][1] - h1[0][0]) * 0.5;
    // Remove final value of the histogram as this is at the upper limit of the range (i.e. count
    // zero)
    h1[0] = Arrays.copyOf(h1[0], h1[0].length - 1);
    h1[1] = Arrays.copyOf(h1[1], h1[0].length);
    h2[1] = Arrays.copyOf(h2[1], h1[0].length);
    h3[1] = Arrays.copyOf(h3[1], h1[0].length);
    // TODO : Fix the smoothing since LOESS sometimes does not work.
    // Perhaps allow configuration of the number of histogram bins and the smoothing bandwidth.
    // Use minimum of 3 points for smoothing
    // Ensure we use at least x% of data
    final double bandwidth = Math.max(3.0 / h1[0].length, 0.15);
    final LoessInterpolator loess = new LoessInterpolator(bandwidth, 1);
    final PolynomialSplineFunction spline1 = loess.interpolate(h1[0], h1[1]);
    final PolynomialSplineFunction spline2 = loess.interpolate(h1[0], h2[1]);
    final PolynomialSplineFunction spline3 = loess.interpolate(h1[0], h3[1]);
    // Use a second interpolator in case the LOESS fails
    final LinearInterpolator lin = new LinearInterpolator();
    final PolynomialSplineFunction spline1b = lin.interpolate(h1[0], h1[1]);
    final PolynomialSplineFunction spline2b = lin.interpolate(h1[0], h2[1]);
    final PolynomialSplineFunction spline3b = lin.interpolate(h1[0], h3[1]);
    // Increase the number of points to show a smooth curve
    final double[] points = new double[bins * 5];
    limits = MathUtils.limits(h1[0]);
    final double interval = (limits[1] - limits[0]) / (points.length - 1);
    final double[] v = new double[points.length];
    final double[] v2 = new double[points.length];
    final double[] v3 = new double[points.length];
    for (int i = 0; i < points.length - 1; i++) {
        points[i] = limits[0] + i * interval;
        v[i] = getSplineValue(spline1, spline1b, points[i]);
        v2[i] = getSplineValue(spline2, spline2b, points[i]);
        v3[i] = getSplineValue(spline3, spline3b, points[i]);
        points[i] += halfBinWidth;
    }
    // Final point on the limit of the spline range
    final int ii = points.length - 1;
    v[ii] = getSplineValue(spline1, spline1b, limits[1]);
    v2[ii] = getSplineValue(spline2, spline2b, limits[1]);
    v3[ii] = getSplineValue(spline3, spline3b, limits[1]);
    points[ii] = limits[1] + halfBinWidth;
    // Calculate recall
    for (int i = 0; i < v.length; i++) {
        v2[i] = v2[i] / v[i];
        v3[i] = v3[i] / v[i];
    }
    final double halfSummaryDepth = settings.summaryDepth * 0.5;
    final String title2 = TITLE + " Depth Histogram (normalised)";
    final Plot plot2 = new Plot(title2, "Depth (nm)", "Recall");
    plot2.setLimits(limits[0] + halfBinWidth, limits[1] + halfBinWidth, 0, MathUtils.min(1, MathUtils.max(v2)));
    plot2.setColor(Color.black);
    plot2.addLabel(0, 0, "Blue = Fitted; Red = Filtered");
    plot2.setColor(Color.blue);
    plot2.addPoints(points, v2, Plot.LINE);
    plot2.setColor(Color.red);
    plot2.addPoints(points, v3, Plot.LINE);
    plot2.setColor(Color.magenta);
    if (-halfSummaryDepth - halfBinWidth >= limits[0]) {
        plot2.drawLine(-halfSummaryDepth, 0, -halfSummaryDepth, getSplineValue(spline3, spline3b, -halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, -halfSummaryDepth - halfBinWidth));
    }
    if (halfSummaryDepth - halfBinWidth <= limits[1]) {
        plot2.drawLine(halfSummaryDepth, 0, halfSummaryDepth, getSplineValue(spline3, spline3b, halfSummaryDepth - halfBinWidth) / getSplineValue(spline1, spline1b, halfSummaryDepth - halfBinWidth));
    }
    ImageJUtils.display(title2, plot2, wo);
    return allAssignments;
}