// Intersect Method for quantitative metallugraphy, 
// by Christopher Coulon, www.gaiag.net, 
// for Ali Soleimani Dorcheh, Max Plank Institute for Intelligent Systems.
//
// This macro is intended to convert a large amount of data from one or more images  
// to a few meaningful quantitative parameters which helps (in this context)
// to understand the three dimensional characteristics of the microstructure of materials.
// A manual systematic image analysis can be very time consuming
// since it requires a person who can characterize the image visually
// via different technics. Computer aided image analysis not only enhances the efficiency
// of the method and avoids biased results due to operating errors,
// but also enables us to collect enough data to provide statistical conclusions.
//
// This work provides a macro which utilizes ImageJ tools and provides
// an automated quantitative method for phase analysis based on counting method.
// The targeted systems are multi-phase alloys which have different phases
// with distinct interphase boundaries. A system composed of two phases in Cr-Cr3Si
// system is used to evaluate the functionality of the macro in such systems. 
//
// For details of usage, please email Chris (chris@gaiag.net)
// or Ali (soleimani.ali@gmail.com).

var step =50;
var angle =0;
var name, title, path, dir, resultsDir;
var areaFraction;
var gridType;
var Width, Height, imageArea = 0;
var voxelW, voxelH;
var lineCount;
var lineTotal = 0, ROItotal = 0;
var invertImage = true;
var timeString;
var ANDlinesPercent = 0, ROIpercent = 0;
var distance, unit;
var totalSum, totalLineCount, noEdgeSum, noEdgeLineCount, edgeParticles;
var intersections, wEdgeLines, noEdgeLines;
var options = false, start = true;

macro "Grid Tool - C00cT1f1aG" {}

macro "Grid Tool Options" {

	path = File.openDialog("Select a File");
	dir = File.getParent(path);
	name = File.getName(path);
	list = getFileList(dir);
	
	start = true;
	
	for (i=0; i<list.length; i++){
		setBatchMode(true);
		if (!endsWith(list[i], "/")) {
			open (list[i]);
			processImage();
			run("Close All");
			closeResults();
		}
		setBatchMode(false);
	}
}

function processImage() {
		
	getVoxelSize(width, height, depth, unit);	
//	write ("width = " + width + "  height = " + height + "  depth = " + depth + "  unit= " + unit);

	distance = 1 / width;
//	write ("distance = " + distance);

	run("Set Scale...", "distance=" + distance + " known=1 pixel=1 unit=" + unit +" global");

	Width = getWidth();
	Height = getHeight();
	imageArea = Width * Height;
	
	w = Width * width;
	h = Height * height;
	
	if (start) {
		gridOptions();
	}
	
	start = false;
	
//	if (options) gridOptions();
	
	setup();
	
	write ("Calibrated Dimensions in " + unit + "\nImage Width = " + w);
	write ("Image Height = " + h);
	write ("Image Area = " + w * h);
	write ("Pixel Dimensions" + "\nWidth = " + Width);
	write ("Height = " + Height);	
	write("Image Area", imageArea);	
	
	write ("grid step = " + step);
	write ("angle = " + angle);
	
	drawLines(); 
	
	processGrid(); 		
	
	selectWindow("Log");

	run("Set Scale...", "distance=" + distance + " known=1 pixel=1 unit=" + unit);

}

macro "Close All Tool - C00cT1f1aC" {}

macro "Close All Tool Options" {
	run("Close All");
}

macro "Close Results Tool - C00cT1f1aR" {}

macro "Close Results Tool Options" {
	closeResults();
}

function drawLines() {
	max = maxOf(getWidth,getHeight);
	a = PI*angle/180;
				
	xc=getWidth/2 + 2*step*(random-0.5);
	yc=getHeight/2 + 2*step*(random-0.5);
	
	Overlay.remove;

	for (i=-1*floor(max/step);i<max/step;i++) {
		x = xc-max*cos(a)+i*step*sin(a);
		y = yc+max*sin(a)+i*step*cos(a);
		Overlay.moveTo(x, y);
		x = xc+max*cos(a)+i*step*sin(a);
		y = yc-max*sin(a)+i*step*cos(a);
		Overlay.lineTo(x, y);		
	}
	Overlay.show;		
}
	
function setup() {
	
	title = getTitle();
	
	dot = lastIndexOf(title, ".");
	name = substring(title, 0, dot);
//	write ("title: " + title + "\nname: " + name);
	
	path = getDirectory("image");

	dir = name + "_" + angle + "_results";
	resultsDir = path+dir+File.separator;
	File.makeDirectory(resultsDir);
	
	write ("********** " + title + " **********");
	
	timeString = getTimeString();
	
	write (timeString);
	
	setThreshold(127, 255);
	
	run("Set Measurements...", "area area_fraction limit add redirect=None decimal=5");
	run("Measure");

	ROItotal = getResult("Area", nResults - 1);
	ROIpercent = getResult("%Area", nResults - 1);
	write ("Total ROI Area: " + ROItotal);
	write ("ROI %Area = " + ROIpercent);
	
	resetThreshold();
	
}

function getTimeString() {
     MonthNames = newArray("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec");
     DayNames = newArray("Sun", "Mon","Tue","Wed","Thu","Fri","Sat");
     getDateAndTime(year, month, dayOfWeek, dayOfMonth, hour, minute, second, msec);
     TimeString ="Date: "+DayNames[dayOfWeek]+" ";
     if (dayOfMonth<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+dayOfMonth+"-"+MonthNames[month]+"-"+year+"\nTime: ";
     if (hour<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+hour+":";
     if (minute<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+minute+":";
     if (second<10) {TimeString = TimeString+"0";}
     TimeString = TimeString+second;	
     return TimeString;
}

function processGrid() {
	
	run("Duplicate...", "title=temp");
	run("Colors...", "foreground=black background=white selection=yellow");
	setBackgroundColor(255, 255, 255);
	run("Select All");
	run("Clear", "slice");
	run("Select None");
	run("Flatten");
	run("8-bit");
	setThreshold(0, 254);
	setOption("BlackBackground", false);
	run("Convert to Mask");
	
	setOption("BlackBackground", false);
	run("Skeletonize");
	
	tls = getImageID();
	
	setThreshold(255, 255);
	
	run("Set Measurements...", "  feret's limit redirect=None decimal=5");
	run("Analyze Particles...", "  clear");
	
	totalLineLength = 0;
	for (i = 0; i < nResults; i++) {
		totalLineLength += getResult("Feret", i);
	}
	
	totalLineCount = nResults;

//	totalLineLength = totalSum;
	write ("Total Grid Length: " + totalLineLength);

	selectImage(tls);
	resetThreshold();
	thisName = name+"_" + angle + "-TotalLineCount.tif";
	saveAs("tiff", resultsDir+thisName); 

	imageCalculator("AND create", title, thisName);

	getMode();
	
	setThreshold(129, 255);

	numberOfIntercepts = sumLines();

	ANDlines = totalSum;
	
	areaFraction = ROIpercent;
	
	resetThreshold();

	saveAs("tiff", resultsDir+name+"_" + angle + "-Line_count.tif"); 
	
	dumpImages();
	
//	LL = ANDlines / totalLineLength;
//	NL = number of intercepts (= non-edge touching + 1/2 *edge touching)/ totalLineLength;
//	PL = 2 * NL
//	SV= 2*PL
//	LA = (PI / 4) * Sv;

	LL = ANDlines / totalLineLength;
	NL = numberOfIntercepts / totalLineLength;
	PL = 2 * NL;
	SV = 2 * PL;
	LA = (PI / 4) * SV;	
	
	write ("\n********** Variable Check **********");
	write ("LL = ANDlines / totalLineLength");
	write ("LL = " + LL);
	write ("ANDlines = " + ANDlines);
	write ("totalLineLength = " + totalLineLength);
	write ("\nNL = number of intercepts");
	write ("(=non-edge touching+1/2*edge touching)/totalLineLength");
	write ("number of intercepts (NL) = " + NL);
	write ("non-edge touching = " + noEdgeLineCount);
	write ("1/2 *edge touching = " + (edgeParticles / 2));	
	write ("**************************************\n");
	
	run("Clear Results");
		
	n = nResults;
	
	setResult("Label", n, name);
	setResult("Phase %", n, areaFraction);
	setResult("LL", n, LL);
	setResult("NL", n, NL);
	setResult("PL", n, PL);
	setResult("SV", n, SV);	
	setResult("LA", n, LA);	
	
	selectWindow("Log");
	saveAs("Text", resultsDir + name + "_" + angle + "_Log.txt");
	
	selectWindow("Results");
	saveAs("Results", resultsDir + name + "_" + angle + "_measurments.xls");
	
	selectWindow(title);
	resetThreshold();
}

function getMode() {
	run("Set Measurements...", "  modal redirect=None decimal=5");
	run("Measure");
	
	mode = getResult("Mode", nResults - 1);

	if (mode > 0) run("Invert");
}

function sumLines() {
	
	run("Set Measurements...", "  feret's limit redirect=None decimal=5");
	run("Analyze Particles...", "  clear");
	
	nr = nResults;
	wEdgeLines = newArray(nr);
	
	totalSum = 0;
	for (i = 0; i < nr; i++) {
		totalSum += getResult("Feret", i);
		wEdgeLines[i] = getResult("Feret", i);
	}
	
	totalMean = totalSum / nr;
	sumDiff = 0;
	
	for (i = 0; i < nr; i++) {
		x = getResult("Feret", i);
		sumDiff += (x - totalMean)*(x - totalMean); 
	}
	
	totalSD = sqrt(sumDiff / (nr - 1));
	
	totalLineCount = nResults;
	
	selectWindow("Results");
//	saveAs("Results", resultsDir + name + "_allLines.xls");
	
	run("Analyze Particles...", "  exclude clear");
	
	nr = nResults;
	noEdgeLines = newArray(nr);
	
	noEdgeSum = 0;
	for (i = 0; i < nResults; i++) {
		noEdgeSum += getResult("Feret", i);
		noEdgeLines[i] = getResult("Feret", i);
	}
	
	noEdgeMean = noEdgeSum / nr;
	sumDiff = 0;
	
	for (i = 0; i < nr; i++) {
		x = getResult("Feret", i);
		sumDiff += (x - noEdgeMean)*(x - noEdgeMean); 
	}
	
	noEdgeSD = sqrt(sumDiff / (nr - 1));

	noEdgeLineCount = nResults;	
	edgeParticles = totalLineCount - noEdgeLineCount;
	
	numberOfIntercepts = noEdgeLineCount + (edgeParticles / 2);
	
	selectWindow("Results");
//	saveAs("Results", resultsDir + name + "_noEdgeLines.xls");
	
	write ("All intercepts = " + totalLineCount);
	write ("No edge intercepts = " + noEdgeLineCount);
	write ("Edge intercepts = " + edgeParticles);
	
	write ("\n********** All Intercepts **********");
	ne = totalLineCount * 2;
	write ("Total Overlap Length = " + totalSum);
	write ("Line Intercepts = " + totalLineCount); 
	write ("Intersections = " + ne);
	
	write ("Average Length All = " + totalMean);
	write ("Standard Deviation All = " + totalSD);
	
	write ("****************************************");	
	
	write ("\n*********** All with Half Edge Intercepts ***********");
	intersections = noEdgeLineCount * 2 + edgeParticles;
	write ("Total Overlap Length = " + totalSum);
	write ("Line Intercepts = " + numberOfIntercepts); 
	write ("Intersections = " + intersections);	
	write ("****************************************");
	
	write ("\n********** Only No Edge Intercepts **********");
	ne = noEdgeLineCount * 2;
	write ("No Edge Overlap Length = " + noEdgeSum);
	write ("Line Intercepts = " + noEdgeLineCount); 
	write ("Intersections = " + ne);
	
	write ("\nAverage Length non-Edge = " + noEdgeMean);
	write ("Standard Deviation non-Edge = " + noEdgeSD);

	write ("****************************************");	
	
	return numberOfIntercepts;
	
}

function gridOptions() {
	Dialog.create("Grid");
	Dialog.addNumber("Step (in pixels):", step);
	Dialog.addNumber("Angle (in degrees):", angle);
//	Dialog.addCheckbox("Set Grid for All", false);
	Dialog.show;
	
	step = Dialog.getNumber();
	angle = Dialog.getNumber();
//	options = Dialog.getCheckbox();
}

function dumpImages() {
	selectWindow(title);
	close("\\Others");
}

function closeResults() {
	if (isOpen("Results")) {
		selectWindow("Results");
		run("Close");
	}
	if (isOpen("ROI Manager")) {
		selectWindow("ROI Manager");
		run("Close");	
	}
	if (isOpen("Log")) {
		selectWindow("Log");
		run("Close");
	}	
	if (isOpen("Threshold")) {
		selectWindow("Threshold");
		run("Close");
	}			
	if (isOpen("Debug")) {
		selectWindow("Debug");
		run("Close");
	}			
}