Examples with PrecisionResultProcedure uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure used on opensource projects

Example 6 with PrecisionResultProcedure

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

the class TraceMolecules method runOptimiser.

private void runOptimiser(TraceManager manager) {
    // Get an estimate of the number of molecules without blinking
    final Statistics stats = new Statistics();
    final double nmPerPixel = this.results.getNmPerPixel();
    final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
    pp.getPrecision();
    stats.add(pp.precisions);
    // Use twice the precision to get the initial distance threshold
    // Use 2.5x sigma as per the PC-PALM protocol in Sengupta, et al (2013) Nature Protocols 8, 345
    final double dEstimate = stats.getMean() * 2.5 / nmPerPixel;
    final int traceCount = manager.traceMolecules(dEstimate, 1);
    if (!getParameters(traceCount, dEstimate)) {
        return;
    }
    // TODO - Convert the distance threshold to use nm instead of pixels?
    final List<double[]> results = runTracing(manager, settings.getMinDistanceThreshold(), settings.getMaxDistanceThreshold(), settings.getMinTimeThreshold(), settings.getMaxTimeThreshold(), settings.getOptimiserSteps());
    // Compute fractional difference from the true value:
    // Use blinking rate directly or the estimated number of molecules
    double reference;
    int statistic;
    if (optimiseBlinkingRate) {
        reference = settings.getBlinkingRate();
        statistic = 3;
        IJ.log(String.format("Estimating blinking rate: %.2f", reference));
    } else {
        reference = traceCount / settings.getBlinkingRate();
        statistic = 2;
        IJ.log(String.format("Estimating number of molecules: %d / %.2f = %.2f", traceCount, settings.getBlinkingRate(), reference));
    }
    for (final double[] result : results) {
        if (optimiseBlinkingRate) {
            result[2] = (reference - result[statistic]) / reference;
        } else {
            result[2] = (result[statistic] - reference) / reference;
        }
    }
    // Locate the optimal parameters with a fit of the zero contour
    final boolean found = findOptimalParameters(results);
    createPlotResults(results);
    if (!found) {
        return;
    }
    // Make fractional difference absolute so that lowest is best
    for (final double[] result : results) {
        result[2] = Math.abs(result[2]);
    }
    // Set the optimal thresholds using the lowest value
    double[] best = new double[] { 0, 0, Double.MAX_VALUE };
    for (final double[] result : results) {
        if (best[2] > result[2]) {
            best = result;
        }
    }
    settings.setDistanceThreshold(best[0]);
    // The optimiser works using frames so convert back to the correct units
    final TypeConverter<TimeUnit> convert = UnitConverterUtils.createConverter(TimeUnit.FRAME, settings.getTimeUnit(), getExposureTimeInMilliSeconds());
    settings.setTimeThreshold(convert.convert(best[1]));
    IJ.log(String.format("Optimal fractional difference @ D-threshold=%g nm, T-threshold=%f s (%d frames)", settings.getDistanceThreshold(), timeThresholdInSeconds(), timeThresholdInFrames()));
    writeSettings();
}

Example 7 with PrecisionResultProcedure

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

the class TrackPopulationAnalysis method getTracks.

/**
 * Gets the tracks. Each track has contiguous frames and the length is enough to fit
 * {@code minTrackLength} overlapping windows of the specified size:
 *
 * <pre>
 * length >= window + minTrackLength - 1
 * </pre>
 *
 * @param combinedResults the combined results
 * @param window the window size
 * @param minTrackLength the minimum track length (assumed to be {@code >= 1})
 * @return the tracks
 */
private static List<Trace> getTracks(List<MemoryPeakResults> combinedResults, int window, int minTrackLength) {
    final LocalList<Trace> tracks = new LocalList<>();
    final Statistics stats = new Statistics();
    final int minSize = window + Math.max(minTrackLength, 1) - 1;
    combinedResults.forEach(results -> {
        final int start = tracks.size();
        // Sort by id then frame
        results = results.copy();
        results.sort(IdFramePeakResultComparator.INSTANCE);
        final int size = results.size();
        // Skip IDs not associated with clustering
        int index = 0;
        while (index < size && results.get(index).getId() < 1) {
            index++;
        }
        // Initialise current id and frame
        int id = results.get(index).getId() - 1;
        int frame = results.get(index).getFrame();
        Trace track = new Trace();
        for (; index < size; index++) {
            final PeakResult result = results.get(index);
            // Same ID and contiguous frames
            if (result.getId() != id || result.getFrame() != frame + 1) {
                addTrack(minSize, tracks, track);
                track = new Trace();
            }
            id = result.getId();
            frame = result.getFrame();
            track.add(result);
        }
        addTrack(minSize, tracks, track);
        stats.reset();
        for (int i = start; i < tracks.size(); i++) {
            stats.add(tracks.unsafeGet(i).size());
        }
        final StringBuilder sb = new StringBuilder(256);
        TextUtils.formatTo(sb, "%s tracks=%d, length=%s +/- %s", results.getName(), stats.getN(), MathUtils.rounded(stats.getMean(), 3), MathUtils.rounded(stats.getStandardDeviation(), 3));
        // Limit of diffusion coefficient from the localisation precision.
        // Just use the entire dataset for simplicity (i.e. not the tracks of min length).
        final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
        try {
            pp.getPrecision();
            final Mean mean = new Mean();
            for (final double p : pp.precisions) {
                mean.add(p);
            }
            // 2nDt = MSD (n = number of dimensions)
            // D = MSD / 2nt
            final CalibrationReader reader = results.getCalibrationReader();
            final double t = reader.getExposureTime() / 1000.0;
            // Assume computed in nm. Convert to um.
            final double x = mean.getMean() / 1000;
            final double d = x * x / (2 * t);
            TextUtils.formatTo(sb, ", precision=%s nm, D limit=%s um^2/s", MathUtils.rounded(x * 1000, 4), MathUtils.rounded(d, 4));
        } catch (final DataException ex) {
        // No precision
        }
        IJ.log(sb.toString());
    });
    return tracks;
}

Example 8 with PrecisionResultProcedure

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

the class SpotInspector method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    results = ResultsManager.loadInputResults(settings.inputOption, false, DistanceUnit.PIXEL, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        IJ.showStatus("");
        return;
    }
    // Check if the original image is open
    ImageSource source = results.getSource();
    if (source == null) {
        IJ.error(TITLE, "Unknown original source image");
        return;
    }
    source = source.getOriginal();
    source.setReadHint(ReadHint.NONSEQUENTIAL);
    if (!source.open()) {
        IJ.error(TITLE, "Cannot open original source image: " + source.toString());
        return;
    }
    final float stdDevMax = getStandardDeviation(results);
    if (stdDevMax < 0) {
        // TODO - Add dialog to get the initial peak width
        IJ.error(TITLE, "Fitting configuration (for initial peak width) is not available");
        return;
    }
    // Rank spots
    rankedResults = new ArrayList<>(results.size());
    // Data for the sorting
    final PrecisionResultProcedure pp;
    if (settings.sortOrderIndex == 1) {
        pp = new PrecisionResultProcedure(results);
        pp.getPrecision();
    } else {
        pp = null;
    }
    // Build procedures to get:
    // Shift = position in pixels - originXY
    final StandardResultProcedure sp;
    if (settings.sortOrderIndex == 9) {
        sp = new StandardResultProcedure(results, DistanceUnit.PIXEL);
        sp.getXyr();
    } else {
        sp = null;
    }
    // SD = gaussian widths only for Gaussian PSFs
    final WidthResultProcedure wp;
    if (settings.sortOrderIndex >= 6 && settings.sortOrderIndex <= 8) {
        wp = new WidthResultProcedure(results, DistanceUnit.PIXEL);
        wp.getWxWy();
    } else {
        wp = null;
    }
    // Amplitude for Gaussian PSFs
    final HeightResultProcedure hp;
    if (settings.sortOrderIndex == 2) {
        hp = new HeightResultProcedure(results, IntensityUnit.PHOTON);
        hp.getH();
    } else {
        hp = null;
    }
    final Counter c = new Counter();
    results.forEach((PeakResultProcedure) result -> {
        final float[] score = getScore(result, c.getAndIncrement(), pp, sp, wp, hp, stdDevMax);
        rankedResults.add(new PeakResultRank(result, score[0], score[1]));
    });
    Collections.sort(rankedResults, PeakResultRank::compare);
    // Prepare results table
    final ImageJTablePeakResults table = new ImageJTablePeakResults(false, results.getName(), true);
    table.copySettings(results);
    table.setTableTitle(TITLE);
    table.setAddCounter(true);
    table.setShowZ(results.is3D());
    // TODO - Add to settings
    table.setShowFittingData(true);
    table.setShowNoiseData(true);
    if (settings.showCalibratedValues) {
        table.setDistanceUnit(DistanceUnit.NM);
        table.setIntensityUnit(IntensityUnit.PHOTON);
    } else {
        table.setDistanceUnit(DistanceUnit.PIXEL);
        table.setIntensityUnit(IntensityUnit.COUNT);
    }
    table.begin();
    // Add a mouse listener to jump to the frame for the clicked line
    textPanel = table.getResultsWindow().getTextPanel();
    // We must ignore old instances of this class from the mouse listeners
    id = currentId.incrementAndGet();
    textPanel.addMouseListener(new MouseAdapter() {

        @Override
        public void mouseClicked(MouseEvent event) {
            SpotInspector.this.mouseClicked(event);
        }
    });
    // Add results to the table
    int count = 0;
    for (final PeakResultRank rank : rankedResults) {
        rank.rank = count++;
        table.add(rank.peakResult);
    }
    table.end();
    if (settings.plotScore || settings.plotHistogram) {
        // Get values for the plots
        float[] xValues = null;
        float[] yValues = null;
        double yMin;
        double yMax;
        int spotNumber = 0;
        xValues = new float[rankedResults.size()];
        yValues = new float[xValues.length];
        for (final PeakResultRank rank : rankedResults) {
            xValues[spotNumber] = spotNumber + 1;
            yValues[spotNumber++] = recoverScore(rank.score);
        }
        // Set the min and max y-values using 1.5 x IQR
        final DescriptiveStatistics stats = new DescriptiveStatistics();
        for (final float v : yValues) {
            stats.addValue(v);
        }
        if (settings.removeOutliers) {
            final double lower = stats.getPercentile(25);
            final double upper = stats.getPercentile(75);
            final double iqr = upper - lower;
            yMin = Math.max(lower - iqr, stats.getMin());
            yMax = Math.min(upper + iqr, stats.getMax());
            IJ.log(String.format("Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax));
        } else {
            yMin = stats.getMin();
            yMax = stats.getMax();
            IJ.log(String.format("Data range: %f - %f", yMin, yMax));
        }
        plotScore(xValues, yValues, yMin, yMax);
        plotHistogram(yValues);
    }
    // Extract spots into a stack
    final int w = source.getWidth();
    final int h = source.getHeight();
    final int size = 2 * settings.radius + 1;
    final ImageStack spots = new ImageStack(size, size, rankedResults.size());
    // To assist the extraction of data from the image source, process them in time order to allow
    // frame caching. Then set the appropriate slice in the result stack
    Collections.sort(rankedResults, (o1, o2) -> Integer.compare(o1.peakResult.getFrame(), o2.peakResult.getFrame()));
    for (final PeakResultRank rank : rankedResults) {
        final PeakResult r = rank.peakResult;
        // Extract image
        // Note that the coordinates are relative to the middle of the pixel (0.5 offset)
        // so do not round but simply convert to int
        final int x = (int) (r.getXPosition());
        final int y = (int) (r.getYPosition());
        // Extract a region but crop to the image bounds
        int minX = x - settings.radius;
        int minY = y - settings.radius;
        final int maxX = Math.min(x + settings.radius + 1, w);
        final int maxY = Math.min(y + settings.radius + 1, h);
        int padX = 0;
        int padY = 0;
        if (minX < 0) {
            padX = -minX;
            minX = 0;
        }
        if (minY < 0) {
            padY = -minY;
            minY = 0;
        }
        final int sizeX = maxX - minX;
        final int sizeY = maxY - minY;
        float[] data = source.get(r.getFrame(), new Rectangle(minX, minY, sizeX, sizeY));
        // Prevent errors with missing data
        if (data == null) {
            data = new float[sizeX * sizeY];
        }
        ImageProcessor spotIp = new FloatProcessor(sizeX, sizeY, data, null);
        // Pad if necessary, i.e. the crop is too small for the stack
        if (padX > 0 || padY > 0 || sizeX < size || sizeY < size) {
            final ImageProcessor spotIp2 = spotIp.createProcessor(size, size);
            spotIp2.insert(spotIp, padX, padY);
            spotIp = spotIp2;
        }
        final int slice = rank.rank + 1;
        spots.setPixels(spotIp.getPixels(), slice);
        spots.setSliceLabel(MathUtils.rounded(rank.originalScore), slice);
    }
    source.close();
    // Reset to the rank order
    Collections.sort(rankedResults, PeakResultRank::compare);
    final ImagePlus imp = ImageJUtils.display(TITLE, spots);
    imp.setRoi((PointRoi) null);
    // Make bigger
    for (int i = 10; i-- > 0; ) {
        imp.getWindow().getCanvas().zoomIn(imp.getWidth() / 2, imp.getHeight() / 2);
    }
}

Example 9 with PrecisionResultProcedure

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

the class TraceLengthAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    if (!showDialog()) {
        return;
    }
    // Load the results
    MemoryPeakResults results = ResultsManager.loadInputResults(settings.inputOption, false, null, null);
    if (MemoryPeakResults.isEmpty(results)) {
        IJ.error(TITLE, "No results could be loaded");
        return;
    }
    try {
        distanceConverter = results.getDistanceConverter(DistanceUnit.UM);
        timeConverter = results.getTimeConverter(TimeUnit.SECOND);
    } catch (final Exception ex) {
        IJ.error(TITLE, "Cannot convert units to um or seconds: " + ex.getMessage());
        return;
    }
    // Get the localisation error (4s^2) in raw units^2
    double precision = 0;
    try {
        final PrecisionResultProcedure p = new PrecisionResultProcedure(results);
        p.getPrecision();
        // Precision in nm using the median
        precision = new Percentile().evaluate(p.precisions, 50);
        // Convert from nm to um to raw units
        final double rawPrecision = distanceConverter.convertBack(precision / 1e3);
        // Get the localisation error (4s^2) in units^2
        error = 4 * rawPrecision * rawPrecision;
    } catch (final Exception ex) {
        ImageJUtils.log(TITLE + " - Unable to compute precision: " + ex.getMessage());
    }
    // Analyse the track lengths
    results = results.copy();
    results.sort(IdFramePeakResultComparator.INSTANCE);
    // Ensure the first result triggers an id change
    lastid = results.getFirst().getId() - 1;
    results.forEach(this::processTrackLength);
    // For the final track
    store();
    msds = msdList.toArray();
    lengths = lengthList.toArray();
    ids = idList.toArray();
    final int[] limits = MathUtils.limits(lengths);
    h1 = new int[limits[1] + 1];
    h2 = new int[h1.length];
    x1 = SimpleArrayUtils.newArray(h1.length, 0, 1f);
    y1 = new float[x1.length];
    y2 = new float[x1.length];
    // Sort by MSD
    final int[] indices = SimpleArrayUtils.natural(msds.length);
    SortUtils.sortIndices(indices, msds, false);
    final double[] msds2 = msds.clone();
    final int[] lengths2 = lengths.clone();
    final int[] ids2 = ids.clone();
    for (int i = 0; i < indices.length; i++) {
        msds[i] = msds2[indices[i]];
        lengths[i] = lengths2[indices[i]];
        ids[i] = ids2[indices[i]];
    }
    // Interactive analysis
    final NonBlockingExtendedGenericDialog gd = new NonBlockingExtendedGenericDialog(TITLE);
    ImageJUtils.addMessage(gd, "Split traces into fixed or moving using the track diffusion coefficient (D).\n" + "Localisation error has been subtracted from jumps (%s nm).", MathUtils.rounded(precision));
    final Statistics s = Statistics.create(msds);
    final double av = s.getMean();
    final String msg = String.format("Average D per track = %s um^2/s", MathUtils.rounded(av));
    gd.addMessage(msg);
    // Histogram the diffusion coefficients
    final WindowOrganiser wo = new WindowOrganiser();
    final HistogramPlot histogramPlot = new HistogramPlotBuilder("Trace diffusion coefficient", StoredData.create(msds), "D (um^2/s)").setRemoveOutliersOption(1).setPlotLabel(msg).build();
    histogramPlot.show(wo);
    final double[] xvalues = histogramPlot.getPlotXValues();
    final double min = xvalues[0];
    final double max = xvalues[xvalues.length - 1];
    // see if we can build a nice slider range from the histogram limits
    if (max - min < 5) {
        // Because sliders are used when the range is <5 and floating point
        gd.addSlider("D_threshold", min, max, settings.msdThreshold);
    } else {
        gd.addNumericField("D_threshold", settings.msdThreshold, 2, 6, "um^2/s");
    }
    gd.addCheckbox("Normalise", settings.normalise);
    gd.addDialogListener((gd1, event) -> {
        settings.msdThreshold = gd1.getNextNumber();
        settings.normalise = gd1.getNextBoolean();
        update();
        return true;
    });
    if (ImageJUtils.isShowGenericDialog()) {
        draw(wo);
        wo.tile();
    }
    gd.setOKLabel("Save datasets");
    gd.setCancelLabel("Close");
    gd.addHelp(HelpUrls.getUrl("trace-length-analysis"));
    gd.showDialog();
    if (gd.wasCanceled()) {
        return;
    }
    // Sort by ID
    final PeakResult[] list = results.toArray();
    Arrays.sort(list, IdFramePeakResultComparator.INSTANCE);
    createResults(results, "Fixed", 0, lastIndex, list);
    createResults(results, "Moving", lastIndex, msds.length, list);
}

Example 10 with PrecisionResultProcedure

use of uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure 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();
}