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 edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D; import edu.emory.mathcs.jtransforms.fft.FloatFFT_2D; import gdsc.smlm.ij.plugins.About; import gdsc.smlm.ij.plugins.Parameters; import gdsc.smlm.ij.utils.Utils; import gdsc.smlm.model.MaskDistribution; import gdsc.smlm.utils.ImageWindow; import gdsc.smlm.utils.Maths; import gdsc.smlm.utils.Statistics; import ij.IJ; import ij.ImagePlus; import ij.WindowManager; import ij.gui.GenericDialog; import ij.gui.Plot; import ij.gui.Roi; import ij.gui.Plot2; import ij.plugin.filter.PlugInFilter; import ij.process.FHT2; import ij.process.FloatProcessor; import ij.process.ImageProcessor; import ij.text.TextWindow; import java.awt.Color; import java.awt.Frame; import java.awt.Rectangle; import java.io.File; import java.io.FileInputStream; import java.io.FileOutputStream; import java.io.FilenameFilter; import java.io.IOException; import java.util.ArrayList; import java.util.Arrays; import java.util.Collections; import org.apache.commons.math3.util.FastMath; import com.thoughtworks.xstream.XStream; import com.thoughtworks.xstream.XStreamException; import com.thoughtworks.xstream.io.xml.DomDriver; /** * Use the PC-PALM protocol to analyse a set of molecules to produce a correlation curve. * <p> * Frequency domain analysis see Sengupta, et al (2013). Quantifying spatial resolution in point-localisation * superresolution images using pair correlation analysis. Nature Protocols 8, pp345-354. * <p> * Spatial domain analysis see Puchnar, et al (2013). Counting molecules in single organelles with superresolution * microscopy allows tracking of the endosome maturation trajectory. PNAS. doi:10.1073/pnas.1309676110 */ public class PCPALMAnalysis implements PlugInFilter { static String TITLE = "PC-PALM Analysis"; private static String resultsDirectory = ""; private static double correlationDistance = 800; // nm private static double correlationInterval = 20; // nm private static boolean binaryImage = false; static double blinkingRate = -1; private static double copiedBlinkingRate = -1; private static double nmPerPixel = -1; private static double copiedNmPerPixel = -1; private static boolean showErrorBars = false; private static boolean applyWindow = false; private static boolean showHighResolutionImage = false; private static boolean showCorrelationImages = false; private static boolean useBorder = true; // Limits for the molecules constructed from the input ROI private double minx, miny, maxx, maxy; private double area, weightedArea, weightedAreaInPx; private int areaInPx, noOfMolecules; private double uniquePoints; // Cache the ImageWindow for faster repeat processing private static ImageWindow imageWindow = new ImageWindow(); // Used for the results table private static TextWindow resultsTable = null; static ArrayList<CorrelationResult> results = new ArrayList<CorrelationResult>(); private boolean spatialDomain; // Area of the region cropped from the PCPALM Molecules list double croppedArea = 0; /* * (non-Javadoc) * * @see ij.plugin.filter.PlugInFilter#setup(java.lang.String, ij.ImagePlus) */ public int setup(String arg, ImagePlus imp) { if ("save".equalsIgnoreCase(arg)) return saveResults(); if ("load".equalsIgnoreCase(arg)) return loadResults(); spatialDomain = "spatial".equalsIgnoreCase(arg); boolean noAreaRoi = (imp.getRoi() == null || !imp.getRoi().isArea()); if (imp == null || (!spatialDomain && noAreaRoi)) { error("Require an input image with an area ROI.\n" + "Please create a binary molecule image using " + PCPALMMolecules.TITLE); return DONE; } if (PCPALMMolecules.molecules == null || PCPALMMolecules.molecules.size() < 2) { error("Require a set of molecules for analysis.\n" + "Please create a binary molecule image using " + PCPALMMolecules.TITLE); return DONE; } if (!showDialog()) return DONE; PCPALMMolecules.logSpacer(); log(TITLE); PCPALMMolecules.logSpacer(); ArrayList<Molecule> molecules = cropToRoi(imp); if (molecules.size() < 2) { error("No results within the crop region"); return DONE; } log("Using %d molecules", molecules.size()); long start = System.currentTimeMillis(); analyse(molecules); double seconds = (System.currentTimeMillis() - start) / 1000.0; String msg = TITLE + " complete : " + seconds + "s"; IJ.showStatus(msg); log(msg); return DONE; } /** * Show a directory selection dialog for the results directory * * @return True if a directory was selected */ private boolean getDirectory() { resultsDirectory = Utils.getDirectory("Results_directory", resultsDirectory); return resultsDirectory != null; } /** * Save all the results to a directory * * @return DONE */ private int saveResults() { if (results.isEmpty()) { error("No results in memory"); } else if (getDirectory()) { XStream xs = new XStream(new DomDriver()); for (CorrelationResult result : results) saveResult(xs, result); } return DONE; } private void saveResult(XStream xs, CorrelationResult result) { String outputFilename = String.format("%s/%s.%d.xml", resultsDirectory, (result.spatialDomain) ? "Spatial" : "Frequency", result.id); FileOutputStream fs = null; try { fs = new FileOutputStream(outputFilename); xs.toXML(result, fs); } catch (XStreamException ex) { //ex.printStackTrace(); IJ.log("Failed to save correlation result to file: " + outputFilename); } catch (Exception e) { IJ.log("Failed to save correlation result to file: " + outputFilename); } finally { if (fs != null) { try { fs.close(); } catch (IOException e) { e.printStackTrace(); } } } } /** * Load all the results from a directory. File must have the XML suffix * * @return DONE */ private int loadResults() { if (getDirectory()) { File[] fileList = (new File(resultsDirectory)).listFiles(new FilenameFilter() { public boolean accept(File arg0, String arg1) { return arg1.endsWith("xml"); } }); if (fileList == null) return DONE; int count = 0; for (int i = 0; i < fileList.length; i++) { XStream xs = new XStream(new DomDriver()); if (fileList[i].isFile()) if (loadResult(xs, fileList[i].getPath())) count++; } if (count > 0) Collections.sort(results); log("Loaded %d results", count); } return DONE; } private boolean loadResult(XStream xs, String path) { FileInputStream fs = null; try { fs = new FileInputStream(path); CorrelationResult result = (CorrelationResult) xs.fromXML(fs); // Replace a result with the same id for (int i = 0; i < results.size(); i++) { if (results.get(i).id == result.id) { results.set(i, result); result = null; break; } } // Add to the results if we did not replace any if (result != null) results.add(result); return true; } catch (XStreamException ex) { //ex.printStackTrace(); IJ.log("Failed to load correlation result from file: " + path); } catch (Exception e) { IJ.log("Failed to load correlation result from file: " + path); } finally { if (fs != null) { try { fs.close(); } catch (IOException e) { e.printStackTrace(); } } } return false; } private void error(String message) { log("ERROR : " + message); IJ.error(TITLE, message); } /* * (non-Javadoc) * * @see ij.plugin.filter.PlugInFilter#run(ij.process.ImageProcessor) */ public void run(ImageProcessor ip) { // Nothing to do } private boolean showDialog() { GenericDialog gd = new GenericDialog(TITLE); gd.addHelp(About.HELP_URL); if (blinkingRate < 1 || copiedBlinkingRate != PCPALMMolecules.blinkingRate) { copiedBlinkingRate = blinkingRate = PCPALMMolecules.blinkingRate; } if (nmPerPixel < 1 || copiedNmPerPixel != PCPALMMolecules.nmPerPixel) { copiedNmPerPixel = nmPerPixel = PCPALMMolecules.nmPerPixel; } gd.addMessage("Analyse clusters using Pair Correlation."); gd.addNumericField("Correlation_distance (nm)", correlationDistance, 0); if (!spatialDomain) { gd.addMessage("-=- Frequency domain analysis -=-"); gd.addCheckbox("Binary_image", binaryImage); gd.addNumericField("Blinking_rate", blinkingRate, 2); gd.addNumericField("nm_per_pixel", nmPerPixel, 2); gd.addCheckbox("Show_error_bars", showErrorBars); gd.addCheckbox("Apply_window", applyWindow); gd.addCheckbox("Show_high_res_image", showHighResolutionImage); gd.addCheckbox("Show_correlation_images", showCorrelationImages); } else { gd.addMessage("-=- Spatial domain analysis -=-"); gd.addCheckbox("Use_border", useBorder); gd.addNumericField("Correlation_interval (nm)", correlationInterval, 0); } gd.showDialog(); if (gd.wasCanceled()) return false; correlationDistance = gd.getNextNumber(); if (!spatialDomain) { binaryImage = gd.getNextBoolean(); blinkingRate = gd.getNextNumber(); nmPerPixel = gd.getNextNumber(); showErrorBars = gd.getNextBoolean(); applyWindow = gd.getNextBoolean(); showHighResolutionImage = gd.getNextBoolean(); showCorrelationImages = gd.getNextBoolean(); } else { useBorder = gd.getNextBoolean(); correlationInterval = gd.getNextNumber(); } // Check arguments try { if (!spatialDomain) { Parameters.isAbove("Correlation distance", correlationDistance, 1); Parameters.isEqualOrAbove("Blinking_rate", blinkingRate, 1); Parameters.isAboveZero("nm per pixel", nmPerPixel); } else { Parameters.isAboveZero("Correlation interval", correlationInterval); } } catch (IllegalArgumentException ex) { error(ex.getMessage()); return false; } return true; } /** * Extract all the PC-PALM molecules that are within the image ROI region. The coordinates bounds are * converted using relative scaling to the limits of the PC-PALM molecules. If a non-rectangular ROI is used then a * mask is extracted and used for the crop. If no image is provided then the full set of molecules is returned. * <p> * Set the area property to the region covered by the molecules. * * @param imp * @return */ ArrayList<Molecule> cropToRoi(ImagePlus imp) { croppedArea = PCPALMMolecules.area; if (PCPALMMolecules.molecules == null || PCPALMMolecules.molecules.isEmpty()) { return PCPALMMolecules.molecules; } double pcw = PCPALMMolecules.maxx - PCPALMMolecules.minx; double pch = PCPALMMolecules.maxy - PCPALMMolecules.miny; Roi roi = (imp == null) ? null : imp.getRoi(); if (roi == null || !roi.isArea()) { log("Roi = %s nm x %s nm = %s um^2", Utils.rounded(pcw, 3), Utils.rounded(pch, 3), Utils.rounded(croppedArea, 3)); minx = PCPALMMolecules.minx; maxx = PCPALMMolecules.maxx; miny = PCPALMMolecules.miny; maxy = PCPALMMolecules.maxy; return PCPALMMolecules.molecules; } int w = imp.getWidth(); int h = imp.getHeight(); Rectangle bounds = roi.getBounds(); // Construct relative coordinates minx = PCPALMMolecules.minx + pcw * ((double) bounds.x / w); miny = PCPALMMolecules.miny + pch * ((double) bounds.y / h); maxx = PCPALMMolecules.minx + pcw * ((double) (bounds.x + bounds.width) / w); maxy = PCPALMMolecules.miny + pch * ((double) (bounds.y + bounds.height) / h); final double roix = maxx - minx; final double roiy = maxy - miny; ArrayList<Molecule> newMolecules = new ArrayList<Molecule>(PCPALMMolecules.molecules.size() / 2); // Support non-rectangular ROIs if (roi.getMask() != null) { MaskDistribution dist = createMaskDistribution(roi.getMask(), roix, roiy); double fraction = (double) dist.getSize() / (roi.getMask().getWidth() * roi.getMask().getHeight()); log("Roi is a mask of %d pixels", dist.getSize()); croppedArea = fraction * roix * roiy / 1e6; log("Roi area %s x %s nm x %s nm = %s um^2", Utils.rounded(fraction), Utils.rounded(roix, 3), Utils.rounded(roiy, 3), Utils.rounded(croppedArea, 3)); double[] xyz = new double[3]; // The mask is 0,0 in the centre so offset by this as well double xoffset = minx + roix / 2; double yoffset = miny + roiy / 2; for (Molecule m : PCPALMMolecules.molecules) { xyz[0] = m.x - xoffset; xyz[1] = m.y - yoffset; if (dist.isWithinXY(xyz)) newMolecules.add(m); } } else { croppedArea = roix * roiy / 1e6; log("Roi = %s nm x %s nm = %s um^2", Utils.rounded(roix, 3), Utils.rounded(roiy, 3), Utils.rounded(croppedArea, 3)); for (Molecule m : PCPALMMolecules.molecules) { if (m.x < minx || m.x >= maxx || m.y < miny || m.y >= maxy) continue; newMolecules.add(m); } } return newMolecules; } private MaskDistribution createMaskDistribution(ImageProcessor ip, double roix, double roiy) { // Calculate the scale of the mask final int w = ip.getWidth(); final int h = ip.getHeight(); final double scaleX = roix / w; final double scaleY = roiy / h; // Use an image for the distribution int[] mask = extractMask(ip); return new MaskDistribution(mask, w, h, 0, scaleX, scaleY); } private int[] extractMask(ImageProcessor ip) { int[] mask = new int[ip.getPixelCount()]; for (int i = 0; i < mask.length; i++) { mask[i] = ip.get(i); } return mask; } /** * Log a message to the IJ log window * * @param format * @param args */ private static void log(String format, Object... args) { IJ.log(String.format(format, args)); } /** * Perform the PC Analysis * * @param molecules */ private void analyse(ArrayList<Molecule> molecules) { // Check if the plots are currently shown String spatialPlotTitle = TITLE + " molecules/um^2"; String frequencyDomainTitle = TITLE + " g(r)"; boolean noPlots; String topPlotTitle; long start = System.nanoTime(); if (spatialDomain) { // ----------------- // Analysis in the spatial domain // ----------------- log("---"); log("Spatial domain analysis"); log("Computing density histogram"); // Compute all-vs-all distance matrix. // Create histogram of distances at different radii. final int nBins = (int) (correlationDistance / correlationInterval) + 1; final double maxDistance2 = correlationDistance * correlationDistance; int[] H = new int[nBins]; // TODO - Update this using a grid with a resolution of maxDistance to increase speed // by only comparing to neighbours within range. // An all-vs-all analysis does not account for a border. // A simple solution is to only process molecules within the border but compare them // to all molecules within the region. Thus every molecule has a complete circle of the max // radius around them to use: // ---------------------- // | | // | -------------- | // | | Within | | // | | Border | | // | | | | // | -------------- | // | Region | // ---------------------- // If the fraction of points within the correlation distance of the edge is low then this // will not make much difference. final double boundaryMinx = (useBorder) ? minx + correlationDistance : minx; final double boundaryMaxx = (useBorder) ? maxx - correlationDistance : maxx; final double boundaryMiny = (useBorder) ? miny + correlationDistance : miny; final double boundaryMaxy = (useBorder) ? maxy - correlationDistance : maxy; int N = 0; if (boundaryMaxx <= boundaryMinx || boundaryMaxy <= boundaryMiny) { log("ERROR: 'Use border' option of %s nm is not possible: Width = %s nm, Height = %s nm", Utils.rounded(correlationDistance, 4), Utils.rounded(maxx - minx, 3), Utils.rounded(maxy - miny, 3)); return; } else { for (int i = molecules.size(); i-- > 0;) { final Molecule m = molecules.get(i); // Optionally ignore molecules that are near the edge of the boundary if (useBorder && (m.x < boundaryMinx || m.x > boundaryMaxx || m.y < boundaryMiny || m.y > boundaryMaxy)) continue; N++; for (int j = molecules.size(); j-- > 0;) { if (i == j) continue; double d = m.distance2(molecules.get(j)); if (d < maxDistance2) { H[(int) (Math.sqrt(d) / correlationInterval)]++; } } } } double[] r = new double[nBins + 1]; for (int i = 0; i <= nBins; i++) r[i] = i * correlationInterval; double[] pcf = new double[nBins]; if (N > 0) { // Note: Convert nm^2 to um^2 final double N_pi = N * Math.PI / 1000000.0; for (int i = 0; i < nBins; i++) { // Pair-correlation is the count at the given distance divided by N and the area at distance ri: // H(r_i) / (N x (pi x (r_i+1)^2 - pi x r_i^2)) pcf[i] = H[i] / (N_pi * (r[i + 1] * r[i + 1] - r[i] * r[i])); } } // The final bin may be empty if the correlation interval was a factor of the correlation distance if (pcf[pcf.length - 1] == 0) { r = Arrays.copyOf(r, nBins - 1); pcf = Arrays.copyOf(pcf, nBins - 1); } else { r = Arrays.copyOf(r, nBins); } double[][] gr = new double[][] { r, pcf, null }; CorrelationResult result = new CorrelationResult(results.size() + 1, PCPALMMolecules.results.getSource(), boundaryMinx, boundaryMiny, boundaryMaxx, boundaryMaxy, N, correlationInterval, 0, false, gr, true); results.add(result); noPlots = WindowManager.getFrame(spatialPlotTitle) == null; topPlotTitle = frequencyDomainTitle; plotCorrelation(gr, 0, spatialPlotTitle, "molecules/um^2", true, false); } else { // ----------------- // Image correlation in the Frequency Domain // ----------------- log("Frequency domain analysis"); // Create a binary image for the molecules ImageProcessor im = PCPALMMolecules.drawImage(molecules, minx, miny, maxx, maxy, nmPerPixel, true, binaryImage); log("Image scale = %.2f nm/pixel : %d x %d pixels", nmPerPixel, im.getWidth(), im.getHeight()); // Apply a window function to the image to reduce FFT edge artifacts. if (applyWindow) { im = applyWindow(im, imageWindow); } if (showHighResolutionImage) { PCPALMMolecules.displayImage(PCPALMMolecules.results.getName() + " " + ((binaryImage) ? "Binary" : "Count") + " Image (high res)", im, nmPerPixel); } // Create weight image (including windowing) ImageProcessor w = createWeightImage(im, applyWindow); // Store the area of the image in um^2 weightedAreaInPx = areaInPx = im.getWidth() * im.getHeight(); if (applyWindow) { weightedAreaInPx *= w.getStatistics().mean; } area = areaInPx * nmPerPixel * nmPerPixel / 1e6; weightedArea = weightedAreaInPx * nmPerPixel * nmPerPixel / 1e6; noOfMolecules = molecules.size(); // Pad the images to the largest scale being investigated by the correlation function. // Original Sengupta paper uses 800nm for the padding size. // Limit to within 80% of the minimum dimension of the image. double maxRadius = correlationDistance / nmPerPixel; int imageSize = FastMath.min(im.getWidth(), im.getHeight()); if (imageSize < 1.25 * maxRadius) maxRadius = imageSize / 1.25; int pad = (int) Math.round(maxRadius); log("Analysing up to %.0f nm = %d pixels", maxRadius * nmPerPixel, pad); im = padImage(im, pad); w = padImage(w, pad); // // Used for debugging // { // ImageProcessor w2 = w.duplicate(); // w2.setMinAndMax(0, 1); // PCPALMMolecules.displayImage(PCPALMMolecules.results.getName() + " Binary Image Mask", w2, nmPerPixel); // } final double peakDensity = getDensity(im); // Create 2D auto-correlation double[][] gr; try { // Use the FFT library as it is multi-threaded. This may not be in the user's path. gr = computeAutoCorrelationCurveFFT(im, w, pad, nmPerPixel, peakDensity); } catch (Exception e) { // Default to the ImageJ built-in FHT gr = computeAutoCorrelationCurveFHT(im, w, pad, nmPerPixel, peakDensity); } if (gr == null) return; // Add the g(r) curve to the results addResult(peakDensity, gr); noPlots = WindowManager.getFrame(frequencyDomainTitle) == null; topPlotTitle = spatialPlotTitle; // Do not plot r=0 value on the curve plotCorrelation(gr, 1, frequencyDomainTitle, "g(r)", false, showErrorBars); } if (noPlots) { // Position the plot underneath the other one Frame f1 = WindowManager.getFrame(topPlotTitle); if (f1 != null) { String bottomPlotTitle = (topPlotTitle.equals(spatialPlotTitle) ? frequencyDomainTitle : spatialPlotTitle); Frame f2 = WindowManager.getFrame(bottomPlotTitle); if (f2 != null) f2.setLocation(f2.getLocation().x, f2.getLocation().y + f1.getHeight()); } } log("%s domain analysis computed in %s ms", (spatialDomain) ? "Spatial" : "Frequency", Utils.rounded((System.nanoTime() - start) * 1e-6, 4)); log("---"); } private ImageProcessor applyWindow(ImageProcessor im, ImageWindow imageWindow) { float[] image = (float[]) im.toFloat(0, null).getPixels(); image = imageWindow.applySeperable(image, im.getWidth(), im.getHeight(), ImageWindow.WindowFunction.TUKEY); return new FloatProcessor(im.getWidth(), im.getHeight(), image, null); } public static Plot2 plotCorrelation(double[][] gr, int offset, String plotTitle, String yAxisTitle, boolean barChart, boolean showErrorBars) { double[] x = new double[gr[1].length - offset]; double[] y = new double[x.length]; System.arraycopy(gr[0], offset, x, 0, x.length); System.arraycopy(gr[1], offset, y, 0, y.length); Plot2 plot = new Plot2(plotTitle, "r (nm)", yAxisTitle); plot.setLimits(0, x[x.length - 1], Maths.min(y) * 0.95, Maths.max(y) * 1.05); plot.addPoints(x, y, (barChart) ? Plot2.BAR : Plot.LINE); Utils.display(plotTitle, plot); if (showErrorBars && !barChart) { plot.setColor(Color.magenta); for (int i = 0; i < x.length; i++) { double sd = gr[2][i + offset]; plot.drawLine(x[i], y[i] - sd, x[i], y[i] + sd); } Utils.display(plotTitle, plot); } return plot; } /** * Pad the image by the specified number of pixels * * @param im * @param pad * @return */ private ImageProcessor padImage(ImageProcessor im, int pad) { // int newW = pad * 2 + im.getWidth(); // int newH = pad * 2 + im.getHeight(); int newW = pad + im.getWidth(); int newH = pad + im.getHeight(); ImageProcessor im2 = im.createProcessor(newW, newH); // im2.insert(im, pad, pad); im2.insert(im, 0, 0); return im2; } /** * Create a weight image of the same size. All pixels corresponding to the original image area * are set to 1. A window function is optionally applied. * * @param im * @param applyWindow * @return The weight image */ private ImageProcessor createWeightImage(ImageProcessor im, boolean applyWindow) { float[] data = new float[im.getWidth() * im.getHeight()]; Arrays.fill(data, 1); ImageProcessor w = new FloatProcessor(im.getWidth(), im.getHeight(), data, null); if (applyWindow) { w = applyWindow(w, imageWindow); } return w; } /** * Compute the auto-correlation curve using FHT (ImageJ built-in). Computes the correlation * image and then samples the image at radii up to the specified length to get the average * correlation at a given radius. * * @param im * @param w * @param maxRadius * @param nmPerPixel * @param density * @return { distances[], gr[], gr_se[] } */ private double[][] computeAutoCorrelationCurveFHT(ImageProcessor im, ImageProcessor w, int maxRadius, double nmPerPixel, double density) { log("Creating Hartley transforms"); FHT2 fht2Im = padToFHT2(im); FHT2 fht2W = padToFHT2(w); if (fht2Im == null || fht2W == null) { error("Unable to perform Hartley transform"); return null; } log("Performing correlation"); FloatProcessor corrIm = computeAutoCorrelationFHT(fht2Im); FloatProcessor corrW = computeAutoCorrelationFHT(fht2W); IJ.showProgress(1); final int centre = corrIm.getHeight() / 2; Rectangle crop = new Rectangle(centre - maxRadius, centre - maxRadius, maxRadius * 2, maxRadius * 2); if (showCorrelationImages) { displayCorrelation(corrIm, "Image correlation", crop); displayCorrelation(corrW, "Window correlation", crop); } log("Normalising correlation"); FloatProcessor correlation = normaliseCorrelation(corrIm, corrW, density); if (showCorrelationImages) displayCorrelation(correlation, "Normalised correlation", crop); return computeRadialAverage(maxRadius, nmPerPixel, correlation); } /** * Gets the density of peaks in the image. The density is in squared pixels * * @param im * @return The density (in pixels^-2) */ private double getDensity(ImageProcessor im) { // PCPALMMolecules.densityPeaks is in nm^-2 double density = PCPALMMolecules.densityPeaks * 1e6; // Alternatively use the density in the sample double sampleDensity = (double) noOfMolecules / area; // Actually count the density in the image uniquePoints = 0; for (int i = im.getPixelCount(); i-- > 0;) { // The image may not be binary so use the number uniquePoints += im.getf(i); } double imageDensity = (double) uniquePoints / weightedArea; log(" %d molecules plotted as %.1f unique points", noOfMolecules, uniquePoints); log(" Total Density = %g um^-2, Sample density = %g (%.2fx), Image density = %g (%.2fx)", density, sampleDensity, sampleDensity / density, imageDensity, imageDensity / density); // This is the method used by the PC-PALM MATLAB code. // The sum of the image divided by the sum of the normalisation window function return (double) uniquePoints / weightedAreaInPx; } /** * Compute the auto-correlation curve using FFT. Computes the correlation image and then samples * the image at radii up to the specified length to get the average correlation at a given * radius. * * @param im * @param w * @param maxRadius * @param nmPerPixel * @return { distances[], gr[], gr_se[] } */ private double[][] computeAutoCorrelationCurveFFT(ImageProcessor im, ImageProcessor w, int maxRadius, double nmPerPixel, double density) { log("Performing FFT correlation"); FloatProcessor corrIm = computeAutoCorrelationFFT(im); FloatProcessor corrW = computeAutoCorrelationFFT(w); if (corrIm == null || corrW == null) { error("Unable to perform Fourier transform"); return null; } final int centre = corrIm.getHeight() / 2; Rectangle crop = new Rectangle(centre - maxRadius, centre - maxRadius, maxRadius * 2, maxRadius * 2); if (showCorrelationImages) { displayCorrelation(corrIm, "Image correlation FFT", crop); displayCorrelation(corrW, "Window correlation FFT", crop); } log(" Normalising correlation"); FloatProcessor correlation = normaliseCorrelation(corrIm, corrW, density); if (showCorrelationImages) displayCorrelation(correlation, "Normalised correlation FFT", crop); return computeRadialAverage(maxRadius, nmPerPixel, correlation); } /** * Compute the radial average correlation function (gr) * * @param maxRadius * the maximum radius to process (in pixels) * @param nmPerPixel * covert pixel distances to nm * @param correlation * auto-correlation * @return { distances[], gr[], gr_se[] } */ private double[][] computeRadialAverage(int maxRadius, double nmPerPixel, FloatProcessor correlation) { // Perform averaging of the correlation function using integer distance bins log(" Computing distance vs correlation curve"); final int centre = correlation.getHeight() / 2; // Count the number of pixels at each distance and sum the correlations Statistics[] gr = new Statistics[maxRadius + 1]; // Cache distances int[] d2 = new int[maxRadius + 1]; for (int dy = 0; dy <= maxRadius; dy++) { gr[dy] = new Statistics(); d2[dy] = dy * dy; } int[][] distance = new int[maxRadius + 1][maxRadius + 1]; for (int dy = 0; dy <= maxRadius; dy++) { for (int dx = dy; dx <= maxRadius; dx++) { distance[dy][dx] = distance[dx][dy] = (int) Math.round(Math.sqrt(d2[dx] + d2[dy])); } } float[] data = (float[]) correlation.getPixels(); for (int dy = -maxRadius; dy <= maxRadius; dy++) { final int absY = Math.abs(dy); int index = (centre + dy) * correlation.getWidth() + centre - maxRadius; for (int dx = -maxRadius; dx <= maxRadius; dx++, index++) { final int d = distance[absY][Math.abs(dx)]; if (d > maxRadius || d == 0) continue; gr[d].add(data[index]); } } // Create the final data: a curve showing distance (in nm) verses the average correlation double[] x = new double[maxRadius + 1]; double[] y = new double[maxRadius + 1]; double[] sd = new double[maxRadius + 1]; for (int i = 0; i < x.length; i++) { x[i] = i * nmPerPixel; y[i] = gr[i].getMean(); sd[i] = gr[i].getStandardError(); } y[0] = correlation.getf(centre, centre); // For debugging // double[] H = new double[x.length]; // for (int i = 0; i < x.length; i++) // H[i] = gr[i].getN(); // Plot2 p = new Plot2("Histogram", "r", "F", x, H); // Utils.display("Histogram", p); return new double[][] { x, y, sd }; } private void displayCorrelation(FloatProcessor correlation, String title, Rectangle crop) { correlation.setRoi(crop); ImageProcessor ip = correlation.crop(); ip.resetMinAndMax(); Utils.display(title, ip); } /** * @param fftIm * in frequency domain * @return */ private FloatProcessor computeAutoCorrelationFHT(FHT2 fftIm) { FHT2 FHT2 = fftIm.conjugateMultiply(fftIm); FHT2.inverseTransform(); FHT2.swapQuadrants(); return FHT2; } private FloatProcessor normaliseCorrelation(FloatProcessor corrIm, FloatProcessor corrW, double density) { float[] data = new float[corrIm.getWidth() * corrIm.getHeight()]; float[] dataIm = (float[]) corrIm.getPixels(); float[] dataW = (float[]) corrW.getPixels(); // Square the density for normalisation density *= density; for (int i = 0; i < data.length; i++) { data[i] = (float) ((double) dataIm[i] / (density * dataW[i])); } FloatProcessor correlation = new FloatProcessor(corrIm.getWidth(), corrIm.getHeight(), data, null); return correlation; } /** * Pads the image to the next power of two and transforms into the frequency domain * * @param ip * @return An FHT2 image in the frequency domain */ private FHT2 padToFHT2(ImageProcessor ip) { FloatProcessor im2 = pad(ip); if (im2 == null) return null; FHT2 FHT2 = new FHT2(im2); FHT2.setShowProgress(false); FHT2.transform(); return FHT2; } /** * Pads the image to the next power of two * * @param ip * @return padded image */ private FloatProcessor pad(ImageProcessor ip) { // Pad to a power of 2 final int size = FastMath.max(ip.getWidth(), ip.getHeight()); int newSize = nextPowerOfTwo(size); if (size > newSize) return null; // Error FloatProcessor im2 = new FloatProcessor(newSize, newSize); // If the binary processor has a min and max this breaks the conversion to float // since values are mapped from 0-255 using the min-max look-up calibration table ip.resetMinAndMax(); im2.insert(ip, 0, 0); return im2; } private int nextPowerOfTwo(final int size) { int newSize = 0; for (int i = 4; i < 15; i++) { newSize = (int) Math.pow(2.0, i); if (size <= newSize) { break; } } return newSize; } /** * Compute the auto-correlation using the JTransforms FFT library * * @param ip * @return */ private FloatProcessor computeAutoCorrelationFFT(ImageProcessor ip) { FloatProcessor paddedIp = pad(ip); if (paddedIp == null) return null; final int size = paddedIp.getWidth(); boolean doubleFFT = false; float[] pixels = (float[]) paddedIp.getPixels(); float[] correlation = new float[size * size]; // JTransform library // ------------------ // The data is stored in 1D array in row-major order. Complex number is // stored as two float values in sequence: the real and imaginary part // Correlation = fft^-1(abs(fft).^2) // The absolute value of a complex number z = x + y*i is the value sqrt(x*x+y*y). if (doubleFFT) { DoubleFFT_2D fft = new DoubleFFT_2D(size, size); double[] data = new double[size * size * 2]; for (int i = 0; i < pixels.length; i++) data[i] = pixels[i]; fft.realForwardFull(data); // Re-use data for (int i = 0, j = 0; i < data.length; i += 2, j++) { data[j] = data[i] * data[i] + data[i + 1] * data[i + 1]; } // Zero fill for (int j = correlation.length; j < data.length; j++) data[j] = 0; DoubleFFT_2D fftCorr = new DoubleFFT_2D(size, size); fftCorr.realInverseFull(data, true); // Get the real part of the data for (int i = 0, j = 0; i < data.length; i += 2, j++) { correlation[j] = (float) data[i]; } } else { FloatFFT_2D fft = new FloatFFT_2D(size, size); float[] data = new float[size * size * 2]; System.arraycopy(pixels, 0, data, 0, pixels.length); fft.realForwardFull(data); // Re-use data for (int i = 0, j = 0; i < data.length; i += 2, j++) { data[j] = data[i] * data[i] + data[i + 1] * data[i + 1]; } // Zero fill for (int j = correlation.length; j < data.length; j++) data[j] = 0; FloatFFT_2D fftCorr = new FloatFFT_2D(size, size); fftCorr.realInverseFull(data, true); // Get the real part of the data for (int i = 0, j = 0; i < data.length; i += 2, j++) { correlation[j] = data[i]; } } // Swap quadrants FloatProcessor fp = new FloatProcessor(size, size, correlation, null); new FHT2().swapQuadrants(fp); return fp; } private void addResult(double peakDensity, double[][] gr) { int id = results.size() + 1; // Convert density from pixel^-2 to um^-2 peakDensity *= 1e6 / (nmPerPixel * nmPerPixel); CorrelationResult result = new CorrelationResult(id, PCPALMMolecules.results.getSource(), minx, miny, maxx, maxy, uniquePoints, nmPerPixel, peakDensity, binaryImage, gr, false); results.add(result); createResultsTable(); final double pcw = (PCPALMMolecules.maxx - PCPALMMolecules.minx) / 100.0; final double pch = (PCPALMMolecules.maxy - PCPALMMolecules.miny) / 100.0; StringBuilder sb = new StringBuilder(); sb.append(id).append("\t"); sb.append(PCPALMMolecules.results.getName()).append("\t"); sb.append(IJ.d2s(minx)).append("\t"); sb.append(IJ.d2s((minx) / pcw)).append("\t"); sb.append(IJ.d2s(miny)).append("\t"); sb.append(IJ.d2s((miny) / pch)).append("\t"); sb.append(IJ.d2s(maxx - minx)).append("\t"); sb.append(IJ.d2s((maxx - minx) / pcw)).append("\t"); sb.append(IJ.d2s(maxy - miny)).append("\t"); sb.append(IJ.d2s((maxy - miny) / pch)).append("\t"); sb.append(Utils.rounded(uniquePoints, 4)).append("\t"); sb.append(Utils.rounded(peakDensity, 4)).append("\t"); sb.append(Utils.rounded(nmPerPixel, 4)).append("\t"); sb.append(binaryImage).append("\t"); resultsTable.append(sb.toString()); } private void createResultsTable() { if (resultsTable == null || !resultsTable.isVisible()) { StringBuilder sb = new StringBuilder(); sb.append("ID\t"); sb.append("Image Source\t"); sb.append("X\t"); sb.append("X %\t"); sb.append("Y\t"); sb.append("Y %\t"); sb.append("Width\t"); sb.append("Width %\t"); sb.append("Height\t"); sb.append("Height %\t"); sb.append("N\t"); sb.append("PeakDensity (um^-2)\t"); sb.append("nm/pixel\t"); sb.append("Binary\t"); resultsTable = new TextWindow(TITLE, sb.toString(), (String) null, 800, 300); } } }