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

Example 1 with StoredData

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

the class BenchmarkSpotFilter method histogramFailures.

/**
 * Histogram the number of negatives preceding each positive.
 *
 * @param benchmarkFilterResult the filter result
 * @return the cumulative histogram
 */
private static double[][] histogramFailures(BenchmarkSpotFilterResult benchmarkFilterResult) {
    final StoredData data = new StoredData();
    benchmarkFilterResult.filterResults.forEachEntry((TIntObjectProcedure<FilterResult>) (peak, filterResult) -> {
        for (final ScoredSpot spot : filterResult.spots) {
            if (spot.match) {
                data.add(spot.fails);
            }
        }
        return true;
    });
    final double[][] h = MathUtils.cumulativeHistogram(data.getValues(), true);
    benchmarkFilterResult.cumul = h;
    benchmarkFilterResult.stats = data;
    return h;
}

Example 2 with StoredData

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

the class PsfCreator method plotCumulativeSignal.

/**
 * Show a plot of the cumulative signal vs distance from the centre.
 *
 * @param z The slice to plot
 * @param normalise normalise the sum to 1
 * @param resetScale Reset the y-axis maximum
 * @param distanceThreshold The distance threshold for the cumulative total shown in the plot
 *        label
 */
private void plotCumulativeSignal(int z, boolean normalise, boolean resetScale, double distanceThreshold) {
    final float[] data = (float[]) psf.getProcessor(z).getPixels();
    final int size = psf.getWidth();
    if (indexLookup == null || indexLookup.length != data.length) {
        // Precompute square distances
        final double[] d2 = new double[size];
        for (int y = 0, y2 = -size / 2; y < size; y++, y2++) {
            d2[y] = y2 * y2;
        }
        // Precompute distances
        final double[] d = new double[data.length];
        for (int y = 0, i = 0; y < size; y++) {
            for (int x = 0; x < size; x++, i++) {
                d[i] = Math.sqrt(d2[y] + d2[x]);
            }
        }
        // Sort
        final int[] indices = SimpleArrayUtils.natural(d.length);
        SortUtils.sortData(indices, d, true, false);
        // Store a unique cumulative index for each distance
        double lastD = d[0];
        int lastI = 0;
        int counter = 0;
        final StoredData distance = new StoredData();
        indexLookup = new int[indices.length];
        for (int i = 0; i < indices.length; i++) {
            if (lastD != d[i]) {
                distance.add(lastD * psfNmPerPixel);
                for (int j = lastI; j < i; j++) {
                    indexLookup[indices[j]] = counter;
                }
                lastD = d[i];
                lastI = i;
                counter++;
            }
        }
        // Do the final distance
        distance.add(lastD * psfNmPerPixel);
        for (int j = lastI; j < indices.length; j++) {
            indexLookup[indices[j]] = counter;
        }
        counter++;
        distances = distance.getValues();
    }
    // Get the signal at each distance
    final double[] signal = new double[distances.length];
    for (int i = 0; i < data.length; i++) {
        if (data[i] > 0) {
            signal[indexLookup[i]] += data[i];
        }
    }
    // Get the cumulative signal
    for (int i = 1; i < signal.length; i++) {
        signal[i] += signal[i - 1];
    }
    // Get the total up to the distance threshold
    double sum = 0;
    for (int i = 0; i < signal.length; i++) {
        if (distances[i] > distanceThreshold) {
            break;
        }
        sum = signal[i];
    }
    if (normalise && distanceThreshold > 0) {
        for (int i = 0; i < signal.length; i++) {
            signal[i] /= sum;
        }
    }
    if (resetScale) {
        maxCumulativeSignal = 0;
    }
    maxCumulativeSignal = MathUtils.maxDefault(maxCumulativeSignal, signal);
    final String title = "Cumulative Signal";
    final boolean alignWindows = (WindowManager.getFrame(title) == null);
    final Plot plot = new Plot(title, "Distance (nm)", "Signal");
    plot.addPoints(distances, signal, Plot.LINE);
    plot.setLimits(0, distances[distances.length - 1], 0, maxCumulativeSignal);
    plot.addLabel(0, 0, String.format("Total = %s (@ %s nm). z = %s nm", MathUtils.rounded(sum), MathUtils.rounded(distanceThreshold), MathUtils.rounded((z - zCentre) * settings.getNmPerSlice())));
    plot.setColor(Color.green);
    plot.drawLine(distanceThreshold, 0, distanceThreshold, maxCumulativeSignal);
    plot.setColor(Color.blue);
    final PlotWindow plotWindow = ImageJUtils.display(title, plot);
    if (alignWindows && plotWindow != null) {
        final PlotWindow otherWindow = getPlot(TITLE_PSF_PARAMETERS);
        if (otherWindow != null) {
            // Put the two plots tiled together so both are visible
            final Point l = plotWindow.getLocation();
            l.x = otherWindow.getLocation().x + otherWindow.getWidth();
            l.y = otherWindow.getLocation().y + otherWindow.getHeight();
            plotWindow.setLocation(l);
        }
    }
    // Update the PSF to the correct slice
    if (psfImp != null) {
        psfImp.setSlice(z);
    }
}

Example 3 with StoredData

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

the class DarkTimeAnalysis method analyse.

private void analyse(MemoryPeakResults results) {
    // Find min and max time frames
    results.sort();
    final int min = results.getFirstFrame();
    final int max = results.getLastFrame();
    // Trace results:
    // TODO - The search distance could have units to avoid assuming the results are in pixels
    final double d = settings.searchDistance / results.getCalibrationReader().getNmPerPixel();
    int range = max - min + 1;
    if (settings.maxDarkTime > 0) {
        range = Math.max(1, (int) Math.round(settings.maxDarkTime * 1000 / msPerFrame));
    }
    final TrackProgress tracker = SimpleImageJTrackProgress.getInstance();
    tracker.status("Analysing ...");
    tracker.log("Analysing (d=%s nm (%s px) t=%s s (%d frames)) ...", MathUtils.rounded(settings.searchDistance), MathUtils.rounded(d), MathUtils.rounded(range * msPerFrame / 1000.0), range);
    Trace[] traces;
    if (settings.method == 0) {
        final TraceManager tm = new TraceManager(results);
        tm.setTracker(tracker);
        tm.traceMolecules(d, range);
        traces = tm.getTraces();
    } else {
        final ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), algorithms[settings.method - 1], tracker);
        final List<Cluster> clusters = engine.findClusters(TraceMolecules.convertToClusterPoints(results), d, range);
        traces = TraceMolecules.convertToTraces(results, clusters);
    }
    tracker.status("Computing histogram ...");
    // Build dark-time histogram
    final int[] times = new int[range];
    final StoredData stats = new StoredData();
    for (final Trace trace : traces) {
        if (trace.getBlinks() > 1) {
            for (final int t : trace.getOffTimes()) {
                times[t]++;
            }
            stats.add(trace.getOffTimes());
        }
    }
    plotDarkTimeHistogram(stats);
    // Cumulative histogram
    for (int i = 1; i < times.length; i++) {
        times[i] += times[i - 1];
    }
    final int total = times[times.length - 1];
    // Plot dark-time up to 100%
    double[] x = new double[range];
    double[] y = new double[range];
    int truncate = 0;
    for (int i = 0; i < x.length; i++) {
        x[i] = i * msPerFrame;
        if (times[i] == total) {
            // Final value at 100%
            y[i] = 100.0;
            truncate = i + 1;
            break;
        }
        y[i] = (100.0 * times[i]) / total;
    }
    if (truncate > 0) {
        x = Arrays.copyOf(x, truncate);
        y = Arrays.copyOf(y, truncate);
    }
    final String title = "Cumulative Dark-time";
    final Plot plot = new Plot(title, "Time (ms)", "Percentile");
    plot.addPoints(x, y, Plot.LINE);
    ImageJUtils.display(title, plot);
    // Report percentile
    for (int i = 0; i < y.length; i++) {
        if (y[i] >= settings.percentile) {
            ImageJUtils.log("Dark-time Percentile %.1f @ %s ms = %s s", settings.percentile, MathUtils.rounded(x[i]), MathUtils.rounded(x[i] / 1000));
            break;
        }
    }
    tracker.status("");
}

Example 4 with StoredData

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

the class TcPalmAnalysis method plotHistogram.

/**
 * Plot a histogram of the extracted statistic.
 *
 * @param show set to true to show the histogram
 * @param wo the window organiser
 * @param clusters the clusters
 * @param name the name
 * @param function the function to extract the plotted statistic
 * @param predicate the predicate to filter the stream of data
 * @param action the action to use on the histogram plot (to set non-standard options)
 */
private static void plotHistogram(boolean show, WindowOrganiser wo, LocalList<ClusterData> clusters, String name, ToDoubleFunction<ClusterData> function, DoublePredicate predicate, Consumer<HistogramPlotBuilder> action) {
    if (!show) {
        return;
    }
    final StoredData data = new StoredData(clusters.size());
    DoubleStream stream = clusters.stream().mapToDouble(function);
    if (predicate != null) {
        stream = stream.filter(predicate);
    }
    stream.forEach(data::add);
    final HistogramPlotBuilder builder = new HistogramPlotBuilder(TITLE, data, name);
    action.accept(builder);
    builder.show(wo);
}

Example 5 with StoredData

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

the class TrackPopulationAnalysis method run.

@Override
public void run(String arg) {
    SmlmUsageTracker.recordPlugin(this.getClass(), arg);
    if (MemoryPeakResults.isMemoryEmpty()) {
        IJ.error(TITLE, "No localisations in memory");
        return;
    }
    settings = Settings.load();
    // Saved by reference so just save now
    settings.save();
    // Read in multiple traced datasets
    // All datasets must have the same pixel pitch and exposure time
    // Get parameters
    // Convert datasets to tracks
    // For each track compute the 4 local track features using the configured window
    // 
    // Optional:
    // Fit a multi-variate Gaussian mixture model to the data
    // (using the configured number of components/populations)
    // Assign each point in the track using the model.
    // Smooth the assignments.
    // 
    // The alternative is to use the localisation category to assign populations.
    // 
    // Plot histograms of each track parameter, coloured by component
    final List<MemoryPeakResults> combinedResults = new LocalList<>();
    if (!showInputDialog(combinedResults)) {
        return;
    }
    final boolean hasCategory = showHasCategoryDialog(combinedResults);
    if (!showDialog(hasCategory)) {
        return;
    }
    ImageJUtils.log(TITLE + "...");
    final List<Trace> tracks = getTracks(combinedResults, settings.window, settings.minTrackLength);
    if (tracks.isEmpty()) {
        IJ.error(TITLE, "No tracks. Please check the input data and min track length setting.");
        return;
    }
    final Calibration cal = combinedResults.get(0).getCalibration();
    final CalibrationReader cr = new CalibrationReader(cal);
    // Use micrometer / second
    final TypeConverter<DistanceUnit> distanceConverter = cr.getDistanceConverter(DistanceUnit.UM);
    final double exposureTime = cr.getExposureTime() / 1000.0;
    final Pair<int[], double[][]> trackData = extractTrackData(tracks, distanceConverter, exposureTime, hasCategory);
    final double[][] data = trackData.getValue();
    // Histogram the raw data.
    final Array2DRowRealMatrix raw = new Array2DRowRealMatrix(data, false);
    final WindowOrganiser wo = new WindowOrganiser();
    // Store the histogram data for plotting the components
    final double[][] columns = new double[FEATURE_NAMES.length][];
    final double[][] limits = new double[FEATURE_NAMES.length][];
    // Get column data
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        columns[i] = raw.getColumn(i);
        if (i == FEATURE_D) {
            // Plot using a logarithmic scale
            SimpleArrayUtils.apply(columns[i], Math::log10);
        }
        limits[i] = MathUtils.limits(columns[i]);
    }
    // Compute histogram bins
    final int[] bins = new int[FEATURE_NAMES.length];
    if (settings.histogramBins > 0) {
        Arrays.fill(bins, settings.histogramBins);
    } else {
        for (int i = 0; i < FEATURE_NAMES.length; i++) {
            bins[i] = HistogramPlot.getBins(StoredData.create(columns[i]), BinMethod.FD);
        }
        // Use the maximum so all histograms look the same
        Arrays.fill(bins, MathUtils.max(bins));
    }
    // Compute plots
    final Plot[] plots = new Plot[FEATURE_NAMES.length];
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        final double[][] hist = HistogramPlot.calcHistogram(columns[i], limits[i][0], limits[i][1], bins[i]);
        plots[i] = new Plot(TITLE + " " + FEATURE_NAMES[i], getFeatureLabel(i, i == FEATURE_D), "Frequency");
        plots[i].addPoints(hist[0], hist[1], Plot.BAR);
        ImageJUtils.display(plots[i].getTitle(), plots[i], ImageJUtils.NO_TO_FRONT, wo);
    }
    wo.tile();
    // The component for each data point
    int[] component;
    // The number of components
    int numComponents;
    // Data used to fit the Gaussian mixture model
    double[][] fitData;
    // The fitted model
    MixtureMultivariateGaussianDistribution model;
    if (hasCategory) {
        // Use the category as the component.
        // No fit data and no output model
        fitData = null;
        model = null;
        // The component is stored at the end of the raw track data.
        final int end = data[0].length - 1;
        component = Arrays.stream(data).mapToInt(d -> (int) d[end]).toArray();
        numComponents = MathUtils.max(component) + 1;
        // In the EM algorithm the probability of each data point is computed and normalised to
        // sum to 1. The normalised probabilities are averaged to create the weights.
        // Note the probability of each data point uses the previous weight and the algorithm
        // iterates.
        // This is not a fitted model but the input model so use
        // zero weights to indicate no fitting was performed.
        final double[] weights = new double[numComponents];
        // Remove the trailing component to show the 'model' in a table.
        createModelTable(Arrays.stream(data).map(d -> Arrays.copyOf(d, end)).toArray(double[][]::new), weights, component);
    } else {
        // Multivariate Gaussian mixture EM
        // Provide option to not use the anomalous exponent in the population mix.
        int sortDimension = SORT_DIMENSION;
        if (settings.ignoreAlpha) {
            // Remove index 0. This shifts the sort dimension.
            sortDimension--;
            fitData = Arrays.stream(data).map(d -> Arrays.copyOfRange(d, 1, d.length)).toArray(double[][]::new);
        } else {
            fitData = SimpleArrayUtils.deepCopy(data);
        }
        final MultivariateGaussianMixtureExpectationMaximization mixed = fitGaussianMixture(fitData, sortDimension);
        if (mixed == null) {
            IJ.error(TITLE, "Failed to fit a mixture model");
            return;
        }
        model = sortComponents(mixed.getFittedModel(), sortDimension);
        // For the best model, assign to the most likely population.
        component = assignData(fitData, model);
        // Table of the final model using the original data (i.e. not normalised)
        final double[] weights = model.getWeights();
        numComponents = weights.length;
        createModelTable(data, weights, component);
    }
    // Output coloured histograms of the populations.
    final LUT lut = LutHelper.createLut(settings.lutIndex);
    IntFunction<Color> colourMap;
    if (LutHelper.getColour(lut, 0).equals(Color.BLACK)) {
        colourMap = i -> LutHelper.getNonZeroColour(lut, i, 0, numComponents - 1);
    } else {
        colourMap = i -> LutHelper.getColour(lut, i, 0, numComponents - 1);
    }
    for (int i = 0; i < FEATURE_NAMES.length; i++) {
        // Extract the data for each component
        final double[] col = columns[i];
        final Plot plot = plots[i];
        for (int n = 0; n < numComponents; n++) {
            final StoredData feature = new StoredData();
            for (int j = 0; j < component.length; j++) {
                if (component[j] == n) {
                    feature.add(col[j]);
                }
            }
            if (feature.size() == 0) {
                continue;
            }
            final double[][] hist = HistogramPlot.calcHistogram(feature.values(), limits[i][0], limits[i][1], bins[i]);
            // Colour the points
            plot.setColor(colourMap.apply(n));
            plot.addPoints(hist[0], hist[1], Plot.BAR);
        }
        plot.updateImage();
    }
    createTrackDataTable(tracks, trackData, fitData, model, component, cal, colourMap);
// Analysis.
// Assign the original localisations to their track component.
// Q. What about the start/end not covered by the window?
// Save tracks as a dataset labelled with the sub-track ID?
// Output for the bound component and free components track parameters.
// Compute dwell times.
// Other ...
// Track analysis plugin:
// Extract all continuous segments of the same component.
// Produce MSD plot with error bars.
// Fit using FBM model.
}