// Written by Cenna van Manen.


// This macros is for analyzing the IFT accumulations per image, based on the user defined variables.

// As input, the macro assumes all images in the chosen folders are single layered tif files for which the ciliary, basal body and IFT
// channels are stored in separate folders, but with the same name.If required, these files and folders can be generated using the channel splitter macro.
// Furthermore, make sure the files and foldername do not contain comma's (,) or semicolons (;). Last, the indicated folder should not contain any other files

// As output a control image and an excel file per photo will be generated, as well as, a summary of all measured IFT accumulations.
// Furthermore, three text files will be generated, one containing the macro, one containing the log of the run, and the other contiainging the used variables.



print ("Starting IFT accumulation analysis");
requires("1.52n")

//Determine if the macro is called upon by another otherwise the vars will be gathered here
//Data format data_from_other_macro: input folder, general output folder, channel name list, channel names, image type, is silent mode
data_from_other_macro = getArgument();
if (lengthOf(data_from_other_macro)>1){
	stand_alone = 0;
	data_from_other_macro = split(data_from_other_macro, ",");
	//Sets silent mode based on previously indicated preferences
	silent_mode = data_from_other_macro[1];
	if(silent_mode){
		setBatchMode(1);
	}
	else {
		setBatchMode(0);
	}
	parent_input_folder = File.getParent(data_from_other_macro[0]);
	input_folder_IFT = parent_input_folder  + File.separator + "Split channels" + File.separator + "IFT marker" + File.separator;
	input_folder_BB = parent_input_folder  + File.separator + "Split channels" + File.separator + "Basal body marker" + File.separator;
	input_folder_cilia = parent_input_folder  + File.separator + "Split channels" + File.separator + "Ciliary marker" + File.separator;
	general_output_folder = data_from_other_macro[0];
	output_folder = CreateOutputDirectory(general_output_folder);
}

else {
	stand_alone = 1;
	print ("Choose input folders");
	input_folder_IFT = getDirectory("Choose input folder for the IFT channel");
	input_folder_BB = getDirectory("Choose input folder for the basal body channel");
	input_folder_cilia = getDirectory("Choose input folder for the ciliary channel");
	silent_mode = RunSilentMode();
	output_folder = CreateOutputDirectory(input_folder_IFT);
	general_output_folder = output_folder;
	print("Input folder IFT:", input_folder_IFT);
	print("Input folder Basal bodies:", input_folder_BB);
	print("Input folder Cilia:", input_folder_cilia);
	print("Output folder:", output_folder);
}


//Save the macro that was used
File.copy(File.directory + File.name, general_output_folder + "Macro_used_to_analyze_IFT_accumulations.txt");

//Gather all variables used in the analysis
automated_variables = CalculateVariables(input_folder_IFT);
user_variables = UserVariables();

PrintVariables(automated_variables, user_variables, general_output_folder);

AnalyzeImages(input_folder_IFT, input_folder_BB, input_folder_cilia, user_variables, automated_variables, output_folder, silent_mode);
SaveLog(output_folder,stand_alone);
print ("Finished analyzing IFT accumulations");


//Create output directory
function CreateOutputDirectory(input_folder){
	parent_folder = File.getParent(input_folder);
	output_folder = parent_folder + File.separator + "IFT analysis" + File.separator;
	File.makeDirectory(output_folder);
	return output_folder;
}

//Calculates the optimal values for analyzing the nuclei
function CalculateVariables(input_folder){
	images = getFileList(input_folder);
	input_path = input_folder + images[0];
	open(input_path);
	getPixelSize(unit, pixel_width, pixel_height);
	
	//Calculate the optimal sigma value for the gaussian blur based on the image resolution. The sigma value can never be smaller then one
	sigma = round(0.25/pixel_width);
	if (sigma<1){
		sigma = 1;
	}
	variables = newArray(3);
	variables[0] = sigma;

	//Calculate the particle size of the cilia based on the image resolution
	min_area_cilia = 0; //Minimal area in um^2
	particle_size_min_cilia = (min_area_cilia/(pixel_width * pixel_height));
	max_area_cilia = 41; //Maximal area in um^2
	particle_size_max_cilia = (max_area_cilia/(pixel_width * pixel_height));
	string_size_cilia = "" + particle_size_min_cilia + "-" + particle_size_max_cilia;
	variables[1] = string_size_cilia;
	
	//Calculate the particle size of the basal bodies based on the image resolution
	min_area_BB = 4; //Minimal area in um^2
	particle_size_min_BB = (min_area_BB/(pixel_width * pixel_height));
	max_area_BB = 31; //Maximal area in um^2
	particle_size_max_BB = (max_area_BB/(pixel_width * pixel_height));
	string_size_BB = "" + particle_size_min_BB + "-" + particle_size_max_BB;
	variables[2] = string_size_BB;
	
	close();
	return variables;
}

//Gather user defined variables
function UserVariables(){
	user_variables = newArray(6);
	user_variables[0] = DetermineVariableByUser("Signal threshold", "IFT", 0);
	user_variables[1] = DetermineVariableByUser("Particle size", "IFT", 0);
	user_variables[2] = DetermineVariableByUser("Signal offset", "Cilia", 0);
	user_variables[3] = DetermineVariableByUser("Signal threshold", "Cilia", 0);
	user_variables[4] = RemoveBB();
	if (user_variables[4]){
		user_variables[5] = DetermineVariableByUser("Signal threshold", "BB", 0);
	}
	else {
		user_variables[5] = 0;
	}
	return user_variables;
}

//Define user variables
function DetermineVariableByUser(value_type, channel_name, number){
	Dialog.create("Define variable");
	if (value_type == "Signal offset"){
		message = "Measure cilia:"
		+ "\n" + "Define the " + value_type + " you want to apply. You can only entre a whole number." 
		+ "\n" + "If you do not want to use this function type 0.";
	}
	if (value_type == "Signal threshold"){
		message = "Define the " + value_type + " you want to apply to the " + channel_name + " channel by"
		+ "\n" + "defining the minimal value. You can only enter a whole number <254. If you want to use" 
		+ "\n" + "automated treasholding per image type 0.";
	}	
	if (value_type == "Particle size"){
		message = "Define the " + value_type + " you want to apply the " + channel_name + " channel. You can only" 
		+ "\n" + "enter a whole number. If the image treshold is sufficient you can leave this at 0, however if there"
		+ "\n" + "are single background pixels that interfere with your analysis you can increase it to 10 or 100.";
	}
	Dialog.addMessage(message);
	Dialog.addNumber(value_type, number);
	Dialog.show();
	value_output = Dialog.getNumber(); 
	return value_output;
}

//USER defines if the macro will exclude BB regions from the measurements
function RemoveBB(){
	BB = getBoolean("Would you like to remove the Basal Body staining from the area that is measured", "Yes, only measure tip accumulations", "No, measure all accumulations");
	return BB;
}


//Define if the macro will be ran in silent mode
function RunSilentMode(){
	silent_mode = getBoolean("Do you want to run the macro in silent mode? This will not open images while running." 
	+ "\n" + "If you select 'no' you will not be able to use the computer while the macro is running. ");
	setBatchMode(silent_mode);
	return silent_mode;
}

//Saves the macro that was used and the user defined variables in separate text files for future reference
function PrintVariables(automated_variables, user_variables, output_folder){
	file = File.open(output_folder + "Variables_used_for_IFT_analysis.txt");
	print(file, "User defined variables:");
	print(file, "Signal threshold IFT channel: " + user_variables[0]+ ", 255 (min, max). 0 indicates the MaxEntropy Auto treshold was used, as indicated in the Log.");
	print(file, "Particle size IFT channel: " + user_variables[1]);
	print(file, "Signal offset cilia channel: " + user_variables[2] + ", 255, 7 (min, max, colour channel)");
	print(file, "Signal threshold cilia channel: " + user_variables[3] + ", 255 (min, max). 0 indicates the MaxEntropy Auto treshold was used, as indicated in the Log.");
	if (user_variables[4]){
		print(file, "Only accumulations in the ciliary tip region were measured. This might include some axonemal accumulations as well.");
		print(file, "Signal threshold basal body channel: " + user_variables[5] + ", 255 (min, max). 0 indicates the Li Auto treshold was used, as indicated in the Log.");
	}
	else {
		print(file, "Accumulations along the whole cilium were measured.");
	}
	print(file, "");
	print(file, "Automated variables:");
	print(file, "Gaussian blur sigma all channels: " + automated_variables[0]);
	print(file, "Particle size cilia channel: " + automated_variables[1]);
	print(file, "Particle size basal body channel: " + automated_variables[2]);
	File.close(file);
}

//Image analysis
function AnalyzeImages(input_IFT, input_BB, input_cilia, user_variables, automated_variables, output_folder, silent_mode){
	images = getFileList(input_IFT);
	print ("Analyzing IFT accumulations:");
	for (i=0; i<images.length; i++) {
		print (images[i]);
		
	//Obtain ciliary ROI
		input_path = input_cilia + images[i];
		open(input_path);
		//Set signal offset if indicated by the user
		if (user_variables[2] > 0){
			run("RGB Color");
			setMinAndMax(user_variables[2], 255, 7);
			run("8-bit");
		}
		//Set signal treshold as indicated by the user (0=automatic tresholding)
		if (user_variables[3] > 0){
			setAutoThreshold("Default dark");
			setThreshold(user_variables[3], 255);
		}
		else {
			print ("Ciliary channel treshold value:");
			run("Auto Threshold", "method=MaxEntropy white show");
		}
		//Crop image to remove all non-ciliary regions
		run("Convert to Mask");
		for (j=0; j<50; j++) {
			setOption("BlackBackground", false);
			run("Dilate");
		}
		//Obtain the ciliary ROIs from the remaining image
		sigma = "sigma=" + automated_variables[0];
		run("Gaussian Blur...", sigma);
		run("Make Binary");
		size = "size=" + automated_variables[1] + " pixel circularity=0.00-1.00 show=Overlay add";
		run("Analyze Particles...", size);
		run("Create Selection");
		run("Make Inverse");
		run("Create Mask");
		rename("Ciliary mask");
		selectWindow(images[i]);
		close();
		roiManager("reset");

	//Obtain basal bodies ROI
		if (user_variables[4]){
			input_path = input_BB + images[i];
			open(input_path);
			run("8-bit");
			//Set signal treshold as indicated by the user (0=automatic tresholding)
			if (user_variables[5] > 0){
				setAutoThreshold("Default dark");
				setThreshold(user_variables[5], 255);
			}
			else {
				print ("Basal body channel treshold value:");
				run("Auto Threshold", "method=Li white show");
			}
			//Obtain the basal body ROIs from the remaining image
			run("Convert to Mask");
			run("Dilate");
			run("Gaussian Blur...", sigma);
			run("Make Binary");
			run("Fill Holes");
			size = "size=" + automated_variables[2] + " pixel circularity=0.00-1.00 show=Overlay add";
			run("Analyze Particles...", size);
			run("Create Selection");
			run("Create Mask");
			rename("Basal Bodies mask");
			selectWindow(images[i]);
			close();
			roiManager("reset");
		}
		
	//Analyze IFT image
		print ("Analyzing IFT accumulations");
		input_path = input_IFT + images[i];
		open(input_path);
		run("8-bit");
		//Remove basal body regions from the images with only ciliary regions
		if (user_variables[4]){
			masks = newArray("Basal Bodies mask", "Ciliary mask");
		}
		else {
			masks = newArray("Ciliary mask");
		}
		setBackgroundColor(0, 0, 0);
		for (k=0; k<masks.length; k++) {
			selectWindow(masks[k]);
			run("Create Selection");
			selectWindow(images[i]);
			run("Restore Selection");
			run("Clear", "slice");
			run("Select None");
		}
		output_path = output_folder + "crop_" + images[i];
		save(output_path);
		close("Basal Bodies mask");
		close("Ciliary mask");
		//Set signal treshold as indicated by the user (0=automatic tresholding)
		if (user_variables[0] > 0){
			setAutoThreshold("Default dark");
			setThreshold(user_variables[0], 255);
		}
		else {
			run("Auto Threshold", "method=MaxEntropy white show");
		}
		//Measure IFT accumulations in the remaining ciliary regions without basal bodies or background
		run("Convert to Mask");
		run("Gaussian Blur...", sigma);
		run("Make Binary");
		run("Fill Holes");
		
		size = "size=" + user_variables[1] + " pixel circularity=0.00-1.00 show=[Overlay Outlines] display summarize record add";
		run("Analyze Particles...", size);
		
	//Creating output		
		output_path = output_folder + images[i] +"RoiSet.zip";
		roiManager("Save", output_path);
		run("Labels...", "color=white font=36 show draw");
		run("Flatten");
		output_path = output_folder + "Labeled_" + images[i];
		saveAs("Tiff", output_path);
		close();
		close();

		//Create a merged image with accumulation overlay
		open(input_path);
		rename("IFT");
		input_path = input_cilia + images[i];
		open(input_path);
		rename("Cilia");
		input_path = input_BB + images[i];
		open(input_path);
		rename("BB");
		run("Merge Channels...", "c1=IFT c2=Cilia c5=BB create");
		run("Stack to RGB");
		run("From ROI Manager");
		roiManager("show all without labels");
		run("Flatten");
		run("Images to Stack", "name=Stack title=[] use");
		run("Remove Overlay");
		output_path = output_folder + "merge_" + images[i];
		saveAs("Tiff", output_path);
		roiManager("reset");
		close();

		//Save results
		selectWindow("Results"); 
	   	name = images[i] + ".xls"; 
	   	saveAs("Measurements", output_folder + name);
		selectWindow("Results"); 
	   	run("Close"); 
	 	print("Silent mode: ", silent_mode);
	 	if (silent_mode == false){
			selectWindow("ROI Manager"); 
			run("Close");
			run("Close");
		}
	}
}

//Saves the summary and the log file of the run
function SaveLog(output_folder,run_type){
	selectWindow("Summary");
	saveAs("Results", output_folder +"Summary.csv");
	run("Close");
	if (run_type == 1){
		selectWindow("Log");
		saveAs("txt", output_folder +"Log.txt");
	}
}