diff --git a/src/main/java/Coloc_2.java b/src/main/java/Coloc_2.java index 0309442..d5a4879 100644 --- a/src/main/java/Coloc_2.java +++ b/src/main/java/Coloc_2.java @@ -1,717 +1,718 @@ import algorithms.Algorithm; import algorithms.AutoThresholdRegression; 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 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.util.Util; import net.imglib2.view.Views; import results.PDFWriter; import results.ResultHandler; import results.SingleWindowDisplay; import results.Warning; /** Copyright 2010, 2011 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/ . */ /** * A plugin which does analysis colocalisation on a pair of images. * * @param */ 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(); // default indices of image, mask and ROI choices protected static int index1 = 0; protected static int index2 = 1; protected static int indexMask = 0; protected static int indexRoi = 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; 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] = ""; } } // set up the users preferences displayImages = Prefs.get(PREF_KEY+"displayImages", false); autoSavePdf = Prefs.get(PREF_KEY+"autoSavePdf", true); boolean displayShuffledCostes = Prefs.get(PREF_KEY+"displayShuffledCostes", false); boolean useLiCh1 = Prefs.get(PREF_KEY+"useLiCh1", true); boolean useLiCh2 = Prefs.get(PREF_KEY+"useLiCh2", true); boolean useLiICQ = Prefs.get(PREF_KEY+"useLiICQ", true); boolean useSpearmanRank = Prefs.get(PREF_KEY+"useSpearmanRank", true); boolean useManders = Prefs.get(PREF_KEY+"useManders", true); boolean useKendallTau = Prefs.get(PREF_KEY+"useKendallTau", true); boolean useScatterplot = Prefs.get(PREF_KEY+"useScatterplot", true); boolean useCostes = Prefs.get(PREF_KEY+"useCostes", true); int psf = (int) Prefs.get(PREF_KEY+"psf", 3); int nrCostesRandomisations = (int) Prefs.get(PREF_KEY+"nrCostesRandomisations", 10); /* make sure the default indices are no bigger * than the amount of images we have */ index1 = clip( index1, 0, titles.length ); index2 = clip( index2, 0, titles.length ); indexMask = clip( indexMask, 0, roisAndMasks.length - 1); gd.addChoice("Channel_1", titles, titles[index1]); gd.addChoice("Channel_2", titles, titles[index2]); gd.addChoice("ROI_or_mask", roisAndMasks, roisAndMasks[indexMask]); //gd.addChoice("Use ROI", roiLabels, roiLabels[indexRoi]); gd.addCheckbox("Show_\"Save_PDF\"_Dialog", autoSavePdf); gd.addCheckbox("Display_Images_in_Result", displayImages); gd.addCheckbox("Display_Shuffled_Images", displayShuffledCostes); final Checkbox shuffleCb = (Checkbox) gd.getCheckboxes().lastElement(); // Add algorithm options gd.addMessage("Algorithms:"); gd.addCheckbox("Li_Histogram_Channel_1", useLiCh1); gd.addCheckbox("Li_Histogram_Channel_2", useLiCh2); gd.addCheckbox("Li_ICQ", useLiICQ); gd.addCheckbox("Spearman's_Rank_Correlation", useSpearmanRank); gd.addCheckbox("Manders'_Correlation", useManders); gd.addCheckbox("Kendall's_Tau_Rank_Correlation", useKendallTau); gd.addCheckbox("2D_Instensity_Histogram", useScatterplot); gd.addCheckbox("Costes'_Significance_Test", useCostes); final Checkbox costesCb = (Checkbox) gd.getCheckboxes().lastElement(); gd.addNumericField("PSF", psf, 1); gd.addNumericField("Costes_randomisations", nrCostesRandomisations, 0); // disable shuffle checkbox if costes checkbox is set to "off" shuffleCb.setEnabled(useCostes); costesCb.addItemListener(new ItemListener() { @Override public void itemStateChanged(ItemEvent e) { shuffleCb.setEnabled(costesCb.getState()); } }); // show the dialog, finally gd.showDialog(); // do nothing if dialog has been canceled if (gd.wasCanceled()) return false; ImagePlus imp1 = WindowManager.getImage(gd.getNextChoiceIndex() + 1); ImagePlus imp2 = WindowManager.getImage(gd.getNextChoiceIndex() + 1); // get image names for output Ch1Name = imp1.getTitle(); Ch2Name = imp2.getTitle(); // make sure both images have the same bit-depth if (imp1.getBitDepth() != imp2.getBitDepth()) { IJ.showMessage("Both images must have the same bit-depth."); return false; } // get information about the mask/ROI to use indexMask = gd.getNextChoiceIndex(); if (indexMask == 0) roiConfig = RoiConfiguration.None; else if (indexMask == 1) roiConfig = RoiConfiguration.Img1; else if (indexMask == 2) roiConfig = RoiConfiguration.Img2; else if (indexMask == 3) roiConfig = RoiConfiguration.RoiManager; else { roiConfig = RoiConfiguration.Mask; /* Make indexMask the reference to the mask image to use. * To do this we reduce it by three for the first three * entries in the combo box. */ indexMask = indexMask - 4; } // save the ImgLib wrapped images as members img1 = ImagePlusAdapter.wrap(imp1); img2 = ImagePlusAdapter.wrap(imp2); /* check if we have a valid ROI for the selected configuration * and if so, get the ROI's bounds. Alternatively, a mask can * be selected (that is basically all, but a rectangle). */ if (roiConfig == RoiConfiguration.Img1 && hasValidRoi(imp1)) { createMasksFromImage(imp1); } else if (roiConfig == RoiConfiguration.Img2 && hasValidRoi(imp2)) { createMasksFromImage(imp2); } else if (roiConfig == RoiConfiguration.RoiManager) { if (!createMasksFromRoiManager(imp1.getWidth(), imp1.getHeight())) return false; } else if (roiConfig == RoiConfiguration.Mask) { // get the image to be used as mask ImagePlus maskImp = WindowManager.getImage(windowList[indexMask]); Img 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)); } // 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+"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)) ); // 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() ); + 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 }; } } diff --git a/src/main/java/algorithms/AutoThresholdRegression.java b/src/main/java/algorithms/AutoThresholdRegression.java index 9d121d0..e4fe9c4 100644 --- a/src/main/java/algorithms/AutoThresholdRegression.java +++ b/src/main/java/algorithms/AutoThresholdRegression.java @@ -1,281 +1,282 @@ 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.util.Util; 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 { /* the threshold for y-intercept to y-max to * raise a warning about it being to high. */ final double warnYInterceptToYMaxRatioThreshold = 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 bToYMaxRatio = 0.0; //This is the Pearson's correlation we will use for further calculations PearsonsCorrelation pearsonsCorrellation; public AutoThresholdRegression(PearsonsCorrelation pc) { super("auto threshold regression"); pearsonsCorrellation = pc; } @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(); - T type1 = cursor.getFirst(); - double ch1 = type1.getRealDouble(); - T type2 = cursor.getSecond(); - double ch2 = type2.getRealDouble(); + 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) ); double m = num/denom; double b = ch2Mean - m*ch1Mean ; // initialize some variables relevant for regression // indicates whether the threshold has been found or not boolean thresholdFound = false; // the maximum number of iterations to look for the threshold final int maxIterations = 100; // the current iteration int iteration = 0; // the initial thresholds double threshold1 = (container.getMaxCh1() + container.getMinCh1()) * 0.5; double threshold2 = container.getMaxCh1(); // Min threshold not yet implemented double ch1ThreshMax = container.getMaxCh1(); double ch2ThreshMax = container.getMaxCh2(); // define some image type specific threshold variables - T thresholdCh1 = img1.randomAccess().get(); - T thresholdCh2 = img2.randomAccess().get(); + 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 = img1.randomAccess().get(); + T dummyT = type.createVariable(); final double minVal = dummyT.getMinValue(); final double maxVal = dummyT.getMaxValue(); // do regression while (!thresholdFound && iteration 0) { // as long as r > 0 we go half the way down threshold1 = threshold1 - thrDiff * 0.5; } // reset the cursor to reuse it cursor.reset(); // increment iteration counter iteration++; } /* 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 = img1.randomAccess().get(); + ch1MinThreshold = type.createVariable(); ch1MinThreshold.setReal(minVal); - ch1MaxThreshold = img1.randomAccess().get(); + ch1MaxThreshold = type.createVariable(); ch1MaxThreshold.setReal(clamp(ch1ThreshMax, minVal, maxVal)); - ch2MinThreshold = img2.randomAccess().get(); + ch2MinThreshold = type.createVariable(); ch2MinThreshold.setReal(minVal); - ch2MaxThreshold = img2.randomAccess().get(); + ch2MaxThreshold = type.createVariable(); ch2MaxThreshold.setReal(clamp(ch2ThreshMax, minVal, maxVal)); autoThresholdSlope = m; autoThresholdIntercept = b; bToYMaxRatio = b / container.getMaxCh2(); // add warnings if values are not in tolerance range if ( Math.abs(bToYMaxRatio) > warnYInterceptToYMaxRatioThreshold ) { addWarning("y-intercept high", "The absolute y-intercept of the auto threshold regression line is high. Maybe you should use a ROI, maybe do a background subtraction in both channels."); } // 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() ) { addWarning("thresholds too low", "The auto threshold method could not find a positive threshold."); } } /** * 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; } 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-max ratio", bToYMaxRatio, 2 ); handler.handleValue( "Ch1 Max Threshold", ch1MaxThreshold.getRealDouble(), 2); handler.handleValue( "Ch2 Max Threshold", ch2MaxThreshold.getRealDouble(), 2); } public double getBToYMaxRatio() { return bToYMaxRatio; } public double getWarnYInterceptToYMaxRatioThreshold() { return warnYInterceptToYMaxRatioThreshold; } 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 6bc8a95..3441974 100644 --- a/src/main/java/algorithms/CostesSignificanceTest.java +++ b/src/main/java/algorithms/CostesSignificanceTest.java @@ -1,401 +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.util.Util; 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(); + 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() ); + 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; } 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/test/java/tests/MandersColocalizationTest.java b/src/test/java/tests/MandersColocalizationTest.java index 64048a1..e48fe66 100644 --- a/src/test/java/tests/MandersColocalizationTest.java +++ b/src/test/java/tests/MandersColocalizationTest.java @@ -1,131 +1,131 @@ package tests; import static org.junit.Assert.assertEquals; import algorithms.MandersColocalization; import algorithms.MandersColocalization.MandersResults; import algorithms.MissingPreconditionException; import net.imglib2.TwinCursor; import net.imglib2.type.numeric.integer.UnsignedByteType; import net.imglib2.util.Util; import net.imglib2.view.Views; import org.junit.Test; public class MandersColocalizationTest extends ColocalisationTest { @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()); + 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); } } \ No newline at end of file diff --git a/src/test/java/tests/TestImageAccessor.java b/src/test/java/tests/TestImageAccessor.java index e8ed819..cf318dd 100644 --- a/src/test/java/tests/TestImageAccessor.java +++ b/src/test/java/tests/TestImageAccessor.java @@ -1,399 +1,399 @@ package tests; import static org.junit.Assume.assumeNotNull; import algorithms.MissingPreconditionException; import gadgets.MaskFactory; import ij.ImagePlus; import ij.gui.NewImage; import ij.gui.Roi; import ij.io.Opener; import ij.process.ImageProcessor; import java.awt.Color; import java.io.BufferedInputStream; import java.io.InputStream; import java.util.Arrays; import net.imglib2.Cursor; import net.imglib2.Interval; import net.imglib2.Localizable; import net.imglib2.Point; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessible; import net.imglib2.RandomAccessibleInterval; import net.imglib2.TwinCursor; import net.imglib2.algorithm.gauss.Gauss; import net.imglib2.algorithm.math.ImageStatistics; import net.imglib2.img.ImagePlusAdapter; 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.type.numeric.real.FloatType; import net.imglib2.util.Util; import net.imglib2.view.Views; /** * A class containing some testing helper methods. It allows * to open Tiffs from within the Jar file and can generate noise * images. * * @author Dan White & Tom Kazimiers */ public class TestImageAccessor { /* a static opener for opening images without the * need for creating every time a new opener */ static Opener opener = new Opener(); /** * Loads a Tiff file from within the jar. The given path is treated * as relative to this tests-package (i.e. "Data/test.tiff" refers * to the test.tiff in sub-folder Data). * * @param The wanted output type. * @param relPath The relative path to the Tiff file. * @return The file as ImgLib image. */ public static & NativeType> RandomAccessibleInterval loadTiffFromJar(String relPath) { InputStream is = TestImageAccessor.class.getResourceAsStream(relPath); BufferedInputStream bis = new BufferedInputStream(is); ImagePlus imp = opener.openTiff(bis, "The Test Image"); assumeNotNull(imp); return ImagePlusAdapter.wrap(imp); } /** * Creates a noisy image that is created by repeatedly adding points * with random intensity to the canvas. That way it tries to mimic the * way a microscope produces images. This convenience method uses the * default values of a point size of 3.0 and produces 5000 points. * After the creation the image is smoothed with a sigma of one in each * direction. * * @param The wanted output type. * @param width The image width. * @param height The image height. * @return The noise image. */ public static & NativeType> RandomAccessibleInterval produceNoiseImageSmoothed(T type, int width, int height) { return produceNoiseImageSmoothed(type, width, height, 3.0f, 5000, new double[] {1.0,1.0}); } /** * Creates a noisy image that is created by repeatedly adding points * with random intensity to the canvas. That way it tries to mimic the * way a microscope produces images. * * @param The wanted output type. * @param width The image width. * @param height The image height. * @param dotSize The size of the dots. * @param numDots The number of dots. * @param smoothingSigma The two dimensional sigma for smoothing. * @return The noise image. */ public static & NativeType> RandomAccessibleInterval produceNoiseImage(int width, int height, float dotSize, int numDots) { /* For now (probably until ImageJ2 is out) we use an * ImageJ image to draw circles. */ int options = NewImage.FILL_BLACK + NewImage.CHECK_AVAILABLE_MEMORY; ImagePlus img = NewImage.createByteImage("Noise", width, height, 1, options); ImageProcessor imp = img.getProcessor(); float dotRadius = dotSize * 0.5f; int dotIntSize = (int) dotSize; for (int i=0; i < numDots; i++) { int x = (int) (Math.random() * width - dotRadius); int y = (int) (Math.random() * height - dotRadius); imp.setColor(Color.WHITE); imp.fillOval(x, y, dotIntSize, dotIntSize); } // we changed the data, so update it img.updateImage(); // create the new image RandomAccessibleInterval noiseImage = ImagePlusAdapter.wrap(img); return noiseImage; } public static & NativeType> RandomAccessibleInterval produceNoiseImageSmoothed(T type, int width, int height, float dotSize, int numDots, double[] smoothingSigma) { RandomAccessibleInterval noiseImage = produceNoiseImage(width, height, dotSize, numDots); return gaussianSmooth(noiseImage, smoothingSigma); } /** * This method creates a noise image that has a specified mean. * Every pixel has a value uniformly distributed around mean with * the maximum spread specified. * * @return a new noise image * @throws MissingPreconditionException if specified means and spreads are not valid */ public static & NativeType> RandomAccessibleInterval produceMeanBasedNoiseImage(T type, int width, int height, double mean, double spread, double[] smoothingSigma) throws MissingPreconditionException { if (mean < spread || (mean + spread) > type.getMaxValue()) { throw new MissingPreconditionException("Mean must be larger than spread, and mean plus spread must be smaller than max of the type"); } // create the new image ImgFactory imgFactory = new ArrayImgFactory(); RandomAccessibleInterval noiseImage = imgFactory.create( new int[] {width, height}, type); // "Noise image"); for (T value : Views.iterable(noiseImage)) { value.setReal( mean + ( (Math.random() - 0.5) * spread ) ); } return gaussianSmooth(noiseImage, smoothingSigma); } /** * This method creates a noise image that is made of many little * sticks oriented in a random direction. How many of them and * what the length of them are can be specified. * * @return a new noise image that is not smoothed */ public static & NativeType> RandomAccessibleInterval produceSticksNoiseImage(int width, int height, int numSticks, int lineWidth, double maxLength) { /* For now (probably until ImageJ2 is out) we use an * ImageJ image to draw lines. */ int options = NewImage.FILL_BLACK + NewImage.CHECK_AVAILABLE_MEMORY; ImagePlus img = NewImage.createByteImage("Noise", width, height, 1, options); ImageProcessor imp = img.getProcessor(); imp.setColor(Color.WHITE); imp.setLineWidth(lineWidth); for (int i=0; i < numSticks; i++) { // find random starting point int x = (int) (Math.random() * width); int y = (int) (Math.random() * height); // create random stick length and direction double length = Math.random() * maxLength; double angle = Math.random() * 2 * Math.PI; // calculate random point on circle, for the direction int destX = x + (int) (length * Math.cos(angle)); int destY = y + (int) (length * Math.sin(angle)); // now draw the line imp.drawLine(x, y, destX, destY); } // we changed the data, so update it img.updateImage(); return ImagePlusAdapter.wrap(img); } /** * This method creates a smoothed noise image that is made of * many little sticks oriented in a random direction. How many * of them and what the length of them are can be specified. * * @return a new noise image that is smoothed */ public static & NativeType> RandomAccessibleInterval produceSticksNoiseImageSmoothed(T type, int width, int height, int numSticks, int lineWidth, double maxLength, double[] smoothingSigma) { RandomAccessibleInterval noiseImage = produceSticksNoiseImage(width, height, numSticks, lineWidth, maxLength); return gaussianSmooth(noiseImage, smoothingSigma); } /** * Generates a Perlin noise image. It is based on Ken Perlin's * reference implementation (ImprovedNoise class) and a small * bit of Kas Thomas' sample code (http://asserttrue.blogspot.com/). */ public static & NativeType> RandomAccessibleInterval producePerlinNoiseImage(T type, int width, int height, double z, double scale) { // create the new image ImgFactory imgFactory = new ArrayImgFactory(); RandomAccessibleInterval noiseImage = imgFactory.create( new int[] {width, height}, type); Cursor noiseCursor = Views.iterable(noiseImage).localizingCursor(); double xOffset = Math.random() * (width*width); double yOffset = Math.random() * (height*height); while (noiseCursor.hasNext()) { noiseCursor.fwd(); double x = (noiseCursor.getDoublePosition(0) + xOffset) * scale; double y = (noiseCursor.getDoublePosition(1) + yOffset) * scale; float t = (float)ImprovedNoise.noise( x, y, z); // ImprovedNoise.noise returns a float in the range [-1..1], // whereas we want a float in the range [0..1], so: t = (1 + t) * 0.5f; noiseCursor.get().setReal(t); } //return gaussianSmooth(noiseImage, imgFactory, smoothingSigma); return noiseImage; } /** * Gaussian Smooth of the input image using intermediate float format. * @param * @param img * @param sigma * @return */ public static & NativeType> RandomAccessibleInterval gaussianSmooth( RandomAccessibleInterval img, double[] sigma) { Interval interval = Views.iterable(img); ImgFactory outputFactory = new ArrayImgFactory(); final long[] dim = new long[ img.numDimensions() ]; img.dimensions(dim); RandomAccessibleInterval output = outputFactory.create( dim, - img.randomAccess().get() ); + img.randomAccess().get().createVariable() ); final long[] pos = new long[ img.numDimensions() ]; Arrays.fill(pos, 0); Localizable origin = new Point(pos); ImgFactory tempFactory = new ArrayImgFactory(); RandomAccessible input = Views.extendMirrorSingle(img); Gauss.inFloat(sigma, input, interval, output, origin, tempFactory); return output; } /** * Inverts an image. * * @param The images data type. * @param image The image to convert. * @return The inverted image. */ public static & NativeType> RandomAccessibleInterval invertImage( RandomAccessibleInterval image) { Cursor imgCursor = Views.iterable(image).localizingCursor(); // invert the image long[] dim = new long[ image.numDimensions() ]; image.dimensions(dim); ArrayImgFactory imgFactory = new ArrayImgFactory(); RandomAccessibleInterval invImg = imgFactory.create( - dim, image.randomAccess().get() ); // "Inverted " + image.getName()); + dim, image.randomAccess().get().createVariable() ); // "Inverted " + image.getName()); RandomAccess invCursor = invImg.randomAccess(); while (imgCursor.hasNext()) { imgCursor.fwd(); invCursor.setPosition(imgCursor); invCursor.get().setReal( imgCursor.get().getMaxValue() - imgCursor.get().getRealDouble() ); } return invImg; } /** * Converts an arbitrary image to a black/white version of it. * All image data lower or equal 0.5 times the maximum value * of the image type will get black, the rest will turn white. */ public static & NativeType> RandomAccessibleInterval makeBinaryImage( RandomAccessibleInterval image) { T binSplitValue = image.randomAccess().get(); binSplitValue.setReal( binSplitValue.getMaxValue() * 0.5 ); return TestImageAccessor.makeBinaryImage(image, binSplitValue); } /** * Converts an arbitrary image to a black/white version of it. * All image data lower or equal the splitValue will get black, * the rest will turn white. */ public static & NativeType> RandomAccessibleInterval makeBinaryImage( RandomAccessibleInterval image, T splitValue) { Cursor imgCursor = Views.iterable(image).localizingCursor(); // make a new image of the same type, but binary long[] dim = new long[ image.numDimensions() ]; image.dimensions(dim); ArrayImgFactory imgFactory = new ArrayImgFactory(); RandomAccessibleInterval binImg = imgFactory.create( dim, - image.randomAccess().get() ); // "Binary image of " + image.getName()); + image.randomAccess().get().createVariable() ); // "Binary image of " + image.getName()); RandomAccess invCursor = binImg.randomAccess(); while (imgCursor.hasNext()) { imgCursor.fwd(); invCursor.setPosition(imgCursor); T currentValue = invCursor.get(); if (currentValue.compareTo(splitValue) > 0) currentValue.setReal( currentValue.getMaxValue() ); else currentValue.setZero(); } return binImg; } /** * A method to combine a foreground image and a background image. * If data on the foreground image is above zero, it will be * placed on the background. While doing that, the image data from * the foreground is scaled to be in range of the background. */ public static > void combineImages(RandomAccessibleInterval background, RandomAccessibleInterval foreground) { final long[] dim = new long[ background.numDimensions() ]; background.dimensions(dim); RandomAccessibleInterval alwaysTrueMask = MaskFactory.createMask(dim, true); TwinCursor cursor = new TwinCursor( background.randomAccess(), foreground.randomAccess(), Views.iterable(alwaysTrueMask).localizingCursor()); // find a scaling factor for scale forground range into background double bgMin = ImageStatistics.getImageMin(background).getRealDouble(); double bgMax = ImageStatistics.getImageMax(background).getRealDouble(); double fgMin = ImageStatistics.getImageMin(foreground).getRealDouble(); double fgMax = ImageStatistics.getImageMax(foreground).getRealDouble(); double scaling = (bgMax - bgMin ) / (fgMax - fgMin); // iterate over both images while (cursor.hasNext()) { cursor.fwd(); T bgData = cursor.getFirst(); double fgData = cursor.getSecond().getRealDouble() * scaling; if (fgData > 0.01) { /* if the foreground data is above zero, copy * it to the background. */ bgData.setReal(fgData); } } } /** * Creates a mask image with a black background and a white * rectangular foreground. * * @param width The width of the result image. * @param height The height of the result image. * @param offset The offset of the rectangular mask. * @param size The size of the rectangular mask. * @return A black image with a white rectangle on it. */ public static & NativeType> RandomAccessibleInterval createRectengularMaskImage( long width, long height, long[] offset, long[] size) { /* For now (probably until ImageJ2 is out) we use an * ImageJ image to draw lines. */ int options = NewImage.FILL_BLACK + NewImage.CHECK_AVAILABLE_MEMORY; ImagePlus img = NewImage.createByteImage("Noise", (int)width, (int)height, 1, options); ImageProcessor imp = img.getProcessor(); imp.setColor(Color.WHITE); Roi rect = new Roi(offset[0], offset[1], size[0], size[1]); imp.fill(rect); // we changed the data, so update it img.updateImage(); return ImagePlusAdapter.wrap(img); } }