/* Macro to quantify foci in nuclei/cells. Works on 2D/3D images, including multiseries files, as long as all series have the same specs * More info on https://imagej.net/plugins/foci-analyzer * * ► Requires the following Fiji update sites: * - 3D ImageJ Suite * - CLIJ * - CLIJ2 * - CLIJx-assistent * - CLIJx-assistent-extensions * - CSBDeep * - IJPB-plugins * - StarDist *(- TensorFlow) In case StarDist gives an error. See https://forum.image.sc/t/stardist-error-since-update/107729 * * ► In order to run Cellpose segmentation you also need: * - A working Cellpose Python environment * - PTBIOP update site, with proper settings. See https://github.com/BIOP/ijl-utilities-wrappers/blob/master/README.md#cellpose * * * Author: Bram van den Broek, The Netherlands Cancer Institute * For questions please use the Image.sc forum (https://forum.image.sc/) with tag @bramvdbroek. * * * * Changelog * --------- * version 0.94, April 2022: * - Invert foci before MorpholibJ marker-controlled watershed -> get rid of holes in the foci. * - Convert foci to 32-bit before DoG filtering. Currently not converted back to 16 bit when determining foci thresholds. * Both changes will substantially impact the measured foci from version 0.93! * * Version 0.95, April 2022: * - Fixed overwriting of output files for multi-series images (e.g. .lif, .czi) * - Fixed mistakes in output metadata (foci size A and B were both A, 2 x 'threshold bias A' in the key) * - Added nucleus 2D area as output column * * Version 0.96, April 2022: * - Added Cellpose as segmentation possibility (requires PTBIOP update site and a working Cellpose installation in Python) * * Version 0.97, May 2022: * - Added possibility to bin the image in XY (all channels) * * Version 0.98, June 2022: * - Fixed a bug: foci detection on max projections with 2 foci channels does not give an error any more * * Version 0.99, June 2022: * - Updated the parameters saved in the metadata of the output images * * Version 1.0, July 2022: * - Added possibility to manually remove segmented nuclei by clicking. Masks are saved and can be loaded for re-analysis. * * Version 1.1, July 2022: * - Fixed a small bug: nuclear intensity in foci channel B is now measured also when detectFociChannelB is false. * - Changes in AreaMaximum detection (minimum foci size) * - Now asks to create output folder if it doesn't exist. * - Classic segmentation works again * - Fixed a bug where the nucleus area was not correctly reported when Stardist nuclei binning factor > 1 * * Version 1.2, October 2022: * - Fixed critical bug where after threshold optimization of the second image Fiji would crash (because an image was released from the GPU) * - Improved log window readability * - Fixed a bug (in debug mode) where Fiji crashed after threshold optimization if the overlay image was not selected * * Version 1.22, October 2022: * - Fixed a bug where manual nuclei removal didn't work on max projections * - Added labels to channels in the overlay image * - Implemented a maximum tile size for Stardist * * Version 1.23, October 2022: * - Added option to crop image borders * * Version 1.24, December 2022: * - Added total nucleus intensity (mean * area) as output column * - Save cropBorders and XYBinning in the output image metadata * - check if Foci channel B actually exists (before it would just take the last channel) and is not the same as Foci channel A * * Version 1.25, December 2022 * - Fixed bug where foci spots were not displayed correctly when threshold was changed * - Improved visualization options * * Version 1.3, December 2022 * - If 'Alt' is pressed before a file is loaded, the user can select a ROI to do a preview analysis on. * - downscaled ROIs are now smoothed using spline fitting * - Improved GUI: 3D image handling radiobuttons instead of separate checkboxes * - Improved logging: Threshold factor * bias is now mentioned in the log window * - Colocalization overlap image is now also 3D and saved as .Gif file * * Version 1.32, April 2023 * - Added the actually used threshold value to the output image metadata * - Included option to detect foci on extended depth of focus projection (Sobel, GPU) * * Version 1.35, June 2023 * - Added option to include foci outside nuclei. Distance range can be specified. * * Version 1.36, February 2024 * - Fixed a bug causing a crash with 32-bit images when producing the overlay. * - Added a cellpose diameter parameter * * Renamed to version 1.4 for release on GitHub: https://github.com/BioImaging-NKI/Foci-analyzer/ * * Version 1.41, May 2024 * - Added possibility to load 2D nuclei segmentations from a .zip file containing ImageJ ROIs (e.g. exported from QuPath) * * Version 1.5, September 2024 * - Added foci overlap count and area/volume for both channels to the result, as well as in the overlap map image (thanks to Mabel Baxter Dalrymple) * - Overlay and overlap images are now saves as .zip files to save storage space (typically 10 times), at the expense of slower processing. * - Possibility for automatic rescaling nuclei for StarDist (to a pixel size 0.4 µm), with non-integer rescaling factor. * Upscaling (downscale factor < 1) is now also allowed. Interpolation is set to bilinear. * - Settings optimization: Channels Tool and B&C window are displayed in upper left corner of the screen. Added some tips in the optimization dialog. * * Version 1.52, September 2024 * - Added a 'Help' button in the optimization dialog, linking to the GitHub site. * - Added the possibility (in the optimization dialog) to use the same fixed absolute threshold for all images. * - Improved starting dialog, with headings * * Version 1.53, November 2024 * - Added possibility to fix the absolute threshold value for all images. Changed the optimization dialog for this. * - Channels and slices can now be flipped in images with multiple slices but only a single channel. * * Version 1.54, December 2024 * - Expanded Cellpose detection options. * * Version 1.59, December 2024 * - Measure foci using MorphoLibJ instead of CLIJ2. Faster and a bit easier to handle. * - A table with all foci statistics is now also saved (thanks to Harry Osborne) * - Added possibility to load settings from a previously analyzed file * * Version 1.6-1.7, January-March 2025 * - Cellpose envPath and envType parameters are automatically retrieved (new Cellpose wrapper) * - Added 3D segmentation with Cellpose 2.5D/3D, and 3D label visualization. Cellpose 3.1.0 is required (for --dP_smooth/flow3D_smooth parameter) * - Changed automatic StarDist downscaling to 0.25 um/pixel (was 0.4) * - Persistent foci detection script parameters are now saved after running the optimization dialog * - Fixed a bug where optimization dialog would be skipped (thanks to Harry Osborne) * - Foci centroids of removed foci due to size restrictions are now removed as well * - Improved filtering in Z -> better foci detection in 3D * - Fixed bug where the median foci volume was given in pixels, thanks to Harry Osborne * - Fixed bug where the foci area was incorrect for 3D images when analyzed in 2D (foci area was multiplied by voxelDepth) * - Input images can now be time-lapse series. Frames are split and processed sequentially. * - Option to load label images instead of performing segmentation. (Label images need to have the same name as the original images, in a different folder.) * - If ROI/label images is not found it skips the image. * - Foci count can be overlayed on the image, as well as the nucleus number * - Added StdDev projection method for handling 3D images * - Use own code for ROI manager to labelmap (instead of SCF MPI CBG plugin) * * Version 1.74: * - Possibility to choose a color for foci count overlay * - Possibility to choose 'none' for nucleus ID / foci count overlay * * Version 1.8: * - Added checkbox to hide the 3D segmentation dialog window for subsequent images. * - Fixed bug where the macro would crash when 3D Cellpose segmentation was chosen in combination with *not* excluding cells on edges. * - Release version for Elmi 2025 * * Version 1.83: * - Added foci centroid coordinates to the 'All foci statistics' table (in 1.82 - bugfix in 1.83) * - Added descriptions to the Scijava script parameters * * * --> TO DO: open timelapse output files, concatenate, and restore overlays. Now via external scripts. * */ #@ String Foci_Analyzer_message (value="
Foci Analyzer (v1.84)
File settings
", visibility="MESSAGE") #@ File[] files (label = "Input files", style="File", description="Here you can specify which files to analyze, by adding them to the list, or drag&drop from a file explorer window.") #@ File outputFolder (label = "Output folder", style = "directory", description="The folder where all the analyzed images and results will be written.") #@ String image_message (value="Image settings
", visibility="MESSAGE") #@ Boolean loadSettingsFromFile (label = "Load settings from previously analyzed image?", description="When checked, a separate dialog will popup where an output .zip file can be selected (overlay or colocalization map).\nAll settings are loaded from the metadata in the saved image. The script parameter entries in this large dialog are ignored.") #@ Integer nucleiChannel (label = "Nuclei channel (-1 if not used)", value = 1, description="The image channel that contains the nuclei. For StarDist nuclei segmentation is performed using this channel, which always happens in 2D (in the case of 3D images on a maximum intensity projection).\nFor Cellpose, there are two possibilities, depending on the value of Cytoplasm/membrane channel below.") #@ Integer cytoChannel (label = "Cytoplasm/membrane channel (-1 if not used)", value = -1, description="The image channel that contains cytoplasm or membrane. If set to -1, segmentation is performed on the nucleus channel alone.\nIf not -1 and Cellpose is chosen as segmentation method, segmentation is performed on this channel. In this case the nuclei channel is the 'additional helper channel'.\nIf not set to -1 and StarDist is chosen, this channel is not used in the segmentation, but instead just displayed in the overlay image. (default: -1)") #@ Integer fociChannelA (label = "Foci channel A", value = 2, description="The first foci channel (default: 2)") #@ Integer fociChannelB (label = "Foci channel B", value = 3, description="The second foci channel (default: 3)") #@ Boolean detectFociChannelB (label = "Also detect foci in channel B and perform colocalization?", persistence=true, value=true, description="If checked, foci in both channels A and B will be analyzed, followed by a simple colocalization analysis. (default: checked)") #@ String ThreeDHandling (label = "2D/3D foci detection", choices={"Detect foci in 3D (or 2D when N/A)", "Detect foci on the Maximum Intensity Projection", "Detect foci on a Extended Depth of Focus Projection", "Detect foci on the Standard Deviation Projection", "Detect foci on the Summed Intensity Projection", "Use quasi-2D foci detection (detect foci separately in every Z-slice)", "Process a single z-slice only (specify below which slice)"}, style="listbox", description="2D/3D foci detection:\nThis parameter determines how foci in 3D images should be analyzed. For 2D input images this setting is ignored. The options are:\n\n- Detect foci in 3D (or 2D when N/A) (default) performs foci analysis using 3D filters and 3D marker-controlled watershed functions. Connected foci in consecutive slices are counted once.\n\n- Detect foci on the Maximum Intensity Projection performs 2D foci analysis on the MIP projection.\n\n- Detect foci on a Extended Depth of Focus Projection performs 2D foci analysis on an EDF projection.\n\n- Detect foci on the Standard Deviation Projection performs 2D foci analysis on the STDEV projection.\n\n- Detect foci on the Summed Intensity Projection performs 2D foci analysis on the SUM projection.\n\n- Use quasi-2D foci detection (detect foci separately in every Z-slice) analyzes every z-slice in a 3D image as a separate 2D image.\n This setting is useful in cases where the z-spacing is very large and each focus is visible in only one z-slice.\n Hence, connected foci in consecutive slices will be counted multiple times.\n\n- Process a single z-slice only (specify below which slice) allows the user to analyze foci only in a particular z-slice.") #@ Integer singleSlice (label = "[single z-slice foci detection only] Slice nr", value=1, min=1, description="The single z-slice used for the option above. For any other choice this parameter is ignored.") #@ Integer cropBorder (label = "Remove image borders (pixels)", value = 0, min=0, description="Possibility to remove edges from the image. This can in particular be useful when the image edges have very different intensities, causing incorrect automatic nuclei segmentation. (default: 0)") #@ Integer XYBinning (label = "Image XY binning before analysis [1-n]", value = 1, min=1, description="Optional pixel binning in case the resolution is very high and the foci consist of many pixels.\nA value of 2 means: 2x2 pixels will be binned into 1 pixel. This reduces noise in the image and speeds up analysis. (default: 1)") #@ String nuclei_message (value = "Nuclei/cell detection settings
", visibility="MESSAGE") #@ String nucleiSegmentationChoice (label = "Nuclei/cell segmentation method", choices={"StarDist nuclei segmentation 2D (or on 2D projection)", "Cellpose segmentation 2D (or on 2D projection)", "Cellpose segmentation 3D", "Classic segmentation", "Load ROIs from file", "Load label images"}, style="listBox", description="Nuclei/cell segmentation method:\nStardist nuclei segmentation 2D (or on 2D projection) (default) uses the pretrained convolutional neural network StarDist to recognize cell nuclei in fluorescence microscopy images.\nIn general this method works very well on a large variety of samples.\n\n- Cellpose segmentation 2D (or on 2D projection) uses the deep learning network Cellpose to recognize whole cells or nuclei.\nUse this option if you want to measure foci in entire cells, or if you prefer Cellpose nuclei segmentation over StarDist.\nN.B. Cellpose requires additional installations (see Installation / Requirements).\n\n- Cellpose segmentation 3D: If this option is chosen a new dialog pops up with extra settings. These are the most important parameters for 3D segmentation.\nMore parameters can be added in the 'Additional Cellpose parameters' field. The Help button takes you to the Cellpose CLI with explanations of all parameters.\n\n- Classic nuclei segmentation allows the user to segment nuclei using manual/automatic thresholding is provided for legacy reasons.\nThe method is almost always outperformed by the two other methods.\n\n- Load ROIs from file: ImageJ ROI .zip files can be loaded instead of performing segmentation. This option is used in the (near future) QuPath-Fiji workflow.\nROI files should have the same name as the input images without extensions, followed by '_ROIs.zip'.\n\n- Load label images allows loading a labelmap, if the segmentation has been done by external programs, or to quickly re-run files with the same segmentation.\nLabel image files should be present in another folder and have the exact same name as the input images.") #@ Double downsampleFactorStarDist (label = "Stardist nuclei rescaling factor [1-n], 0 for automatic, 1 for no rescaling", value = 0, min=0, description="Stardist is trained on medium resolution images, and generally performs well on images with pixel sizes around 0.2-0.5 µm.\nSet to 0 for automatic rescaling the nuclei to an optimal pixel size of 0.25 µm, or put any other number for manual control of the rescaling.") #@ Double probabilityThreshold (label = "Probability/flow threshold [0.0-1.0] (StarDist/Cellpose)", value = 0.5, min=0, max=1, style="format:0.0", description="Lower values will accept more nuclei/cells; higher values will be more stringent. For Cellpose this is actually the flow_threshold parameter.") #@ String CellposeModel (label = "Cellpose model", choices={"cyto3","nuclei","tissuenet_cp3","cpsam","custom"}, style="listBox", value="cyto3", description="The model (built-in or custom) used for segmentation.") #@ Integer CellposeDiameter (label = "Cellpose cell diameter (pixels), 0 for automatic", value = 0, min=0, description="Estimated diameter of the cells, in pixels. Setting this parameter to 0 will trigger Cellpose to estimate it.") #@ Integer minNucleusSize_setting (label = "Remove nuclei/cells with diameter smaller than (µm)", value = 4, description="Objects smaller than circles having an area corresponding to this diameter will be removed.") #@ Integer maxNucleusSize_setting (label = "Remove nuclei/cells with diameter larger than (µm)", value = 50, description="Likewise, but for large objects.") #@ Boolean excludeOnEdges (label = "Exclude nuclei/cells on image edges", value = true, description="When checked, nuclei that touch the image edge will not be analyzed. (default: checked).") #@ String manualRemoveNuclei (label = "Manually remove segmented nuclei/cells", choices={"No thanks","Manually remove nuclei", "Load previously saved removals (from output folder)"}, value = "No thanks", description="Manual nuclei removal: allows the user to erase ill-segmented nuclei before analysis. (default: No thanks)\nOptions:\n\n- No thanks means no manual nuclei editing\n\n- Manually remove nuclei : Remove nuclei by leftclicking them in the image with the mouse. Editings will be saved to a small text file in the output folder.\n\n- Load previously saved removals (from output folder) : If you have edited the segmented nuclei on this image before, it will load the previous nuclei removals\n from the file in the specified output folder. (Hence, if you change the output folder parameter this option will not work.)") #@ String foci_message1 (value="Foci detection settings
", visibility="MESSAGE") #@ Boolean optimizationModeSetting (label = "Preview foci detection (for parameter fine-tuning)?", value=true, description="Checking this will allow the user to adapt the foci detection settings on a preview analysis before quantifying.") #@ String fociSizeA (label = "Foci size channel A (after XY binning)", choices={"tiny","small","average","large","huge","other"}, style="listBox", value="average", description="This parameter controls several foci image filtering steps and steers the macro towards detecting smaller or larger foci.") #@ String fociSizeB (label = "Foci size channel B (after XY binning)", choices={"tiny","small","average","large","huge","other"}, style="listBox", value="average", description="This parameter controls several foci image filtering steps and steers the macro towards detecting smaller or larger foci.") #@ String fociDetectionMethod (label = "Foci detection method", choices={"Marker-controlled watershed (recommended)","AreaMaxima local maximum detection"}, style="listBox", description="Select the method for foci detection.") #@ Double thresholdFactorA (label = "Foci intensity threshold bias channel A", value = 0.0, min=-2.5, max=2.5, style="scroll bar", stepSize=0.01, description="The macro will automatically estimate the intensity threshold for foci detection. This default threshold can be biased with the slider.") #@ Double thresholdFactorB (label = "Foci intensity threshold bias channel B", value = 0.0, min=-2.5, max=2.5, style="scroll bar", stepSize=0.01, description="The macro will automatically estimate the intensity threshold for foci detection. This default threshold can be biased with the slider.") #@ Integer minFociSize (label = "Minimum foci size (area) (pixels/voxels)", value = 3, description="Foci occupying an area/volume smaller than this value will be deleted.") #@ Integer maxFociSize (label = "Maximum foci size (area) (pixels/voxels)", value = 9999, description="The upper limit for the foci size, in pixels/voxels.") #@ Double maxFociDistanceOutsideNuclei_setting (label = "Max distance of foci outside nuclei/cells (µm); -1 for full image", value = 0, min=-1, style="format:0.0", description="This controls how far outside the cell/nucleus segmentation foci should still be counted. (default: 0)") #@ Integer minOverlapSize (label = "Minimum overlap of foci to colocalize (pixels/voxels)", min = 1, value = 1, description="Foci in channel A and B will be counted as colocalizing only if they overlap with at least this area/volume.") #@ String visualization_message (value = "Visualization options
", visibility="MESSAGE") #@ String overlayChoice (label = "Nuclei/cell overlay choice", choices={"nucleus/cell ID","foci count","both","none"}, description="Select which numbers are added to the overlay.") #@ ColorRGB fontColorCells (label = "Nuclei/cell label color", value="orange", description="The color of the nucleus/cell ID text overlay.") #@ ColorRGB fontColorFoci (label = "Foci count label color", value="red", description="The color of the foci count text overlay.") #@ Integer labelFontSize (label = "Nuclei/cell label font size", value=12, min=4, description="The size of the nuclei/cell ID text overlay.") #@ String overlayBrightness (label = "Nuclei/cell outline brightness", choices={"bright","dim"}, description="The brightness of the nuclei outlines overlay.") #@ String outlineColor (label = "Nuclei/cell outline color", choices={"White","Red","Green","Blue","Cyan","Magenta","Yellow"}, value = "Cyan", description="The color of the nuclei outlines overlay.") #@ Boolean debugMode (label = "Debug mode (show intermediate images)", value=false, description="Used for development and bug fixing: checking this option will trigger displaying many intermediate results during the processing. It will also slow down the analysis.") #@ String file_and_image_message0 (value = "Need help? Visit the Foci Analyzer website on ImageJ.net
", visibility="MESSAGE") version = 1.84; requires("1.54i"); //Minimum required ImageJ version /* KNOWN ISSUES * * ! Maka a possibility to skip the 3D Cellpose parameters dialog * ! Foci outside nuclei: Make isotropic - dilate - make non-isotropic * ! 3D outlines also in colocalization image * ! Make filter before classic segmentation pixelsize-dependent * * * TO DO | IDEAS * * * Include 3D Cellpose settings when loading settings from previously analyzed image * * * Output: * (Centroids of all foci? (MorphoLibJ Analyze regions?)) * (Distance from edge for all foci?) * * * Include Voronoi-Otsu segmentation! * * * Use Roi.setMinStrokeWidth() for a better scaling of outlines? * * * Do not remove lines in log window when areaMaximum detection is used. * * * Use ALT/right mouse click to pop up options menu (background subtraction, nuclear segmentation, foci size, etc. * and allow the user to determine the right settings! See MosaicExplorerJ macro. * Currently, pressing ALT during opening the image allows for making a selection. * * * Add more: description="..." to the options (onmouseover) * * * Background intensity is now very crude (just the mean outside the nuclei). So: mask the background, set (percentile) threshold, measure mean/median. * */ //Check if all dependency Update Sites are installed: // * - 3D ImageJ Suite // * - CLIJ // * - CLIJ2 // * - CLIJx-assistent // * - CLIJx-assistent-extensions // * - CSBDeep // * - IJPB-plugins // * - PT-BIOP // * - StarDist missingPlugin = ""; List.setCommands; if (List.get("CLIJ2 Macro Extensions")=="") missingPlugin += "CLIJ2, "; if (List.get("3D Manager")=="") missingPlugin += "3D Image Suite, "; if (List.get("MorphoLibJ Marker-controlled Watershed (CLIJx, experimental)")=="") missingPlugin += "CLIJx, "; if (List.get("Intensity Measurements 2D/3D")=="") missingPlugin += "MorphoLibJ (IJPB-Plugins), "; if (List.get("Run your network")=="") missingPlugin += "CDBDeep, "; if (List.get("Command From Macro")=="") missingPlugin += "StarDist, "; if (missingPlugin != "") { print("\\Clear"); missingPlugin = missingPlugin.substring(0, missingPlugin.length-2); print("Error: Required plugin(s) not found:\n"+missingPlugin+"\n \nFoci Analyzer requires the following Fiji Update Sites to be activated:\n* 3D ImageJ Suite\n* CLIJ\n* CLIJ2\n* CLIJx-assistent\n* CLIJx-assistent-extensions\n* CSBDeep\n* IJPB-plugins\n* StarDist\n \nGo to Help -> Update... -> Manage Update Sites and check the relevant boxes."); exit("Error: Required plugin(s) not found:\n"+missingPlugin+"\n \nFoci Analyzer requires the following Fiji Update Sites to be activated:\n* 3D ImageJ Suite\n* CLIJ\n* CLIJ2\n* CLIJx-assistent\n* CLIJx-assistent-extensions\n* CSBDeep\n* IJPB-plugins\n* StarDist\n \nGo to Help -> Update... -> Manage Update Sites and check the relevant boxes.\nThis info is also printed to the Log Window."); } if(loadSettingsFromFile == true) { //TO DO: Add 3D Cellpose settings setBatchMode(true); loadSettingsPath = File.openDialog("Select a previously analyzed output file (.zip) to load settings from"); if(endsWith(loadSettingsPath, ".zip")) open(loadSettingsPath); else exit("Please select an output .zip file from a previous analysis. Exiting macro."); print("Loading settings from "+loadSettingsPath); loadVersion = Property.getNumber("00. Version "); if(loadVersion <= version) showMessage("Warning: The file you have selected has been analyzed with version "+loadVersion+", while you are currently running version "+version+".\nUnexpected things may happen."); nucleiChannel = Property.getNumber("02. Nuclei channel "); cytoChannel = Property.getNumber("03. Cytoplasm/membrane channel (-1 if not used) "); fociChannelA = Property.getNumber("04. Foci channel A "); fociChannelB = Property.getNumber("05. Foci channel B "); detectFociChannelB = Property.getNumber("06. Also detect foci channel B? "); ThreeDHandling = Property.get("07. 3D handling? "); singleSlice = Property.getNumber("08. [single z-slice foci detection only] Slice nr "); cropBorder = Property.getNumber("09. Remove image borders (pixels) "); XYBinning = Property.getNumber("10. Image XY binning before analysis [1-n] "); nucleiSegmentationChoice = Property.get("11. Nuclei/cell segmentation method "); downsampleFactorStarDist = Property.getNumber("12. Stardist nuclei rescaling factor [1-n], 0 for automatic, 1 for no rescaling "); probabilityThreshold = Property.getNumber("13. Probablility/flow threshold [0.0-1.0] (StarDist/Cellpose) "); CellposeModel = Property.getNumber("14. Cellpose model "); CellposeDiameter = Property.getNumber("15. Cellpose cell diameter (pixels), 0 for automatic "); minNucleusSize_setting = Property.getNumber("16. Remove nulei/cells with diameter smaller than (µm) "); maxNucleusSize_setting = Property.getNumber("17. Remove nulei/cells with diameter larger than (µm) "); manualRemoveNuclei = Property.get("18. Manually remove segmentated nuclei/cells"); excludeOnEdges = Property.getNumber("19. Exclude nulei on image edges "); //optimizationMode = Property.getNumber("20. Enable foci parameters optimization mode? "); fociSizeA = Property.get("21. Foci size channel A "); fociSizeB = Property.get("22. Foci size channel B "); fociDetectionMethod = Property.get("23. Foci detection method "); thresholdFactorA = Property.getNumber("24. Foci intensity threshold bias channel A "); thresholdFactorB = Property.getNumber("25. Foci intensity threshold bias channel B "); //thresholdA = Property.getNumber("26. Actual threshold value channel A "); //thresholdB = Property.getNumber("27. Actual threshold value channel B "); minFociSize = Property.getNumber("28. Minimum foci size (area)"); maxFociSize = Property.getNumber("29. Maximum foci size (area)"); maxFociDistanceOutsideNuclei_setting = Property.getNumber("30. Max distance of foci outside nuclei/cells (µm); -1 for full image"); minOverlapSize = Property.getNumber("31. Minimum overlap of foci to colocalize (pixels/voxels) "); overlayChoice = Property.get("32. Nuclei overlay choice "); fontColorCells = Property.get("33. Nuclei label color "); fontColorFoci = Property.get("34. Nuclei label color "); labelFontSize = Property.getNumber("35. Nuclei label size "); overlayBrightness = Property.get("36. Nuclei/cell outline bightness "); outlineColor = Property.get("37. Nuclei outline color "); close(); setBatchMode(false); } //nuclei detection nucleiBlurRadiusXY = 2; nucleiBlurRadiusZ = 2; nucleiMedian3DradiusXY = 2; nucleiMedian3DradiusZ = 2; //minNucleusSize = 4; //maxNucleusSize = 40; excludeOnEdges = excludeOnEdges; StarDistDownsampleInterpolation = "Bilinear"; //None, Bilinear or Bicubic if(nucleiSegmentationChoice == "Cellpose segmentation 2D (or on 2D projection)" || nucleiSegmentationChoice=="Cellpose segmentation 3D") cellpose = true; else cellpose=false; //foci detection var fociFilterSizeXY; //foci radius for DifferenceOfGaussians var fociFilterSizeZ; var fociSizeXY; //foci radius for detectMaximaBox var fociSizeZ; maxTileSize = 2048; // Maximum StarDist tile size thresholdMultiplier = 3; //Default threshold multiplier - foci should be this number of times more intense than the standard deviation of the background (average of all nuclei) minFociSize = minFociSize; maxFociSize = maxFociSize; var flipSlicesAndChannels = false; //H_max_peak_flooding = 100; //Percentage of the peak thickOutlines = true; //Width of nuclei outlines, 1 or 2 pixels thickOutlinesThresholdSize = 1200; //Above this image size nuclei outlines are always thick (2 pixels) useLargeSpots = false; //Make spots in overlay larger LABELOPACITY = 100; //opacity of outlines addNumbersOverlay = true; saveOverlayImage = true; var hideCellpose3DDialog = false; var optimizationMode = optimizationModeSetting; var processchoice; var doneOptimizing; var processAllOtherImages = false; var processAllOtherImagesFixedThreshold = false; var useROI = false; var reAnalyzeFullImage = false; var currentFileNr = 0; var displayMode = "composite"; //Default display settings at optimization phase var activeChannels = "11111"; var maxDisplaySetting; var processTime = 0; //var threshold = 0; var thresholdA = 0; var thresholdB = 0; //Create the azure and orange LUTs for the nuclei and cells (TO DO: maybe do this in a more elegant way than globals) b_reds = newArray(256); b_greens = newArray(256); b_blues = newArray(256); create_azure_lut(b_reds, b_greens, b_blues); o_reds = newArray(256); o_greens = newArray(256); o_blues = newArray(256); create_orange_lut(o_reds, o_greens, o_blues); saveSettings(); run("Set Measurements...", "area mean standard integrated median redirect=None decimal=3"); run("Conversions...", " "); setOption("BlackBackground", true); setForegroundColor(255, 255, 255); run("Colors...", "foreground=white background=black selection=cyan"); roiManager("Reset"); run("Clear Results"); print("\\Clear"); close("\\Others"); run("Clear Results"); run("CLIJ2 Macro Extensions", "cl_device="); Ext.CLIJ2_clear(); // Start of workflow //------------------------------ setBatchMode(true); run("Close All"); nrOfImages = files.length; if(!File.exists(outputFolder)) { createOutputFolder = getBoolean("The output folder '"+outputFolder+"' does not exist. Create?", "Of course, go ahead!", "See if I care!"); if(createOutputFolder) File.makeDirectory(outputFolder); else { formatHardDrive = getBoolean("Allright, try this:\nFormat the hard drive?", "Yes, goodbye forever", "No! Please mr. Foci Analyzer, I'll do anything you ask!"); if(formatHardDrive) { showMessage("Ok, you wished for it!"); exit("Oh, wait. I'll erase myself as well.\nCall it your lucky day then!"); } else exit("That's more like it. Now, run the macro and try again."); } } if (nucleiSegmentationChoice == "Load ROIs from file") { // ROIsFolder = call("ij.Prefs.get", "ROIs.Folder", File.getParent(files[0])); ROIsFolder = File.getParent(files[0]); Dialog.createNonBlocking("Select a ROIs folder"); Dialog.addMessage("ROI files should have the same name as the input images without extensions, followed by '_ROIs.zip'."); Dialog.addDirectory("Folder containing ROI .zip files", ROIsFolder); Dialog.show(); ROIsFolder = Dialog.getString(); call("ij.Prefs.set", "ROIs.Folder", ROIsFolder); print("Getting segmentations from ROIs in "+ROIsFolder); } if (nucleiSegmentationChoice == "Load label images") { labelImageFolder = call("ij.Prefs.get", "label.Image.Folder", File.getParent(files[0])); Dialog.createNonBlocking("Select a label image folder"); Dialog.addMessage("Label image files should have the exact same name as the input images."); Dialog.addDirectory("Folder containing label image files", labelImageFolder); Dialog.show(); labelImageFolder = Dialog.getString(); call("ij.Prefs.set", "label.Image.Folder", labelImageFolder); print("Getting segmentations from label images in "+labelImageFolder); } if (cellpose == true && CellposeModel == "custom") { CellposeModelPath = call("ij.Prefs.get", "Cellpose.custom.model.path", File.getParent(files[0])); CellposeDiameterPref = call("ij.Prefs.get", "Cellpose.diameter", CellposeDiameter); Dialog.createNonBlocking("Select a custom Cellpose model"); Dialog.addFile("Select or drag&drop custom Cellpose model", CellposeModelPath, 80); if(CellposeDiameter==0) Dialog.addNumber("Automatic diameter estimation does not work with custom models. Enter a cell diameter", CellposeDiameterPref, 0, 5, "pixels"); if(cytoChannel>0 && nucleiChannel>0) Dialog.addMessage("N.B. The custom model should be trained on images with a cytoplasm/membrane channel *and* a nucleus channel!\nIf not, no cells will be detected, but Fiji will not crash."); else if(cytoChannel>0 && nucleiChannel<=0) Dialog.addMessage("N.B. The custom model should be trained on images with a cytoplasm/membrane channel *only*!\nIf not, no cells will be detected, but Fiji will not crash."); else if(cytoChannel<=0 && nucleiChannel>0) Dialog.addMessage("N.B. The custom model should be trained on images with a nucleus channel *only*!\nIf not, no cells will be detected, but Fiji will not crash."); Dialog.show(); CellposeModelPath = Dialog.getString(); if(CellposeDiameter==0) CellposeDiameter = Dialog.getNumber(); CellposeModel = CellposeModelPath; //Both need to be the same call("ij.Prefs.set", "Cellpose.custom.model.path", CellposeModelPath); call("ij.Prefs.set", "Cellpose.diameter", CellposeDiameter); print("Using Cellpose model "+CellposeModelPath); } else CellposeModelPath = "path\\to\\own_cellpose_model"; //dummy name for (currentFileNr = 0; currentFileNr < nrOfImages; currentFileNr++) { print("\nProcessing file "+currentFileNr+1+"/"+nrOfImages+": "+files[currentFileNr] + "\n"); processFile(currentFileNr, files[currentFileNr], outputFolder); } close("Results"); logWindowContents = getInfo("log"); File.saveString(logWindowContents, outputFolder + File.separator + "Log.txt"); print("\n-------------------------------------------------------------------------"); print("Finished processing "+nrOfImages+" files in "+processTime*60+" seconds ("+d2s(processTime,1)+" minutes)."); print("Average speed: "+d2s((nrOfImages)/processTime,1)+" images per minute."); print("\\Update3:Finished processing "+nrOfImages+" files in "+processTime*60+" seconds ("+d2s(processTime,1)+" minutes)."); restoreSettings(); function processFile(current_image_nr, file, outputFolder) { run("Close All"); startTimeSeries = getTime(); print("\\Update1:Processing file "+current_image_nr+1+"/"+nrOfImages+": " + file); print("\\Update2:Average speed: "+d2s((current_image_nr)/processTime,1)+" images per minute."); time_to_run = (nrOfImages-(current_image_nr)) * processTime/(current_image_nr); if(time_to_run<5) print("\\Update3:Projected run time: "+d2s(time_to_run*60,0)+" seconds ("+d2s(time_to_run,1)+" minutes)."); else if(time_to_run<60) print("\\Update3:Projected run time: "+d2s(time_to_run,1)+" minutes. You'd better get some coffee."); else if(time_to_run<480) print("\\Update3:Projected run time: "+d2s(time_to_run,1)+" minutes ("+d2s(time_to_run/60,1)+" hours). You'd better go and do something useful."); else if(time_to_run<1440) print("\\Update3:Projected run time: "+d2s(time_to_run,1)+" minutes. ("+d2s(time_to_run/60,1)+" hours). You'd better come back tomorrow."); else if(time_to_run>1440) print("\\Update3:Projected run time: "+d2s(time_to_run,1)+" minutes. This is never going to work. Give it up!"); print("\\Update4:-------------------------------------------------------------------------"); run("Bio-Formats Macro Extensions"); //Necessary to do this here, because you can only activate one Macro Extension at the time Ext.setId(file); Ext.getSeriesCount(nr_series); run("CLIJ2 Macro Extensions", "cl_device="); //Necessary to do this here, because you can only activate one Macro Extension at the time Ext.CLIJ2_clear(); filename = File.getName(file); fileExtension = substring(filename, lastIndexOf(filename, "."), lengthOf(filename)); if(endsWith(fileExtension, "tif") || endsWith(fileExtension, "jpg") || endsWith(fileExtension, "png")) { //Use standard opener open(file); if(XYBinning > 1) run("Bin...", "x="+XYBinning+" y="+XYBinning+" z=1 bin=Average"); if(reAnalyzeFullImage == false && optimizationModeSetting == true && processAllOtherImages == false) optimizationMode = true; process_current_series(file, true); if(useROI) { reAnalyzeFullImage = true; currentFileNr--; //ugly change of global variable, but I can't see another easy way if(optimizationModeSetting == true && processAllOtherImages == false) optimizationMode = false; } else reAnalyzeFullImage = false; } else { //Use Bio-Formats for(s = 0; s < nr_series; s++) { run("Close All"); run("Bio-Formats Importer", "open=["+file+"] autoscale color_mode=Default view=Hyperstack stack_order=XYCZT series_"+s+1); seriesName = getTitle(); seriesName = replace(seriesName,"\\/","-"); //replace slashes by dashes in the seriesName // print(seriesName); // outputPath = output + File.separator + substring(seriesName) if(XYBinning > 1) run("Bin...", "x="+XYBinning+" y="+XYBinning+" z=1 bin=Average"); process_current_series(seriesName, false); if(useROI) s--; //Analyze this series again without ROI //may work, but only for files opened by bioformats } } } function process_current_series(image, nameHasExtension) { if(ThreeDHandling == "Process a single z-slice only (specify below which slice)") run("Duplicate...", "duplicate slices="+singleSlice); if(nameHasExtension) imageName = File.getNameWithoutExtension(image); else imageName = image; //Initialize image and table getDimensions(gwidth, gheight, gchannels, gslices, gframes); // global variables getVoxelSize(pixelWidth, pixelHeight, pixelDepth, unit); if(bitDepth() == 24) { print("[WARNING] The input image is RGB. It will be converted to composite multichannel image with channels 1:Red, 2:Green; 3:Blue."); run("Make Composite"); } if(gslices>1) imageIs3D = true; else imageIs3D = false; //timelapse handling - save individual frames and call processFile recursively [smiley with sunglasses] if(gframes > 1){ original = getTitle; for(t=1; t<=gframes; t++) { showStatus("Splitting time frames... "+t+"/"+gframes); showProgress(t, gframes); run("Duplicate...", "duplicate frames="+t); saveAs("tiff", outputFolder + File.separator + imageName + "__t="+IJ.pad(t, 3)); close(); } if (nucleiSegmentationChoice == "Load label images") { open(labelImageFolder + File.separator + File.getNameWithoutExtension(image) + ".tif"); for(t=1; t<=gframes; t++) { showStatus("Splitting label image time frames... "+t+"/"+gframes); showProgress(t, gframes); run("Duplicate...", "duplicate range="+t+"-"+t); saveAs("tiff", labelImageFolder + File.separator + imageName + "__t="+IJ.pad(t, 3)); close(); } } for(t=1; t<=gframes; t++) { processFile(t-1, outputFolder + File.separator + imageName + "__t="+IJ.pad(t, 3)+".tif", outputFolder); } //overlay_image = "Foci_overlay_ch3"; // for(t=1; t<=gframes; t++) { // open(outputFolder + File.separator + imageName + "__t="+IJ.pad(t, 3) + "__" + overlay_image + ".zip"); //// run("To ROI Manager"); //// Overlay.copy; // if(t>1) run("Concatenate...", "image1=all_frames image2=" + imageName + "__t="+IJ.pad(t, 3) + "__" + overlay_image + ".tif"); // rename("all_frames"); // } //setBatchMode("exit and display"); return; } original = getTitle(); print("Image path: "+image+"\n"); if(gchannels == 1 && gslices > 1) { if(flipSlicesAndChannels == false) { setBatchMode("show"); Dialog.createNonBlocking("Single channel detected"); Dialog.addMessage("Warning: ["+imageName+"] has only 1 Channel, but "+gslices+" Slices."); Dialog.addRadioButtonGroup("Do you want to flip Channels and Slices?", newArray("Yes, and do this for all subsequent images", "Yes, only for this image", "No (exit macro)"), 3, 1, "Yes, and do this for all subsequent images"); Dialog.show(); answer = Dialog.getRadioButton(); if(answer == "Yes, and do this for all subsequent images") flipSlicesAndChannels = true; if(answer == "Yes, only for this image") run("Re-order Hyperstack ...", "channels=[Slices (z)] slices=[Channels (c)] frames=[Frames (t)]"); else if(answer == "No (exit macro)") exit(); } if(flipSlicesAndChannels == true) run("Re-order Hyperstack ...", "channels=[Slices (z)] slices=[Channels (c)] frames=[Frames (t)]"); setBatchMode("hide"); } Stack.setDisplayMode("grayscale"); getDimensions(gwidth, gheight, gchannels, gslices, gframes); Stack.setChannel(nucleiChannel); if(cropBorder>0) { makeRectangle(cropBorder, cropBorder, gwidth-2*cropBorder, gheight-2*cropBorder); run("Crop"); getDimensions(gwidth, gheight, gchannels, gslices, gframes); // global variables } if(gwidth > thickOutlinesThresholdSize && gheight > thickOutlinesThresholdSize) { thickOutlines = true; } if(gslices > 1) anisotropyFactor = pixelDepth / pixelWidth; else anisotropyFactor = 0; bits = bitDepth(); //TO DO: FIND A BETTER WAY FOR THIS: if(bits == 8) run("16-bit"); //Convert to 16-bit, because 8-bit restricts the foci labelmap to 255 foci Stack.setSlice(gslices/2); setBatchMode("show"); run("Enhance Contrast", "saturated=0.35"); resultTable = "Foci results per cell"; if(isOpen(resultTable)) Table.reset(resultTable); else { Table.create(resultTable); Table.setLocationAndSize(0, 0, 1000, 500); } allFociResultsTable = "All foci statistics"; if(isOpen(allFociResultsTable)) close(allFociResultsTable); Table.create(allFociResultsTable); Table.setLocationAndSize(0, 500, 1000, 500); if(gslices > 1 && ThreeDHandling == "Detect foci on the Maximum Intensity Projection") { selectWindow(original); run("Z Project...", "projection=[Max Intensity]"); setBatchMode("show"); original = getTitle(); gslices = 1; } if(gslices > 1 && ThreeDHandling == "Detect foci on the Standard Deviation Projection") { selectWindow(original); run("Z Project...", "projection=[Standard Deviation]"); setBatchMode("show"); original = getTitle(); gslices = 1; } if(gslices > 1 && ThreeDHandling == "Detect foci on the Summed Intensity Projection") { selectWindow(original); run("Z Project...", "projection=[Sum Slices]"); setBatchMode("show"); original = getTitle(); gslices = 1; run("Conversions...", "scale"); } if(gslices > 1 && ThreeDHandling == "Detect foci on a Extended Depth of Focus Projection") { //TO DO: Doesn't work yet, because multichannel is not compatible with CLIJ2 selectWindow(original); getDimensions(width, height, channels, slices, frames); mergeString = ""; for(c=1; c<=channels; c++) { selectWindow(original); Stack.setChannel(c); Ext.CLIJ2_push(original); EDF = "EDF"; Ext.CLIJ2_extendedDepthOfFocusSobelProjection(original, EDF, 10); Ext.CLIJ2_pull(EDF); rename("EDF_channel_"+c); mergeString += " c"+c+"=EDF_channel_"+c; } Ext.CLIJ2_release(EDF); run("Merge Channels...", mergeString+" create"); rename(original+"_EDF"); Stack.setDisplayMode("grayscale"); original = getTitle(); setBatchMode("show"); gslices = 1; } if(gslices == 1) pixelDepth = 1; //Run analysis only on a ROI wait(50); x_ROI = 0; y_ROI = 0; useROI = false; if(isKeyDown("alt")) { setTool("rectangle"); waitForUser("Select a (rectangular) ROI to optimize the analysis on"); if(selectionType==0) { useROI = true; getSelectionBounds(x_ROI, y_ROI, width_ROI, height_ROI); run("Duplicate...", "duplicate title=selection"); process_image = getTitle(); } else { print("Ignoring non-rectangular selection..."); run("Select None"); } } if(useROI == false) process_image = original; //Segment and label nuclei if (nucleiSegmentationChoice == "StarDist nuclei segmentation 2D (or on 2D projection)") nuclei_info = segmentNucleiStarDist(process_image, nucleiChannel, probabilityThreshold, pixelWidth, unit, resultTable); else if (nucleiSegmentationChoice == "Classic segmentation") nuclei_info = segmentNucleiClassic(process_image, nucleiChannel, nucleiBlurRadiusXY, nucleiBlurRadiusZ, nucleiMedian3DradiusXY, nucleiMedian3DradiusZ); else if (cellpose == true) nuclei_info = segmentCellsCellpose(process_image, nucleiChannel, cytoChannel, probabilityThreshold, pixelWidth, pixelDepth, unit, resultTable); else if (nucleiSegmentationChoice == "Load ROIs from file") nuclei_info = loadROIs(original, nucleiChannel, ROIsFolder); else if (nucleiSegmentationChoice == "Load label images") nuclei_info = loadROIs(original, nucleiChannel, labelImageFolder); if(nuclei_info[0] == "FileNotFound") continue; labelmap_nuclei = nuclei_info[0]; nuclei_outlines = nuclei_info[1]; nrNuclei = nuclei_info[2]; if(manualRemoveNuclei != "No thanks" && nrNuclei > 0) { nuclei_info = manually_remove_labels(labelmap_nuclei, nuclei_outlines, nrNuclei, process_image, imageName); labelmap_nuclei = nuclei_info[0]; nuclei_outlines = nuclei_info[1]; nrNuclei = nuclei_info[2]; labelmap_nuclei = "Labelmap_nuclei_edited"; } if(nrNuclei > 0) { //Write ID and area to foci results table run("Clear Results"); Ext.CLIJ2_statisticsOfLabelledPixels(labelmap_nuclei, labelmap_nuclei); nucleus_id_ = Table.getColumn("IDENTIFIER", "Results"); nucleus_area_ = Table.getColumn("PIXEL_COUNT", "Results"); nucleus_area_ = multiplyArraywithScalar(nucleus_area_, Math.sqr(pixelWidth)); Table.setColumn("Cell ID", nucleus_id_, resultTable); if (nucleiSegmentationChoice == "Load ROIs from file") { for (i = 0; i < nrNuclei; i++) { roiManager("select", i); Table.set("Cell UUID", i, Roi.getName, resultTable); } roiManager("deselect"); } Table.setColumn("Cell area 2D ("+unit+"^2)", nucleus_area_, resultTable); //Create a 3D version of the 2D nuclei labelmap, if required Ext.CLIJ2_getDimensions(labelmap_nuclei, labelmap_width, labelmap_height, labelmap_depth); if(gslices>1 && labelmap_depth==1) Ext.CLIJ2_imageToStack(labelmap_nuclei, labelmap_nuclei_3D, gslices); else labelmap_nuclei_3D = labelmap_nuclei; //Foci filtering and detection - in a loop to enable parameter optimization firstTimeProcessing = true; zoom = 1; //Add selected ROI to the ROI Manager and add as overlay - must do this after StarDist if(useROI) { selectWindow(original); roiManager("reset"); roiManager("add"); roiManager("Select",0); roiManager("Set Color", "white"); roiManager("Set Line Width", 1); roiManager("add"); } do { if(processAllOtherImages == false) doneOptimizing = false; //Reset this parameter from the previous round //Create dilated labelmap for including foci outside the nuclei labelmap_nuclei_3D_dilated = "labelmap_nuclei_3D_dilated"; if(maxFociDistanceOutsideNuclei_setting != 0) { maxFociDistanceOutsideNuclei = round(maxFociDistanceOutsideNuclei_setting / pixelWidth); if(maxFociDistanceOutsideNuclei_setting < 0) Ext.CLIJ2_dilateLabels(labelmap_nuclei_3D, labelmap_nuclei_3D_dilated, maxOf(gwidth, gheight)); else Ext.CLIJ2_dilateLabels(labelmap_nuclei_3D, labelmap_nuclei_3D_dilated, maxFociDistanceOutsideNuclei); if(ThreeDHandling == "Detect foci in 3D (or 2D when N/A)" && nucleiSegmentationChoice == "Cellpose segmentation 3D" && pixelDepth/pixelWidth > 1.33) print("[WARNING] Due to non-isotropic pixels the foci region will be expanded from the nuclei more in Z than in X and Y! (by a factor of "+d2s(pixelDepth/pixelWidth,1)+")"); //N.B. For non-isotropic 3D data this 3D dilation is not fair, but making the labelmap isotropic creates intermediate (non-integer) values. Oh well.. } else Ext.CLIJ2_copy(labelmap_nuclei_3D, labelmap_nuclei_3D_dilated); if(debugMode) showImagefromGPU(labelmap_nuclei_3D_dilated); //Create outlines from dilated nuclei labelmap nuclei_dilated_outlines = "nuclei_dilated_outlines"; if(isOpen("nuclei_dilated_outlines")) close(nuclei_dilated_outlines); if(maxFociDistanceOutsideNuclei_setting != 0 && nucleiSegmentationChoice != "Cellpose segmentation 3D") { Ext.CLIJ2_copySlice(labelmap_nuclei_3D_dilated, labelmap_nuclei_2D_dilated, 0); Ext.CLIJ2_detectLabelEdges(labelmap_nuclei_2D_dilated, nuclei_dilated_edges); Ext.CLIJ2_mask(labelmap_nuclei_2D_dilated, nuclei_dilated_edges, nuclei_dilated_outlines); Ext.CLIJ2_release(nuclei_dilated_edges); Ext.CLIJ2_release(labelmap_nuclei_2D_dilated); Ext.CLIJ2_pullBinary(nuclei_dilated_outlines); if(outlineColor == "White") outlineColor = "Grays"; run(outlineColor); } else if(maxFociDistanceOutsideNuclei_setting != 0 && nucleiSegmentationChoice == "Cellpose segmentation 3D") { Ext.CLIJ2_maximumSliceBySliceSphere(labelmap_nuclei_3D_dilated, labelmap_maximum, 1, 1); Ext.CLIJ2_minimumSliceBySliceSphere(labelmap_nuclei_3D_dilated, labelmap_minimum, 1, 1); Ext.CLIJ2_subtractImages(labelmap_maximum, labelmap_minimum, label_edges); Ext.CLIJ2_release(labelmap_maximum); Ext.CLIJ2_release(labelmap_minimum); Ext.CLIJ2_threshold(label_edges, nuclei_dilated_outlines, 1); Ext.CLIJ2_pull(nuclei_dilated_outlines); Ext.CLIJ2_release(label_edges); Ext.CLIJ2_release(nuclei_dilated_outlines); if(outlineColor == "White") outlineColor = "Grays"; run(outlineColor); } //Check if fociChannelB actually exists and is different from Foci channel A if(detectFociChannelB == true && fociChannelB > gchannels) exit("Error: The selected Foci channel B ("+fociChannelB+") is higher than the number of channels of the image ("+gchannels+")."); if(detectFociChannelB == true && fociChannelB == fociChannelA) { showMessageWithCancel("Selected Foci channels are the same", "Warning: The selected Foci channel B ("+fociChannelB+") is the same as Foci channel A!\nIf you continue only this channel will be analyzed."); detectFociChannelB = false; } //Foci channel A detections_fociA = detect_foci(process_image, fociChannelA, fociSizeA, anisotropyFactor, firstTimeProcessing, labelmap_nuclei, labelmap_nuclei_3D, thresholdFactorA, thresholdA); labelmap_fociA = detections_fociA[0]; mask_fociA = detections_fociA[1]; spots_fociA = detections_fociA[2]; thresholdA = detections_fociA[3]; if(firstTimeProcessing == false) call("ij.gui.ImageWindow.setNextLocation", x_image, y_image); overlayA = mergeOriginalAndDetection(original, nrNuclei, nuclei_outlines, mask_fociA, spots_fociA, fociChannelA, useROI, x_ROI, y_ROI); overlay_image = overlayA; //Will be overwritten if channel B is also used if(firstTimeProcessing == false && detectFociChannelB == false) close("Processing..."); //Foci channel B if(detectFociChannelB == true) { detections_fociB = detect_foci(process_image, fociChannelB, fociSizeB, anisotropyFactor, firstTimeProcessing, labelmap_nuclei, labelmap_nuclei_3D, thresholdFactorB, thresholdB); labelmap_fociB = detections_fociB[0]; mask_fociB = detections_fociB[1]; spots_fociB = detections_fociB[2]; thresholdB = detections_fociB[3]; overlayB = mergeOriginalAndDetection(original, nrNuclei, nuclei_outlines, mask_fociB, spots_fociB, fociChannelB, useROI, x_ROI, y_ROI); if(firstTimeProcessing == false) close("Processing..."); Overlay.copy(); //Preserve nuclei outline overlays overlay_image = "Foci_overlay_ch"+fociChannelA+"_and_ch"+fociChannelB; run("Concatenate...", " title="+overlay_image+" image1="+overlayA+" image2="+overlayB+" image3=[-- None --]"); run("Stack to Hyperstack...", "order=xyczt(default) channels="+nSlices/(gslices*2)+" slices="+gslices+" frames=2 display=Composite"); Overlay.paste(); Stack.setFrame(2); run("Label...", "format=Text starting=0 interval=0 x=20 y=10 font=18 text=[Channel "+fociChannelB+"] range=[2-2] use"); Stack.setFrame(1); run("Label...", "format=Text starting=0 interval=0 x=20 y=10 font=18 text=[Channel "+fociChannelA+"] range=[1-1] use"); } //Set display settings and location if(firstTimeProcessing == false) { Stack.setChannel(1); setMinAndMax(minDisplayNuclei, maxDisplayNuclei); Stack.setChannel(2); setMinAndMax(minDisplayFoci, maxDisplayFoci); if(cytoChannel>0) { Stack.setChannel(5); setMinAndMax(minDisplayCells, maxDisplayCells); } // Stack.setDisplayMode(displayMode); // Stack.setActiveChannels(activeChannels); } // Property.set("CompositeProjection", "Max"); //Use 'Composite Max' display setting Stack.setDisplayMode(displayMode); Stack.setActiveChannels(activeChannels); if(firstTimeProcessing == false) Stack.setPosition(currentChannel, currentSlice, currentFrame); else if(gslices > 1) Stack.setSlice(gslices/2); setBatchMode("show"); if(firstTimeProcessing == false) run("Set... ", "zoom="+zoom*100+" x="+displayX + displayWidth/2+" y="+displayY + displayHeight/2); run("Channels Tool..."); selectWindow("Channels"); setLocation(0, 0); run("Brightness/Contrast..."); selectWindow("B&C"); setLocation(-5, 210); if(firstTimeProcessing == true) getLocationAndSize(x_image, y_image, imageWidth, imageHeight); //if(slices>1) run("Animation Options...", "speed=4 loop start"); //Parameter optimization dialog if(optimizationMode == true && processAllOtherImages == false) { Dialog.createNonBlocking("Optimize settings for foci detection"); Dialog.addMessage("Inspect the detected foci and optimize the settings. Some tips:\n• Hide and show the foci with the checkbox 'channel 3' in the Channels Tool (upper left corner of the screen).\n• Change the display brighness of the 'channels' (1:nuclei, 2:foci, 3:detected foci, 4:foci centers, [5:cells]) in the B&C window.\n• Zoom [+/- or up/down keys] and pan [hold space & drag the image].\n• Use the sliders below the image to change the active channel [arrow keys] and displayed z-slice [Ctrl + arrow keys],\n and to switch between Foci channels A and B [Alt + arrow keys].\n\n ", 12, "#000080"); Dialog.addChoice("Foci size channel "+fociChannelA, newArray("tiny","small","average","large","huge","other"), fociSizeA); if(detectFociChannelB) Dialog.addChoice("Foci size channel "+fociChannelB, newArray("tiny","small","average","large","huge","other"), fociSizeB); Dialog.addChoice("Detection method", newArray("Marker-controlled watershed (recommended)","AreaMaxima local maximum detection"), fociDetectionMethod); Dialog.addSlider("Threshold bias channel "+fociChannelA+" (higher is more strict)", -2.5, 2.5, thresholdFactorA); if(detectFociChannelB) Dialog.addSlider("Threshold bias channel "+fociChannelB+" (higher is more strict)", -2.5, 2.5, thresholdFactorB); if(gslices > 1) { Dialog.addNumber("Minimum foci size", minFociSize, 1, 4, "voxels"); Dialog.addNumber("Maximum foci size", maxFociSize, 1, 4, "voxels"); } else { Dialog.addNumber("Minimum foci size", minFociSize, 1, 4, "pixels"); Dialog.addNumber("Maximum foci size", maxFociSize, 1, 4, "pixels"); } Dialog.addNumber("Extend foci detection outside nuclei/cells with", maxFociDistanceOutsideNuclei_setting, 1, 4, unit+" (-1 for full image)") if(detectFociChannelB) { if(gslices > 1) Dialog.addNumber("Minimum overlap of foci", minOverlapSize, 1, 4, "voxels"); else Dialog.addNumber("Minimum overlap of foci", minOverlapSize, 1, 4, "pixels"); } Dialog.addChoice("Nuclei outline visualization", newArray("bright","dim"), overlayBrightness); items = newArray("Recalculate with these settings", "Done optimizing | Process and optimize the next image", "Done optimizing | Process all other images with these settings (calculate thresholds for each image)", "Done optimizing | Process all other images with these settings (fix current threshold levels)"); Dialog.addChoice("Action", items, "Recalculate with these settings"); Dialog.addMessage("Click Help for more info (Foci Analyzer ImageJ site)", 12, "#000080"); Dialog.addHelp("https://imagej.net/plugins/foci-analyzer"); //determine dialog location if(x_image + imageWidth + 500 < screenWidth) Dialog.setLocation(x_image+imageWidth-15, y_image); else if(x_image > 500) Dialog.setLocation(x_image-530, y_image); // else if(y_image + gheight + 300 < screenHeight) Dialog.setLocation(x_image, y_image+gheight+100); else Dialog.setLocation(x_image, y_image); Dialog.show(); //Get dialog entries fociSizeA = Dialog.getChoice(); if(detectFociChannelB == true) fociSizeB = Dialog.getChoice(); fociDetectionMethod = Dialog.getChoice(); thresholdFactorA = Dialog.getNumber(); if(detectFociChannelB) thresholdFactorB = Dialog.getNumber(); minFociSize = Dialog.getNumber(); maxFociSize = Dialog.getNumber(); maxFociDistanceOutsideNuclei_setting = Dialog.getNumber(); if(detectFociChannelB) minOverlapSize = Dialog.getNumber(); overlayBrightness = Dialog.getChoice(); processChoice = Dialog.getChoice(); if(processChoice == "Done optimizing | Process and optimize the next image") doneOptimizing = true; processAllOtherImages = false; if(processChoice == "Done optimizing | Process all other images with these settings (calculate thresholds for each image)") { doneOptimizing = true; processAllOtherImages = true; } if(processChoice == "Done optimizing | Process all other images with these settings (fix current threshold levels)") { doneOptimizing = true; processAllOtherImages = true; processAllOtherImagesFixedThreshold = true; } //Write persistent Script parameters for the next time the macro is run. setScriptParameterValue("fociSizeA", fociSizeA); setScriptParameterValue("fociSizeB", fociSizeB); setScriptParameterValue("thresholdFactorA", thresholdFactorA); setScriptParameterValue("thresholdFactorB", thresholdFactorB); setScriptParameterValue("fociDetectionMethod",fociDetectionMethod); setScriptParameterValue("minFociSize",minFociSize); setScriptParameterValue("maxFociSize",maxFociSize); setScriptParameterValue("maxFociDistanceOutsideNuclei_setting",maxFociDistanceOutsideNuclei_setting); setScriptParameterValue("minOverlapSize",minOverlapSize); setScriptParameterValue("overlayBrightness",overlayBrightness); //get image display properties selectWindow(overlay_image); getLocationAndSize(x_image, y_image, imageWidth, imageHeight); Stack.getActiveChannels(activeChannels); Stack.getDisplayMode(displayMode); Stack.getPosition(currentChannel, currentSlice, currentFrame); Stack.setChannel(1); getMinAndMax(minDisplayNuclei, maxDisplayNuclei); Stack.setChannel(2); getMinAndMax(minDisplayFoci, maxDisplayFoci); if(cytoChannel>0) { Stack.setChannel(5); getMinAndMax(minDisplayCells, maxDisplayCells); } zoom = getZoom(); getDisplayedArea(displayX, displayY, displayWidth, displayHeight); if(doneOptimizing == false) { close("foci_ch"+fociChannelA); //close("MAX_Foci_ch"+channel+"_filtered_and_masked"); close("Labelmap_detected_foci_filtered_ch"+fociChannelA); close(mask_fociA); // close("foci_RAW"+fociChannelA); if(detectFociChannelB == false) { selectWindow("Foci_overlay_ch"+fociChannelA); rename("Processing..."); } else if(detectFociChannelB == true) { close("foci_ch"+fociChannelB); //close("MAX_Foci_ch"+channel+"_filtered_and_masked"); close("Labelmap_detected_foci_filtered_ch"+fociChannelB); close(mask_fociB); // close("foci_RAW_ch"+fociChannelB); selectWindow("Foci_overlay_ch"+fociChannelA+"_and_ch"+fociChannelB); rename("Processing..."); } } else useROI = false; // N.B. REMOVING THE PREVIOUS LINE WILL RESULT IN THE FULL IMAGE BEING OPENED AND ANALYZED! THIS HAPPENS ANYWAY WHEN OPTIMIZATION IS OFF. } firstTimeProcessing = false; } while(optimizationMode == true && doneOptimizing == false); if(useROI == false) { //Get nuclei positions if(maxFociDistanceOutsideNuclei_setting != 0) Ext.CLIJ2_statisticsOfLabelledPixels(labelmap_nuclei_3D_dilated, labelmap_nuclei_3D_dilated); selectWindow("Results"); boundingBox_X = Table.getColumn("BOUNDING_BOX_X"); boundingBox_Y = Table.getColumn("BOUNDING_BOX_Y"); boundingBox_Z = Table.getColumn("BOUNDING_BOX_Z"); boundingBox_width = Table.getColumn("BOUNDING_BOX_WIDTH"); boundingBox_height = Table.getColumn("BOUNDING_BOX_HEIGHT"); boundingBox_depth = Table.getColumn("BOUNDING_BOX_DEPTH"); Array.getStatistics(boundingBox_width, minWidth, maxWidth); Array.getStatistics(boundingBox_height, minHeight, maxHeight); Array.getStatistics(boundingBox_depth, minDepth, maxDepth); if(debugMode) print("\nMaximum nucleus bounding box: "+maxWidth+", "+maxHeight); //Measure the foci nrFoci = measureFoci(original, fociChannelA, nrNuclei, labelmap_nuclei_3D, labelmap_fociA, boundingBox_X, boundingBox_Y, boundingBox_Z, maxWidth, maxHeight, maxDepth); if(detectFociChannelB) nrFoci = measureFoci(original, fociChannelB, nrNuclei, labelmap_nuclei_3D, labelmap_fociB, boundingBox_X, boundingBox_Y, boundingBox_Z, maxWidth, maxHeight, maxDepth); if(gslices > 1 && nrFoci>0) Table.renameColumn("Volume", "Volume (voxels)", allFociResultsTable); else if (gslices == 1 && nrFoci>0) Table.renameColumn("Volume", "Area (pixels)", allFociResultsTable); Table.update(allFociResultsTable); Table.deleteRows(0, 0, "Results"); //remove background label overlay_numbers_on_image(overlay_image); // print("\n\nGPU Memory after channel "+fociChannelB); // Ext.CLIJ2_reportMemory(); //Compute A->B foci colocalization if(detectFociChannelB) { foci_overlap_map = computeOverlap("Mask_foci_ch" + fociChannelA, "Mask_foci_ch" + fociChannelB, nrNuclei, labelmap_nuclei_3D, boundingBox_X, boundingBox_Y, maxWidth, maxHeight); //These images are still open in RAM/GPU // print("\n\nGPU Memory after computing overlap"); // Ext.CLIJ2_reportMemory(); } //Measure intensity in non-foci channels measure_nuclear_intensities(original, nrNuclei, labelmap_nuclei_3D, gchannels, fociChannelA, fociChannelB, resultTable); //Save data //Create parameter list getDateAndTime(year, month, dayOfWeek, dayOfMonth, hour, minute, second, msec); MonthNames = newArray("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"); DayNames = newArray("Sun", "Mon","Tue","Wed","Thu","Fri","Sat"); List.set("00. ==Foci Analyzer settings=", ""); List.set("00. Date ", " " + DayNames[dayOfWeek] + " " + dayOfMonth + " " + MonthNames[month] + " " + year); List.set("00. Time ", " " + hour +":"+IJ.pad(minute,2)+":"+IJ.pad(second,2)); List.set("00. Version ", version); if(loadSettingsFromFile == true) List.set("01. Analysis settings were loaded from ", loadSettingsPath); else List.set("01. Settins were entered from dialog", 1); List.set("02. Nuclei channel ", nucleiChannel); List.set("03. Cytoplasm/membrane channel (-1 if not used) ", cytoChannel); List.set("04. Foci channel A ", fociChannelA); List.set("05. Foci channel B ", fociChannelB); List.set("06. Also detect foci channel B? ", detectFociChannelB); List.set("07. 3D handling? ", ThreeDHandling); List.set("08. [single z-slice foci detection only] Slice nr ", singleSlice); List.set("09. Remove image borders (pixels) ", cropBorder); List.set("10. Image XY binning before analysis [1-n] ", XYBinning); List.set("11. Nuclei/cell segmentation method ", nucleiSegmentationChoice); List.set("12. Stardist nuclei rescaling factor [1-n], 0 for automatic, 1 for no rescaling ", downsampleFactorStarDist); List.set("13. Probablility/flow threshold [0.0-1.0] (StarDist/Cellpose) ", probabilityThreshold); List.set("14. Cellpose model ", CellposeModel); List.set("15. Cellpose cell diameter (pixels), 0 for automatic ", CellposeDiameter); List.set("16. Remove nulei/cells with diameter smaller than (µm) ", minNucleusSize_setting); List.set("17. Remove nulei/cells with diameter larger than (µm) ", maxNucleusSize_setting); List.set("18. Manually remove segmentated nuclei/cells", manualRemoveNuclei); List.set("19. Exclude nulei on image edges ", excludeOnEdges); List.set("20. Enable foci parameters optimization mode? ", optimizationMode); List.set("21. Foci size channel A ", fociSizeA); List.set("22. Foci size channel B ", fociSizeB); List.set("23. Foci detection method ", fociDetectionMethod); List.set("24. Foci intensity threshold bias channel A ", thresholdFactorA); List.set("25. Foci intensity threshold bias channel B ", thresholdFactorB); List.set("26. Actual threshold value channel A ", thresholdA); List.set("27. Actual threshold value channel B ", thresholdB); List.set("28. Minimum foci size (area) ", minFociSize); List.set("29. Maximum foci size (area) ", maxFociSize); List.set("30. Max distance of foci outside nuclei/cells (µm); -1 for full image ", maxFociDistanceOutsideNuclei_setting); List.set("31. Minimum overlap of foci to colocalize (pixels/voxels) ", minOverlapSize); List.set("32. Nuclei overlay choice ", overlayChoice); List.set("33. Nuclei/cell label color ", fontColorCells); List.set("34. Nuclei/cell label color ", fontColorFoci); List.set("35. Nuclei/cell label size ", labelFontSize); List.set("36. Nuclei/cell outline bightness ", overlayBrightness); List.set("37. Nuclei/cell outline color ", outlineColor); list = List.getList(); selectWindow(resultTable); Table.save(outputFolder + File.separator + imageName + "__Foci_results.tsv"); selectWindow(allFociResultsTable); Table.save(outputFolder + File.separator + imageName + "__All_Foci_statistics.tsv"); if(saveOverlayImage == true) { if(detectFociChannelB == true) { selectImage(foci_overlap_map); Stack.setSlice(gslices/2); setBatchMode("show"); setMetadata("info", list); saveAs("zip", outputFolder + File.separator + imageName + "__foci_coloc_map"); } } selectImage(overlay_image); setBatchMode("hide"); setBatchMode("show"); //Lame trick to move this window to the front setMetadata("info", list); saveAs("zip", outputFolder + File.separator + imageName + "__" + overlay_image); overlay_image = getTitle(); //Name has changed run("Clear Results"); // Ext.CLIJ2_release(labelmap_nuclei); // Ext.CLIJ2_release(labelmap_nuclei_3D); // Ext.CLIJ2_release(labelmap_nuclei_3D_dilated); // Ext.CLIJ2_release("foci_ch"+fociChannelA); // if(detectFociChannelB == true) Ext.CLIJ2_release("foci_ch"+fociChannelB); // Ext.CLIJ2_release("Labelmap_detected_foci_ch"+fociChannelA); // if(detectFociChannelB == true) Ext.CLIJ2_release("Labelmap_detected_foci_ch"+fociChannelB); // Ext.CLIJ2_reportMemory(); Ext.CLIJ2_clear(); } } else print("[WARNING] 0 nuclei detected in "+ image +" !"); //Ext.CLIJx_organiseWindows(Number startX, Number startY, Number tilesX, Number tilesY, Number tileWidth, Number tileHeight); cleanup(); //Ext.CLIJ2_reportMemory(); endTime = getTime(); processTime = processTime+(endTime-startTimeSeries)/60000; //------------------------------ // End of workflow } function segmentNucleiClassic(image, channel, nucleiBlurRadiusXY, nucleiBlurRadiusZ, nucleiMedian3DradiusXY, nucleiMedian3DradiusZ) { selectWindow(image); getDimensions(width, height, channels, slices, frames); minNucleusSize = PI*Math.sqr((minNucleusSize_setting / pixelWidth / 2)); //Calculate the nucleus area as if it were a circle maxNucleusSize = PI*Math.sqr((maxNucleusSize_setting / pixelWidth / 2)); nuclei = "nuclei"; run("Duplicate...", "title=nuclei duplicate channels=&channel"); run("32-bit"); showStatus("Segmenting nuclei..."); Ext.CLIJ2_push(nuclei); Ext.CLIJ2_gaussianBlur3D(nuclei, nuclei_filtered1, nucleiBlurRadiusXY, nucleiBlurRadiusXY, nucleiBlurRadiusZ); //run("Median 3D...", "x=" + nucleiMedian3DradiusXY + " y=" + nucleiMedian3DradiusXY + " z=" + nucleiMedian3DradiusZ); // CLIJ2_Median3D gives problems, at least on my GPU // Ext.CLIJ2_median3DSphere(nuclei, nuclei_filtered1, nucleiMedian3DradiusXY/2, nucleiMedian3DradiusXY/2, nucleiMedian3DradiusZ/2); Ext.CLIJ2_maximumZProjection(nuclei_filtered1, nuclei_filtered_MAX); Ext.CLIJ2_automaticThreshold(nuclei_filtered_MAX, thresholded, "Otsu"); Ext.CLIJ2_pullBinary(thresholded); run("Watershed"); //Not on GPU because results are not so good. run("Properties...", "unit=&unit pixel_width=&pixelWidth pixel_height=&pixelHeight voxel_depth=1.0000"); if(excludeOnEdges) run("Analyze Particles...", "size=" + minNucleusSize + "-" + maxNucleusSize + " circularity=0.25-1.00 show=[Count Masks] exclude include add"); else run("Analyze Particles...", "size=" + minNucleusSize + "-" + maxNucleusSize + " pixel circularity=0.25-1.00 show=[Count Masks] include add"); rename("Labelmap_nuclei"); labelmap_nuclei = "Labelmap_nuclei"; Ext.CLIJ2_push(labelmap_nuclei); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei, nrNuclei); //Count the number of nuclei //Create nuclei outlines from nuclei labelmap Ext.CLIJ2_detectLabelEdges(labelmap_nuclei, labelmap_edges); if(thickOutlines == false) { Ext.CLIJ2_mask(labelmap_nuclei, labelmap_edges, labelmap_outlines); Ext.CLIJ2_release(labelmap_edges); } else labelmap_outlines = labelmap_edges; Ext.CLIJ2_pullBinary(labelmap_outlines); Ext.CLIJ2_release(labelmap_outlines); rename("nuclei_outlines"); run(outlineColor); nuclei_outlines = "nuclei_outlines"; labelmap_nuclei = "Labelmap_nuclei"; return newArray(labelmap_nuclei, nuclei_outlines, nrNuclei); } function segmentNucleiStarDist (image, channel, probabilityThreshold, pixelWidth, unit, resultTable) { selectWindow(image); getDimensions(width, height, channels, slices, frames); minNucleusSize = PI*Math.sqr((minNucleusSize_setting / pixelWidth / 2)); //Calculate the nucleus area as if it were a circle maxNucleusSize = PI*Math.sqr((maxNucleusSize_setting / pixelWidth / 2)); run("Duplicate...", "title=nuclei duplicate channels=" +channel); //Get nucleus channel starDist_input_image = getTitle(); //Downsample and Z-project.**** TO DO (maybe): Perform this on the GPU (but uses more GPU RAM) and requires splitting and combining the channels. if(downsampleFactorStarDist == 0 && unit == "µm" || unit == "um" || unit == "microns" || unit == "micron") { downsampleFactorStarDist = 0.25/pixelWidth; //scale to 0.25 um/pixel if(downsampleFactorStarDist < 1 && downsampleFactorStarDist > 0.5) downsampleFactorStarDist = 1; //Do not upsample unless the pixel size is > 0.5 um. else if(downsampleFactorStarDist < 0.25 && XYBinning>1) print ("[WARNING] The pixel size is very large ("+pixelWidth+" "+unit+"). Stardist may have difficulties segmenting the nuclei. Try setting XY binning to a lower number (currently "+XYBinning); else if(downsampleFactorStarDist < 0.25 && XYBinning==1) print ("[WARNING] The pixel size is very large ("+pixelWidth+" "+unit+"). Stardist may have difficulties segmenting the nuclei."); else if(downsampleFactorStarDist > 0.8 && downsampleFactorStarDist < 1.2) downsampleFactorStarDist = 1; //Too small change - skip rescaling print("Pixel size (after "+XYBinning+"x"+XYBinning+" XY binning): "+pixelWidth+" µm\nStarDist downsample factor (automatic): "+downsampleFactorStarDist+" (effective pixel size: "+pixelWidth*downsampleFactorStarDist+" µm)\n"); } else if (downsampleFactorStarDist == 0) { print("[WARNING] Pixel size seems incorrect ("+pixelWidth+" "+unit+")"+". Cannot determine the nuclei downsample factor for Stardist. Stardist may have difficulties segmenting the nuclei."); downsampleFactorStarDist = 1; } else { if(downsampleFactorStarDist < 0.15) { print("[WARNING] Overruling StarDist downsample factor ("+downsampleFactorStarDist+"). Will be set to the minimum value of 0.15"); downsampleFactorStarDist = 0.15; //Allow upscaling up to 6.667x } print("Pixel size (after "+XYBinning+"x"+XYBinning+" XY binning): "+pixelWidth+" µm\nStarDist downsample factor set to: "+downsampleFactorStarDist+" (effective pixel size: "+pixelWidth*downsampleFactorStarDist+" µm)\n"); } if(downsampleFactorStarDist != 1) { run("Scale...", "x="+1/downsampleFactorStarDist+" y="+1/downsampleFactorStarDist+" interpolation="+StarDistDownsampleInterpolation+" average process create"); rename("nuclei_downscaled"); starDist_input_image = getTitle(); } if(slices>1) { run("Z Project...", "projection=[Max Intensity]"); rename("nuclei_Zprojected"); starDist_input_image = getTitle(); if(downsampleFactorStarDist != 1) close("nuclei_downscaled"); } //Run StarDist getDimensions(dswidth, dsheight, dschannels, dsslices, dsframes); starDistTiles = pow(floor((maxOf(dswidth, dsheight)/maxTileSize)-1)+1,2); //Determine the nr. of tiles run("Command From Macro", "command=[de.csbdresden.stardist.StarDist2D], args=['input':'"+starDist_input_image+"', 'modelChoice':'Versatile (fluorescent nuclei)', 'normalizeInput':'true', 'percentileBottom':'1.0', 'percentileTop':'99.60000000000001', 'probThresh':'"+probabilityThreshold+"', 'nmsThresh':'0.3', 'outputType':'ROI Manager', 'nTiles':'"+starDistTiles+"', 'excludeBoundary':'2', 'roiPosition':'Stack', 'verbose':'false', 'showCsbdeepProgress':'false', 'showProbAndDist':'false'], process=[false]"); //Scale up ROIs if(downsampleFactorStarDist != 1) RoiManager.scale(downsampleFactorStarDist, downsampleFactorStarDist, false); //Spline fit ROIs selectWindow(image); //Have to select an image with the correct size for (i = 0; i < roiManager("count"); i++) { roiManager("select", i); run("Fit Spline"); // Roi.getSplineAnchors(x, y); // Roi.setPolygonSplineAnchors(x, y); roiManager("Update"); } //Convert ROIs to label map run("Select None"); labelmap_nuclei = "Labelmap_nuclei_unfiltered"; // run("ROI Manager to LabelMap(2D)"); labelmap_nuclei = ROI_Manager_to_labelmap(image); setBatchMode("show"); run("glasbey_on_dark"); setMinAndMax(0, 255); run("Clear Results"); selectWindow(labelmap_nuclei); getDimensions(uswidth, usheight, uschannels, usslices, usframes); if(uswidth != width || usheight != height) run("Canvas Size...", "width="+width+" height="+height+" position=Center zero"); // Make sure that the size of the upscaled labelmap is correct - upscaling can cause rounding errors //Close unused images close("Label Image"); if(slices>1) close("nuclei_Zprojected"); //Exclude nuclei on edges and count nr of nuclei Ext.CLIJ2_push(labelmap_nuclei); close(labelmap_nuclei); //Ext.CLIJ2_excludeLabelsWithValuesOutOfRange(area, labelmap_nuclei, labelmap_nuclei_filtered, minNucleusSize, maxNucleusSize); if(excludeOnEdges) { Ext.CLIJ2_excludeLabelsOnEdges(labelmap_nuclei, labelmap_nuclei_edges_excluded); Ext.CLIJ2_release(labelmap_nuclei); } else labelmap_nuclei_edges_excluded = labelmap_nuclei; Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei_edges_excluded, nrNucleiBeforeSizeFiltering); Ext.CLIJ2_excludeLabelsOutsideSizeRange(labelmap_nuclei_edges_excluded, labelmap_nuclei_filtered, minNucleusSize, maxNucleusSize); Ext.CLIJ2_release(labelmap_nuclei_edges_excluded); labelmap_nuclei_final = "Labelmap_nuclei"; Ext.CLIJ2_closeIndexGapsInLabelMap(labelmap_nuclei_filtered, labelmap_nuclei_gapsClosed); // if(downsampleFactorStarDist > 1) { Ext.CLIJ2_greyscaleOpeningSphere(labelmap_nuclei_gapsClosed, labelmap_nuclei_final, round(downsampleFactorStarDist + 1), round(downsampleFactorStarDist + 1), 0); //Smooth labels a bit Ext.CLIJ2_release(labelmap_nuclei_gapsClosed); // } // else labelmap_nuclei_final = labelmap_nuclei_gapsClosed; Ext.CLIJ2_release(labelmap_nuclei_filtered); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei_final, nrNuclei); //Count the number of nuclei print(nrNucleiBeforeSizeFiltering - nrNuclei + " nuclei were removed due to size restrictions (" + minNucleusSize_setting + " - " + maxNucleusSize_setting + " " + unit + " | "+d2s(minNucleusSize,0)+" - "+d2s(maxNucleusSize,0)+" pixels)."); //Create nuclei outlines from nuclei labelmap Ext.CLIJ2_detectLabelEdges(labelmap_nuclei_final, labelmap_edges); if(thickOutlines == false) { Ext.CLIJ2_mask(labelmap_nuclei_final, labelmap_edges, labelmap_outlines); Ext.CLIJ2_release(labelmap_edges); } else labelmap_outlines = labelmap_edges; Ext.CLIJ2_pullBinary(labelmap_outlines); Ext.CLIJ2_release(labelmap_outlines); rename("nuclei_outlines"); run(outlineColor); nuclei_outlines = "nuclei_outlines"; labelmap_nuclei_final = "Labelmap_nuclei"; return newArray(labelmap_nuclei_final, nuclei_outlines, nrNuclei); } function loadROIs(image, channel, folder) { // File.setDefaultDir(outputFolder); run("Duplicate...", "title=nuclei duplicate channels=" +channel); //Get nucleus channel roiManager("reset"); if(nucleiSegmentationChoice == "Load ROIs from file") { if(File.exists(folder + File.separator + File.getNameWithoutExtension(image) + "_ROIs.zip")) roiManager("open", folder + File.separator + File.getNameWithoutExtension(image) + "_ROIs.zip"); else { print("WARNING: ROI file not found! "+folder + File.separator + File.getNameWithoutExtension(image) + "_ROIs.zip\nSkipping this image"); return newArray("FileNotFound"); } run("ROI Manager to LabelMap(2D)"); } else if(nucleiSegmentationChoice == "Load label images") { if(File.exists(folder + File.separator + File.getNameWithoutExtension(image) + ".tif")) open(folder + File.separator + File.getNameWithoutExtension(image) + ".tif"); //Label images must have the same name as the original image else { print("WARNING: Labelmap file not found! "+folder + File.separator + File.getNameWithoutExtension(image) + ".tif\nSkipping this image"); return newArray("FileNotFound"); } } labelmap_nuclei_final = "Labelmap_nuclei"; rename(labelmap_nuclei_final); Ext.CLIJ2_push(labelmap_nuclei_final); // close(labelmap_nuclei); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei_final, nrNuclei); //Count the number of nuclei //Create nuclei outlines from nuclei labelmap Ext.CLIJ2_detectLabelEdges(labelmap_nuclei_final, labelmap_edges); if(thickOutlines == false) { Ext.CLIJ2_mask(labelmap_nuclei_final, labelmap_edges, labelmap_outlines); Ext.CLIJ2_release(labelmap_edges); } else labelmap_outlines = labelmap_edges; Ext.CLIJ2_pullBinary(labelmap_outlines); Ext.CLIJ2_release(labelmap_outlines); rename("nuclei_outlines"); run(outlineColor); nuclei_outlines = "nuclei_outlines"; labelmap_nuclei_final = "Labelmap_nuclei"; return newArray(labelmap_nuclei_final, nuclei_outlines, nrNuclei); } function segmentCellsCellpose (image, nucleiChannel, cytoChannel, probabilityThreshold, pixelWidth, pixelDepth, unit, resultTable) { //Check which version of the BIOP Cellpose wrapper is present, if at all. List.setCommands; if (List.get("Cellpose ...")!="") { //Get Cellpose settings //This works for the new wrapper envPath = getPref("Packages.ch.epfl.biop.wrappers.cellpose.ij2commands.Cellpose", "env_path"); envType = getPref("Packages.ch.epfl.biop.wrappers.cellpose.ij2commands.Cellpose", "env_type"); if(envType == "