Mercurial > repos > bgruening > woundhealing_scratch_assay
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; +} +