Raw File
Tip revision: 5b1a3920afa8e85132c94bcc6dfce94575f939ce authored by lauracmurphy on 08 April 2022, 11:27:04 UTC
Tip revision: 5b1a392
* Erosion Territories                                                                                                                                                                                                                                                                                                                                      
* Expected Input:
* The expected input for this script are 2D non-timelapse images in tiff/tif/ND2 format that are either 3 or 4 channel 
* One channel is expect to be a DAPI or other nuclear stain
* At start of macro, user selects folder with images in (nested folders dealth with fine), sets channel order and chooses desired number of bins
* Description of script workflow:
* The DAPI channel is used to segment the nuclei (using gaussian blur and huang threshold and the 'watershed irregular features' from BioVoxxel - see installation notes)
* Each nucleus area is duplicated from original image to isolate nuclei for further processing
* For each nuclei image, the area is measured and the expected area of the bins is calculated based on the number requsted
* Using the expected areas of these bins and iteratively thresholding on a distance map of the nucleus ROI bins are produced
* Mean intensity of these ROIs is measured in all the channels                                                        
* Output:
* Each nucleus found in each image will be saved as an RGB image with the outline of created bins overlaid on the image
* A .csv file will be saved with the areas of the nuclei, the bins and the mean intensity of each bin for each channel in the image
* Bin/ territory naming in the results go from outside in so bin 1 is the outermost one.
* Installation Notes: 
* Download FIJI here -
* Details of the Biovoxxel Update site -
* How to add an update site -                                                        
* How to use an ImageJ macro from Github -
* This script is hosted here with an example image -       
* Written by Laura Murphy (                                                                                                                                                                                                                                                                                                                             
* IGC Advanced Imaging Resource (
* First written: November 2019  Last updated January 2022                                                                                                        

//-- Part 0: Preparation steps: get directories from user, set up arrays to store results that are being added to and some other setup

// Get user to select folder with images for input and the output folder where they want results saved
input_folder = getDirectory("Parent directory of your image folders");
output_folder = getDirectory("Choose the folder where you want to save your results");

// Set up empty arrays that will be filled during script and displayed and saved at the end
Image_Name = newArray();
Nucleus_No = newArray();
Nucleus_Area = newArray();
Bin_No = newArray()
Bin_Area = newArray()
DAPI_MeanIntensity = newArray();
FITC_MeanIntensity = newArray();
TxRd_MeanIntensity = newArray();
Cy5_MeanIntensity = newArray();

// Magic numbers for nuclei segmentation
filter_sigma = 2; // Sigma for Gaussian blur
erosion_cycle = 1; // Erosion cycle for watershed
watershed_convexity = 0.99; // Convexity threshold for watershed
min_particle_size = 50; // Minimum particle size for Analyze Particles
max_particle_size = 700; // Maxmimum particle size for Analyze Particles
min_particle_circularity = 0.7; // Minimum circularity for Analyze Particles
max_particle_circularity = 1.00; // Maxmimum circularity for Analyze Particles

// Array of channel colours
colour_array = newArray("Blue", "Green", "Red", "Cyan");

// Get list of images in nested folders with certain extensions  - function at bottom of macro
directory_list = newArray();
directory_list = get_file_tree(input_folder, directory_list);

//-- Part 1: Image loop starts, each image goes through nuclear segmentation and each nuclei is duplicated from original image

// Start looping through images in input folders
for (img = 0; img < directory_list.length; img++) {
	path = directory_list[img];
	run("Bio-Formats Importer", "open=[" + path + "] autoscale color_mode=Grayscale rois_import=[ROI manager] view=Hyperstack stack_order=XYCZT");
	getDimensions(width, height, channels, slices, frames);

	original_image = getTitle();
	image_name = File.nameWithoutExtension;

	// When first image is open, get user to input desired number of bins and assign channels in order of the stack
	if (img == 0) {
		Dialog.create("User input");
		Dialog.addMessage("Enter the number of concentric territories you want the nucleus split into");
		Dialog.addNumber("Bins", 4);
		Dialog.addMessage("Enter the order of your channels below");
		Dialog.addNumber("DAPI", 4);
		Dialog.addNumber("FITC", 1);
		Dialog.addNumber("TxRd", 2);
		Dialog.addNumber("Cy5", 3);;
		//Store user input
		bins = Dialog.getNumber;
		channel_order = newArray(channels);
		for (ch = 0; ch < channels; ch++){
			channel_order[ch] = Dialog.getNumber;

	// Get pixel sizing, if image isn't scaled then skip image
	getVoxelSize(width, height, depth, unit);
	if (unit == "pixels"){
		waitForUser("Image has no scaling and will be skipped");

	// Set channel colours 
	for (lut = 0; lut < channel_order.length; lut++){

	// Duplicate image to preserve for display use 
	run("Duplicate...", "title=Display duplicate");

	// Segment nuclei in image
	run("Duplicate...", "title=Nuclear_Mask");
	run("Gaussian Blur...", "sigma=" + filter_sigma);
	setAutoThreshold("Huang dark");
	run("Convert to Mask"); 
	run("Watershed Irregular Features", "erosion=" + erosion_cycle + " convexity_threshold=" + watershed_convexity + " separator_size=0-Infinity");
	run("Analyze Particles...", "size=" + min_particle_size + "-" + max_particle_size + " circularity=" + min_particle_circularity + "-" + max_particle_circularity + " exclude include add");

	// Store number of nuclei
	number_of_nuclei = roiManager("Count");

	// Loop through nuclei, duplicate image of each one in image, naming them so they can be processed further	
	for (nuclei = 0; nuclei < number_of_nuclei; nuclei++){
		roiManager("Select", nuclei);
		title = "Nuclei_" + nuclei + 1;
		run("Duplicate...", "title=" + title);
		roiManager("Select", nuclei);
		run("Duplicate...", "duplicate");
		rename("Measure" + title);
		roiManager("Select", nuclei);
		run("Duplicate...", "duplicate");
		rename("Display" + title);

//-- Part 2: Looping through isolated nuclei image starts, creation of bins

	// Start of loop for going through the duplicated images of the solo nuclei
	for (image_of_nucleus = 0; image_of_nucleus < number_of_nuclei; image_of_nucleus++){
		title = "Nuclei_" + image_of_nucleus + 1;

		// Create distance map to use for erosion
		run("Distance Map");
		run("Clear Outside");

		// Get area of whole nuclei and maximum value of distance map (further distance in nuclei mask from edge)
		getStatistics(whole_area, mean, min, whole_max, std, histogram);
		nuclear_area = whole_area;
		run("Select None");

		// Work out the desired area each bin should be by dividing total nuclear area by number of bins
		area_fraction = whole_area/bins;
		prebins = newArray(bins);

		// Work out the expected areas of the shells (i.e. the smaller circles that share centroid with nucleus)
		for (areas_of_prebins = 0; areas_of_prebins < prebins.length; areas_of_prebins++){
			addtoArray = whole_area;
			prebins[areas_of_prebins] = addtoArray;
			whole_area = whole_area-area_fraction;

		// For each of the chosen bin number, iteratively reduce threshold of distance map until area of thresholded area is lower/equal to expected
		// Once that is reached, add to ROI Manager.
		for (bin = 0; bin < prebins.length; bin++) {
			current_area_aim = prebins[bin];
			for (threshold = 1; threshold < whole_max; threshold++) {
				getStatistics(area, mean, min, max, std, histogram);
				run("Select None");
				setThreshold(threshold, whole_max);
				run("Create Selection");
				getStatistics(area, mean, min, max, std, histogram);
				if (area <= current_area_aim){
		// Use prebins to create bins by using XOR function between the pre-bins in the ROI Manager
		number_of_rois = roiManager("Count");
		to_be_deleted = newArray(number_of_rois);

		// Add current ROIs to array to delete later
		for (bins = 0; bins < number_of_rois; bins++){
			to_be_deleted[bins] = bins;

		// Loop through rois, selecting pairs to get XOR area between them, add to ROI and name them.
		for (roi = 0; roi < number_of_rois-1; roi++){
			roiManager("select", newArray(roi, roi+1));
			Roi.setName("BIN_" + roi+1);
			roiManager("select", newArray(roi, roi+1));
		roiManager("select", number_of_rois-1);
		Roi.setName("BIN_" + number_of_rois);

		// Delete the shells from earlier
		roiManager("select", to_be_deleted);

//-- Part 3: Time to save results, for each bin, the image name, nucleus and nucleus area is saved to keep arrays same length
		//   Then for each of the bins (also called territories), the area is measured and then the intensity of each of the channels
		//   All things are saved to arrays which are displayed once complete. On the last image, the arrays are saved to the output folder

		// Loop through final shapes, adding nuclei results to arrays to keep arrays even
		for (territories = 0; territories < roiManager("Count"); territories++){                        	
			Image_Name = Array.concat(Image_Name, image_name);
			Nucleus_No = Array.concat(Nucleus_No, image_of_nucleus+1);
			Nucleus_Area = Array.concat(Nucleus_Area, nuclear_area);
			Bin_No = Array.concat(Bin_No, territories+1);    
			selectWindow("Measure" + title);
			roiManager("Select", territories);
			getStatistics(bin_area, mean, min, max, std, histogram);
			Bin_Area = Array.concat(Bin_Area, bin_area);

			roiManager("Select", territories);
			getStatistics(bin_area, mean_DAPI, min, max, std, histogram);
			DAPI_MeanIntensity = Array.concat(DAPI_MeanIntensity, mean_DAPI);
			roiManager("Select", territories);
			getStatistics(bin_area, mean_FITC, min, max, std, histogram);
			FITC_MeanIntensity = Array.concat(FITC_MeanIntensity, mean_FITC);
			roiManager("Select", territories);
			getStatistics(binArea, mean_TXRD, min, max, std, histogram);
			TxRd_MeanIntensity = Array.concat(TxRd_MeanIntensity, mean_TXRD);

			// If the images are 4 channel, then save the Cy5 channel
			if (channel_order.length > 3){
				roiManager("Select", territories);
				getStatistics(binArea, mean_Cy5, min, max, std, histogram);
				Cy5_MeanIntensity = Array.concat(Cy5_MeanIntensity, mean_Cy5);

		// Take original image, convert to RGB and save with bin boundaries overlayed on image for checking.
		selectWindow("Display" + title);
		run("Make Composite");
		run("RGB Color");
		roiManager("Show all without labels");                  
		saveAs("Tiff", output_folder + image_name + "_" + title);

	// Set up for next image coming up by closing open images and reset ROI manager and results table
	run("Close All");
	run("Clear Results");
	// Show updated array of results each time an image is completed. If image was 3 channel, the Cy5 array will be empty., Nucleus_No, Nucleus_Area, Bin_No, Bin_Area, DAPI_MeanIntensity, FITC_MeanIntensity, TxRd_MeanIntensity, Cy5_MeanIntensity);

// Save table of arrays to user directed output folder now all images are processed
saveAs("Results", output_folder + File.separator + "Results.csv");

//-- Part 4: Let user know the macro is complete

Dialog.addMessage("Macro Complete!");;

//-- Functions

function get_file_tree(directory, file_tree){
file_list = getFileList(directory);

for(f = 0; f < file_list.length; f++){
	if (matches(file_list[f], "(?i).*\\.(tif|tiff|nd2)$"))
		file_tree = Array.concat(file_tree, directory + file_list[f]);
	if(File.isDirectory(directory + File.separator + file_list[f]))
		file_tree = get_file_tree(directory + file_list[f], file_tree);
return file_tree;
back to top