use of uk.ac.sussex.gdsc.smlm.results.Trace in project GDSC-SMLM by aherbert.
the class TraceDiffusion method calculatePrecision.
/**
* Calculate the average precision of localisation in the traces.
*
* @param traces the traces
* @param multi true if the traces were from multiple input results
*/
private void calculatePrecision(Trace[] traces, boolean multi) {
// Check the diffusion simulation for a precision
if (DiffusionRateTest.isSimulated(results.getName()) && !multi) {
precision = DiffusionRateTest.getLastSimulationPrecision();
} else {
precision = 999;
try {
final Gaussian2DPeakResultCalculator calculator = Gaussian2DPeakResultHelper.create(results.getPsf(), results.getCalibration(), Gaussian2DPeakResultHelper.LSE_PRECISION);
// Get the average precision of the localisations
precision = 0;
int count = 0;
for (final Trace trace : traces) {
for (int k = 0; k < trace.size(); k++) {
final PeakResult r = trace.get(k);
precision += calculator.getLsePrecision(r.getParameters(), r.getNoise());
}
count += trace.size();
}
precision /= count;
} catch (final ConfigurationException ex) {
// Ignore this and we will ask the user for the precision
}
}
if (precision > 100) {
final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("The average precision of the traced results is " + MathUtils.rounded(precision, 4) + " nm.\nPlease verify the precision.");
gd.addSlider("Precision (nm)", 5, 100, precision);
gd.showDialog();
if (!(gd.wasCanceled() || gd.invalidNumber())) {
precision = Math.abs(gd.getNextNumber());
}
}
}
use of uk.ac.sussex.gdsc.smlm.results.Trace 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);
}
use of uk.ac.sussex.gdsc.smlm.results.Trace in project GDSC-SMLM by aherbert.
the class TraceDiffusion method createTraceFunction.
/**
* Creates the trace function for the configured trace diffusion mode.
*
* @return the function
*/
private Function<MemoryPeakResults, Trace[]> createTraceFunction() {
if (clusteringSettings.getTraceDiffusionMode() == 1) {
final DmttConfiguration config = createDmttConfiguration();
return r -> new DynamicMultipleTargetTracing(r).traceMolecules(config).toArray(new Trace[0]);
}
// Nearest neighbour
// Convert from NM to the native units of the results
final Converter c = CalibrationHelper.getDistanceConverter(results.getCalibration(), DistanceUnit.NM);
final double distanceThreshold = c.convertBack(clusteringSettings.getDistanceThreshold());
final double distanceExclusion = c.convertBack(clusteringSettings.getDistanceExclusion());
return r -> {
final TraceManager manager = new TraceManager(r);
// Run the tracing
manager.setTracker(SimpleImageJTrackProgress.getInstance());
manager.setDistanceExclusion(distanceExclusion);
manager.traceMolecules(distanceThreshold, 1);
return manager.getTraces();
};
}
use of uk.ac.sussex.gdsc.smlm.results.Trace in project GDSC-SMLM by aherbert.
the class TraceMolecules method convertToTraces.
/**
* Convert the clusters from the clustering engine into traces composed of the original list of
* peak results.
*
* @param results the results
* @param clusters the clusters
* @return the traces
*/
public static Trace[] convertToTraces(MemoryPeakResults results, List<Cluster> clusters) {
final Trace[] traces = new Trace[clusters.size()];
int index = 0;
for (final Cluster cluster : clusters) {
final Trace trace = new Trace();
trace.setId(index + 1);
for (ClusterPoint point = cluster.getHeadClusterPoint(); point != null; point = point.getNext()) {
// The point Id was the position in the original results array
trace.add(results.get(point.getId()));
}
trace.sort();
traces[index++] = trace;
}
return traces;
}
use of uk.ac.sussex.gdsc.smlm.results.Trace 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);
}
Aggregations