diff --git a/src/main/java/ColocImgLibGadgets.java b/src/main/java/ColocImgLibGadgets.java index 4d52124..2ba7eb4 100644 --- a/src/main/java/ColocImgLibGadgets.java +++ b/src/main/java/ColocImgLibGadgets.java @@ -1,131 +1,132 @@ import ij.IJ; import ij.ImagePlus; import ij.plugin.PlugIn; import java.util.ArrayList; import java.util.Collections; import java.util.List; import java.util.Random; import net.imglib2.Cursor; import net.imglib2.algorithm.math.ImageStatistics; import net.imglib2.img.ImagePlusAdapter; import net.imglib2.img.Img; import net.imglib2.img.ImgFactory; import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.type.NativeType; import net.imglib2.type.numeric.RealType; public class ColocImgLibGadgets & NativeType> implements PlugIn { protected Img img1, img2; + @Override public void run(String arg) { ImagePlus imp1 = IJ.openImage("/Users/dan/Documents/Dresden/ipf/colocPluginDesign/red.tif"); img1 = ImagePlusAdapter.wrap(imp1); ImagePlus imp2 = IJ.openImage("/Users/dan/Documents/Dresden/ipf/colocPluginDesign/green.tif"); img2 = ImagePlusAdapter.wrap(imp2); double pearson = calculatePearson(); Img ranImg = generateRandomImageStack(img1, new int[] {2,2,1}); } /** * To randomize blockwise we enumerate the blocks, shuffle that list and * write the data to their new position based on the shuffled list. */ protected Img generateRandomImageStack(Img img, int[] blockDimensions) { int numberOfDimensions = Math.min(img.numDimensions(), blockDimensions.length); int numberOfBlocks = 0; long[] numberOfBlocksPerDimension = new long[numberOfDimensions]; for (int i = 0 ; i allTheBlocks = new ArrayList(numberOfBlocks); for (int i = 0; i cursor = img.cursor(); // create factories for new image stack //ContainerFactory containerFactory = new ImagePlusContainerFactory(); ImgFactory imgFactory = new ArrayImgFactory(); //new ImageFactory(cursor.getType(), containerFactory); // create a new stack for the random images final long[] dim = new long[ img.numDimensions() ]; img.dimensions(dim); Img randomStack = imgFactory.create(dim, img.firstElement().createVariable()); // iterate over image data while (cursor.hasNext()) { cursor.fwd(); T type = cursor.get(); // type.getRealDouble(); } return randomStack; } protected double calculatePearson() { Cursor cursor1 = img1.cursor(); Cursor cursor2 = img2.cursor(); double mean1 = getImageMean(img1); double mean2 = getImageMean(img2); // Do some rather simple performance testing long startTime = System.currentTimeMillis(); double pearson = calculatePearson(cursor1, mean1, cursor2, mean2); // End performance testing long finishTime = System.currentTimeMillis(); long elapsed = finishTime - startTime; // print some output to IJ log IJ.log("mean of ch1: " + mean1 + " " + "mean of ch2: " + mean2); IJ.log("Pearson's Coefficient " + pearson); IJ.log("That took: " + elapsed + " ms"); return pearson; } protected double calculatePearson(Cursor cursor1, double mean1, Cursor cursor2, double mean2) { double pearsonDenominator = 0; double ch1diffSquaredSum = 0; double ch2diffSquaredSum = 0; while (cursor1.hasNext() && cursor2.hasNext()) { cursor1.fwd(); cursor2.fwd(); T type1 = cursor1.get(); double ch1diff = type1.getRealDouble() - mean1; T type2 = cursor2.get(); double ch2diff = type2.getRealDouble() - mean2; pearsonDenominator += ch1diff*ch2diff; ch1diffSquaredSum += (ch1diff*ch1diff); ch2diffSquaredSum += (ch2diff*ch2diff); } double pearsonNumerator = Math.sqrt(ch1diffSquaredSum * ch2diffSquaredSum); return pearsonDenominator / pearsonNumerator; } protected double getImageMean(Img img) { double sum = 0; Cursor cursor = img.cursor(); while (cursor.hasNext()) { cursor.fwd(); T type = cursor.get(); sum += type.getRealDouble(); } return sum / ImageStatistics.getNumPixels(img); } } diff --git a/src/main/java/Coloc_2.java b/src/main/java/Coloc_2.java index 817f1a4..1ee843a 100644 --- a/src/main/java/Coloc_2.java +++ b/src/main/java/Coloc_2.java @@ -1,753 +1,754 @@ import algorithms.Algorithm; import algorithms.AutoThresholdRegression; import algorithms.AutoThresholdRegression.Implementation; import algorithms.CostesSignificanceTest; import algorithms.Histogram2D; import algorithms.InputCheck; import algorithms.KendallTauRankCorrelation; import algorithms.LiHistogram2D; import algorithms.LiICQ; import algorithms.MandersColocalization; import algorithms.MissingPreconditionException; import algorithms.PearsonsCorrelation; import algorithms.SpearmanRankCorrelation; import fiji.Debug; import gadgets.DataContainer; import ij.IJ; import ij.ImagePlus; import ij.Prefs; import ij.WindowManager; import ij.gui.GenericDialog; import ij.gui.Roi; import ij.gui.ShapeRoi; import ij.plugin.PlugIn; import ij.plugin.frame.RoiManager; import ij.process.Blitter; import ij.process.ImageProcessor; import java.awt.Checkbox; import java.awt.Frame; import java.awt.Rectangle; import java.awt.event.ItemEvent; import java.awt.event.ItemListener; import java.awt.event.WindowAdapter; import java.awt.event.WindowEvent; import java.util.ArrayList; import java.util.Arrays; import java.util.List; import net.imglib2.Cursor; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.img.ImagePlusAdapter; import net.imglib2.img.Img; import net.imglib2.img.ImgFactory; import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.type.NativeType; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; import results.PDFWriter; import results.ResultHandler; import results.SingleWindowDisplay; import results.Warning; /** Copyright 2010-2015, Daniel J. White, Tom Kazimiers, Johannes Schindelin and the Fiji project. Fiji is just imageJ - batteries included. 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. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/ . */ /** * An ImageJ plugin which does pixel intensity correlation based * colocalisation analysis on a pair of images, * with optional Mask or ROI. * * @param * @author Daniel J. White * @author Tom Kazimiers * @author Johannes Schindelin */ public class Coloc_2 & NativeType< T >> implements PlugIn { // a small bounding box container protected class BoundingBox { public long[] offset; public long[] size; public BoundingBox(long [] offset, long[] size) { this.offset = offset.clone(); this.size = size.clone(); } } // a storage class for ROI information protected class MaskInfo { BoundingBox roi; public Img mask; // constructors public MaskInfo(BoundingBox roi, Img mask) { this.roi = roi; this.mask = mask; } public MaskInfo() { } } // the storage key for Fiji preferences protected final static String PREF_KEY = "Coloc_2."; // Allowed types of ROI configuration protected enum RoiConfiguration { None, Img1, Img2, Mask, RoiManager }; // the ROI configuration to use protected RoiConfiguration roiConfig = RoiConfiguration.Img1; // A list of all ROIs/masks found protected ArrayList masks = new ArrayList(); // A list of auto threshold implementations protected String[] regressions = new String[ AutoThresholdRegression.Implementation.values().length]; // default indices of image, mask, ROI and regression choices protected static int index1 = 0; protected static int index2 = 1; protected static int indexMask = 0; protected static int indexRoi = 0; protected static int indexRegr = 0; // the images to work on protected Img img1, img2; // names of the images working on protected String Ch1Name = ""; protected String Ch2Name = ""; // the channels of the images to use protected int img1Channel = 1, img2Channel = 1; /* The different algorithms this plug-in provides. * If a reference is null it will not get run. */ protected PearsonsCorrelation pearsonsCorrelation; protected LiHistogram2D liHistogramCh1; protected LiHistogram2D liHistogramCh2; protected LiICQ liICQ; protected SpearmanRankCorrelation SpearmanRankCorrelation; protected MandersColocalization mandersCorrelation; protected KendallTauRankCorrelation kendallTau; protected Histogram2D histogram2D; protected CostesSignificanceTest costesSignificance; // indicates if images should be printed in result protected boolean displayImages; // indicates if a PDF should be saved automatically protected boolean autoSavePdf; + @Override public void run(String arg0) { if (showDialog()) { try { for (MaskInfo mi : masks) { colocalise(img1, img2, mi.roi, mi.mask); } } catch (MissingPreconditionException e) { IJ.handleException(e); IJ.showMessage("An error occured, could not colocalize!"); return; } } } public boolean showDialog() { // get IDs of open windows int[] windowList = WindowManager.getIDList(); // if there are less than 2 windows open, cancel if (windowList == null || windowList.length < 2) { IJ.showMessage("At least 2 images must be open!"); return false; } /* create a new generic dialog for the * display of various options. */ final GenericDialog gd = new GenericDialog("Coloc 2"); String[] titles = new String[windowList.length]; /* the masks and ROIs array needs three more entries than * windows to contain "none", "ROI ch 1" and "ROI ch 2" */ String[] roisAndMasks= new String[windowList.length + 4]; roisAndMasks[0]=""; roisAndMasks[1]="ROI(s) in channel 1"; roisAndMasks[2]="ROI(s) in channel 2"; roisAndMasks[3]="ROI Manager"; // go through all open images and add them to GUI for (int i=0; i < windowList.length; i++) { ImagePlus imp = WindowManager.getImage(windowList[i]); if (imp != null) { titles[i] = imp.getTitle(); roisAndMasks[i + 4] =imp.getTitle(); } else { titles[i] = ""; } } // find all available regression strategies Implementation[] regressionImplementations = AutoThresholdRegression.Implementation.values(); for (int i=0; i maskImg = ImagePlusAdapter.wrap( maskImp ); // get a valid mask info for the image MaskInfo mi = getBoundingBoxOfMask(maskImg); masks.add( mi ) ; } else { /* if no ROI/mask is selected, just add an empty MaskInfo * to colocalise both images without constraints. */ masks.add(new MaskInfo(null, null)); } // get information about the mask/ROI to use indexRegr = gd.getNextChoiceIndex(); // read out GUI data autoSavePdf = gd.getNextBoolean(); displayImages = gd.getNextBoolean(); displayShuffledCostes = gd.getNextBoolean(); useLiCh1 = gd.getNextBoolean(); useLiCh2 = gd.getNextBoolean(); useLiICQ = gd.getNextBoolean(); useSpearmanRank = gd.getNextBoolean(); useManders = gd.getNextBoolean(); useKendallTau = gd.getNextBoolean(); useScatterplot = gd.getNextBoolean(); useCostes = gd.getNextBoolean(); psf = (int) gd.getNextNumber(); nrCostesRandomisations = (int) gd.getNextNumber(); // save user preferences Prefs.set(PREF_KEY+"regressionImplementation", indexRegr); Prefs.set(PREF_KEY+"autoSavePdf", autoSavePdf); Prefs.set(PREF_KEY+"displayImages", displayImages); Prefs.set(PREF_KEY+"displayShuffledCostes", displayShuffledCostes); Prefs.set(PREF_KEY+"useLiCh1", useLiCh1); Prefs.set(PREF_KEY+"useLiCh2", useLiCh2); Prefs.set(PREF_KEY+"useLiICQ", useLiICQ); Prefs.set(PREF_KEY+"useSpearmanRank", useSpearmanRank); Prefs.set(PREF_KEY+"useManders", useManders); Prefs.set(PREF_KEY+"useKendallTau", useKendallTau); Prefs.set(PREF_KEY+"useScatterplot", useScatterplot); Prefs.set(PREF_KEY+"useCostes", useCostes); Prefs.set(PREF_KEY+"psf", psf); Prefs.set(PREF_KEY+"nrCostesRandomisations", nrCostesRandomisations); // Parse algorithm options pearsonsCorrelation = new PearsonsCorrelation(PearsonsCorrelation.Implementation.Fast); if (useLiCh1) liHistogramCh1 = new LiHistogram2D("Li - Ch1", true); if (useLiCh2) liHistogramCh2 = new LiHistogram2D("Li - Ch2", false); if (useLiICQ) liICQ = new LiICQ(); if (useSpearmanRank) SpearmanRankCorrelation = new SpearmanRankCorrelation(); if (useManders) mandersCorrelation = new MandersColocalization(); if (useKendallTau) kendallTau = new KendallTauRankCorrelation(); if (useScatterplot) histogram2D = new Histogram2D("2D intensity histogram"); if (useCostes) { costesSignificance = new CostesSignificanceTest(pearsonsCorrelation, psf, nrCostesRandomisations, displayShuffledCostes); } return true; } /** * Call this method to run a whole colocalisation configuration, * all selected algorithms get run on the supplied images. You * can specify the data further by supplying appropriate * information in the mask structure. * * @param img1 * @param img2 * @param roi * @param mask * @param maskBB * @throws MissingPreconditionException */ public void colocalise(Img img1, Img img2, BoundingBox roi, Img mask) throws MissingPreconditionException { // create a new container for the selected images and channels DataContainer container; if (mask != null) { container = new DataContainer(img1, img2, img1Channel, img2Channel, Ch1Name, Ch2Name, mask, roi.offset, roi.size); } else if (roi != null) { // we have no mask, but a regular ROI in use container = new DataContainer(img1, img2, img1Channel, img2Channel, Ch1Name, Ch2Name, roi.offset, roi.size); } else { // no mask and no ROI is present container = new DataContainer(img1, img2, img1Channel, img2Channel, Ch1Name, Ch2Name); } // create a results handler final List> listOfResultHandlers = new ArrayList>(); final PDFWriter pdfWriter = new PDFWriter(container); final SingleWindowDisplay swDisplay = new SingleWindowDisplay(container, pdfWriter); listOfResultHandlers.add(swDisplay); listOfResultHandlers.add(pdfWriter); //ResultHandler resultHandler = new EasyDisplay(container); // this list contains the algorithms that will be run when the user clicks ok List> userSelectedJobs = new ArrayList>(); // add some pre-processing jobs: userSelectedJobs.add( container.setInputCheck( new InputCheck()) ); userSelectedJobs.add( container.setAutoThreshold( new AutoThresholdRegression(pearsonsCorrelation, AutoThresholdRegression.Implementation.values()[indexRegr]))); // add user selected algorithms addIfValid(pearsonsCorrelation, userSelectedJobs); addIfValid(liHistogramCh1, userSelectedJobs); addIfValid(liHistogramCh2, userSelectedJobs); addIfValid(liICQ, userSelectedJobs); addIfValid(SpearmanRankCorrelation, userSelectedJobs); addIfValid(mandersCorrelation, userSelectedJobs); addIfValid(kendallTau, userSelectedJobs); addIfValid(histogram2D, userSelectedJobs); addIfValid(costesSignificance, userSelectedJobs); // execute all algorithms int count = 0; int jobs = userSelectedJobs.size(); for (Algorithm a : userSelectedJobs){ try { count++; IJ.showStatus(count + "/" + jobs + ": Running " + a.getName()); a.execute(container); } catch (MissingPreconditionException e){ for (ResultHandler r : listOfResultHandlers){ r.handleWarning( new Warning( "Probem with input data", a.getName() + ": " + e.getMessage() ) ); } } } // clear status IJ.showStatus(""); // let the algorithms feed their results to the handler for (Algorithm a : userSelectedJobs){ for (ResultHandler r : listOfResultHandlers) a.processResults(r); } // if we have ROIs/masks, add them to results if (displayImages) { RandomAccessibleInterval channel1, channel2; if (mask != null || roi != null) { long[] offset = container.getMaskBBOffset(); long[] size = container.getMaskBBSize(); channel1 = createMaskImage( container.getSourceImage1(), container.getMask(), offset, size ); channel2 = createMaskImage( container.getSourceImage2(), container.getMask(), offset, size ); } else { channel1 = container.getSourceImage1(); channel2 = container.getSourceImage2(); } for (ResultHandler r : listOfResultHandlers) { r.handleImage (channel1, "Channel 1"); r.handleImage (channel2, "Channel 2"); } } // do the actual results processing swDisplay.process(); // add window to the IJ window manager swDisplay.addWindowListener(new WindowAdapter() { @Override public void windowClosed(WindowEvent e) { WindowManager.removeWindow((Frame) swDisplay); } }); WindowManager.addWindow(swDisplay); // show PDF saving dialog if requested if (autoSavePdf) pdfWriter.process(); } /** * A method to get the bounding box from the data in the given * image that is above zero. Those values are interpreted as a * mask. It will return null if no mask information was found. * * @param mask The image to look for "on" values in * @return a new MaskInfo object or null */ protected MaskInfo getBoundingBoxOfMask(Img mask) { Cursor cursor = mask.localizingCursor(); int numMaskDims = mask.numDimensions(); // the "off type" of the mask T offType = mask.firstElement().createVariable(); offType.setZero(); // the corners of the bounding box long[] min = null; long[] max = null; // indicates if mask data has been found boolean maskFound = false; // a container for temporary position information long[] pos = new long[numMaskDims]; // walk over the mask while (cursor.hasNext() ) { cursor.fwd(); T data = cursor.get(); // test if the current mask data represents on or off if (data.compareTo(offType) > 0) { // get current position cursor.localize(pos); if (!maskFound) { // we found mask data, first time maskFound = true; // init min and max with the current position min = Arrays.copyOf(pos, numMaskDims); max = Arrays.copyOf(pos, numMaskDims); } else { /* Is is at least second hit, compare if it * has new "extreme" positions, i.e. does * is make the BB bigger? */ for (int d=0; d max[d]) { // is it larger than max max[d] = pos[d]; } } } } } if (!maskFound) { return null; } else { // calculate size long[] size = new long[numMaskDims]; for (int d=0; d a, List> list) { if (a != null) list.add(a); } /** * Returns true if a custom ROI has been selected, i.e if the current * ROI does not have the extent of the whole image. * @return true if custom ROI selected, false otherwise */ protected boolean hasValidRoi(ImagePlus imp) { Roi roi = imp.getRoi(); if (roi == null) return false; Rectangle theROI = roi.getBounds(); // if the ROI is the same size as the image (default ROI), return false return (theROI.height != imp.getHeight() || theROI.width != imp.getWidth()); } /** * Clips a value to the specified bounds. */ protected static int clip(int val, int min, int max) { return Math.max( Math.min( val, max ), min ); } /** * This method checks if the given ImagePlus contains any * masks or ROIs. If so, the appropriate date structures * are created and filled. */ protected void createMasksFromImage(ImagePlus imp) { // get ROIs from current image in Fiji Roi[] impRois = split(imp.getRoi()); // create the ROIs createMasksAndRois(impRois, imp.getWidth(), imp.getHeight()); } /** * A method to fill the masks array with data based on the ROI manager. */ protected boolean createMasksFromRoiManager(int width, int height) { RoiManager roiManager = RoiManager.getInstance(); if (roiManager == null) { IJ.error("Could not get ROI Manager instance."); return false; } Roi[] selectedRois = roiManager.getSelectedRoisAsArray(); // create the ROIs createMasksAndRois(selectedRois, width, height); return true; } /** * Creates appropriate data structures from the ROI information * passed. If an irregular ROI is found, it will be put into a * frame of its bounding box size and put into an Image. * * In the end the members ROIs, masks and maskBBs will be * filled if ROIs or masks were found. They will be null * otherwise. */ protected void createMasksAndRois(Roi[] rois, int width, int height) { // create empty list masks.clear(); for (Roi r : rois ){ MaskInfo mi = new MaskInfo(); // add it to the list of masks/ROIs masks.add(mi); // get the ROIs/masks bounding box Rectangle rect = r.getBounds(); mi.roi = new BoundingBox( new long[] {rect.x, rect.y} , new long[] {rect.width, rect.height}); ImageProcessor ipMask = r.getMask(); // check if we got a regular ROI and return if so if (ipMask == null) { continue; } // create a mask processor of the same size as a slice ImageProcessor ipSlice = ipMask.createProcessor(width, height); // fill the new slice with black ipSlice.setValue(0.0); ipSlice.fill(); // position the mask on the new mask processor ipSlice.copyBits(ipMask, (int)mi.roi.offset[0], (int)mi.roi.offset[1], Blitter.COPY); // create an Image out of it ImagePlus maskImp = new ImagePlus("Mask", ipSlice); // and remember it and the masks bounding box mi.mask = ImagePlusAdapter.wrap( maskImp ); } } /** * This method duplicates the given images, but respects * ROIs if present. Meaning, a sub-picture will be created when * source images are ROI/MaskImages. * @throws MissingPreconditionException */ protected RandomAccessibleInterval createMaskImage( RandomAccessibleInterval image, RandomAccessibleInterval mask, long[] offset, long[] size) throws MissingPreconditionException { long[] pos = new long[ image.numDimensions() ]; // sanity check if (pos.length != offset.length || pos.length != size.length) { throw new MissingPreconditionException("Mask offset and size must be of same dimensionality like image."); } // use twin cursor for only one image TwinCursor cursor = new TwinCursor( image.randomAccess(), image.randomAccess(), Views.iterable(mask).localizingCursor()); // prepare output image ImgFactory maskFactory = new ArrayImgFactory(); //Img maskImage = maskFactory.create( size, name ); RandomAccessibleInterval maskImage = maskFactory.create( size, image.randomAccess().get().createVariable() ); RandomAccess maskCursor = maskImage.randomAccess(); // go through the visible data and copy it to the output while (cursor.hasNext()) { cursor.fwd(); cursor.localize(pos); // shift coordinates by offset for (int i=0; i < pos.length; ++i) { pos[i] = pos[i] - offset[i]; } // write out to correct position maskCursor.setPosition( pos ); maskCursor.get().set( cursor.getFirst() ); } return maskImage; } /** * Splits a non overlapping composite ROI into its sub ROIs. * * @param roi The ROI to split * @return A list of one or more ROIs */ public static Roi[] split(Roi roi) { if (roi instanceof ShapeRoi) return ((ShapeRoi)roi).getRois(); return new Roi[] { roi }; } /** * Main method for easier development. To run this plugin with Maven, use: * mvn exec:java -Dexec.mainClass="Coloc_2" */ public static void main(String[] args) { Debug.run(null, null); } } diff --git a/src/main/java/Colocalisation_Test.java b/src/main/java/Colocalisation_Test.java index 4e7fcdc..c75eab0 100644 --- a/src/main/java/Colocalisation_Test.java +++ b/src/main/java/Colocalisation_Test.java @@ -1,741 +1,742 @@ //version 21 4 05 //added van Steensel CCF analysis //JCellSci v109.p787 //version 29/4/05 //Costes randomisation uses pixels in ch2 only once per random image import ij.IJ; import ij.ImagePlus; import ij.ImageStack; import ij.Prefs; import ij.WindowManager; import ij.gui.GenericDialog; import ij.gui.PlotWindow; import ij.gui.Roi; import ij.measure.Calibration; import ij.plugin.PlugIn; import ij.plugin.filter.GaussianBlur; import ij.process.ByteProcessor; import ij.process.ImageProcessor; import ij.process.ImageStatistics; import ij.process.ShortProcessor; import ij.text.TextWindow; import java.awt.Rectangle; import java.text.DecimalFormat; public class Colocalisation_Test implements PlugIn {static boolean headingsSet2; private static int index1=0; private static int index2=1; private ImagePlus imp1, imp2,impmask, imp3; private ImageProcessor ipmask; private int indexRoi= (int)Prefs.get("Rand_indexRoi.int",0); private int indexRand= (int)Prefs.get("Rand_indexRand.int",0); private boolean useROI, useMask ; private Roi roi, roi1, roi2; private static boolean randZ= Prefs.get("Rand_randZ.boolean", false); private static boolean ignoreZeroZero= Prefs.get("Rand_ignore.boolean", true); private static boolean smooth= Prefs.get("Rand_smooth.boolean", true); private static boolean keep = Prefs.get("Rand_keep.boolean", false); private static boolean currentSlice= Prefs.get("Rand_currentSlice.boolean", true); private static boolean useManPSF= Prefs.get("Rand_useManPSF.boolean", false); private static boolean showR= Prefs.get("Rand_showR.boolean", false); private static double psf =0; private static int manPSF= (int)Prefs.get("Rand_manPSF.int",10); private static int iterations= (int)Prefs.get("Rand_iterations.int",0); private static double ch2Lambda = Prefs.get("Rand_ch2L.double",520); private static double NA = Prefs.get("Rand_NA.double",1.4); private static double pixelSize = Prefs.get("Rand_pixelSize.double",0.10); DecimalFormat df3 = new DecimalFormat("##0.000"); DecimalFormat df2 = new DecimalFormat("##0.00"); DecimalFormat df1 = new DecimalFormat("##0.0"); DecimalFormat df0 = new DecimalFormat("##0"); String[] chooseRand= { "Fay (x,y,z translation)","Costes approximation (smoothed noise)", "van Steensel (x translation)"}; private int width, height, rwidth, rheight, xOffset, yOffset, mask; String randMethod = "Fay"; private long startTime; private long endTime; StringBuffer rVals = new StringBuffer(); boolean Costes = false; boolean Fay = false; boolean vanS = false; boolean rBlocks= false; protected static TextWindow textWindow; + @Override public void run(String arg) {startTime = System.currentTimeMillis(); if (showDialog()) correlate(imp1, imp2, imp3); } public boolean showDialog() { int[] wList = WindowManager.getIDList(); if (wList==null) { IJ.noImage(); return false; } String[] titles = new String[wList.length]; String[] chooseROI= new String[wList.length+3]; chooseROI[0] = "None"; chooseROI[1] = "ROI in channel 1 "; chooseROI[2] = "ROI in channel 2"; if (indexRoi>wList.length+3) indexRoi=0; for (int i=0; i=titles.length)index1 = 0; if (index2>=titles.length)index2 = 0; GenericDialog gd = new GenericDialog("Colocalisation Test"); gd.addChoice("Channel 1", titles, titles[index1]); gd.addChoice("Channel 2", titles, titles[index2]); gd.addChoice("ROI or Mask", chooseROI, chooseROI[indexRoi]); gd.addChoice("Randomization method", chooseRand,chooseRand[indexRand]); //gd.addCheckbox("Ignore zero-zero pixels", ignoreZeroZero); gd.addCheckbox("Current slice only (Ch1)", currentSlice); gd.addCheckbox("Keep example randomized image", keep); gd.addCheckbox("Show all R values from Ch1 vs Ch2(rand)", showR); gd.addMessage("See: http://uhnresearch.ca/wcif/imagej"); gd.showDialog(); if (gd.wasCanceled()) return false; index1 = gd.getNextChoiceIndex(); index2 = gd.getNextChoiceIndex(); indexRoi = gd.getNextChoiceIndex(); indexRand=gd.getNextChoiceIndex(); ignoreZeroZero = true; currentSlice= gd.getNextBoolean(); keep = gd.getNextBoolean(); showR=gd.getNextBoolean(); String title1 = titles[index1]; String title2 = titles[index2]; imp1 = WindowManager.getImage(wList[index1]); imp2 = WindowManager.getImage(wList[index2]); if (imp1.getType()==imp1.COLOR_RGB || imp2.getType()==imp1.COLOR_RGB) { IJ.showMessage("Colocalisation Test", "Both images must be grayscale."); return false; } useMask=false; if (indexRoi >=3) { imp3 = WindowManager.getImage(indexRoi-2); useMask=true; } else imp3 = WindowManager.getImage(wList[index2]); Calibration cal = imp2.getCalibration(); pixelSize = cal.pixelWidth; if(indexRand==0) {Fay=true; randMethod = "Fay"; } if(indexRand==1) {Costes = true; randMethod = "Costes X, Y"; } if(indexRand==2) {vanS=true; randMethod = "van Steensel"; } //test to ensure all images match up. boolean matchWidth=false; boolean matchHeight=false; boolean matchSlice = false; if (imp1.getWidth()==imp2.getWidth()&&imp1.getWidth()==imp3.getWidth()) matchWidth = true; if (imp1.getHeight()==imp2.getHeight()&&imp1.getHeight()==imp3.getHeight()) matchHeight = true; if (imp1.getStackSize()==imp2.getStackSize()&&imp1.getStackSize()==imp3.getStackSize()) matchSlice = true; if (!(matchWidth&&matchHeight&&matchSlice)) {IJ.showMessage("Image mismatch","Images do not match. Exiting"); return false; } if (Costes||rBlocks) { GenericDialog gd2 = new GenericDialog("PSF details"); gd2.addCheckbox("Randomize pixels in z-axis", randZ); gd2.addNumericField("Pixel Size (µm)", pixelSize,3); gd2.addNumericField("Channel 2 wavelength (nm)", ch2Lambda,0); gd2.addNumericField("NA of objective", NA,2); gd2.addNumericField("Iterations",iterations,0); gd2.addMessage(""); gd2.addCheckbox("Use manual PSF", useManPSF); gd2.addNumericField("PSF radius in pixels", manPSF, 0); gd2.showDialog(); if (gd2.wasCanceled()) return false; randZ = gd2.getNextBoolean(); if (randZ) randMethod+=", Z"; pixelSize =gd2.getNextNumber(); ch2Lambda = gd2.getNextNumber(); NA = gd2.getNextNumber(); iterations = (int)gd2.getNextNumber(); useManPSF = gd2.getNextBoolean(); manPSF = (int)gd2.getNextNumber(); psf = (0.61*ch2Lambda)/NA; psf = (psf)/(pixelSize*1000); if (useManPSF) psf = manPSF; //IJ.showMessage("PSF radius = "+df3.format(psf)); } return true; } public void correlate(ImagePlus imp1, ImagePlus imp2, ImagePlus imp3) { String Ch1fileName = imp1.getTitle(); String Ch2fileName = imp2.getTitle(); ImageStack img1 = imp1.getStack(); ImageStack img2 = imp2.getStack(); ImageStack img3=imp3.getStack(); int width = imp1.getWidth(); int height = imp1.getHeight(); String fileName = Ch1fileName + " and " + Ch2fileName; ImageProcessor ip1 = imp1.getProcessor(); ImageProcessor ip2 = imp2.getProcessor(); ImageProcessor ip3= null; if (useMask) ip3 = imp3.getProcessor(); ImageStack stackRand = new ImageStack(rwidth,rheight); ImageProcessor ipRand= img2.getProcessor(1); int currentSliceNo = imp1.getCurrentSlice(); if(currentSlice) fileName = fileName+ ". slice: "+currentSliceNo ; double pearsons1 = 0; double pearsons2 = 0; double pearsons3 = 0; double r2= 1; double r=1; double ch1Max=0; double ch1Min = ip1.getMax(); double ch2Max=0; double ch2Min = ip1.getMax(); int nslices = imp1.getStackSize(); int ch1, ch2, count; double sumX = 0; double sumXY = 0; double sumXX = 0; double sumYY = 0; double sumY = 0; double sumXtotal=0; double sumYtotal=0; double colocX=0; double colocY=0; int N = 0; int N2 = 0; double r2min=1; double r2max=-1; sumX = 0; sumXY = 0; sumXX = 0; sumYY = 0; sumY = 0; int i=1; double coloc2M1 = 0; double coloc2M2 = 0; int colocCount=0; int colocCount1=0; int colocCount2=0; double r2sd=0; double sumr2sqrd=0; double sumr2=0; double ICQ2mean=0; double sumICQ2sqrd=0; double sumICQ2=0; double ICQobs=0; int countICQ=0; int Nr=0; int Ng=0; //get stack2 histogram ImageStatistics stats = imp2.getStatistics(); if (imp2.getType() == imp2.GRAY16) stats.nBins = 1<<16; int [] histogram = new int[stats.nBins]; //roi code if (indexRoi==1||indexRoi==2) useROI = true; else useROI=false; ip1 = imp1.getProcessor(); ip2 = imp2.getProcessor(); ip3 = imp2.getProcessor(); if (useMask) ip3 = imp3.getProcessor(); roi1 = imp1.getRoi(); roi2 = imp2.getRoi(); Rectangle rect =ip1.getRoi(); if (indexRoi==1) {if(roi1==null) useROI=false; else { ipmask = imp1.getMask(); rect = ip1.getRoi(); } } if (indexRoi==2) {if(roi2==null) useROI=false; else{ipmask = imp2.getMask(); rect = ip2.getRoi();} } if (useROI==false) {xOffset = 0;yOffset = 0; rwidth=width; rheight =height;} else {xOffset = rect.x; yOffset = rect.y; rwidth=rect.width; rheight =rect.height;} int g1=0;int g2=0; int histCount=0; //calulate pearsons for existing image; for (int s=1; s<=nslices;s++) {if (currentSlice) {s=currentSliceNo; nslices=s; } ip1 = img1.getProcessor(s); ip2 = img2.getProcessor(s); ip3 = img3.getProcessor(s); for (int y=0; ych1) ch1Min = ch1; if (ch2Maxch2) ch2Min = ch2; N++; if(ch1+ch2!=0) N2++; sumXtotal += ch1; sumYtotal += ch2; if(ch2>0) colocX += ch1; if(ch1>0) colocY += ch2; sumX +=ch1; sumXY += (ch1 * ch2); sumXX += (ch1 * ch1); sumYY += (ch2 *ch2); sumY += ch2; if(ch1>0) Nr++; if(ch2>0)Ng++; //add ch2 value to histogram histCount = histogram[ch2]; histCount++; histogram[ch2]=histCount; } } } } if (ignoreZeroZero) N = N2; //N = N2; //double ch1Mean = sumX/Nr; //double ch2Mean = sumY/Ng; double ch1Mean = sumX/N; double ch2Mean = sumY/N; // IJ.showMessage("Ch1: "+ch1Mean +" Ch2: "+ch2Mean +" count nonzerozero: "+N); pearsons1 = sumXY - (sumX*sumY/N); pearsons2 = sumXX - (sumX*sumX/N); pearsons3 = sumYY - (sumY*sumY/N); //IJ.showMessage("p1: "+pearsons1+" p2: "+pearsons2+" p3: "+pearsons3); r= pearsons1/(Math.sqrt(pearsons2*pearsons3)); double colocM1 = (double)(colocX/sumXtotal); double colocM2 = (double)(colocY/sumYtotal); //calucalte ICQ int countAll=0; int countPos=0; double PDMobs=0; double PDM=0; double ICQ2; int countAll2=0; for (int s=1; s<=nslices;s++) {if (currentSlice) {s=currentSliceNo; nslices=s; } ip1 = img1.getProcessor(s); ip2 = img2.getProcessor(s); ip3 = img3.getProcessor(s); for (int y=0; y<=rheight; y++) {IJ.showStatus("Calculating r for original images. Press 'Esc' to abort"); if (IJ.escapePressed()) {IJ.beep(); return;} for (int x=0; x<=rwidth; x++) {mask = (int)ip3.getPixelValue(x,y); if (indexRoi==0) mask=1; if((useROI)&&(ipmask!=null)) mask = (int)ipmask.getPixelValue(x,y); if (mask!=0) { ch1 = (int)ip1.getPixel(x+xOffset,y+yOffset); ch2 = (int)ip2.getPixel(x+xOffset,y+yOffset); if (ch1+ch2!=0) { PDMobs = ((double)ch1-(double)ch1Mean)*((double)ch2-(double)ch2Mean); if (PDMobs>0) countPos++; countAll2++; } } } } } //IJ.showMessage("count+ = "+countPos +" CountNonZeroPair= "+countAll2); ICQobs = ((double)countPos/(double)countAll2)-0.5; boolean ch3found= false; //do random localisations int rx=0; int ry = 0; int rz=0; int ch3; GaussianBlur gb = new GaussianBlur(); double r2mean=0; int slicesDone = 0; int xCount=0; int xOffset2 = -15; int yOffset2 = -10; int zOffset2=0; int startSlice=1; if(Costes) { xOffset2=0; yOffset2=0; } if (Fay) iterations = 25; if (nslices>=2&&Fay) zOffset2=-1; if (Fay&&nslices>=2) {startSlice=2; nslices-=2; iterations=75; } if (vanS) {xOffset2=-21; startSlice=1; iterations=41; } int blockNumberX = (int)(width/(psf*2)); int blockNumberY = (int)(height/(psf*2)); ImageProcessor blockIndex = new ByteProcessor(blockNumberX,blockNumberY); double [] vSx = new double [41]; double [] vSr = new double [41] ; int ch4=0; int [] zUsed = new int [nslices]; //stackRand = new ImageStack(rwidth,rheight); int blockCount=0; boolean rBlock=true; int vacant; //start randomisations and calculation or Rrands for (int c=1; c<=iterations; c++) { stackRand = new ImageStack(rwidth,rheight); if(Fay) {if (c==26||c==51) {zOffset2 += 1; xOffset2=-15; yOffset2=-10; } if (xOffset2<10) xOffset2+=5; else {xOffset2 = -15; yOffset2+=5;} } if(vanS) { //IJ.showMessage("xOffset: "+xOffset2); xOffset2+=1; } for (int s=startSlice; s<=nslices; s++) { ipRand = new ShortProcessor(rwidth,rheight); slicesDone++; if (currentSlice) {s=currentSliceNo; nslices=s; } IJ.showStatus("Iteration "+c+ "/"+iterations+" Slice: "+s +"/" +nslices+" 'Esc' to abort"); if (IJ.escapePressed()) {IJ.beep(); return;} ip1= img1.getProcessor(s); ip2 = img2.getProcessor(s+zOffset2); ip3 = img3.getProcessor(s); for (int y=0; y1)&&ch2!=0) { ch3 = (int)((Math.random()*(ch2Max-ch2Min))+ch2Min) ; ipRand.putPixel(x,y,ch3); } } if (IJ.escapePressed()) {IJ.beep(); return;} //add to random image } } } if (Costes) gb.blur(ipRand, psf); stackRand.addSlice("Correlation Plot", ipRand); } //random image created now calculate r //reset values for r sumXX=0; sumX=0; sumXY = 0; sumYY = 0; sumY = 0; N = 0; N2=0; int s2=0; sumXtotal = 0; sumYtotal = 0; colocX = 0; colocY = 0; double ICQrand=0; int countPos2=0; countAll=0; //xOffset2=-21; if (IJ.escapePressed()) {IJ.beep(); return;} for (int s=startSlice; s<=nslices;s++) { s2=s; if (Fay&&nslices>=2) s2-=1; if (currentSlice) {s=currentSliceNo; nslices=s; s2=1; } ip1= img1.getProcessor(s); ip2 = stackRand.getProcessor(s2); for (int y=0; ych1) ch1Min = ch1; if (ch2Maxch2) ch2Min = ch2; N++; //Mander calc sumXtotal = sumXtotal+ch1; sumYtotal = sumYtotal+ch2; if(ch2>0) colocX = colocX + ch1; if(ch1>0) colocY = colocY + ch2; if((ch1+ch2!=0)) N2++; sumX = sumX+ch1; sumXY = sumXY + (ch1 * ch2); sumXX = sumXX + (ch1 * ch1); sumYY = sumYY + (ch2 *ch2); sumY = sumY + ch2; if (ch1+ch2!=0) {PDM = ((double)ch1-(double)ch1Mean)*((double)ch2-(double)ch2Mean); if (PDM>0) countPos2++; countAll++; } } } } } if (ignoreZeroZero) N = N2; ICQ2 = ((double)countPos2/(double)countAll)-0.5; ICQ2mean+=ICQ2; if (ICQobs>ICQ2) countICQ++; pearsons1 = sumXY - (sumX*sumY/N); pearsons2 = sumXX - (sumX*sumX/N); pearsons3 = sumYY - (sumY*sumY/N); r2= pearsons1/(Math.sqrt(pearsons2*pearsons3)); if(vanS) { vSx[c-1] = (double)xOffset2; vSr[c-1] = (double)r2 ; } if (r2r2max) r2max = r2; //IJ.write("Random "+ c + "\t"+df3.format(r2)+ "\t"+ df3.format(coloc2M1) + "\t"+df3.format(coloc2M2)); //IJ.write(df3.format(r2)); rVals.append(r2+"\n"); r2mean = r2mean+r2; if (r>r2) colocCount++; sumr2sqrd =sumr2sqrd +(r2*r2); sumr2 = sumr2+r2; sumICQ2sqrd +=(ICQ2*ICQ2); sumICQ2 +=ICQ2; //IJ.write(IJ.d2s(ICQ2,3)); } //done randomisations //calcualte mean Rrand r2mean = r2mean/iterations; r2sd = Math.sqrt(((iterations*(sumr2sqrd))-(sumr2*sumr2))/(iterations*(iterations-1))); ICQ2mean=ICQ2mean/iterations; double ICQ2sd =Math.sqrt(((iterations*(sumICQ2sqrd))-(sumICQ2*sumICQ2))/(iterations*(iterations-1))); double Zscore =(r-r2mean)/r2sd; double ZscoreICQ = (ICQobs-ICQ2mean)/ICQ2sd; String icqPercentile= "<50%"; String Percentile = ""+(iterations-colocCount)+"/"+iterations; //calculate percentage of Rrand that is less than Robs //code from: //http://www.cs.princeton.edu/introcs/26function/MyMath.java.html //Thanks to Bob Dougherty //50*{1 + erf[(V -mean)/(sqrt(2)*sdev)] double fx = 0.5*(1+erf(r-r2mean)/(Math.sqrt(2)*r2sd)); if (fx>=1) fx=1; if (fx<=0) fx=0; String Percentile2 = IJ.d2s(fx,3)+""; if(keep) new ImagePlus("Example random image", stackRand).show(); double percColoc = ((double)colocCount/(double)iterations)*100; double percICQ = ((double)countICQ/(double)iterations)*100; String Headings2 = "Image" +"\tR(obs)" +"\tR(rand) mean±sd" +"\tP-value" +"\tR(rand)>R(obs)" +"\tIterations" +" \tRandomisation" +"\tPSF width\n"; String strPSF = "na"; if (Costes||rBlocks) strPSF = df3.format(psf*pixelSize*2)+ " µm ("+df0.format(psf*2)+" pix.)" ; String str = fileName + "\t"+df3.format(r)+ "\t" +df3.format(r2mean) + "±"+ df3.format(r2sd)+ "\t"+Percentile2+ "\t" + (Percentile )+ "\t" +df0.format(iterations)+ "\t" + randMethod+ "\t" +strPSF; if (textWindow == null) textWindow = new TextWindow("Results", Headings2, str, 400, 250); else { textWindow.getTextPanel().setColumnHeadings(Headings2); textWindow.getTextPanel().appendLine(str); } IJ.selectWindow("Results"); if (showR) new TextWindow( "Random R values", "R(rand)", rVals.toString(),300, 400); if(vanS) {PlotWindow plot = new PlotWindow("CCF","x-translation","Pearsons",vSx,vSr); //r2min = (1.05*r2min); //r2max= (r2max*1.05); //plot.setLimits(-20, 20, r2min, r2max); plot.draw(); } Prefs.set("Rand_ignore.boolean", ignoreZeroZero); Prefs.set("Rand_keep.boolean", keep); Prefs.set("Rand_manPSF.int", manPSF); Prefs.set("Rand_smooth.boolean", smooth); if (Costes) Prefs.set("Rand_iterations.int", (int)iterations); Prefs.set("Rand_ch2L.double",ch2Lambda); Prefs.set("Rand_NA.double",NA); Prefs.set("Rand_pixelSize.double",pixelSize); Prefs.set("Rand_currentSlice.boolean", currentSlice); Prefs.set("Rand_useManPSF.boolean", useManPSF); Prefs.set("Rand_showR.boolean", showR); Prefs.set("Rand_indexRoi.int",indexRoi); Prefs.set("Rand_indexRand.int",indexRand); Prefs.set("Rand_randZ.boolean",randZ); long elapsedTime = (System.currentTimeMillis()-startTime)/1000; String units = "secs"; if (elapsedTime>90) {elapsedTime /= 60; units = "mins";} IJ.showStatus("Done. "+ elapsedTime+ " "+ units); } //code from: //http://www.cs.princeton.edu/introcs/26function/MyMath.java.html public static double erf(double z) { double t = 1.0 / (1.0 + 0.5 * Math.abs(z)); // use Horner's method double ans = 1 - t * Math.exp( -z*z - 1.26551223 + t * ( 1.00002368 + t * ( 0.37409196 + t * ( 0.09678418 + t * (-0.18628806 + t * ( 0.27886807 + t * (-1.13520398 + t * ( 1.48851587 + t * (-0.82215223 + t * ( 0.17087277)))))))))); if (z >= 0) return ans; else return -ans; } } diff --git a/src/main/java/Colocalisation_Threshold.java b/src/main/java/Colocalisation_Threshold.java index e17a7ac..94cd773 100644 --- a/src/main/java/Colocalisation_Threshold.java +++ b/src/main/java/Colocalisation_Threshold.java @@ -1,945 +1,946 @@ //22/4/5 import ij.IJ; import ij.ImagePlus; import ij.ImageStack; import ij.Prefs; import ij.WindowManager; import ij.gui.GenericDialog; import ij.gui.Roi; import ij.measure.Calibration; import ij.plugin.PlugIn; import ij.process.ColorProcessor; import ij.process.FloatProcessor; import ij.process.ImageConverter; import ij.process.ImageProcessor; import ij.process.ImageStatistics; import ij.process.ShortProcessor; import ij.text.TextWindow; import java.awt.Rectangle; import java.text.DecimalFormat; public class Colocalisation_Threshold implements PlugIn { boolean headingsSetCTC; private static int index1=0; private static int index2=1; private static int indexMask; private static boolean displayCounts; private static boolean useMask; private static boolean useRoi; private ImagePlus imp1, imp2; private static boolean threshold; private int ch1, ch2, ch3, nslices, width, height; private int indexRoi= (int)Prefs.get("CTC_indexRoi.int",0); private DecimalFormat df4 = new DecimalFormat("##0.0000"); private DecimalFormat df3 = new DecimalFormat("##0.000"); private DecimalFormat df2 = new DecimalFormat("##0.00"); private DecimalFormat df1 = new DecimalFormat("##0.0"); private DecimalFormat df0 = new DecimalFormat("##0"); private static boolean colocValConst= Prefs.get("CTC_colocConst.boolean", false); private int dualChannelIndex = (int)Prefs.get("CTC_channels.int",0); private static boolean bScatter= Prefs.get("CTC_bScatter.boolean", false); private static boolean bShowLocalisation=Prefs.get("CTC_show.boolean", false); boolean opt0 = Prefs.get("CTC_opt0.boolean", true); boolean opt1 = Prefs.get("CTC_opt1.boolean", true); boolean opt1a = Prefs.get("CTC_opt1a.boolean", true); boolean opt2 = Prefs.get("CTC_opt2.boolean", true); boolean opt3a = Prefs.get("CTC_opt3a.boolean", true); boolean opt3b = Prefs.get("CTC_opt3b.boolean", true); boolean opt4 = Prefs.get("CTC_opt4.boolean", true); boolean opt5 = Prefs.get("CTC_opt5.boolean", true); boolean opt6 = Prefs.get("CTC_opt6.boolean", true); boolean opt7 = Prefs.get("CTC_opt7.boolean", true); boolean opt8 = Prefs.get("CTC_opt8.boolean", true); boolean opt9 = Prefs.get("CTC_opt9.boolean", true); boolean opt10 = Prefs.get("CTC_opt10.boolean", true); String[] dualChannels= { "Red : Green","Red : Blue", "Green : Blue",}; private int colIndex1 = 0; private int colIndex2 = 1; private int colIndex3 = 2; ImageProcessor ip1, ip2, ipmask; ColorProcessor ipColoc; ImagePlus colocPix; private int rwidth, rheight, xOffset, yOffset; String[] chooseROI= { "None","Channel 1", "Channel 2",}; protected static TextWindow textWindow; + @Override public void run(String arg) { if (showDialog()) correlate(imp1, imp2); } public boolean showDialog() { int[] wList = WindowManager.getIDList(); if (wList==null) { IJ.noImage(); return false; } String[] titles = new String[wList.length]; String[] chooseMask= new String[wList.length+1]; chooseMask[0]=""; for (int i=0; i=titles.length)index1 = 0; if (index2>=titles.length)index2 = 0; if (indexMask>=titles.length)indexMask = 0; displayCounts = false; threshold = false; GenericDialog gd = new GenericDialog("Colocalisation Thresholds"); gd.addChoice("Channel_1", titles, titles[index1]); gd.addChoice("Channel_2", titles, titles[index2]); gd.addChoice("Use ROI", chooseROI, chooseROI[indexRoi]); // gd.addChoice("Mask channel", chooseMask, chooseMask[indexMask]); gd.addChoice("Channel Combination", dualChannels, dualChannels[dualChannelIndex]); gd.addCheckbox("Show Colocalized Pixel Map",bShowLocalisation); gd.addCheckbox("Use constant intensity for colocalized pixels",colocValConst); gd.addCheckbox("Show Scatter plot",bScatter); gd.addCheckbox("Include zero-zero pixels in threshold calculation",opt0); gd.addCheckbox("Set options",false); gd.addMessage("See: http://uhnresearch.ca/wcif/imagej"); gd.showDialog(); if (gd.wasCanceled()) return false; index1 = gd.getNextChoiceIndex(); index2 = gd.getNextChoiceIndex(); indexRoi = gd.getNextChoiceIndex(); // indexMask = gd.getNextChoiceIndex(); //IJ.showMessage(""+indexMask); dualChannelIndex = gd.getNextChoiceIndex(); bShowLocalisation = gd.getNextBoolean(); colocValConst = gd.getNextBoolean(); bScatter= gd.getNextBoolean(); opt0 =gd.getNextBoolean(); boolean options=gd.getNextBoolean(); imp1 = WindowManager.getImage(wList[index1]); imp2 = WindowManager.getImage(wList[index2]); useMask=false; //IJ.showMessage(""+indexMask); imp1 = WindowManager.getImage(wList[index1]); imp2 = WindowManager.getImage(wList[index2]); if (imp1.getType()!=ImagePlus.GRAY8&&imp1.getType()!=ImagePlus.GRAY16&&imp2.getType()!=ImagePlus.GRAY16 &&imp2.getType()!=ImagePlus.GRAY8) { IJ.showMessage("Image Correlator", "Both images must be 8-bit or 16-bit grayscale."); return false; } ip1 = imp1.getProcessor(); ip2 = imp2.getProcessor(); Roi roi1 = imp1.getRoi(); Roi roi2= imp2.getRoi(); width = imp1.getWidth(); height = imp1.getHeight(); useRoi=true; if (indexRoi== 0) useRoi = false; Rectangle rect =ip1.getRoi(); //IJ.showMessage("index"+rect.width); if ((indexRoi==1)) { if (roi1==null) { useRoi=false; } else { if (roi1.getType()==Roi.RECTANGLE) { IJ.showMessage("Does not work with rectangular ROIs"); return false; } ipmask = imp1.getMask(); //if (keepROIimage) new ImagePlus("Mask",ipmask).show(); rect = ip1.getRoi(); } } if ((indexRoi==2)) { if (roi2==null) { useRoi=false; } else { if (roi2.getType()==Roi.RECTANGLE) { IJ.showMessage("Does not work with rectangular ROIs"); return false; } ipmask = imp2.getMask(); //if (keepROIimage) new ImagePlus("Mask",ipmask).show(); rect = ip2.getRoi(); } } if (indexRoi==0) { xOffset = 0; yOffset = 0; rwidth=width; rheight =height; } else { xOffset = rect.x; yOffset = rect.y; rwidth=rect.width; rheight =rect.height; } //if red-blue if (dualChannelIndex==1) { colIndex2 = 2; colIndex3=1; }; //if blue-green if (dualChannelIndex==2) { colIndex1 = 1; colIndex2 =2 ; colIndex3=0; } if (options) { GenericDialog gd2 = new GenericDialog("Set Results Options"); gd2.addMessage("See online manual for detailed description of these values"); gd2.addCheckbox("Show linear regression solution",opt1a); gd2.addCheckbox("Show thresholds",opt1); gd2.addCheckbox("Pearson's for whole image",opt2); gd2.addCheckbox("Pearson's for image above thresholds",opt3a); gd2.addCheckbox("Pearson's for image below thresholds (should be ~0)",opt3b); gd2.addCheckbox("Mander's original coefficients (threshold = 0)",opt4); gd2.addCheckbox("Mander's using thresholds",opt5); gd2.addCheckbox("Number of colocalized voxels",opt6); gd2.addCheckbox("% Image volume colocalized",opt7); gd2.addCheckbox("% Voxels colocalized",opt8); gd2.addCheckbox("% Intensity colocalized",opt9); gd2.addCheckbox("% Intensity above threshold colocalized",opt10); gd2.showDialog(); if (gd2.wasCanceled()) return false; opt1=gd2.getNextBoolean(); opt1a=gd2.getNextBoolean(); opt2=gd2.getNextBoolean(); opt3a=gd2.getNextBoolean(); opt3b=gd2.getNextBoolean(); opt4=gd2.getNextBoolean(); opt5=gd2.getNextBoolean(); opt6=gd2.getNextBoolean(); opt7=gd2.getNextBoolean(); opt8=gd2.getNextBoolean(); opt9=gd2.getNextBoolean(); opt10=gd2.getNextBoolean(); headingsSetCTC = false; } //IJ.showMessage(""+indexMask); return true; } public void correlate(ImagePlus imp1, ImagePlus imp2) {//IJ.showMessage("mask? "+useMask); String ch1fileName = imp1.getTitle(); String ch2fileName = imp2.getTitle(); //String maskName = impMask.getTitle(); String fileName = ch1fileName + " & " + ch2fileName; ImageProcessor ip1 = imp1.getProcessor(); ImageProcessor ip2 = imp2.getProcessor(); Calibration spatialCalibration = imp1.getCalibration(); ImageProcessor ipMask = imp1.getMask(); if (indexRoi>1) ipMask = imp2.getMask(); // ImageStack imgMask = impMask.getStack(); ImageStack img1 = imp1.getStack(); ImageStack img2 = imp2.getStack(); if (indexRoi== 0) useRoi = false; Rectangle rect1 =ip1.getRoi(); Rectangle rect2 =ip2.getRoi(); Roi roi1 = imp1.getRoi(); Roi roi2= imp2.getRoi(); nslices = imp1.getStackSize(); width = imp1.getWidth(); height = imp1.getHeight(); ipColoc = new ColorProcessor(rwidth,rheight); ImageStack stackColoc = new ImageStack(rwidth,rheight); MinMaxContainer minMax1 = getMinMax(ip1); MinMaxContainer minMax2 = getMinMax(ip2); double ch1threshmin=0; double ch1threshmax=minMax1.max; double ch2threshmin=0; double ch2threshmax=minMax2.max; double pearsons1 = 0; double pearsons2 = 0; double pearsons3 = 0; double r = 1; pearsons1 =0; pearsons2 = 0; pearsons3 = 0; double r2= 1; boolean thresholdFound=false; boolean unsigned = true; int count =0; double sumX = 0; double sumXY = 0; double sumXX = 0; double sumYY = 0; double sumY = 0; double colocX = 0; double colocY = 0; double countX = 0; double countY = 0; double sumXYm=0; int Nch1=0,Nch2=0; double oldMax=0; int sumCh2gtT =0; int sumCh1gtT =0; int N = 0; int N2 = 0; int Nzero=0; int Nch1gtT=0; int Nch2gtT=0; double oldMax2=0; int ch1Max=(int)minMax1.max; int ch2Max=(int)minMax2.max; int ch1Min = (int)minMax1.min; int ch2Min = (int)minMax2.min; int ch1ROIMax=0; int ch2ROIMax=0; String Headings = "\t \t \t \t \t \t \t \n"; ImageProcessor plot32 = new FloatProcessor(256, 256); ImageProcessor plot16 = new ShortProcessor(256, 256); int scaledXvalue =0; int scaledYvalue=0; if (ch1Max<255) ch1Max=255; if (ch2Max<255) ch2Max=255; double ch1Scaling = (double)255/(double)ch1Max; double ch2Scaling = (double)255/(double)ch2Max; // scaling for both channels should be the same so the scatterplot is not squewed double chScaling = 1; if (ch1Scaling>ch2Scaling) chScaling = ch1Scaling; else chScaling = ch2Scaling; boolean divByZero=false; StringBuffer sb= new StringBuffer(); String str=""; int i = imp1.getCurrentSlice(); double bBest=0; double mBest = 0; double bestr2=1; double ch1BestThresh=0; double ch2BestThresh=0; String mString; //start regression IJ.showStatus("1/4: Performing regression. Press 'Esc' to abort"); int ch1Sum=0; int ch2Sum=0; int ch3Sum=0; double ch1mch1MeanSqSum=0; double ch2mch2MeanSqSum= 0; double ch3mch3MeanSqSum= 0; ImageProcessor ipPlot = new ShortProcessor (256,256); int mask=0; if (indexRoi== 0) useRoi = false; Rectangle rect =ip1.getRoi(); if ((indexRoi==1)) { if (roi1==null) { useRoi=false; } else { if (roi1.getType()==Roi.RECTANGLE) { IJ.showMessage("Does not work with rectangular ROIs"); return; } ipMask = imp1.getMask(); //if (keepROIimage) new ImagePlus("Mask",ipmask).show(); rect = ip1.getRoi(); } } if ((indexRoi==2)) { if (roi2==null) { useRoi=false; } else { if (roi2.getType()==Roi.RECTANGLE) { IJ.showMessage("Does not work with rectangular ROIs"); return ; } ipMask = imp2.getMask(); //if (keepROIimage) new ImagePlus("Mask",ipmask).show(); rect = ip2.getRoi(); } } if (useRoi==false) { xOffset = 0; yOffset = 0; rwidth=width; rheight =height; } else { xOffset = rect.x; yOffset = rect.y; rwidth=rect.width; rheight =rect.height; } //new ImagePlus("Mask",ipMask).show(); //get means for (int s=1; s<=nslices; s++) { if (IJ.escapePressed()) { IJ.beep(); return; } ip1 = img1.getProcessor(s); ip2 = img2.getProcessor(s); //ipMask = imgMask.getProcessor(s); for (int y=0; ych1ROIMax) ch1ROIMax=ch1; if (ch2>ch2ROIMax) ch2ROIMax=ch2; ch3 = ch1+ch2; ch1Sum+=ch1; ch2Sum+=ch2; ch3Sum+=ch3; if (ch1+ch2!=0) N++; } } } } double ch1Mean = ch1Sum/N; double ch2Mean = ch2Sum/N; double ch3Mean = ch3Sum/N; N=0; // Do some rather simple performance testing long startTime = System.currentTimeMillis(); //calulate variances for (int s=1; s<=nslices; s++) { if (IJ.escapePressed()) { IJ.beep(); return; } ip1 = img1.getProcessor(s); ip2 = img2.getProcessor(s); //ipMask = imgMask.getProcessor(s); for (int y=0; yr2*r2)) { ch1BestThresh=ch1threshmax; bestr2=r2; } //if our r is close to our level of tolerance then set threshold has been found if ((r2-tolerance) )thresholdFound = true; //if we've reached ch1 =1 then we've exhausted posibilities if (Math.round(ch1threshmax)==0) thresholdFound =true; oldMax= newMax; //change threshold max if (r2>=0) { if ((r2>=r2Prev)&&(!divByZero)) newMax = newMax/2; if ((r20) { Nch1++; mCh2coloc = mCh2coloc+ch2; } if (ch2>0) { Nch2++; mCh1coloc = mCh1coloc+ch1; } if ((double)ch2>=ch2threshmax) { Nch2gtT++; sumCh2gtT = sumCh2gtT+ch2; colocX=colocX+ch1; } if ((double)ch1>=ch1threshmax) { Nch1gtT++; sumCh1gtT = sumCh1gtT+ch1; colocY=colocY+ch2; } if (((double)ch1>ch1threshmax)&&((double)ch2>ch2threshmax)) { sumColocCh1 = sumColocCh1+ch1; sumColocCh2 = sumColocCh2+ch2; Ncoloc++; //calc pearsons sumX = sumX+ch1; sumXY = sumXY + (ch1 * ch2); sumXX = sumXX + (ch1 * ch1); sumYY = sumYY + (ch2 *ch2); sumY = sumY + ch2; } } } } } //IJ.showMessage("Totoal"+N+" N0:"+Nzero+" Nc :"+ Ncoloc); pearsons1 = sumXY - (sumX*sumY/Ncoloc); pearsons2 = sumXX - (sumX*sumX/Ncoloc); pearsons3 = sumYY - (sumY*sumY/Ncoloc); //Pearsons for coloclaised volume double Rcoloc= pearsons1/(Math.sqrt(pearsons2*pearsons3)); //Mander's original //[i.e. E(ch1if ch2>0) ÷ E(ch1total)] double M1 = mCh1coloc /sumCh1total; double M2 = mCh2coloc /sumCh2total; //Manders using threshold //[i.e. E(ch1 if ch2>ch2threshold) ÷ (Ech1total)] double colocM1 = (double) colocX/(double)sumCh1total; double colocM2 = (double) colocY/(double)sumCh2total; //as in Coste's paper //[i.e. E(ch1>ch1threshold) ÷ E(ch1total)] double colocC1 = (double)sumCh1gtT/ (double)sumCh1total; double colocC2 = (double)sumCh2gtT/(double)sumCh2total; //Imaris percentage volume double percVolCh1 = (double)Ncoloc/ (double)Nch1gtT; double percVolCh2 = (double)Ncoloc/(double)Nch2gtT; double percTotCh1 = (double) sumColocCh1/ (double)sumCh1total; double percTotCh2 = (double) sumColocCh2/ (double)sumCh2total; //Imaris percentage material double percMatCh1 = (double) sumColocCh1/(double)sumCh1gtT; double percMatCh2 = (double)sumColocCh2/(double)sumCh2gtT; //IJ.showMessage("Totoal"+N+" N0:"+Nzero+" Nc :"+ Ncoloc); sb.append(fileName+"\n"); //if (!useMask) maskName = ""; str = fileName +"\t"+"ROI" + indexRoi+"\t"; str += opt0 ? "incl.\t" : "excl.\t"; if (opt2) str+= df3.format(rTotal)+ "\t"; if (opt1a) str+= df3.format(m)+ "\t "+df1.format(b)+ "\t"; if (opt1) str+= IJ.d2s(ch1threshmax,0)+"\t"+IJ.d2s(ch2threshmax,0)+"\t"; if (opt3a) str+= df4.format(Rcoloc) +"\t"; if (opt3b) str+= df3.format(bestr2)+"\t"; if (opt4) str+= df4.format(M1)+ "\t "+df4.format(M2)+"\t"; if (opt5) str+= df4.format(colocM1)+ "\t"+df4.format(colocM2)+"\t"; if (opt6) str+= Ncoloc+ "\t"; if (opt7) str+= df2.format(((double)Ncoloc*(double)100)/((double)width*(double)height*(double)nslices))+"%\t"; if (opt8) str+= df2.format(percVolCh1*100 )+ "%\t"; if (opt8) str+= df2.format(percVolCh2*100 )+ "%\t"; if (opt9) str+= df2.format(percTotCh1*100 )+ "%\t"; if (opt9) str+= df2.format(percTotCh2*100 )+ "%\t"; if (opt10) str+= df2.format(percMatCh1*100 )+ "%\t"; if (opt10) str+= df2.format(percMatCh2*100 )+ "%\t"; String heading = "Images\tMask\tZeroZero\t"; if (opt2) heading += "Rtotal\t"; if (opt1a) heading += "m\tb\t"; if (opt1) heading += "Ch1 thresh\tCh2 thresh\t"; if (opt3a) heading += "Rcoloc\t"; if (opt3b) heading += "RcolocPixelsImageThresh1)&&((double)ch2>colocPixelsImageThresh2)) { colocInt=255; if (!colocValConst) { colocInt = (int)Math.sqrt(ch1*ch2); } color[colIndex1 ]=(int)colocInt; color[colIndex2 ]=(int)colocInt; color[colIndex3 ]=(int)colocInt; ipColoc.putPixel(x,y,color ); } } } } //IJ.showMessage(stackColoc.getWidth()+ " - " + ipColoc.getWidth()); stackColocPix.addSlice("Colocalized Pixel Map Image", ipColoc); } return stackColocPix; } private class MinMaxContainer { public double min; public double max; public MinMaxContainer(double min, double max){ this.min = min; this.max = max; } } } diff --git a/src/main/java/algorithms/AutoThresholdRegression.java b/src/main/java/algorithms/AutoThresholdRegression.java index 2b14874..8aac127 100644 --- a/src/main/java/algorithms/AutoThresholdRegression.java +++ b/src/main/java/algorithms/AutoThresholdRegression.java @@ -1,326 +1,327 @@ package algorithms; import gadgets.DataContainer; import gadgets.ThresholdMode; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; import results.ResultHandler; /** * A class implementing the automatic finding of a threshold * used for Person colocalisation calculation. */ public class AutoThresholdRegression> extends Algorithm { // Identifiers for choosing which implementation to use public enum Implementation {Costes, Bisection}; Implementation implementation = Implementation.Bisection; /* The threshold for ration of y-intercept : y-mean to raise a warning about * it being to high or low, meaning far from zero. Don't use y-max as before, * since this could be a very high value outlier. Mean is probably more * reliable. */ final double warnYInterceptToYMeanRatioThreshold = 0.01; // the slope and and intercept of the regression line double autoThresholdSlope = 0.0, autoThresholdIntercept = 0.0; /* The thresholds for both image channels. Pixels below a lower * threshold do NOT include the threshold and pixels above an upper * one will NOT either. Pixels "in between (and including)" thresholds * do include the threshold values. */ T ch1MinThreshold, ch1MaxThreshold, ch2MinThreshold, ch2MaxThreshold; // additional information double bToYMeanRatio = 0.0; //This is the Pearson's correlation we will use for further calculations PearsonsCorrelation pearsonsCorrellation; public AutoThresholdRegression(PearsonsCorrelation pc) { this(pc, Implementation.Costes); } public AutoThresholdRegression(PearsonsCorrelation pc, Implementation impl) { super("auto threshold regression"); pearsonsCorrellation = pc; implementation = impl; } @Override public void execute(DataContainer container) throws MissingPreconditionException { // get the 2 images for the calculation of Pearson's final RandomAccessibleInterval img1 = container.getSourceImage1(); final RandomAccessibleInterval img2 = container.getSourceImage2(); final RandomAccessibleInterval mask = container.getMask(); double ch1Mean = container.getMeanCh1(); double ch2Mean = container.getMeanCh2(); double combinedMean = ch1Mean + ch2Mean; // get the cursors for iterating through pixels in images TwinCursor cursor = new TwinCursor( img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); // variables for summing up the double ch1MeanDiffSum = 0.0, ch2MeanDiffSum = 0.0, combinedMeanDiffSum = 0.0; double combinedSum = 0.0; int N = 0, NZero = 0; // reference image data type final T type = cursor.getFirst(); while (cursor.hasNext()) { cursor.fwd(); double ch1 = cursor.getFirst().getRealDouble(); double ch2 = cursor.getSecond().getRealDouble(); combinedSum = ch1 + ch2; // TODO: Shouldn't the whole calculation take only pixels // into account that are combined above zero? And not just // the denominator (like it is done now)? // calculate the numerators for the variances ch1MeanDiffSum += (ch1 - ch1Mean) * (ch1 - ch1Mean); ch2MeanDiffSum += (ch2 - ch2Mean) * (ch2 - ch2Mean); combinedMeanDiffSum += (combinedSum - combinedMean) * (combinedSum - combinedMean); // count only pixels that are above zero if ( (ch1 + ch2) > 0.00001) NZero++; N++; } double ch1Variance = ch1MeanDiffSum / (N - 1); double ch2Variance = ch2MeanDiffSum / (N - 1); double combinedVariance = combinedMeanDiffSum / (N - 1.0); //http://mathworld.wolfram.com/Covariance.html //?2 = X2?(X)2 // = E[X2]?(E[X])2 //var (x+y) = var(x)+var(y)+2(covar(x,y)); //2(covar(x,y)) = var(x+y) - var(x)-var(y); double ch1ch2Covariance = 0.5*(combinedVariance - (ch1Variance + ch2Variance)); // calculate regression parameters double denom = 2*ch1ch2Covariance; double num = ch2Variance - ch1Variance + Math.sqrt( (ch2Variance - ch1Variance) * (ch2Variance - ch1Variance) + (4 * ch1ch2Covariance *ch1ch2Covariance) ); final double m = num/denom; final double b = ch2Mean - m*ch1Mean ; // A stepper that walks thresholds Stepper stepper; // to map working thresholds to channels ChannelMapper mapper; // let working threshold walk on channel one if the regression line // leans more towards the abscissa (which represents channel one) for // positive and negative correlation. if (m > -1 && m < 1.0) { // Map working threshold to channel one (because channel one has a // larger maximum value. mapper = new ChannelMapper() { @Override public double getCh1Threshold(double t) { return t; } @Override public double getCh2Threshold(double t) { return (t * m) + b; } }; // Select a stepper if (implementation == Implementation.Bisection) { // Start at the midpoint of channel one stepper = new BisectionStepper( Math.abs(container.getMaxCh1() + container.getMinCh1()) * 0.5, container.getMaxCh1()); } else { stepper = new SimpleStepper(container.getMaxCh1()); } } else { // Map working threshold to channel two (because channel two has a // larger maximum value. mapper = new ChannelMapper() { @Override public double getCh1Threshold(double t) { return (t - b) / m; } @Override public double getCh2Threshold(double t) { return t; } }; // Select a stepper if (implementation == Implementation.Bisection) { // Start at the midpoint of channel two stepper = new BisectionStepper( Math.abs(container.getMaxCh2() + container.getMinCh2()) * 0.5, container.getMaxCh2()); } else { stepper = new SimpleStepper(container.getMaxCh2()); } } // Min threshold not yet implemented double ch1ThreshMax = container.getMaxCh1(); double ch2ThreshMax = container.getMaxCh2(); // define some image type specific threshold variables T thresholdCh1 = type.createVariable(); T thresholdCh2 = type.createVariable(); // reset the previously created cursor cursor.reset(); /* Get min and max value of image data type. Since type of image * one and two are the same, we dont't need to distinguish them. */ T dummyT = type.createVariable(); final double minVal = dummyT.getMinValue(); final double maxVal = dummyT.getMaxValue(); // do regression while (!stepper.isFinished()) { // round ch1 threshold and compute ch2 threshold ch1ThreshMax = Math.round(mapper.getCh1Threshold(stepper.getValue())); ch2ThreshMax = Math.round(mapper.getCh2Threshold(stepper.getValue())); /* Make sure we don't get overflow the image type specific threshold variables * if the image data type doesn't support this value. */ thresholdCh1.setReal(clamp(ch1ThreshMax, minVal, maxVal)); thresholdCh2.setReal(clamp(ch2ThreshMax, minVal, maxVal)); try { // do persons calculation within the limits final double currentPersonsR = pearsonsCorrellation.calculatePearsons(cursor, ch1Mean, ch2Mean, thresholdCh1, thresholdCh2, ThresholdMode.Below); stepper.update(currentPersonsR); } catch (MissingPreconditionException e) { /* the exception that could occur is due to numerical * problems within the Pearsons calculation. */ stepper.update(Double.NaN); } // reset the cursor to reuse it cursor.reset(); } /* Store the new results. The lower thresholds are the types * min value for now. For the max threshold we do a clipping * to make it fit into the image type. */ ch1MinThreshold = type.createVariable(); ch1MinThreshold.setReal(minVal); ch1MaxThreshold = type.createVariable(); ch1MaxThreshold.setReal(clamp(ch1ThreshMax, minVal, maxVal)); ch2MinThreshold = type.createVariable(); ch2MinThreshold.setReal(minVal); ch2MaxThreshold = type.createVariable(); ch2MaxThreshold.setReal(clamp(ch2ThreshMax, minVal, maxVal)); autoThresholdSlope = m; autoThresholdIntercept = b; bToYMeanRatio = b / container.getMeanCh2(); // add warnings if values are not in tolerance range if ( Math.abs(bToYMeanRatio) > warnYInterceptToYMeanRatioThreshold ) { addWarning("y-intercept far from zero", "The ratio of the y-intercept of the auto threshold regression " + "line to the mean value of Channel 2 is high. This means the " + "y-intercept is far from zero, implying a significant positive " + "or negative zero offset in the image data intensities. Maybe " + "you should use a ROI. Maybe do a background subtraction in " + "both channels. Make sure you didn't clip off the low " + "intensities to zero. This might not affect Pearson's " + "correlation values very much, but might harm other results."); } // add warning if threshold is above the image mean if (ch1ThreshMax > ch1Mean) { addWarning("Threshold of ch. 1 too high", "Too few pixels are taken into account for above-threshold calculations. The threshold is above the channel's mean."); } if (ch2ThreshMax > ch2Mean) { addWarning("Threshold of ch. 2 too high", "Too few pixels are taken into account for above-threshold calculations. The threshold is above the channel's mean."); } // add warnings if values are below lowest pixel value of images if ( ch1ThreshMax < container.getMinCh1() || ch2ThreshMax < container.getMinCh2() ) { String msg = "The auto threshold method could not find a positive " + "threshold, so thresholded results are meaningless."; msg += implementation == Implementation.Costes ? "" : " Maybe you should try classic thresholding."; addWarning("thresholds too low", msg); } } /** * Clamp a value to a min or max value. If the value is below min, min is * returned. Accordingly, max is returned if the value is larger. If it is * neither, the value itself is returned. */ public static double clamp(double val, double min, double max) { return min > val ? min : max < val ? max : val; } + @Override public void processResults(ResultHandler handler) { super.processResults(handler); handler.handleValue( "m (slope)", autoThresholdSlope , 2 ); handler.handleValue( "b (y-intercept)", autoThresholdIntercept, 2 ); handler.handleValue( "b to y-mean ratio", bToYMeanRatio, 2 ); handler.handleValue( "Ch1 Max Threshold", ch1MaxThreshold.getRealDouble(), 2); handler.handleValue( "Ch2 Max Threshold", ch2MaxThreshold.getRealDouble(), 2); handler.handleValue( "Threshold regression", implementation.toString()); } public double getBToYMeanRatio() { return bToYMeanRatio; } public double getWarnYInterceptToYMaxRatioThreshold() { return warnYInterceptToYMeanRatioThreshold; } public double getAutoThresholdSlope() { return autoThresholdSlope; } public double getAutoThresholdIntercept() { return autoThresholdIntercept; } public T getCh1MinThreshold() { return ch1MinThreshold; } public T getCh1MaxThreshold() { return ch1MaxThreshold; } public T getCh2MinThreshold() { return ch2MinThreshold; } public T getCh2MaxThreshold() { return ch2MaxThreshold; } } diff --git a/src/main/java/algorithms/CostesSignificanceTest.java b/src/main/java/algorithms/CostesSignificanceTest.java index 2c2f307..888394a 100644 --- a/src/main/java/algorithms/CostesSignificanceTest.java +++ b/src/main/java/algorithms/CostesSignificanceTest.java @@ -1,400 +1,401 @@ package algorithms; import gadgets.DataContainer; import gadgets.DataContainer.MaskType; import gadgets.Statistics; import java.util.ArrayList; import java.util.Arrays; import java.util.Collections; import java.util.List; import net.imglib2.Cursor; import net.imglib2.IterableInterval; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; import net.imglib2.algorithm.gauss.Gauss; import net.imglib2.img.Img; import net.imglib2.img.ImgFactory; import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.roi.RectangleRegionOfInterest; import net.imglib2.type.NativeType; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; import results.ResultHandler; public class CostesSignificanceTest & NativeType> extends Algorithm { // radius of the PSF in pixels, its size *must* for now be three protected double[] psfRadius = new double[3]; // indicates if the shuffled images should be shown as a result boolean showShuffledImages = false; // the number of randomization tests int nrRandomizations; // the shuffled image last worked on Img smoothedShuffledImage; // the Pearson's algorithm (that should have been run before) PearsonsCorrelation pearsonsCorrelation; // a list of resulting Pearsons values from the randomized images List shuffledPearsonsResults; /* the amount of Pearson's values with shuffled data * that has the value of the original one or is larger. */ int shuffledPearsonsNotLessOriginal = 0; // The mean of the shuffled Pearson values double shuffledMean = 0.0; // The standard derivation of the shuffled Pearson values double shuffledStdDerivation = 0.0; /* The Costes P-Value which is the probability that * Pearsons r is different from the mean of the randomized * r values. */ double costesPValue; // the maximum retries in case of Pearson numerical errors protected final int maxErrorRetries = 3; /** * Creates a new Costes significance test object by using a * cube block with the given edge length. * * @param psfRadiusInPixels The edge width of the 3D cube block. */ public CostesSignificanceTest(PearsonsCorrelation pc, int psfRadiusInPixels, int nrRandomizations, boolean showShuffledImages) { super("Costes significance test"); this.pearsonsCorrelation = pc; Arrays.fill(psfRadius, psfRadiusInPixels); this.nrRandomizations = nrRandomizations; this.showShuffledImages = showShuffledImages; } /** * Builds a list of blocks that represent the images. To * do so we create a list image ROI cursors. If a block * does not fit into the image it will get a out-of-bounds * strategy. */ @Override public void execute(DataContainer container) throws MissingPreconditionException { final RandomAccessibleInterval img1 = container.getSourceImage1(); final RandomAccessibleInterval img2 = container.getSourceImage2(); final RandomAccessibleInterval mask = container.getMask(); /* To determine the number of needed blocks, we need * the effective dimensions of the image. Since the * mask is responsible for this, we ask for its size. */ long[] dimensions = container.getMaskBBSize(); int nrDimensions = dimensions.length; // calculate the needed number of blocks per image int nrBlocksPerImage = 1; long[] nrBlocksPerDimension = new long[3]; for (int i = 0; i < nrDimensions; i++) { // add the amount of full fitting blocks to the counter nrBlocksPerDimension[i] = (long) (dimensions[i] / psfRadius[i]); // if there is the need for a out-of-bounds block, increase count if ( dimensions[i] % psfRadius[i] != 0 ) nrBlocksPerDimension[i]++; // increase total count nrBlocksPerImage *= nrBlocksPerDimension[i]; } /* For creating the input and output blocks we need * offset and size as floating point array. */ double[] floatOffset = new double[ img1.numDimensions() ]; long[] longOffset = container.getMaskBBOffset(); for (int i=0; i< longOffset.length; ++i ) floatOffset[i] = longOffset[i]; double[] floatDimensions = new double[ nrDimensions ]; for (int i=0; i< nrDimensions; ++i ) floatDimensions[i] = dimensions[i]; /* Create the ROI blocks. The image dimensions might not be * divided cleanly by the block size. Therefore we need to * have an out of bounds strategy -- a mirror. */ List> blockIntervals; blockIntervals = new ArrayList>( nrBlocksPerImage ); RandomAccessible< T> infiniteImg = Views.extendMirrorSingle( img1 ); generateBlocks( infiniteImg, blockIntervals, floatOffset, floatDimensions); // create input and output cursors and store them along their offset List> inputBlocks = new ArrayList>(nrBlocksPerImage); List> outputBlocks = new ArrayList>(nrBlocksPerImage); for (IterableInterval roiIt : blockIntervals) { inputBlocks.add(roiIt.localizingCursor()); outputBlocks.add(roiIt.localizingCursor()); } // we will need a zero variable final T zero = img1.randomAccess().get().createVariable(); zero.setZero(); /* Create a new image to contain the shuffled data and with * same dimensions as the original data. */ final long[] dims = new long[img1.numDimensions()]; img1.dimensions(dims); ImgFactory factory = new ArrayImgFactory(); Img shuffledImage = factory.create( dims, img1.randomAccess().get().createVariable() ); RandomAccessible< T> infiniteShuffledImage = Views.extendValue(shuffledImage, zero ); // create a double version of the PSF for the smoothing double[] smoothingPsfRadius = new double[nrDimensions]; for (int i = 0; i < nrDimensions; i++) { smoothingPsfRadius[i] = (double) psfRadius[i]; } // the retry count for error cases int retries = 0; shuffledPearsonsResults = new ArrayList(); for (int i=0; i < nrRandomizations; i++) { // shuffle the list Collections.shuffle( inputBlocks ); // get an output random access RandomAccess output = infiniteShuffledImage.randomAccess(); // check if a mask is in use and further actions are needed if (container.getMaskType() == MaskType.Irregular) { Cursor siCursor = shuffledImage.cursor(); // black the whole intermediate image, just in case we have irr. masks while (siCursor.hasNext()) { siCursor.fwd(); output.setPosition(siCursor); output.get().setZero(); } } // write out the shuffled input blocks into the output blocks for (int j=0; j inputCursor = inputBlocks.get(j); Cursor outputCursor = outputBlocks.get(j); /* Iterate over both blocks. Theoretically the iteration * order could be different. Because we are dealing with * randomized data anyway, this is not a problem here. */ while (inputCursor.hasNext() && outputCursor.hasNext()) { inputCursor.fwd(); outputCursor.fwd(); output.setPosition(outputCursor); // write the data output.get().set( inputCursor.get() ); } /* Reset both cursors. If we wouldn't do that, the * image contents would not change on the next pass. */ inputCursor.reset(); outputCursor.reset(); } smoothedShuffledImage = Gauss.inFloat( smoothingPsfRadius, shuffledImage); try { // calculate correlation value... double pValue = pearsonsCorrelation.calculatePearsons( smoothedShuffledImage, img2, mask); // ...and add it to the results list shuffledPearsonsResults.add( pValue ); } catch (MissingPreconditionException e) { /* if the randomized input data does not suit due to numerical * problems, try it three times again and then fail. */ if (retries < maxErrorRetries) { // increase retry count and the number of randomizations retries++; nrRandomizations++; } else { throw new MissingPreconditionException("Maximum retries have been made (" + + retries + "), but errors keep on coming: " + e.getMessage(), e); } } } // calculate statistics on the randomized values and the original one double originalVal = pearsonsCorrelation.getPearsonsCorrelationValue(); calculateStatistics(shuffledPearsonsResults, originalVal); } /** * This method drives the creation of RegionOfInterest-Cursors on the given image. * It does not matter if those generated blocks are used for reading and/or * writing. The resulting blocks are put into the given list and are in the * responsibility of the caller, i.e. he or she must make sure the cursors get * closed on some point in time. * * @param img The image to create cursors on. * @param blockList The list to put newly created cursors into * @param outOfBoundsFactory Defines what to do if a block has parts out of image bounds. */ protected void generateBlocks(RandomAccessible img, List> blockList, double[] offset, double[] size) throws MissingPreconditionException { // get the number of dimensions int nrDimensions = img.numDimensions(); if (nrDimensions == 2) { // for a 2D image... generateBlocksXY(img, blockList, offset, size); } else if (nrDimensions == 3) { // for a 3D image... final double depth = size[2]; double z; double originalZ = offset[2]; // go through the depth in steps of block depth for ( z = psfRadius[2]; z <= depth; z += psfRadius[2] ) { offset[2] = originalZ + z - psfRadius[2]; generateBlocksXY(img, blockList, offset, size); } // check is we need to add a out of bounds strategy cursor if (z > depth) { offset[2] = originalZ + z - psfRadius[2]; generateBlocksXY(img, blockList, offset, size); } offset[2] = originalZ; } else throw new MissingPreconditionException("Currently only 2D and 3D images are supported."); } /** * Goes stepwise through the y-dimensions of the image data and adds cursors * for each row to the given list. The method does not check if there is a * y-dimensions, so this should be made sure before. you can enforce to * create all cursors as out-of-bounds one. * * @param img The image to get the data and cursors from. * @param blockList The list to put the blocks into. * @param offset The current offset configuration. Only [0] and [1] will be changed. * @param outOfBoundsFactory The factory to create out-of-bounds-cursors with. * @param forceOutOfBounds Indicates if all cursors created should be out-of-bounds ones. */ protected void generateBlocksXY(RandomAccessible img, List> blockList, double[] offset, double[] size) { // potentially masked image height double height = size[1]; final double originalY = offset[1]; // go through the height in steps of block width double y; for ( y = psfRadius[1]; y <= height; y += psfRadius[1] ) { offset[1] = originalY + y - psfRadius[1]; generateBlocksX(img, blockList, offset, size); } // check is we need to add a out of bounds strategy cursor if (y > height) { offset[1] = originalY + y - psfRadius[1]; generateBlocksX(img, blockList, offset, size); } offset[1] = originalY; } /** * Goes stepwise through a row of image data and adds cursors to the given list. * If there is not enough image data for a whole block, an out-of-bounds cursor * is generated. The creation of out-of-bound cursors could be enforced as well. * * @param img The image to get the data and cursors from. * @param blockList The list to put the blocks into. * @param offset The current offset configuration. Only [0] of it will be changed. * @param outOfBoundsFactory The factory to create out-of-bounds-cursors with. * @param forceOutOfBounds Indicates if all cursors created should be out-of-bounds ones. */ protected void generateBlocksX(RandomAccessible img, List> blockList, double[] offset, double[] size) { // potentially masked image width double width = size[0]; final double originalX = offset[0]; // go through the width in steps of block width double x; for ( x = psfRadius[0]; x <= width; x += psfRadius[0] ) { offset[0] = originalX + x - psfRadius[0]; RectangleRegionOfInterest roi = new RectangleRegionOfInterest(offset.clone(), psfRadius.clone()); IterableInterval roiInterval = roi.getIterableIntervalOverROI(img); blockList.add(roiInterval); } // check is we need to add a out of bounds strategy cursor if (x > width) { offset[0] = originalX + x - psfRadius[0]; RectangleRegionOfInterest roi = new RectangleRegionOfInterest(offset.clone(), psfRadius.clone()); IterableInterval roiInterval = roi.getIterableIntervalOverROI(img); blockList.add(roiInterval); } offset[0] = originalX; } protected void calculateStatistics(List compareValues, double originalVal) { shuffledPearsonsNotLessOriginal = 0; int iterations = shuffledPearsonsResults.size(); double compareSum = 0.0; for( Double shuffledVal : shuffledPearsonsResults ) { double diff = shuffledVal - originalVal; /* check if the randomized Pearsons value is equal * or larger than the original one. */ if( diff > -0.00001 ) { shuffledPearsonsNotLessOriginal++; } compareSum += shuffledVal; } shuffledMean = compareSum / iterations; shuffledStdDerivation = Statistics.stdDeviation(compareValues); // get the quantile of the original value in the shuffle distribution costesPValue = Statistics.phi(originalVal, shuffledMean, shuffledStdDerivation); if (costesPValue > 1.0) costesPValue = 1.0; else if (costesPValue < 0.0) costesPValue = 0.0; } + @Override public void processResults(ResultHandler handler) { super.processResults(handler); // if desired, show the last shuffled image available if ( showShuffledImages ) { handler.handleImage( smoothedShuffledImage, "Smoothed & shuffled channel 1" ); } handler.handleValue("Costes P-Value", costesPValue, 2); handler.handleValue("Costes Shuffled Mean", shuffledMean, 2); handler.handleValue("Costes Shuffled Std.D.", shuffledStdDerivation, 2); /* give the ratio of results at least as large as the * original value. */ double ratio = 0.0; if (shuffledPearsonsNotLessOriginal > 0) { ratio = (double)shuffledPearsonsResults.size() / (double)shuffledPearsonsNotLessOriginal; } handler.handleValue("Ratio of rand. Pearsons >= actual Pearsons value ", ratio, 2); } public double getCostesPValue() { return costesPValue; } public double getShuffledMean() { return shuffledMean; } public double getShuffledStdDerivation() { return shuffledStdDerivation; } public double getShuffledPearsonsNotLessOriginal() { return shuffledPearsonsNotLessOriginal; } } diff --git a/src/main/java/algorithms/Histogram2D.java b/src/main/java/algorithms/Histogram2D.java index 9ae3b7c..9a5a04b 100644 --- a/src/main/java/algorithms/Histogram2D.java +++ b/src/main/java/algorithms/Histogram2D.java @@ -1,378 +1,380 @@ package algorithms; import gadgets.DataContainer; import ij.measure.ResultsTable; import java.util.EnumSet; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.img.ImgFactory; import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.type.numeric.integer.LongType; import net.imglib2.view.Views; import results.ResultHandler; /** * Represents the creation of a 2D histogram between two images. * Channel 1 is set out in x direction, while channel 2 in y direction. * @param The source images value type */ public class Histogram2D> extends Algorithm { // An enumeration of possible drawings public enum DrawingFlags { Plot, RegressionLine, Axes } // the drawing configuration EnumSet drawingSettings; // The width of the scatter-plot protected int xBins = 256; // The height of the scatter-plot protected int yBins = 256; // The name of the result 2D histogram to pass elsewhere protected String title = ""; // Swap or not swap ch1 and ch2 protected boolean swapChannels = false; // member variables for labeling protected String ch1Label = "Channel 1"; protected String ch2Label = "Channel 2"; // Result keeping members // the generated plot image private RandomAccessibleInterval plotImage; // the bin widths for each channel private double xBinWidth = 0.0, yBinWidth = 0.0; // labels for the axes private String xLabel = "", yLabel = ""; // ranges for the axes private double xMin = 0.0, xMax = 0.0, yMin = 0.0, yMax = 0.0; public Histogram2D(){ this("2D Histogram"); } public Histogram2D(String title){ this(title, false); } public Histogram2D(String title, boolean swapChannels){ this(title, swapChannels, EnumSet.of( DrawingFlags.Plot, DrawingFlags.RegressionLine )); } public Histogram2D(String title, boolean swapChannels, EnumSet drawingSettings){ super(title); this.title = title; this.swapChannels = swapChannels; if (swapChannels) { int xBins = this.xBins; this.xBins = this.yBins; this.yBins = xBins; } this.drawingSettings = drawingSettings; } /** * Gets the minimum of channel one. Takes channel * swapping into consideration and will return min * of channel two if swapped. * * @return The minimum of what is seen as channel one. */ protected double getMinCh1(DataContainer container) { return swapChannels ? container.getMinCh2() : container.getMinCh1(); } /** * Gets the minimum of channel two. Takes channel * swapping into consideration and will return min * of channel one if swapped. * * @return The minimum of what is seen as channel two. */ protected double getMinCh2(DataContainer container) { return swapChannels ? container.getMinCh1() : container.getMinCh2(); } /** * Gets the maximum of channel one. Takes channel * swapping into consideration and will return max * of channel two if swapped. * * @return The maximum of what is seen as channel one. */ protected double getMaxCh1(DataContainer container) { return swapChannels ? container.getMaxCh2() : container.getMaxCh1(); } /** * Gets the maximum of channel two. Takes channel * swapping into consideration and will return max * of channel one if swapped. * * @return The maximum of what is seen as channel two. */ protected double getMaxCh2(DataContainer container) { return swapChannels ? container.getMaxCh1() : container.getMaxCh2(); } /** * Gets the image of channel one. Takes channel * swapping into consideration and will return image * of channel two if swapped. * * @return The image of what is seen as channel one. */ protected RandomAccessibleInterval getImageCh1(DataContainer container) { return swapChannels ? container.getSourceImage2() : container.getSourceImage1(); } /** * Gets the image of channel two. Takes channel * swapping into consideration and will return image * of channel one if swapped. * * @return The image of what is seen as channel two. */ protected RandomAccessibleInterval getImageCh2(DataContainer container) { return swapChannels ? container.getSourceImage1() : container.getSourceImage2(); } /** * Gets the label of channel one. Takes channel * swapping into consideration and will return label * of channel two if swapped. * * @return The label of what is seen as channel one. */ protected String getLabelCh1() { return swapChannels ? ch2Label : ch1Label; } /** * Gets the label of channel two. Takes channel * swapping into consideration and will return label * of channel one if swapped. * * @return The label of what is seen as channel two. */ protected String getLabelCh2() { return swapChannels ? ch1Label : ch2Label; } + @Override public void execute(DataContainer container) throws MissingPreconditionException { generateHistogramData(container); } protected void generateHistogramData(DataContainer container) { double ch1BinWidth = getXBinWidth(container); double ch2BinWidth = getYBinWidth(container); // get the 2 images for the calculation of Pearson's final RandomAccessibleInterval img1 = getImageCh1(container); final RandomAccessibleInterval img2 = getImageCh2(container); final RandomAccessibleInterval mask = container.getMask(); // get the cursors for iterating through pixels in images TwinCursor cursor = new TwinCursor(img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); // create new image to put the scatter-plot in final ImgFactory scatterFactory = new ArrayImgFactory< LongType >(); plotImage = scatterFactory.create(new int[] {xBins, yBins}, new LongType() ); // create access cursors final RandomAccess histogram2DCursor = plotImage.randomAccess(); // iterate over images long[] pos = new long[ plotImage.numDimensions() ]; while (cursor.hasNext()) { cursor.fwd(); double ch1 = cursor.getFirst().getRealDouble(); double ch2 = cursor.getSecond().getRealDouble(); /* Scale values for both channels to fit in the range. * Moreover mirror the y value on the x axis. */ pos[0] = getXValue(ch1, ch1BinWidth, ch2, ch2BinWidth); pos[1] = getYValue(ch1, ch1BinWidth, ch2, ch2BinWidth); // set position of input/output cursor histogram2DCursor.setPosition( pos ); // get current value at position and increment it long count = histogram2DCursor.get().getIntegerLong(); count++; histogram2DCursor.get().set(count); } xBinWidth = ch1BinWidth; yBinWidth = ch2BinWidth; xLabel = getLabelCh1(); yLabel = getLabelCh2(); xMin = getXMin(container); xMax = getXMax(container); yMin = getYMin(container); yMax = getYMax(container); } /** * A table of x-values, y-values and the counts is generated and * returned as a string. The single fields in one row (X Y Count) * are separated by tabs. * * @return A String representation of the histogram data. */ public String getData() { StringBuffer sb = new StringBuffer(); double xBinWidth = 1.0 / getXBinWidth(); double yBinWidth = 1.0 / getYBinWidth(); double xMin = getXMin(); double yMin = getYMin(); // check if we have bins of size one or other ones boolean xBinWidthIsOne = Math.abs(xBinWidth - 1.0) < 0.00001; boolean yBinWidthIsOne = Math.abs(yBinWidth - 1.0) < 0.00001; // configure decimal places accordingly int xDecimalPlaces = xBinWidthIsOne ? 0 : 3; int yDecimalPlaces = yBinWidthIsOne ? 0 : 3; // create a cursor to access the histogram data RandomAccess cursor = plotImage.randomAccess(); // loop over 2D histogram for (int i=0; i < plotImage.dimension(0); ++i) { for (int j=0; j < plotImage.dimension(1); ++j) { cursor.setPosition(i, 0); cursor.setPosition(j, 1); sb.append( ResultsTable.d2s(xMin + (i * xBinWidth), xDecimalPlaces) + "\t" + ResultsTable.d2s(yMin + (j * yBinWidth), yDecimalPlaces) + "\t" + ResultsTable.d2s(cursor.get().getRealDouble(), 0) + "\n"); } } return sb.toString(); } + @Override public void processResults(ResultHandler handler) { super.processResults(handler); handler.handleHistogram( this, title ); } /** * Calculates the bin width of one bin in x/ch1 direction. * @param container The container with images to work on * @return The width of one bin in x direction */ protected double getXBinWidth(DataContainer container) { double ch1Max = getMaxCh1(container); if (ch1Max < yBins) { // bin widths must not exceed 1 return 1; } // we need (ch1Max * width + 0.5) < yBins, but just so, i.e. // ch1Max * width + 0.5 == yBins - eps // width = (yBins - 0.5 - eps) / ch1Max return (yBins - 0.50001) / ch1Max; } /** * Calculates the bin width of one bin in y/ch2 direction. * @param container The container with images to work on * @return The width of one bin in y direction */ protected double getYBinWidth(DataContainer container) { double ch2Max = getMaxCh2(container); if (ch2Max < yBins) { // bin widths must not exceed 1 return 1; } return (yBins - 0.50001) / ch2Max; } /** * Calculates the locations x value. * @param ch1Val The intensity of channel one * @param ch1BinWidt The bin width for channel one * @return The x value of the data point location */ protected int getXValue(double ch1Val, double ch1BinWidth, double ch2Val, double ch2BinWidth) { return (int)(ch1Val * ch1BinWidth + 0.5); } /** * Calculates the locations y value. * @param ch2Val The intensity of channel one * @param ch2BinWidt The bin width for channel one * @return The x value of the data point location */ protected int getYValue(double ch1Val, double ch1BinWidth, double ch2Val, double ch2BinWidth) { return (yBins - 1) - (int)(ch2Val * ch2BinWidth + 0.5); } protected double getXMin(DataContainer container) { return 0; } protected double getXMax(DataContainer container) { return swapChannels ? getMaxCh2(container) : getMaxCh1(container); } protected double getYMin(DataContainer container) { return 0; } protected double getYMax(DataContainer container) { return swapChannels ? getMaxCh1(container) : getMaxCh2(container); } // Result access methods public RandomAccessibleInterval getPlotImage() { return plotImage; } public double getXBinWidth() { return xBinWidth; } public double getYBinWidth() { return yBinWidth; } public String getXLabel() { return xLabel; } public String getYLabel() { return yLabel; } public double getXMin() { return xMin; } public double getXMax() { return xMax; } public double getYMin() { return yMin; } public double getYMax() { return yMax; } public String getTitle() { return title; } public EnumSet getDrawingSettings() { return drawingSettings; } } diff --git a/src/main/java/algorithms/LiHistogram2D.java b/src/main/java/algorithms/LiHistogram2D.java index f20c93e..dd1a595 100644 --- a/src/main/java/algorithms/LiHistogram2D.java +++ b/src/main/java/algorithms/LiHistogram2D.java @@ -1,148 +1,152 @@ package algorithms; import gadgets.DataContainer; import java.util.EnumSet; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; /** * Represents the creation of a 2D histogram between two images. * Channel 1 is set out in x direction, while channel 2 in y direction. * The value calculation is done after Li. * @param */ public class LiHistogram2D> extends Histogram2D { // On execution these variables hold the images means double ch1Mean, ch2Mean; // On execution these variables hold the images min and max double ch1Min, ch1Max, ch2Min, ch2Max; // On execution these variables hold the Li's value min and max and their difference double liMin, liMax, liDiff; // On execution these variables hold the scaling factors double ch1Scaling, ch2Scaling; // boolean to test which channel we are using for eg. Li 2D histogram y axis boolean useCh1 = true; public LiHistogram2D(boolean useCh1) { this("Histogram 2D (Li)", useCh1); } public LiHistogram2D(String title, boolean useCh1) { this(title, false, useCh1); } public LiHistogram2D(String title, boolean swapChannels, boolean useCh1) { this(title, swapChannels, useCh1, EnumSet.of(DrawingFlags.Plot)); } public LiHistogram2D(String title, boolean swapChannels, boolean useCh1, EnumSet drawingSettings) { super(title, swapChannels, drawingSettings); this.useCh1 = useCh1; } + @Override public void execute(DataContainer container) throws MissingPreconditionException { ch1Mean = swapChannels ? container.getMeanCh2() : container.getMeanCh1(); ch2Mean = swapChannels ? container.getMeanCh1() : container.getMeanCh2(); ch1Min = getMinCh1(container); ch1Max = getMaxCh1(container); ch2Min = getMinCh2(container); ch2Max = getMaxCh2(container); /* A scaling to the x bins has to be made: * For that to work we need the min and the * max value that could occur. */ // get the 2 images and the mask final RandomAccessibleInterval img1 = getImageCh1(container); final RandomAccessibleInterval img2 = getImageCh2(container); final RandomAccessibleInterval mask = container.getMask(); // get the cursors for iterating through pixels in images TwinCursor cursor = new TwinCursor(img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); // give liMin and liMax appropriate starting values at the top and bottom of the range liMin = Double.MAX_VALUE; liMax = Double.MIN_VALUE; // iterate over images while (cursor.hasNext()) { cursor.fwd(); double ch1 = cursor.getFirst().getRealDouble(); double ch2 = cursor.getSecond().getRealDouble(); double productOfDifferenceOfMeans = (ch1Mean - ch1) * (ch2Mean - ch2); if (productOfDifferenceOfMeans < liMin) liMin = productOfDifferenceOfMeans; if (productOfDifferenceOfMeans > liMax) liMax = productOfDifferenceOfMeans; } liDiff = Math.abs(liMax - liMin); generateHistogramData(container); } @Override protected double getXBinWidth(DataContainer container) { return (double) xBins / (double)(liDiff + 1); } @Override protected double getYBinWidth(DataContainer container) { double max; if (useCh1) { max = getMaxCh1(container); } else { max = getMaxCh2(container); } return (double) yBins / (double)(max + 1); } @Override protected int getXValue(double ch1Val, double ch1BinWidth, double ch2Val, double ch2BinWidth) { /* We want the values to be scaled and shifted by and * offset in a way that the smallest (possibly negative) * value is in first bin and highest value in largest bin. */ return (int)( (( (ch1Mean - ch1Val) * (ch2Mean - ch2Val)) - liMin) * ch1BinWidth); } @Override protected int getYValue(double ch1Val, double ch1BinWidth, double ch2Val, double ch2BinWidth) { if (useCh1) return (yBins - 1) - (int)(ch1Val * ch2BinWidth); else return (yBins - 1) - (int)(ch2Val * ch2BinWidth); } @Override protected double getXMin(DataContainer container) { return swapChannels ? (useCh1 ? container.getMinCh1(): container.getMinCh2()) : liMin; } + @Override protected double getXMax(DataContainer container) { return swapChannels ? (useCh1 ? container.getMaxCh1(): container.getMaxCh2()) : liMax; } + @Override protected double getYMin(DataContainer container) { return swapChannels ? liMin : (useCh1 ? container.getMinCh1(): container.getMinCh2()); } + @Override protected double getYMax(DataContainer container) { return swapChannels ? liMax : (useCh1 ? container.getMaxCh1(): container.getMaxCh2()); } } diff --git a/src/main/java/algorithms/MandersColocalization.java b/src/main/java/algorithms/MandersColocalization.java index 533c25f..51acefa 100644 --- a/src/main/java/algorithms/MandersColocalization.java +++ b/src/main/java/algorithms/MandersColocalization.java @@ -1,262 +1,268 @@ package algorithms; import gadgets.DataContainer; import gadgets.ThresholdMode; import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; import results.ResultHandler; /** *

This algorithm calculates Manders et al.'s split colocalization * coefficients, M1 and M2. These are independent of signal intensities, * but are directly proportional to the amount of * fluorescence in the colocalized objects in each colour channel of the * image, relative to the total amount of fluorescence in that channel. * See "Manders, Verbeek, Aten - Measurement of colocalization * of objects in dual-colour confocal images. J. Microscopy, vol. 169 * pt 3, March 1993, pp 375-382".

* *

M1 = sum of Channel 1 intensity in pixels over the channel 2 threshold / total Channel 1 intensity. * M2 is vice versa. * The threshold may be everything > 0 in the other channel, which we call M1 and M2: without thresholds * or everything above some thresholds in the opposite channels 1 or 2, called tM1 and tM2: with thresholds * The result is a fraction (range 0-1, but often misrepresented as a %. We wont do that here.

* *

TODO: Further, it could/should/will calculate other split colocalization coefficients, * such as fraction of pixels (voxels) colocalized, * or fraction of intensity colocalized, as described at: * WCIF * copy pasted here - credits to Tony Collins.

* *

Number of colocalised voxels - Ncoloc * This is the number of voxels which have both channel 1 and channel 2 intensities above threshold * (i.e., the number of pixels in the yellow area of the scatterplot).

* *

%Image volume colocalised - %Volume * This is the percentage of voxels which have both channel 1 and channel 2 intensities above threshold, * expressed as a percentage of the total number of pixels in the image (including zero-zero pixels); * in other words, the number of pixels in the scatterplot's yellow area / total number of pixels in the scatter plot (the Red + Green + Blue + Yellow areas).

* *

%Voxels Colocalised - %Ch1 Vol; %Ch2 Vol * This generates a value for each channel. This is the number of voxels for each channel * which have both channel 1 and channel 2 intensities above threshold, * expressed as a percentage of the total number of voxels for each channel * above their respective thresholds; in other words, for channel 1 (along the x-axis), * this equals the (the number of pixels in the Yellow area) / (the number of pixels in the Blue + Yellow areas). * For channel 2 this is calculated as follows: * (the number of pixels in the Yellow area) / (the number of pixels in the Red + Yellow areas).

* *

%Intensity Colocalised - %Ch1 Int; %Ch2 Int * This generates a value for each channel. For channel 1, this value is equal to * the sum of the pixel intensities, with intensities above both channel 1 and channel 2 thresholds * expressed as a percentage of the sum of all channel 1 intensities; * in other words, it is calculated as follows: * (the sum of channel 1 pixel intensities in the Yellow area) / (the sum of channel 1 pixels intensities in the Red + Green + Blue + Yellow areas).

* *

%Intensities above threshold colocalised - %Ch1 Int > thresh; %Ch2 Int > thresh * This generates a value for each channel. For channel 1, * this value is equal to the sum of the pixel intensities * with intensities above both channel 1 and channel 2 thresholds * expressed as a percentage of the sum of all channel 1 intensities above the threshold for channel 1. * In other words, it is calculated as follows: * (the sum of channel 1 pixel intensities in the Yellow area) / (sum of channel 1 pixels intensities in the Blue + Yellow area)

* *

The results are often represented as % values, but to make them consistent with Manders' * split coefficients, we will also report them as fractions (range 0-1). * Perhaps this helps prevent the confusion in comparisons of results * caused by one person's %coloc being a totally * different measurement than another person's %coloc.

* * @param */ public class MandersColocalization> extends Algorithm { // Manders' split coefficients, M1 and M2: without thresholds // fraction of intensity of a channel, in pixels above zero in the other channel. private double mandersM1, mandersM2; // thresholded Manders M1 and M2 values, // Manders' split coefficients, tM1 and tM2: with thresholds // fraction of intensity of a channel, in pixels above threshold in the other channel. private double mandersThresholdedM1, mandersThresholdedM2; // Number of colocalized voxels (pixels) - Ncoloc private long numberOfPixelsAboveBothThresholds; // Fraction of Image volume colocalized - Fraction of Volume private double fractionOfPixelsAboveBothThresholds; // Fraction Voxels (pixels) Colocalized - Fraction of Ch1 Vol; Fraction of Ch2 Vol private double fractionOfColocCh1Pixels, fractionOfColocCh2Pixels; // Fraction Intensity Colocalized - Fraction of Ch1 Int; Fraction of Ch2 Int private double fractionOfColocCh1Intensity, fractionOfColocCh2Intensity; // Fraction of Intensities above threshold, colocalized - // Fraction of Ch1 Int > thresh; Fraction of Ch2 Int > thresh private double fractionOfColocCh1IntensityAboveCh1Thresh, fractionOfColocCh2IntensityAboveCh2Thresh; /** * A result container for Manders' calculations. */ public static class MandersResults { public double m1; public double m2; } public MandersColocalization() { super("Manders correlation"); } @Override public void execute(DataContainer container) throws MissingPreconditionException { // get the two images for the calculation of Manders' split coefficients RandomAccessible img1 = container.getSourceImage1(); RandomAccessible img2 = container.getSourceImage2(); RandomAccessibleInterval mask = container.getMask(); TwinCursor cursor = new TwinCursor(img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); // calculate Manders' split coefficients without threshold, M1 and M2. MandersResults results = calculateMandersCorrelation(cursor, img1.randomAccess().get().createVariable()); // save the results mandersM1 = results.m1; mandersM2 = results.m2; // calculate the thresholded Manders' split coefficients, tM1 and tM2, if possible AutoThresholdRegression autoThreshold = container.getAutoThreshold(); if (autoThreshold != null ) { // thresholded Manders' split coefficients, tM1 and tM2 cursor.reset(); results = calculateMandersCorrelation(cursor, autoThreshold.getCh1MaxThreshold(), autoThreshold.getCh2MaxThreshold(), ThresholdMode.Above); // save the results mandersThresholdedM1 = results.m1; mandersThresholdedM2 = results.m2; } } /** * Calculates Manders' split coefficients, M1 and M2: without thresholds * * @param cursor A TwinCursor that walks over two images * @param type A type instance, its value is not relevant * @return Both Manders' split coefficients, M1 and M2. */ public MandersResults calculateMandersCorrelation(TwinCursor cursor, T type) { return calculateMandersCorrelation(cursor, type, type, ThresholdMode.None); } /** * Calculates Manders' split coefficients, tM1 and tM2: with thresholds * * @param cursor A TwinCursor that walks over two images * @param type A type instance, its value is not relevant * @param thresholdCh1 type T * @param thresholdCh2 type T * @param tmode A ThresholdMode the threshold mode * @return Both thresholded Manders' split coefficients, tM1 and tM2. */ public MandersResults calculateMandersCorrelation(TwinCursor cursor, final T thresholdCh1, final T thresholdCh2, ThresholdMode tMode) { SplitCoeffAccumulator mandersAccum; // create a zero-values variable to compare to later on final T zero = thresholdCh1.createVariable(); zero.setZero(); // iterate over images - set the boolean value for if a pixel is thresholded // without thresholds: M1 and M1 if (tMode == ThresholdMode.None) { mandersAccum = new SplitCoeffAccumulator(cursor) { + @Override final boolean acceptMandersCh1(T type1, T type2) { return (type2.compareTo(zero) > 0); } + @Override final boolean acceptMandersCh2(T type1, T type2) { return (type1.compareTo(zero) > 0); } }; // with thresholds - below thresholds } else if (tMode == ThresholdMode.Below) { mandersAccum = new SplitCoeffAccumulator(cursor) { + @Override final boolean acceptMandersCh1(T type1, T type2) { return (type2.compareTo(zero) > 0) && (type2.compareTo(thresholdCh2) <= 0); } + @Override final boolean acceptMandersCh2(T type1, T type2) { return (type1.compareTo(zero) > 0) && (type1.compareTo(thresholdCh1) <= 0); } }; // with thresholds - above thresholds: tM1 and tM2 } else if (tMode == ThresholdMode.Above) { mandersAccum = new SplitCoeffAccumulator(cursor) { + @Override final boolean acceptMandersCh1(T type1, T type2) { return (type2.compareTo(zero) > 0) && (type2.compareTo(thresholdCh2) >= 0); } + @Override final boolean acceptMandersCh2(T type1, T type2) { return (type1.compareTo(zero) > 0) && (type1.compareTo(thresholdCh1) >= 0); } }; } else { throw new UnsupportedOperationException(); } MandersResults results = new MandersResults(); // calculate the results, see description above, as a fraction. results.m1 = mandersAccum.mandersSumCh1 / mandersAccum.sumCh1; results.m2 = mandersAccum.mandersSumCh2 / mandersAccum.sumCh2; return results; } @Override public void processResults(ResultHandler handler) { super.processResults(handler); handler.handleValue( "Manders' M1 (without thresholds)", mandersM1 ); handler.handleValue( "Manders' M2 (without thresholds)", mandersM2 ); handler.handleValue( "Manders' tM1 (with thresholds)", mandersThresholdedM1 ); handler.handleValue( "Manders' tM2 (with thresholds)", mandersThresholdedM2 ); } /** * A class similar to the Accumulator class, but more specific * to the Manders' split and other split channel coefficient calculations. */ protected abstract class SplitCoeffAccumulator { double sumCh1, sumCh2, mandersSumCh1, mandersSumCh2; public SplitCoeffAccumulator(TwinCursor cursor) { while (cursor.hasNext()) { cursor.fwd(); T type1 = cursor.getFirst(); T type2 = cursor.getSecond(); double ch1 = type1.getRealDouble(); double ch2 = type2.getRealDouble(); // boolean logics for adding or not adding to the different value counters for a pixel. if (acceptMandersCh1(type1, type2)) mandersSumCh1 += ch1; if (acceptMandersCh2(type1, type2)) mandersSumCh2 += ch2; // add this pixel's two intensity values to the ch1 and ch2 sum counters sumCh1 += ch1; sumCh2 += ch2; } } abstract boolean acceptMandersCh1(T type1, T type2); abstract boolean acceptMandersCh2(T type1, T type2); } } diff --git a/src/main/java/algorithms/PearsonsCorrelation.java b/src/main/java/algorithms/PearsonsCorrelation.java index c7a7203..4a750ba 100644 --- a/src/main/java/algorithms/PearsonsCorrelation.java +++ b/src/main/java/algorithms/PearsonsCorrelation.java @@ -1,382 +1,389 @@ package algorithms; import gadgets.DataContainer; import gadgets.MaskFactory; import gadgets.ThresholdMode; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.algorithm.math.ImageStatistics; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.RealType; import net.imglib2.view.Views; import results.ResultHandler; /** * A class that represents the mean calculation of the two source * images in the data container. * * @param */ public class PearsonsCorrelation> extends Algorithm { // Identifiers for choosing which implementation to use public enum Implementation {Classic, Fast}; // The member variable to store the implementation of the Pearson's Coefficient calculation used. Implementation theImplementation = Implementation.Fast; // resulting Pearsons value without thresholds double pearsonsCorrelationValue; // resulting Pearsons value below threshold double pearsonsCorrelationValueBelowThr; // resulting Pearsons value above threshold double pearsonsCorrelationValueAboveThr; /** * Creates a new Pearson's Correlation and allows us to define * which implementation of the calculation to use. * @param theImplementation The implementation of Pearson's Coefficient calculation to use. */ public PearsonsCorrelation(Implementation implementation) { super("Pearson correlation"); this.theImplementation = implementation; } /** * Creates a new Pearson's Correlation with default (fast) implementation parameter. */ public PearsonsCorrelation() { this(Implementation.Fast); } + @Override public void execute(DataContainer container) throws MissingPreconditionException { // get the 2 images for the calculation of Pearson's RandomAccessibleInterval img1 = container.getSourceImage1(); RandomAccessibleInterval img2 = container.getSourceImage2(); RandomAccessibleInterval mask = container.getMask(); // get the thresholds of the images AutoThresholdRegression autoThreshold = container.getAutoThreshold(); if (autoThreshold == null ) { throw new MissingPreconditionException("Pearsons calculation need AutoThresholdRegression to be run before it."); } T threshold1 = autoThreshold.getCh1MaxThreshold(); T threshold2 = autoThreshold.getCh2MaxThreshold(); if (threshold1 == null || threshold2 == null ) { throw new MissingPreconditionException("Pearsons calculation needs valid (not null) thresholds."); } /* Create cursors to walk over the images. First go over the * images without a mask. */ TwinCursor cursor = new TwinCursor( img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); MissingPreconditionException error = null; if (theImplementation == Implementation.Classic) { // get the means from the DataContainer double ch1Mean = container.getMeanCh1(); double ch2Mean = container.getMeanCh2(); try { cursor.reset(); pearsonsCorrelationValue = classicPearsons(cursor, ch1Mean, ch2Mean); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValue = Double.NaN; error = e; } try { cursor.reset(); pearsonsCorrelationValueBelowThr = classicPearsons(cursor, ch1Mean, ch2Mean, threshold1, threshold2, ThresholdMode.Below); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValueBelowThr = Double.NaN; error = e; } try { cursor.reset(); pearsonsCorrelationValueAboveThr = classicPearsons(cursor, ch1Mean, ch2Mean, threshold1, threshold2, ThresholdMode.Above); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValueAboveThr = Double.NaN; error = e; } } else if (theImplementation == Implementation.Fast) { try { cursor.reset(); pearsonsCorrelationValue = fastPearsons(cursor); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValue = Double.NaN; error = e; } try { cursor.reset(); pearsonsCorrelationValueBelowThr = fastPearsons(cursor, threshold1, threshold2, ThresholdMode.Below); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValueBelowThr = Double.NaN; error = e; } try { cursor.reset(); pearsonsCorrelationValueAboveThr = fastPearsons(cursor, threshold1, threshold2, ThresholdMode.Above); } catch (MissingPreconditionException e) { // probably a numerical error occurred pearsonsCorrelationValueAboveThr = Double.NaN; error = e; } } // if an error occurred, throw it one level up if (error != null) throw error; } /** * Calculates Pearson's R value without any constraint in values, thus it uses no thresholds. * If additional data like the images mean is needed, it is calculated. * * @param The images base type. * @param img1 The first image to walk over. * @param img2 The second image to walk over. * @return Pearson's R value. * @throws MissingPreconditionException */ public > double calculatePearsons( RandomAccessibleInterval img1, RandomAccessibleInterval img2) throws MissingPreconditionException { // create an "always true" mask to walk over the images final long[] dims = new long[img1.numDimensions()]; img1.dimensions(dims); RandomAccessibleInterval alwaysTrueMask = MaskFactory.createMask(dims, true); return calculatePearsons(img1, img2, alwaysTrueMask); } /** * Calculates Pearson's R value without any constraint in values, thus it uses no * thresholds. A mask is required to mark which data points should be visited. If * additional data like the images mean is needed, it is calculated. * * @param The images base type. * @param img1 The first image to walk over. * @param img2 The second image to walk over. * @param mask A mask for the images. * @return Pearson's R value. * @throws MissingPreconditionException */ public > double calculatePearsons( RandomAccessibleInterval img1, RandomAccessibleInterval img2, RandomAccessibleInterval mask) throws MissingPreconditionException { TwinCursor cursor = new TwinCursor( img1.randomAccess(), img2.randomAccess(), Views.iterable(mask).localizingCursor()); double r; if (theImplementation == Implementation.Classic) { /* since we need the means and apparently don't have them, * calculate them. */ double mean1 = ImageStatistics.getImageMean(img1); double mean2 = ImageStatistics.getImageMean(img2); // do the actual calculation r = classicPearsons(cursor, mean1, mean2); } else { r = fastPearsons(cursor); } return r; } /** * Calculates Pearson's R value with the possibility to constraint in values. * This could be useful of one wants to apply thresholds. You need to provide * the images means, albeit not used by all implementations. * * @param The images base type. * @param cursor The cursor to walk over both images. * @return Pearson's R value. * @throws MissingPreconditionException */ public > double calculatePearsons(TwinCursor cursor, double mean1, double mean2, S thresholdCh1, S thresholdCh2, ThresholdMode tMode) throws MissingPreconditionException { if (theImplementation == Implementation.Classic) { // do the actual calculation return classicPearsons(cursor, mean1, mean2, thresholdCh1, thresholdCh2, tMode); } else { return fastPearsons(cursor, thresholdCh1, thresholdCh2, tMode); } } /** * Calculates Person's R value by using a Classic implementation of the * algorithm. This method allows the specification of a TwinValueRangeCursor. * With such a cursor one for instance can combine different thresholding * conditions for each channel. The cursor is not closed in here. * * @param The image base type * @param img1 Channel one * @param img2 Channel two * @param cursor The cursor that defines the walk over both images. * @return Person's R value */ public static > double classicPearsons(TwinCursor cursor, double meanCh1, double meanCh2) throws MissingPreconditionException { return classicPearsons(cursor, meanCh1, meanCh2, null, null, ThresholdMode.None); } public static > double classicPearsons(TwinCursor cursor, double meanCh1, double meanCh2, final T thresholdCh1, final T thresholdCh2, ThresholdMode tMode) throws MissingPreconditionException { // the actual accumulation of the image values is done in a separate object Accumulator acc; if (tMode == ThresholdMode.None) { acc = new Accumulator(cursor, meanCh1, meanCh2) { + @Override final public boolean accept(T type1, T type2) { return true; } }; } else if (tMode == ThresholdMode.Below) { acc = new Accumulator(cursor, meanCh1, meanCh2) { + @Override final public boolean accept(T type1, T type2) { return type1.compareTo(thresholdCh1) < 0 || type2.compareTo(thresholdCh2) < 0; } }; } else if (tMode == ThresholdMode.Above) { acc = new Accumulator(cursor, meanCh1, meanCh2) { + @Override final public boolean accept(T type1, T type2) { return type1.compareTo(thresholdCh1) > 0 || type2.compareTo(thresholdCh2) > 0; } }; } else { throw new UnsupportedOperationException(); } double pearsonsR = acc.xy / Math.sqrt(acc.xx * acc.yy); checkForSanity(pearsonsR, acc.count); return pearsonsR; } /** * Calculates Person's R value by using a fast implementation of the * algorithm. This method allows the specification of a TwinValueRangeCursor. * With such a cursor one for instance can combine different thresholding * conditions for each channel. The cursor is not closed in here. * * @param The image base type * @param img1 Channel one * @param img2 Channel two * @param cursor The cursor that defines the walk over both images. * @return Person's R value */ public static > double fastPearsons(TwinCursor cursor) throws MissingPreconditionException { return fastPearsons(cursor, null, null, ThresholdMode.None); } public static > double fastPearsons(TwinCursor cursor, final T thresholdCh1, final T thresholdCh2, ThresholdMode tMode) throws MissingPreconditionException { // the actual accumulation of the image values is done in a separate object Accumulator acc; if (tMode == ThresholdMode.None) { acc = new Accumulator(cursor) { + @Override final public boolean accept(T type1, T type2) { return true; } }; } else if (tMode == ThresholdMode.Below) { acc = new Accumulator(cursor) { + @Override final public boolean accept(T type1, T type2) { return type1.compareTo(thresholdCh1) < 0 || type2.compareTo(thresholdCh2) < 0; } }; } else if (tMode == ThresholdMode.Above) { acc = new Accumulator(cursor) { + @Override final public boolean accept(T type1, T type2) { return type1.compareTo(thresholdCh1) > 0 || type2.compareTo(thresholdCh2) > 0; } }; } else { throw new UnsupportedOperationException(); } // for faster computation, have the inverse of N available double invCount = 1.0 / acc.count; double pearsons1 = acc.xy - (acc.x * acc.y * invCount); double pearsons2 = acc.xx - (acc.x * acc.x * invCount); double pearsons3 = acc.yy - (acc.y * acc.y * invCount); double pearsonsR = pearsons1 / (Math.sqrt(pearsons2 * pearsons3)); checkForSanity(pearsonsR, acc.count); return pearsonsR; } /** * Does a sanity check for calculated Pearsons values. Wrong * values can happen for fast and classic implementation. * * @param val The value to check. */ private static void checkForSanity(double value, int iterations) throws MissingPreconditionException { if ( Double.isNaN(value) || Double.isInfinite(value)) { /* For the _fast_ implementation this could happen: * Infinity could happen if only the numerator is 0, i.e.: * sum1squared == sum1 * sum1 * invN * and * sum2squared == sum2 * sum2 * invN * If the denominator is also zero, one will get NaN, i.e: * sumProduct1_2 == sum1 * sum2 * invN * * For the classic implementation it could happen, too: * Infinity happens if one channels sum of value-mean-differences * is zero. If it is negative for one image you will get NaN. * Additionally, if is zero for both channels at once you * could get NaN. NaN */ throw new MissingPreconditionException("A numerical problem occured: the input data is unsuitable for this algorithm. Possibly too few pixels (in range were: " + iterations + ")."); } } @Override public void processResults(ResultHandler handler) { super.processResults(handler); handler.handleValue("Pearson's R value (no threshold)", pearsonsCorrelationValue, 2); handler.handleValue("Pearson's R value (below threshold)", pearsonsCorrelationValueBelowThr, 2); handler.handleValue("Pearson's R value (above threshold)", pearsonsCorrelationValueAboveThr, 2); } public double getPearsonsCorrelationValue() { return pearsonsCorrelationValue; } public double getPearsonsCorrelationBelowThreshold() { return pearsonsCorrelationValueBelowThr; } public double getPearsonsCorrelationAboveThreshold() { return pearsonsCorrelationValueAboveThr; } } diff --git a/src/main/java/net/imglib2/TwinCursor.java b/src/main/java/net/imglib2/TwinCursor.java index 25ab861..465d080 100644 --- a/src/main/java/net/imglib2/TwinCursor.java +++ b/src/main/java/net/imglib2/TwinCursor.java @@ -1,149 +1,152 @@ package net.imglib2; import net.imglib2.predicate.MaskPredicate; import net.imglib2.predicate.Predicate; import net.imglib2.type.Type; import net.imglib2.type.logic.BitType; /** * The TwinCursor moves over two images with respect to a mask. The mask * has to be of the same dimensionality as the images. Position information * obtained from this class comes from the mask. * * @author Johannes Schindelin and Tom Kazimiers */ public class TwinCursor> implements Cursor, PairIterator { final protected PredicateCursor mask; final protected RandomAccess channel1; final protected RandomAccess channel2; /* * For performance, we keep one position array (to avoid * having to create a new array in every single step). */ final protected long[] position; /* To avoid calling next() too often */ protected boolean gotNext; public TwinCursor(final RandomAccess channel1, final RandomAccess channel2, final Cursor mask) { final Predicate predicate = new MaskPredicate(); this.mask = new PredicateCursor(mask, predicate); this.channel1 = channel1; this.channel2 = channel2; position = new long[mask.numDimensions()]; mask.localize(position); } + @Override final public boolean hasNext() { gotNext = false; return mask.hasNext(); } final public void getNext() { if (gotNext) return; mask.next(); mask.localize(position); channel1.setPosition(position); channel2.setPosition(position); gotNext = true; } + @Override final public T getFirst() { getNext(); return channel1.get(); } + @Override final public T getSecond() { getNext(); return channel2.get(); } @Override public void reset() { gotNext = false; mask.reset(); } @Override public void fwd() { if (hasNext()) getNext(); } @Override public void jumpFwd(long arg0) { throw new UnsupportedOperationException("This method has not been implemented, yet."); } @Override public T next() { throw new UnsupportedOperationException("This method has not been implemented, yet."); } @Override public void remove() { throw new UnsupportedOperationException("This method has not been implemented, yet."); } @Override public double getDoublePosition(int arg0) { return mask.getDoublePosition(arg0); } @Override public float getFloatPosition(int arg0) { return mask.getFloatPosition(arg0); } @Override public void localize(float[] arg0) { mask.localize(arg0); } @Override public void localize(double[] arg0) { mask.localize(arg0); } @Override public int numDimensions() { return mask.numDimensions(); } @Override public Sampler copy() { throw new UnsupportedOperationException("This method has not been implemented, yet."); } @Override public T get() { throw new UnsupportedOperationException("This method has not been implemented, yet."); } @Override public int getIntPosition(int arg0) { return mask.getIntPosition(arg0); } @Override public long getLongPosition(int arg0) { return mask.getLongPosition(arg0); } @Override public void localize(int[] arg0) { mask.localize(arg0); } @Override public void localize(long[] arg0) { mask.localize(arg0); } @Override public Cursor copyCursor() { throw new UnsupportedOperationException("This method has not been implemented, yet."); } } \ No newline at end of file diff --git a/src/main/java/net/imglib2/predicate/MaskPredicate.java b/src/main/java/net/imglib2/predicate/MaskPredicate.java index 9a5abb8..9c4e95e 100644 --- a/src/main/java/net/imglib2/predicate/MaskPredicate.java +++ b/src/main/java/net/imglib2/predicate/MaskPredicate.java @@ -1,11 +1,12 @@ package net.imglib2.predicate; import net.imglib2.Cursor; import net.imglib2.type.logic.BitType; public class MaskPredicate implements Predicate { + @Override public boolean test(final Cursor cursor) { return cursor.get().get(); } -} \ No newline at end of file +} diff --git a/src/main/java/results/EasyDisplay.java b/src/main/java/results/EasyDisplay.java index cd4c2c5..e1b7e56 100644 --- a/src/main/java/results/EasyDisplay.java +++ b/src/main/java/results/EasyDisplay.java @@ -1,90 +1,96 @@ package results; import algorithms.Histogram2D; import gadgets.DataContainer; import ij.IJ; import ij.ImagePlus; import ij.text.TextWindow; import net.imglib2.RandomAccessibleInterval; import net.imglib2.algorithm.math.ImageStatistics; import net.imglib2.img.display.imagej.ImageJFunctions; import net.imglib2.type.numeric.RealType; public class EasyDisplay> implements ResultHandler { // the text window to present value and text results protected static TextWindow textWindow; /* the data container with general information about the * source images that were processed by the algorithms. */ protected DataContainer container; public EasyDisplay(DataContainer container) { final int twWidth = 170; final int twHeight = 250; //test if the results windows is already there, if so use it. if (textWindow == null || !textWindow.isVisible()) textWindow = new TextWindow("Results", "Result\tValue\n", "", twWidth, twHeight); else { // set dimensions textWindow.setSize(twWidth, twHeight); } // deactivate the window for now textWindow.setVisible(false); // save a reference to the data container this.container = container; } + @Override public void handleImage(RandomAccessibleInterval image, String name) { ImagePlus imp = ImageJFunctions.wrapFloat( image, name ); double max = ImageStatistics.getImageMax( image ).getRealDouble(); showImage( imp, max ); } + @Override public void handleHistogram(Histogram2D histogram, String name) { ImagePlus imp = ImageJFunctions.wrapFloat( histogram.getPlotImage(), name ); double max = ImageStatistics.getImageMax( histogram.getPlotImage() ).getRealDouble(); showImage( imp, max ); } protected void showImage(ImagePlus imp, double max) { // set the display range imp.setDisplayRange(0.0, max); imp.show(); } + @Override public void handleWarning(Warning warning) { // no warnings are shown in easy display } @Override public void handleValue(String name, String value) { textWindow.getTextPanel().appendLine(name + "\t" + value + "\n"); } + @Override public void handleValue(String name, double value) { handleValue(name, value, 3); } + @Override public void handleValue(String name, double value, int decimals) { handleValue(name, IJ.d2s(value, decimals)); } protected void printTextStatistics(DataContainer container){ textWindow.getTextPanel().appendLine("Ch1 Mean\t" + container.getMeanCh1() + "\n"); textWindow.getTextPanel().appendLine("Ch2 Mean\t" + container.getMeanCh2() + "\n"); textWindow.getTextPanel().appendLine("Ch1 Min\t" + container.getMinCh1() + "\n"); textWindow.getTextPanel().appendLine("Ch2 Min\t" + container.getMinCh2() + "\n"); textWindow.getTextPanel().appendLine("Ch1 Max\t" + container.getMaxCh1() + "\n"); textWindow.getTextPanel().appendLine("Ch2 Max\t" + container.getMaxCh2() + "\n"); } + @Override public void process() { // print some general information about images printTextStatistics(container); // show the results textWindow.setVisible(true); IJ.selectWindow("Results"); } } diff --git a/src/main/java/results/NamedContainer.java b/src/main/java/results/NamedContainer.java index ffcc8cc..cd1727b 100644 --- a/src/main/java/results/NamedContainer.java +++ b/src/main/java/results/NamedContainer.java @@ -1,24 +1,25 @@ package results; /** * A small container to name objects by overriding * the toString method. * */ public class NamedContainer { T object; String name; public NamedContainer(T object, String name) { this.object = object; this.name = name; } public T getObject() { return object; } + @Override public String toString() { return name; } } diff --git a/src/main/java/results/PDFWriter.java b/src/main/java/results/PDFWriter.java index 749124e..acba2b4 100644 --- a/src/main/java/results/PDFWriter.java +++ b/src/main/java/results/PDFWriter.java @@ -1,229 +1,235 @@ package results; import algorithms.Histogram2D; import com.itextpdf.text.BadElementException; import com.itextpdf.text.Document; import com.itextpdf.text.DocumentException; import com.itextpdf.text.PageSize; import com.itextpdf.text.Paragraph; import com.itextpdf.text.pdf.PdfContentByte; import com.itextpdf.text.pdf.PdfWriter; import gadgets.DataContainer; import gadgets.DataContainer.MaskType; import ij.IJ; import ij.ImagePlus; import ij.io.SaveDialog; import java.io.FileOutputStream; import java.io.IOException; import java.util.ArrayList; import java.util.List; import net.imglib2.RandomAccessibleInterval; import net.imglib2.algorithm.math.ImageStatistics; import net.imglib2.img.display.imagej.ImageJFunctions; import net.imglib2.type.numeric.RealType; import net.imglib2.type.numeric.integer.LongType; public class PDFWriter> implements ResultHandler { // indicates if we want to produce US letter or A4 size boolean isLetter = false; // indicates if the content is the first item on the page boolean isFirst = true; // show the name of the image static boolean showName=true; // a static counter for this sessions created PDFs static int succeededPrints = 0; // show the size in pixels of the image static boolean showSize=true; // a reference to the data container DataContainer container; PdfWriter writer; Document document; // a list of the available result images, no matter what specific kinds protected List listOfPDFImages = new ArrayList(); protected List listOfPDFTexts = new ArrayList(); /** * Creates a new PDFWriter that can access the container. * * @param container The data container for source image data */ public PDFWriter(DataContainer container) { this.container = container; } + @Override public void handleImage(RandomAccessibleInterval image, String name) { ImagePlus imp = ImageJFunctions.wrapFloat( image, name ); // set the display range double max = ImageStatistics.getImageMax(image).getRealDouble(); imp.setDisplayRange(0.0, max); addImageToList(imp, name); } /** * Handles a histogram the following way: create snapshot, log data, reset the * display range, apply the Fire LUT and finally store it as an iText PDF image. * Afterwards the image is reset to its orignal state again */ + @Override public void handleHistogram(Histogram2D histogram, String name) { RandomAccessibleInterval image = histogram.getPlotImage(); ImagePlus imp = ImageJFunctions.wrapFloat( image, name ); // make a snapshot to be able to reset after modifications imp.getProcessor().snapshot(); imp.getProcessor().log(); imp.updateAndDraw(); imp.getProcessor().resetMinAndMax(); IJ.run(imp,"Fire", null); addImageToList(imp, name); // reset the imp from the log scaling we applied earlier imp.getProcessor().reset(); } protected void addImageToList(ImagePlus imp, String name) { java.awt.Image awtImage = imp.getImage(); try { com.itextpdf.text.Image pdfImage = com.itextpdf.text.Image.getInstance(awtImage, null); pdfImage.setAlt(name); // iText-1.3 setMarkupAttribute("name", name); listOfPDFImages.add(pdfImage); } catch (BadElementException e) { IJ.log("Could not convert image to correct format for PDF generation"); IJ.handleException(e); } catch (IOException e) { IJ.log("Could not convert image to correct format for PDF generation"); IJ.handleException(e); } } + @Override public void handleWarning(Warning warning) { listOfPDFTexts.add(new Paragraph("Warning! " + warning.getShortMessage() + " - " + warning.getLongMessage())); } @Override public void handleValue(String name, String value) { // send (output parameter name, value) to IJ.log for scraping batch results IJ.log(name + ", " + value); listOfPDFTexts.add(new Paragraph(name + ": " + value)); } + @Override public void handleValue(String name, double value) { handleValue(name, value, 3); } + @Override public void handleValue(String name, double value, int decimals) { listOfPDFTexts.add(new Paragraph(name + ": " + IJ.d2s(value, decimals))); } /** * Prints an image into the opened PDF. * @param img The image to print. * @param printName The name to print under the image. */ protected void addImage(com.itextpdf.text.Image image) throws DocumentException, IOException { if (! isFirst) { document.add(new Paragraph("\n")); float vertPos = writer.getVerticalPosition(true); if (vertPos - document.bottom() < image.getHeight()) { document.newPage(); } else { PdfContentByte cb = writer.getDirectContent(); cb.setLineWidth(1f); if (isLetter) { cb.moveTo(PageSize.LETTER.getLeft(50), vertPos); cb.lineTo(PageSize.LETTER.getRight(50), vertPos); } else { cb.moveTo(PageSize.A4.getLeft(50), vertPos); cb.lineTo(PageSize.A4.getRight(50), vertPos); } cb.stroke(); } } if (showName) { Paragraph paragraph = new Paragraph(image.getAlt()); // iText-1.3: getMarkupAttribute("name")); paragraph.setAlignment(Paragraph.ALIGN_CENTER); document.add(paragraph); //spcNm = 40; } if (showSize) { Paragraph paragraph = new Paragraph(image.getWidth() + " x " + image.getHeight()); paragraph.setAlignment(Paragraph.ALIGN_CENTER); document.add(paragraph); //spcSz = 40; } image.setAlignment(com.itextpdf.text.Image.ALIGN_CENTER); document.add(image); isFirst = false; } + @Override public void process() { try { // Use the getJobName() in DataContainer for the job name. String jobName = container.getJobName(); /* If a mask is in use, add a counter * information to the jobName. */ if (container.getMaskType() != MaskType.None) { // maskHash is now used as the mask or ROI unique ID in the // jobName but we can still increment and use succeededPrints at // the end of the filename for PDFs when there is a mask. jobName += (succeededPrints + 1); } // get the path to the file we are about to create SaveDialog sd = new SaveDialog("Save as PDF", jobName, ".pdf"); // update jobName if the user changes it in the save file dialog. jobName = sd.getFileName(); String directory = sd.getDirectory(); // make sure we have what we need next if ((jobName == null) || (directory == null)) { return; } String path = directory+jobName; // create a new iText Document and add date and title document = new Document(isLetter ? PageSize.LETTER : PageSize.A4); document.addCreationDate(); document.addTitle(jobName); // get a writer object to do the actual output writer = PdfWriter.getInstance(document, new FileOutputStream(path)); document.open(); // write job name at the top of the PDF file as a title Paragraph titlePara = new Paragraph(jobName); document.add(titlePara); // iterate over all produced images for (com.itextpdf.text.Image img : listOfPDFImages) { addImage(img); } //iterate over all produced text objects for (Paragraph p : listOfPDFTexts) { document.add(p); } } catch(DocumentException de) { IJ.showMessage("PDF Writer", de.getMessage()); } catch(IOException ioe) { IJ.showMessage("PDF Writer", ioe.getMessage()); } finally { if (document !=null) { document.close(); succeededPrints++; } } } } diff --git a/src/main/java/results/SingleWindowDisplay.java b/src/main/java/results/SingleWindowDisplay.java index d7214e3..c1ba9a2 100644 --- a/src/main/java/results/SingleWindowDisplay.java +++ b/src/main/java/results/SingleWindowDisplay.java @@ -1,635 +1,647 @@ package results; import algorithms.AutoThresholdRegression; import algorithms.Histogram2D; import fiji.util.gui.JImagePanel; import gadgets.DataContainer; import ij.IJ; import ij.ImageJ; import ij.ImagePlus; import ij.gui.Line; import ij.gui.Overlay; import ij.process.ImageProcessor; import ij.text.TextWindow; import java.awt.Container; import java.awt.Dimension; import java.awt.FlowLayout; import java.awt.GridBagConstraints; import java.awt.GridBagLayout; import java.awt.Panel; import java.awt.datatransfer.Clipboard; import java.awt.datatransfer.ClipboardOwner; import java.awt.datatransfer.StringSelection; import java.awt.datatransfer.Transferable; import java.awt.event.ActionEvent; import java.awt.event.ActionListener; import java.awt.event.ItemEvent; import java.awt.event.ItemListener; import java.awt.event.MouseEvent; import java.awt.event.MouseMotionListener; import java.io.PrintWriter; import java.io.StringWriter; import java.util.ArrayList; import java.util.HashMap; import java.util.List; import java.util.Map; import javax.swing.JButton; import javax.swing.JCheckBox; import javax.swing.JComboBox; import javax.swing.JEditorPane; import javax.swing.JFrame; import javax.swing.JScrollPane; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessibleInterval; import net.imglib2.img.display.imagej.ImageJFunctions; import net.imglib2.type.numeric.RealType; import net.imglib2.type.numeric.integer.LongType; /** * This class displays the container contents in one single window * and offers features like the use of different LUTs. * */ public class SingleWindowDisplay> extends JFrame implements ResultHandler, ItemListener, ActionListener, ClipboardOwner, MouseMotionListener { private static final long serialVersionUID = -5642321584354176878L; protected static final int WIN_WIDTH = 350; protected static final int WIN_HEIGHT = 600; // a static list for keeping track of all other SingleWindowDisplays protected static ArrayList> displays = new ArrayList>(); // indicates if original images should be displayed or not protected boolean displayOriginalImages = false; // this is the image currently selected by the drop down menu protected RandomAccessibleInterval> currentlyDisplayedImageResult; // a list of the available result images, no matter what specific kinds protected List< NamedContainer< RandomAccessibleInterval>> > listOfImages = new ArrayList< NamedContainer< RandomAccessibleInterval>> >(); protected Map, Histogram2D> mapOf2DHistograms = new HashMap, Histogram2D>(); // a list of warnings protected List warnings = new ArrayList(); // a list of named values, collected from algorithms protected List valueResults = new ArrayList(); /* a map of images and corresponding LUTs. When an image is not in * there no LUT should be applied. */ protected Map listOfLUTs = new HashMap(); //make a cursor so we can get pixel values from the image protected RandomAccess> pixelAccessCursor; // A PDF writer to call if user wants PDF print protected PDFWriter pdfWriter; // The current image protected ImagePlus imp; // GUI elements protected JImagePanel imagePanel; protected JButton listButton, copyButton; protected JCheckBox log; /* The data container with general information about * source images */ protected DataContainer dataContainer = null; public SingleWindowDisplay(DataContainer container, PDFWriter pdfWriter){ // Show job name in title bar super(container.getJobName()); setPreferredSize(new Dimension(WIN_WIDTH, WIN_HEIGHT)); // save a reference to the container dataContainer = container; this.pdfWriter = pdfWriter; // don't show ourself on instantiation this.setVisible(false); // add ourself to the list of single window displays displays.add(this); } public void setup() { JComboBox dropDownList = new JComboBox(); for(NamedContainer< RandomAccessibleInterval> > img : listOfImages) { dropDownList.addItem( new NamedContainer> >( img.object, img.name)); } dropDownList.addItemListener(this); imagePanel = new JImagePanel(ij.IJ.createImage("dummy", "8-bit", 10, 10, 1)); imagePanel.addMouseMotionListener(this); // Create something to display it in final JEditorPane editor = new JEditorPane(); editor.setEditable(false); // we're browsing not editing editor.setContentType("text/html"); // must specify HTML text editor.setText(makeHtmlText()); // specify the text to display // Put the JEditorPane in a scrolling window and add it JScrollPane scrollPane = new JScrollPane(editor); scrollPane.setPreferredSize(new Dimension(256, 150)); Panel buttons = new Panel(); buttons.setLayout(new FlowLayout(FlowLayout.RIGHT)); // add button for data display of histograms listButton = new JButton("List"); listButton.addActionListener(new ActionListener() { + @Override public void actionPerformed(ActionEvent e) { showList(); } }); buttons.add(listButton); // add button for data copy of histograms copyButton = new JButton("Copy"); copyButton.addActionListener(new ActionListener() { + @Override public void actionPerformed(ActionEvent arg0) { copyToClipboard(); } }); buttons.add(copyButton); // add button for PDF printing JButton pdfButten = new JButton("PDF"); pdfButten.addActionListener(new ActionListener() { + @Override public void actionPerformed(ActionEvent arg0) { pdfWriter.process(); } }); buttons.add(pdfButten); /* We want the image to be log scale by default * so the user can see something. */ log = new JCheckBox("Log"); log.setSelected(true); log.addActionListener(this); buttons.add(log); final GridBagLayout layout = new GridBagLayout(); final Container pane = getContentPane(); getContentPane().setLayout(layout); final GridBagConstraints c = new GridBagConstraints(); c.fill = GridBagConstraints.BOTH; c.weightx = 1; c.gridwidth = GridBagConstraints.BOTH; c.gridy++; pane.add(dropDownList, c); c.gridy++; c.weighty = 1; pane.add(imagePanel, c); c.gridy++; c.weighty = 1; pane.add(scrollPane, c); c.weighty = 0; c.gridy++; pane.add(buttons, c); pack(); } + @Override public void process() { // if wanted, display source images if ( displayOriginalImages ) { listOfImages.add( new NamedContainer> >( dataContainer.getSourceImage1(), dataContainer.getSourceImage1Name() ) ); listOfImages.add( new NamedContainer> >( dataContainer.getSourceImage2(), dataContainer.getSourceImage2Name() ) ); } // set up the GUI, which runs makeHtmlText() for the value results // formatting. setup(); // display the first image available, if any if (listOfImages.size() > 0) { adjustDisplayedImage(listOfImages.get(0).object); } // show the GUI this.setVisible(true); } + @Override public void handleImage(RandomAccessibleInterval image, String name) { listOfImages.add( new NamedContainer> >( image, name ) ); } + @Override public void handleHistogram(Histogram2D histogram, String name) { listOfImages.add( new NamedContainer> >( histogram.getPlotImage(), name ) ); mapOf2DHistograms.put(histogram.getPlotImage(), histogram); // link the histogram to a LUT listOfLUTs.put(histogram.getPlotImage(), "Fire"); } + @Override public void handleWarning(Warning warning) { warnings.add( warning ); } @Override public void handleValue(String name, String value) { valueResults.add( new ValueResult(name, value)); } + @Override public void handleValue(String name, double value) { handleValue(name, value, 3); } + @Override public void handleValue(String name, double value, int decimals) { valueResults.add( new ValueResult(name, value, decimals)); } /** * Prints an HTML table entry onto the stream. */ protected void printTableRow(PrintWriter out, String name, String text) { out.print("" + name + "" + text + ""); // add some style information out.print(""); out.print(""); // print out warnings, if any if ( warnings.size() > 0 ) { out.print("

Warnings

"); // Print out the table out.print(""); out.print(""); for (Warning w : warnings) { printTableRow(out, w.getShortMessage(), w.getLongMessage()); } out.println("
TypeMessage
"); } else { out.print("

No warnings occured

"); } // print out simple value results out.print("

Results

"); // Print out the table out.print(""); out.print(""); // Print table rows and spit results to the IJ log. for ( ValueResult vr : valueResults) { if (vr.isNumber) { printTableRow(out, vr.name, vr.number, vr.decimals); IJ.log(vr.name + ", " + IJ.d2s(vr.number, vr.decimals)); } else { printTableRow(out, vr.name, vr.value); IJ.log(vr.name + ", " + vr.value); } } out.println("
NameResult
"); out.print(""); out.close(); // Get the string of HTML from the StringWriter and return it. return sout.toString(); } /** * If the currently selected ImageResult is an HistrogramResult, * a table of x-values, y-values and the counts. */ protected void showList() { /* check if we are dealing with an histogram result * or a generic image result */ if (isHistogram(currentlyDisplayedImageResult)) { Histogram2D hr = mapOf2DHistograms.get(currentlyDisplayedImageResult); double xBinWidth = 1.0 / hr.getXBinWidth(); double yBinWidth = 1.0 / hr.getYBinWidth(); // check if we have bins of size one or other ones boolean xBinWidthIsOne = Math.abs(xBinWidth - 1.0) < 0.00001; boolean yBinWidthIsOne = Math.abs(yBinWidth - 1.0) < 0.00001; // configure table headings accordingly String vHeadingX = xBinWidthIsOne ? "X value" : "X bin start"; String vHeadingY = yBinWidthIsOne ? "Y value" : "Y bin start"; // get the actual histogram data String histogramData = hr.getData(); TextWindow tw = new TextWindow(getTitle(), vHeadingX + "\t" + vHeadingY + "\tcount", histogramData, 250, 400); tw.setVisible(true); } } /** * If the currently selected ImageResult is an HistogramRestult, * this method copies its data into to the clipboard. */ protected void copyToClipboard() { /* check if we are dealing with an histogram result * or a generic image result */ if (isHistogram(currentlyDisplayedImageResult)) { /* try to get the system clipboard and return * if we can't get it */ Clipboard systemClipboard = null; try { systemClipboard = getToolkit().getSystemClipboard(); } catch (Exception e) { systemClipboard = null; } if (systemClipboard==null) { IJ.error("Unable to copy to Clipboard."); return; } // copy histogram values IJ.showStatus("Copying histogram values..."); String text = mapOf2DHistograms.get(currentlyDisplayedImageResult).getData(); StringSelection contents = new StringSelection( text ); systemClipboard.setContents(contents, this); IJ.showStatus(text.length() + " characters copied to Clipboard"); } } @Override public void mouseDragged(MouseEvent e) { // nothing to do here } @Override public void mouseMoved(MouseEvent e) { if (e.getSource().equals(imagePanel)) { /* * calculate the mouse position relative to the upper left * corner of the displayed image. */ final int imgWidth = imagePanel.getSrcRect().width; final int imgHeight = imagePanel.getSrcRect().height; int displayWidth = (int)(imgWidth * imagePanel.getMagnification()); int displayHeight = (int)(imgHeight * imagePanel.getMagnification()); int offsetX = (imagePanel.getWidth() - displayWidth) / 2; int offsetY = (imagePanel.getHeight() - displayHeight) / 2; int onImageX = imagePanel.screenX(e.getX() - offsetX); int onImageY = imagePanel.screenY(e.getY() - offsetY); // make sure we stay within the image boundaries if (onImageX >= 0 && onImageX < imgWidth && onImageY >= 0 && onImageY < imgHeight ) { mouseMoved(onImageX, onImageY); } else { IJ.showStatus(""); } } } /** * Displays information about the pixel below the mouse cursor of * the currently displayed image result. The coordinates passed are * expected to be within the image boundaries. * * @param x * @param y */ public void mouseMoved( int x, int y) { final ImageJ ij = IJ.getInstance(); if (ij != null && currentlyDisplayedImageResult != null) { /* If Alt key is not pressed, display the calibrated data. * If not, display image positions and data. * Non log image intensity from original image or 2D histogram result is always shown in status bar, * not the log intensity that might actually be displayed in the image. */ if (!IJ.altKeyDown()) { // the alt key is not pressed use x and y values that are bin widths or calibrated intensities not the x y image coordinates. if (isHistogram(currentlyDisplayedImageResult)) { Histogram2D histogram = mapOf2DHistograms.get(currentlyDisplayedImageResult); synchronized( pixelAccessCursor ) { // set position of output cursor pixelAccessCursor.setPosition(x, 0); pixelAccessCursor.setPosition(y, 1); // for a histogram coordinate display we need to invert the Y axis y = (int)currentlyDisplayedImageResult.dimension(1) - 1 - y; // get current value at position RandomAccess cursor = (RandomAccess)pixelAccessCursor; long val = cursor.get().getIntegerLong(); double calibratedXBinBottom = histogram.getXMin() + x / histogram.getXBinWidth(); double calibratedXBinTop = histogram.getXMin() + (x + 1) / histogram.getXBinWidth(); double calibratedYBinBottom = histogram.getYMin() + y / histogram.getYBinWidth(); double calibratedYBinTop = histogram.getYMin() + (y + 1) / histogram.getYBinWidth(); IJ.showStatus("x = " + IJ.d2s(calibratedXBinBottom) + " to " + IJ.d2s(calibratedXBinTop) + ", y = " + IJ.d2s(calibratedYBinBottom) + " to " + IJ.d2s(calibratedYBinTop) + ", value = " + val ); } } else { RandomAccessibleInterval img = (RandomAccessibleInterval) currentlyDisplayedImageResult; ImagePlus imp = ImageJFunctions.wrapFloat( img, "TODO" ); imp.mouseMoved(x, y); } } else { // alt key is down, so show the image coordinates for x y in status bar. RandomAccessibleInterval img = (RandomAccessibleInterval) currentlyDisplayedImageResult; ImagePlus imp = ImageJFunctions.wrapFloat( img, "TODO" ); imp.mouseMoved(x, y); } } } /** * Draws the passed ImageResult on the ImagePlus of this class. * If the image is part of a CompositeImageResult then contained * lines will also be drawn */ protected void drawImage(RandomAccessibleInterval> img) { // get ImgLib image as ImageJ image imp = ImageJFunctions.wrapFloat( (RandomAccessibleInterval) img, "TODO" ); imagePanel.updateImage(imp); // set the display range // check if a LUT should be applied if ( listOfLUTs.containsKey(img) ) { // select linked look up table IJ.run(imp, listOfLUTs.get(img), null); } imp.getProcessor().resetMinAndMax(); boolean overlayModified = false; Overlay overlay = new Overlay(); // if it is the 2d histogram, we want to show the regression line if (isHistogram(img)) { Histogram2D histogram = mapOf2DHistograms.get(img); /* check if we should draw a regression line for the * current histogram. */ if ( histogram.getDrawingSettings().contains(Histogram2D.DrawingFlags.RegressionLine) ) { AutoThresholdRegression autoThreshold = dataContainer.getAutoThreshold(); if (histogram != null && autoThreshold != null) { if (img == histogram.getPlotImage()) { drawLine(overlay, img, autoThreshold.getAutoThresholdSlope(), autoThreshold.getAutoThresholdIntercept()); overlayModified = true; } } } } if (overlayModified) { overlay.setStrokeColor(java.awt.Color.WHITE); imp.setOverlay(overlay); } imagePanel.repaint(); } /** * Tests whether the given image is a histogram or not. * @param img The image to test * @return true if histogram, false otherwise */ protected boolean isHistogram(RandomAccessibleInterval> img) { return mapOf2DHistograms.containsKey(img); } /** * Draws the line on the overlay. */ protected void drawLine(Overlay overlay, RandomAccessibleInterval> img, double slope, double intercept) { double startX, startY, endX, endY; long imgWidth = img.dimension(0); long imgHeight = img.dimension(1); /* since we want to draw the line over the whole image * we can directly use screen coordinates for x values. */ startX = 0.0; endX = imgWidth; // check if we can get some exta information for drawing if (isHistogram(img)) { Histogram2D histogram = mapOf2DHistograms.get(img); // get calibrated start y coordinates double calibratedStartY = slope * histogram.getXMin() + intercept; double calibratedEndY = slope * histogram.getXMax() + intercept; // convert calibrated coordinates to screen coordinates startY = calibratedStartY * histogram.getYBinWidth(); endY = calibratedEndY * histogram.getYBinWidth(); } else { startY = slope * startX + intercept; endY = slope * endX + intercept; } /* since the screen origin is in the top left * of the image, we need to x-mirror our line */ startY = ( imgHeight - 1 ) - startY; endY = ( imgHeight - 1 ) - endY; // create the line ROI and add it to the overlay Line lineROI = new Line(startX, startY, endX, endY); /* Set drawing width of line to one, in case it has * been changed globally. */ lineROI.setStrokeWidth(1.0f); overlay.add(lineROI); } protected void adjustDisplayedImage (RandomAccessibleInterval> img) { /* when changing the result image to display * need to set the image we were looking at * back to not log scale, * so we don't log it twice if its reselected. */ if (log.isSelected()) toggleLogarithmic(false); currentlyDisplayedImageResult = img; pixelAccessCursor = img.randomAccess(); // Currently disabled, due to lag of non-histograms :-) // disable list and copy button if it is no histogram result listButton.setEnabled( isHistogram(img) ); copyButton.setEnabled( isHistogram(img) ); drawImage(img); toggleLogarithmic(log.isSelected()); // ensure a valid layout, we changed the image getContentPane().validate(); getContentPane().repaint(); } + @Override public void itemStateChanged(ItemEvent e) { if (e.getStateChange() == ItemEvent.SELECTED) { RandomAccessibleInterval> img = ((NamedContainer> >)(e.getItem())).getObject(); adjustDisplayedImage(img); } } protected void toggleLogarithmic(boolean enabled){ if (imp == null) return; ImageProcessor ip = imp.getProcessor(); if (enabled) { ip.snapshot(); ip.log(); ip.resetMinAndMax(); } else ip.reset(); imagePanel.repaint(); } + @Override public void actionPerformed(ActionEvent e) { if (e.getSource() == log) { toggleLogarithmic(log.isSelected()); } } + @Override public void lostOwnership(Clipboard clipboard, Transferable contents) { // nothing to do here } } diff --git a/src/test/java/tests/MandersColocalizationTest.java b/src/test/java/tests/MandersColocalizationTest.java index 6c97778..7d31a5a 100644 --- a/src/test/java/tests/MandersColocalizationTest.java +++ b/src/test/java/tests/MandersColocalizationTest.java @@ -1,198 +1,198 @@ package tests; import static org.junit.Assert.assertEquals; import algorithms.MandersColocalization; import algorithms.MandersColocalization.MandersResults; import algorithms.MissingPreconditionException; import net.imglib2.Cursor; import net.imglib2.TwinCursor; import net.imglib2.converter.Converter; import net.imglib2.converter.Converters; import net.imglib2.type.logic.BitType; import net.imglib2.type.numeric.integer.UnsignedByteType; import net.imglib2.view.Views; import gadgets.ThresholdMode; import org.junit.Test; /** * This class contains JUnit 4 test cases for the calculation * of Manders' split colocalization coefficients * * @author Dan White & Tom Kazimiers */ public class MandersColocalizationTest extends ColocalisationTest { /** * This method tests artificial test images as detailed in * the Manders et al. paper, using above zero threshold (none). * Note: It is not sensitive to choosing the wrong channel to test for * above threshold, because threshold is same for both channels: above zero, * and also that the blobs overlap perfectly or not at all. */ @Test public void mandersPaperImagesTest() throws MissingPreconditionException { MandersColocalization mc = new MandersColocalization(); TwinCursor cursor; MandersResults r; // test A-A combination cursor = new TwinCursor( mandersA.randomAccess(), mandersA.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get().createVariable()); assertEquals(1.0d, r.m1, 0.0001); assertEquals(1.0d, r.m2, 0.0001); // test A-B combination cursor = new TwinCursor( mandersA.randomAccess(), mandersB.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.75d, r.m1, 0.0001); assertEquals(0.75d, r.m2, 0.0001); // test A-C combination cursor = new TwinCursor( mandersA.randomAccess(), mandersC.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.5d, r.m1, 0.0001); assertEquals(0.5d, r.m2, 0.0001); // test A-D combination cursor = new TwinCursor( mandersA.randomAccess(), mandersD.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.25d, r.m1, 0.0001); assertEquals(0.25d, r.m2, 0.0001); // test A-E combination cursor = new TwinCursor( mandersA.randomAccess(), mandersE.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.0d, r.m1, 0.0001); assertEquals(0.0d, r.m2, 0.0001); // test A-F combination cursor = new TwinCursor( mandersA.randomAccess(), mandersF.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.25d, r.m1, 0.0001); assertEquals(0.3333d, r.m2, 0.0001); // test A-G combination.firstElement( cursor = new TwinCursor( mandersA.randomAccess(), mandersG.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.25d, r.m1, 0.0001); assertEquals(0.50d, r.m2, 0.0001); // test A-H combination cursor = new TwinCursor( mandersA.randomAccess(), mandersH.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.25d, r.m1, 0.0001); assertEquals(1.00d, r.m2, 0.0001); // test A-I combination cursor = new TwinCursor( mandersA.randomAccess(), mandersI.randomAccess(), Views.iterable(mandersAlwaysTrueMask).localizingCursor()); r = mc.calculateMandersCorrelation(cursor, mandersA.randomAccess().get()); assertEquals(0.083d, r.m1, 0.001); assertEquals(0.75d, r.m2, 0.0001); } /** * This method tests real experimental noisy but * biologically perfectly colocalized test images, * using previously calculated autothresholds (.above mode) * Amongst other things, hopefully it is sensitive to * choosing the wrong channel to test for above threshold */ @Test public void mandersRealNoisyImagesTest() throws MissingPreconditionException { MandersColocalization mrnc = new MandersColocalization(); // test biologically perfect but noisy image coloc combination // this cast is bad, so use Views.iterable instead. //Cursor mask = Converters.convert((IterableInterval) positiveCorrelationMaskImage, Cursor mask = Converters.convert(Views.iterable(positiveCorrelationMaskImage), new Converter() { @Override public void convert(UnsignedByteType arg0, BitType arg1) { arg1.set(arg0.get() > 0); } }, new BitType()).cursor(); TwinCursor twinCursor; MandersResults r; // Manually set the thresholds for ch1 and ch2 with the results from a // Costes Autothreshold using bisection implementation of regression, of the images used UnsignedByteType thresholdCh1 = new UnsignedByteType(); thresholdCh1.setInteger(70); UnsignedByteType thresholdCh2 = new UnsignedByteType(); thresholdCh2.setInteger(53); //Set the threshold mode ThresholdMode tMode; tMode = ThresholdMode.Above; // Set the TwinCursor to have the mask image channel, and 2 images. twinCursor = new TwinCursor( positiveCorrelationImageCh1.randomAccess(), positiveCorrelationImageCh2.randomAccess(), mask); // Use the constructor that takes ch1 and ch2 autothresholds and threshold mode. r = mrnc.calculateMandersCorrelation(twinCursor, thresholdCh1, thresholdCh2, tMode); assertEquals(0.705665d, r.m1, 0.000001); assertEquals(0.724752d, r.m2, 0.000001); } -} \ No newline at end of file +}