use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.
the class FIRE method cropToRoi.
private MemoryPeakResults cropToRoi(MemoryPeakResults results) {
if (roiBounds == null)
return results;
// Adjust bounds relative to input results image
//Rectangle2D.Float bounds = results.getDataBounds();
Rectangle bounds = results.getBounds(true);
double xscale = roiImageWidth / bounds.width;
double yscale = roiImageHeight / bounds.height;
float minX = (float) (bounds.x + roiBounds.x / xscale);
float maxX = (float) (minX + roiBounds.width / xscale);
float minY = (float) (bounds.y + (roiBounds.y / yscale));
float maxY = (float) (minY + roiBounds.height / yscale);
// Create a new set of results within the bounds
MemoryPeakResults newResults = new MemoryPeakResults();
newResults.begin();
for (PeakResult peakResult : results.getResults()) {
float x = peakResult.params[Gaussian2DFunction.X_POSITION];
float y = peakResult.params[Gaussian2DFunction.Y_POSITION];
if (x < minX || x > maxX || y < minY || y > maxY)
continue;
newResults.add(peakResult);
}
newResults.end();
newResults.copySettings(results);
newResults.setBounds(new Rectangle((int) minX, (int) minY, (int) (maxX - minX), (int) (maxY - minY)));
return newResults;
}
use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.
the class FIRE method runQEstimation.
private void runQEstimation() {
IJ.showStatus(TITLE + " ...");
if (!showQEstimationInputDialog())
return;
MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
return;
}
if (results.getCalibration() == null) {
IJ.error(TITLE, "The results are not calibrated");
return;
}
results = cropToRoi(results);
if (results.size() < 2) {
IJ.error(TITLE, "No results within the crop region");
return;
}
initialise(results, null);
// We need localisation precision.
// Build a histogram of the localisation precision.
// Get the initial mean and SD and plot as a Gaussian.
PrecisionHistogram histogram = calculatePrecisionHistogram();
if (histogram == null) {
IJ.error(TITLE, "No localisation precision available.\n \nPlease choose " + PrecisionMethod.FIXED + " and enter a precision mean and SD.");
return;
}
StoredDataStatistics precision = histogram.precision;
//String name = results.getName();
double fourierImageScale = SCALE_VALUES[imageScaleIndex];
int imageSize = IMAGE_SIZE_VALUES[imageSizeIndex];
// Create the image and compute the numerator of FRC.
// Do not use the signal so results.size() is the number of localisations.
IJ.showStatus("Computing FRC curve ...");
FireImages images = createImages(fourierImageScale, imageSize, false);
// DEBUGGING - Save the two images to disk. Load the images into the Matlab
// code that calculates the Q-estimation and make this plugin match the functionality.
//IJ.save(new ImagePlus("i1", images.ip1), "/scratch/i1.tif");
//IJ.save(new ImagePlus("i2", images.ip2), "/scratch/i2.tif");
FRC frc = new FRC();
frc.progress = progress;
frc.setFourierMethod(fourierMethod);
frc.setSamplingMethod(samplingMethod);
frc.setPerimeterSamplingFactor(perimeterSamplingFactor);
FRCCurve frcCurve = frc.calculateFrcCurve(images.ip1, images.ip2, images.nmPerPixel);
if (frcCurve == null) {
IJ.error(TITLE, "Failed to compute FRC curve");
return;
}
IJ.showStatus("Running Q-estimation ...");
// Note:
// The method implemented here is based on Matlab code provided by Bernd Rieger.
// The idea is to compute the spurious correlation component of the FRC Numerator
// using an initial estimate of distribution of the localisation precision (assumed
// to be Gaussian). This component is the contribution of repeat localisations of
// the same molecule to the numerator and is modelled as an exponential decay
// (exp_decay). The component is scaled by the Q-value which
// is the average number of times a molecule is seen in addition to the first time.
// At large spatial frequencies the scaled component should match the numerator,
// i.e. at high resolution (low FIRE number) the numerator is made up of repeat
// localisations of the same molecule and not actual structure in the image.
// The best fit is where the numerator equals the scaled component, i.e. num / (q*exp_decay) == 1.
// The FRC Numerator is plotted and Q can be determined by
// adjusting Q and the precision mean and SD to maximise the cost function.
// This can be done interactively by the user with the effect on the FRC curve
// dynamically updated and displayed.
// Compute the scaled FRC numerator
double qNorm = (1 / frcCurve.mean1 + 1 / frcCurve.mean2);
double[] frcnum = new double[frcCurve.getSize()];
for (int i = 0; i < frcnum.length; i++) {
FRCCurveResult r = frcCurve.get(i);
frcnum[i] = qNorm * r.getNumerator() / r.getNumberOfSamples();
}
// Compute the spatial frequency and the region for curve fitting
double[] q = FRC.computeQ(frcCurve, false);
int low = 0, high = q.length;
while (high > 0 && q[high - 1] > maxQ) high--;
while (low < q.length && q[low] < minQ) low++;
// Require we fit at least 10% of the curve
if (high - low < q.length * 0.1) {
IJ.error(TITLE, "Not enough points for Q estimation");
return;
}
// Obtain initial estimate of Q plateau height and decay.
// This can be done by fitting the precision histogram and then fixing the mean and sigma.
// Or it can be done by allowing the precision to be sampled and the mean and sigma
// become parameters for fitting.
// Check if we can sample precision values
boolean sampleDecay = precision != null && FIRE.sampleDecay;
double[] exp_decay;
if (sampleDecay) {
// Random sample of precision values from the distribution is used to
// construct the decay curve
int[] sample = Random.sample(10000, precision.getN(), new Well19937c());
final double four_pi2 = 4 * Math.PI * Math.PI;
double[] pre = new double[q.length];
for (int i = 1; i < q.length; i++) pre[i] = -four_pi2 * q[i] * q[i];
// Sample
final int n = sample.length;
double[] hq = new double[n];
for (int j = 0; j < n; j++) {
// Scale to SR pixels
double s2 = precision.getValue(sample[j]) / images.nmPerPixel;
s2 *= s2;
for (int i = 1; i < q.length; i++) hq[i] += FastMath.exp(pre[i] * s2);
}
for (int i = 1; i < q.length; i++) hq[i] /= n;
exp_decay = new double[q.length];
exp_decay[0] = 1;
for (int i = 1; i < q.length; i++) {
double sinc_q = sinc(Math.PI * q[i]);
exp_decay[i] = sinc_q * sinc_q * hq[i];
}
} else {
// Note: The sigma mean and std should be in the units of super-resolution
// pixels so scale to SR pixels
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Smoothing
double[] smooth;
if (loessSmoothing) {
// Note: This computes the log then smooths it
double bandwidth = 0.1;
int robustness = 0;
double[] l = new double[exp_decay.length];
for (int i = 0; i < l.length; i++) {
// Original Matlab code computes the log for each array.
// This is equivalent to a single log on the fraction of the two.
// Perhaps the two log method is more numerically stable.
//l[i] = Math.log(Math.abs(frcnum[i])) - Math.log(exp_decay[i]);
l[i] = Math.log(Math.abs(frcnum[i] / exp_decay[i]));
}
try {
LoessInterpolator loess = new LoessInterpolator(bandwidth, robustness);
smooth = loess.smooth(q, l);
} catch (Exception e) {
IJ.error(TITLE, "LOESS smoothing failed");
return;
}
} else {
// Note: This smooths the curve before computing the log
double[] norm = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
norm[i] = frcnum[i] / exp_decay[i];
}
// Median window of 5 == radius of 2
MedianWindow mw = new MedianWindow(norm, 2);
smooth = new double[exp_decay.length];
for (int i = 0; i < norm.length; i++) {
smooth[i] = Math.log(Math.abs(mw.getMedian()));
mw.increment();
}
}
// Fit with quadratic to find the initial guess.
// Note: example Matlab code frc_Qcorrection7.m identifies regions of the
// smoothed log curve with low derivative and only fits those. The fit is
// used for the final estimate. Fitting a subset with low derivative is not
// implemented here since the initial estimate is subsequently optimised
// to maximise a cost function.
Quadratic curve = new Quadratic();
SimpleCurveFitter fit = SimpleCurveFitter.create(curve, new double[2]);
WeightedObservedPoints points = new WeightedObservedPoints();
for (int i = low; i < high; i++) points.add(q[i], smooth[i]);
double[] estimate = fit.fit(points.toList());
double qValue = FastMath.exp(estimate[0]);
//System.out.printf("Initial q-estimate = %s => %.3f\n", Arrays.toString(estimate), qValue);
// This could be made an option. Just use for debugging
boolean debug = false;
if (debug) {
// Plot the initial fit and the fit curve
double[] qScaled = FRC.computeQ(frcCurve, true);
double[] line = new double[q.length];
for (int i = 0; i < q.length; i++) line[i] = curve.value(q[i], estimate);
String title = TITLE + " Initial fit";
Plot2 plot = new Plot2(title, "Spatial Frequency (nm^-1)", "FRC Numerator");
String label = String.format("Q = %.3f", qValue);
plot.addPoints(qScaled, smooth, Plot.LINE);
plot.setColor(Color.red);
plot.addPoints(qScaled, line, Plot.LINE);
plot.setColor(Color.black);
plot.addLabel(0, 0, label);
Utils.display(title, plot, Utils.NO_TO_FRONT);
}
if (fitPrecision) {
// Q - Should this be optional?
if (sampleDecay) {
// If a sample of the precision was used to construct the data for the initial fit
// then update the estimate using the fit result since it will be a better start point.
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
// Do a multivariate fit ...
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
PointValuePair p = null;
MultiPlateauness f = new MultiPlateauness(frcnum, q, low, high);
double[] initial = new double[] { histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, qValue };
p = findMin(p, opt, f, scale(initial, 0.1));
p = findMin(p, opt, f, scale(initial, 0.5));
p = findMin(p, opt, f, initial);
p = findMin(p, opt, f, scale(initial, 2));
p = findMin(p, opt, f, scale(initial, 10));
if (p != null) {
double[] point = p.getPointRef();
histogram.mean = point[0] * images.nmPerPixel;
histogram.sigma = point[1] * images.nmPerPixel;
qValue = point[2];
}
} else {
// If so then this should be optional.
if (sampleDecay) {
if (precisionMethod != PrecisionMethod.FIXED) {
histogram.sigma = precision.getStandardDeviation();
// Normalise sum-of-squares to the SR pixel size
double meanSumOfSquares = (precision.getSumOfSquares() / (images.nmPerPixel * images.nmPerPixel)) / precision.getN();
histogram.mean = images.nmPerPixel * Math.sqrt(meanSumOfSquares - estimate[1] / (4 * Math.PI * Math.PI));
}
exp_decay = computeExpDecay(histogram.mean / images.nmPerPixel, histogram.sigma / images.nmPerPixel, q);
}
// Estimate spurious component by promoting plateauness.
// The Matlab code used random initial points for a Simplex optimiser.
// A Brent line search should be pretty deterministic so do simple repeats.
// However it will proceed downhill so if the initial point is wrong then
// it will find a sub-optimal result.
UnivariateOptimizer o = new BrentOptimizer(1e-3, 1e-6);
Plateauness f = new Plateauness(frcnum, exp_decay, low, high);
UnivariatePointValuePair p = null;
p = findMin(p, o, f, qValue, 0.1);
p = findMin(p, o, f, qValue, 0.2);
p = findMin(p, o, f, qValue, 0.333);
p = findMin(p, o, f, qValue, 0.5);
// Do some Simplex repeats as well
SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10);
p = findMin(p, opt, f, qValue * 0.1);
p = findMin(p, opt, f, qValue * 0.5);
p = findMin(p, opt, f, qValue);
p = findMin(p, opt, f, qValue * 2);
p = findMin(p, opt, f, qValue * 10);
if (p != null)
qValue = p.getPoint();
}
QPlot qplot = new QPlot(frcCurve, qValue, low, high);
// Interactive dialog to estimate Q (blinking events per flourophore) using
// sliders for the mean and standard deviation of the localisation precision.
showQEstimationDialog(histogram, qplot, frcCurve, images.nmPerPixel);
IJ.showStatus(TITLE + " complete");
}
use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.
the class FreeFilterResults method run.
/*
* (non-)
*
* @see ij.plugin.PlugIn#run(java.lang.String)
*/
public void run(String arg) {
SMLMUsageTracker.recordPlugin(this.getClass(), arg);
if (MemoryPeakResults.isMemoryEmpty()) {
// Ask user if they want to show the demo filters
ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
gd.enableYesNoCancel();
gd.hideCancelButton();
gd.addMessage("No results in memory. Show the demo filters?");
gd.showDialog();
if (gd.wasOKed())
logDemoFilters(TITLE);
return;
}
if (!showDialog())
return;
results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
IJ.showStatus("");
return;
}
// Filter results
Filter filter = Filter.fromXML(filterSettings.freeFilter);
if (filter != null) {
MemoryPeakResults newResults = filter.filter(results);
if (newResults.size() > 0) {
newResults.setName(results.getName() + " Free Filtered");
MemoryPeakResults.addResults(newResults);
}
IJ.showStatus(String.format("Filtered %d results to %d", results.size(), newResults.size()));
} else {
IJ.showStatus("ERROR: Unable to create filter");
}
}
use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.
the class FIRE method run.
/*
* (non-Javadoc)
*
* @see ij.plugin.PlugIn#run(java.lang.String)
*/
public void run(String arg) {
extraOptions = Utils.isExtraOptions();
SMLMUsageTracker.recordPlugin(this.getClass(), arg);
// Require some fit results and selected regions
int size = MemoryPeakResults.countMemorySize();
if (size == 0) {
IJ.error(TITLE, "There are no fitting results in memory");
return;
}
if ("q".equals(arg)) {
TITLE += " Q estimation";
runQEstimation();
return;
}
IJ.showStatus(TITLE + " ...");
if (!showInputDialog())
return;
MemoryPeakResults results = ResultsManager.loadInputResults(inputOption, false);
if (results == null || results.size() == 0) {
IJ.error(TITLE, "No results could be loaded");
return;
}
MemoryPeakResults results2 = ResultsManager.loadInputResults(inputOption2, false);
results = cropToRoi(results);
if (results.size() < 2) {
IJ.error(TITLE, "No results within the crop region");
return;
}
if (results2 != null) {
results2 = cropToRoi(results2);
if (results2.size() < 2) {
IJ.error(TITLE, "No results2 within the crop region");
return;
}
}
initialise(results, results2);
if (!showDialog())
return;
long start = System.currentTimeMillis();
// Compute FIRE
String name = results.getName();
double fourierImageScale = SCALE_VALUES[imageScaleIndex];
int imageSize = IMAGE_SIZE_VALUES[imageSizeIndex];
if (this.results2 != null) {
name += " vs " + results2.getName();
FireResult result = calculateFireNumber(fourierMethod, samplingMethod, thresholdMethod, fourierImageScale, imageSize);
if (result != null) {
logResult(name, result);
if (showFRCCurve)
showFrcCurve(name, result, thresholdMethod);
}
} else {
FireResult result = null;
int repeats = (randomSplit) ? Math.max(1, FIRE.repeats) : 1;
if (repeats == 1) {
result = calculateFireNumber(fourierMethod, samplingMethod, thresholdMethod, fourierImageScale, imageSize);
if (result != null) {
logResult(name, result);
if (showFRCCurve)
showFrcCurve(name, result, thresholdMethod);
}
} else {
// Multi-thread this ...
int nThreads = Maths.min(repeats, getThreads());
ExecutorService executor = Executors.newFixedThreadPool(nThreads);
TurboList<Future<?>> futures = new TurboList<Future<?>>(repeats);
TurboList<FIREWorker> workers = new TurboList<FIREWorker>(repeats);
setProgress(repeats);
IJ.showProgress(0);
IJ.showStatus(TITLE + " computing ...");
for (int i = 1; i <= repeats; i++) {
FIREWorker w = new FIREWorker(i, fourierImageScale, imageSize);
workers.add(w);
futures.add(executor.submit(w));
}
// Wait for all to finish
for (int t = futures.size(); t-- > 0; ) {
try {
// The future .get() method will block until completed
futures.get(t).get();
} catch (Exception e) {
// This should not happen.
// Ignore it and allow processing to continue (the number of neighbour samples will just be smaller).
e.printStackTrace();
}
}
IJ.showProgress(1);
executor.shutdown();
// Show a combined FRC curve plot of all the smoothed curves if we have multiples.
LUT valuesLUT = LUTHelper.createLUT(LutColour.FIRE_GLOW);
@SuppressWarnings("unused") LUT // Black at max value
noSmoothLUT = LUTHelper.createLUT(LutColour.GRAYS).createInvertedLut();
LUTHelper.DefaultLUTMapper mapper = new LUTHelper.DefaultLUTMapper(0, repeats);
FrcCurve curve = new FrcCurve();
Statistics stats = new Statistics();
WindowOrganiser wo = new WindowOrganiser();
boolean oom = false;
for (int i = 0; i < repeats; i++) {
FIREWorker w = workers.get(i);
if (w.oom)
oom = true;
if (w.result == null)
continue;
result = w.result;
if (!Double.isNaN(result.fireNumber))
stats.add(result.fireNumber);
if (showFRCCurveRepeats) {
// Output each FRC curve using a suffix.
logResult(w.name, result);
wo.add(Utils.display(w.plot.getTitle(), w.plot));
}
if (showFRCCurve) {
int index = mapper.map(i + 1);
//@formatter:off
curve.add(name, result, thresholdMethod, LUTHelper.getColour(valuesLUT, index), Color.blue, //LUTHelper.getColour(noSmoothLUT, index)
null);
//@formatter:on
}
}
if (result != null) {
wo.cascade();
double mean = stats.getMean();
logResult(name, result, mean, stats);
if (showFRCCurve) {
curve.addResolution(mean);
Plot2 plot = curve.getPlot();
Utils.display(plot.getTitle(), plot);
}
}
if (oom) {
//@formatter:off
IJ.error(TITLE, "ERROR - Parallel computation out-of-memory.\n \n" + TextUtils.wrap("The number of results will be reduced. " + "Please reduce the size of the Fourier image " + "or change the number of threads " + "using the extra options (hold down the 'Shift' " + "key when running the plugin).", 80));
//@formatter:on
}
}
// Only do this once
if (showFRCTimeEvolution && result != null && !Double.isNaN(result.fireNumber))
showFrcTimeEvolution(name, result.fireNumber, thresholdMethod, nmPerPixel / result.getNmPerPixel(), imageSize);
}
IJ.showStatus(TITLE + " complete : " + Utils.timeToString(System.currentTimeMillis() - start));
}
use of gdsc.smlm.results.MemoryPeakResults in project GDSC-SMLM by aherbert.
the class FilterAnalysis method showDialog.
private boolean showDialog(List<MemoryPeakResults> resultsList, boolean fileInput) {
GenericDialog gd = new GenericDialog(TITLE);
gd.addHelp(About.HELP_URL);
int total = 0;
int tp = 0;
for (MemoryPeakResults r : resultsList) {
total += r.size();
for (PeakResult p : r.getResults()) if (p.origValue != 0)
tp++;
}
gd.addMessage(String.format("%d files, %d results, %d True-Positives", resultsList.size(), total, tp));
if (!fileInput) {
gd.addCheckbox("SNR_filter", snrFilter);
gd.addNumericField("Min_SNR", minSnr, 0);
gd.addNumericField("Max_SNR", maxSnr, 0);
gd.addNumericField("Min_Width", minWidth, 2);
gd.addNumericField("Max_Width", maxWidth, 2);
gd.addNumericField("Increment_Width", incWidth, 2);
gd.addCheckbox("Precision_filter", precisionFilter);
gd.addNumericField("Min_Precision", minPrecision, 0);
gd.addNumericField("Max_Precision", maxPrecision, 0);
gd.addCheckbox("Trace_filter", traceFilter);
gd.addNumericField("Min_distance", minDistance, 2);
gd.addNumericField("Max_distance", maxDistance, 2);
gd.addNumericField("Increment_distance", incDistance, 2);
gd.addNumericField("Min_time", minTime, 0);
gd.addNumericField("Max_time", maxTime, 0);
gd.addNumericField("Increment_time", incTime, 0);
gd.addCheckbox("Hysteresis_SNR_filter", hysteresisSnrFilter);
gd.addNumericField("Min_SNR_gap", minSnrGap, 0);
gd.addNumericField("Max_SNR_gap", maxSnrGap, 0);
gd.addNumericField("Increment_SNR_gap", incSnrGap, 0);
gd.addCheckbox("Hysteresis_Precision_filter", hysteresisPrecisionFilter);
gd.addNumericField("Min_Precision_gap", minPrecisionGap, 0);
gd.addNumericField("Max_Precision_gap", maxPrecisionGap, 0);
gd.addNumericField("Increment_Precision_gap", incPrecisionGap, 0);
gd.addCheckbox("Save_filters", saveFilterSets);
}
gd.addCheckbox("Show_table", showResultsTable);
gd.addSlider("Plot_top_n", 0, 20, plotTopN);
gd.addCheckbox("Calculate_sensitivity", calculateSensitivity);
gd.addSlider("Delta", 0.01, 1, delta);
if (!fileInput) {
// Re-arrange to 2 columns
if (gd.getLayout() != null) {
GridBagLayout grid = (GridBagLayout) gd.getLayout();
Component splitLabel = (Component) gd.getCheckboxes().get(3);
int xOffset = 0, yOffset = 0;
int lastY = -1, rowCount = 0;
for (Component comp : gd.getComponents()) {
// Check if this should be the second major column
if (comp == splitLabel) {
xOffset += 2;
yOffset -= rowCount;
rowCount = 0;
}
// Reposition the field
GridBagConstraints c = grid.getConstraints(comp);
if (lastY != c.gridy)
rowCount++;
lastY = c.gridy;
c.gridx = c.gridx + xOffset;
c.gridy = c.gridy + yOffset;
c.insets.left = c.insets.left + 10 * xOffset;
c.insets.top = 0;
c.insets.bottom = 0;
grid.setConstraints(comp, c);
}
if (IJ.isLinux())
gd.setBackground(new Color(238, 238, 238));
}
}
gd.showDialog();
if (gd.wasCanceled() || !readDialog(gd, fileInput))
return false;
return true;
}
Aggregations