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

Example 11 with Statistics

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

the class PeakResultsReader method readNStorm.

private MemoryPeakResults readNStorm() {
    final MemoryPeakResults results = createResults();
    results.setName(FileUtils.getName(filename));
    try (FileInputStream fis = new FileInputStream(filename);
        BufferedReader input = new BufferedReader(new UnicodeReader(fis, null))) {
        final ProgressReporter reporter = createProgressReporter(fis);
        String line;
        int errors = 0;
        // The single line header
        final String header = input.readLine();
        if (header == null) {
            throw new IOException("NStorm header missing");
        }
        // NStorm files added more column fields for later formats.
        // If the header contains 'Photons' then this can be used to determine the gain
        boolean readPhotons = header.contains("\tPhotons\t");
        while ((line = input.readLine()) != null) {
            if (line.isEmpty()) {
                continue;
            }
            final PeakResult result = createNStormResult(line, readPhotons);
            if (result != null) {
                results.add(result);
                // Just read the photons from the first 100
                if (readPhotons) {
                    readPhotons = results.size() < 100;
                }
            } else if (++errors >= 10) {
                break;
            }
            reporter.showProgress();
        }
    } catch (final IOException ex) {
        logError(ex);
    }
    // The following relationship holds when length == 1:
    // intensity = height * 2 * pi * sd0 * sd1 / pixel_pitch^2
    // => Pixel_pitch = sqrt(height * 2 * pi * sd0 * sd1 / intensity)
    // Try and create a calibration
    final Statistics pixelPitch = new Statistics();
    results.forEach(new PeakResultProcedureX() {

        static final double TWO_PI = 2 * Math.PI;

        @Override
        public boolean execute(PeakResult peakResult) {
            if (peakResult.getFrame() == peakResult.getEndFrame()) {
                final float height = peakResult.getOrigValue();
                final float intensity = peakResult.getParameter(PeakResult.INTENSITY);
                final float sd0 = peakResult.getParameter(INDEX_SX);
                final float sd1 = peakResult.getParameter(INDEX_SY);
                pixelPitch.add(Math.sqrt(height * TWO_PI * sd0 * sd1 / intensity));
                // Stop when we have enough for a good guess
                return (pixelPitch.getN() > 100);
            }
            return false;
        }
    });
    // Determine the gain using the photons column
    final Statistics gain = new Statistics();
    results.forEach((PeakResultProcedureX) peakResult -> {
        double photons = peakResult.getError();
        if (photons != 0) {
            peakResult.setError(0);
            gain.add(peakResult.getIntensity() / photons);
            return false;
        }
        return true;
    });
    // TODO - Support all the NSTORM formats: one-axis, two-axis, rotated, 3D.
    // Is this information in the header?
    // We could support setting the PSF as a Gaussian2D with one/two axis SD.
    // This would mean updating all the result params if it is a one axis PSF.
    // For now just record it as a 2 axis PSF.
    // Create a calibration
    calibration = new CalibrationWriter();
    // NSTORM data is in counts when the Photons column is present.
    // Q. Is it in counts when this column is not present?
    calibration.setIntensityUnit(IntensityUnit.COUNT);
    calibration.setDistanceUnit(DistanceUnit.NM);
    if (pixelPitch.getN() > 0) {
        final double nmPerPixel = pixelPitch.getMean();
        calibration.setNmPerPixel(nmPerPixel);
    }
    if (gain.getN() > 0 && gain.getStandardError() < 1e-3) {
        calibration.setCountPerPhoton(gain.getMean());
    }
    results.setCalibration(calibration.getCalibration());
    return results;
}

Example 12 with Statistics

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

the class TraceDiffusion method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    jumpDistanceParametersRef.set(null);
    extraOptions = ImageJUtils.isExtraOptions();
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    settings = Settings.load();
    // Saved by reference so just save now
    settings.save();
    final ArrayList<MemoryPeakResults> allResults = new ArrayList<>();
    // Option to pick multiple input datasets together using a list box.
    if ("multi".equals(arg) && !showMultiDialog(allResults)) {
        return;
    }
    // This shows the dialog for selecting trace options
    if (!showTraceDialog(allResults)) {
        return;
    }
    if (allResults.isEmpty()) {
        return;
    }
    ImageJUtils.log(TITLE + "...");
    // - Trace each single dataset (and store in memory)
    // - Combine trace results held in memory
    final Trace[] traces = getTraces(allResults);
    // This still allows a zero entry in the results table.
    if (traces.length > 0 && !showDialog()) {
        return;
    }
    final int count = traces.length;
    double[] fitMsdResult = null;
    int numberOfDataPoints = 0;
    double[][] jdParams = null;
    if (count > 0) {
        calculatePrecision(traces, allResults.size() > 1);
        // --- MSD Analysis ---
        // Conversion constants
        final double px2ToUm2 = MathUtils.pow2(results.getCalibrationReader().getNmPerPixel()) / 1e6;
        final double px2ToUm2PerSecond = px2ToUm2 / exposureTime;
        // Get the maximum trace length
        int length = clusteringSettings.getMinimumTraceLength();
        if (!clusteringSettings.getTruncate()) {
            for (final Trace trace : traces) {
                if (length < trace.size()) {
                    length = trace.size();
                }
            }
        }
        // Get the localisation error (4s^2) in um^2
        final double error = (clusteringSettings.getPrecisionCorrection()) ? 4 * precision * precision / 1e6 : 0;
        // Pre-calculate MSD correction factors. This accounts for the fact that the distance moved
        // in the start/end frames is reduced due to the averaging of the particle location over the
        // entire frame into a single point. The true MSD may be restored by applying a factor.
        // Note: These are used for the calculation of the diffusion coefficients per molecule and
        // the MSD passed to the Jump Distance analysis. However the error is not included in the
        // jump distance analysis so will be subtracted from the fitted D coefficients later.
        final double[] factors;
        if (clusteringSettings.getMsdCorrection()) {
            factors = new double[length];
            for (int t = 1; t < length; t++) {
                factors[t] = JumpDistanceAnalysis.getConversionfactor(t);
            }
        } else {
            factors = SimpleArrayUtils.newArray(length, 0.0, 1.0);
        }
        // Extract the mean-squared distance statistics
        final Statistics[] stats = new Statistics[length];
        for (int i = 0; i < stats.length; i++) {
            stats[i] = new Statistics();
        }
        final ArrayList<double[]> distances = (settings.saveTraceDistances || settings.displayTraceLength) ? new ArrayList<>(traces.length) : null;
        // Store all the jump distances at the specified interval
        final StoredDataStatistics jumpDistances = new StoredDataStatistics();
        final int jumpDistanceInterval = clusteringSettings.getJumpDistance();
        // Compute squared distances
        final StoredDataStatistics msdPerMoleculeAllVsAll = new StoredDataStatistics();
        final StoredDataStatistics msdPerMoleculeAdjacent = new StoredDataStatistics();
        for (final Trace trace : traces) {
            final PeakResultStoreList results = trace.getPoints();
            // Sum the MSD and the time
            final int traceLength = (clusteringSettings.getTruncate()) ? clusteringSettings.getMinimumTraceLength() : trace.size();
            // Get the mean for each time separation
            final double[] sumDistance = new double[traceLength + 1];
            final double[] sumTime = new double[sumDistance.length];
            // Do the distances to the origin (saving if necessary)
            final float x0 = results.get(0).getXPosition();
            final float y0 = results.get(0).getYPosition();
            if (distances != null) {
                final double[] msd = new double[traceLength - 1];
                for (int j = 1; j < traceLength; j++) {
                    final int t = j;
                    final double d = distance2(x0, y0, results.get(j));
                    msd[j - 1] = px2ToUm2 * d;
                    if (t == jumpDistanceInterval) {
                        jumpDistances.add(msd[j - 1]);
                    }
                    sumDistance[t] += d;
                    sumTime[t] += t;
                }
                distances.add(msd);
            } else {
                for (int j = 1; j < traceLength; j++) {
                    final int t = j;
                    final double d = distance2(x0, y0, results.get(j));
                    if (t == jumpDistanceInterval) {
                        jumpDistances.add(px2ToUm2 * d);
                    }
                    sumDistance[t] += d;
                    sumTime[t] += t;
                }
            }
            if (clusteringSettings.getInternalDistances()) {
                // Do the internal distances
                for (int i = 1; i < traceLength; i++) {
                    final float x = results.get(i).getXPosition();
                    final float y = results.get(i).getYPosition();
                    for (int j = i + 1; j < traceLength; j++) {
                        final int t = j - i;
                        final double d = distance2(x, y, results.get(j));
                        if (t == jumpDistanceInterval) {
                            jumpDistances.add(px2ToUm2 * d);
                        }
                        sumDistance[t] += d;
                        sumTime[t] += t;
                    }
                }
                // Add the average distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    // Note: (traceLength - t) == count
                    stats[t].add(sumDistance[t] / (traceLength - t));
                }
            } else {
                // Add the distance per time separation to the population
                for (int t = 1; t < traceLength; t++) {
                    stats[t].add(sumDistance[t]);
                }
            }
            // Fix this for the precision and MSD adjustment.
            // It may be necessary to:
            // - sum the raw distances for each time interval (this is sumDistance[t])
            // - subtract the precision error
            // - apply correction factor for the n-frames to get actual MSD
            // - sum the actual MSD
            double sumD = 0;
            final double sumD_adjacent = Math.max(0, sumDistance[1] - error) * factors[1];
            double sumT = 0;
            final double sumT_adjacent = sumTime[1];
            for (int t = 1; t < traceLength; t++) {
                sumD += Math.max(0, sumDistance[t] - error) * factors[t];
                sumT += sumTime[t];
            }
            // Calculate the average displacement for the trace (do not simply use the largest
            // time separation since this will miss moving molecules that end up at the origin)
            msdPerMoleculeAllVsAll.add(px2ToUm2PerSecond * sumD / sumT);
            msdPerMoleculeAdjacent.add(px2ToUm2PerSecond * sumD_adjacent / sumT_adjacent);
        }
        StoredDataStatistics dperMoleculeAllVsAll = null;
        StoredDataStatistics dperMoleculeAdjacent = null;
        if (settings.saveTraceDistances || (clusteringSettings.getShowHistograms() && settings.displayDHistogram)) {
            dperMoleculeAllVsAll = calculateDiffusionCoefficient(msdPerMoleculeAllVsAll);
            dperMoleculeAdjacent = calculateDiffusionCoefficient(msdPerMoleculeAdjacent);
        }
        if (settings.saveTraceDistances) {
            saveTraceDistances(traces.length, distances, msdPerMoleculeAllVsAll, msdPerMoleculeAdjacent, dperMoleculeAllVsAll, dperMoleculeAdjacent);
        }
        if (settings.displayTraceLength) {
            final StoredDataStatistics lengths = calculateTraceLengths(distances);
            showHistogram(lengths, "Trace length (um)");
        }
        if (settings.displayTraceSize) {
            final StoredDataStatistics sizes = calculateTraceSizes(traces);
            showHistogram(sizes, "Trace size", true);
        }
        // Plot the per-trace histogram of MSD and D
        if (clusteringSettings.getShowHistograms()) {
            if (settings.displayMsdHistogram) {
                showHistogram(msdPerMoleculeAllVsAll, "MSD/Molecule (all-vs-all)");
                showHistogram(msdPerMoleculeAdjacent, "MSD/Molecule (adjacent)");
            }
            if (settings.displayDHistogram) {
                showHistogram(dperMoleculeAllVsAll, "D/Molecule (all-vs-all)");
                showHistogram(dperMoleculeAdjacent, "D/Molecule (adjacent)");
            }
        }
        // Calculate the mean squared distance (MSD)
        final double[] x = new double[stats.length];
        final double[] y = new double[x.length];
        final double[] sd = new double[x.length];
        // Intercept is the 4s^2 (in um^2)
        y[0] = 4 * precision * precision / 1e6;
        for (int i = 1; i < stats.length; i++) {
            x[i] = i * exposureTime;
            y[i] = stats[i].getMean() * px2ToUm2;
            // sd[i] = stats[i].getStandardDeviation() * px2ToUm2;
            sd[i] = stats[i].getStandardError() * px2ToUm2;
        }
        final String title = TITLE + " MSD";
        final Plot plot = plotMsd(x, y, sd, title);
        // Fit the MSD using a linear fit
        fitMsdResult = fitMsd(x, y, title, plot);
        // Jump Distance analysis
        if (settings.saveRawData) {
            saveStatistics(jumpDistances, "Jump Distance", "Distance (um^2)", false);
        }
        // Calculate the cumulative jump-distance histogram
        final double[][] jdHistogram = JumpDistanceAnalysis.cumulativeHistogram(jumpDistances.getValues());
        // Always show the jump distance histogram
        jdTitle = TITLE + " Jump Distance";
        jdPlot = new Plot(jdTitle, "Distance (um^2)", "Cumulative Probability");
        jdPlot.addPoints(jdHistogram[0], jdHistogram[1], Plot.LINE);
        display(jdTitle, jdPlot);
        // Fit Jump Distance cumulative probability
        numberOfDataPoints = jumpDistances.getN();
        jdParams = fitJumpDistance(jumpDistances, jdHistogram);
        jumpDistanceParametersRef.set(jdParams);
    }
    summarise(traces, fitMsdResult, numberOfDataPoints, jdParams);
}

Example 13 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics 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 14 with Statistics

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

the class TrackPopulationAnalysis method extractTrackData.

/**
 * Extract the track data. This extracts different descriptors of the track using a rolling local
 * window.
 *
 * <p>Distances are converted to {@code unit} using the provided converter and time units are
 * converted from frame to seconds (s). The diffusion coefficients is in unit^2/s.
 *
 * <p>If categories are added they are remapped to be a natural sequence starting from 0.
 *
 * @param tracks the tracks
 * @param distanceConverter the distance converter
 * @param deltaT the time step of each frame in seconds (delta T)
 * @param hasCategory if true add the category from the result to the end of the data
 * @return the track data (track lengths and descriptors)
 */
private Pair<int[], double[][]> extractTrackData(List<Trace> tracks, TypeConverter<DistanceUnit> distanceConverter, double deltaT, boolean hasCategory) {
    final List<double[]> data = new LocalList<>(tracks.size());
    double[] x = new double[0];
    double[] y = new double[0];
    final int wm1 = settings.window - 1;
    final int valuesLength = hasCategory ? 5 : 4;
    final int mid = settings.window / 2;
    final MsdFitter msdFitter = new MsdFitter(settings, deltaT);
    final double significance = settings.significance;
    final int[] fitResult = new int[4];
    // Factor for the diffusion coefficient: 1/N * 1/(2*dimensions*deltaT) = 1 / 4Nt
    // with N the number of points to average.
    final double diffusionCoefficientFactor = 1.0 / (4 * wm1 * deltaT);
    // Used for the standard deviations
    final Statistics statsX = new Statistics();
    final Statistics statsY = new Statistics();
    final Ticker ticker = ImageJUtils.createTicker(tracks.size(), 1, "Computing track features...");
    // Collect the track categories. Later these are renumbered to a natural sequence from 0.
    final TIntHashSet catSet = new TIntHashSet();
    // final StoredDataStatistics statsAlpha = new StoredDataStatistics(tracks.size());
    // Process each track
    final TIntArrayList lengths = new TIntArrayList(tracks.size());
    for (final Trace track : tracks) {
        // Get xy coordinates
        final int size = track.size();
        if (x.length < size) {
            x = new double[size];
            y = new double[size];
        }
        for (int i = 0; i < size; i++) {
            final PeakResult peak = track.get(i);
            x[i] = distanceConverter.convert(peak.getXPosition());
            y[i] = distanceConverter.convert(peak.getYPosition());
        }
        final int smwm1 = size - wm1;
        final int previousSize = data.size();
        for (int k = 0; k < smwm1; k++) {
            final double[] values = new double[valuesLength];
            data.add(values);
            // middle of even sized windows is between two localisations.
            if (hasCategory) {
                final int cat = track.get(k + mid).getCategory();
                values[4] = cat;
                catSet.add(cat);
            }
            // First point in window = k
            // Last point in window = k + w - 1 = k + wm1
            final int end = k + wm1;
            // 1. Anomalous exponent.
            msdFitter.setData(x, y);
            try {
                msdFitter.fit(k, null);
                // statsAlpha.add(msdFitter.alpha);
                int fitIndex = msdFitter.positiveSlope ? 0 : 2;
                // If better then this is anomalous diffusion
                final double alpha = msdFitter.lvmSolution2.getPoint().getEntry(2);
                values[0] = alpha;
                if (msdFitter.pValue > significance) {
                    fitIndex++;
                }
                fitResult[fitIndex]++;
            // values[0] = msdFitter.alpha;
            // // Debug
            // if (
            // // msdFitter.pValue < 0.2
            // msdFitter.pValue < 0.2 && values[0] < 0
            // // !msdFitter.positiveSlope
            // ) {
            // final RealVector p = msdFitter.lvmSolution2.getPoint();
            // final String title = "anomalous exponent";
            // final Plot plot = new Plot(title, "time (s)", "MSD (um^2)");
            // final double[] t = SimpleArrayUtils.newArray(msdFitter.s.length, deltaT, deltaT);
            // plot.addLabel(0, 0, msdFitter.lvmSolution2.getPoint().toString() + " p="
            // + msdFitter.pValue + ". " + msdFitter.lvmSolution1.getPoint().toString());
            // plot.addPoints(t, msdFitter.s, Plot.CROSS);
            // plot.addPoints(t, msdFitter.model2.value(p).getFirst().toArray(), Plot.LINE);
            // plot.setColor(Color.BLUE);
            // plot.addPoints(t,
            // msdFitter.model1.value(msdFitter.lvmSolution1.getPoint()).getFirst().toArray(),
            // Plot.LINE);
            // plot.setColor(Color.RED);
            // final double[] yy = Arrays.stream(t).map(msdFitter.reg::predict).toArray();
            // plot.addPoints(t, yy, Plot.CIRCLE);
            // System.out.printf("%s : %s", msdFitter.lvmSolution2.getPoint(), values[0]);
            // ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT);
            // System.out.println();
            // }
            } catch (TooManyIterationsException | ConvergenceException ex) {
                if (settings.debug) {
                    ImageJUtils.log("Failed to fit anomalous exponent: " + ex.getMessage());
                }
                // Ignore this and leave as Brownian motion
                values[0] = 1.0;
            }
            // Referenced papers:
            // Hozé, N. H., D. (2017) Statistical methods for large ensembles of super-resolution
            // stochastic single particle trajectories in cell biology.
            // Annual Review of Statistics and Its Application 4, 189-223
            // 
            // Amitai, A., Seeber, A., Gasser, S. M. & Holcman, D. (2017) Visualization of Chromatin
            // Decompaction and Break Site Extrusion as Predicted by Statistical Polymer
            // Modeling of Single-Locus Trajectories. Cell reports 18, 1200-1214
            // 2. Effective diffusion coefficient (Hozé, eq 10).
            // This is the average squared jump distance between successive points
            // divided by 1 / (2 * dimensions * deltaT), i.e. 1 / 4t.
            double sum = 0;
            for (int i = k; i < end; i++) {
                sum += MathUtils.distance2(x[i], y[i], x[i + 1], y[i + 1]);
            }
            values[1] = sum * diffusionCoefficientFactor;
            // 3. Length of confinement (Amitai et al, eq 1).
            // Compute the average of the standard deviation of the position in each dimension.
            statsX.reset();
            statsY.reset();
            for (int i = k; i <= end; i++) {
                statsX.add(x[i]);
                statsY.add(y[i]);
            }
            values[2] = (statsX.getStandardDeviation() + statsY.getStandardDeviation()) / 2;
            // 4. Magnitude of drift vector (Hozé, eq 9).
            // Note: The drift field is given as the expected distance between successive points, i.e.
            // the average step. Since all track windows are the same length this is the same
            // as the distance between the first and last point divided by the number of points.
            // The drift field in each dimension is combined to create a drift norm, i.e. Euclidean
            // distance.
            values[3] = MathUtils.distance(x[k], y[k], x[end], y[end]) / wm1;
        }
        lengths.add(data.size() - previousSize);
        ticker.tick();
    }
    ImageJUtils.finished();
    if (settings.debug) {
        ImageJUtils.log("  +Slope, significant:   %d", fitResult[0]);
        ImageJUtils.log("  +Slope, insignificant: %d", fitResult[1]);
        ImageJUtils.log("  -Slope, significant:   %d", fitResult[2]);
        ImageJUtils.log("  -Slope, insignificant: %d", fitResult[3]);
    }
    ImageJUtils.log("Insignificant anomalous exponents: %d / %d", fitResult[1] + fitResult[3], data.size());
    // System.out.println(statsAlpha.getStatistics().toString());
    final double[][] trackData = data.toArray(new double[0][0]);
    if (hasCategory) {
        final int[] categories = catSet.toArray();
        Arrays.sort(categories);
        // Only remap if non-compact (i.e. not 0 to n)
        final int max = categories[categories.length - 1];
        if (categories[0] != 0 || max + 1 != categories.length) {
            final int[] remap = new int[max + 1];
            for (int i = 0; i < categories.length; i++) {
                remap[categories[i]] = i;
            }
            final int end = valuesLength - 1;
            for (final double[] values : trackData) {
                values[end] = remap[(int) values[end]];
            }
        }
    }
    return Pair.create(lengths.toArray(), trackData);
}

Example 15 with Statistics

use of uk.ac.sussex.gdsc.core.utils.Statistics 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;
}