//Graph Assistant 5-Aug-2017 
//05-Aug-2017 	error bars were not correct or missing when columns were only filled partially, this is fixed
//06-Mar-2016    error bars were not compatible with new ImageJ version > 1.48, this is fixed
//   <https://sils.fnwi.uva.nl/bcb/objectj/examples/Coli-Inspector/utils/>

macro "Unused Tool -" {}
var menuItems = getTitles("Axis Dia FluorTotal *Center");

var dCmds = newMenu("Graph Assistant Menu Tool", menuItems);
var k90 = 1.645, k95 = 1.9600, k98 = 2.3264, k99 = 2.5758;
var xAxis = "", yAxis = "";
var myList= "Axis-Min=0.0\nAxis-Max=3.5\nAxis-Bin=0.1\nDia-Min=0\nDia-Max=1.5";
var binnedXvalues, binnedYvalues, binnedErrorBars;

menuItems = getTitles("Dia FluorTotal *Center");

//return array of selected column titles plus chooser
function getTitles(titleSelection){
	patterns = split(titleSelection);
	titles = "";
	for (col = 1; col <= ojNColumns(); col++){
		title = ojGetColumnTitle(col);
		found = false;
		for (jj = 0; jj < patterns.length; jj++){
			if (ojMatches(title, patterns[jj]))
				found = true;
		}
		if (found)
			titles = titles + " " + title;
	}
	titles = titles + " - Axes_and_Ranges Qualify_ROI Disqualify_ROI Qualify_All - Custom_Columns";
	return split(titles);
}



function doPlot(xTitle, yTitle){
	requires("1.50");
	confid = k95;                              // 95% confidence
	List.setList(myList);
	minX = parseFloat(List.get(xTitle + "-Min"));
	maxX = parseFloat(List.get(xTitle + "-Max"));
	minY = parseFloat(List.get(yTitle + "-Min"));
	maxY = parseFloat(List.get(yTitle + "-Max"));
	binWidth = parseFloat(List.get(xTitle + "-Bin"));
	if (isNaN(minX)) minX = 0;
	if (isNaN(maxX)) maxX = 10;
	if (isNaN(minY)) minY = minOf(0, 1.05*ojGetStatistics(yTitle, "min"));
	if (isNaN(maxY)) maxY = 1.05 * ojGetStatistics(yTitle, "max");
	if (isNaN(binWidth)) binWidth = 0;
	plotTitle = yTitle + " vs " + xTitle;
	count = ojGetStatistics(xTitle, "count");
	if (ojGetStatistics(yTitle, "count") <1 ) {
		showMessageWithCancel("Zero qualified entries - Qualify All?");
		for (obj = 1; obj<= ojNObjects(); obj++) 
			ojQualify(obj, true);
		exit;
	}
	xValues = newArray(count);
	yValues	 = newArray(count);
	nQualified = 0;//num qualified objects
	for (obj = 1; obj<= ojNObjects(); obj++){
		if (ojQualified(obj)){
			valx = ojResult(xTitle, obj);
			valy = ojResult(yTitle, obj);
			if (!isNaN(valx + valy)){
				xValues[nQualified] = ojResult(xTitle, obj);
				yValues[nQualified++] = ojResult(yTitle, obj);
			}
		}
	}
	if (nQualified == 0) exit("zero objects");
	
	xValues = Array.trim(xValues, nQualified);//5.8.2017
	yValues = Array.trim(yValues, nQualified);

	if (binWidth > 0)
		getConfidences(xValues, yValues, binWidth, minX);
	
	while(isOpen(plotTitle)){
		selectImage(plotTitle); 
		close;
	}
	//if (binWidth > 0) rearranged 6-mar-2016
	//	Plot.create(plotTitle, xTitle, yTitle, binnedXvalues, binnedYvalues);
	//else
	Plot.create(plotTitle, xTitle, yTitle);
	Plot.setLimits(minX, maxX, minY, maxY);
	Plot.setColor("blue"); 
	Plot.add("dots", xValues, yValues);
	Plot.setColor("red");
	Plot.setLineWidth(2); 
	if (binWidth > 0){
		Plot.add("line", binnedXvalues, binnedYvalues);
		Plot.add("error bars", binnedErrorBars);
	}
	Plot.show;
}  
  

function qualifyRoi(inside){
	name = getTitle;
	if (endsWith(name, ".tif")){
		name = substring(name, 0, lengthOf(name) - 4);
	}
	titles = split(name);//e.g. Axis vs Dia
	xTitle = titles[titles.length-1];
	yTitle = titles[0];
	ojSetColumnProperty(xTitle, "sort", 0);//otherwise it crashes

	for (obj = 1; obj <= ojNObjects(); obj++){
		xval = ojResult(xTitle, obj);
		yval = ojResult(yTitle, obj);
		toUnscaled(xval, yval);
		flag = selectionContains(xval, yval)^!inside;
		ojQualify(obj, flag);
	}
	if (!endsWith(getTitle, ".tif")){
		showMessageWithCancel("Rebuild " + getTitle + " Plot ?");
		columns = split(getTitle);
		doPlot(columns[columns.length - 1], columns[0]);
	}
}


macro "Graph Assistant Menu Tool - C037T0c12GT8c12A" {

	if(xAxis == "")
		xAxis = ojGetColumnTitle(1);
	if(yAxis == "")
		yAxis = ojGetColumnTitle(2);
	ojRequires("1.01w");
	title = getArgument();
	setKeyDown("");
	if (ojColumnNumber(title) > 0){	
		doPlot(xAxis, title);
	}
	else if (title == "Axes_and_Ranges"){
		axesAndRangesDialog();
	}
	else if (title == "Qualify_ROI")
		qualifyRoi(true);
	else if (title == "Disqualify_ROI")
		qualifyRoi(false);
	else if (title == "Qualify_All"){
		for (obj = 1; obj <= ojNObjects(); obj++)
			ojQualify(obj, true);
		ojSetColumnProperty(xAxis, "sort", 0);//updates statistics
	}
	else if (title == "Custom_Columns"){
		makeCustomColumns();
	}
}

function axesAndRangesDialog(){
	nCols = ojNColumns();
	colTitles = newArray(nCols);
	for (jj = 0; jj < nCols; jj++){
		title = ojGetColumnTitle(jj+1);
		colTitles[jj] = title;
	}
	Dialog.create("Choose Columns for Scatter Plot");
	Dialog.addChoice("Use for X-Axis:", colTitles, xAxis);
	Dialog.addChoice("Use for Y-Axis:", colTitles, yAxis);
	Dialog.show();


	xAxis = Dialog.getChoice();
	yAxis = Dialog.getChoice();
	
	
	List.setList(myList);
	Dialog.create("Choose Plot Parameters");
	Dialog.addMessage("-- Horizontal --");
	key1 = xAxis+ "-Min";
	value = parseFloat(List.get(key1));
	Dialog.addNumber(key1, value);
	
	key2 = xAxis+ "-Max";
	value = parseFloat(List.get(key2));
	Dialog.addNumber(key2, value);

	key3 = xAxis+ "-Bin";
	value = parseFloat(List.get(key3));
	Dialog.addNumber(key3, value);

	Dialog.addMessage("-- Vertical --");
	key4 = yAxis+ "-Min";
	value = parseFloat(List.get(key4));
	Dialog.addNumber(key4, value);

	key5 = yAxis+ "-Max";
	value = parseFloat(List.get(key5));
	Dialog.addNumber(key5, value);
	
	Dialog.addCheckbox("Create Plot Now", true);
	Dialog.show();	
	val1 = Dialog.getNumber();
	val2 = Dialog.getNumber();
	val3 = Dialog.getNumber();
	val4 = Dialog.getNumber();
	val5 = Dialog.getNumber();
	List.set(key1, val1);
	List.set(key2, val2);
	List.set(key3, val3);
	List.set(key4, val4);
	List.set(key5, val5);
	myList=List.getList;
	plotNow = Dialog.getCheckbox;	
	if (plotNow) 
		doPlot(xAxis, yAxis);
}


//calculates globals binnedXvalues, binnedYvalues, binnedErrorBars;

function getConfidences(xValues, yValues, binWidth, xStart){//eg binWidth = 0.5
	rankPArr = Array.rankPositions(xValues);
	count = 0;
	group = newArray (9000);
	
	means = newArray(500);
	stDevs = newArray(500);
	e95 =  newArray(500);
	counts = newArray(500);
	midBins = newArray(500);
	nBins = 0;//need not be contiguous			
	for (jj = 0; jj < nQualified; jj++){
		index = rankPArr[jj];//start with smallest value
		xVal = xValues[index];
		yVal =  yValues[index];
		bin = floor((xVal - xStart)/binWidth);//0, 1, 2...
		full =	(jj == (nQualified -1));
		if (!full){
			nextIndex = rankPArr[jj + 1];
			nextBin = floor((xValues[nextIndex] - xStart)/binWidth);
			full = bin != nextBin;
		}
		group[count++] = yVal;
		if (full && count > 1){
			showProgress(jj/nQualified);
			group = Array.trim(group, count);
			Array.getStatistics(group, min, max, mean, stdDev);
			midBin = xStart + binWidth * bin + binWidth/2;
			midBins[nBins] = midBin;//0, 0.5, 1
			means[nBins] = mean;
			stDevs[nBins] = stdDev;
			e95[nBins] = confid*stdDev/sqrt(count); 
			counts[nBins] = count;	
			group = newArray (9000);
			count = 0;
			nBins++;
		}
		if (full) count =0;

	}
	midBins = Array.trim(midBins, nBins);
	means = Array.trim(means, nBins);
	stDevs = Array.trim(stDevs, nBins);
	e95 = Array.trim(e95, nBins);
	counts = Array.trim(counts, nBins);
	binnedXvalues = midBins;
	binnedYvalues = means;
	binnedErrorBars = e95;


}

//age = Td * log(1 - 0.5 * rank/nCells)/log(0.5);
function makeCustomColumns(){
	showProgress(0.1);

	ojHideResults();
	ojSetColumnProperty("Axis", "sort", 0);//updates statistics
	nCells = ojNObjects();
	createAgeColumn("Axis", "Age");
	createAgeColumn("Dia", "AgeDia");
	showProgress(0.2);
	ojInitColumn("RingFraction");
	for(obj = 1; obj <= nCells; obj++){
		ringFract = 100 * ojResult("FCPlus", obj) /ojResult("FluorTotal", obj);
		ojSetResult("RingFraction", obj, ringFract);
	}
	
	
	ojInitColumn("ConcTotal");
	for(obj = 1; obj <= nCells; obj++){
		conc = ojResult("FluorTotal", obj) /ojResult("Volume", obj);
		ojSetResult("ConcTotal", obj, conc);
	}
	showProgress(0.6);

	ojInitColumn("WallArea");
	for(obj = 1; obj <= nCells; obj++){
		axis = ojResult("Axis", obj);
		dia = ojResult("Dia", obj);
		cylinderArea =  (axis- dia) * dia * PI;
		ballArea =  (dia/2) *(dia/2) * 4 * PI;
		cellArea = cylinderArea + ballArea;
		ojSetResult("WallArea", obj, cellArea);
	}
	showProgress(0.9);

	ojInitColumn("ConcWall");
	for(obj = 1; obj <= nCells; obj++){
		conc = ojResult("FluorTotal", obj) /ojResult("WallArea", obj);
		ojSetResult("ConcWall", obj, conc);
	}
	ojShowResults();
}
 
function createAgeColumn(fromTitle, toTitle){
	nCells = ojNObjects();
	ojInitColumn(toTitle);
	axes = newArray(nCells);
	for (obj = 1; obj <= nCells; obj++)
		axes[obj-1] = ojResult(fromTitle, obj);
	rankPosArr = Array.rankPositions(axes);
	ranks = Array.rankPositions(rankPosArr);
	Td = 100;
	for (obj = 1; obj <= nCells; obj++){
		relRank = ranks[obj-1] / nCells;
		age = Td * log(1 - 0.5 * relRank) / log(0.5);
		ojSetResult(toTitle, obj, age);
	}
}