// Written by Cenna van Manen.

// This macro will measure the cilia per image, based on the user defined variables.

// As input, the macro assumes all images in the chosen folder are single layered tif files. If required, these files 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 cilia.
// Furthermore, three text files will be generated, one containing the macro, one containing the log of the run, and the other containing the used variables.



print ("Start measuring cilia length");
requires("1.52n")

//Determine if the macro is called upon by another macro, otherwise the variables 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 = 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 folder");
	input_folder = getDirectory("Choose input folder");
	output_folder = CreateOutputDirectory(input_folder);
	general_output_folder = output_folder;
	print("Input folder:", input_folder);
	print("Output folder:", output_folder);
	silent_mode = RunSilentMode();
}

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

//Gather all variables used in the analysis
minmax_threshold = DetermineVariableByUser("Signal offset", 0);
greyvalue_threshold = DetermineVariableByUser("Signal threshold", 0);
gaussian_blur = CalculateSigma(input_folder);

PrintVariables(minmax_threshold, greyvalue_threshold, gaussian_blur, general_output_folder);

AnalyzeImages(input_folder, minmax_threshold, greyvalue_threshold, gaussian_blur, output_folder, silent_mode);
SaveLog(output_folder,stand_alone);

print ("Finished measuring cilia");


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

//Define user variables
function DetermineVariableByUser(value_type, number){
	Dialog.create("Define variables to measure the cilia");
	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 = "Measure cilia:"
		+ "\n" + "Define the " + value_type + " you want to apply by defining the minimal value. You can only enter a" 
		+ "\n" + " whole number <254. If you want to use automated treasholding per image type 0.";
	}	
	Dialog.addMessage(message);
	Dialog.addNumber(value_type, number);
	Dialog.show();
	value_output = Dialog.getNumber(); 
	return value_output;
}

function CalculateSigma(input_folder){
	images = getFileList(input_folder);
	input_path = input_folder + images[0];
	open(input_path);
	getPixelSize(unit, pixel_width, pixel_height);
	sigma = round(0.25/pixel_width);
	if (sigma<1){
		sigma = 1;
	}
	if (sigma>5){
		sigma = 5;
	}
	close();
	print ("Optimal gaussian blur sigma: ", sigma);
	return sigma
}

//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;
}

//Define dilate cycles
function DetermineDilateRepeatsAndCutOff(){
	getPixelSize(unit, pixel_width, pixel_height);
	gap_width = 0.5; //Set as a constant of 0.5 um
	dilate_correction = (0.5*gap_width*(1/pixel_width));
	dilate_repeat = floor(dilate_correction)+1; //+1 to make sure the gap is covert
	cut_off = 3 * pixel_width; //The minimal size that can be measured. Based on the skeleton width of 1 to 2 pixels, all correctly measured values have to be smaller then 3 pixels.
	return_values = newArray(dilate_repeat, cut_off);
	return return_values;
}

//Saves the user defined variables in separate text files for future reference
function PrintVariables(minmax_threshold, greyvalue_threshold, gaussian_blur, output_folder){
	file = File.open(output_folder + "Variables_used_to_measure_cilia.txt");
	print(file, "User defined variables:");
	print(file, "Signal offset: " + minmax_threshold + ", 255, 7 (min, max, colour channel)");
	print(file, "Signal threshold: " + greyvalue_threshold + ", 255 (min, max). 0 indicates the MaxEntropy Auto treshold was used, as indicated in the Log.");
	print(file, "");
	print(file, "Automated variables:");
	print(file, "Gaussian blur sigma: " + gaussian_blur);
	File.close(file);
}

//Image analysis
function AnalyzeImages(input_folder, minmax_threshold, greyvalue_threshold, gaussian_blur, output_folder, silent_mode){
	images = getFileList(input_folder);
	for (i=0; i<images.length; i++) {
		input_path = input_folder + images[i];
		print ("Measuring cilia in ", images[i]);
		open(input_path);
		if (i==0){
			determined_values = DetermineDilateRepeatsAndCutOff();
			dilate_repeat = determined_values[0];
			cut_off = determined_values[1];
		}
		if (minmax_threshold > 0){
			run("RGB Color");
			setMinAndMax(minmax_threshold, 255, 7);
			run("8-bit");
		}
		if (greyvalue_threshold > 0){
			setAutoThreshold("Default dark");
			setThreshold(greyvalue_threshold, 255);
		}
		else {
			run("Auto Threshold", "method=MaxEntropy white show");
		}
		run("Convert to Mask");
		for (j=0; j<dilate_repeat; j++) {
			setOption("BlackBackground", false);
			run("Dilate");
		}
		sigma = "sigma=" + gaussian_blur;
		run("Gaussian Blur...", sigma);
		run("Make Binary");
		getPixelSize(unit, pixel_width, pixel_height);
		area = ((2*(cut_off+(dilate_repeat*pixel_width)))^2);
		size = "size=" + area + "-Infinity circularity=0.00-1.00 show=[Overlay Outlines] display summarize record add";
		run("Analyze Particles...", size);
		name = replace(images[i], ".tif", ".csv");
		selectWindow("Results");
		saveAs("Measurements", output_folder + "counts_" + name); 
		run("Close");
		//Remove background
		setForegroundColor(0, 0, 0);
		setBackgroundColor(255,255,255);
		
		nROIs = roiManager("count");
		if (nROIs > 1){
			roiManager("Combine");
		}
		if (nROIs == 1) {
			roiManager("select", 0);
		}
		run("Make Inverse");
		run("Clear", "slice");
		run("Select None");
		
		roiManager("Set Color", "cyan");
		roiManager("Set Line Width", 1);
		run("Skeletonize (2D/3D)");

		//Creating output
		saveAs("Tiff", output_folder + "skeleton_" + images[i]);
		size = "size=" + 3*cut_off + "-Infinity pixel show=Overlay add";
		run("Analyze Particles...", size);
		run("Analyze Skeleton (2D/3D)", "prune=none");
		//saveAs("Tiff", output_folder + "length_" + images[i]);
	   	name = replace(images[i], ".tif", ".csv");
		selectWindow("Results");
		saveAs("Measurements", output_folder + "length_" + name); 
		run("Close");
		close();
		close();
		
		//Create image with skeleton overlay
		open(input_path);
		run("From ROI Manager");
		roiManager("show all without labels");
		run("Flatten");
		run("Images to Stack", "name=Stack title=[] use");
		run("Remove Overlay");
		saveAs("Tiff", output_folder + "merge_" + images[i]);
		roiManager("reset");
		close();
	}
	if (silent_mode == false){
		selectWindow("ROI Manager"); 
		run("Close");
	}		
	print ("Cut off value: ", cut_off);
	SummarizeData(input_folder, output_folder, cut_off);
}

//Create the summary file
function SummarizeData(input_folder, output_folder, cut_off){
	names = getFileList(input_folder);
	file = File.open(output_folder + "Summary_length.csv");
	print(file, "Label, Uncorrected measurement, Corrected cilium length");
	for (i=0; i<names.length; i++) {
		file_name = replace(names[i], ".tif", ".csv");
		file_name = "length_" + file_name;
		input_path = output_folder + file_name;
		open(input_path);
		selectWindow(file_name);
		IJ.renameResults(file_name,"Results");
		measurements = nResults;
		for (j = 0; j < measurements; j++) {
			branches = getResultString("# Branches",j);
			endpoint_voxels = getResultString("# End-point voxels",j);
			measured_length = getResult("Maximum Branch Length",j);
			//Values that are smaller then the cut off value or do not have exactly one branch with two end points are excluded from the summary
			if (branches==1){
				if (endpoint_voxels==2){
					if (cut_off<measured_length){
						//Measurements were corrected for the error introduced by the skeletonization algeorithm.
						//The corrected length of the measured values is based on the formula y = ax+b, for which y is the measured length,
						//x is the corrected length and a & b were determined using a large test set, see manual.
						a = 0.964857;
						b = -0.572969;
						
						//The corrected length can then be determine as x = (y-b)/a.
						y_b = (measured_length-b);
						corrected_length = ((y_b)/a);
						
						//Summarize the resultes
						file_name = replace(names[i], ".tif", "");
						correct_measurement = file_name + "," + measured_length + "," + corrected_length;
						print (file, correct_measurement);
					}
				}
			}
		}
	}
	close("Results");
}

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