// @int(label = "Min Diffraction spot size (pix^2)", value=7, description="Minimum number of pixels that consitutes a diffraction spot.") minSize
// @String(label = "Beamstop Status:", choices={"Retracted", "Inserted"}, style="radioButtonHorizontal", description="If a beamstop was used, the center of the pattern is calculated differently than if no beamstop was used.") beamStopStatus
// @String(label = "Binning Size:", choices={"Small", "Medium", "Large"}, style="radioButtonHorizontal", description="Try changing this if you are not getting good binning to your individual d-spacing measurements.") stepSize
// @String(label = "Display Modes:", choices={"Measurements", "Summary", "All"}, style="radioButtonHorizontal", description="Which windows do you want to keep open at the end?") displayModes
// @String(label = "Include d-spacing list?", choices={"Yes", "No"}, style="radioButtonHorizontal", description="Click YES if you want to enable this feature.") txtWindow

/* 
 *  This plugin finds the center of every spot in your diffraction pattern above a threshold   
 *  number of pixels, chosen by the user, and then finds the distance between each spot and the 
 *  direct beam. This scalar quantity represents the length of the g vector in reciprocal space. 
 *  It is the inverse of d, the interplanar spacing in real space.
 *  
 *  Your diffraction pattern should be collected using Gatan's Digital Micrograph 
 *  software (.dm3 file type), and your images should be properly calibrated in reciprocal space, 
 *  i.e. scale bar units display as 1/nm in Digital Micrograph. 
 *  
 *  Written by: Mak Koten, Ph.D. (2018)
 */

// Duplicate original image (this could be optional)...
run("Duplicate...", " ");
diffDuplicate=getImageID;

// Set up the thresholding.
thresholdSpots();

// Use analyze particles method to measure the position of each diffraction spot's center.
measureSpots();

/*
 * The user must help the plugin calculate the center of the direct beam. This is critical since
 * it is the position from which all d-spacing values are calculated. If the beamstop was retracted, 
 * then it is just the center spot. If the beamstop was inserted then it is the midpoint between any
 * two symmetric spots chosen by the user.
 */

if (beamStopStatus == "Retracted") {
	Dialog.create("Direct Beam Center");
	Dialog.addMessage("Enter the label that corre-\nsponds to the direct beam.\n ");
	Dialog.addNumber("             (000):", 10);
	Dialog.show();
	nextRow = Dialog.getNumber();
	dbRow = nextRow-1;
	directBeamX = getResult("X", dbRow);
	directBeamY = getResult("Y", dbRow);
} else if (beamStopStatus == "Inserted") {
	Dialog.create("Direct Beam Center");
	Dialog.addMessage("Enter two spots that are mirrored about the direct \nbeam, and are equidistant from its center.\n ");
	Dialog.addNumber("                    Spot 1:", 10);
	Dialog.addNumber("                    Spot 2:", 30);
	Dialog.addMessage("");
	Dialog.addNumber("                  No. of Artifacts:", 2);
	Dialog.addMessage("");
	Dialog.show();
	
	nextRowSpotOne = Dialog.getNumber();
	nextRowSpotTwo = Dialog.getNumber();
	ignoreSpots = Dialog.getNumber();
	rowSpotOne = nextRowSpotOne-1;
	rowSpotTwo = nextRowSpotTwo-1;

	spotOneX = getResult("X", rowSpotOne);
	spotOneY = getResult("Y", rowSpotOne);
	spotTwoX = getResult("X", rowSpotTwo);
	spotTwoY = getResult("Y", rowSpotTwo);

// Taking the average makes the order entered irrelevant
	directBeamX = (spotOneX+spotTwoX)/2;	
	directBeamY = (spotOneY+spotTwoY)/2;
} else {
	print("Must select BeamStop status.");
}

// Read the measurements made into arrays and calculate G and d for each spot.
label = newArray();
area = newArray();
xPos = newArray();
yPos = newArray();
G = newArray();
dSpace = newArray();

for (row=0; row<nResults; row++) {
	indexRow = row+1;
	a = getResult("Area",row);
	spotX = getResult("X",row);
	spotY = getResult("Y",row);
	g = sqrt((spotX-directBeamX)*(spotX-directBeamX)+(spotY-directBeamY)*(spotY-directBeamY));	// Thanks Pythagoras!
	d = 10/g;
	G = Array.concat(G, g);
	dSpace = Array.concat(dSpace, d);
	label = Array.concat(label, indexRow);
	area = Array.concat(area, a);
	xPos = Array.concat(xPos, spotX);
	yPos = Array.concat(yPos, spotY);
	setResult("G (1/nm)", row, g);
	setResult("d (A)", row, d);
}

// Update results table with new columns for G and d. ALL spots included here even artifacts.
updateResults();

// The lists are trimmed to exclude measurement artifacts.
if (beamStopStatus == "Retracted") {
	cut = nResults-1;
	}
	else if (beamStopStatus == "Inserted") {
	cut = nResults-ignoreSpots;
	}
	else {
	print("Forgot to enter number of spot to ignore");
}

// Which step size does the user want to use for binning?
if (stepSize == "Small") {
		stepSize = 0.05;
	}
	else if (stepSize == "Medium") {
		stepSize = 0.1;
	}
	else if (stepSize == "Large") {
		stepSize = 0.3;
	} else {
	print("Error: No bin size entered.");
}

// Close undesired windows or exit for selected display mode.
if (displayModes == "Summary") {
	selectWindow("Results");
	run("Close");
	selectImage(diffDuplicate);
	run("Close");
}
if (displayModes == "Measurements") {
	exit;
}

// Trim measurement artifacts from dSpace Array
dSpace = Array.sort(dSpace);
dSpace = Array.trim(dSpace, cut);

// Automatically assign and average individual measurements of d to a single value for plane families.
dSpaceAve = aveByFamily(dSpace, stepSize);
d_Spacing = sortReverse(dSpaceAve);

// Calculate the standard deviations for the average d spacings.
dSpaceSTDEV = dStandDev(dSpace, stepSize);
Three_Sigma = sortReverse(dSpaceSTDEV);

// Calculate the number of spots used for each average value.
numSpots = dAveNumber(dSpace, stepSize);
Spot_Count = sortReverse(numSpots);

// Plot the averaged values and the measured values on the same graph for visualization.
Plot.create("Comparison Plot", "Measured Spots", "Plane Spacings");
Plot.setColor("blue", "blue");
Plot.add("circle", dSpace);
Plot.setColor("red", "red");
Plot.add("boxes", dSpaceAve);

// Create Results Summary table (Index, # spots, d ave, d stdev).
Array.show("Results Summary in Angstroms (row numbers)", Spot_Count, d_Spacing, Three_Sigma);

// Create text window list?
if (txtWindow == "Yes") {
	printToText("d (A) = ", d_Spacing);
}

/*
 *************************
 * * * * FUNCTIONS * * * *
 *************************
 */

// This calls a sequence of commands that allows the user to optimize thresholding before analysis.
function thresholdSpots() {
	run("8-bit");
	run("Clear Results");
	resetThreshold();
	setAutoThreshold("Default dark");
	//setThreshold(130, 255);
	run("Threshold...");
	waitForUser("Adjust Threshold, and click OK when done.");
	if (isOpen("Threshold") == true) {
		selectWindow("Threshold");
		run("Close");
	}
}

// This calls a sequence of commands that will measure the position of the thresholded spots.
function measureSpots() {
	run("Set Measurements...", "area centroid redirect=None decimal=4");
	run("Analyze Particles...", "size=" + minSize + "-Infinity pixel circularity=0.0-1.00 show=[Overlay Masks] display exclude include");
	run("Labels...", "color=red font=28 show"); 
}

// This will average the measurements from a family of planes and return a single array with the average values
function aveByFamily(d0, step) {
	d1=newArray();
	d2=newArray();
	for (i=0; i<d0.length; i++) {
		d1=Array.concat(d1, d0[i]);
		dDiff=d0[i+1]-d0[i];
		if (dDiff>=step) {
			Array.getStatistics(d1, min, max, mean, std);
			d2=Array.concat(d2, mean);
			d1=newArray();
		}
		if (i==d0.length-2) {
			last=d0.length-1;
			d1=Array.concat(d1, d0[last]); 
			Array.getStatistics(d1, min, max, mean, std);
			d2=Array.concat(d2, mean);
			return d2;
			exit;
		}
	}
}

// This will return an array of 3*standard deviations for the averaged planes. 
function dStandDev(d0, step) {
	d1=newArray();
	d2=newArray();
	for (i=0; i<d0.length; i++) {
		d1=Array.concat(d1, d0[i]);
		dDiff=d0[i+1]-d0[i];
		if (dDiff>=step) {
			Array.getStatistics(d1, min, max, mean, std);
			d2=Array.concat(d2, 3*std);
			d1=newArray();
		}
		if (i==d0.length-2) {
			last=d0.length-1;
			d1=Array.concat(d1, d0[last]); 
			Array.getStatistics(d1, min, max, mean, std);
			d2=Array.concat(d2, 3*std);
			return d2;
			exit;
		}
	}
}

// This will return an array containing the number of spots averaged for each value. A psudo intensity... 
function dAveNumber(d0, step) {
	d1=newArray();
	n=0;
	for (i=0; i<d0.length; i++) {
		n++;
		dDiff=d0[i+1]-d0[i];
		if (dDiff>=step) {
			d1=Array.concat(d1, n);
			n=0;
		}
		if (i==d0.length-2) {
			d1=Array.concat(d1, n+1); 
			return d1;
			exit;
	}
}

// Print to Text Window
function printToText(name, a) {
	title1 = "Text Window";
	title2 = "["+title1+"]";
	f = title2;
	run("Text Window...", "name="+title2+" width=72 height=4 menu");
	print(f, "\\Update:"); // clears the window
	print(f, name);
	n=a.length;
	for (i=0; i<n; i++) {
		if (i==n-1) {
			print(f, d2s(a[i],2)+".");
			exit;
		}
		print(f, d2s(a[i],2)+", ");
	}
}

// This is just like Array.sort(), but it sorts the array by descending values rather than ascending.
function sortDescending(array1) {
	array2=newArray();
	array3=newArray();
	for (i=0; i<array1.length; i++) {
		negValue=-1*array1[i];
		array2=Array.concat(array2, negValue);
		array2=Array.sort(array2);
	}
	for (i=0; i<array2.length; i++) {
		posValue=-1*array2[i];
		array3=Array.concat(array3, posValue);
	}
	return array3;
}

// This will sort an array in reverse without the sort function.
function sortReverse(array1) {
	array2=newArray();
	l=array1.length;
	for (i=0; i<l; i++) {
		n=-1*i+l-1;
		array2=Array.concat(array2, array1[n]);
	}
	return array2;
}

// Sometimes useful for troubleshooting.
function list(name, a) {
    print(name); 
    Array.print(a);
}