wrap(maskImp);
- // get a valid mask info for the image
- final 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 = gdIndexRegr;
-
- // read out GUI data
- autoSavePdf = gdAutoSavePdf;
- displayImages = gdDisplayImages;
-
- // Parse algorithm options
- pearsonsCorrelation = new PearsonsCorrelation<>(
- PearsonsCorrelation.Implementation.Fast);
-
- if (gdUseLiCh1) liHistogramCh1 = new LiHistogram2D<>("Li - Ch1", true);
- if (gdUseLiCh2) liHistogramCh2 = new LiHistogram2D<>("Li - Ch2", false);
- if (gdUseLiICQ) liICQ = new LiICQ<>();
- if (gdUseSpearmanRank) {
- SpearmanRankCorrelation = new SpearmanRankCorrelation<>();
- }
- if (gdUseManders) mandersCorrelation = new MandersColocalization<>();
- if (gdUseKendallTau) kendallTau = new KendallTauRankCorrelation<>();
- if (gdUseScatterplot) histogram2D = new Histogram2D<>(
- "2D intensity histogram");
- if (gdUseCostes) {
- costesSignificance = new CostesSignificanceTest<>(pearsonsCorrelation,
- gdPsf, gdNrCostesRandomisations, gdDisplayShuffledCostes);
- }
-
- 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.
- *
- * NB: This method returns {@code void} for binary backwards compatibility
- * with old code which might invoke it. If you want access to the
- * {@link AnalysisResults} directly, call
- * {@link #colocalise(Img, Img, Coloc_2.BoundingBox, Img, List)} with
- * {@code null} for {@code extraHandlers}.
- *
- *
- * @param image1
- * @param image2
- * @param roi
- * @param mask
- * @throws MissingPreconditionException
- */
- public void colocalise(final Img image1, final Img image2,
- final BoundingBox roi, final Img mask)
- throws MissingPreconditionException
- {
- colocalise(image1, image2, roi, mask, null);
- }
-
- /**
- * 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 image1 First image.
- * @param image2 Second image.
- * @param roi Region of interest to which analysis is confined.
- * @param mask Mask to which analysis is confined.
- * @param extraHandlers additional objects to be notified of analysis results.
- * @return Data structure housing the results.
- * @throws MissingPreconditionException
- */
- public AnalysisResults colocalise(final Img image1, final Img image2,
- final BoundingBox roi, final Img mask,
- final List> extraHandlers)
- throws MissingPreconditionException
- {
- // create a new container for the selected images and channels
- DataContainer container;
- if (mask != null) {
- container = new DataContainer<>(image1, image2, 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<>(image1, image2, img1Channel, img2Channel,
- Ch1Name, Ch2Name, roi.offset, roi.size);
- }
- else {
- // no mask and no ROI is present
- container = new DataContainer<>(image1, image2, img1Channel, img2Channel,
- Ch1Name, Ch2Name);
- }
-
- // create a results handler
- final List> listOfResultHandlers =
- new ArrayList<>();
- final AnalysisResults analysisResults = new AnalysisResults<>();
- listOfResultHandlers.add(analysisResults);
- final PDFWriter pdfWriter = new PDFWriter<>(container);
- final boolean headless = Boolean.getBoolean("java.awt.headless");
- final SingleWindowDisplay swDisplay;
- if (headless) swDisplay = null;
- else {
- swDisplay = new SingleWindowDisplay<>(container, pdfWriter);
- listOfResultHandlers.add(swDisplay);
- }
- listOfResultHandlers.add(pdfWriter);
- if (extraHandlers != null) listOfResultHandlers.addAll(extraHandlers);
- // ResultHandler resultHandler = new EasyDisplay(container);
-
- // this contains the algorithms that will be run when the user clicks ok
- final 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;
- final int jobs = userSelectedJobs.size();
- for (final Algorithm a : userSelectedJobs) {
- try {
- count++;
- IJ.showStatus(count + "/" + jobs + ": Running " + a.getName());
- a.execute(container);
- }
- catch (final MissingPreconditionException e) {
- for (final 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 (final Algorithm a : userSelectedJobs) {
- for (final 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) {
- final long[] offset = container.getMaskBBOffset();
- final 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();
- }
- channel1 = project(channel1);
- channel2 = project(channel2);
-
- for (final ResultHandler r : listOfResultHandlers) {
- r.handleImage(channel1, "Channel 1 (Max Projection)");
- r.handleImage(channel2, "Channel 2 (Max Projection)");
- }
- }
- if (swDisplay != null) {
- // do the actual results processing
- swDisplay.process();
- // add window to the IJ window manager
- swDisplay.addWindowListener(new WindowAdapter() {
-
- @Override
- public void windowClosing(final WindowEvent e) {
- WindowManager.removeWindow(swDisplay);
- swDisplay.dispose();
- // NB: For some reason, garbage collection of this bundle of objects
- // does not occur when this window listener reference remains in
- // place. As such, we explicitly unregister ourself here.
- swDisplay.removeWindowListener(this);
- }
- });
- WindowManager.addWindow(swDisplay);
- }
-
- // show PDF saving dialog if requested
- if (autoSavePdf) pdfWriter.process();
-
- return analysisResults;
- }
-
- private RandomAccessibleInterval project(
- final RandomAccessibleInterval image)
- {
- if (image.numDimensions() < 2) {
- throw new IllegalArgumentException("Dimensionality too small: " + //
- image.numDimensions());
- }
-
- final IterableInterval input = Views.iterable(image);
- final T type = input.firstElement(); // e.g. unsigned 8-bit
- final long xLen = image.dimension(0);
- final long yLen = image.dimension(1);
-
- // initialize output image with minimum value of the pixel type
- final long[] outputDims = { xLen, yLen };
- final Img output = new ArrayImgFactory().create(outputDims, type);
- for (final T sample : output) {
- sample.setReal(type.getMinValue());
- }
-
- // loop over the input image, performing the max projection
- final Cursor inPos = input.localizingCursor();
- final RandomAccess outPos = output.randomAccess();
- while (inPos.hasNext()) {
- final T inPix = inPos.next();
- final long xPos = inPos.getLongPosition(0);
- final long yPos = inPos.getLongPosition(1);
- outPos.setPosition(xPos, 0);
- outPos.setPosition(yPos, 1);
- final T outPix = outPos.get();
- if (outPix.compareTo(inPix) < 0) {
- outPix.set(inPix);
- }
- }
- return output;
- }
-
- /**
- * 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(final Img mask) {
- final Cursor cursor = mask.localizingCursor();
-
- final int numMaskDims = mask.numDimensions();
- // the "off type" of the mask
- final T offType = mask.firstElement().createVariable();
- offType.setZero();
- // the corners of the bounding box
- final long[][] minMax = new long[2][];
- // indicates if mask data has been found
- boolean maskFound = false;
- // a container for temporary position information
- final long[] pos = new long[numMaskDims];
- // walk over the mask
- while (cursor.hasNext()) {
- cursor.fwd();
- final 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
- minMax[0] = Arrays.copyOf(pos, numMaskDims);
- minMax[1] = 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 < numMaskDims; d++) {
- if (pos[d] < minMax[0][d]) {
- // is it smaller than min
- minMax[0][d] = pos[d];
- }
- else if (pos[d] > minMax[1][d]) {
- // is it larger than max
- minMax[1][d] = pos[d];
- }
- }
- }
- }
- }
- if (!maskFound) return null;
-
- // calculate size
- final long[] size = new long[numMaskDims];
- for (int d = 0; d < numMaskDims; d++)
- size[d] = minMax[1][d] - minMax[0][d] + 1;
- // create and add bounding box
- final BoundingBox bb = new BoundingBox(minMax[0], size);
- return new MaskInfo(bb, mask);
- }
-
- /**
- * Adds the provided Algorithm to the list if it is not null.
- */
- protected void addIfValid(final Algorithm a,
- final 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(final ImagePlus imp) {
- final Roi roi = imp.getRoi();
- if (roi == null) return false;
-
- final 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(final int val, final int min, final 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(final ImagePlus imp) {
- // get ROIs from current image in Fiji
- final 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(final int width,
- final int height)
- {
- final RoiManager roiManager = RoiManager.getInstance();
- if (roiManager == null) {
- IJ.error("Could not get ROI Manager instance.");
- return false;
- }
- final 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 {@code 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(final Roi[] rois, final int width,
- final int height)
- {
- // create empty list
- masks.clear();
-
- for (final Roi r : rois) {
- final MaskInfo mi = new MaskInfo();
- // add it to the list of masks/ROIs
- masks.add(mi);
- // get the ROIs/masks bounding box
- final Rectangle rect = r.getBounds();
- mi.roi = new BoundingBox(new long[] { rect.x, rect.y }, new long[] {
- rect.width, rect.height });
- final 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
- final 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
- final 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(
- final RandomAccessibleInterval image,
- final RandomAccessibleInterval mask, final long[] offset,
- final long[] size) throws MissingPreconditionException
- {
- final 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
- final TwinCursor cursor = new TwinCursor<>(image.randomAccess(), //
- image.randomAccess(), Views.iterable(mask).localizingCursor());
- // prepare output image
- final ImgFactory maskFactory = new ArrayImgFactory<>();
- // Img maskImage = maskFactory.create( size, name );
- final RandomAccessibleInterval maskImage = maskFactory.create(size, //
- image.randomAccess().get().createVariable());
- final 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(final Roi roi) {
- if (roi instanceof ShapeRoi) return ((ShapeRoi) roi).getRois();
- return new Roi[] { roi };
- }
-}
diff --git a/src/main/java/sc/fiji/coloc/Colocalisation_Test.java b/src/main/java/sc/fiji/coloc/Colocalisation_Test.java
deleted file mode 100644
index 7b46ed6..0000000
--- a/src/main/java/sc/fiji/coloc/Colocalisation_Test.java
+++ /dev/null
@@ -1,765 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc;
-
-//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/sc/fiji/coloc/Colocalisation_Threshold.java b/src/main/java/sc/fiji/coloc/Colocalisation_Threshold.java
deleted file mode 100644
index 1e44ca4..0000000
--- a/src/main/java/sc/fiji/coloc/Colocalisation_Threshold.java
+++ /dev/null
@@ -1,969 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc;
-
-//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/sc/fiji/coloc/algorithms/Accumulator.java b/src/main/java/sc/fiji/coloc/algorithms/Accumulator.java
deleted file mode 100644
index 66eb70d..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/Accumulator.java
+++ /dev/null
@@ -1,109 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-import net.imglib2.TwinCursor;
-import net.imglib2.type.numeric.RealType;
-
-/**
- * A class allowing an easy accumulation of values visited by a
- * TwinCursor. After instantiation the sum of channel one,
- * channel two, products with them self and a product of both of
- * them will be available. It additionally provides the possibility
- * to subtract values from the data before the adding them to the
- * sum.
- *
- * @author Johannes Schindelin and Tom Kazimiers
- */
-public abstract class Accumulator> {
- protected double x, y, xx, xy, yy;
- protected int count;
-
- /**
- * The two values x and y from each cursor iteration to get
- * summed up as single values and their combinations.
- */
- public Accumulator(final TwinCursor cursor) {
- this(cursor, false, 0.0d, 0.0d);
- }
-
- /**
- * The two values (x - xDiff) and (y - yDiff) from each cursor
- * iteration to get summed up as single values and their combinations.
- */
- public Accumulator(final TwinCursor cursor, double xDiff, double yDiff) {
- this(cursor, true, xDiff, yDiff);
- }
-
- protected Accumulator(final TwinCursor cursor, boolean substract, double xDiff, double yDiff) {
- while (cursor.hasNext()) {
- cursor.fwd();
-
- T type1 = cursor.getFirst();
- T type2 = cursor.getSecond();
-
- if (!accept(type1, type2))
- continue;
-
- double value1 = type1.getRealDouble();
- double value2 = type2.getRealDouble();
-
- if (substract) {
- value1 -= xDiff;
- value2 -= yDiff;
- }
-
- x += value1;
- y += value2;
- xx += value1 * value1;
- xy += value1 * value2;
- yy += value2 * value2;
- count++;
- }
- }
-
- public abstract boolean accept(T type1, T type2);
-
- public double getX() {
- return x;
- }
-
- public double getY() {
- return y;
- }
-
- public double getXX() {
- return xx;
- }
-
- public double getXY() {
- return xy;
- }
-
- public double getYY() {
- return yy;
- }
-
- public int getCount() {
- return count;
- }
-}
diff --git a/src/main/java/sc/fiji/coloc/algorithms/Algorithm.java b/src/main/java/sc/fiji/coloc/algorithms/Algorithm.java
deleted file mode 100644
index 46ac322..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/Algorithm.java
+++ /dev/null
@@ -1,92 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import net.imglib2.type.numeric.RealType;
-
-import sc.fiji.coloc.gadgets.DataContainer;
-import sc.fiji.coloc.results.ResultHandler;
-import sc.fiji.coloc.results.Warning;
-
-/**
- * An algorithm is an abstraction of techniques like the
- * calculation of the Persons coefficient or Li'S ICQ. It
- * allows to separate initialization and execution of
- * such an algorithm.
- */
-public abstract class Algorithm> {
- // a name for the algorithm
- protected String name;
- /* a list of warnings that can be filled by the
- * execute method
- */
- List warnings = new ArrayList();
-
- public Algorithm(String name) {
- this.name = name;
- }
-
- /**
- * Executes the previously initialized {@link Algorithm}.
- */
- public abstract void execute(DataContainer container) throws MissingPreconditionException;
-
- public String getName() {
- return name;
- }
-
- /**
- * A method to give the algorithm the opportunity to let
- * its results being processed by the passed handler.
- * By default this methods passes the collected warnings to
- * the handler and sub-classes should make use of this by
- * adding custom behavior and call the super class.
- *
- * @param handler The ResultHandler to process the results.
- */
- public void processResults(ResultHandler handler) {
- for (Warning w : warnings)
- handler.handleWarning( w );
- }
-
- /**
- * Gets a reference to the warnings.
- *
- * @return A reference to the warnings list
- */
- public List getWarnings() {
- return warnings;
- }
-
- /**
- * Adds a warning to the list of warnings.
- *
- * @param shortMsg A short descriptive message
- * @param longMsg A long message
- */
- protected void addWarning(String shortMsg, String longMsg) {
- warnings.add( new Warning(shortMsg, longMsg) );
- }
-}
diff --git a/src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java b/src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java
deleted file mode 100644
index b91b356..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java
+++ /dev/null
@@ -1,349 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-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 sc.fiji.coloc.gadgets.DataContainer;
-import sc.fiji.coloc.gadgets.ThresholdMode;
-import sc.fiji.coloc.results.ResultHandler;
-
-/**
- * A class implementing the automatic finding of a threshold
- * used for Pearson 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 ratio 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/sc/fiji/coloc/algorithms/BisectionStepper.java b/src/main/java/sc/fiji/coloc/algorithms/BisectionStepper.java
deleted file mode 100644
index de98904..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/BisectionStepper.java
+++ /dev/null
@@ -1,91 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-/**
- * Try to converge a threshold based on an update value condition. If the
- * update value is larger zero, the threshold is lowered by half the distance
- * between the last thresholds. If the update value falls below zero or is not
- * a number, the threshold is increased by such a half step.
- *
- * @author Tom Kazimiers
- *
- */
-public class BisectionStepper extends Stepper {
- protected double threshold1;
- protected double threshold2;
- protected double thrDiff = Double.NaN;
- protected int iterations = 0;
- protected int maxIterations = 100;
-
- /**
- * Initialize the bisection stepper with a start threshold and its
- * last threshold
- *
- * @param threshold The current threshold
- * @param lastThreshold The last threshold
- */
- public BisectionStepper(double threshold, double lastThreshold) {
- threshold1 = threshold;
- threshold2 = lastThreshold;
- thrDiff = Math.abs(threshold1 - threshold2);
- }
-
- /**
- * Update threshold by a bisection step. If {@code value} is below zero or
- * not a number, the step is made upwards. If it is above zero, the stoep is
- * downwards.
- */
- @Override
- public void update(double value) {
- // update working thresholds for next iteration
- threshold2 = threshold1;
- if (Double.NaN == value || value < 0) {
- // we went too far, increase by the absolute half
- threshold1 = threshold1 + thrDiff * 0.5;
- } else if (value > 0) {
- // as long as r > 0 we go half the way down
- threshold1 = threshold1 - thrDiff * 0.5;
- }
- // Update difference to last threshold
- thrDiff = Math.abs(threshold1 - threshold2);
- // Update update counter
- iterations++;
- }
-
- /**
- * Get current threshold.
- */
- @Override
- public double getValue() {
- return threshold1;
- }
-
- /**
- * If the difference between both thresholds is < 1, we consider
- * that as reasonable close to abort the regression.
- */
- @Override
- public boolean isFinished() {
- return iterations > maxIterations || thrDiff < 1.0;
- }
-}
diff --git a/src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java b/src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java
deleted file mode 100644
index 61d70af..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java
+++ /dev/null
@@ -1,33 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-/**
- * A channel mapper should map an input value to either channel one or
- * channel two.
- *
- * @author Tom Kazimiers
- */
-public interface ChannelMapper {
- double getCh1Threshold(double t);
- double getCh2Threshold(double t);
-}
diff --git a/src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java b/src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java
deleted file mode 100644
index f092315..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java
+++ /dev/null
@@ -1,421 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-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 sc.fiji.coloc.gadgets.DataContainer;
-import sc.fiji.coloc.gadgets.DataContainer.MaskType;
-import sc.fiji.coloc.gadgets.Statistics;
-import sc.fiji.coloc.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 offset
- * @param size
- */
- 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 size
- */
- 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 size
- */
- 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/sc/fiji/coloc/algorithms/Histogram2D.java b/src/main/java/sc/fiji/coloc/algorithms/Histogram2D.java
deleted file mode 100644
index 79a8b0a..0000000
--- a/src/main/java/sc/fiji/coloc/algorithms/Histogram2D.java
+++ /dev/null
@@ -1,415 +0,0 @@
-/*-
- * #%L
- * Fiji's plugin for colocalization analysis.
- * %%
- * Copyright (C) 2009 - 2017 Fiji developers.
- * %%
- * 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
- * .
- * #L%
- */
-package sc.fiji.coloc.algorithms;
-
-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 sc.fiji.coloc.gadgets.DataContainer;
-import sc.fiji.coloc.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