/*
 * Roughness measurement macro
 * v1.17.10
 * Developed by M Foley, Sydney Microscopy & Microanalysis
 *
 * function processFolder loops through to find files for analysis & then produces the output
 * function processImage actually runs the analysis on individual files
 * function updateLogBar produces a scrolling 'progress' bar for the log window
 * 
 * Assumption: slices start at the bottom of the sample, highest slice = furtherest from the base
 */

startTime = getTime();																						// Timekeeping

setBatchMode(true);

var fileNames = newArray();																					// Arrays to hold results as Global Variables
var filezRq = newArray();
var filezRa = newArray();
var filezRv = newArray();
var filezRp = newArray();
var filezRt = newArray();
var filezRsk = newArray();
var filezRku = newArray();
var fileStdDevRA = newArray();

dir1 = getDirectory("Choose Source Directory ");															// Location of datasets
extension = ".am";																							// Format - .am or .tiff

if (getBoolean("Process folder "+dir1, "Yes please", "No thanks") ==1){  processFolder(dir1);}				// Last chance to check things



  function processFolder(dir1) {
  	print("\\Clear");
  	print("Processing folder "+dir1);																		// Line 0;
	print("\n");																							// Reserved for *****
	print("\n");																							// Update of what we're doing
	print("------------------------------------------");													// Spacer
     list = getFileList(dir1);
     for (i=0; i<list.length; i++) {
          if (endsWith(list[i], "/"))																		// if directory, redirect
              processFolder(dir1+list[i]);
          else if (endsWith(list[i], extension)&&indexOf(list[i], "Composite")<0)							// Ignore previously processed files
             processImage(dir1, list[i], i, list.length);
      }
///////////////////// processFolder finished loop, cleaning up /////////////////////
Files 			= fileNames;																				// Rename arrays for Results window
Sq 				= filezRq;
Min_Sv 			= filezRv;
Sa 				= filezRa;
Sa_StdDev		= fileStdDevRA;
Max_Sp 			= filezRp;
Max_St 			= filezRt;
Mean_Ssk 		= filezRsk;
Mean_Sku 		= filezRku;

totalTime 		= getTime()-startTime;																		// Time in ms
totalMinutes	= (totalTime-totalTime%60000)/60000;
totalSeconds	= round((totalTime-totalMinutes*60000)/1000);
if(totalSeconds<10){totalSeconds="0"+totalSeconds;}															// Formatting fix

Array.show("Results", Files, Sa, Sa_StdDev, Sq, Min_Sv, Max_Sp, Max_St, Mean_Ssk, Mean_Sku);
	saveAs("Results", dir1+"Analysis Results.csv");
   	print("\\Update1: ");
 	print("\\Update2: ");
	print("------------------------------------------");													// Spacer
  	print("All files in "+dir1+" have been processed");
  	print("Total time: "+totalMinutes+":"+totalSeconds);													// Time for run
  }
////////////////////////////////////////////////////////////////////////////////////

function processImage(dir1, name, i_curr, i_max) {
	print("\\Update2:Opening file "+name);
	if(extension == ".am")
		{run("Amira...", "amirafile=["+dir1+name+"]"); }													// If .am, use importer, else open as normal
	else
		{open(dir1+name);}
	
	run("Select None");																						// Ensure that no selections exist
	
	run("8-bit");																							// Ensure 8-bit compatability
	
	setAutoThreshold("Default dark");																		// Convert to binary
	setThreshold(1, 255);
	run("Convert to Mask", "method=Default background=Default black");
	
	getVoxelSize(vox_width, vox_height, vox_depth, vox_unit);												// Store the image dimension details
	
	originalName = getTitle();																				// Get image name
	originalName = substring(originalName,0, indexOf(originalName,extension));								// Trim image name of extension
	rename(originalName);
	
	print("\\Update2:Reslicing "+originalName);																// Update log window
	run("Reslice [/]...", "output=1 start=Bottom avoid");													// Reorient dataset
	rename("temporary");																					// Close the original dataset
	selectWindow(originalName);
	close();
	selectWindow("temporary");
	rename(originalName);
		
	run("Select None");																						// Ensure that no selections exist
	run("Duplicate...", "title=[ROI_Poly_"+originalName+"] duplicate");										// Duplicate window to create the Poly-filled version
	selectWindow(originalName);																				// Return to original window
	
	print("\\Update2:Finding edges of "+originalName);														// Update log window
	run("Find Edges", "stack");																				// Find edges
	setOption("BlackBackground", true);
	setAutoThreshold("Default dark");																		
	run("Convert to Mask", "method=Default background=Default black");										// Re-thresh to clear any 'fuzziness'
	print("\\Update2:Skeletonization of "+originalName+" (may take some time)");							// Update log window
	run("Skeletonize", "stack");																			// Skeletonize
	
	getDimensions(width, height, channels, slices, frames);													// Get dimensions for array storage. Channels & frames irrelevant
	
	xPos = Array.getSequence(width);																		// Prepopulate the xPosition array with values
	yPos = newArray(width);																					// All other arrays are 0
	zRa = newArray(slices);
	zRq = newArray(slices);
	zRv = newArray(slices);
	zRp = newArray(slices);
	zRt = newArray(slices);
	zRsk = newArray(slices);
	zRku = newArray(slices);
	
	print("\\Update2:Processing slices of "+originalName);													// Update log window

	for (z_Image=1; z_Image<=slices; z_Image++){															// Begin loop through slices
		print("\\Update1:"+updateLogBar(z_Image));															// Update log window
		
		///////////////// Status bar information /////////////////
		
			image_complete = round(z_Image*100/slices);
			analysis_complete = round(i_curr*100/i_max);
			image_complete = image_complete-image_complete%5;
			analysis_complete = analysis_complete-analysis_complete%5;
			showStatus("Approximately "+image_complete+" / "+analysis_complete+"% processed (file / overall)");
			showProgress(z_Image/slices);
		
		///////////////// Status bar information /////////////////
		
		selectWindow(originalName);																			// Ensure we're on the right image
		Overlay.remove;																						// Remove any overlay
		setSlice(z_Image);																					// Set slice
		
		for (x=0; x<width; x++){																			// Loop through each y value (from 'top' to bottom)
			for (y=height; y>0; y--){																		// to find where the first non-zero pixel exists
				val = getPixel(x,y);																		// when found, dump to the xPos & yPos arrays
				if (val==255){
					xPos[x] = x;
					yPos[x] = y;
					y = 0;
				}
				if(x==width-1)
				{yPos[x] = yPos[x-1];}																		// Fix for edge - skeletonisation erodes the final pixel
				
			}
		}
		
		yPos_Rq = 0;																						// Reset the Rn value arrays for this slice
		yPos_Ra = 0;
		yPos_Rv = height;
		yPos_Rp = 0;
		yPos_Rt = 0;
		yPos_Rsk = 0;
		yPos_Rku = 0;
		Array.getStatistics(yPos, min, max, mean, std);														// Find the mean yPos value
	
		
		for (i=0; i<lengthOf(yPos); i++)																	// loop through each column
		{
			yi = yPos[i] - mean;																			// yi used repeatedly - cleaner to calculate first

			yPos_Ra = yPos_Ra + abs(yi)*vox_height;
			yPos_Rq = yPos_Rq + pow(yi*vox_height,2);
			// Rsk
			yPos_Rsk = yPos_Rsk + pow((yi)*vox_height, 3);
			// Rku
			yPos_Rku = yPos_Rku + pow((yi)*vox_height, 4);
		}
		
		zRa[z_Image-1] = yPos_Ra/lengthOf(yPos);															// Append values against _this_ slice position
		zRq[z_Image-1] = sqrt(yPos_Rq/lengthOf(yPos));
		zRv[z_Image-1] = (min-mean)*vox_height;
		zRp[z_Image-1] = (max-mean)*vox_height;
		zRt[z_Image-1] = zRp[z_Image-1]-zRv[z_Image-1];
		zRsk[z_Image-1] = yPos_Rsk/(lengthOf(yPos)*pow(zRq[z_Image-1],3));
		zRku[z_Image-1] = yPos_Rku/(lengthOf(yPos)*pow(zRq[z_Image-1],4));
		
		selectWindow("ROI_Poly_"+originalName);																// Prepare to generate the filled polygon
		setSlice(z_Image);																					// Ensure we're on the right slice
		
		xPolygon = Array.concat(0, xPos, width, width);														// Additional points required to complete polygon
		yPolygon = Array.concat(0, yPos, yPos[i-1], 0);
		makeSelection("polygon",xPolygon, yPolygon);														// Make & fill polygon with white pixels
		setColor("white");
		run("Fill", "slice");
	}																										// Slice loop complete
	
	print("\\Update2:Calculating results");																	// Update log window
																											// Calculate overall results for sample
	Array.getStatistics(zRa, Ramin,Ramax,Ramean,Rastd);	
	Array.getStatistics(zRq, Rqmin,Rqmax,Rqmean,Rqstd);
	Array.getStatistics(zRv, Rvmin,Rvmax,Rvmean,Rvstd);
	Array.getStatistics(zRp, Rpmin,Rpmax,Rpmean,Rpstd);
	Rtmax 		= Rpmax - Rvmin;																			// Array.getStatistics(zRt, Rtmin,Rtmax,Rtmean,Rtstd) not used
																											// Rt should be between max & min over WHOLE surface
	Array.getStatistics(zRsk, Rskmin,Rskmax,Rskmean,Rskstd);
	Array.getStatistics(zRku, Rkumin,Rkumax,Rkumean,Rkustd);
	
	fileNames 	= Array.concat(fileNames, originalName);													// Append results to global arrays
	filezRa 	= Array.concat(filezRa, Ramean);
	filezRq 	= Array.concat(filezRq, Rqmean);
	filezRv 	= Array.concat(filezRv, Rvmin);
	filezRp 	= Array.concat(filezRp, Rpmax);
	filezRt 	= Array.concat(filezRt, Rtmax);
	filezRsk 	= Array.concat(filezRsk, Rskmean);
	filezRku 	= Array.concat(filezRku, Rkumean);
	//
	fileStdDevRA= Array.concat(fileStdDevRA, Rastd);
	//
	
	print("\\Update2:Preparing final image of "+originalName);												// Update log window
	run("Merge Channels...", "c1=[ROI_Poly_"+originalName+"] c4=["+originalName+"] create keep ignore");	// Merge composite image of original outline & filled regions
	rename(originalName+"-Composite");																		// Rename accordingly
	print("\\Update2:Converting image of "+originalName+" for saving");										// Update log window
	run("Stack to RGB", "slices");																			// Convert from composite to RGB
	run("8-bit");																							// Convert from RGB to 8-bit greyscale
	print("\\Update2:Saving image of "+originalName);														// Update log window
	
	if (extension ==".am")																					// Save in original extension
		{run("AmiraMesh ...", "amirafile=["+dir1+originalName+"-Composite.am]");}
	else
		{saveAs(".tif", dir1+originalName+"-Composite");}

	totalTime 		= getTime()-startTime;																	// Time in ms
	totalMinutes	= (totalTime-totalTime%60000)/60000;
	totalSeconds	= round((totalTime-totalMinutes*60000)/1000);
	if(totalSeconds<10){totalSeconds="0"+totalSeconds;}														// Formatting fix

	
	print(totalMinutes+":"+totalSeconds+"     Completed "+originalName);									// Print to log
	run("Close All");																						// Clean up

///////////////////// Updating results window /////////////////////
Files = fileNames;																							// Rename arrays for Results window
Sq = filezRq;																								// Updates on each loop
Min_Sv = filezRv;
Sa = filezRa;
Sa_StdDev		= fileStdDevRA;
Max_Sp = filezRp;
Max_St = filezRt;
Mean_Ssk = filezRsk;
Mean_Sku = filezRku;
	
Array.show("Results", Files, Sa, Sa_StdDev, Sq, Min_Sv, Max_Sp, Max_St, Mean_Ssk, Mean_Sku);
///////////////////// ---------------------- /////////////////////

}


function updateLogBar(i) {
	a_= fromCharCode(9633,9632,9632,9632,9632);
	b_= fromCharCode(9632,9633,9632,9632,9632);
	c_ = fromCharCode(9632,9632,9633,9632,9632);
	d_ = fromCharCode(9632,9632,9632,9633,9632);
	e_ = fromCharCode(9632,9632,9632,9632,9633);
	
	bars = newArray(a_,b_,c_,d_,e_,d_,c_,b_);
	n = bars.length;

	return bars[i/n%n];
}