/*  Enigma+ Version 19.0404
 *  Matt Foley @ SMM (University of Sydney)
 *
 *  This macro is two parts
 *
 *  CraterExtract - generates a collection of craters from a thresholded stack and a CSV list
 *
 *  Enigma+ - analyses the crater stack
 *
	
 *	19.0404	update
 *	Line #165: 	
 *			run("Remove Outliers...", "radius=1 threshold=0 which=Dark"); // remove the thin 1 pixel border introduced by the hull
 *		changed to 	
 *			run("Remove Outliers...", "radius=1 threshold=0 which=Bright"); // remove the thin 1 pixel border introduced by the hull
 *		owing to changes in default settings for FIJI
 *	
 */

macro "CraterExtract & Enigma...[j]" {

res= 17.66; // resolution of scan in um


  if (isOpen("ROI Manager")) {
     selectWindow("ROI Manager");
     run("Close");
  }
  setBatchMode(true);

// Crater extract
// Check if dataset is appropriate
if(nSlices<2) exit("Error - no stack selected");				// Checks to ensure a stack is being analysed
getHistogram(values, counts, 256);								// Generates a list of all pixel values from 0-255
total=0;
for (i=1; i<255; i++) total+=counts[i];							// Checks if there are any pixels that are no black (0) or white (255)
if(total>0) exit("8-bit binary image (0 and 255) required.");	//Exits if there are "grey" pixels or if the image is not set correctly

origWindow = getTitle();										// Store original filename
run("Select None");

// Loads tab delimited file (TDF) containing positions

/*
 * Note: The format of the TDF MUST MATCH THE FOLLOWING EXAMPLE
 *
 * Crater#	Position	SliceStart	SliceEnd	X1	Y1	X2	Y2
 */

pathfile=File.openDialog("Choose the crater positions file:");	// Opens the tab delimited file
filestring=File.openAsString(pathfile);

rows=split(filestring, "\n");									// Splits the lines
craterID = newArray(rows.length);
cratLoc = newArray(rows.length);
SliceStart = newArray(rows.length);
SliceEnd = newArray(rows.length);
X1 = newArray(rows.length);
Y1 = newArray(rows.length);
X2 = newArray(rows.length);
Y2 = newArray(rows.length);


for(i=0; i<rows.length; i++){
	columns=split(rows[i], "\t");
	craterID[i] = columns[0];
	cratLoc[i] = columns[1];
	SliceStart[i] = columns[2];
	SliceEnd[i] = columns[3];
	X1[i] = columns[4];
	Y1[i] = columns[5];
	X2[i] = columns[6];
	Y2[i] = columns[7];
	}


// Safety checks
getDimensions(width,height,channels,slices,frames);
Array.getStatistics(X1, X1min, X1max);
Array.getStatistics(Y1, Y1min, Y1max);
Array.getStatistics(X2, X2min, X2max);
Array.getStatistics(Y2, Y2min, Y2max);
Array.getStatistics(SliceStart, SliceStartmin, SliceStartmax);
Array.getStatistics(SliceEnd, SliceEndmin, SliceEndmax);

if ((X1min <0) || (Y1min <0) || (X1max > width) || (Y1max > height)){exit("Check your X1/Y1 values");};
if ((X2min <0) || (Y2min <0) || (X2max > width) || (Y2max > height)){exit("Check your X2/Y2 values");};
if ((SliceStartmin < 0) || (SliceStartmax > slices) || (SliceEndmin < 0) || (SliceEndmax > slices)){exit("Check your Slice values");}


// generate results table
// Builds results table
tablename1 = "Crater results";
f= "["+tablename1+"]";
if(isOpen(tablename1))
	null=null;
	else{
run("Table...", "name="+f+" width=600 height=250");
print(f, "\\Headings:File\tCrater Volume (pixels)\tCrater Volume (um^3)");
	};


// Prepare to save craters and analyse individually
saveDir = getDirectory("Please select the saving directory");

for (i=1; i< rows.length; i++){
	extractCrater(i);
	enigma(i);
	saveCrater(i);
};


// return to original
selectWindow(origWindow);
run("Select None");
showMessage("Extract & Enigma of \n \n"+origWindow+"\n \nis complete.");

setBatchMode(false);





/*
 *
 *  FUNCTIONS
 *
 */

function extractCrater(n){
	selectWindow(origWindow);
	makeRectangle(X1[n], Y1[n], X2[n]-X1[n], Y2[n]-Y1[n]);
		run("Duplicate...", "title="+cratLoc[n]+"_"+craterID[n]+"_ duplicate range="+SliceStart[n]+"-"+SliceEnd[n]);
	saveExtract(n);
};





function saveExtract(n){
	folder_suff = "Extract_"+cratLoc[n]+"_"+craterID[n];
	saveDir_long = saveDir+folder_suff;
	File.makeDirectory(saveDir_long);
	fs = File.separator;
	run("Image Sequence... ", "format=TIFF save=["+saveDir_long+"]");
};

function saveCrater(n){
	folder_suff = "Crater_"+cratLoc[n]+"_"+craterID[n];
	saveDir_long = saveDir+folder_suff;
	File.makeDirectory(saveDir_long);
	fs = File.separator;
	run("Image Sequence... ", "format=TIFF save=["+saveDir_long+"]");
};


function enigma(n){

	selectWindow(cratLoc[n]+"_"+craterID[n]+"_");
run("Analyze Particles...", "size=0-Infinity circularity=0.00-1.00 show=Nothing clear include add stack");
TotROI=nResults;
roiManager("Show None");				// These lines find all the "tooth" pixels

for (i=0; i<TotROI; i++){
//	showStatus("Enigma currently running on ROI:"+i+" of "+TotROI+" for "+cratLoc[n]+"_"+craterID[n]);	// Status bar
//	showProgress(i/TotROI);
	roiManager("Select", i);			// Select the tooth pixels on relevant slice
	run("Convex Hull");				// Run the convex hull
	run("Invert", "slice");				// Invert this selection - effectively creating crater
	run("Select None");				// deselect everything
	run("Remove Outliers...", "radius=1 threshold=0 which=Bright"); // remove the thin 1 pixel border introduced by the hull
};							// repeat until all slices are complete
Stack.getStatistics(area,mean);				// gets the total "area" of stack (widthxheight x number of slices) along with "average" grey value
crater = (area*mean/255);				// multiplying A*M and dividing by 255 gives the % coverage by the tooth pixels
	crater_px = crater*res*res*res;				// allows for voxel conversion
	print(f, cratLoc[n]+"_"+craterID[n]+"\t"+crater+"\t"+crater_px);
//	showStatus("Enigma is complete for "+cratLoc[n]+"_"+craterID[n]);
};



}