diff measureWoundClosing.groovy @ 0:8948cc562b7c draft default tip

planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/tools commit 4d4d0b10eb6be3b1f13b2becc8f057ac41d3f0de
author bgruening
date Thu, 29 Feb 2024 17:43:50 +0000
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/measureWoundClosing.groovy	Thu Feb 29 17:43:50 2024 +0000
@@ -0,0 +1,232 @@
+/**
+ * This script runs in Fiji
+ *
+ * It needs the following Fiji update sites:
+ * - IJPB-Plugins (MorpholibJ)
+ */
+
+
+import fiji.threshold.Auto_Threshold
+import ij.IJ
+import ij.ImagePlus
+import ij.gui.Roi
+import ij.measure.Measurements
+import ij.measure.ResultsTable
+import ij.plugin.FolderOpener
+import ij.plugin.ImageCalculator
+import ij.plugin.filter.Analyzer
+import ij.plugin.frame.RoiManager
+import inra.ijpb.binary.BinaryImages
+import inra.ijpb.morphology.Morphology
+import inra.ijpb.morphology.Reconstruction
+import inra.ijpb.morphology.strel.SquareStrel
+import inra.ijpb.segment.Threshold
+import net.imagej.ImageJ
+
+// INPUT UI AND CLI
+//
+#@ File (label="Input directory", style="directory") inputDir
+#@ String (label="Dataset ID") datasetId
+#@ Double (label="CoV threshold (-1: auto)", default="-1") threshold
+#@ Boolean (label="Run headless", default="false") headless
+#@ Boolean (label="Save results", default="true") saveResults
+#@ String (label="Output directory name (will be created next to input directory)", default="analysis") outDirName
+
+// INPUT FOR TESTING WITHIN IDE
+//
+//def inputDir = new File("/Users/tischer/Documents/wound-healing-htm-screen/data/input")
+//def datasetId = "A3ROI2_Slow"; // C4ROI1_Fast A3ROI2_Slow
+//def outDirName = "analysis"
+//def threshold = (double) -1.0 // auto
+//def headless = false;
+//def saveResults = false;
+//new ij.ImageJ().setVisible(true)
+
+// FIXED PARAMETERS
+//
+def cellDiameter = 20
+def scratchDiameter = 500
+def binningFactor = 2
+
+// DERIVED PARAMETERS
+//
+def cellFilterRadius = cellDiameter/binningFactor
+def scratchFilterRadius = scratchDiameter/binningFactor
+
+println("Cell filter radius: " + cellFilterRadius)
+println("Scratch filter radius: " + scratchFilterRadius)
+
+// CODE
+//
+
+// open the images
+//
+IJ.run("Close All");
+println("Opening: " + datasetId)
+def imp = FolderOpener.open(inputDir.toString(), " filter=(.*"+datasetId+".*)")
+println("Number of slices: " + imp.getNSlices())
+
+if ( imp == null  || imp.getNSlices() == 0 ) {
+    println("Could not find any files matching the pattern!")
+    System.exit(1)
+}
+
+// process the images to enhance regions with cells
+// using the fact the the cell regions have a higher local variance
+//
+println("Process images to enhance regions containing cells...")
+// remove scaling to work in pixel units
+IJ.run(imp,"Properties...", "pixel_width=1 pixel_height=1 voxel_depth=1");
+// bin to save compute time
+IJ.run(imp, "Bin...", "x=" + binningFactor + " y=" + binningFactor + " z=1 bin=Average");
+def binnedImp = imp.duplicate() // keep for saving
+// enhance cells
+IJ.run(imp, "32-bit", "");
+def sdevImp = imp.duplicate()
+sdevImp.setTitle(datasetId + " sdev" )
+IJ.run(sdevImp, "Find Edges", "stack"); // removes larger structures, such as dirt in the background
+IJ.run(sdevImp, "Variance...", "radius=" + cellFilterRadius + " stack");
+IJ.run(sdevImp, "Square Root", "stack");
+// mean
+def meanImp = imp.duplicate()
+meanImp.setTitle(datasetId + " mean")
+IJ.run(meanImp, "Mean...", "radius=" + cellFilterRadius + " stack");
+// cov
+def covImp = ImageCalculator.run(sdevImp, meanImp, "Divide create 32-bit stack");
+IJ.run(covImp, "Enhance Contrast", "saturated=0.35");
+IJ.run(covImp, "8-bit", ""); // otherwise the thresholding does not seem to work
+covImp.setTitle(datasetId + " cov" )
+if (!headless) covImp.duplicate().show();
+
+// create binary image (cell-free regions are foreground)
+//
+println("Creating binary image of cell-free regions...")
+// configure black background
+IJ.run("Options...", "iterations=1 count=1 black");
+// determine threshold in first frame, because there we are
+// closest to a 50/50 occupancy of the image with signal,
+// which is best for most auto-thresholding algorithms
+covImp.setPosition(1)
+def histogram = covImp.getProcessor().getHistogram()
+// Auto threshold with Huang method and
+// multiply the threshold with a fixed factor (as is done in CellProfiler),
+// based on the observation that the threshold is consistently
+// a bit too high, which may be due to
+// the fact that the majority of the image is foreground
+if ( threshold == -1 )
+  threshold = Auto_Threshold.Huang(histogram) * 0.8
+println("Threshold: " + threshold)
+// create binary image of whole movie,
+// using the threshold of the first image
+// defining the cell free regions as foreground
+def binaryImp = (ImagePlus) Threshold.threshold(covImp, 0.0, threshold)
+// dilate the cell free regions, because due to the cell filter radius
+// the cell sizes are over estimated (blurred into cell free regions)
+binaryImp = Morphology.dilation(binaryImp, SquareStrel.fromRadius((int) cellFilterRadius))
+binaryImp.setTitle(datasetId + " binary")
+if(!headless) binaryImp.duplicate().show()
+
+// create scratch ROIs
+//
+println("Creating scratch ROI...")
+binaryImp.setPosition(1) // scratch is most visible in first frame
+def scratchIp = binaryImp.crop("whole-slice").getProcessor().duplicate();
+// identify largest cell free region as scratch region
+scratchIp = BinaryImages.keepLargestRegion(scratchIp)
+// remove cells inside scratch region
+scratchIp = Reconstruction.fillHoles(scratchIp)
+if(!headless) new ImagePlus("Scratch", scratchIp.duplicate()).show()
+// disconnect from cell free regions outside scratch
+scratchIp = Morphology.opening(scratchIp, SquareStrel.fromRadius((int)(scratchFilterRadius/20)))
+// in case the morphological opening cut off some cell free
+// areas outside the scratch we again only keep the largest region
+scratchIp = BinaryImages.keepLargestRegion(scratchIp)
+// smoothen scratch edges
+scratchIp = Morphology.closing(scratchIp, SquareStrel.fromRadius((int)(scratchFilterRadius/2)))
+
+// convert binary image to ROI, which is handy for measurements
+def scratchImp = new ImagePlus("Finale scratch", scratchIp)
+if(!headless) scratchImp.show()
+IJ.run(scratchImp, "Create Selection", "");
+def scratchROI = scratchImp.getRoi()
+
+// measure occupancy of scratch ROI
+// `area_fraction` returns the fraction of foreground pixels
+// (cell free area) within the measurement ROI
+println("Performing measurements...")
+IJ.run("Set Measurements...", "area bounding area_fraction redirect=None decimal=2");
+def rt = RoiManager.multiMeasure(binaryImp, new Roi[]{scratchROI}, false)
+
+// show results
+//
+if ( !headless ) {
+    rt.show("Results")
+    binnedImp.show()
+    binnedImp.setTitle(datasetId + " binned")
+    binnedImp.setRoi(scratchROI, true)
+    binaryImp.setPosition(1)
+    binaryImp.show()
+    binaryImp.setTitle(datasetId + " binary")
+    binaryImp.setRoi(scratchROI, true)
+}
+
+// save results
+//
+if ( saveResults ) {
+    // create output directory next to input directory
+    def outputDir = new File(inputDir.getParent(), outDirName);
+    println("Ensuring existence of output directory: " + outputDir)
+    outputDir.mkdir()
+    // save table
+    rt.save(new File(outputDir, datasetId + ".csv").toString());
+    // save binned image with ROI
+    binnedImp.setRoi(scratchROI, false)
+    IJ.save(binnedImp, new File(outputDir, datasetId + ".tif").toString());
+}
+
+println("Analysis of "+datasetId+" is done!")
+if ( headless ) System.exit(0)
+
+// FUNCTIONS
+//
+
+// copied from ImageJ RoiManager because
+// https://forum.image.sc/t/make-multimeasure-public-in-roimanager/69273
+private static ResultsTable multiMeasure(ImagePlus imp, Roi[] rois) {
+    int nSlices = imp.getStackSize();
+    Analyzer aSys = new Analyzer(imp); // System Analyzer
+    ResultsTable rtSys = Analyzer.getResultsTable();
+    ResultsTable rtMulti = new ResultsTable();
+    rtMulti.showRowNumbers(true);
+    rtSys.reset();
+    int currentSlice = imp.getCurrentSlice();
+    for (int slice=1; slice<=nSlices; slice++) {
+        int sliceUse = slice;
+        if (nSlices==1) sliceUse = currentSlice;
+        imp.setSliceWithoutUpdate(sliceUse);
+        rtMulti.incrementCounter();
+        if ((Analyzer.getMeasurements()& Measurements.LABELS)!=0)
+            rtMulti.addLabel("Label", imp.getTitle());
+        int roiIndex = 0;
+        for (int i=0; i<rois.length; i++) {
+            imp.setRoi(rois[i]);
+            roiIndex++;
+            aSys.measure();
+            for (int j=0; j<=rtSys.getLastColumn(); j++){
+                float[] col = rtSys.getColumn(j);
+                String head = rtSys.getColumnHeading(j);
+                String suffix = ""+roiIndex;
+                Roi roi = imp.getRoi();
+                if (roi!=null) {
+                    String name = roi.getName();
+                    if (name!=null && name.length()>0 && (name.length()<9||!Character.isDigit(name.charAt(0))))
+                        suffix = "("+name+")";
+                }
+                if (head!=null && col!=null && !head.equals("Slice"))
+                    rtMulti.addValue(head+suffix, rtSys.getValue(j,rtSys.getCounter()-1));
+            }
+        }
+    }
+    return rtMulti;
+}
+