//Umis_Array
//Read layout file
//Determine position of array within image
//Draw an ROI at each point
//Adjust position, store adjustment and apply to next position
//log and recall error values for each position
//Classify tissue at each point
//Measure greyscale at each position
//GoTo any indent within an array
//Go forward or backwards through the array
//Handles arrays in any orientation
//Correlates qBSE backscatter coefficient
//to nanoindentation elastic modulus

//Umis_Array v0.53 ImageJ macro for combined qBSE and nanoindentation analysis
//Copyright (C) 2007 Michael Doube
//m.doube at qmul.ac.uk
//2007-04-02 18:03 GMT

//This program is free software; you can redistribute it and/or
//modify it under the terms of the GNU General Public License
//as published by the Free Software Foundation (version 2).
//http://www.gnu.org/copyleft/gpl.html

//This program is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.

//You should have received a copy of the GNU General Public License
//along with this program; if not, write to the Free Software
//Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

//declare some global variables

var xyf;
var row;
var pos;
var width;
var height;
var nindents;
var radius;
var startx;
var starty;
var newx;
var newy;
var scalex;
var scaley;
var thetax;
var thetay;
var lengthy;
var lengthx;
var lines;
var imagelist;
var xerror;
var yerror;
var goto = -1;
var xerrors;
var yerrors;
var errorlogged;
var xfactor = 1;
var yfactor = 1;
var batch = -1;
var moduli;
var force_macro = -1;
var nimages = 0;
var xrefx;
var xrefy;
var yrefx;
var yrefy;
var win_size = 300;
var moduli;
var modulus_min = 0;
var modulus_max = 30;
var mean_greys;


//set the ROI radius - usually 15 pixels
var radius = 15;

//first macro starts here

macro "Setup Array... [b]" {
nimages = 0;
imagelist = newArray(3);
//open the qBSE image
open();
imagelist[0]=getImageID();
nimages++;
run("Set Scale...", "distance=1 known=1 pixel=1 unit=pixel global");
//open the array map
open();
imagelist[1]=getImageID();
nimages++;
run("Set Scale...", "distance=1 known=1 pixel=1 unit=pixel global");
sync = isOpen("Sync Windows 1.6");
if (sync == false) run("Sync Windows");

//Open the array file which must be a csv list like x,y,f
string = File.openAsString("");

//Initialise the results window
if(isOpen("Results")==false) row = 0;
else if (pos < nResults()) row =pos;
else row = nResults();

//initialise the position stepper
pos=-1;

//Split the CSV file into an array of lines
lines = split(string, "\n");
nindents = lengthOf(lines);

//initialise the error log arrays
xerrors = newArray(nindents);
yerrors = newArray(nindents);
errorlogged = newArray(nindents);

//find the dimensions of the indent array
dimension = lines[nindents-1];
dimension = split(dimension, "\,");
width = parseFloat(dimension[0]);
height = parseFloat(dimension[1]);

//Clear old selections
run("Select None");

//Locate important array landmarks on the image
//bottom left (start), bottom right (xref), top left (yref)

leftButton=16;
shift=1;
alt=8;
rightButton=23;
ctrl=2;
x2=-1;
y2=-1;
z2=-1;
flags2=-1;

//set the last value to be determined to be the loop escape value
yrefy = -1;

while (yrefy==-1) {
          getCursorLoc(x, y, z, flags);
          if (x!=x2 || y!=y2 || z!=z2 || flags!=flags2) {

// the "-radius" bits are to shift the ROI's to the centre of the cursor position
//since the x,y position of each ROI is determined by the top left corner of its bounding box

//capture the start position with left click plus ctrl
		if (flags == 18){
			startx = x-radius;
			starty = y-radius;
			print("start position captured",x, y,"\n");
		}		

//capture the bottom right indent point with left click plus alt
		if (flags == 24){
			xrefx = x-radius;
			xrefy = y-radius;
			print("x reference point captured",x,y,"\n");
		}		
//capture the top left indent point with left click plus shift
		if (flags == 17){
			yrefx = x-radius;
			yrefy = y-radius;
			print("y reference point captured",x,y,"\n");
		}
	}

	x2=x; y2=y; z2=z; flags2=flags;
	wait(10);
}

//work out the baseline angle
thetax = atan((starty-xrefy)/(xrefx-startx));
print ("baseline angle", thetax,"\n");

//work out the baseline length
lengthx = sqrt((xrefx-startx)*(xrefx-startx)+(xrefy-starty)*(xrefy-starty));
print ("Baseline length", lengthx, width, "\n");

//work out the mast angle
thetay = atan((starty-yrefy)/(yrefx-startx));
print ("Mast angle", thetay,"\n");

//work out the mast length
lengthy = cos((thetax)-(thetay)+(PI/2))*sqrt((yrefx-startx)*(yrefx-startx)+(yrefy-starty)*(yrefy-starty));
print("Mast length", abs(lengthy), height, "\n");

//work out a scale factor in x
scalex = (lengthx)/(width);

//work out a scale factor in y
scaley = abs((lengthy)/(height));

//handle top-right start positions
if (starty < yrefy){
	scaley = -1*scaley;
	scalex = -1*scalex;
}

//handle reflected axes
if ((starty > yrefy && startx > xrefx)||(starty < yrefy && startx < xrefx)){
	scalex = -1*scalex;
}
}
}

//macro - adjust the scale factors
macro "Adjust Scale... [u]" {
Dialog.create("Adjust Zoom");
Dialog.addNumber("Window Size", win_size);
Dialog.show();
win_size = Dialog.getNumber();
drawarray();
}

//macro - draw the whole array
macro "Draw array... [y]" {
Dialog.create("Draw Options");
array_styles = newArray("Outline", "Numbered", "Numbered Outline", "Filled");
Dialog.addChoice("Style", array_styles, array_styles[0]);
Dialog.show();
array_style = Dialog.getChoice();
batch = 1;
setBatchMode(true);
for (pos=0; pos<nindents;pos++){
	drawarray();
	if (array_style == array_styles[0]) run("Draw");
	else if (array_style == array_styles[1]){
		//print the number
		drawString(pos+1,startx+newx+xerror+radius,starty-newy-yerror+radius);
	}
	else if (array_style == array_styles[2]){
		run("Draw");
		drawString(pos+1,startx+newx+xerror,starty-newy-yerror);	
	}
	else fill();
}
batch = -1;
setBatchMode(false);
run("Select None");
pos = 0;
}

//macro - advance to the next indent
macro "Next Indent... [l]" {
setBatchMode(true);
poserror();
if (pos < nindents-1) pos++; else exit("End of indentation array!");
drawarray();
setBatchMode(false);
}

//macro - go to the previous indent
macro "Previous Indent... [j]" {
setBatchMode(true);
poserror();
if (pos > 0) pos--; else exit("Beginning of indentation array!");
drawarray();
setBatchMode(false);
}

//macro: make a dialog box where the user can enter the indent they want to go to
macro "Go To Indent... [g]" {
	poserror();
	Dialog.create("Go To Indent");
	Dialog.addString("Indent No.", "");
	Dialog.addMessage("Indent range 1-"+nindents+"\n");
	Dialog.show();
	goto = parseFloat(Dialog.getString());
	if (goto > nindents || goto < 1) showMessage("Oops!", "Please enter a value between 1 and "+nindents+"\n");
	else {
		goto = goto - 1;
		drawarray();
	}
}

//macro - set tissue type, measure pixels and record results
macro "Tissue type... [t]" {

//draw a plot
if (!isOpen("Scatter Plot")){
parse_moduli();
mean_greys = newArray(nindents);
}

if(isOpen("Results")==false) row = 0;
//update the row in the results table that relates to the indent
if (pos <= nResults) row = pos;
	xyf= split(lines[pos], "\,");
	for (p=0; p<lengthOf(xyf); p++){
		setResult("N", row, pos+1);
		setResult("X", row, xyf[0]);
		setResult("Y", row, xyf[1]);
		setResult("Force", row, xyf[2]);
	}
selectImage(imagelist[0]);
//Make a dialog box where the user can enter tissue type
//and go to the next indent
	Dialog.create("Classify Tissue");
	Dialog.addMessage("Indent no.: "+(pos+1)+"\n\n0 Lacuna\n1 Hyaline articular cartilage\n2 Articular calcified cartilage\n3 Cement line\n4 Bone\n5 Marrow space");
	Dialog.addString("Tissue Type", "");
	Dialog.show();
	type = Dialog.getString();
	setResult("Tissue", row, type);

//collect data on pixels in ROI
	getRawStatistics(count, mean, min, max, std, histogram);

//add the mean value to the array for a plot
mean_greys[pos] = mean;

//set the threshold (qBSE images use lots of 0's for non-mineralised tissue)
        	lower = 1;
	upper = 255;

//work out the non-zero area fraction
	background = 0;
       	for (k=0; k<lower; k++) {
            		background = background + histogram[k];
            		fraction = (count - background) / count;
		setResult("AF", row, fraction);
       	}

//work out the mean mineral density of the mineralised fraction
	cumulative = 0;
	av= 0;
	for (j=lower; j<upper; j++) {
		cumulative = cumulative + histogram[j];
		av = av + j * histogram[j];
		bmmd = av / cumulative;
		setResult("Mean Matrix", row, bmmd);
	}

//work out the modal mineral density of the mineralised fraction
	mode = 0;
	for (p=lower; p<upper; p++) {
		if (mode < histogram[p]) {
			mode = histogram[p];
       			setResult("Mode Matrix", row, p);
			modep = p;
		}
	setResult("Mean ROI", row, mean);
	}
print ('row', row, 'pos', pos);
	updateResults;

//plot the result
Plot.create("Scatter Plot", "Mean Backscatter Coefficient", "Elastic Modulus");
Plot.setLimits(0, 255, 0, modulus_max);
Plot.setColor("green");
graph_moduli = newArray(nindents);
graph_greys = newArray(nindents);
for (n = 0; n<=pos; n++){
graph_moduli[n] = moduli[n];
graph_greys[n] = mean_greys[n];
}
Plot.add("boxes", graph_greys, graph_moduli);
Plot.update();


poserror();
if (pos < nindents-1) pos++;
else {exit("End of indentation array!");}
drawarray();
}
}

//Macro: set the fill colour based on the indentation modulus
//requires a CSV list in (n,mean,sd) format

macro "Colour Code Array... [f]"{

force_macro = 1;

//parse the moduli file
parse_moduli();

//ask the user for the array spacing
Dialog.create("Array Spacing");
Dialog.addMessage("Please enter\nthe array spacing");
Dialog.addNumber("x (microns)", 20);
Dialog.addNumber("y (microns)", 20);
Dialog.show;
xspacing  = parseFloat(Dialog.getNumber());
yspacing = parseFloat(Dialog.getNumber());


selectImage(imagelist[0]);
name = getTitle();
w = getWidth();
h=getHeight();
newImage("moduli_"+name, "8-bit Black", w, h, 1);
imagelist[2] = getImageID();
nimages++;
scale = 250/(modulus_max-modulus_min);

//Draw the array
batch = 1;
setBatchMode(true);
for (pos=0; pos<nindents;pos++){
    if (20+(moduli[pos]-modulus_min)*scale < 255 ) setColor(20+(moduli[pos]-modulus_min)*scale);
    else setColor(255);
	drawarray();
	fill();
}
batch = -1;
setBatchMode(false);
pos = 0;

run("Select None");
run("Interactive 3D Surface Plot");

//only plot the values from the array area, not the background black
makePolygon(startx-radius, starty+radius, xrefx+radius, xrefy+radius, yrefx+xrefx-startx+radius, yrefy+xrefy-starty-radius, yrefx-radius, yrefy-radius);
run("Interactive 3D Surface Plot");

force_macro = -1;
}


//poserror() - determine the error at each position
function poserror() {
if (pos > -1){

//measure the manual adjustment to add to the auto adjust value
getSelectionBounds(adjposx,adjposy,wi, he);

//recall the logged error values for that position if they exist
if (errorlogged[pos] == -8){
	xerror = xerrors[pos];
	yerror = yerrors[pos];
}

xerroradj = adjposx-(startx+floor(newx)+xerror);
yerroradj = adjposy-(starty-floor(newy)-yerror);

//if the selection is accidentally deselected, fix the way-out error value
if (adjposx==0 && adjposy==0){
	xerroradj = 0;
	yerroradj = -1;
}

xerror = xerror+xerroradj;
//yerror needs to be adjusted by -1 to correct for float to integer rounding errors
yerror = yerror-yerroradj-1;
if (yerror == -1) yerror = 0;

//log the x and y errors
xerrors[pos] = xerror;
yerrors[pos] = yerror;
errorlogged[pos] = -8;
print(pos, xerrors[pos], yerrors[pos]);
}

}

//drawarray() function

function drawarray() {

setBatchMode(true);

//Get each x,y coordinate and display it
	if (goto != -1) pos = goto;
	xyf= split(lines[pos], "\,");

//set the x scale
	x = (xyf[0])*(scalex);

//set the y scale
	y = (xyf[1])*(scaley);
	
//work out the rotated positions
	newx = (x-y*(tan(thetax)))*(cos(thetax));
	newy = y/cos(thetax)+sin(thetax)*(x-y*tan(thetax));

//work out the skew
	skewx = cos(thetax)*y*tan(thetax+(PI/2)-thetay);
	skewy = sin(thetax)*y*tan(thetax+(PI/2)-thetay);

//adjust the positions
	newx = (newx + skewx)*xfactor;
	newy = (newy + skewy)*yfactor;

//reset the error values when returning to the origin (1st indent)
if (pos == 0){
	xerroradj = 0;
	yerroradj = -1;
	xerror = 0;
	yerror = 0;
}

//recall the logged error values for the position if they exist
if (errorlogged[pos] == -8){
	xerror = xerrors[pos];
	yerror = yerrors[pos];
}

//Zoom and centre the window
for (n = 0; n<nimages; n++){
	if (batch < 0){
		selectImage(imagelist[n]);
		makeRectangle(startx+newx+xerror-win_size/2+radius, starty-newy-yerror-win_size/2+radius, win_size, win_size);
		run("To Selection");
}

//if we are running the modulus macro then draw a polygon
if (force_macro > 0){

xp = startx+newx+radius+xerror;
yp = starty-newy+radius-yerror;

h = abs(yspacing * lengthy / height);
l = xspacing * lengthx / width;

x0 = xp - (h * cos(thetay) / 2) - (l * cos(thetax) / 2);
y0 = yp - (h * sin(thetay) / 2) - (l * sin(thetax) / 2); 

x1 = x0 - h*cos(thetay);
y1 = y0 + h*sin(thetay);

x3 = x0 + l*cos(thetax);
y3 = y0 - l*sin(thetax);

x2 = x1 + x3 - x0;
y2 = y1 + y3 - y0;

makePolygon(x0,y0, x1, y1, x2, y2, x3, y3);
} else {

//draw a circular ROI at the rotated, skewed and error adjusted position
	makeOval(startx+newx+xerror,starty-newy-yerror,2*radius,2*radius);
}

}

//Tell us which indent we are looking at
if (nimages == 3) showStatus("Indent: "+(pos+1)+", x="+xyf[0]+", y="+xyf[1]+", E="+round(10*moduli[pos])/10);
else showStatus("Indent: "+(pos+1)+", x="+xyf[0]+", y="+xyf[1]);
	
//reset the GoTo value
goto = -1;

}

function parse_moduli() {
//parse the moduli file
	modulus_min = 99999999;
	modulus_max = 0;
	moduli_file = File.openAsString("");
	moduli_lines = split(moduli_file, "\n");
	moduli = newArray(nindents);
	for (m=0; m<lengthOf(moduli_lines); m++){
		moduli_n= split(moduli_lines[m], "\,");
		if (isNaN(moduli_n[1]) == true) moduli[m] = 0; else moduli[m] = parseFloat(moduli_n[1]);
		if (moduli[m] < modulus_min) modulus_min = moduli[m];
		if (moduli[m] > modulus_max) modulus_max = moduli[m];
	}
}