/* 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://github.com/BioImaging-NKI/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, b.vd.broek@nki.nl * * * * 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 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) * * --> TO DO: open timelapse output files, concatenate, and restore overlays. Now via external scripts. * */ #@ String file_message (value="

File settings

", visibility="MESSAGE") #@ File[] files (label = "Input files", style="File") #@ File outputFolder (label = "Output folder", style = "directory") #@ String image_message (value="

Image settings

", visibility="MESSAGE") #@ Boolean loadSettingsFromFile (label = "Load settings from analyzed image? (will not update values in this dialog)") #@ Integer nucleiChannel (label = "Nuclei channel (-1 if not used)", value = 1) #@ Integer cytoChannel (label = "Cytoplasm/membrane channel (-1 if not used)", value = -1) #@ Integer fociChannelA (label = "Foci channel A", value = 2) #@ Integer fociChannelB (label = "Foci channel B", value = 3) #@ Boolean detectFociChannelB (label = "Also detect foci in channel B and perform colocalization?", persistence=true, value=true, description="Check this box to calculate colocalization between two foci channels.") #@ String ThreeDHandling (label = "3D image handling", choices={"Detect foci in 3D (or 2D when not applicable)", "Detect foci on the Maximum Intensity Projection", "Detect foci on the Standard Deviation Projection", "Detect foci on the Summed Intensity 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") #@ Integer singleSlice (label = "[single z-slice foci detection only] Slice nr", value=1, min=1) #@ Integer cropBorder (label = "Remove image borders (pixels)", value = 0, min=0) #@ Integer XYBinning (label = "Image XY binning before analysis [1-n]", value = 1, min=1, description="Image binning will greatly shorten the analysis time, but can be less accurate.") #@ 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 image"}, style="listBox") #@ Double downsampleFactorStarDist(label = "Stardist nuclei rescaling factor [1-n], 0 for automatic, 1 for no rescaling", value = 0, min=0, description="A binning >1 for high-resolution input images can prevent 'oversegmentation', and will also somewhat speed up processing.") #@ Double probabilityThreshold (label = "Probablility/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") #@ String CellposeModel (label = "Cellpose model", choices={"cyto3","cyto2_cp3","nuclei","custom"}, style="listBox", value="cyto3") #@ Integer CellposeDiameter (label = "Cellpose cell diameter (pixels), 0 for automatic", value = 0, min=0, description="Estimation of the cell diameter") #@ Integer minNucleusSize_setting (label = "Remove nulei with diameter smaller than (µm)", value = 4) #@ Integer maxNucleusSize_setting (label = "Remove nulei with diameter larger than (µm)", value = 50) #@ Boolean excludeOnEdges (label = "Exclude nulei on image edges", value = true) #@ String manualRemoveNuclei (label = "Manually remove segmentated nuclei/cells", choices={"No thanks","Manually remove nuclei", "load previously saved removals (from output folder)"}, value = "No thanks", description="Allow nuclei/cell removal by clicking, or load previous removals (Masks should be in the output folder)") #@ String foci_message1 (value="

Foci detection settings

", visibility="MESSAGE") #@ Boolean optimizationModeSetting (label = "Show preview for foci parameters fine-tuning?", value=true) #@ String fociSizeA (label = "Foci size channel A (after XY binning)", choices={"tiny","small","average","large","huge","other"}, style="listBox", value="average", description="Simplified foci detection parameter. This parameter determines the Difference of Gaussians filter size for foci detection.") #@ String fociSizeB (label = "Foci size channel B (after XY binning)", choices={"tiny","small","average","large","huge","other"}, style="listBox", value="average", description="Simplified foci detection parameter. This parameter determines the Difference of Gaussians filter size for foci detection.") #@ String fociDetectionMethod (label = "Foci detection method", choices={"Marker-controlled watershed (recommended)","AreaMaxima local maximum detection"}, style="listBox") #@ Double thresholdFactorA (label = "Foci intensity threshold bias channel A", value = 0.0, min=-2.5, max=2.5, style="scroll bar", stepSize=0.01) #@ Double thresholdFactorB (label = "Foci intensity threshold bias channel B", value = 0.0, min=-2.5, max=2.5, style="scroll bar", stepSize=0.01) #@ Integer minFociSize (label = "Minimum foci size (area) (pixels/voxels)", value = 3) #@ Integer maxFociSize (label = "Maximum foci size (area) (pixels/voxels)", value = 9999) #@ Double maxFociDistanceOutsideNuclei_setting (label = "Max distance of foci outside nuclei/cells (µm); -1 for full image", value = 0, min=-1, style="format:0.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"}) #@ ColorRGB fontColor (label = "Nuclei/cell label color", value="orange") #@ Integer labelFontSize (label = "Nuclei/cell label font size", value=12, min=4) #@ String overlayBrightness (label = "Nuclei/cell outline bightness", choices={"bright","dim"}) #@ String outlineColor (label = "Nuclei/cell outline color", choices={"White","Red","Green","Blue","Cyan","Magenta","Yellow"}, value = "Cyan") #@ Boolean debugMode (label = "Debug mode (show intermediate images)", value=false) #@ String file_and_image_message0 (value = "
Need help? Visit Foci Analyzer on GitHub", visibility="MESSAGE") version = 1.7; requires("1.54a"); //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 * * 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) * * Open and split timelapse image; save data as individual .tsv files? (in order to prevent too much bookkeeping) * * 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 != "") exit("Error: "+missingPlugin+" not found!\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."); if(loadSettingsFromFile == true) { 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 "); fontColor = Property.get("33. Nuclei label color "); labelFontSize = Property.getNumber("34. Nuclei label size "); overlayBrightness = Property.get("35. Nuclei/cell outline bightness "); outlineColor = Property.get("36. 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 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.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 image") { labelImageFolder = call("ij.Prefs.get", "label.Image.Folder", File.getParent(files[0])); Dialog.createNonBlocking("Select a label image folder"); 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"); } //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 image") { 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 Standard Deviation Projection") { selectWindow(original); run("Z Project...", "projection=[Max Intensity]"); setBatchMode("show"); original = getTitle(); gslices = 1; } if(gslices > 1 && ThreeDHandling == "Detect foci on the Maximum Intensity 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 image") 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 not applicable)" && 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=4 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("Max distance of foci outside nuclei/cells", 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 (relative thresholds)", "Done optimizing | Process all other images with these settings (fixed absolute thresholds)"); Dialog.addChoice("Action", items, "Recalculate with these settings"); Dialog.addMessage("Click Help for more info (Foci Analyzer GitHub site)", 12, "#000080"); Dialog.addHelp("https://github.com/BioImaging-NKI/Foci-analyzer/blob/main/README.md"); //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 (relative thresholds)") { doneOptimizing = true; processAllOtherImages = true; } if(processChoice == "Done optimizing | Process all other images with these settings (fixed absolute thresholds)") { 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_width = Table.getColumn("BOUNDING_BOX_WIDTH"); boundingBox_height = Table.getColumn("BOUNDING_BOX_HEIGHT"); Array.getStatistics(boundingBox_width, minWidth, maxWidth); Array.getStatistics(boundingBox_height, minHeight, maxHeight); 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, maxWidth, maxHeight); if(detectFociChannelB) nrFoci = measureFoci(original, fociChannelB, nrNuclei, labelmap_nuclei_3D, labelmap_fociB, boundingBox_X, boundingBox_Y, maxWidth, maxHeight); 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 ", fontColor); List.set("34. Nuclei/cell label size ", labelFontSize); List.set("35. Nuclei/cell outline bightness ", overlayBrightness); List.set("36. 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_overlap_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 if(nrNucleiBeforeSizeFiltering - nrNuclei > 0) print(nrNucleiBeforeSizeFiltering - nrNuclei + " nuclei were removed due to size restrictions (" + minNucleusSize_setting + " - " + maxNucleusSize_setting + " " + unit + ")"); //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 image") { 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 == "") envType = "conda"; //Default is conda, but returns oldCellposeWrapper = false; } else if (List.get("Cellpose Advanced")!="") { oldCellposeWrapper = true; // This works for the old wrapper, but not for the new one... // prefs = split(File.openAsString(getDir("preferences") + "IJ_Prefs.txt"),"\n"); // for (i = 0; i < prefs.length; i++) { // data = split(prefs[i],"="); // if(data.length>1) List.set(data[0], data[1]); // } // envPath = List.get(".ch.epfl.biop.wrappers.cellpose.Cellpose.envDirPath"); // envPath = envPath.replace("\\\\:", ":"); //This line is needed because sometimes the path has an extra '\' after the drive letter. // envType = List.get(".ch.epfl.biop.wrappers.cellpose.Cellpose.envType"); print("[WARNING] Old version of the Cellpose wrapper detected. It may work now, but may be broken in the future. Please Update Fiji via Help -> Update..." } else exit("ERROR: Cellpose wrapper not found! Update Fiji and activate the PT-BIOP Update site."); List.clear(); if(nucleiSegmentationChoice == "Cellpose segmentation 3D" && ThreeDHandling != "Detect foci in 3D (or 2D when not applicable)" && gslices > 1) { print("[WARNING] 3D Cellpose segmentation is selected, but '3D image handling' is set to '"+ThreeDHandling+"'. Proceeding with 2D Cellpose segmentation."); nucleiSegmentationChoice = "Cellpose segmentation 2D (or on 2D projection)"; } if(nucleiSegmentationChoice == "Cellpose segmentation 3D" && ThreeDHandling == "Detect foci in 3D (or 2D when not applicable)" && gslices == 1) { print("[WARNING] 3D Cellpose segmentation is selected, but the image has only 1 slice. Proceeding with 2D Cellpose segmentation."); nucleiSegmentationChoice = "Cellpose segmentation 2D (or on 2D projection)"; } 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)); cellpose_input_image = getTitle(); if(cytoChannel<=0) { run("Duplicate...", "title=nuclei duplicate channels=" +nucleiChannel); //Get nuclei for segmentation cellpose_input_image = getTitle(); } else run("Duplicate...", "title=nuclei duplicate channels=" +cytoChannel); //Get cells (and pretend they are nuclei) Cellpose_minCellIntensity = 0; //To do: maybe make this a parameter for 2D Cellpose as well if(slices>1) { if(nucleiSegmentationChoice == "Cellpose segmentation 2D (or on 2D projection)") { run("Z Project...", "projection=[Max Intensity]"); rename("forCellpose_Zprojected"); cellpose_input_image = getTitle(); } else if(nucleiSegmentationChoice == "Cellpose segmentation 3D") { if(oldCellposeWrapper == true) exit("Sorry, 3D segmentation with Foci Analyzer is not compatible with the the old BIOP Cellpose wrapper.\nRun Help -> Update... and update the PTBIOP plugins to get the new one."); //Get Prefs Cellpose_use_GPU = call("ij.Prefs.get", "Cellpose.use.GPU", true); Cellpose_3Dseg_type = call("ij.Prefs.get", "Cellpose.3Dseg.type", "2.5D"); if(CellposeDiameter == 0) CellposeDiameter = call("ij.Prefs.get", "Cellpose.diameter", 30); Cellpose_cellprob_threshold = call("ij.Prefs.get", "Cellpose.cellprob.threshold", 0.0); Cellpose_stitch_threshold = call("ij.Prefs.get", "Cellpose.stitch.threshold", 0.0); Cellpose_dP_smooth = call("ij.Prefs.get", "Cellpose.dP.smooth", 0.0); Cellpose_min_size = call("ij.Prefs.get", "Cellpose.min.size", 15); Cellpose_additional_parameters = call("ij.Prefs.get", "Cellpose.additional.parameters", ""); Cellpose_normalize3D = call("ij.Prefs.get", "Cellpose.normalize.3D", false); Cellpose_minCellIntensity = call("ij.Prefs.get", "Cellpose.min.cell.intensity", 0); Dialog.createNonBlocking("Cellpose 3D advanced settings"); Dialog.addMessage("Cellpose 3D segmentation parameters", 14, "#000080"); Dialog.addMessage("If unsure, leave as is or click the Help button below for more info.", 12, "#000080"); Dialog.addChoice("Segmentation type", newArray("2.5D", "3D"), Cellpose_3Dseg_type); // Dialog.setInsets(-25, 250, 0); Dialog.setInsets(-25, 350, 0); Dialog.addMessage("'2.5D' means: 2D segmentation on every slice + stitching together"); Dialog.addNumber("--diameter", CellposeDiameter, 1, 6, "pixels (N.B. '0 = automatic' does not work for '3D')"); Dialog.addNumber("--flow_threshold (2.5D only)", probabilityThreshold, 1, 6, "[0-3], default: 0.4"); Dialog.addNumber("--cellprob_threshold", Cellpose_cellprob_threshold, 1, 6, "[(-6)-6], default: 0.0"); Dialog.addNumber("--anisotropy", d2s(pixelDepth/pixelWidth,2), 2, 6, "obtained from image dimensions"); Dialog.addNumber("--stitch_threshold", Cellpose_stitch_threshold, 2, 6, "[0-1], 2D + stitching across planes; only for '2.5D'"); Dialog.addNumber("--dP_smooth (3D only)", Cellpose_dP_smooth, 1, 6, "stddev of Gaussian filter for smoothing flows"); Dialog.addNumber("--min_size", Cellpose_min_size, 0, 6, "pixels, default: 15"); Dialog.addString("Additional Cellpose parameters", Cellpose_additional_parameters, 60); // Dialog.setInsets(0, 258, 0); //probably depends on display scaling and/or font size... Dialog.setInsets(0, 290, 0); //probably depends on display scaling and/or font size... Dialog.addCheckbox("--use_gpu", Cellpose_use_GPU); Dialog.setInsets(0, 290, 0); Dialog.addCheckbox("Normalize the image in 3D? (Otherwise slice by slice, possibly leading to falsely detected cells in dim slices.)", Cellpose_normalize3D); Dialog.setInsets(10, 0, 0); Dialog.addNumber("Remove cells with mean intensity lower than", Cellpose_minCellIntensity, 1, 6, "gray values"); Dialog.addMessage("Additional Cellpose parameters are e.g. '--no_norm, --restore_type, nuclei, --niter, 200'. Click the Help button for a list of all parameters."); Dialog.addHelp("https://cellpose.readthedocs.io/en/latest/cli.html"); Dialog.show(); Cellpose_use_GPU = Dialog.getCheckbox(); Cellpose_3Dseg_type = Dialog.getChoice(); CellposeDiameter = Dialog.getNumber(); probabilityThreshold = Dialog.getNumber(); Cellpose_cellprob_threshold = Dialog.getNumber(); Cellpose_anisotropy = Dialog.getNumber(); Cellpose_stitch_threshold = Dialog.getNumber(); Cellpose_dP_smooth = Dialog.getNumber(); Cellpose_min_size = Dialog.getNumber(); Cellpose_additional_parameters = Dialog.getString(); Cellpose_normalize3D = Dialog.getCheckbox(); Cellpose_minCellIntensity = Dialog.getNumber(); //Set Prefs call("ij.Prefs.set", "Cellpose.use.GPU", Cellpose_use_GPU); call("ij.Prefs.set", "Cellpose.3Dseg.type", Cellpose_3Dseg_type); call("ij.Prefs.set", "Cellpose.diameter", CellposeDiameter); setScriptParameterValue(CellposeDiameter, CellposeDiameter); call("ij.Prefs.set", "Cellpose.cellprob.threshold", Cellpose_cellprob_threshold); setScriptParameterValue(Cellpose_cellprob_threshold, Cellpose_cellprob_threshold); call("ij.Prefs.set", "Cellpose.stitch.threshold", Cellpose_stitch_threshold); call("ij.Prefs.set", "Cellpose.dP.smooth", Cellpose_dP_smooth); call("ij.Prefs.set", "Cellpose.min.size", Cellpose_min_size); call("ij.Prefs.set", "Cellpose.additional.parameters", Cellpose_additional_parameters); call("ij.Prefs.set", "Cellpose.normalize.3D", Cellpose_normalize3D); call("ij.Prefs.set", "Cellpose.min.cell.intensity", Cellpose_minCellIntensity); if(Cellpose_use_GPU == true) Cellpose_use_GPU = "--use_gpu"; if(Cellpose_3Dseg_type == "3D") Cellpose_3Dseg_type = "--do_3D"; else Cellpose_3Dseg_type = ""; if(Cellpose_normalize3D == true) { Ext.CLIJ2_push(cellpose_input_image); Ext.CLIJx_normalize(cellpose_input_image, cellpose_input_image_normalized); Ext.CLIJ2_pull(cellpose_input_image_normalized); Ext.CLIJ2_release(cellpose_input_image); Ext.CLIJ2_release(cellpose_input_image_normalized); cellpose_input_image = cellpose_input_image_normalized; Cellpose_additional_parameters += ", --no_norm"; } } } setBatchMode("exit and display"); //Cellpose doesn't work in batch mode selectWindow(cellpose_input_image); showStatus("Performing Cellpose segmentation... (check the Console for live info)"); if(nucleiSegmentationChoice == "Cellpose segmentation 2D (or on 2D projection)") { if(cytoChannel > 0) run("Cellpose ...", "env_path="+envPath+" env_type="+envType+" model=["+CellposeModel+"] model_path=["+CellposeModelPath+"] diameter="+CellposeDiameter+" ch1="+cytoChannel+" ch2="+nucleiChannel+" additional_flags=[--use_gpu, --flow_threshold, "+probabilityThreshold+", --cellprob_threshold, 0.0]"); else run("Cellpose ...", "env_path="+envPath+" env_type="+envType+" model=["+CellposeModel+"] model_path=["+CellposeModelPath+"] diameter="+CellposeDiameter+" ch1="+cytoChannel+" ch2=0 additional_flags=[--use_gpu, --flow_threshold, "+probabilityThreshold+", --cellprob_threshold, 0.0]"); } else if(nucleiSegmentationChoice == "Cellpose segmentation 3D") { if(cytoChannel > 0) run("Cellpose ...", "env_path="+envPath+" env_type="+envType+" model=["+CellposeModel+"] model_path=["+CellposeModelPath+"] diameter="+CellposeDiameter+" ch1="+cytoChannel+" ch2="+nucleiChannel+" additional_flags=["+Cellpose_use_GPU+", "+Cellpose_3Dseg_type+", --flow_threshold, "+probabilityThreshold+", --cellprob_threshold, "+Cellpose_cellprob_threshold+", --stitch_threshold, "+Cellpose_stitch_threshold+", --anisotropy, "+Cellpose_anisotropy+", --dP_smooth, "+Cellpose_dP_smooth+", --min_size, "+Cellpose_min_size+", "+Cellpose_additional_parameters+"]"); else run("Cellpose ...", "env_path="+envPath+" env_type="+envType+" model=["+CellposeModel+"] model_path=["+CellposeModelPath+"] diameter="+CellposeDiameter+" ch1="+nucleiChannel+" ch2=0 additional_flags=["+Cellpose_use_GPU+", "+Cellpose_3Dseg_type+", --flow_threshold, "+probabilityThreshold+", --cellprob_threshold, "+Cellpose_cellprob_threshold+", --stitch_threshold, "+Cellpose_stitch_threshold+", --anisotropy, "+Cellpose_anisotropy+", --dP_smooth, "+Cellpose_dP_smooth+", --min_size, "+Cellpose_min_size+", "+Cellpose_additional_parameters+"]"); } else if(oldCellposeWrapper == true) run("Cellpose Advanced", "diameter="+CellposeDiameter+" cellproba_threshold=0 flow_threshold="+probabilityThreshold+" anisotropy=1.0 diam_threshold=12.0 model="+CellposeModel+" nuclei_channel=0 cyto_channel=1 dimensionmode=2D stitch_threshold=-1.0 omni=false cluster=false additional_flags="); //open("C:/Users/b.vd.broek/surfdrive/DATA/Luka Wolf/Labelmap_nuclei.tif"); //TEMP FAST PROCESSING labelmap_nuclei = getTitle(); // selectWindow(image); selectWindow(labelmap_nuclei); setBatchMode("hide"); // selectWindow("nuclei"); // setBatchMode("hide"); //Close unused images if(slices>1 && nucleiSegmentationChoice == "Cellpose segmentation 2D (or on 2D projection)") close("forCellpose_Zprojected"); //Exclude nuclei on edges and count nr of nuclei Ext.CLIJ2_push(labelmap_nuclei); close(labelmap_nuclei); 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); //Count the number of nuclei 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_final); Ext.CLIJ2_release(labelmap_nuclei_filtered); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei_final, nrNuclei); //Count the number of nuclei Ext.CLIJ2_pull(labelmap_nuclei_final); if(nrNucleiBeforeSizeFiltering - nrNuclei > 0) print(nrNucleiBeforeSizeFiltering - nrNuclei + " cells were removed due to size restrictions (" + minNucleusSize_setting + " - " + maxNucleusSize_setting + " " + unit + ")"); if(nrNuclei > 0 && Cellpose_minCellIntensity > 0) { if(debugMode) { showImagefromGPU(labelmap_nuclei_final); run("glasbey on dark"); } run("Intensity Measurements 2D/3D", "input=nuclei labels=Labelmap_nuclei mean volume"); //OF TOCH BETER CLIJ2_GETSTATISTICSOFBACKGROUNDANDLABELEDPIXELS?? labels_meanIntensity = Table.getColumn("Mean", "nuclei-intensity-measurements"); close("nuclei-intensity-measurements"); //labels_meanIntensity = Table.getColumn("MEAN_INTENSITY", "Results"); labels_meanIntensity = prependToArray(0, labels_meanIntensity); Ext.CLIJ2_pushArray(labels_meanIntensity_vector, labels_meanIntensity, labels_meanIntensity.length, 1, 1); Ext.CLIJ2_generateParametricImage(labelmap_nuclei_final, labels_meanIntensity_vector, meanIntensity_map); Ext.CLIJ2_withinIntensityRange(meanIntensity_map, meanIntensity_map_filtered, Cellpose_minCellIntensity, 1E10); Ext.CLIJ2_release(meanIntensity_map); //labelmap_nuclei_filtered = "labelmap_nuclei_filtered"; //labelmap_outlines_filtered = "labelmap_nuclei_filtered_gapsclosed"; Ext.CLIJ2_multiplyImages(labelmap_nuclei_final, meanIntensity_map_filtered, labelmap_nuclei_filtered); Ext.CLIJ2_closeIndexGapsInLabelMap(labelmap_nuclei_filtered, labelmap_nuclei_filtered_gapsclosed); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_nuclei_filtered_gapsclosed, nrNuclei); Ext.CLIJ2_release(labelmap_nuclei_filtered); //Ext.CLIJ2_multiplyImages(nuclei_outlines_org, meanIntensity_map_filtered, nuclei_outlines_filtered); //Ext.CLIJ2_pull(nuclei_outlines_filtered); //nuclei_outlines = nuclei_outlines_filtered; Ext.CLIJ2_release(meanIntensity_map_filtered); labelmap_nuclei_final = labelmap_nuclei_filtered_gapsclosed; if(debugMode) { showImagefromGPU(labelmap_nuclei_filtered_gapsclosed); run("glasbey on dark"); } } //Create nuclei outlines from nuclei labelmap if (nucleiSegmentationChoice == "Cellpose segmentation 2D (or on 2D projection)") { 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); } else if(nucleiSegmentationChoice == "Cellpose segmentation 3D") { //Create a 3D image with label edges, without the ugly 'filled planes' that Ext.CLIJ2_detectLabelEdges() outputs. Ext.CLIJ2_maximumSliceBySliceSphere(labelmap_nuclei_final, labelmap_maximum, 1, 1); Ext.CLIJ2_minimumSliceBySliceSphere(labelmap_nuclei_final, 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, labelmap_outlines, 1); Ext.CLIJ2_pull(labelmap_outlines); Ext.CLIJ2_release(label_edges); Ext.CLIJ2_release(labelmap_outlines); rename("nuclei_outlines"); run(outlineColor); } nuclei_outlines = "nuclei_outlines"; return newArray(labelmap_nuclei_final, nuclei_outlines, nrNuclei); } //Manually remove labels by creating point selections. Masks are saved as PNG files (compressed) for later use. function manually_remove_labels(labelmap, label_edges, nrNuclei, windowName, imageName) { selectWindow(windowName); Stack.setDisplayMode("composite"); Stack.setChannel(nucleiChannel); setLut(b_reds, b_greens, b_blues); Stack.setChannel(fociChannelA); run("Green"); if(detectFociChannelB) { Stack.setChannel(fociChannelB); run("Red"); } run("Add Image...", "image="+label_edges+" x=0 y=0 opacity=100 zero"); setBatchMode("show"); if(manualRemoveNuclei == "load previously saved removals (from output folder)" && File.exists(outputFolder + File.separator + imageName + "__removalMask.png")) { open(outputFolder + File.separator + imageName + "__removalMask.png"); if(is("Inverting LUT")) run("Invert LUT"); run("Points from Mask"); } else if(manualRemoveNuclei == "Manually remove nuclei") { setTool("multipoint"); run("Select None"); run("Point Tool...", "type=Dot color=White size=[Extra Large] show counter=0"); waitForUser("Select nuclei to remove, then press OK"); if(selectionType == 10) { //Only create mask if there is a selection run("Create Mask"); saveAs("png", outputFolder + File.separator + imageName+"__removalMask"); } else if(File.exists(outputFolder + File.separator + imageName + "__removalMask.png")) File.delete(outputFolder + File.separator + imageName + "__removalMask.png"); //Remove Mask } selectWindow(windowName); run("Select None"); close(imageName + "__removalMask"); Ext.CLIJ2_getMaximumOfAllPixels(labelmap, nrLabels); Ext.CLIJ2_pull(labelmap); selectWindow(labelmap); run("Restore Selection"); run("Clear Results"); run("Set Measurements...", "mean redirect=None decimal=3"); run("Measure"); run("Select None"); setTool("rectangle"); nuclei_to_be_removed = Table.getColumn("Mean", "Results"); //List of nuclei to be removed flaglist = "flaglist"; newImage(flaglist, "32-bit black", nrLabels+1, 1, 1); //+1 because background = 0 for (i = 0; i < nuclei_to_be_removed.length; i++) setPixel(nuclei_to_be_removed[i], 0, 1); Ext.CLIJ2_push(flaglist); setBatchMode("show"); close(flaglist); labelmap_edited = "Labelmap_nuclei_edited"; Ext.CLIJ2_excludeLabels(flaglist, labelmap, labelmap_edited); Ext.CLIJ2_release(labelmap); labelmap = "Labelmap_nuclei_edited"; Ext.CLIJ2_getMaximumOfAllPixels(labelmap_edited, nrNuclei); label_edges_edited = "label_edges_edited"; Ext.CLIJ2_detectLabelEdges(labelmap_edited, label_edges_edited); Ext.CLIJ2_pull(label_edges_edited); run("Cyan"); resetMinAndMax; //Subtract new label_edges from original label_edges imageCalculator("Subtract", label_edges, label_edges_edited); run("Red"); setMinAndMax(0, 2); //Add nuclei outlines as overlay to the original image selectWindow(original); Overlay.remove; run("Add Image...", "image="+label_edges+" x=0 y=0 opacity=100 zero"); run("Add Image...", "image="+label_edges_edited+" x=0 y=0 opacity=100 zero"); close("label_edges"); // draw_label_numbers(labelmap_edited, labelFontSize, fontColor); updateDisplay(); close("label_edges"); return newArray(labelmap_edited, label_edges_edited, nrNuclei); } function ROI_Manager_to_labelmap(image) { setBatchMode(true); getDimensions(width, height, channels, slices, frames); newImage("labelmap", "16-bit", width, height, 1); for (i = 0; i < roiManager("count"); i++) { roiManager("select", i); changeValues(0, 65535, i+1); } run("Select None"); return "labelmap"; } function detect_foci(image, channel, fociSize, anisotropyFactor, firstTimeProcessing, labelmap_nuclei, labelmap_nuclei_3D, thresholdFactor, threshold) { //Determine threshold on the maximum projection (more robust), only where nuclei are present // Small FLAW: the area outside the nuclei is set to zero, but is included in the threshold calculation. (Could be fixed by setting 0 to NaN, but how to do that on the GPU?) // On the other hand: there is almost no difference, it seems. selectWindow(image); foci = "foci_ch"+channel; if(firstTimeProcessing == true) { rename(foci); Stack.setChannel(channel); Ext.CLIJ2_pushCurrentZStack(foci); rename(image); //Ext.CLIJ2_release(foci); if(debugMode) showImagefromGPU(foci); } else if(debugMode) showImagefromGPU(foci); //Determine filter sizes if(fociSize == "tiny") { fociFilterSizeXY = 0.25; fociSizeXY = 0.5; } if(fociSize == "small") { fociFilterSizeXY = 0.5; fociSizeXY = 1; } if(fociSize == "average"){fociFilterSizeXY = 1.0; fociSizeXY = 1.5; } if(fociSize == "large") { fociFilterSizeXY = 2; fociSizeXY = 2; } if(fociSize == "huge") { fociFilterSizeXY = 4; fociSizeXY = 3; } if(fociSize == "other") { fociSize = getNumber("Enter foci size in pixels", 2); fociFilterSizeXY = fociSize; fociSizeXY = fociSize; } //pull from GPU, remove outliers, push to GPU Ext.CLIJ2_convertFloat(foci, foci_32bit); Ext.CLIJ2_pull(foci_32bit); rename(foci+"_outliers_removed"); foci_outliers_removed = foci+"_outliers_removed"; selectWindow(foci+"_outliers_removed"); run("Remove Outliers", "block_radius_x="+minOf(25*fociSizeXY, 50)+" block_radius_y="+minOf(25*fociSizeXY, 50)+" standard_deviations=2 stack"); if(debugMode) setBatchMode("show"); Ext.CLIJ2_push(foci_outliers_removed); // if(gslices>1) Ext.CLIJ2_mask(foci_outliers_removed, labelmap_nuclei_3D, foci_outliers_removed_masked); // else Ext.CLIJ2_mask(foci_outliers_removed, labelmap_nuclei, foci_outliers_removed_masked); // showImagefromGPU(foci_outliers_removed_masked); //Determine the median of standard deviation signals in the oulier-removed nuclei run("Clear Results"); Ext.CLIJ2_statisticsOfLabelledPixels(foci_outliers_removed, labelmap_nuclei); nucleiStddev_ = Table.getColumn("STANDARD_DEVIATION_INTENSITY", "Results"); medianNucleiStddev = medianOfArray(nucleiStddev_); print("Detecting foci in channel "+channel); print("Median of the Standard Deviations of all nuclei: "+medianNucleiStddev); fociFilterSizeZ = 0.5 * anisotropyFactor * fociFilterSizeXY; //semi-arbitrary setting fociSizeZ = 0.5 * anisotropyFactor * fociSizeXY; //Filter the foci - Difference of Gaussians, but then subtracting a blurred outlier-removed image if(gslices>1) Ext.CLIJ2_gaussianBlur3D(foci_32bit, foci_blurred, fociFilterSizeXY/2, fociFilterSizeXY/2, fociFilterSizeZ/2); else Ext.CLIJ2_gaussianBlur2D(foci_32bit, foci_blurred, fociFilterSizeXY/2, fociFilterSizeXY/2); if(gslices>1) Ext.CLIJ2_gaussianBlur3D(foci_outliers_removed, foci_outliers_removed_blurred, fociFilterSizeXY*4, fociFilterSizeXY*4, fociFilterSizeZ*4); else Ext.CLIJ2_gaussianBlur2D(foci_outliers_removed, foci_outliers_removed_blurred, fociFilterSizeXY*4, fociFilterSizeXY*4); if(debugMode) showImagefromGPU(foci_blurred); if(debugMode) showImagefromGPU(foci_outliers_removed_blurred); Ext.CLIJ2_release(foci_outliers_removed); foci_filtered = "foci_ch"+channel+"_filtered"; Ext.CLIJ2_subtractImages(foci_blurred, foci_outliers_removed_blurred, foci_filtered); Ext.CLIJ2_release(foci_blurred); Ext.CLIJ2_release(foci_outliers_removed_blurred); if(debugMode) { print("Foci filter size XY: "+fociFilterSizeXY+", Z: "+fociFilterSizeZ); showImagefromGPU(foci_filtered); //if(bits == 8) run("16-bit"); } //TO DO: Check Labelmap_detected_foci for 8-bit images? (no values >255) if(processAllOtherImagesFixedThreshold == false) { threshold = medianNucleiStddev * thresholdMultiplier; //Set default threshold at n times the stddev (of the outlier-removed foci, which is a bit lower than the actual stddev, but less dependent on many foci being present) threshold = threshold*exp(thresholdFactor); print("Threshold set at "+d2s(thresholdMultiplier*exp(thresholdFactor),1)+" times above background stddev (bias set at "+thresholdFactor+"): Threshold used: "+threshold); } else print("Threshold fixed at "+threshold); maxDisplaySetting = minOf(pow(2,bits), threshold * 5); //logWindow = split(getInfo("log"),"\n"); //waitForUser(logWindow.length); if(fociDetectionMethod == "AreaMaxima local maximum detection") { Ext.CLIJ2_pull(foci_filtered); selectWindow(foci_filtered); // setBatchMode("show"); // hMin = 100; //setBatchMode(false); //H-watershed does not work in batch mode if this image is not shown. // run("H_Watershed", "impin=["+foci+"] hmin="+hMin+" thresh="+threshold+" peakflooding="+H_max_peak_flooding + " outputmask=false allowsplitting=false"); // print(getTitle); run("AreaMaxima local maximum detection (2D, 3D)", "minimum="+minFociSize+" threshold="+threshold/2); rename("Labelmap_detected_foci_ch"+channel); Ext.CLIJ2_push("Labelmap_detected_foci_ch"+channel); labeledSpots = "labeledSpots_ch"+channel; Ext.CLIJ2_reduceLabelsToCentroids("Labelmap_detected_foci_ch"+channel, labeledSpots); Ext.CLIJ2_mask(labeledSpots, labelmap_nuclei_3D_dilated, labeledSpots_masked); if(isOpen("foci_spots_ch"+channel)) close("foci_spots_ch"+channel); Ext.CLIJ2_pullBinary(labeledSpots_masked); rename("foci_spots_ch"+channel); Ext.CLIJ2_release(labeledSpots); Ext.CLIJ2_release(labeledSpots_masked); close(foci_filtered); close("Labelmap_detected_foci_ch"+channel); } else if(fociDetectionMethod == "Marker-controlled watershed (recommended)") { // Alternative function - could be faster: // Ext.CLIJx_detectAndLabelMaximaAboveThreshold(foci, detectedMaxima, 0, 0, 0, threshold, false) // Ext.CLIJ2_pull(detectedMaxima); // rename("detectedMaxima"); //Detect all maxima allMaxima = "allMaxima"; if(ThreeDHandling == "Use quasi-2D foci detection (detect foci separately in every Z-slice)" && gslices>1) Ext.CLIJ2_detectMaximaSliceBySliceBox(foci_filtered, allMaxima, fociSizeXY, fociSizeXY); else Ext.CLIJ2_detectMaxima3DBox(foci_filtered, allMaxima, fociSizeXY, fociSizeXY, fociSizeZ); //remove spots lower than threshold MaskFociAboveThreshold = "MaskFociAboveThreshold_ch"+channel; Ext.CLIJ2_threshold(foci_filtered, MaskFociAboveThreshold, threshold); Ext.CLIJ2_mask(allMaxima, MaskFociAboveThreshold, maskedSpots); Ext.CLIJ2_release(allMaxima); if(debugMode) showImagefromGPU(MaskFociAboveThreshold); labeledSpots = "labeledSpots_ch"+channel; Ext.CLIJ2_connectedComponentsLabelingBox(maskedSpots, labeledSpots); //Create label image with single pixels Ext.CLIJ2_release(maskedSpots); Ext.CLIJ2_mask(labeledSpots, labelmap_nuclei_3D_dilated, labeledSpots_masked); Ext.CLIJ2_release(labeledSpots); // if(bits == 8) Ext.CLIJ2_release(maskedSpots_16bit); if(debugMode) { showImagefromGPU(labeledSpots_masked); //run("16-bit"); run("glasbey on dark"); run("Merge Channels...", "c1="+foci_filtered+" c2="+labeledSpots_masked+" create keep"); rename("Overlay_spots_ch"+channel); Stack.setChannel(2); run("Magenta"); Stack.setChannel(1); run("Green"); setBatchMode("show"); } //run("16-bit"); //run("glasbey on dark"); //run("Merge Channels...", "c1="+foci+" c2="+labeledSpots+" create keep"); //setBatchMode("show"); //Ext.CLIJx_detectAndLabelMaximaAboveThreshold(foci, detectedMaxima, 0, 0, 0, threshold, false) //Ext.CLIJx_seededWatershed(detectedMaxima, foci, Image_label_map_destination, threshold); //Doesn't quite cut it - too many maxima detected? //Use CLIJx-MorphoLibJ marker-controlled watershed instead labelmap_foci_filled = "Labelmap_detected_foci_ch"+channel; //Invert foci before marker-controlled watershed! Ext.CLIJ2_invert(foci_filtered, foci_filtered_inverted); Ext.CLIJx_morphoLibJMarkerControlledWatershed(foci_filtered_inverted, labeledSpots_masked, MaskFociAboveThreshold, labelmap_foci); Ext.CLIJx_morphoLibJFillHoles(labelmap_foci, labelmap_foci_filled); Ext.CLIJ2_release(labelmap_foci); Ext.CLIJ2_release(foci_filtered); Ext.CLIJ2_release(foci_filtered_inverted); Ext.CLIJ2_release(MaskFociAboveThreshold); } labelmapDetectedFoci = "Labelmap_detected_foci_ch"+channel; //Exclude foci too small or too large if(debugMode) { showImagefromGPU(labelmapDetectedFoci); run("glasbey on dark"); resetMinAndMax; } Ext.CLIJ2_excludeLabelsOutsideSizeRange(labelmapDetectedFoci, labelmapDetectedFociFiltered, minFociSize, maxFociSize); if(debugMode) showImagefromGPU(labelmapDetectedFociFiltered); //Exclude foci outside of nuclei, or nuclei+range Ext.CLIJ2_closeIndexGapsInLabelMap(labelmapDetectedFociFiltered, labelmap_foci_final); Ext.CLIJ2_release(labelmapDetectedFociFiltered); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_foci_final, nrFoci); //Mask the spots image with the remaining foci Ext.CLIJ2_mask(labeledSpots_masked, labelmap_foci_final, labeledSpots_masked_filtered); Ext.CLIJ2_release(labeledSpots_masked); if(isOpen("foci_spots_ch"+channel)) close("foci_spots_ch"+channel); if(useLargeSpots == true) Ext.CLIJ2_dilateBox(labeledSpots_masked_filtered, labeledSpots_masked_filtered_dilated); else labeledSpots_masked_filtered_dilated = labeledSpots_masked_filtered; Ext.CLIJ2_pullBinary(labeledSpots_masked_filtered_dilated); rename("foci_spots_ch"+channel); Ext.CLIJ2_release(labeledSpots_masked_filtered); //Remove some lines in the log window (unnecessary info by MorphoLibJ) logWindowContents = getInfo("log"); logWindowLines = split(logWindowContents, "\n"); logWindowLines = Array.trim(logWindowLines, logWindowLines.length - 7); logWindowLines = arrayToString(logWindowLines); print("\\Clear"); print(logWindowLines); print("\\Update:Channel "+channel+": "+nrFoci+" foci detected in "+nrNuclei+" nuclei ("+d2s(nrFoci/nrNuclei,1)+" foci per nucleus)" + "\n"); Ext.CLIJ2_pull(labelmap_foci_final); if(debugMode) { run("glasbey on dark"); resetMinAndMax; setBatchMode("show"); } rename("Labelmap_detected_foci_filtered_ch"+channel); //setMinAndMax(0, max); //run("glasbey_on_dark"); //setBatchMode("show"); //create binary from the labelmap if(isOpen("foci_spots_ch"+channel)) close("Mask_foci_ch"+channel); Ext.CLIJ2_pullBinary(labelmap_foci_final); rename("Mask_foci_ch"+channel); Ext.CLIJ2_release(labelmap_foci_final); if(!debugMode) close(foci+"_outliers_removed"); return newArray("Labelmap_detected_foci_filtered_ch"+channel, "Mask_foci_ch"+channel, "foci_spots_ch"+channel, threshold); } function mergeOriginalAndDetection(original, nrNuclei, nuclei_outlines, foci_mask, foci_spots, channel, useROI, x_ROI, y_ROI) { // Merge the original foci channel with the detected foci if(useROI == false) { //prepare nuclei selectWindow("nuclei"); run(bits+"-bit"); //convert to allow merging if(cytoChannel>0 && cellpose==true) { rename("cells"); selectWindow(original); run("Duplicate...", "title=nuclei duplicate channels="+nucleiChannel); run(bits+"-bit"); //convert to allow merging } else if(cytoChannel>0) { selectWindow(original); run("Duplicate...", "title=cells duplicate channels="+cytoChannel); run(bits+"-bit"); //convert to allow merging } //prepare foci raw foci_RAW = "foci_ch"+channel; if(bits == 16) Ext.CLIJ2_convertUInt16(foci_RAW, foci_RAW_8or16bit); else if(bits == 8) Ext.CLIJ2_convertUInt8(foci_RAW, foci_RAW_8or16bit); else foci_RAW_8or16bit = foci_RAW; Ext.CLIJ2_pull(foci_RAW_8or16bit); if(bits == 8 || bits == 16) Ext.CLIJ2_release(foci_RAW_8or16bit); } if(useROI == true) { //prepare nuclei close("nuclei"); selectWindow(original); run("Select None"); run("Duplicate...", "title=nuclei duplicate channels="+nucleiChannel); run(bits+"-bit"); //convert to allow merging if(cytoChannel>0) { selectWindow(original); run("Duplicate...", "title=cells duplicate channels="+cytoChannel); run(bits+"-bit"); } //prepare foci raw selectWindow(original); foci_RAW_8or16bit = "foci_ch"+channel; run("Duplicate...", "title="+foci_RAW_8or16bit+ " duplicate channels="+channel); run(bits+"-bit"); //convert to allow merging } //prepare foci_mask - channel 3 selectWindow(foci_mask); run(bits+"-bit"); //convert to allow merging if(useROI == true) { Ext.CLIJ2_push(foci_mask); Ext.CLIJ2_create3D(foci_mask_large, gwidth, gheight, gslices, bits); Ext.CLIJ2_set(foci_mask_large, 0); Ext.CLIJ2_paste3D(foci_mask, foci_mask_large, x_ROI, y_ROI, 0); Ext.CLIJ2_pull(foci_mask_large); setMinAndMax(0, 255); } else foci_mask_large = foci_mask; //prepare foci_spots - channel 4 selectWindow(foci_spots); run(bits+"-bit"); //convert to allow merging if(useROI == true) { Ext.CLIJ2_push(foci_spots); Ext.CLIJ2_create3D(foci_spots_large, gwidth, gheight, gslices, bits); Ext.CLIJ2_set(foci_spots_large, 0); Ext.CLIJ2_paste3D(foci_spots, foci_spots_large, x_ROI, y_ROI, 0); Ext.CLIJ2_pull(foci_spots_large); setMinAndMax(0, 255); } else foci_spots_large = foci_spots; if(debugMode) { selectWindow("nuclei"); setBatchMode("show"); // selectWindow(foci_RAW_8or16bit); // setBatchMode("show"); selectWindow(foci_mask_large); setBatchMode("show"); selectWindow(foci_spots_large); setBatchMode("show"); if(isOpen("Cells")) selectWindow("cells"); setBatchMode("show"); } if(cytoChannel<=0) run("Merge Channels...", "c1=nuclei c2=" + foci_RAW_8or16bit + " c3=" + foci_mask_large + " c4=" + foci_spots_large + " create keep"); else run("Merge Channels...", "c1=nuclei c2=" + foci_RAW_8or16bit + " c3=" + foci_mask_large + " c4=" + foci_spots_large + " c5=cells create keep"); rename("Foci_overlay_ch"+channel); getDimensions(width, height, channels, slices, frames); Stack.setSlice(slices/2); Stack.setChannel(1); setLabel("Foci_overlay_ch"+channel, "NUCLEI"); setLut(b_reds, b_greens, b_blues); run("Enhance Contrast...", "saturated=0.02"); Stack.setChannel(2); setLabel("Foci_overlay_ch"+channel, "FOCI"); run("Green"); run("Enhance Contrast...", "saturated=0.02"); getMinAndMax(min, max); setMinAndMax(min, maxOf(max, min + maxDisplaySetting)); Stack.setChannel(3); setLabel("Foci_overlay_ch"+channel, "DETECTED FOCI"); resetMinAndMax; run("Magenta"); Stack.setChannel(4); setLabel("Foci_overlay_ch"+channel, "CENTROIDS OF DETECTED FOCI"); resetMinAndMax; run("Magenta"); if(cytoChannel>0) { Stack.setChannel(5); setLabel("Foci_overlay_ch"+channel, "CELLS"); run("Enhance Contrast...", "saturated=0.35"); setLut(o_reds, o_greens, o_blues); } //Add nuclei outlines as overlay to the original image selectWindow(nuclei_outlines); if(overlayBrightness == "dim") setMinAndMax(0, 512); //Silly way to set nuclei outline overlay to dim, but it works if(overlayBrightness == "bright") resetMinAndMax; //Silly way to set nuclei outline overlay to dim, but it works if(debugMode) setBatchMode("show"); selectWindow("Foci_overlay_ch"+channel); for (i = 1; i <= gslices; i++) { if(gslices>1 && nucleiSegmentationChoice == "Cellpose segmentation 3D") { selectImage(nuclei_outlines); Stack.setSlice(i); if(maxFociDistanceOutsideNuclei_setting != 0) { selectImage(nuclei_dilated_outlines); Stack.setSlice(i); } } if(gslices>1) { selectImage("Foci_overlay_ch"+channel); Stack.setSlice(i); } run("Add Image...", "image="+nuclei_outlines+" x=0 y=0 opacity="+LABELOPACITY+" zero"); //Add labelmap to image as overlay Overlay.setPosition(0, i, 0); if(maxFociDistanceOutsideNuclei_setting != 0) run("Add Image...", "image="+nuclei_dilated_outlines+" x=0 y=0 opacity="+maxOf(LABELOPACITY/4, 20)+" zero"); Overlay.setPosition(0, i, 0); } /* if(addNumbersOverlay) { run("Clear Results"); Ext.CLIJ2_statisticsOfLabelledPixels(labelmap_nuclei, labelmap_nuclei); //Somehow if you put (image, labelmap) as arguments the pixel count is wrong setFont("SansSerif", labelFontSize, "bold antialiased"); color = color_to_hex(fontColor); setColor(color); for (i = 0; i < nrNuclei; i++) { x = getResult("MASS_CENTER_X", i); y = getResult("MASS_CENTER_Y", i); Overlay.drawString(i+1, x - labelFontSize/2, y + labelFontSize/2); } } */ if(useROI) { run("Translate...", "x="+x_ROI+" y="+y_ROI+" interpolation=None overlay"); roiManager("Remove Slice Info"); run("From ROI Manager"); //Add the selected ROI as overlay } return "Foci_overlay_ch"+channel; } function setLabel(image, label) { selectWindow(image); getDimensions(width, height, channels, slices, frames); Stack.getPosition(ch, sl, fr); if(slices>1) { for (i=1; i<=slices; i++) { Stack.setSlice(i); //run("Set Label...", "label="+label.replace(" ","_")); setMetadata("Label", label.replace(" ","_")); } } else setMetadata("Label", label.replace(" ","_")); Stack.setPosition(ch, sl, fr); } function measureFoci(original, channel, nrNuclei, labelmap_nuclei_3D, labelmap_foci, boundingBox_X, boundingBox_Y, boundingBox_width, boundingBox_height) { Ext.CLIJ2_push(labelmap_foci); // if(debugMode) showImagefromGPU(labelmap_foci); selectWindow(original); run("Duplicate...", "title=foci_RAW_ch"+channel + " duplicate channels=" + channel); // nucleus_id_ = newArray(nrNuclei); // nucleus_area_ = newArray(nrNuclei); foci_count_ = newArray(nrNuclei); foci_mean_int_ = newArray(nrNuclei); foci_median_int_ = newArray(nrNuclei); foci_size_ = newArray(nrNuclei); foci_median_size_ = newArray(nrNuclei); foci_sum_size_ = newArray(nrNuclei); foci_sum_int_ = newArray(nrNuclei); nucleus_mean_int_ = newArray(nrNuclei); nucleus_sum_intensity_ = newArray(nrNuclei); background_ = newArray(nrNuclei); //Loop over all nuclei startTime = getTime(); foci_raw = "foci_RAW_ch"+channel; Ext.CLIJ2_push(foci_raw); for (i = 0; i < nrNuclei; i++) { showStatus("Analyzing foci in channel "+channel+", nucleus "+i+1 + "/" + nrNuclei); showProgress(i, nrNuclei); nrResults = nResults; labelmap_foci_nucleus = "labelmap_foci_nucleus"; labelmap_foci_nucleus_relabeled = "labelmap_foci_nucleus_relabeled"; foci_raw_cropped = "foci_nucleus"; // Crop labelmaps and raw images before measuring - much faster for large images Ext.CLIJ2_crop3D(foci_raw, foci_raw_cropped, boundingBox_X[i], boundingBox_Y[i], 0, boundingBox_width, boundingBox_height, gslices); if(debugMode) { // showImagefromGPU(foci_raw_cropped); // getDimensions(width, height, channels, slices, frames); // print("dimensions of nucleus "+i+1+" image: "+width+ ", "+height); } Ext.CLIJ2_crop3D(labelmap_nuclei_3D_dilated, labelmap_nuclei_cropped, boundingBox_X[i], boundingBox_Y[i], 0, boundingBox_width, boundingBox_height, gslices); Ext.CLIJ2_crop3D(labelmap_foci, labelmap_foci_cropped, boundingBox_X[i], boundingBox_Y[i], 0, boundingBox_width, boundingBox_height, gslices); Ext.CLIJ2_maskLabel(labelmap_foci_cropped, labelmap_nuclei_cropped, labelmap_foci_nucleus, i+1); //This is roughly equally fast as maskLabel: //Ext.CLIJ2_labelToMask(labelmap_nuclei_3D, mask_nucleus, i+1); //Ext.CLIJ2_mask(labelmap_foci, mask_nucleus, labelmap_foci_nucleus); //Also try: //Ext.CLIJ2_maskStackWithPlane(Image_source, Image_mask, Image_destination); //relabel foci in the current nucleus label and measure using MprpholibJ [faster than Ext.CLIJ2_statisticsOfLabelledPixels(foci_raw_cropped, labelmap_foci_nucleus_relabeled);] Ext.CLIJ2_closeIndexGapsInLabelMap(labelmap_foci_nucleus, labelmap_foci_nucleus_relabeled); Ext.CLIJ2_pull(foci_raw_cropped); Ext.CLIJ2_pull(labelmap_foci_nucleus_relabeled); run("Intensity Measurements 2D/3D", "input="+foci_raw_cropped+" labels="+labelmap_foci_nucleus_relabeled+" mean stddev max min median skewness volume"); fociTable = "foci_nucleus-intensity-measurements"; headings = split(Table.headings(fociTable), "\t"); //Get column headers. Ok, it is the same for every nucleus, but it's fast anyway. close(foci_raw_cropped); close(labelmap_foci_nucleus_relabeled); //Count the number of foci Ext.CLIJ2_getMaximumOfAllPixels(labelmap_foci_nucleus_relabeled, nrFoci); foci_count_[i] = nrFoci; //Copy measured data to allFociResults table totalNrFoci = Table.size(allFociResultsTable); for (k = totalNrFoci; k < totalNrFoci + nrFoci; k++) { Table.set("Cell ID", k, i+1, allFociResultsTable); Table.set("channel", k, channel, allFociResultsTable); } for (col=0; col0) { int_ = Table.getColumn("Mean", fociTable); foci_mean_int_[i] = meanOfArray(int_); foci_median_int_[i] = medianOfArray(int_); foci_sum_int_[i] = sumArray(int_); size_ = Table.getColumn("Volume", fociTable); foci_size_[i] = meanOfArray(size_) * pixelWidth * pixelHeight * pixelDepth; foci_median_size_[i] = medianOfArray(size_) * pixelWidth * pixelHeight * pixelDepth; foci_sum_size_[i] = foci_size_[i] * nrFoci; //average size * nr } else { foci_mean_int_[i] = NaN; foci_median_int_[i] = NaN; foci_sum_int_[i] = NaN; foci_size_[i] = NaN; foci_median_size_[i] = NaN; foci_sum_size_[i] = NaN; } } Table.update(allFociResultsTable); Ext.CLIJ2_release(foci_raw_cropped); Ext.CLIJ2_release(labelmap_nuclei_cropped); Ext.CLIJ2_release(labelmap_foci_cropped); Ext.CLIJ2_release(labelmap_foci_nucleus_relabeled); Ext.CLIJ2_release(labelmap_foci_nucleus); if(debugMode) print("Measuring all foci in channel "+channel+" took "+ getTime() - startTime +" ms."); //Measure overall nuclear intensity run("Clear Results"); Ext.CLIJ2_statisticsOfBackgroundAndLabelledPixels(foci_raw, labelmap_nuclei_3D); background = getResult("MEAN_INTENSITY", 0); Array.fill(background_, background); for (i = 0; i < nResults-1; i++) { nucleus_mean_int_[i] = Table.get("MEAN_INTENSITY", i+1, "Results"); nucleus_sum_intensity_[i] = Table.get("SUM_INTENSITY", i+1, "Results"); } Ext.CLIJ2_release(foci_raw); //Add relevant results to the resultTable Table.setColumn("Background intensity ch"+channel, background_, resultTable); Table.setColumn("Mean cell intensity ch"+channel, nucleus_mean_int_, resultTable); Table.setColumn("Sum cell intensity ch"+channel, nucleus_sum_intensity_, resultTable); Table.setColumn("Foci count ch"+channel, foci_count_, resultTable); Table.setColumn("Mean foci intensity ch"+channel, foci_mean_int_, resultTable); Table.setColumn("Median foci intensity ch"+channel, foci_median_int_, resultTable); if(gslices > 1) Table.setColumn("Mean foci volume ("+unit+"^3) ch"+channel, foci_size_, resultTable); else Table.setColumn("Mean foci area ("+unit+"^2) ch"+channel, foci_size_, resultTable); if(gslices > 1) Table.setColumn("Median foci volume ("+unit+"^3) ch"+channel, foci_median_size_, resultTable); else Table.setColumn("Median foci area ("+unit+"^2) ch"+channel, foci_median_size_, resultTable); if(gslices > 1) Table.setColumn("Total foci volume ("+unit+"^3) ch"+channel, foci_sum_size_, resultTable); else Table.setColumn("Total foci area ("+unit+"^2) ch"+channel, foci_sum_size_, resultTable); Table.setColumn("Total foci intensity ch"+channel, foci_sum_int_, resultTable); Table.update(); close(fociTable); return nrFoci; } function measure_nuclear_intensities(original, nrNuclei, labelmap_nuclei_3D, gchannels, fociChannelA, fociChannelB, resultTable) { background_ = newArray(nrNuclei); for(c=1; c<=gchannels; c++) { if(c == fociChannelA || detectFociChannelB == true && c == fociChannelB) continue; else { selectWindow(original); run("Clear Results"); Stack.setChannel(c); channel_to_measure = original+"_ch"+c; run("Duplicate...", "title=[" + channel_to_measure + "] duplicate channels=" + c); Ext.CLIJ2_push(channel_to_measure); Ext.CLIJ2_statisticsOfBackgroundAndLabelledPixels(channel_to_measure, labelmap_nuclei_3D); mean_intensity_ = Table.getColumn("MEAN_INTENSITY", "Results"); sum_intensity_ = Table.getColumn("SUM_INTENSITY", "Results"); background = mean_intensity_[0]; Array.fill(background_, background); mean_intensity_ = Array.slice(mean_intensity_, 1, mean_intensity_.length); sum_intensity_ = Array.slice(sum_intensity_, 1, sum_intensity_.length); Table.setColumn("Background intensity ch"+c, background_, resultTable); Table.setColumn("Mean cell intensity ch"+c, mean_intensity_, resultTable); Table.setColumn("Sum cell intensity ch"+c, sum_intensity_, resultTable); Table.update; Ext.CLIJ2_release(channel_to_measure); } } } //Draw nuclei and/or foci numbers as overlay function overlay_numbers_on_image(overlay_image) { selectImage(overlay_image); if(addNumbersOverlay) { setFont("SansSerif", labelFontSize, "bold antialiased"); color1 = color_to_hex(fontColor); color2 = "#ff0051"; setColor(color1); if(overlayChoice == "nucleus/cell ID") { for (i = 0; i < nrNuclei; i++) { x = getResult("CENTROID_X", i); y = getResult("CENTROID_Y", i); z_start = getResult("BOUNDING_BOX_Z", i); z_end = getResult("BOUNDING_BOX_END_Z", i); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(i+1, x - labelFontSize/2 + x_ROI, y + labelFontSize/2 + y_ROI); Overlay.setPosition(0, k+1, 0); } } } else if(overlayChoice == "foci count") { for (i = 0; i < nrNuclei; i++) { x = getResult("CENTROID_X", i); y = getResult("CENTROID_Y", i); z_start = getResult("BOUNDING_BOX_Z", i); z_end = getResult("BOUNDING_BOX_END_Z", i); nrf = Table.get("Foci count ch"+fociChannelA, i, resultTable); setColor(color2); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(nrf, x - labelFontSize/2, y + labelFontSize/2); Overlay.setPosition(0, k+1, 1); } if(detectFociChannelB) { nrf = Table.get("Foci count ch"+fociChannelB, i, resultTable); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(nrf, x - labelFontSize/2, y + labelFontSize/2); Overlay.setPosition(0, k+1, 2); } } } } else if(overlayChoice == "both") { for (i = 0; i < nrNuclei; i++) { x = getResult("CENTROID_X", i); y = getResult("CENTROID_Y", i) - labelFontSize/2 - 2; z_start = getResult("BOUNDING_BOX_Z", i); z_end = getResult("BOUNDING_BOX_END_Z", i); setColor(color1); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(i+1, x - labelFontSize/2, y + labelFontSize/2); Overlay.setPosition(0, k+1, 0); } y = getResult("CENTROID_Y", i) + labelFontSize/2 + 2; nrf = Table.get("Foci count ch"+fociChannelA, i, resultTable); setColor(color2); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(nrf, x - labelFontSize/2, y + labelFontSize/2); Overlay.setPosition(0, k+1, 1); } if(detectFociChannelB) { nrf = Table.get("Foci count ch"+fociChannelB, i, resultTable); for(k = z_start; k <= z_end ; k++) { Overlay.drawString(nrf, x - labelFontSize/2, y + labelFontSize/2); Overlay.setPosition(0, k+1, 2); } } } } updateDisplay(); } } function computeOverlap(mask_fociA, mask_fociB, nrNuclei, labelmap_nuclei_3D, boundingBox_X, boundingBox_Y, boundingBox_width, boundingBox_height) { showStatus("Computing foci colocalization..."); foci_overlap_map = "Foci overlap map"; startTime = getTime(); run("Clear Results"); Ext.CLIJ2_push(mask_fociA); Ext.CLIJ2_push(mask_fociB); Ext.CLIJ2_binaryIntersection(mask_fociA, mask_fociB, overlap); Ext.CLIJ2_connectedComponentsLabelingBox(overlap, labelmap_overlap); //Ext.CLIJ2_generateBinaryOverlapMatrix(labelmap_fociA, labelmap_fociB, overlap_matrix); //Ext.CLIJ2_pull(overlap_matrix); //setBatchMode("show"); Ext.CLIJ2_getJaccardIndex(mask_fociA, mask_fociB, jaccardIndex); //print("Jaccard Index (before overlap size exclusion): "+jaccardIndex); Ext.CLIJ2_release(mask_fociA); Ext.CLIJ2_release(mask_fociB); Ext.CLIJ2_release(overlap); //Filter on overlap size labelmap_overlap_filtered = "labelmap_overlap_filtered"; Ext.CLIJ2_excludeLabelsOutsideSizeRange(labelmap_overlap, labelmap_overlap_filtered, minOverlapSize, maxFociSize); Ext.CLIJ2_release(labelmap_overlap); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_overlap_filtered, max); print("\n"+max+" colocalizing foci with overlap > "+ minOverlapSize +" voxels detected in "+nrNuclei+" nuclei ("+d2s(max/nrNuclei,1)+" per nucleus)."); //Count overlapping foci for all nuclei overlap_count_ = newArray(nrNuclei); overlap_volume_ = newArray(nrNuclei); for (i = 0; i < nrNuclei; i++) { showStatus("Computing foci overlap in nucleus "+i+1 + "/" + nrNuclei); showProgress(i, nrNuclei); labelmap_overlap_nucleus = "labelmap_overlap_nucleus_"+i+1; // Crop labelmaps (nuclei and overlap) before measuring - much faster for large images Ext.CLIJ2_crop3D(labelmap_nuclei_3D, labelmap_nuclei_cropped, boundingBox_X[i], boundingBox_Y[i], 0, boundingBox_width, boundingBox_height, gslices); Ext.CLIJ2_crop3D(labelmap_overlap_filtered, labelmap_overlap_cropped, boundingBox_X[i], boundingBox_Y[i], 0, boundingBox_width, boundingBox_height, gslices); if(debugMode) { // showImagefromGPU(labelmap_nuclei_cropped); // getDimensions(width, height, channels, slices, frames); // print("dimensions of nucleus "+i+1+" image: "+width+ ", "+height); } Ext.CLIJ2_maskLabel(labelmap_overlap_cropped, labelmap_nuclei_cropped, labelmap_overlap_nucleus, i+1); Ext.CLIJ2_closeIndexGapsInLabelMap(labelmap_overlap_nucleus, labelmap_overlap_nucleus_relabeled); Ext.CLIJ2_getMaximumOfAllPixels(labelmap_overlap_nucleus_relabeled, nrOverlappingFoci); Ext.CLIJ2_countNonZeroPixels(labelmap_overlap_nucleus_relabeled); overlap_count_[i] = nrOverlappingFoci; overlap_volume_[i] = getResult("CountNonZero", i) * pixelWidth * pixelHeight * pixelDepth; //Cannot re-use the memory, because the shape is different for every nucleus! Ext.CLIJ2_release(labelmap_nuclei_cropped); Ext.CLIJ2_release(labelmap_overlap_cropped); Ext.CLIJ2_release(labelmap_overlap_nucleus); Ext.CLIJ2_release(labelmap_overlap_nucleus_relabeled); } Ext.CLIJ2_release(labelmap_overlap_filtered); //Add overlap results to the resultTable Table.setColumn("Overlapping foci count", overlap_count_, resultTable); if(gslices > 1) Table.setColumn("Overlapping foci volume ("+unit+"^3)", overlap_volume_, resultTable); else Table.setColumn("Overlapping foci area ("+unit+"^2)", overlap_volume_, resultTable); //Channel A fociCountChannelA_ = Table.getColumn("Foci count ch"+fociChannelA, resultTable); overlap_count_fraction_A_ = divideArrays(overlap_count_, fociCountChannelA_); overlap_count_percentage_A_ = multiplyArraywithScalar(overlap_count_fraction_A_, 100); Table.setColumn("Overlap count % (ch"+fociChannelA+")", overlap_count_percentage_A_, resultTable); if(gslices > 1) fociVolumeChannelA_ = Table.getColumn("Total foci volume ("+unit+"^3) ch"+fociChannelA, resultTable); else fociVolumeChannelA_ = Table.getColumn("Total foci area ("+unit+"^2) ch"+fociChannelA, resultTable); overlap_volume_fraction_A_ = divideArrays(overlap_volume_, fociVolumeChannelA_); overlap_volume_percentage_A_ = multiplyArraywithScalar(overlap_volume_fraction_A_, 100); if(gslices > 1) Table.setColumn("Overlap volume % (ch"+fociChannelA+")", overlap_volume_percentage_A_, resultTable); else Table.setColumn("Overlap area % (ch"+fociChannelA+")", overlap_volume_percentage_A_, resultTable); //Channel B fociCountChannelB_ = Table.getColumn("Foci count ch"+fociChannelB, resultTable); overlap_count_fraction_B_ = divideArrays(overlap_count_, fociCountChannelB_); overlap_count_percentage_B_ = multiplyArraywithScalar(overlap_count_fraction_B_, 100); Table.setColumn("Overlap count % (ch"+fociChannelB+")", overlap_count_percentage_B_, resultTable); if(gslices > 1) fociVolumeChannelB_ = Table.getColumn("Total foci volume ("+unit+"^3) ch"+fociChannelB, resultTable); else fociVolumeChannelB_ = Table.getColumn("Total foci area ("+unit+"^2) ch"+fociChannelB, resultTable); overlap_volume_fraction_B_ = divideArrays(overlap_volume_, fociVolumeChannelB_); overlap_volume_percentage_B_ = multiplyArraywithScalar(overlap_volume_fraction_B_, 100); if(gslices > 1) Table.setColumn("Overlap volume % (ch"+fociChannelB+")", overlap_volume_percentage_B_, resultTable); else Table.setColumn("Overlap area % (ch"+fociChannelB+")", overlap_volume_percentage_B_, resultTable); // TO DO (maybe): Use CLIJx_labelOverlapCountMap to calculate the number of overlapping foci. // It will replace the above code, except for the overlap size filtering, which may not be possible. // Also, you can use it to count foci in every cell by passing the nuclei labelmap and a foci labelmap. //Create Results table column for the overlap visualization, including the background (entry 0) // run("Clear Results"); // setResult("Overlap percentage", 0, 0); // for (i = 0; i < nrNuclei; i++) { // setResult("Overlap percentage", i+1, overlap_count_percentage_[i]); // } // updateResults(); //Visualize the amount of overlapping foci in a map //The following line doesn't work. Why not? // Ext.CLIJ2_generateParametricImageFromResultsTableColumn(labelmap_nuclei_3D, overlap_count_map, "Overlap percentage"); // Ext.CLIJ2_pushResultsTableColumn(overlap_fill, "Overlap percentage ("); overlap_count_map_A = create_overlap_map(overlap_count_percentage_A_, labelmap_nuclei_3D, "overlap_count_percentage_map_A"); overlap_volume_map_A = create_overlap_map(overlap_volume_percentage_A_, labelmap_nuclei_3D, "overlap_volume_percentage_map_A"); overlap_count_map_B = create_overlap_map(overlap_count_percentage_B_, labelmap_nuclei_3D, "overlap_count_percentage_map_B"); overlap_volume_map_B = create_overlap_map(overlap_volume_percentage_B_, labelmap_nuclei_3D, "overlap_volume_percentage_map_B"); if(debugMode) print("Measuring colocalization took "+ getTime() - startTime +" ms."); //Create foci mask overlap image run("Merge Channels...", "c1="+mask_fociA+" c2="+mask_fociB+" c3="+overlap_count_map_A+" c4="+overlap_volume_map_A+" c5="+overlap_count_map_B+" c6="+overlap_volume_map_B+" create"); rename(foci_overlap_map); Stack.setChannel(1); setLabel(foci_overlap_map, "FOCI ch"+fociChannelA); run("Green"); Stack.setChannel(2); setLabel(foci_overlap_map, "FOCI ch"+fociChannelB); run("Red"); Stack.setChannel(3); setLabel(foci_overlap_map, "FOCI COUNT OVERLAP PERCENTAGE WITH ch"+fociChannelA); setLut(b_reds, b_greens, b_blues); setMinAndMax(0, 100); Stack.setChannel(4); setLabel(foci_overlap_map, "FOCI VOLUME OVERLAP PERCENTAGE WITH ch"+fociChannelA); run("Magenta"); setMinAndMax(0, 100); Stack.setChannel(5); setLabel(foci_overlap_map, "FOCI COUNT OVERLAP PERCENTAGE WITH ch"+fociChannelB); setLut(b_reds, b_greens, b_blues); setMinAndMax(0, 100); Stack.setChannel(6); setLabel(foci_overlap_map, "FOCI VOLUME OVERLAP PERCENTAGE WITH ch"+fociChannelB); run("Magenta"); setMinAndMax(0, 100); Stack.setDisplayMode("color"); setBatchMode("show"); //Add nuclei outlines as overlay to the overlap image selectImage(foci_overlap_map); for (c = 1; c <= 6; c++) { Stack.setChannel(c); for (i = 1; i <= gslices; i++) { if(gslices > 1) Stack.setSlice(i); run("Add Image...", "image="+nuclei_outlines+" x=0 y=0 opacity="+LABELOPACITY+" zero"); //Add labelmap to image as overlay } } if(addNumbersOverlay) { run("Clear Results"); Ext.CLIJ2_statisticsOfLabelledPixels(labelmap_nuclei, labelmap_nuclei); setFont("SansSerif", labelFontSize, "bold antialiased"); color = color_to_hex(fontColor); setColor(color); if(overlayChoice == "nucleus/cell ID") { for (i = 0; i < nrNuclei; i++) { x = getResult("MASS_CENTER_X", i); y = getResult("MASS_CENTER_Y", i); Overlay.drawString(i+1, x - labelFontSize/2, y + labelFontSize/2); } } else if(overlayChoice == "foci count") { for (i = 0; i < nrNuclei; i++) { x = getResult("MASS_CENTER_X", i); y = getResult("MASS_CENTER_Y", i); Overlay.drawString(Table.get("Foci count ch"+fociChannelA, i, resultTable), x - labelFontSize/2, y + labelFontSize/2); } } //TO DO: FIX THIS ALSO FOR BOTH< AND POSSIBLY USE THE FUNCTION overlay_numbers_on_image() } if(gslices > 1) Stack.setSlice(gslices/2); setBatchMode("show"); return foci_overlap_map; } function create_overlap_map(overlap_percentage_, labelmap_nuclei_3D, overlapMapName) { overlap_percentage_ = insertElementIntoArrayAtPosition(0, overlap_percentage_, 0); //Insert element 0 Ext.CLIJ2_pushArray(overlap_fill, overlap_percentage_, overlap_percentage_.length, 1, 1); if(debugMode == true) showImagefromGPU(overlap_fill); overlap_map = "overlap_map"; Ext.CLIJ2_replaceIntensities(labelmap_nuclei_3D, overlap_fill, overlap_map); Ext.CLIJ2_release(overlap_fill); if(bits == 16) Ext.CLIJ2_convertUInt16(overlap_map, overlap_map_8or16bit); else if(bits == 8) Ext.CLIJ2_convertUInt8(overlap_map, overlap_map_8or16bit); Ext.CLIJ2_release(overlap_map); if(bits == 16) Ext.CLIJ2_replaceIntensity(overlap_map_8or16bit, overlapMapName, 65535, 0); else if(bits == 8) Ext.CLIJ2_replaceIntensity(overlap_map_8or16bit, overlapMapName, 255, 0); Ext.CLIJ2_pull(overlapMapName); rename(overlapMapName); if(debugMode == true) { run("Duplicate...", "title="+overlapMapName+"_debug"); setBatchMode("show"); } Ext.CLIJ2_release(overlapMapName); Ext.CLIJ2_release(overlap_map_8or16bit); return overlapMapName; } //Create a nice azure blue LUT for the nuclei function create_azure_lut(b_reds, b_greens, b_blues) { for (i=0; i<256; i++) { b_reds[i] = 0; b_greens[i] = 0.5*i; b_blues[i] = i; } } //Create an orange LUT for the cells function create_orange_lut(o_reds, o_greens, o_blues) { for (i=0; i<256; i++) { o_reds[i] = i; o_greens[i] = 0.5*i; o_blues[i] = 0; } } //Convert a color into a hexadecimal code function color_to_hex(color) { colorArray = split(color,",,"); hexcolor = "#" + IJ.pad(toHex(colorArray[0]),2) + IJ.pad(toHex(colorArray[1]),2) + IJ.pad(toHex(colorArray[2]),2); return hexcolor; } //Converts an array into a string, elements separated by "\n" function arrayToString(array) { outputString = ""; for (i = 0; i < array.length; i++) { outputString += array[i] + "\n"; } return outputString; } //Prepends the value to the array function prependToArray(value, array) { temparray=newArray(lengthOf(array)+1); for (i=0; i\");" + param + ";" ); } //Gets a persistent Script Parameter value. N.B. This returns 'null' when the parameter is set to the default value! function getScriptParameterValue(string_parameter) { return eval("js", "var ctx = Packages.ij.IJ.runPlugIn(\"org.scijava.Context\", \"\");" + "var ps = ctx.service(Packages.org.scijava.prefs.PrefService.class);" + "var " + string_parameter + " = ps.get(org.scijava.script.ScriptModule.class, \"" + string_parameter + "\", \"\");" + string_parameter + ";" ); } //Sets a persistent Script Parameter value. N.B. This returns 'null' when the parameter is set to the default value! function setScriptParameterValue(string_parameter, string_value) { eval("js", "var ctx = Packages.ij.IJ.runPlugIn(\"org.scijava.Context\", \"\");" + "var ps = ctx.service(Packages.org.scijava.prefs.PrefService.class);" + "var " + string_parameter + " = ps.put(org.scijava.script.ScriptModule.class, \"" + string_parameter + "\", \"" + string_value + "\");" ); } function cleanup() { // close("nuclei"); close("foci"); }