use of ij.gui.Plot 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 ij.gui.Plot in project GDSC-SMLM by aherbert.
the class SpotAnalysis method showProfile.
private void showProfile(String title, String yTitle, double[] xValues, double[] yValues, double[] yValues2) {
final Plot plot = new Plot(title, "Frame", yTitle);
plot.addPoints(xValues, yValues, Plot.LINE);
final double[] limits = MathUtils.limits(yValues);
plot.setLimits(xValues[0], xValues[xValues.length - 1], limits[0], limits[1]);
plot.draw();
plot.setColor(Color.red);
plot.addPoints(xValues, yValues2, Plot.LINE);
plot.setColor(Color.magenta);
// Add the on-frames
if (!onFrames.isEmpty()) {
final double[] onx = new double[onFrames.size()];
final double[] ony = new double[onx.length];
int count = 0;
for (final Spot s : onFrames) {
onx[count] = s.frame;
ony[count] = yValues[s.frame - 1];
count++;
}
plot.addPoints(onx, ony, Plot.CIRCLE);
}
// Add the candidate frames
if (!candidateFrames.isEmpty()) {
plot.setColor(Color.cyan);
final double[] onx = new double[candidateFrames.size()];
final double[] ony = new double[onx.length];
int count = 0;
for (int i = 0; i < candidateFrames.size(); i++) {
final int frame = candidateFrames.getQuick(i);
onx[count] = frame;
ony[count] = yValues[frame - 1];
count++;
}
plot.addPoints(onx, ony, Plot.BOX);
plot.setColor(Color.magenta);
}
// Overlay current position
plot.addPoints(new double[] { rawImp.getCurrentSlice(), rawImp.getCurrentSlice() }, limits, Plot.LINE);
plot.setColor(Color.blue);
ImageJUtils.display(title, plot);
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class DensityImage method computeRipleysPlot.
/**
* Compute the Ripley's L-function for user selected radii and show it on a plot.
*
* @param results the results
*/
private void computeRipleysPlot(MemoryPeakResults results) {
final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.addMessage("Compute Ripley's L(r) - r plot");
gd.addNumericField("Min_radius", settings.minR, 2);
gd.addNumericField("Max_radius", settings.maxR, 2);
gd.addNumericField("Increment", settings.incrementR, 2);
gd.addCheckbox("Confidence_intervals", settings.confidenceIntervals);
gd.showDialog();
if (gd.wasCanceled()) {
return;
}
settings.minR = gd.getNextNumber();
settings.maxR = gd.getNextNumber();
settings.incrementR = gd.getNextNumber();
settings.confidenceIntervals = gd.getNextBoolean();
if (settings.minR > settings.maxR || settings.incrementR < 0 || gd.invalidNumber()) {
IJ.error(TITLE, "Invalid radius parameters");
return;
}
final DensityManager dm = createDensityManager(results);
final double[][] values = calculateLScores(dm);
// 99% confidence intervals
final int iterations = (settings.confidenceIntervals) ? 99 : 0;
double[] upper = null;
double[] lower = null;
final Rectangle bounds = results.getBounds();
final double area = (double) bounds.width * bounds.height;
// Use a uniform distribution for the coordinates
final HaltonSequenceGenerator dist = new HaltonSequenceGenerator(2);
dist.skipTo(SeedFactory.createInt());
for (int i = 0; i < iterations; i++) {
IJ.showProgress(i, iterations);
IJ.showStatus(String.format("L-score confidence interval %d / %d", i + 1, iterations));
// Randomise coordinates
final float[] x = new float[results.size()];
final float[] y = new float[x.length];
for (int j = x.length; j-- > 0; ) {
final double[] d = dist.nextVector();
x[j] = (float) (d[0] * bounds.width);
y[j] = (float) (d[1] * bounds.height);
}
final double[][] values2 = calculateLScores(new DensityManager(x, y, area));
if (upper == null || lower == null) {
upper = values2[1];
lower = upper.clone();
} else {
for (int m = upper.length; m-- > 0; ) {
if (upper[m] < values2[1][m]) {
upper[m] = values2[1][m];
}
if (lower[m] > values2[1][m]) {
lower[m] = values2[1][m];
}
}
}
}
final String title = results.getName() + " Ripley's (L(r) - r) / r";
final Plot plot = new Plot(title, "Radius", "(L(r) - r) / r");
plot.addPoints(values[0], values[1], Plot.LINE);
// Get the limits
double yMin = min(0, values[1]);
double yMax = max(0, values[1]);
if (iterations > 0) {
yMin = min(yMin, lower);
yMax = max(yMax, upper);
}
plot.setLimits(0, values[0][values[0].length - 1], yMin, yMax);
if (iterations > 0) {
plot.setColor(Color.BLUE);
plot.addPoints(values[0], upper, 1);
plot.setColor(Color.RED);
plot.addPoints(values[0], lower, 1);
plot.setColor(Color.BLACK);
}
ImageJUtils.display(title, plot);
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class DiffusionRateTest method showExample.
private void showExample(int totalSteps, double diffusionSigma, UniformRandomProvider rng) {
final MoleculeModel m = new MoleculeModel(0, new double[3]);
final float[] xValues = new float[totalSteps];
final float[] x = new float[totalSteps];
final float[] y = new float[totalSteps];
final DiffusionType diffusionType = CreateDataSettingsHelper.getDiffusionType(settings.getDiffusionType());
double[] axis;
if (diffusionType == DiffusionType.LINEAR_WALK) {
axis = nextVector(SamplerUtils.createNormalizedGaussianSampler(rng));
} else {
axis = null;
}
for (int j = 0; j < totalSteps; j++) {
if (diffusionType == DiffusionType.GRID_WALK) {
m.walk(diffusionSigma, rng);
} else if (diffusionType == DiffusionType.LINEAR_WALK) {
m.slide(diffusionSigma, axis, rng);
} else {
m.move(diffusionSigma, rng);
}
x[j] = (float) (m.getX());
y[j] = (float) (m.getY());
xValues[j] = (float) ((j + 1) / settings.getStepsPerSecond());
}
// Plot x and y coords on a timeline
final String title = TITLE + " example coordinates";
final Plot plot = new Plot(title, "Time (seconds)", "Distance (um)");
final float[] xUm = convertToUm(x);
final float[] yUm = convertToUm(y);
float[] limits = MathUtils.limits(xUm);
limits = MathUtils.limits(limits, yUm);
plot.setLimits(0, totalSteps / settings.getStepsPerSecond(), limits[0], limits[1]);
plot.setColor(Color.red);
plot.addPoints(xValues, xUm, Plot.LINE);
plot.setColor(Color.blue);
plot.addPoints(xValues, yUm, Plot.LINE);
ImageJUtils.display(title, plot);
// Scale up and draw 2D position
for (int j = 0; j < totalSteps; j++) {
x[j] *= pluginSettings.magnification;
y[j] *= pluginSettings.magnification;
}
final float[] limitsx = getLimits(x);
final float[] limitsy = getLimits(y);
int width = (int) (limitsx[1] - limitsx[0]);
int height = (int) (limitsy[1] - limitsy[0]);
// Ensure we draw something, even it is a simple dot at the centre for no diffusion
if (width == 0) {
width = (int) (32 * pluginSettings.magnification);
limitsx[0] = -width / 2.0f;
}
if (height == 0) {
height = (int) (32 * pluginSettings.magnification);
limitsy[0] = -height / 2.0f;
}
final ImageProcessor ip = new ByteProcessor(width, height);
// Adjust x and y using the minimum to centre
x[0] -= limitsx[0];
y[0] -= limitsy[0];
for (int j = 1; j < totalSteps; j++) {
// Adjust x and y using the minimum to centre
x[j] -= limitsx[0];
y[j] -= limitsy[0];
// Draw a line
ip.setColor(32 + (223 * j) / (totalSteps - 1));
ip.drawLine(round(x[j - 1]), round(y[j - 1]), round(x[j]), round(y[j]));
}
// Draw the final position
ip.putPixel(round(x[totalSteps - 1]), round(y[totalSteps - 1]), 255);
final ImagePlus imp = ImageJUtils.display(TITLE + " example", ip);
// Apply the fire lookup table
WindowManager.setTempCurrentImage(imp);
final LutLoader lut = new LutLoader();
lut.run("fire");
WindowManager.setTempCurrentImage(null);
}
use of ij.gui.Plot in project GDSC-SMLM by aherbert.
the class CameraModelAnalysis method execute.
/**
* Execute the analysis.
*/
private boolean execute() {
dirty = false;
final CameraModelAnalysisSettings settings = this.settings.build();
if (!(getGain(settings) > 0)) {
ImageJUtils.log(TITLE + "Error: No total gain");
return false;
}
if (!(settings.getPhotons() > 0)) {
ImageJUtils.log(TITLE + "Error: No photons");
return false;
}
// Avoid repeating the same analysis
if (settings.equals(lastSettings)) {
return true;
}
lastSettings = settings;
final IntHistogram h = getHistogram(settings);
// Build cumulative distribution
final double[][] cdf1 = cumulativeHistogram(h);
final double[] x1 = cdf1[0];
final double[] y1 = cdf1[1];
// Interpolate to 300 steps faster evaluation?
// Get likelihood function
final LikelihoodFunction f = getLikelihoodFunction(settings);
// Create likelihood cumulative distribution
final double[][] cdf2 = cumulativeDistribution(settings, cdf1, f);
// Compute Komolgorov distance
final double[] distanceAndValue = getDistance(cdf1, cdf2);
final double distance = distanceAndValue[0];
final double value = distanceAndValue[1];
final double area = distanceAndValue[2];
final double[] x2 = cdf2[0];
final double[] y2 = cdf2[1];
// Fill y1
int offset = 0;
while (x2[offset] < x1[0]) {
offset++;
}
final double[] y1b = new double[y2.length];
System.arraycopy(y1, 0, y1b, offset, y1.length);
Arrays.fill(y1b, offset + y1.length, y2.length, y1[y1.length - 1]);
// KolmogorovSmirnovTest
// n is the number of samples used to build the probability distribution.
final int n = (int) MathUtils.sum(h.histogramCounts);
// From KolmogorovSmirnovTest.kolmogorovSmirnovTest(RealDistribution distribution, double[]
// data, boolean exact):
// Returns the p-value associated with the null hypothesis that data is a sample from
// distribution.
// E.g. If p<0.05 then the null hypothesis is rejected and the data do not match the
// distribution.
double pvalue = Double.NaN;
try {
pvalue = 1d - kolmogorovSmirnovTest.cdf(distance, n);
} catch (final MathArithmeticException ex) {
// Cannot be computed to leave at NaN
}
// Plot
final WindowOrganiser wo = new WindowOrganiser();
String title = TITLE + " CDF";
Plot plot = new Plot(title, "Count", "CDF");
final double max = 1.05 * MathUtils.maxDefault(1, y2);
plot.setLimits(x2[0], x2[x2.length - 1], 0, max);
plot.setColor(Color.blue);
plot.addPoints(x2, y1b, Plot.BAR);
plot.setColor(Color.red);
plot.addPoints(x2, y2, Plot.BAR);
plot.setColor(Color.magenta);
plot.drawLine(value, 0, value, max);
plot.setColor(Color.black);
plot.addLegend("CDF\nModel");
plot.addLabel(0, 0, String.format("Distance=%s @ %.0f (Mean=%s) : p=%s", MathUtils.rounded(distance), value, MathUtils.rounded(area), MathUtils.rounded(pvalue)));
ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
// Show the histogram
title = TITLE + " Histogram";
plot = new Plot(title, "Count", "Frequency");
plot.setLimits(x1[0] - 0.5, x1[x1.length - 1] + 1.5, 0, MathUtils.max(h.histogramCounts) * 1.05);
plot.setColor(Color.blue);
plot.addPoints(x1, SimpleArrayUtils.toDouble(h.histogramCounts), Plot.BAR);
plot.setColor(Color.red);
final double[] x = floatHistogram[0].clone();
final double[] y = floatHistogram[1].clone();
final double scale = n / (MathUtils.sum(y) * (x[1] - x[0]));
for (int i = 0; i < y.length; i++) {
y[i] *= scale;
}
plot.addPoints(x, y, Plot.LINE);
plot.setColor(Color.black);
plot.addLegend("Sample\nExpected");
ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
wo.tile();
return true;
}
Aggregations