/* 
 * THRESHOLD ANALYSER
 * Script to assist in the selection of a threshold value for an image stack
 * Matt Foley, SMM, USYD
 * Version 0.1b (Beta)
 * 2018-02-28--2018-03-01
 */

/*
 * Structure of script
 * 
 * 	Check for suitable input
 * 		Must be an image
 * 		Currently must be a stack
 * 		
 * 	Determine analysis parameters
 * 		Image Name
 * 		BitDepth
 * 		Accuracy of analysis
 * 		
 * 	Initialise variables & arrays
 * 		Arrays
 * 		Threshhold variables
 * 		Status variables
 * 		Timing variables
 * 		Object analysis variables
 * 		 		
 * 	Main body of code
 * 	
 * 	Functions
 */

//	Check for suitable input and exit if failure
if ((getInfo("window.type") != "Image")|nSlices<2)
	{exit("Inappropriate input (Non-image or stack selected)");}

//	Deselect any selection to prevent ROI errors
run("Select None");

//	Determine analysis parameters
//		Determine image name and image bitdepth
var_original 			= getTitle();
var_bitDepth 			= bitDepth();

//		Accuracy of analysis
var_sensitivity 		= getNumber("How accurate do you want to be? (Between 3 and "+var_bitDepth+" where smaller = less accurate)", 5);

//	Enable batch mode for speed
setBatchMode(true);

//	Reposition and resize image window to fit future windows
setLocation(375,436, 553, 620);
			
//	Initialise variables & arrays
//		Arrays
array_Thresh 			= newArray();
array_Measure 			= newArray();
array_ObjectCount 		= newArray();
array_ObjectTotalArea 	= newArray();
array_ObjectPctArea 	= newArray();

//		Threshold variables
var_ThreshMin 			= 0;
var_ThreshMax 			= pow(2,var_bitDepth) -1;
var_ThreshStep 			= (var_ThreshMax+1)/(pow(2,var_sensitivity));	

//		Status variables
status_TitleShort 		= "Progress";
status_Title 			= "["+status_TitleShort+"]";
status_ProgressText		= "\\Update:";

//		Timing variables
status_StartTime 		= getTime();
status_PreviousTime 	= 0;
status_ComplCurrTime	= 0;
status_TimePerStep		= 0;
status_CurrentTime		= 0;

//		Metrics for no reason other than they're there
//		Values chosen as they will be max'd/min'd appropriately
status_MaxRunTime 		= 0;
status_MinRunTime 		= 1e13;

//		Object analysis variables
var_ObjectCount 		= 0;
var_ObjectTotalArea 	= 0;
var_ObjectPctArea 		= 0;

/***************************************************************/

//	Set up progress window in upper left corner
//	If it is open, clear it
//	If it isn't open, create it

if(isOpen(status_TitleShort)){print(status_Title, status_ProgressText);}
else
{run("Text Window...", "name="+ status_Title +" width=50 height=4");
setLocation(0, 0);
}

//	Close previous result windows if open
fn_IfOpenClose("Overall Results");
fn_IfOpenClose("Plotted Results");

//	Loop through threshhold values
for (i=var_ThreshMin; i<=var_ThreshMax+1; i+=var_ThreshStep){

	//	Status updates - time
			status_CurrentTime 						= (getTime()-status_StartTime)/1000;
			status_TimePerStep 						= status_CurrentTime-status_PreviousTime;
			status_CompleteTime 					= round((var_ThreshMax - i)/var_ThreshStep * status_TimePerStep);
			status_PreviousTime						= status_CurrentTime;
			status_ComplCurrTime					= round(status_CompleteTime+status_CurrentTime);
	//	Update max & min run time estimates		
		if(status_ComplCurrTime> status_MaxRunTime)
			status_MaxRunTime=status_ComplCurrTime;
		if(status_ComplCurrTime< status_MinRunTime&&i>0)		//	&&i>0 deals with a Min value existing before any measurements made
			status_MinRunTime=status_ComplCurrTime;
	//	Update the progress text
			status_ProgressText = "Processing threshold "+i+
				" of "+var_ThreshMax+1+" ("+i/var_ThreshStep+1+
				"/"+((var_ThreshMax-var_ThreshMin+1)/var_ThreshStep)+1+")";
	//	Update Progress window				
		print(status_Title, "\\Update:"+status_ProgressText+
			"\nRuntime so far: "+status_CurrentTime+
			" s\nEstimated time remaining: "+status_CompleteTime+
			" s ("+status_ComplCurrTime+" s total)");

	//	Select original image
			selectWindow(var_original);
	//	Duplicate and rename to threshold value
			 run("Duplicate...", "title="+i+" duplicate");
			 selectWindow(i);
	//	Set threshhold range and convert to mask
			setThreshold(0, i);
			run("Convert to Mask", "method=Default background=Dark black");	 
	//	Get stack statistics and drop into arrays
			Stack.getStatistics(voxelCount, mean, min, max, std);
			array_Thresh 						= Array.concat(array_Thresh, i);
			array_Measure 						= Array.concat(array_Measure, mean/255);
	//	Analyse the threshholded image for object details	
		run("Analyze Particles...", "exclude clear summarize");
	//	Reset object values
			var_ObjectCount 					= 0;
			var_ObjectTotalArea 				= 0;
			var_ObjectPctArea 					= 0;
	//	Extract details from summary window
			selectWindow("Summary of "+i); 
			lines 								= split(getInfo(), "\n"); 
			headings 							= split(lines[0], "\t"); 
			values 								= split(lines[1], "\t"); 
			for (j=0; j<headings.length; j++){ 
				if (headings[j] == "Count")
					{var_ObjectCount 			= var_ObjectCount+ values[j];}
				if (headings[j] == "Total Area")
					{var_ObjectTotalArea 		= var_ObjectTotalArea+ values[j];}
				if (headings[j] == "%Area")
					{var_ObjectPctArea 			= var_ObjectPctArea+ values[j];}
			}
	//	Close summary window
			fn_IfOpenClose("Summary of "+i);
	//	Return to threshholded image window and close it
			fn_IfOpenClose(i);

	//	Drop results into arrays
			array_ObjectCount 					= Array.concat(array_ObjectCount, var_ObjectCount);
			array_ObjectTotalArea 				= Array.concat(array_ObjectTotalArea, var_ObjectTotalArea);


/* Loop ends */		
}
//	Disable batch mode
setBatchMode(false);
wait(100);

//	Update Progress window
print(status_Title, "\\Update:All thresholds processed\nTotal runtime: "+status_CurrentTime+
	" seconds\nMin/Max time estimates: "+status_MinRunTime+
	"/"+status_MaxRunTime+" seconds");

//	Output arrays
//		Rename arrays for column titles
Threshhold 		= fn_ArrayRename(array_Thresh);
Measured		= fn_ArrayRename(array_Measure);
ObjectCount		= fn_ArrayRename(array_ObjectCount);
ObjectArea		= fn_ArrayRename(array_ObjectTotalArea);
//	Show arrays and position
Array.show("Overall Results", Threshhold, Measured, ObjectCount, ObjectArea);
selectWindow("Overall Results");
setLocation(0, 112);

//	Scaling of results for concurrent plotting
	Array.getStatistics(array_ObjectCount, min_aOC, max_aOC); 
	Array.getStatistics(array_ObjectTotalArea, min_aOTA, max_aOTA);
	Array.getStatistics(array_Measure, min_aM, max_aM);
	
	var_ScaleFactor = fn_MaxValue(max_aOC, fn_MaxValue(max_aOTA, max_aM));
	
	array_ObjectCountScaled 			= fn_ArrayMultiply(array_ObjectCount, var_ScaleFactor/max_aOC);
	array_MeasureScaled 				= fn_ArrayMultiply(array_Measure, var_ScaleFactor/max_aM);
	array_ObjectTotalAreaScaled 		= fn_ArrayMultiply(array_ObjectTotalArea, var_ScaleFactor/max_aOTA);
// Scaling end

//	Set position for Plot window
call("ij.gui.ImageWindow.setNextLocation", 375,0);

//	Create plot
Plot.create("Plotted Results", "Threshold", "");
Plot.setColor("blue");
Plot.add("line", array_Thresh, array_MeasureScaled);
Plot.setColor("red");
Plot.add("line", array_Thresh, array_ObjectCountScaled);
Plot.setColor("#10b900");
Plot.add("line", array_Thresh, array_ObjectTotalAreaScaled);
Plot.setLimitsToFit();
Plot.setFormatFlags("11001100111101");
Plot.setLegend("Measured\tObject Count\tObject Area");

//	Return to original window
selectWindow(var_original);


/*
 * Function block
 */
 
function fn_IfOpenClose(window){		// Closes named window if open
	//	Checks if window variable is an integer or a string and converts as required
	if(isNaN(parseInt(window)))
		window = window;
		else
		window = d2s(window, 0);
	if(isOpen(window)){selectWindow(window); run("Close");}
}

function fn_ArrayRename(array){
	array_out = newArray();
	for (i=0; i<lengthOf(array); i++){
		array_out = Array.concat(array_out, array[i]);
	}
	return array_out;
}

function fn_ArrayMultiply(array, var_mult){
	array_out = newArray();
	for (i=0; i<lengthOf(array); i++){
		array_out = Array.concat(array_out, array[i]*var_mult);
	}
	return array_out;
}

function fn_MaxValue(val1, val2){
	if (val1>=val2)
		return val1;
	else 
		return val2;
}