Search in sources :

Example 1 with PrecisionMethod

use of uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod in project GDSC-SMLM by aherbert.

the class FitConfiguration method getFilterPrecisionMethod.

/**
 * Gets the precision method that will be used to produce the precision value for filtering.
 *
 * <p>This checks first the direct filter and then the current value for the precision method.
 *
 * <p>The value is written into the current calibration to allow the calibration to be used in
 * saved results.
 *
 * @return the filter precision method
 */
public PrecisionMethod getFilterPrecisionMethod() {
    PrecisionMethod method = getFilterPrecisionMethodInternal();
    switch(method) {
        case MORTENSEN:
            varianceSelector = new VarianceSelector();
            break;
        case MORTENSEN_LOCAL_BACKGROUND:
            varianceSelector = new VarianceSelector2();
            break;
        case POISSON_CRLB:
            varianceSelector = new VarianceSelectorCrlb();
            break;
        default:
            method = PrecisionMethod.PRECISION_METHOD_NA;
            varianceSelector = new BaseVarianceSelector();
            break;
    }
    calibration.setPrecisionMethod(method);
    return method;
}
Also used : PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod)

Example 2 with PrecisionMethod

use of uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod in project GDSC-SMLM by aherbert.

the class ImageJ3DResultsViewer method createSphereSizeFromPrecision.

@Nullable
private static Point3f[] createSphereSizeFromPrecision(MemoryPeakResults results) {
    final PrecisionResultProcedure p = new PrecisionResultProcedure(results);
    try {
        final PrecisionMethod m = p.getPrecision();
        IJ.log("Using precision method " + FitProtosHelper.getName(m));
        final Point3f[] size = new Point3f[results.size()];
        for (int i = 0, j = 0; i < p.precisions.length; i++) {
            // Precision is in NM which matches the rendering
            final float v = (float) p.precisions[i];
            size[j++] = new Point3f(v, v, v);
        }
        return size;
    } catch (final DataException ex) {
        IJ.error(TITLE, "The results have no precision: " + ex.getMessage());
        return null;
    }
}
Also used : DataException(uk.ac.sussex.gdsc.core.data.DataException) Point3f(org.scijava.vecmath.Point3f) PrecisionResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure) PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 3 with PrecisionMethod

use of uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod in project GDSC-SMLM by aherbert.

the class SummariseResults method createSummary.

private static String createSummary(StringBuilder sb, MemoryPeakResults result, int[] removeNullResults) {
    sb.setLength(0);
    final DescriptiveStatistics[] stats = new DescriptiveStatistics[2];
    for (int i = 0; i < stats.length; i++) {
        stats[i] = new DescriptiveStatistics();
    }
    if (result.hasNullResults()) {
        IJ.log("Null results in dataset: " + result.getName());
        if (removeNullResults[0] == UNKNOWN) {
            final GenericDialog gd = new GenericDialog(TITLE);
            gd.addMessage("There are invalid results in memory.\n \nClean these results?");
            gd.enableYesNoCancel();
            gd.hideCancelButton();
            gd.showDialog();
            removeNullResults[0] = (gd.wasOKed()) ? YES : NO;
        }
        if (removeNullResults[0] == NO) {
            result = result.copy();
        }
        result.removeNullResults();
    }
    final CalibrationReader calibration = result.getCalibrationReader();
    PrecisionMethod precisionMethod = PrecisionMethod.PRECISION_METHOD_NA;
    boolean stored = false;
    final int size = result.size();
    if (size > 0) {
        // Precision
        try {
            final PrecisionResultProcedure p = new PrecisionResultProcedure(result);
            // Use stored precision if possible
            stored = result.hasPrecision();
            precisionMethod = p.getPrecision(stored);
            for (final double v : p.precisions) {
                stats[0].addValue(v);
            }
        } catch (final DataException ex) {
        // Ignore
        }
        // SNR
        try {
            final SnrResultProcedure p = new SnrResultProcedure(result);
            p.getSnr();
            for (final double v : p.snr) {
                stats[1].addValue(v);
            }
        } catch (final DataException ex) {
        // Ignore
        }
    }
    sb.append(result.getName());
    int maxT = 0;
    if (result.size() == 0) {
        sb.append("\t0\t0");
    } else {
        sb.append('\t').append(result.size());
        maxT = result.getMaxFrame();
        sb.append('\t').append(maxT);
    }
    if (calibration != null && calibration.hasExposureTime()) {
        sb.append('\t').append(TextUtils.millisToString((long) Math.ceil(maxT * calibration.getExposureTime())));
    } else {
        sb.append("\t-");
    }
    if (size > 0) {
        final boolean includeDeviations = result.hasDeviations();
        final long memorySize = MemoryPeakResults.estimateMemorySize(size, includeDeviations);
        final String memory = TextUtils.bytesToString(memorySize);
        sb.append('\t').append(memory);
    } else {
        sb.append("\t-");
    }
    final Rectangle bounds = result.getBounds(true);
    TextUtils.formatTo(sb, "\t%d,%d,%d,%d", bounds.x, bounds.y, bounds.x + bounds.width, bounds.y + bounds.height);
    if (calibration != null) {
        sb.append('\t').append(calibration.hasNmPerPixel() ? MathUtils.rounded(calibration.getNmPerPixel()) : '-');
        sb.append('\t').append(calibration.hasExposureTime() ? MathUtils.rounded(calibration.getExposureTime()) : '-');
        if (calibration.hasCameraType()) {
            sb.append('\t').append(CalibrationProtosHelper.getName(calibration.getCameraType()));
            if (calibration.isCcdCamera()) {
                sb.append(" bias=").append(calibration.getBias());
                sb.append(" gain=").append(calibration.getCountPerPhoton());
            }
        } else {
            sb.append("\t-");
        }
        sb.append('\t').append(calibration.hasDistanceUnit() ? UnitHelper.getShortName(calibration.getDistanceUnit()) : '-');
        sb.append('\t').append(calibration.hasIntensityUnit() ? UnitHelper.getShortName(calibration.getIntensityUnit()) : '-');
    } else {
        sb.append("\t\t\t\t\t");
    }
    if (result.is3D()) {
        sb.append("\tY");
    } else {
        sb.append("\tN");
    }
    sb.append("\t").append(FitProtosHelper.getName(precisionMethod));
    if (stored) {
        sb.append(" (Stored)");
    }
    for (int i = 0; i < stats.length; i++) {
        if (Double.isNaN(stats[i].getMean())) {
            sb.append("\t-\t-\t-\t-");
        } else {
            sb.append('\t').append(IJ.d2s(stats[i].getMean(), 3));
            sb.append('\t').append(IJ.d2s(stats[i].getPercentile(50), 3));
            sb.append('\t').append(IJ.d2s(stats[i].getMin(), 3));
            sb.append('\t').append(IJ.d2s(stats[i].getMax(), 3));
        }
    }
    return sb.toString();
}
Also used : DescriptiveStatistics(org.apache.commons.math3.stat.descriptive.DescriptiveStatistics) Rectangle(java.awt.Rectangle) CalibrationReader(uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader) PrecisionResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure) PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod) DataException(uk.ac.sussex.gdsc.core.data.DataException) GenericDialog(ij.gui.GenericDialog) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) SnrResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.SnrResultProcedure)

Example 4 with PrecisionMethod

use of uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod in project GDSC-SMLM by aherbert.

the class BenchmarkSpotFit method summariseResults.

private void summariseResults(BenchmarkSpotFitResult spotFitResults, long runTime, final PreprocessedPeakResult[] preprocessedPeakResults, int uniqueIdCount, CandidateData candidateData, TIntObjectHashMap<List<Coordinate>> actualCoordinates) {
    // Summarise the fitting results. N fits, N failures.
    // Optimal match statistics if filtering is perfect (since fitting is not perfect).
    final StoredDataStatistics distanceStats = new StoredDataStatistics();
    final StoredDataStatistics depthStats = new StoredDataStatistics();
    // Get stats for all fitted results and those that match
    // Signal, SNR, Width, xShift, yShift, Precision
    createFilterCriteria();
    final StoredDataStatistics[][] stats = new StoredDataStatistics[3][filterCriteria.length];
    for (int i = 0; i < stats.length; i++) {
        for (int j = 0; j < stats[i].length; j++) {
            stats[i][j] = new StoredDataStatistics();
        }
    }
    final double nmPerPixel = simulationParameters.pixelPitch;
    double tp = 0;
    double fp = 0;
    int failCtp = 0;
    int failCfp = 0;
    int ctp = 0;
    int cfp = 0;
    final int[] singleStatus = new int[FitStatus.values().length];
    final int[] multiStatus = new int[singleStatus.length];
    final int[] doubletStatus = new int[singleStatus.length];
    final int[] multiDoubletStatus = new int[singleStatus.length];
    // Easier to materialise the values since we have a lot of non final variables to manipulate
    final TIntObjectHashMap<FilterCandidates> fitResults = spotFitResults.fitResults;
    final int[] frames = new int[fitResults.size()];
    final FilterCandidates[] candidates = new FilterCandidates[fitResults.size()];
    final int[] counter = new int[1];
    fitResults.forEachEntry((frame, candidate) -> {
        frames[counter[0]] = frame;
        candidates[counter[0]] = candidate;
        counter[0]++;
        return true;
    });
    for (final FilterCandidates result : candidates) {
        // Count the number of fit results that matched (tp) and did not match (fp)
        tp += result.tp;
        fp += result.fp;
        for (int i = 0; i < result.fitResult.length; i++) {
            if (result.spots[i].match) {
                ctp++;
            } else {
                cfp++;
            }
            final MultiPathFitResult fitResult = result.fitResult[i];
            if (singleStatus != null && result.spots[i].match) {
                // Debugging reasons for fit failure
                addStatus(singleStatus, fitResult.getSingleFitResult());
                addStatus(multiStatus, fitResult.getMultiFitResult());
                addStatus(doubletStatus, fitResult.getDoubletFitResult());
                addStatus(multiDoubletStatus, fitResult.getMultiDoubletFitResult());
            }
            if (noMatch(fitResult)) {
                if (result.spots[i].match) {
                    failCtp++;
                } else {
                    failCfp++;
                }
            }
            // We have multi-path results.
            // We want statistics for:
            // [0] all fitted spots
            // [1] fitted spots that match a result
            // [2] fitted spots that do not match a result
            addToStats(fitResult.getSingleFitResult(), stats);
            addToStats(fitResult.getMultiFitResult(), stats);
            addToStats(fitResult.getDoubletFitResult(), stats);
            addToStats(fitResult.getMultiDoubletFitResult(), stats);
        }
        // Statistics on spots that fit an actual result
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            distanceStats.add(fitMatch.distance * nmPerPixel);
            depthStats.add(fitMatch.zdepth * nmPerPixel);
        }
    }
    if (tp == 0) {
        IJ.error(TITLE, "No fit results matched the simulation actual results");
        return;
    }
    // Store data for computing correlation
    final double[] i1 = new double[depthStats.getN()];
    final double[] i2 = new double[i1.length];
    final double[] is = new double[i1.length];
    int ci = 0;
    for (final FilterCandidates result : candidates) {
        for (int i = 0; i < result.match.length; i++) {
            if (!result.match[i].isFitResult()) {
                // For now just ignore the candidates that matched
                continue;
            }
            final FitMatch fitMatch = (FitMatch) result.match[i];
            final ScoredSpot spot = result.spots[fitMatch.index];
            i1[ci] = fitMatch.predictedSignal;
            i2[ci] = fitMatch.actualSignal;
            is[ci] = spot.spot.intensity;
            ci++;
        }
    }
    // We want to compute the Jaccard against the spot metric
    // Filter the results using the multi-path filter
    final ArrayList<MultiPathFitResults> multiPathResults = new ArrayList<>(fitResults.size());
    for (int i = 0; i < frames.length; i++) {
        final int frame = frames[i];
        final MultiPathFitResult[] multiPathFitResults = candidates[i].fitResult;
        final int totalCandidates = candidates[i].spots.length;
        final List<Coordinate> list = actualCoordinates.get(frame);
        final int nActual = (list == null) ? 0 : list.size();
        multiPathResults.add(new MultiPathFitResults(frame, multiPathFitResults, totalCandidates, nActual));
    }
    // Score the results and count the number returned
    final List<FractionalAssignment[]> assignments = new ArrayList<>();
    final TIntHashSet set = new TIntHashSet(uniqueIdCount);
    final FractionScoreStore scoreStore = set::add;
    final MultiPathFitResults[] multiResults = multiPathResults.toArray(new MultiPathFitResults[0]);
    // Filter with no filter
    final MultiPathFilter mpf = new MultiPathFilter(new SignalFilter(0), null, multiFilter.residualsThreshold);
    mpf.fractionScoreSubset(multiResults, NullFailCounter.INSTANCE, this.results.size(), assignments, scoreStore, CoordinateStoreFactory.create(0, 0, imp.getWidth(), imp.getHeight(), config.convertUsingHwhMax(config.getDuplicateDistanceParameter())));
    final double[][] matchScores = new double[set.size()][];
    int count = 0;
    for (int i = 0; i < assignments.size(); i++) {
        final FractionalAssignment[] a = assignments.get(i);
        if (a == null) {
            continue;
        }
        for (int j = 0; j < a.length; j++) {
            final PreprocessedPeakResult r = ((PeakFractionalAssignment) a[j]).peakResult;
            set.remove(r.getUniqueId());
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            score[FILTER_PRECISION + 1] = a[j].getScore();
            matchScores[count++] = score;
        }
    }
    // Add the rest
    set.forEach(new CustomTIntProcedure(count) {

        @Override
        public boolean execute(int uniqueId) {
            // This should not be null or something has gone wrong
            final PreprocessedPeakResult r = preprocessedPeakResults[uniqueId];
            if (r == null) {
                throw new IllegalArgumentException("Missing result: " + uniqueId);
            }
            final double precision = Math.sqrt(r.getLocationVariance());
            final double signal = r.getSignal();
            final double snr = r.getSnr();
            final double width = r.getXSdFactor();
            final double xShift = r.getXRelativeShift2();
            final double yShift = r.getYRelativeShift2();
            // Since these two are combined for filtering and the max is what matters.
            final double shift = (xShift > yShift) ? Math.sqrt(xShift) : Math.sqrt(yShift);
            final double eshift = Math.sqrt(xShift + yShift);
            final double[] score = new double[8];
            score[FILTER_SIGNAL] = signal;
            score[FILTER_SNR] = snr;
            score[FILTER_MIN_WIDTH] = width;
            score[FILTER_MAX_WIDTH] = width;
            score[FILTER_SHIFT] = shift;
            score[FILTER_ESHIFT] = eshift;
            score[FILTER_PRECISION] = precision;
            matchScores[count++] = score;
            return true;
        }
    });
    final FitConfiguration fitConfig = config.getFitConfiguration();
    // Debug the reasons the fit failed
    if (singleStatus != null) {
        String name = PeakFit.getSolverName(fitConfig);
        if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
            name += " Camera";
        }
        IJ.log("Failure counts: " + name);
        printFailures("Single", singleStatus);
        printFailures("Multi", multiStatus);
        printFailures("Doublet", doubletStatus);
        printFailures("Multi doublet", multiDoubletStatus);
    }
    final StringBuilder sb = new StringBuilder(300);
    // Add information about the simulation
    final double signal = simulationParameters.averageSignal;
    final int n = results.size();
    sb.append(imp.getStackSize()).append('\t');
    final int w = imp.getWidth();
    final int h = imp.getHeight();
    sb.append(w).append('\t');
    sb.append(h).append('\t');
    sb.append(n).append('\t');
    final double density = ((double) n / imp.getStackSize()) / (w * h) / (simulationParameters.pixelPitch * simulationParameters.pixelPitch / 1e6);
    sb.append(MathUtils.rounded(density)).append('\t');
    sb.append(MathUtils.rounded(signal)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.pixelPitch)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.depth)).append('\t');
    sb.append(simulationParameters.fixedDepth).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.gain)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.readNoise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.background)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.noise)).append('\t');
    if (simulationParameters.fullSimulation) {
    // The total signal is spread over frames
    }
    sb.append(MathUtils.rounded(signal / simulationParameters.noise)).append('\t');
    sb.append(MathUtils.rounded(simulationParameters.sd / simulationParameters.pixelPitch)).append('\t');
    sb.append(spotFilter.getDescription());
    // nP and nN is the fractional score of the spot candidates
    addCount(sb, (double) candidateData.countPositive + candidateData.countNegative);
    addCount(sb, candidateData.countPositive);
    addCount(sb, candidateData.countNegative);
    addCount(sb, candidateData.fractionPositive);
    addCount(sb, candidateData.fractionNegative);
    String name = PeakFit.getSolverName(fitConfig);
    if (fitConfig.getFitSolver() == FitSolver.MLE && fitConfig.isModelCamera()) {
        name += " Camera";
    }
    add(sb, name);
    add(sb, config.getFitting());
    spotFitResults.resultPrefix = sb.toString();
    // Q. Should I add other fit configuration here?
    // The fraction of positive and negative candidates that were included
    add(sb, (100.0 * ctp) / candidateData.countPositive);
    add(sb, (100.0 * cfp) / candidateData.countNegative);
    // Score the fitting results compared to the original simulation.
    // Score the candidate selection:
    add(sb, ctp + cfp);
    add(sb, ctp);
    add(sb, cfp);
    // TP are all candidates that can be matched to a spot
    // FP are all candidates that cannot be matched to a spot
    // FN = The number of missed spots
    FractionClassificationResult match = new FractionClassificationResult(ctp, cfp, 0, simulationParameters.molecules - ctp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Score the fitting results:
    add(sb, failCtp);
    add(sb, failCfp);
    // TP are all fit results that can be matched to a spot
    // FP are all fit results that cannot be matched to a spot
    // FN = The number of missed spots
    add(sb, tp);
    add(sb, fp);
    match = new FractionClassificationResult(tp, fp, 0, simulationParameters.molecules - tp);
    add(sb, match.getRecall());
    add(sb, match.getPrecision());
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // Do it again but pretend we can perfectly filter all the false positives
    // add(sb, tp);
    match = new FractionClassificationResult(tp, 0, 0, simulationParameters.molecules - tp);
    // Recall is unchanged
    // Precision will be 100%
    add(sb, match.getF1Score());
    add(sb, match.getJaccard());
    // The mean may be subject to extreme outliers so use the median
    double median = distanceStats.getMedian();
    add(sb, median);
    final WindowOrganiser wo = new WindowOrganiser();
    String label = String.format("Recall = %s. n = %d. Median = %s nm. SD = %s nm", MathUtils.rounded(match.getRecall()), distanceStats.getN(), MathUtils.rounded(median), MathUtils.rounded(distanceStats.getStandardDeviation()));
    new HistogramPlotBuilder(TITLE, distanceStats, "Match Distance (nm)").setPlotLabel(label).show(wo);
    median = depthStats.getMedian();
    add(sb, median);
    // Sort by spot intensity and produce correlation
    double[] correlation = null;
    double[] rankCorrelation = null;
    double[] rank = null;
    final FastCorrelator fastCorrelator = new FastCorrelator();
    final ArrayList<Ranking> pc1 = new ArrayList<>();
    final ArrayList<Ranking> pc2 = new ArrayList<>();
    ci = 0;
    if (settings.showCorrelation) {
        final int[] indices = SimpleArrayUtils.natural(i1.length);
        SortUtils.sortData(indices, is, settings.rankByIntensity, true);
        correlation = new double[i1.length];
        rankCorrelation = new double[i1.length];
        rank = new double[i1.length];
        for (final int ci2 : indices) {
            fastCorrelator.add(Math.round(i1[ci2]), Math.round(i2[ci2]));
            pc1.add(new Ranking(i1[ci2], ci));
            pc2.add(new Ranking(i2[ci2], ci));
            correlation[ci] = fastCorrelator.getCorrelation();
            rankCorrelation[ci] = Correlator.correlation(rank(pc1), rank(pc2));
            if (settings.rankByIntensity) {
                rank[ci] = is[0] - is[ci];
            } else {
                rank[ci] = ci;
            }
            ci++;
        }
    } else {
        for (int i = 0; i < i1.length; i++) {
            fastCorrelator.add(Math.round(i1[i]), Math.round(i2[i]));
            pc1.add(new Ranking(i1[i], i));
            pc2.add(new Ranking(i2[i], i));
        }
    }
    final double pearsonCorr = fastCorrelator.getCorrelation();
    final double rankedCorr = Correlator.correlation(rank(pc1), rank(pc2));
    // Get the regression
    final SimpleRegression regression = new SimpleRegression(false);
    for (int i = 0; i < pc1.size(); i++) {
        regression.addData(pc1.get(i).value, pc2.get(i).value);
    }
    // final double intercept = regression.getIntercept();
    final double slope = regression.getSlope();
    if (settings.showCorrelation) {
        String title = TITLE + " Intensity";
        Plot plot = new Plot(title, "Candidate", "Spot");
        final double[] limits1 = MathUtils.limits(i1);
        final double[] limits2 = MathUtils.limits(i2);
        plot.setLimits(limits1[0], limits1[1], limits2[0], limits2[1]);
        label = String.format("Correlation=%s; Ranked=%s; Slope=%s", MathUtils.rounded(pearsonCorr), MathUtils.rounded(rankedCorr), MathUtils.rounded(slope));
        plot.addLabel(0, 0, label);
        plot.setColor(Color.red);
        plot.addPoints(i1, i2, Plot.DOT);
        if (slope > 1) {
            plot.drawLine(limits1[0], limits1[0] * slope, limits1[1], limits1[1] * slope);
        } else {
            plot.drawLine(limits2[0] / slope, limits2[0], limits2[1] / slope, limits2[1]);
        }
        ImageJUtils.display(title, plot, wo);
        title = TITLE + " Correlation";
        plot = new Plot(title, "Spot Rank", "Correlation");
        final double[] xlimits = MathUtils.limits(rank);
        double[] ylimits = MathUtils.limits(correlation);
        ylimits = MathUtils.limits(ylimits, rankCorrelation);
        plot.setLimits(xlimits[0], xlimits[1], ylimits[0], ylimits[1]);
        plot.setColor(Color.red);
        plot.addPoints(rank, correlation, Plot.LINE);
        plot.setColor(Color.blue);
        plot.addPoints(rank, rankCorrelation, Plot.LINE);
        plot.setColor(Color.black);
        plot.addLabel(0, 0, label);
        ImageJUtils.display(title, plot, wo);
    }
    add(sb, pearsonCorr);
    add(sb, rankedCorr);
    add(sb, slope);
    label = String.format("n = %d. Median = %s nm", depthStats.getN(), MathUtils.rounded(median));
    new HistogramPlotBuilder(TITLE, depthStats, "Match Depth (nm)").setRemoveOutliersOption(1).setPlotLabel(label).show(wo);
    // Plot histograms of the stats on the same window
    final double[] lower = new double[filterCriteria.length];
    final double[] upper = new double[lower.length];
    final double[] min = new double[lower.length];
    final double[] max = new double[lower.length];
    for (int i = 0; i < stats[0].length; i++) {
        final double[] limits = showDoubleHistogram(stats, i, wo, matchScores);
        lower[i] = limits[0];
        upper[i] = limits[1];
        min[i] = limits[2];
        max[i] = limits[3];
    }
    // Reconfigure some of the range limits
    // Make this a bit bigger
    upper[FILTER_SIGNAL] *= 2;
    // Make this a bit bigger
    upper[FILTER_SNR] *= 2;
    final double factor = 0.25;
    if (lower[FILTER_MIN_WIDTH] != 0) {
        // (assuming lower is less than 1)
        upper[FILTER_MIN_WIDTH] = 1 - Math.max(0, factor * (1 - lower[FILTER_MIN_WIDTH]));
    }
    if (upper[FILTER_MIN_WIDTH] != 0) {
        // (assuming upper is more than 1)
        lower[FILTER_MAX_WIDTH] = 1 + Math.max(0, factor * (upper[FILTER_MAX_WIDTH] - 1));
    }
    // Round the ranges
    final double[] interval = new double[stats[0].length];
    interval[FILTER_SIGNAL] = SignalFilter.DEFAULT_INCREMENT;
    interval[FILTER_SNR] = SnrFilter.DEFAULT_INCREMENT;
    interval[FILTER_MIN_WIDTH] = WidthFilter2.DEFAULT_MIN_INCREMENT;
    interval[FILTER_MAX_WIDTH] = WidthFilter.DEFAULT_INCREMENT;
    interval[FILTER_SHIFT] = ShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_ESHIFT] = EShiftFilter.DEFAULT_INCREMENT;
    interval[FILTER_PRECISION] = PrecisionFilter.DEFAULT_INCREMENT;
    interval[FILTER_ITERATIONS] = 0.1;
    interval[FILTER_EVALUATIONS] = 0.1;
    // Create a range increment
    final double[] increment = new double[lower.length];
    for (int i = 0; i < increment.length; i++) {
        lower[i] = MathUtils.floor(lower[i], interval[i]);
        upper[i] = MathUtils.ceil(upper[i], interval[i]);
        final double range = upper[i] - lower[i];
        // Allow clipping if the range is small compared to the min increment
        double multiples = range / interval[i];
        // Use 8 multiples for the equivalent of +/- 4 steps around the centre
        if (multiples < 8) {
            multiples = Math.ceil(multiples);
        } else {
            multiples = 8;
        }
        increment[i] = MathUtils.ceil(range / multiples, interval[i]);
        if (i == FILTER_MIN_WIDTH) {
            // Requires clipping based on the upper limit
            lower[i] = upper[i] - increment[i] * multiples;
        } else {
            upper[i] = lower[i] + increment[i] * multiples;
        }
    }
    for (int i = 0; i < stats[0].length; i++) {
        lower[i] = MathUtils.round(lower[i]);
        upper[i] = MathUtils.round(upper[i]);
        min[i] = MathUtils.round(min[i]);
        max[i] = MathUtils.round(max[i]);
        increment[i] = MathUtils.round(increment[i]);
        sb.append('\t').append(min[i]).append(':').append(lower[i]).append('-').append(upper[i]).append(':').append(max[i]);
    }
    // Disable some filters
    increment[FILTER_SIGNAL] = Double.POSITIVE_INFINITY;
    // increment[FILTER_SHIFT] = Double.POSITIVE_INFINITY;
    increment[FILTER_ESHIFT] = Double.POSITIVE_INFINITY;
    wo.tile();
    sb.append('\t').append(TextUtils.nanosToString(runTime));
    createTable().append(sb.toString());
    if (settings.saveFilterRange) {
        GUIFilterSettings filterSettings = SettingsManager.readGuiFilterSettings(0);
        String filename = (silent) ? filterSettings.getFilterSetFilename() : ImageJUtils.getFilename("Filter_range_file", filterSettings.getFilterSetFilename());
        if (filename == null) {
            return;
        }
        // Remove extension to store the filename
        filename = FileUtils.replaceExtension(filename, ".xml");
        filterSettings = filterSettings.toBuilder().setFilterSetFilename(filename).build();
        // Create a filter set using the ranges
        final ArrayList<Filter> filters = new ArrayList<>(4);
        // Create the multi-filter using the same precision type as that used during fitting.
        // Currently no support for z-filter as 3D astigmatism fitting is experimental.
        final PrecisionMethod precisionMethod = getPrecisionMethod((DirectFilter) multiFilter.getFilter());
        Function<double[], Filter> generator;
        if (precisionMethod == PrecisionMethod.POISSON_CRLB) {
            generator = parameters -> new MultiFilterCrlb(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else if (precisionMethod == PrecisionMethod.MORTENSEN) {
            generator = parameters -> new MultiFilter(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        } else {
            // Default
            generator = parameters -> new MultiFilter2(parameters[FILTER_SIGNAL], (float) parameters[FILTER_SNR], parameters[FILTER_MIN_WIDTH], parameters[FILTER_MAX_WIDTH], parameters[FILTER_SHIFT], parameters[FILTER_ESHIFT], parameters[FILTER_PRECISION], 0f, 0f);
        }
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(increment));
        if (saveFilters(filename, filters)) {
            SettingsManager.writeSettings(filterSettings);
        }
        // Create a filter set using the min/max and the initial bounds.
        // Set sensible limits
        min[FILTER_SIGNAL] = Math.max(min[FILTER_SIGNAL], 30);
        max[FILTER_SNR] = Math.min(max[FILTER_SNR], 10000);
        max[FILTER_PRECISION] = Math.min(max[FILTER_PRECISION], 100);
        // Make the 4-set filters the same as the 3-set filters.
        filters.clear();
        filters.add(generator.apply(min));
        filters.add(generator.apply(lower));
        filters.add(generator.apply(upper));
        filters.add(generator.apply(max));
        saveFilters(FileUtils.replaceExtension(filename, ".4.xml"), filters);
    }
    spotFitResults.min = min;
    spotFitResults.max = max;
}
Also used : Color(java.awt.Color) PeakResultPoint(uk.ac.sussex.gdsc.smlm.results.PeakResultPoint) Arrays(java.util.Arrays) CoordinateStoreFactory(uk.ac.sussex.gdsc.smlm.results.filter.CoordinateStoreFactory) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) BasePreprocessedPeakResult(uk.ac.sussex.gdsc.smlm.results.filter.BasePreprocessedPeakResult) MultiPathFitResults(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResults) Filter(uk.ac.sussex.gdsc.smlm.results.filter.Filter) HelpUrls(uk.ac.sussex.gdsc.smlm.ij.plugins.HelpUrls) Pair(org.apache.commons.lang3.tuple.Pair) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults) FilterValidationFlag(uk.ac.sussex.gdsc.smlm.results.filter.FilterValidationFlag) FitProtosHelper(uk.ac.sussex.gdsc.smlm.data.config.FitProtosHelper) ImageJImageConverter(uk.ac.sussex.gdsc.smlm.ij.utils.ImageJImageConverter) PrecisionFilter(uk.ac.sussex.gdsc.smlm.results.filter.PrecisionFilter) WidthFilter2(uk.ac.sussex.gdsc.smlm.results.filter.WidthFilter2) LinearInterpolator(org.apache.commons.math3.analysis.interpolation.LinearInterpolator) BlockingQueue(java.util.concurrent.BlockingQueue) StopWatch(org.apache.commons.lang3.time.StopWatch) ConcurrencyUtils(uk.ac.sussex.gdsc.core.utils.concurrent.ConcurrencyUtils) MultiPathFitResult(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResult) TextUtils(uk.ac.sussex.gdsc.core.utils.TextUtils) Plot(ij.gui.Plot) PeakFit(uk.ac.sussex.gdsc.smlm.ij.plugins.PeakFit) TIntHashSet(gnu.trove.set.hash.TIntHashSet) ImagePlus(ij.ImagePlus) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue) TextArea(java.awt.TextArea) PeakFractionalAssignment(uk.ac.sussex.gdsc.smlm.results.filter.PeakFractionalAssignment) XmlUtils(uk.ac.sussex.gdsc.core.utils.XmlUtils) ShiftFilter(uk.ac.sussex.gdsc.smlm.results.filter.ShiftFilter) FileUtils(uk.ac.sussex.gdsc.core.utils.FileUtils) PlugIn(ij.plugin.PlugIn) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) PolynomialSplineFunction(org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction) Prefs(ij.Prefs) TIntProcedure(gnu.trove.procedure.TIntProcedure) ArrayList(java.util.ArrayList) SortUtils(uk.ac.sussex.gdsc.core.utils.SortUtils) SignalFilter(uk.ac.sussex.gdsc.smlm.results.filter.SignalFilter) FitConfiguration(uk.ac.sussex.gdsc.smlm.engine.FitConfiguration) BenchmarkSpotFilterResult(uk.ac.sussex.gdsc.smlm.ij.plugins.benchmark.BenchmarkSpotFilter.BenchmarkSpotFilterResult) FitEngineConfigurationProvider(uk.ac.sussex.gdsc.smlm.ij.plugins.PeakFit.FitEngineConfigurationProvider) Assignment(uk.ac.sussex.gdsc.core.match.Assignment) FitWorker(uk.ac.sussex.gdsc.smlm.engine.FitWorker) MultiFilterCrlb(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilterCrlb) FileOutputStream(java.io.FileOutputStream) FilterResult(uk.ac.sussex.gdsc.smlm.ij.plugins.benchmark.BenchmarkSpotFilter.FilterResult) FractionClassificationResult(uk.ac.sussex.gdsc.core.match.FractionClassificationResult) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) ResultsMatchCalculator(uk.ac.sussex.gdsc.smlm.ij.plugins.ResultsMatchCalculator) PreprocessedPeakResult(uk.ac.sussex.gdsc.smlm.results.filter.PreprocessedPeakResult) Percentile(org.apache.commons.math3.stat.descriptive.rank.Percentile) EShiftFilter(uk.ac.sussex.gdsc.smlm.results.filter.EShiftFilter) DescriptiveStatistics(org.apache.commons.math3.stat.descriptive.DescriptiveStatistics) ImageStack(ij.ImageStack) MaximaSpotFilter(uk.ac.sussex.gdsc.smlm.filters.MaximaSpotFilter) FitTask(uk.ac.sussex.gdsc.smlm.engine.FitParameters.FitTask) NullFailCounter(uk.ac.sussex.gdsc.smlm.results.count.NullFailCounter) FractionScoreStore(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter.FractionScoreStore) TIntObjectHashMap(gnu.trove.map.hash.TIntObjectHashMap) TextWindow(ij.text.TextWindow) Spot(uk.ac.sussex.gdsc.smlm.filters.Spot) ItemListener(java.awt.event.ItemListener) FitSolver(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.FitSolver) ParameterType(uk.ac.sussex.gdsc.smlm.results.filter.ParameterType) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) OutOfRangeException(org.apache.commons.math3.exception.OutOfRangeException) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) AssignmentComparator(uk.ac.sussex.gdsc.core.match.AssignmentComparator) PeakResults(uk.ac.sussex.gdsc.smlm.results.PeakResults) ResultAssignment(uk.ac.sussex.gdsc.smlm.results.filter.ResultAssignment) PsfCalculator(uk.ac.sussex.gdsc.smlm.ij.plugins.PsfCalculator) PlotWindow(ij.gui.PlotWindow) MathUtils(uk.ac.sussex.gdsc.core.utils.MathUtils) FitParameters(uk.ac.sussex.gdsc.smlm.engine.FitParameters) SettingsManager(uk.ac.sussex.gdsc.smlm.ij.settings.SettingsManager) ItemEvent(java.awt.event.ItemEvent) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) FitEngineConfiguration(uk.ac.sussex.gdsc.smlm.engine.FitEngineConfiguration) Coordinate(uk.ac.sussex.gdsc.core.match.Coordinate) FractionalAssignment(uk.ac.sussex.gdsc.core.match.FractionalAssignment) FilterXStreamUtils(uk.ac.sussex.gdsc.smlm.results.filter.FilterXStreamUtils) FitStatus(uk.ac.sussex.gdsc.smlm.fitting.FitStatus) List(java.util.List) PointPair(uk.ac.sussex.gdsc.core.match.PointPair) SimpleArrayUtils(uk.ac.sussex.gdsc.core.utils.SimpleArrayUtils) ParameterisedFitJob(uk.ac.sussex.gdsc.smlm.engine.ParameterisedFitJob) Rectangle(java.awt.Rectangle) FastCorrelator(uk.ac.sussex.gdsc.core.utils.FastCorrelator) PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod) ImmutableFractionalAssignment(uk.ac.sussex.gdsc.core.match.ImmutableFractionalAssignment) MultiFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter) WindowOrganiser(uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser) ScoredSpot(uk.ac.sussex.gdsc.smlm.ij.plugins.benchmark.BenchmarkSpotFilter.ScoredSpot) AtomicReference(java.util.concurrent.atomic.AtomicReference) Function(java.util.function.Function) TextField(java.awt.TextField) IJImageSource(uk.ac.sussex.gdsc.smlm.ij.IJImageSource) Correlator(uk.ac.sussex.gdsc.core.utils.Correlator) NoiseEstimatorMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.NoiseEstimatorMethod) LinkedList(java.util.LinkedList) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) WidthFilter(uk.ac.sussex.gdsc.smlm.results.filter.WidthFilter) GUIFilterSettings(uk.ac.sussex.gdsc.smlm.ij.settings.GUIProtos.GUIFilterSettings) RampedScore(uk.ac.sussex.gdsc.core.utils.RampedScore) FitResult(uk.ac.sussex.gdsc.smlm.fitting.FitResult) Checkbox(java.awt.Checkbox) FilterSet(uk.ac.sussex.gdsc.smlm.results.filter.FilterSet) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) SnrFilter(uk.ac.sussex.gdsc.smlm.results.filter.SnrFilter) MultiFilter2(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter2) SimpleRegression(org.apache.commons.math3.stat.regression.SimpleRegression) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) ImageJUtils(uk.ac.sussex.gdsc.core.ij.ImageJUtils) SynchronizedPeakResults(uk.ac.sussex.gdsc.smlm.results.SynchronizedPeakResults) IJ(ij.IJ) SmlmUsageTracker(uk.ac.sussex.gdsc.smlm.ij.plugins.SmlmUsageTracker) Collections(java.util.Collections) ArrayList(java.util.ArrayList) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) PrecisionMethod(uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod) MultiPathFitResult(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResult) MultiFilterCrlb(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilterCrlb) PeakFractionalAssignment(uk.ac.sussex.gdsc.smlm.results.filter.PeakFractionalAssignment) FractionalAssignment(uk.ac.sussex.gdsc.core.match.FractionalAssignment) ImmutableFractionalAssignment(uk.ac.sussex.gdsc.core.match.ImmutableFractionalAssignment) FractionClassificationResult(uk.ac.sussex.gdsc.core.match.FractionClassificationResult) BasePreprocessedPeakResult(uk.ac.sussex.gdsc.smlm.results.filter.BasePreprocessedPeakResult) PreprocessedPeakResult(uk.ac.sussex.gdsc.smlm.results.filter.PreprocessedPeakResult) ScoredSpot(uk.ac.sussex.gdsc.smlm.ij.plugins.benchmark.BenchmarkSpotFilter.ScoredSpot) FastCorrelator(uk.ac.sussex.gdsc.core.utils.FastCorrelator) Plot(ij.gui.Plot) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) StoredDataStatistics(uk.ac.sussex.gdsc.core.utils.StoredDataStatistics) Coordinate(uk.ac.sussex.gdsc.core.match.Coordinate) FitConfiguration(uk.ac.sussex.gdsc.smlm.engine.FitConfiguration) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) TIntHashSet(gnu.trove.set.hash.TIntHashSet) GUIFilterSettings(uk.ac.sussex.gdsc.smlm.ij.settings.GUIProtos.GUIFilterSettings) SignalFilter(uk.ac.sussex.gdsc.smlm.results.filter.SignalFilter) WindowOrganiser(uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser) MultiFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter) PeakResultPoint(uk.ac.sussex.gdsc.smlm.results.PeakResultPoint) BasePoint(uk.ac.sussex.gdsc.core.match.BasePoint) PeakFractionalAssignment(uk.ac.sussex.gdsc.smlm.results.filter.PeakFractionalAssignment) FractionScoreStore(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter.FractionScoreStore) SimpleRegression(org.apache.commons.math3.stat.regression.SimpleRegression) Filter(uk.ac.sussex.gdsc.smlm.results.filter.Filter) PrecisionFilter(uk.ac.sussex.gdsc.smlm.results.filter.PrecisionFilter) ShiftFilter(uk.ac.sussex.gdsc.smlm.results.filter.ShiftFilter) MultiPathFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFilter) SignalFilter(uk.ac.sussex.gdsc.smlm.results.filter.SignalFilter) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) EShiftFilter(uk.ac.sussex.gdsc.smlm.results.filter.EShiftFilter) MaximaSpotFilter(uk.ac.sussex.gdsc.smlm.filters.MaximaSpotFilter) MultiFilter(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter) WidthFilter(uk.ac.sussex.gdsc.smlm.results.filter.WidthFilter) SnrFilter(uk.ac.sussex.gdsc.smlm.results.filter.SnrFilter) MultiFilter2(uk.ac.sussex.gdsc.smlm.results.filter.MultiFilter2) MultiPathFitResults(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResults)

Aggregations

PrecisionMethod (uk.ac.sussex.gdsc.smlm.data.config.FitProtos.PrecisionMethod)4 Rectangle (java.awt.Rectangle)2 DescriptiveStatistics (org.apache.commons.math3.stat.descriptive.DescriptiveStatistics)2 ExtendedGenericDialog (uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog)2 TIntObjectHashMap (gnu.trove.map.hash.TIntObjectHashMap)1 TIntProcedure (gnu.trove.procedure.TIntProcedure)1 TIntHashSet (gnu.trove.set.hash.TIntHashSet)1 IJ (ij.IJ)1 ImagePlus (ij.ImagePlus)1 ImageStack (ij.ImageStack)1 Prefs (ij.Prefs)1 GenericDialog (ij.gui.GenericDialog)1 Plot (ij.gui.Plot)1 PlotWindow (ij.gui.PlotWindow)1 PlugIn (ij.plugin.PlugIn)1 TextWindow (ij.text.TextWindow)1 Checkbox (java.awt.Checkbox)1 Color (java.awt.Color)1 TextArea (java.awt.TextArea)1 TextField (java.awt.TextField)1