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, 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/Colocalisation_Test.java b/src/main/java/sc/fiji/coloc/Colocalisation_Test.java
similarity index 99%
rename from src/main/java/Colocalisation_Test.java
rename to src/main/java/sc/fiji/coloc/Colocalisation_Test.java
index 0ed491f..7b46ed6 100644
--- a/src/main/java/Colocalisation_Test.java
+++ b/src/main/java/sc/fiji/coloc/Colocalisation_Test.java
@@ -1,763 +1,765 @@
/*-
* #%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/Colocalisation_Threshold.java b/src/main/java/sc/fiji/coloc/Colocalisation_Threshold.java
similarity index 99%
rename from src/main/java/Colocalisation_Threshold.java
rename to src/main/java/sc/fiji/coloc/Colocalisation_Threshold.java
index 95ff21c..1e44ca4 100644
--- a/src/main/java/Colocalisation_Threshold.java
+++ b/src/main/java/sc/fiji/coloc/Colocalisation_Threshold.java
@@ -1,967 +1,969 @@
/*-
* #%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/algorithms/Accumulator.java b/src/main/java/sc/fiji/coloc/algorithms/Accumulator.java
similarity index 98%
rename from src/main/java/algorithms/Accumulator.java
rename to src/main/java/sc/fiji/coloc/algorithms/Accumulator.java
index bf7923b..66eb70d 100644
--- a/src/main/java/algorithms/Accumulator.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/Accumulator.java
@@ -1,109 +1,109 @@
/*-
* #%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 algorithms;
+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/algorithms/Algorithm.java b/src/main/java/sc/fiji/coloc/algorithms/Algorithm.java
similarity index 93%
rename from src/main/java/algorithms/Algorithm.java
rename to src/main/java/sc/fiji/coloc/algorithms/Algorithm.java
index 5b17c5e..46ac322 100644
--- a/src/main/java/algorithms/Algorithm.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/Algorithm.java
@@ -1,92 +1,92 @@
/*-
* #%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 algorithms;
-
-import gadgets.DataContainer;
+package sc.fiji.coloc.algorithms;
import java.util.ArrayList;
import java.util.List;
import net.imglib2.type.numeric.RealType;
-import results.ResultHandler;
-import results.Warning;
+
+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/algorithms/AutoThresholdRegression.java b/src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java
similarity index 98%
rename from src/main/java/algorithms/AutoThresholdRegression.java
rename to src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java
index 532c016..b91b356 100644
--- a/src/main/java/algorithms/AutoThresholdRegression.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/AutoThresholdRegression.java
@@ -1,348 +1,349 @@
/*-
* #%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 algorithms;
+package sc.fiji.coloc.algorithms;
-import gadgets.DataContainer;
-import gadgets.ThresholdMode;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.TwinCursor;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.view.Views;
-import results.ResultHandler;
+
+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/algorithms/BisectionStepper.java b/src/main/java/sc/fiji/coloc/algorithms/BisectionStepper.java
similarity index 98%
rename from src/main/java/algorithms/BisectionStepper.java
rename to src/main/java/sc/fiji/coloc/algorithms/BisectionStepper.java
index 8970c25..de98904 100644
--- a/src/main/java/algorithms/BisectionStepper.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/BisectionStepper.java
@@ -1,91 +1,91 @@
/*-
* #%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 algorithms;
+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/algorithms/ChannelMapper.java b/src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java
similarity index 96%
rename from src/main/java/algorithms/ChannelMapper.java
rename to src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java
index f26c564..61d70af 100644
--- a/src/main/java/algorithms/ChannelMapper.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/ChannelMapper.java
@@ -1,33 +1,33 @@
/*-
* #%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 algorithms;
+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/algorithms/CostesSignificanceTest.java b/src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java
similarity index 98%
rename from src/main/java/algorithms/CostesSignificanceTest.java
rename to src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java
index e2363ab..f092315 100644
--- a/src/main/java/algorithms/CostesSignificanceTest.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/CostesSignificanceTest.java
@@ -1,421 +1,421 @@
/*-
* #%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 algorithms;
-
-import gadgets.DataContainer;
-import gadgets.DataContainer.MaskType;
-import gadgets.Statistics;
+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 results.ResultHandler;
+
+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/algorithms/Histogram2D.java b/src/main/java/sc/fiji/coloc/algorithms/Histogram2D.java
similarity index 99%
rename from src/main/java/algorithms/Histogram2D.java
rename to src/main/java/sc/fiji/coloc/algorithms/Histogram2D.java
index 33f749b..79a8b0a 100644
--- a/src/main/java/algorithms/Histogram2D.java
+++ b/src/main/java/sc/fiji/coloc/algorithms/Histogram2D.java
@@ -1,414 +1,415 @@
/*-
* #%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 algorithms;
+package sc.fiji.coloc.algorithms;
-import gadgets.DataContainer;
import ij.measure.ResultsTable;
import java.util.EnumSet;
import net.imglib2.RandomAccess;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.TwinCursor;
import net.imglib2.img.ImgFactory;
import net.imglib2.img.array.ArrayImgFactory;
import net.imglib2.type.logic.BitType;
import net.imglib2.type.numeric.RealType;
import net.imglib2.type.numeric.integer.LongType;
import net.imglib2.view.Views;
-import results.ResultHandler;
+
+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