Java tutorial
package gdsc.smlm.ij.plugins.pcpalm; /*----------------------------------------------------------------------------- * GDSC SMLM Software * * Copyright (C) 2013 Alex Herbert * Genome Damage and Stability Centre * University of Sussex, UK * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 3 of the License, or * (at your option) any later version. *---------------------------------------------------------------------------*/ import gdsc.smlm.function.SkewNormalFunction; import gdsc.smlm.function.gaussian.Gaussian2DFunction; import gdsc.smlm.ij.IJTrackProgress; import gdsc.smlm.ij.plugins.About; import gdsc.smlm.ij.plugins.Parameters; import gdsc.smlm.ij.plugins.ResultsManager; import gdsc.smlm.ij.plugins.ResultsManager.InputSource; import gdsc.smlm.ij.utils.Utils; import gdsc.smlm.model.MaskDistribution; import gdsc.smlm.model.StandardFluorophoreSequenceModel; import gdsc.smlm.model.UniformDistribution; import gdsc.smlm.results.MemoryPeakResults; import gdsc.smlm.results.NullSource; import gdsc.smlm.results.PeakResult; import gdsc.smlm.results.Trace; import gdsc.smlm.results.TraceManager; import gdsc.smlm.results.clustering.Cluster; import gdsc.smlm.results.clustering.ClusterPoint; import gdsc.smlm.results.clustering.ClusteringAlgorithm; import gdsc.smlm.results.clustering.ClusteringEngine; import gdsc.smlm.utils.Maths; import gdsc.smlm.utils.Statistics; import gdsc.smlm.utils.StoredDataStatistics; import ij.IJ; import ij.ImagePlus; import ij.Prefs; import ij.WindowManager; import ij.gui.GenericDialog; import ij.gui.Plot2; import ij.measure.Calibration; import ij.plugin.PlugIn; import ij.plugin.frame.Recorder; import ij.process.ByteProcessor; import ij.process.ImageProcessor; import ij.process.ShortProcessor; import java.awt.Color; import java.awt.Rectangle; import java.util.ArrayList; import java.util.Arrays; import java.util.LinkedList; import java.util.List; import org.apache.commons.math3.analysis.MultivariateFunction; import org.apache.commons.math3.analysis.MultivariateMatrixFunction; import org.apache.commons.math3.analysis.MultivariateVectorFunction; import org.apache.commons.math3.analysis.function.Gaussian; import org.apache.commons.math3.exception.NotStrictlyPositiveException; import org.apache.commons.math3.exception.TooManyEvaluationsException; import org.apache.commons.math3.fitting.GaussianFitter; import org.apache.commons.math3.optim.InitialGuess; import org.apache.commons.math3.optim.MaxEval; import org.apache.commons.math3.optim.MaxIter; import org.apache.commons.math3.optim.PointValuePair; import org.apache.commons.math3.optim.PointVectorValuePair; import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; import org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction; import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.NelderMeadSimplex; import org.apache.commons.math3.optim.nonlinear.scalar.noderiv.SimplexOptimizer; import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; import org.apache.commons.math3.optim.nonlinear.vector.Target; import org.apache.commons.math3.optim.nonlinear.vector.Weight; import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer; import org.apache.commons.math3.random.RandomDataGenerator; import org.apache.commons.math3.random.RandomGenerator; import org.apache.commons.math3.random.Well19937c; import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics; import org.apache.commons.math3.util.FastMath; /** * Use the PC-PALM protocol to prepare a set of localisations into molecules. This can be used for for clustering * analysis. * <p> * See Sengupta, et al (2013). Quantifying spatial resolution in point-localisation superresolution images using pair * correlation analysis. Nature Protocols 8, pp345-354. * <p> * See also Veatch, et al (2012). Correlation Functions Quantify Super-Resolution Images and Estimate Apparent * Clustering Due to Over-Counting. PLoS One 7, Issue 2, e31457 */ public class PCPALMMolecules implements PlugIn { static String TITLE = "PC-PALM Molecules"; private static String inputOption = ""; private static boolean chooseRoi = false; private static double nmPerPixelLimit = 0; private static String[] RUN_MODE = { "PC-PALM", "Manual Tracing", "In-memory results", "Simulation" }; private static int runMode = 0; // Mode 0: PC-PALM protocol for estimating localisation precision and then tracing molecules private static String roiImage = ""; private static int histogramBins = 50; private static String[] singlesMode = new String[] { "Ignore", "Include in molecules histogram", "Include in final filtering" }; private static int singlesModeIndex = 1; private static boolean simplexFitting = false; private static boolean showHistograms = true; private static boolean binaryImage = true; static double blinkingRate = 2; private static double p = 0.6; private static String[] BLINKING_DISTRIBUTION = new String[] { "Poisson", "Geometric", "None", "Binomial" }; private static int blinkingDistribution = 0; private static boolean clearResults = false; // Mode 1. Manual tracing of molecules private static double dThreshold = 150; private static double tThreshold = 1; // Mode 2. Direct use of in-memory results // - No parameters needed // Mode 3. Random simulation of molecules private static int nMolecules = 2000; private static double simulationSize = 16; private static boolean distanceAnalysis = false; private static String[] CLUSTER_SIMULATION = new String[] { "None", "Circles", "Non-overlapping circles", "Circles Mask" }; private static int clusterSimulation = 0; private static double clusterNumber = 3; private static double clusterNumberSD = 0; private static double clusterRadius = 50; private static boolean showClusterMask = false; // Low resolution image construction private static int lowResolutionImageSize = 1024; private static double roiSizeInUm = 4; private static boolean showHighResolutionImage = false; private Rectangle roiBounds; private int roiImageWidth, roiImageHeight; private long start; // These package level variables are used by the PCPALMAnalysis plugin. static MemoryPeakResults results; static double minx, miny, maxx, maxy; static double nmPerPixel; static ArrayList<Molecule> molecules = null; static double sigmaS = 20; static double densityPeaks; static double densityProtein; static double seconds; static double area; /* * (non-Javadoc) * * @see ij.plugin.PlugIn#run(java.lang.String) */ public void run(String arg) { // Require some fit results and selected regions boolean resultsAvailable = MemoryPeakResults.countMemorySize() > 0; if (!getRunMode(resultsAvailable)) return; if (runMode != 3) { results = ResultsManager.loadInputResults(inputOption, true); if (results == null || results.size() == 0) { IJ.error(TITLE, "No results could be loaded"); return; } if (results.getCalibration() == null) { IJ.error(TITLE, "Results are not calibrated"); return; } // Get the lifetime before cropping as this is the true representation of the number of frames. // This may truncate the lifetime if the first/last localisation are not near the end of the // acquisition lifetime getLifetime(); results = cropToRoi(results); if (results.size() == 0) { IJ.error(TITLE, "No results within the crop region"); return; } } // Clear cached results molecules = null; // Different run-modes for generating the set of molecules for analysis switch (runMode) { case 0: runPCPALM(); break; case 1: runManualTracing(); break; case 2: runInMemoryResults(); break; case 3: runSimulation(resultsAvailable); area = simulationSize * simulationSize; seconds = 100; // Use an arbitrary lifetime break; } if (molecules == null) return; if (molecules.size() < 2) { IJ.error(TITLE, "Not enough molecules to construct a binary image"); return; } // Generate binary PALM image if (!createImage(molecules)) return; // Density is required for the PC analysis densityPeaks = calculatePeakDensity(); if (runMode == 0 || runMode == 3) { // Blinking rate is mentioned in the PC-PALM protocol and so we include it here. // TODO - Add automated estimation of the blinking rate from the data using the method of // Annibale, et al (2011), Quantitative photo activated localization microscopy: unraveling the // effects of photoblinking. PLoS One, 6(7): e22678 (http://dx.doi.org/10.1371%2Fjournal.pone.0022678) densityProtein = densityPeaks / blinkingRate; log("Peak Density = %s (um^-2). Protein Density = %s (um^-2)", Utils.rounded(densityPeaks * 1e6), Utils.rounded(densityProtein * 1e6)); } else { // No blinking rate for non PC-PALM methods. This can be configured in later plugins if required. blinkingRate = 1; densityProtein = densityPeaks; log("Molecule Density = %s (um^-2)", Utils.rounded(densityPeaks * 1e6)); } log("Results lifetime = %s s", Utils.rounded(seconds)); // Use a second plugin filter that will work on a region drawn on the binary image // and compute the PALM analysis double seconds = (System.currentTimeMillis() - start) / 1000.0; String msg = TITLE + " complete : " + seconds + "s"; IJ.showStatus(msg); log(msg); } private boolean getRunMode(boolean resultsAvailable) { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); // Build a list of all images with a region ROI List<String> titles = new LinkedList<String>(); if (WindowManager.getWindowCount() > 0) { for (int imageID : WindowManager.getIDList()) { ImagePlus imp = WindowManager.getImage(imageID); if (imp != null && imp.getRoi() != null && imp.getRoi().isArea()) titles.add(imp.getTitle()); } } if (!resultsAvailable) { runMode = 3; gd.addMessage( "Simulate molecules for cluster analysis.\nComputes a binary image from localisation data"); gd.addNumericField("Molecules", nMolecules, 0); gd.addNumericField("Simulation_size (um)", simulationSize, 2); gd.addNumericField("Blinking_rate", blinkingRate, 2); gd.addChoice("Blinking_distribution", BLINKING_DISTRIBUTION, BLINKING_DISTRIBUTION[blinkingDistribution]); gd.addNumericField("Average_precision (nm)", sigmaS, 2); gd.addCheckbox("Show_histograms", showHistograms); gd.addCheckbox("Distance_analysis", distanceAnalysis); gd.addChoice("Cluster_simulation", CLUSTER_SIMULATION, CLUSTER_SIMULATION[clusterSimulation]); gd.addNumericField("Cluster_number", clusterNumber, 2); gd.addNumericField("Cluster_variation (SD)", clusterNumberSD, 2); gd.addNumericField("Cluster_radius", clusterRadius, 2); gd.addCheckbox("Show_cluster_mask", showClusterMask); Recorder.recordOption("Run_mode", RUN_MODE[runMode]); } else { gd.addMessage( "Prepare molecules for cluster analysis.\nComputes a binary image from raw localisation data"); ResultsManager.addInput(gd, inputOption, InputSource.MEMORY); if (!titles.isEmpty()) gd.addCheckbox((titles.size() == 1) ? "Use_ROI" : "Choose_ROI", chooseRoi); gd.addChoice("Run_mode", RUN_MODE, RUN_MODE[runMode]); } gd.addMessage("Select options for low resolution image:"); gd.addSlider("Image_size (px)", 512, 2048, lowResolutionImageSize); gd.addSlider("ROI_size (um)", 1.5, 4, roiSizeInUm); gd.addMessage("Select options for high resolution image:"); gd.addCheckbox("Show_high_res_image", showHighResolutionImage); gd.addSlider("nm_per_pixel_limit", 0, 20, nmPerPixelLimit); gd.addMessage("Optionally remove all analysis results from memory"); gd.addCheckbox("Clear_results", clearResults); gd.showDialog(); if (gd.wasCanceled()) return false; if (!resultsAvailable) { nMolecules = (int) Math.abs(gd.getNextNumber()); simulationSize = Math.abs(gd.getNextNumber()); blinkingRate = Math.abs(gd.getNextNumber()); blinkingDistribution = gd.getNextChoiceIndex(); sigmaS = Math.abs(gd.getNextNumber()); showHistograms = gd.getNextBoolean(); distanceAnalysis = gd.getNextBoolean(); clusterSimulation = gd.getNextChoiceIndex(); clusterNumber = Math.abs(gd.getNextNumber()); clusterNumberSD = Math.abs(gd.getNextNumber()); clusterRadius = Math.abs(gd.getNextNumber()); showClusterMask = gd.getNextBoolean(); } else { inputOption = ResultsManager.getInputSource(gd); if (!titles.isEmpty()) chooseRoi = gd.getNextBoolean(); runMode = gd.getNextChoiceIndex(); } lowResolutionImageSize = (int) gd.getNextNumber(); roiSizeInUm = gd.getNextNumber(); showHighResolutionImage = gd.getNextBoolean(); nmPerPixelLimit = Math.abs(gd.getNextNumber()); clearResults = gd.getNextBoolean(); // Check arguments try { if (!resultsAvailable) { Parameters.isAboveZero("Molecules", nMolecules); Parameters.isAboveZero("Simulation size", simulationSize); Parameters.isEqualOrAbove("Blinking rate", blinkingRate, 1); Parameters.isEqualOrAbove("Cluster number", clusterNumber, 1); } Parameters.isAbove("Image scale", lowResolutionImageSize, 1); Parameters.isAboveZero("ROI size", roiSizeInUm); } catch (IllegalArgumentException ex) { IJ.error(TITLE, ex.getMessage()); return false; } if (!titles.isEmpty() && chooseRoi && resultsAvailable) { if (titles.size() == 1) { roiImage = titles.get(0); Recorder.recordOption("Image", roiImage); } else { String[] items = titles.toArray(new String[titles.size()]); gd = new GenericDialog(TITLE); gd.addMessage("Select the source image for the ROI"); gd.addChoice("Image", items, roiImage); gd.showDialog(); if (gd.wasCanceled()) return false; roiImage = gd.getNextChoice(); } ImagePlus imp = WindowManager.getImage(roiImage); roiBounds = imp.getRoi().getBounds(); roiImageWidth = imp.getWidth(); roiImageHeight = imp.getHeight(); } else { roiBounds = null; } if (!resultsAvailable) { if (!getPValue()) return false; } if (clearResults) { PCPALMAnalysis.results.clear(); PCPALMFitting.previous_gr = null; } return true; } private MemoryPeakResults cropToRoi(MemoryPeakResults results) { Rectangle bounds = results.getBounds(true); area = (bounds.width * bounds.height * results.getNmPerPixel() * results.getNmPerPixel()) / 1e6; if (roiBounds == null) { return results; } // Adjust bounds relative to input results image double xscale = (double) roiImageWidth / bounds.width; double yscale = (double) roiImageHeight / bounds.height; roiBounds.x /= xscale; roiBounds.width /= xscale; roiBounds.y /= yscale; roiBounds.height /= yscale; float minX = (int) (roiBounds.x); float maxX = (int) Math.ceil(roiBounds.x + roiBounds.width); float minY = (int) (roiBounds.y); float maxY = (int) Math.ceil(roiBounds.y + roiBounds.height); // Update the area with the cropped region area *= (maxX - minX) / bounds.width; area *= (maxY - minY) / bounds.height; // 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; } private void runPCPALM() { if (!showPCPALMDialog()) return; startLog(); // Follow the PC-PALM protocol log("Fitting localisation precision..."); ArrayList<Molecule> localisations = extractLocalisations(results); double sigmaRaw = calculateAveragePrecision(localisations, "Localisations"); log("%d localisations with an average precision of %.2f", results.size(), sigmaRaw); log("Fitting molecule precision..."); ArrayList<Molecule> singles = new ArrayList<Molecule>(); molecules = extractMolecules(results, sigmaRaw, singles); if (singlesModeIndex == 1) molecules.addAll(singles); sigmaS = calculateAveragePrecision(molecules, "Molecules"); log("%d molecules with an average precision of %.2f", molecules.size(), sigmaS); // Q. Should this filter the original localisations or just the grouped peaks? if (singlesModeIndex == 2) molecules.addAll(singles); molecules = filterMolecules(molecules, sigmaS); log("%d molecules within precision %.2f", molecules.size(), 3 * sigmaS); } private void startLog() { logSpacer(); log(TITLE); logSpacer(); start = System.currentTimeMillis(); } private boolean showPCPALMDialog() { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); gd.addMessage( "Estimate the average localisation precision by fitting histograms.\nUse the precision to trace localisations into molecule pulses."); gd.addNumericField("Histogram_bins", histogramBins, 0); gd.addChoice("Singles_mode", singlesMode, singlesMode[singlesModeIndex]); gd.addCheckbox("Simplex_fit", simplexFitting); gd.addCheckbox("Show_histograms", showHistograms); gd.addCheckbox("Binary_image", binaryImage); gd.addNumericField("Blinking_rate", blinkingRate, 2); gd.showDialog(); if (gd.wasCanceled()) return false; histogramBins = (int) gd.getNextNumber(); singlesModeIndex = gd.getNextChoiceIndex(); simplexFitting = gd.getNextBoolean(); showHistograms = gd.getNextBoolean(); binaryImage = gd.getNextBoolean(); blinkingRate = gd.getNextNumber(); // Check arguments try { Parameters.isAbove("Histogram bins", histogramBins, 1); Parameters.isEqualOrAbove("Blinking rate", blinkingRate, 1); } catch (IllegalArgumentException ex) { IJ.error(TITLE, ex.getMessage()); return false; } return true; } /** * Extract molecules for the PC-PALM analysis. * <p> * Estimate the localisation uncertainty (precision) of each molecule using the formula of Mortensen, et al (2010), * Nature Methods 7, 377-381. Store distance in nm and signal in photons using the calibration * * @param results * @return */ public ArrayList<Molecule> extractLocalisations(MemoryPeakResults results) { ArrayList<Molecule> molecules = new ArrayList<Molecule>(results.size()); final double nmPerPixel = results.getNmPerPixel(); final double gain = results.getGain(); final boolean emCCD = results.isEMCCD(); for (PeakResult r : results.getResults()) { double p = r.getPrecision(nmPerPixel, gain, emCCD); // Remove EMCCD adjustment //p /= Math.sqrt(PeakResult.F); molecules.add(new Molecule(r.getXPosition() * nmPerPixel, r.getYPosition() * nmPerPixel, p, r.getSignal() / gain)); } return molecules; } /** * Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision distribution. * * @param molecules * @param subTitle * @return The average precision */ private double calculateAveragePrecision(ArrayList<Molecule> molecules, String subTitle) { String title = (showHistograms) ? TITLE + " Histogram " + subTitle : null; return calculateAveragePrecision(molecules, title, histogramBins, true, true); } /** * Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision distribution. * <p> * A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does not fit within 3 SDs * of the simple mean then the simple mean is returned. * * @param molecules * @param title * the plot title (null if no plot should be displayed) * @param histogramBins * @param logFitParameters * Record the fit parameters to the ImageJ log * @param removeOutliers * The distribution is created using all values within 1.5x the inter-quartile range (IQR) of the data * @return The average precision */ public double calculateAveragePrecision(ArrayList<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) { // Plot histogram of the precision float[] data = new float[molecules.size()]; DescriptiveStatistics stats = new DescriptiveStatistics(); double yMin = Double.NEGATIVE_INFINITY, yMax = 0; for (int i = 0; i < data.length; i++) { data[i] = (float) molecules.get(i).precision; stats.addValue(data[i]); } // Set the min and max y-values using 1.5 x IQR if (removeOutliers) { double lower = stats.getPercentile(25); double upper = stats.getPercentile(75); if (Double.isNaN(lower) || Double.isNaN(upper)) { if (logFitParameters) Utils.log("Error computing IQR: %f - %f", lower, upper); } else { double iqr = upper - lower; yMin = FastMath.max(lower - iqr, stats.getMin()); yMax = FastMath.min(upper + iqr, stats.getMax()); if (logFitParameters) Utils.log(" Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax); } } if (yMin == Double.NEGATIVE_INFINITY) { yMin = stats.getMin(); yMax = stats.getMax(); if (logFitParameters) Utils.log(" Data range: %f - %f", yMin, yMax); } float[][] hist = Utils.calcHistogram(data, yMin, yMax, histogramBins); Plot2 plot = null; if (title != null) { plot = new Plot2(title, "Precision", "Frequency"); float[] xValues = hist[0]; float[] yValues = hist[1]; if (xValues.length > 0) { double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]); plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.max(yValues) * 1.05); } plot.addPoints(xValues, yValues, Plot2.BAR); Utils.display(title, plot); } // Extract non-zero data float[] x = Arrays.copyOf(hist[0], hist[0].length); float[] y = hist[1]; int count = 0; float dx = (x[1] - x[0]) * 0.5f; for (int i = 0; i < y.length; i++) if (y[i] > 0) { x[count] = x[i] + dx; y[count] = y[i]; count++; } x = Arrays.copyOf(x, count); y = Arrays.copyOf(y, count); // Sense check to fitted data. Get mean and SD of histogram double[] stats2 = Utils.getHistogramStatistics(x, y); double mean = stats2[0]; if (logFitParameters) log(" Initial Statistics: %f +/- %f", stats2[0], stats2[1]); // Standard Gaussian fit double[] parameters = fitGaussian(x, y); if (parameters == null) { log(" Failed to fit initial Gaussian"); return mean; } double newMean = parameters[1]; double error = Math.abs(stats2[0] - newMean) / stats2[1]; if (error > 3) { log(" Failed to fit Gaussian: %f standard deviations from histogram mean", error); return mean; } if (newMean < yMin || newMean > yMax) { log(" Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax); return mean; } mean = newMean; if (logFitParameters) log(" Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]); double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 }; // Fit to a skewed Gaussian (or appropriate function) double[] skewParameters = fitSkewGaussian(x, y, initialSolution); if (skewParameters == null) { log(" Failed to fit Skewed Gaussian"); return mean; } SkewNormalFunction sn = new SkewNormalFunction(skewParameters); if (logFitParameters) log(" Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance())); newMean = sn.getMean(); error = Math.abs(stats2[0] - newMean) / stats2[1]; if (error > 3) { log(" Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error); return mean; } if (newMean < yMin || newMean > yMax) { log(" Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax); return mean; } // Use original histogram x-axis to maintain all the bins if (plot != null) { x = hist[0]; for (int i = 0; i < y.length; i++) x[i] += dx; plot.setColor(Color.red); addToPlot(plot, x, skewParameters, Plot2.LINE); plot.setColor(Color.black); Utils.display(title, plot); } // Return the average precision from the fitted curve return newMean; } private double[] fitGaussian(float[] x, float[] y) { MyGaussianFitter fitter = new MyGaussianFitter( new org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer(), 2000); for (int i = 0; i < x.length; i++) fitter.addObservedPoint(x[i], y[i]); double[] parameters; try { parameters = fitter.fit(); } catch (TooManyEvaluationsException e) { // Use the initial estimate parameters = fitter.guess; } catch (Exception e) { // Just in case there is another exception type, or the initial estimate failed return null; } return parameters; } private double[] fitSkewGaussian(float[] x, float[] y, double[] initialSolution) { try { double[] skewParameters = (simplexFitting) ? optimiseSimplex(x, y, initialSolution) : optimiseLeastSquares(x, y, initialSolution); return skewParameters; } catch (TooManyEvaluationsException e) { return null; } } private double[] optimiseLeastSquares(float[] x, float[] y, double[] initialSolution) { // Least-squares optimisation using numerical gradients final SkewNormalDifferentiableFunction myFunction = new SkewNormalDifferentiableFunction(initialSolution); myFunction.addData(x, y); LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(); PointVectorValuePair optimum = optimizer.optimize(new MaxIter(3000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return myFunction.jacobian(point); } }), new ModelFunction(myFunction), new Target(myFunction.calculateTarget()), new Weight(myFunction.calculateWeights()), new InitialGuess(initialSolution)); double[] skewParameters = optimum.getPoint(); return skewParameters; } private double[] optimiseSimplex(float[] x, float[] y, double[] initialSolution) { // Simplex optimisation SkewNormalMultivariateFunction sn2 = new SkewNormalMultivariateFunction(initialSolution); sn2.addData(x, y); NelderMeadSimplex simplex = new NelderMeadSimplex(4); SimplexOptimizer opt = new SimplexOptimizer(1e-6, 1e-10); PointValuePair solution = opt.optimize(new MaxEval(1000), new InitialGuess(initialSolution), simplex, new ObjectiveFunction(sn2), GoalType.MINIMIZE); double[] skewParameters2 = solution.getPointRef(); return skewParameters2; } /** * Add the skewed gaussian to the histogram plot * * @param plot * @param x * @param parameters * Gaussian parameters * @param alpha * @param shape */ private void addToPlot(Plot2 plot, float[] x, double[] parameters, int shape) { SkewNormalFunction sn = new SkewNormalFunction(parameters); float[] y = new float[x.length]; for (int i = 0; i < x.length; i++) y[i] = (float) sn.evaluate(x[i]); plot.addPoints(x, y, shape); } /** * Group all localisations in successive frames within 2.5x of the initial precision estimate into a single molecule * * @param results * The results * @param sigmaRaw * The initial precision estimate * @param singles * a list of the singles (not grouped into molecules) * @return a list of molecules */ private ArrayList<Molecule> extractMolecules(MemoryPeakResults results, double sigmaRaw, ArrayList<Molecule> singles) { return traceMolecules(results, sigmaRaw * 2.5, 1, singles); } /** * Trace localisations * * @param results * The results * @param distance * The distance threshold (nm) * @param time * The time threshold (frames) * @param singles * a list of the singles (not grouped into molecules) * @return a list of molecules */ private ArrayList<Molecule> traceMolecules(MemoryPeakResults results, double distance, int time, ArrayList<Molecule> singles) { TraceManager tm = new TraceManager(results); double distanceThreshold = distance / results.getNmPerPixel(); tm.traceMolecules(distanceThreshold, time); Trace[] traces = tm.getTraces(); ArrayList<Molecule> molecules = new ArrayList<Molecule>(traces.length); final double nmPerPixel = results.getNmPerPixel(); final double gain = results.getGain(); final boolean emCCD = results.isEMCCD(); for (Trace t : traces) { double p = t.getLocalisationPrecision(nmPerPixel, gain, emCCD); if (t.size() == 1) { float[] centroid = t.getCentroid(); singles.add( new Molecule(centroid[0] * nmPerPixel, centroid[1] * nmPerPixel, p, t.getSignal() / gain)); } else { float[] centroid = t.getCentroid(); molecules.add( new Molecule(centroid[0] * nmPerPixel, centroid[1] * nmPerPixel, p, t.getSignal() / gain)); } } log(" %d localisations traced to %d molecules (%d singles, %d traces) using d=%.2f nm, t=%d frames (%s s)", results.size(), molecules.size() + singles.size(), singles.size(), molecules.size(), distance, time, Utils.rounded(time * results.getCalibration().exposureTime / 1000.0)); return molecules; } /** * Calculate the density of peaks in the original data * * @return The peak density */ private double calculatePeakDensity() { //double pcw = PCPALMMolecules.maxx - PCPALMMolecules.minx; //double pch = PCPALMMolecules.maxy - PCPALMMolecules.miny; //double area = pcw * pch; // Use the area from the source of the molecules return molecules.size() / (area * 1E6); } /** * Return a new list, removing all molecules with a precision over 3x of the precision estimate * * @param molecules * @param sigmaS * The precision estimate * @return */ private ArrayList<Molecule> filterMolecules(ArrayList<Molecule> molecules, double sigmaS) { ArrayList<Molecule> newMolecules = new ArrayList<Molecule>(molecules.size()); final double limit = 3 * sigmaS; for (Molecule m : molecules) { if (m.precision <= limit) newMolecules.add(m); } return newMolecules; } private void runManualTracing() { if (!showManualTracingDialog()) return; startLog(); // Convert seconds to frames int timeInFrames = FastMath.max(1, (int) Math.round(tThreshold * 1000.0 / results.getCalibration().exposureTime)); ArrayList<Molecule> singles = new ArrayList<Molecule>(); molecules = traceMolecules(results, dThreshold, timeInFrames, singles); molecules.addAll(singles); } private boolean showManualTracingDialog() { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); gd.addMessage("Use distance and time thresholds to trace localisations into molecules."); gd.addNumericField("Distance (nm)", dThreshold, 0); gd.addNumericField("Time (seconds)", tThreshold, 2); gd.showDialog(); if (gd.wasCanceled()) return false; dThreshold = Math.abs(gd.getNextNumber()); tThreshold = Math.abs(gd.getNextNumber()); // Check arguments try { Parameters.isAboveZero("Distance threshold", dThreshold); Parameters.isAboveZero("Time threshold", tThreshold); } catch (IllegalArgumentException ex) { IJ.error(TITLE, ex.getMessage()); return false; } return true; } private void runInMemoryResults() { startLog(); molecules = extractLocalisations(results); } private void runSimulation(boolean resultsAvailable) { if (resultsAvailable && !showSimulationDialog()) return; startLog(); log("Simulation parameters"); if (blinkingDistribution == 3) { log(" - Clusters = %d", nMolecules); log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4)); log(" - Molecules/cluster = %s", Utils.rounded(blinkingRate, 4)); log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]); log(" - p-Value = %s", Utils.rounded(p, 4)); } else { log(" - Molecules = %d", nMolecules); log(" - Simulation size = %s um", Utils.rounded(simulationSize, 4)); log(" - Blinking rate = %s", Utils.rounded(blinkingRate, 4)); log(" - Blinking distribution = %s", BLINKING_DISTRIBUTION[blinkingDistribution]); } log(" - Average precision = %s nm", Utils.rounded(sigmaS, 4)); log(" - Clusters simulation = " + CLUSTER_SIMULATION[clusterSimulation]); if (clusterSimulation > 0) { log(" - Cluster number = %s +/- %s", Utils.rounded(clusterNumber, 4), Utils.rounded(clusterNumberSD, 4)); log(" - Cluster radius = %s nm", Utils.rounded(clusterRadius, 4)); } final double nmPerPixel = 100; double width = simulationSize * 1000.0; // Allow a border of 3 x sigma for +/- precision //if (blinkingRate > 1) width -= 3 * sigmaS; RandomGenerator randomGenerator = new Well19937c( System.currentTimeMillis() + System.identityHashCode(this)); RandomDataGenerator dataGenerator = new RandomDataGenerator(randomGenerator); UniformDistribution dist = new UniformDistribution(null, new double[] { width, width, 0 }, randomGenerator.nextInt()); molecules = new ArrayList<Molecule>(nMolecules); // Create some dummy results since the calibration is required for later analysis results = new MemoryPeakResults(); results.setCalibration(new gdsc.smlm.results.Calibration(nmPerPixel, 1, 100)); results.setSource(new NullSource("Molecule Simulation")); results.begin(); int count = 0; // Generate a sequence of coordinates ArrayList<double[]> xyz = new ArrayList<double[]>((int) (nMolecules * 1.1)); Statistics statsRadius = new Statistics(); Statistics statsSize = new Statistics(); String maskTitle = TITLE + " Cluster Mask"; ByteProcessor bp = null; double maskScale = 0; // TODO - Add a fluctuations model to this. if (clusterSimulation > 0) { // Simulate clusters. // Note: In the Veatch et al. paper (Plos 1, e31457) correlation functions are built using circles // with small radii of 4-8 Arbitrary Units (AU) or large radii of 10-30 AU. A fluctuations model is // created at T = 1.075 Tc. It is not clear exactly how the particles are distributed. // It may be that a mask is created first using the model. The particles are placed on the mask using // a specified density. This simulation produces a figure to show either a damped cosine function // (circles) or an exponential (fluctuations). The number of particles in each circle may be randomly // determined just by density. The figure does not discuss the derivation of the cluster size // statistic. // // If this plugin simulation is run with a uniform distribution and blinking rate of 1 then the damped // cosine function is reproduced. The curve crosses g(r)=1 at a value equivalent to the average // distance to the centre-of-mass of each drawn cluster, not the input cluster radius parameter (which // is a hard upper limit on the distance to centre). final int maskSize = lowResolutionImageSize; int[] mask = null; maskScale = width / maskSize; // scale is in nm/pixel ArrayList<double[]> clusterCentres = new ArrayList<double[]>(); int totalSteps = 1 + (int) Math.ceil(nMolecules / clusterNumber); if (clusterSimulation == 2 || clusterSimulation == 3) { // Clusters are non-overlapping circles // Ensure the circles do not overlap by using an exclusion mask that accumulates // out-of-bounds pixels by drawing the last cluster (plus some border) on an image. When no // more pixels are available then stop generating molecules. // This is done by cumulatively filling a mask and using the MaskDistribution to select // a new point. This may be slow but it works. // TODO - Allow clusters of different sizes... mask = new int[maskSize * maskSize]; Arrays.fill(mask, 255); MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); double[] centre; IJ.showStatus("Computing clusters mask"); int roiRadius = (int) Math.round((clusterRadius * 2) / maskScale); if (clusterSimulation == 3) { // Generate a mask of circles then sample from that. // If we want to fill the mask completely then adjust the total steps to be the number of // circles that can fit inside the mask. totalSteps = (int) (maskSize * maskSize / (Math.PI * Math.pow(clusterRadius / maskScale, 2))); } while ((centre = maskDistribution.next()) != null && clusterCentres.size() < totalSteps) { IJ.showProgress(clusterCentres.size(), totalSteps); // The mask returns the coordinates with the centre of the image at 0,0 centre[0] += width / 2; centre[1] += width / 2; clusterCentres.add(centre); // Fill in the mask around the centre to exclude any more circles that could overlap double cx = centre[0] / maskScale; double cy = centre[1] / maskScale; fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 0); //log("[%.1f,%.1f] @ [%.1f,%.1f]", centre[0], centre[1], cx, cy); //Utils.display("Mask", new ColorProcessor(maskSize, maskSize, mask)); try { maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); } catch (IllegalArgumentException e) { // This can happen when there are no more non-zero pixels log("WARNING: No more room for clusters on the mask area (created %d of estimated %d)", clusterCentres.size(), totalSteps); break; } } IJ.showProgress(1); IJ.showStatus(""); } else { // Clusters are overlapping circles // Pick centres randomly from the distribution while (clusterCentres.size() < totalSteps) clusterCentres.add(dist.next()); } if (showClusterMask || clusterSimulation == 3) { // Show the mask for the clusters if (mask == null) mask = new int[maskSize * maskSize]; else Arrays.fill(mask, 0); int roiRadius = (int) Math.round((clusterRadius) / maskScale); for (double[] c : clusterCentres) { double cx = c[0] / maskScale; double cy = c[1] / maskScale; fillMask(mask, maskSize, (int) cx, (int) cy, roiRadius, 1); } if (clusterSimulation == 3) { // We have the mask. Now pick points at random from the mask. MaskDistribution maskDistribution = new MaskDistribution(mask, maskSize, maskSize, 0, maskScale, maskScale, randomGenerator); // Allocate each molecule position to a parent circle so defining clusters. int[][] clusters = new int[clusterCentres.size()][]; int[] clusterSize = new int[clusters.length]; for (int i = 0; i < nMolecules; i++) { double[] centre = maskDistribution.next(); // The mask returns the coordinates with the centre of the image at 0,0 centre[0] += width / 2; centre[1] += width / 2; xyz.add(centre); // Output statistics on cluster size and number. // TODO - Finding the closest cluster could be done better than an all-vs-all comparison double max = distance2(centre, clusterCentres.get(0)); int cluster = 0; for (int j = 1; j < clusterCentres.size(); j++) { double d2 = distance2(centre, clusterCentres.get(j)); if (d2 < max) { max = d2; cluster = j; } } // Assign point i to cluster centre[2] = cluster; if (clusterSize[cluster] == 0) { clusters[cluster] = new int[10]; } if (clusters[cluster].length <= clusterSize[cluster]) { clusters[cluster] = Arrays.copyOf(clusters[cluster], (int) (clusters[cluster].length * 1.5)); } clusters[cluster][clusterSize[cluster]++] = i; } // Generate real cluster size statistics for (int j = 0; j < clusterSize.length; j++) { final int size = clusterSize[j]; if (size == 0) continue; statsSize.add(size); if (size == 1) { statsRadius.add(0); continue; } // Find centre of cluster and add the distance to each point double[] com = new double[2]; for (int n = 0; n < size; n++) { double[] xy = xyz.get(clusters[j][n]); for (int k = 0; k < 2; k++) com[k] += xy[k]; } for (int k = 0; k < 2; k++) com[k] /= size; for (int n = 0; n < size; n++) { double dx = xyz.get(clusters[j][n])[0] - com[0]; double dy = xyz.get(clusters[j][n])[1] - com[1]; statsRadius.add(Math.sqrt(dx * dx + dy * dy)); } } } if (showClusterMask) { bp = new ByteProcessor(maskSize, maskSize); for (int i = 0; i < mask.length; i++) if (mask[i] != 0) bp.set(i, 128); Utils.display(maskTitle, bp); } } // Use the simulated cluster centres to create clusters of the desired size if (clusterSimulation == 1 || clusterSimulation == 2) { for (double[] clusterCentre : clusterCentres) { int clusterN = (int) Math.round( (clusterNumberSD > 0) ? dataGenerator.nextGaussian(clusterNumber, clusterNumberSD) : clusterNumber); if (clusterN < 1) continue; //double[] clusterCentre = dist.next(); if (clusterN == 1) { // No need for a cluster around a point xyz.add(clusterCentre); statsRadius.add(0); statsSize.add(1); } else { // Generate N random points within a circle of the chosen cluster radius. // Locate the centre-of-mass and the average distance to the centre. double[] com = new double[3]; int j = 0; while (j < clusterN) { // Generate a random point within a circle uniformly // http://stackoverflow.com/questions/5837572/generate-a-random-point-within-a-circle-uniformly double t = 2.0 * Math.PI * randomGenerator.nextDouble(); double u = randomGenerator.nextDouble() + randomGenerator.nextDouble(); double r = clusterRadius * ((u > 1) ? 2 - u : u); double x = r * Math.cos(t); double y = r * Math.sin(t); double[] xy = new double[] { clusterCentre[0] + x, clusterCentre[1] + y }; xyz.add(xy); for (int k = 0; k < 2; k++) com[k] += xy[k]; j++; } // Add the distance of the points from the centre of the cluster. // Note this does not account for the movement due to precision. statsSize.add(j); if (j == 1) { statsRadius.add(0); } else { for (int k = 0; k < 2; k++) com[k] /= j; while (j > 0) { double dx = xyz.get(xyz.size() - j)[0] - com[0]; double dy = xyz.get(xyz.size() - j)[1] - com[1]; statsRadius.add(Math.sqrt(dx * dx + dy * dy)); j--; } } } } } } else { // Random distribution for (int i = 0; i < nMolecules; i++) xyz.add(dist.next()); } // The Gaussian sigma should be applied so the overall distance from the centre // ( sqrt(x^2+y^2) ) has a standard deviation of sigmaS? final double sigma1D = sigmaS / Math.sqrt(2); // Show optional histograms StoredDataStatistics intraDistances = null; StoredDataStatistics blinks = null; if (showHistograms) { int capacity = (int) (xyz.size() * blinkingRate); intraDistances = new StoredDataStatistics(capacity); blinks = new StoredDataStatistics(capacity); } Statistics statsSigma = new Statistics(); for (int i = 0; i < xyz.size(); i++) { int nOccurrences = getBlinks(dataGenerator, blinkingRate); if (showHistograms) blinks.add(nOccurrences); final int size = molecules.size(); // Get coordinates in nm final double[] moleculeXyz = xyz.get(i); if (bp != null && nOccurrences > 0) { bp.putPixel((int) Math.round(moleculeXyz[0] / maskScale), (int) Math.round(moleculeXyz[1] / maskScale), 255); } while (nOccurrences-- > 0) { final double[] localisationXy = Arrays.copyOf(moleculeXyz, 2); // Add random precision if (sigma1D > 0) { final double dx = dataGenerator.nextGaussian(0, sigma1D); final double dy = dataGenerator.nextGaussian(0, sigma1D); localisationXy[0] += dx; localisationXy[1] += dy; if (!dist.isWithinXY(localisationXy)) continue; // Calculate mean-squared displacement statsSigma.add(dx * dx + dy * dy); } final double x = localisationXy[0]; final double y = localisationXy[1]; molecules.add(new Molecule(x, y, i, 1)); // Store in pixels float[] params = new float[7]; params[Gaussian2DFunction.X_POSITION] = (float) (x / nmPerPixel); params[Gaussian2DFunction.Y_POSITION] = (float) (y / nmPerPixel); results.add(i + 1, (int) x, (int) y, 0, 0, 0, params, null); } if (molecules.size() > size) { count++; if (showHistograms) { int newCount = molecules.size() - size; if (newCount == 1) { // No intra-molecule distances //intraDistances.add(0); continue; } // Get the distance matrix between these molecules double[][] matrix = new double[newCount][newCount]; for (int ii = size, x = 0; ii < molecules.size(); ii++, x++) { for (int jj = size + 1, y = 1; jj < molecules.size(); jj++, y++) { final double d2 = molecules.get(ii).distance2(molecules.get(jj)); matrix[x][y] = matrix[y][x] = d2; } } // Get the maximum distance for particle linkage clustering of this molecule double max = 0; for (int x = 0; x < newCount; x++) { // Compare to all-other molecules and get the minimum distance // needed to join at least one double linkDistance = Double.POSITIVE_INFINITY; for (int y = 0; y < newCount; y++) { if (x == y) continue; if (matrix[x][y] < linkDistance) linkDistance = matrix[x][y]; } // Check if this is larger if (max < linkDistance) max = linkDistance; } intraDistances.add(Math.sqrt(max)); } } } results.end(); if (bp != null) Utils.display(maskTitle, bp); // Used for debugging //System.out.printf(" * Molecules = %d (%d activated)\n", xyz.size(), count); //if (clusterSimulation > 0) // System.out.printf(" * Cluster number = %s +/- %s. Radius = %s +/- %s\n", // Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4), // Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4)); log("Simulation results"); log(" * Molecules = %d (%d activated)", xyz.size(), count); log(" * Blinking rate = %s", Utils.rounded((double) molecules.size() / xyz.size(), 4)); log(" * Precision (Mean-displacement) = %s nm", (statsSigma.getN() > 0) ? Utils.rounded(Math.sqrt(statsSigma.getMean()), 4) : "0"); if (showHistograms) { if (intraDistances.getN() == 0) { log(" * Mean Intra-Molecule particle linkage distance = 0 nm"); log(" * Fraction of inter-molecule particle linkage @ 0 nm = 0 %%"); } else { plot(blinks, "Blinks/Molecule", true); double[][] intraHist = plot(intraDistances, "Intra-molecule particle linkage distance", false); // Determine 95th and 99th percentile int p99 = intraHist[0].length - 1; double limit1 = 0.99 * intraHist[1][p99]; double limit2 = 0.95 * intraHist[1][p99]; while (intraHist[1][p99] > limit1 && p99 > 0) p99--; int p95 = p99; while (intraHist[1][p95] > limit2 && p95 > 0) p95--; log(" * Mean Intra-Molecule particle linkage distance = %s nm (95%% = %s, 99%% = %s, 100%% = %s)", Utils.rounded(intraDistances.getMean(), 4), Utils.rounded(intraHist[0][p95], 4), Utils.rounded(intraHist[0][p99], 4), Utils.rounded(intraHist[0][intraHist[0].length - 1], 4)); if (distanceAnalysis) { performDistanceAnalysis(intraHist, p99); } } } if (clusterSimulation > 0) { log(" * Cluster number = %s +/- %s", Utils.rounded(statsSize.getMean(), 4), Utils.rounded(statsSize.getStandardDeviation(), 4)); log(" * Cluster radius = %s +/- %s nm (mean distance to centre-of-mass)", Utils.rounded(statsRadius.getMean(), 4), Utils.rounded(statsRadius.getStandardDeviation(), 4)); } } private double[][] plot(StoredDataStatistics stats, String label, boolean integerBins) { String title = TITLE + " " + label; Plot2 plot; double[][] hist = null; if (integerBins) { // The histogram is not need for the return statement Utils.showHistogram(title, stats, label, 1, 0, 0); } else { // Show a cumulative histogram so that the bin size is not relevant hist = Maths.cumulativeHistogram(stats.getValues(), false); // Create the axes double[] xValues = hist[0]; double[] yValues = hist[1]; // Plot plot = new Plot2(title, label, "Frequency", xValues, yValues); Utils.display(title, plot); } return hist; } private void performDistanceAnalysis(double[][] intraHist, int p99) { // We want to know the fraction of distances between molecules at the 99th percentile // that are intra- rather than inter-molecule. // Do single linkage clustering of closest pair at this distance and count the number of // links that are inter and intra. // Convert molecules for clustering ArrayList<ClusterPoint> points = new ArrayList<ClusterPoint>(molecules.size()); for (Molecule m : molecules) // Precision was used to store the molecule ID points.add(ClusterPoint.newClusterPoint((int) m.precision, m.x, m.y, m.photons)); ClusteringEngine engine = new ClusteringEngine(Prefs.getThreads(), ClusteringAlgorithm.PARTICLE_SINGLE_LINKAGE, new IJTrackProgress()); IJ.showStatus("Clustering to check inter-molecule distances"); engine.setTrackJoins(true); ArrayList<Cluster> clusters = engine.findClusters(points, intraHist[0][p99]); IJ.showStatus(""); if (clusters != null) { double[] intraIdDistances = engine.getIntraIdDistances(); double[] interIdDistances = engine.getInterIdDistances(); int all = interIdDistances.length + intraIdDistances.length; log(" * Fraction of inter-molecule particle linkage @ %s nm = %s %%", Utils.rounded(intraHist[0][p99], 4), (all > 0) ? Utils.rounded(100.0 * interIdDistances.length / all, 4) : "0"); // Show a double cumulative histogram plot double[][] intraIdHist = Maths.cumulativeHistogram(intraIdDistances, false); double[][] interIdHist = Maths.cumulativeHistogram(interIdDistances, false); // Plot String title = TITLE + " molecule linkage distance"; Plot2 plot = new Plot2(title, "Distance", "Frequency", intraIdHist[0], intraIdHist[1]); double max = (intraIdHist[1].length > 0) ? intraIdHist[1][intraIdHist[1].length - 1] : 0; if (interIdHist[1].length > 0) max = FastMath.max(max, interIdHist[1][interIdHist[1].length - 1]); plot.setLimits(0, intraIdHist[0][intraIdHist[0].length - 1], 0, max); plot.setColor(Color.blue); plot.addPoints(interIdHist[0], interIdHist[1], Plot2.LINE); plot.setColor(Color.black); Utils.display(title, plot); } else { log("Aborted clustering to check inter-molecule distances"); } } private double distance2(double[] centre1, double[] centre2) { final double dx = centre1[0] - centre2[0]; final double dy = centre1[1] - centre2[1]; return dx * dx + dy * dy; } /** * Fill the given mask with the fill value using a circle with the specified centre and radius * * @param mask * @param maskSize * @param cx * @param cy * @param roiRadius * @param fill */ private void fillMask(int[] mask, final int maskSize, final int cx, final int cy, final int roiRadius, final int fill) { int minx = cx - roiRadius; int maxx = cx + roiRadius; int miny = cy - roiRadius; int maxy = cy + roiRadius; final int r2 = roiRadius * roiRadius; // Pre-calculate x range if (minx < 0) minx = 0; if (maxx >= maskSize) maxx = maskSize - 1; if (minx > maxx) return; int n = 0; int[] dx2 = new int[roiRadius * 2 + 1]; for (int x = minx; x <= maxx; x++) { dx2[n++] = (cx - x) * (cx - x); } if (miny < 0 || maxy >= maskSize) { for (int y = miny, dy = -roiRadius; y <= maxy; y++, dy++) { if (y < 0) continue; if (y >= maskSize) break; final int limit = r2 - (dy * dy); for (int i = y * maskSize + minx, nn = 0; nn < n; i++, nn++) { if (dx2[nn] <= limit) mask[i] = fill; } } } else { for (int y = miny, dy = -roiRadius; y <= maxy; y++, dy++) { final int limit = r2 - (dy * dy); for (int i = y * maskSize + minx, nn = 0; nn < n; i++, nn++) { if (dx2[nn] <= limit) mask[i] = fill; } } } } private int getBlinks(RandomDataGenerator dataGenerator, double averageBlinks) { switch (blinkingDistribution) { case 3: // Binomial distribution return dataGenerator.nextBinomial((int) Math.round(averageBlinks), p); case 2: return (int) Math.round(averageBlinks); case 1: return StandardFluorophoreSequenceModel.getBlinks(true, dataGenerator, averageBlinks); default: return StandardFluorophoreSequenceModel.getBlinks(false, dataGenerator, averageBlinks); } } private boolean showSimulationDialog() { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); gd.addMessage("Simulate a random distribution of molecules."); gd.addNumericField("Molecules", nMolecules, 0); gd.addNumericField("Simulation_size (um)", simulationSize, 2); gd.addNumericField("Blinking_rate", blinkingRate, 2); gd.addChoice("Blinking_distribution", BLINKING_DISTRIBUTION, BLINKING_DISTRIBUTION[blinkingDistribution]); gd.addNumericField("Average_precision (nm)", sigmaS, 2); gd.addCheckbox("Show_histograms", showHistograms); gd.addCheckbox("Distance_analysis", distanceAnalysis); gd.addChoice("Cluster_simulation", CLUSTER_SIMULATION, CLUSTER_SIMULATION[clusterSimulation]); gd.addNumericField("Cluster_number", clusterNumber, 2); gd.addNumericField("Cluster_variation (SD)", clusterNumberSD, 2); gd.addNumericField("Cluster_radius", clusterRadius, 2); gd.addCheckbox("Show_cluster_mask", showClusterMask); gd.showDialog(); if (gd.wasCanceled()) return false; nMolecules = (int) Math.abs(gd.getNextNumber()); simulationSize = Math.abs(gd.getNextNumber()); blinkingRate = Math.abs(gd.getNextNumber()); blinkingDistribution = gd.getNextChoiceIndex(); sigmaS = Math.abs(gd.getNextNumber()); showHistograms = gd.getNextBoolean(); distanceAnalysis = gd.getNextBoolean(); clusterSimulation = gd.getNextChoiceIndex(); clusterNumber = Math.abs(gd.getNextNumber()); clusterNumberSD = Math.abs(gd.getNextNumber()); clusterRadius = Math.abs(gd.getNextNumber()); showClusterMask = gd.getNextBoolean(); // Check arguments try { Parameters.isAboveZero("Molecules", nMolecules); Parameters.isAboveZero("Simulation size", simulationSize); Parameters.isEqualOrAbove("Blinking rate", blinkingRate, 1); Parameters.isEqualOrAbove("Cluster number", clusterNumber, 1); } catch (IllegalArgumentException ex) { IJ.error(TITLE, ex.getMessage()); return false; } return getPValue(); } private boolean getPValue() { if (blinkingDistribution == 3) { GenericDialog gd = new GenericDialog(TITLE); gd.addMessage("Binomial distribution requires a p-value"); gd.addSlider("p-Value (%)", 0, 100, 100 * p); gd.showDialog(); if (gd.wasCanceled()) return false; p = FastMath.max(FastMath.min(gd.getNextNumber(), 100), 0) / 100; } return true; } /** * Get the lifetime of the results using the earliest and latest frames and the calibrated exposure time. */ private void getLifetime() { int start; int end; List<PeakResult> peakResults = results.getResults(); if (peakResults.isEmpty()) { seconds = 0; return; } start = end = peakResults.get(0).peak; for (PeakResult r : peakResults) { if (start > r.peak) start = r.peak; if (end < r.getEndFrame()) end = r.getEndFrame(); } seconds = (end - start + 1) * results.getCalibration().exposureTime / 1000; } private boolean createImage(ArrayList<Molecule> molecules) { if (molecules.isEmpty()) return false; // Find the limits of the image minx = maxx = molecules.get(0).x; miny = maxy = molecules.get(0).y; // Compute limits for (int i = molecules.size(); i-- > 0;) { final Molecule m1 = molecules.get(i); if (minx > m1.x) minx = m1.x; else if (maxx < m1.x) maxx = m1.x; if (miny > m1.y) miny = m1.y; else if (maxy < m1.y) maxy = m1.y; } // // Debug all-vs-all comparison // long start = System.nanoTime(); // double dMin2 = Double.POSITIVE_INFINITY; // for (int i = molecules.size(); i-- > 0;) // { // final Molecule m1 = molecules.get(i); // for (int j = i; j-- > 0;) // { // final Molecule m2 = molecules.get(j); // if (dMin2 > m1.distance2(m2)) // dMin2 = m1.distance2(m2); // } // } // long t1 = System.nanoTime() - start; // Assign to a grid final int gridSize = 500; final double xBinSize = (maxx - minx) / gridSize; final double yBinSize = (maxy - miny) / gridSize; final int nXBins = 1 + (int) ((maxx - minx) / xBinSize); final int nYBins = 1 + (int) ((maxy - miny) / yBinSize); Molecule[][] grid = new Molecule[nXBins][nYBins]; for (Molecule m : molecules) { final int xBin = (int) ((m.x - minx) / xBinSize); final int yBin = (int) ((m.y - miny) / yBinSize); // Build a single linked list m.next = grid[xBin][yBin]; grid[xBin][yBin] = m; } // Find the minimum distance between molecules. double dMin = Double.POSITIVE_INFINITY; IJ.showStatus("Computing minimum distance ..."); IJ.showProgress(0); Molecule[] neighbours = new Molecule[5]; for (int yBin = 0, currentIndex = 0, finalIndex = nXBins * nYBins; yBin < nYBins; yBin++) { for (int xBin = 0; xBin < nXBins; xBin++) { IJ.showProgress(currentIndex, finalIndex); for (Molecule m1 = grid[xBin][yBin]; m1 != null; m1 = m1.next) { // Build a list of which cells to compare up to a maximum of 4 // | 0,0 | 1,0 // ------------+----- // -1,1 | 0,1 | 1,1 int count = 0; neighbours[count++] = m1.next; if (yBin < nYBins - 1) { neighbours[count++] = grid[xBin][yBin + 1]; if (xBin > 0) neighbours[count++] = grid[xBin - 1][yBin + 1]; } if (xBin < nXBins - 1) { neighbours[count++] = grid[xBin + 1][yBin]; if (yBin < nYBins - 1) neighbours[count++] = grid[xBin + 1][yBin + 1]; } // Compare to neighbours while (count-- > 0) { for (Molecule m2 = neighbours[count]; m2 != null; m2 = m2.next) { if (dMin > m1.distance2(m2)) dMin = m1.distance2(m2); } } } } } IJ.showStatus(""); IJ.showProgress(1); nmPerPixel = Math.sqrt(dMin); log("Minimum distance between molecules = %g nm", nmPerPixel); if (nmPerPixel == 0 && nmPerPixelLimit == 0) { IJ.error(TITLE, "Zero minimum distance between molecules - please enter a nm/pixel limit for image reconstruction"); return false; } if (nmPerPixel < nmPerPixelLimit) { log("Minimum distance adjusted to user defined limit %.2f nm", nmPerPixelLimit); nmPerPixel = nmPerPixelLimit; } // Compute the minimum size we can use and stay within memory. // Assume a 4um x 4um analysis section with 800nm feature radius. double limit = (4000.0 + 800.0 + 1) / 4096; if (nmPerPixel < limit) { log("Minimum distance adjusted to %.2f nm to fit in memory", limit); nmPerPixel = limit; } log("X-range %.2f - %.2f : Y-range %.2f - %.2f (nm)", minx, maxx, miny, maxy); // Construct a binary representation String namePrefix = results.getName() + " " + ((binaryImage) ? "Binary" : "Count") + " Image"; double lowResNmPerPixel; double xrange = maxx - minx; double yrange = maxy - miny; if (xrange > 0 || yrange > 0) { lowResNmPerPixel = FastMath.max(xrange, yrange) / lowResolutionImageSize; } else { // The resolution does not matter lowResNmPerPixel = 100; } ImagePlus imp = displayImage(namePrefix + " (low res)", molecules, minx, miny, maxx, maxy, lowResNmPerPixel, false, binaryImage); // Add an ROI to allow the user to select regions. PC-PALM recommends 2x2 to 4x4 um^2 int size = (int) (roiSizeInUm * 1000.0 / lowResNmPerPixel); imp.setRoi(new Rectangle(0, 0, size, size)); if (showHighResolutionImage) displayImage(namePrefix + " (high res)", molecules, minx, miny, maxx, maxy, nmPerPixel, false, binaryImage); // Store the molecules, the data range and the dMin. // This will be used by a filter plugin that crops sections from the image for PC analysis PCPALMMolecules.molecules = molecules; return true; } static ImageProcessor drawImage(ArrayList<Molecule> molecules, double minx, double miny, double maxx, double maxy, double nmPerPixel, boolean checkBounds, boolean binary) { double scalex = maxx - minx; double scaley = maxy - miny; int width = (int) Math.round(scalex / nmPerPixel) + 1; int height = (int) Math.round(scaley / nmPerPixel) + 1; // *** // The PC-PALM + PLoS One papers describe using a binary image. // However both papers provide MatLab code where the number of particles is // calculated using sum(sum(I)). This indicates a non binary image could be input // to the routine to calculate the correlation function g(r). // *** if (binary) { byte[] data = new byte[width * height]; for (Molecule m : molecules) { if (checkBounds) { if (m.x < minx || m.x >= maxx || m.y < miny || m.y >= maxy) continue; } // Shift to the origin. This makes the image more memory efficient. int x = (int) Math.round((m.x - minx) / nmPerPixel); int y = (int) Math.round((m.y - miny) / nmPerPixel); int index = y * width + x; // Construct a binary image data[index] = (byte) 1; } ByteProcessor ip = new ByteProcessor(width, height, data, null); ip.setMinAndMax(0, 1); return ip; } else { short[] data = new short[width * height]; for (Molecule m : molecules) { if (checkBounds) { if (m.x < minx || m.x >= maxx || m.y < miny || m.y >= maxy) continue; } // Shift to the origin. This makes the image more memory efficient. int x = (int) Math.round((m.x - minx) / nmPerPixel); int y = (int) Math.round((m.y - miny) / nmPerPixel); int index = y * width + x; // Construct a count image data[index]++; } ShortProcessor ip = new ShortProcessor(width, height, data, null); ip.setMinAndMax(0, Maths.max(data)); return ip; } } static ImagePlus displayImage(String title, ImageProcessor ip, double nmPerPixel) { ImagePlus imp = Utils.display(title, ip); Calibration cal = new Calibration(); cal.setUnit("um"); cal.pixelWidth = cal.pixelHeight = nmPerPixel / 1000; imp.setCalibration(cal); return imp; } static ImagePlus displayImage(String title, ArrayList<Molecule> molecules, double minx, double miny, double maxx, double maxy, double nmPerPixel, boolean checkBounds, boolean binary) { ImageProcessor ip = drawImage(molecules, minx, miny, maxx, maxy, nmPerPixel, checkBounds, binary); return displayImage(title, ip, nmPerPixel); } /** * Allow optimisation using Apache Commons Math 3 Optimiser */ public abstract class SkewNormalOptimiserFunction extends SkewNormalFunction { public SkewNormalOptimiserFunction(double[] parameters) { super(parameters); } protected List<Double> x = null; protected List<Double> y = null; public void addPoint(double x, double y) { if (this.x == null) { this.x = new ArrayList<Double>(); this.y = new ArrayList<Double>(); } this.x.add(x); this.y.add(y); } public void addData(float[] x, float[] y) { this.x = new ArrayList<Double>(); this.y = new ArrayList<Double>(); for (int i = 0; i < x.length; i++) { this.x.add((double) x[i]); this.y.add((double) y[i]); } } public double[] calculateTarget() { double[] target = new double[y.size()]; for (int i = 0; i < y.size(); i++) { target[i] = y.get(i).doubleValue(); } return target; } public double[] calculateWeights() { double[] w = new double[y.size()]; for (int i = 0; i < y.size(); i++) { w[i] = 1; } return w; } } /** * Allow optimisation using Apache Commons Math 3 Gradient Optimiser */ public class SkewNormalDifferentiableFunction extends SkewNormalOptimiserFunction implements MultivariateVectorFunction { // Adapted from http://commons.apache.org/proper/commons-math/userguide/optimization.html // Use the deprecated API since the new one is not yet documented. public SkewNormalDifferentiableFunction(double[] parameters) { super(parameters); } private double[][] jacobian(double[] variables) { // Compute the gradients using numerical differentiation double[][] jacobian = new double[x.size()][4]; final double delta = 0.001; double[][] d = new double[variables.length][variables.length]; for (int i = 0; i < variables.length; i++) d[i][i] = delta * Math.abs(variables[i]); // Should the delta be changed for each parameter ? for (int i = 0; i < jacobian.length; ++i) { double x = this.x.get(i); double value = evaluate(x, variables); for (int j = 0; j < variables.length; j++) { double value2 = evaluate(x, variables[0] + d[0][j], variables[1] + d[1][j], variables[2] + d[2][j], variables[3] + d[3][j]); jacobian[i][j] = (value2 - value) / d[j][j]; } } return jacobian; } /* * (non-Javadoc) * * @see org.apache.commons.math3.analysis.MultivariateVectorFunction#value(double[]) */ public double[] value(double[] variables) { double[] values = new double[x.size()]; for (int i = 0; i < values.length; i++) values[i] = evaluate(x.get(i), variables); return values; } } /** * Allow optimisation using Apache Commons Math 3 Simplex */ public class SkewNormalMultivariateFunction extends SkewNormalOptimiserFunction implements MultivariateFunction { public SkewNormalMultivariateFunction(double[] parameters) { super(parameters); } /* * (non-Javadoc) * * @see org.apache.commons.math3.analysis.MultivariateFunction#value(double[]) */ public double value(double[] point) { // Objective function is to minimise sum-of-squares double ss = 0; for (int i = x.size(); i-- > 0;) { double dx = y.get(i) - evaluate(x.get(i), point); ss += dx * dx; } return ss; } } /** * Apache GaussianFitter extended to allow maxEvaluations to be set * */ public class MyGaussianFitter extends GaussianFitter { public int maxEvaluations = 1000; public double[] guess = null; public MyGaussianFitter(MultivariateVectorOptimizer optimizer, int maxEvaluations) { super(optimizer); this.maxEvaluations = maxEvaluations; } /* * (non-Javadoc) * * @see org.apache.commons.math3.fitting.GaussianFitter#fit(double[]) */ public double[] fit(double[] initialGuess) { final Gaussian.Parametric f = new Gaussian.Parametric() { @Override public double value(double x, double... p) { double v = Double.POSITIVE_INFINITY; try { v = super.value(x, p); } catch (NotStrictlyPositiveException e) { // NOPMD // Do nothing. } return v; } @Override public double[] gradient(double x, double... p) { double[] v = { Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY }; try { v = super.gradient(x, p); } catch (NotStrictlyPositiveException e) { // NOPMD // Do nothing. } return v; } }; return fit(maxEvaluations, f, initialGuess); } /** * Fits a Gaussian function to the observed points. * * @return the parameters of the Gaussian function that best fits the * observed points (in the same order as above). */ public double[] fit() { guess = (new ParameterGuesser(getObservations())).guess(); return fit(guess); } } /** * Log a message to the IJ log window * * @param format * @param args */ private static void log(String format, Object... args) { Utils.log(format, args); } /** * Output a spacer to the ImageJ log */ static void logSpacer() { log("-=-=-=-=-=-=-"); } }