//Author: N. Vischer
//18.04.24, 8:26

var 
dia = 54, //diameter of a well in pixels
thresh=70, //suggested, will be changed by macro
background = 0,// will be changed by macro
plotMinY=70, //lower plot limit 
plotMaxY=160, //upper plot limit
show3Points = true, //paint t20, t50, t80 points plus text in plot?
includeArea = false, //show area in tableA?
flagColor = "magenta",
;

var //internally used:
baseName ="", basePath = "", tabA = "", tabB = "", labelFlag = 0,
;

// Sets 4 circles close to the corners so the user can drag 
//them to the correct place
macro "Suggest 4 Circles[1]"{
	requires("1.54j");
	if(baseName != "")
		showMessageWithCancel("Forget previous register");
	baseName = "";	
	title = getTitle;
	if(!endsWith(title, ".tif"))
		exit("title should end with 'tif'");
	if(nSlices  == 1)
		exit("Timelapse stack expected");
	if(Overlay.size > 0)
		showMessageWithCancel("Forget current overlay?");
	Overlay.clear;
	roiManager("reset");
	off = 0.15;
	makeOval(getWidth * off, getHeight * off, dia, dia);
	Overlay.addSelection;
	makeOval(getWidth * (1-off), getHeight * off, dia, dia);
	Overlay.addSelection;
	makeOval(getWidth * off, getHeight * (1-off), dia, dia);
	Overlay.addSelection;
	makeOval(getWidth * (1-off), getHeight * (1-off), dia, dia);
	Overlay.addSelection;
	run("Overlay Options...", "stroke=green width=0 fill=none set show");//show numbers
	run("Select None");
	Overlay.useNamesAsLabels(false);
	run("Labels...", "color=white font=12 show draw");
	Overlay.show;
}

//if 4 circles exist, the remaining 92 will be created by 2D interpolation
//if >= 96 circles exist, corner circles remain, rest is recalculated
macro "Create 96 Circles [2]"{
	resetThreshold;
	if(Overlay.size > 96)
		showMessageWithCancel("OK to delete overlays > #96?");
	xx = newArray(4);
	yy = newArray(4);
	if(Overlay.size == 4)
		corners = newArray(0, 1, 2, 3);//use these 4
	else if(Overlay.size >= 96)
		corners = newArray(0, 11, 84, 95);//use 4 corners and recalculate the rest
	else
		exit("Expected: 4 or min. 96 circles");
	for(jj=0; jj < 4; jj++){
		Overlay.activateSelection(corners[jj]);
		getSelectionBounds(x, y, width, height);
		xx[jj] = x + width/2;//centers
		yy[jj] = y + height/2;
	}
	
	Overlay.clear;
	//Calculate well positions via 2D-interpolation
	leftX = newArray(xx[0], xx[2]);
	leftY = newArray(yy[0], yy[2]);
	rightX = newArray(xx[1], xx[3]);
	rightY = newArray(yy[1], yy[3]);
	
	leftX = Array.resample(leftX, 8);
	leftY = Array.resample(leftY, 8);
	rightX = Array.resample(rightX, 8);
	rightY = Array.resample(rightY, 8);
	
	for(row = 0; row < 8; row++){
		xRow = newArray(leftX[row], rightX[row]);
		xRow = Array.resample(xRow, 12);
		yRow = newArray(leftY[row], rightY[row]);	
		yRow = Array.resample(yRow, 12);
		
		for(col = 0; col < 12; col++){
			makeOval(xRow[col] -dia/2, yRow[col] - dia/2, dia, dia);
			Roi.setStrokeColor("green");
			run("Add Selection...");
			Overlay.setPosition(0);//show on all slices
		}
	}
	run("Select None");	
	Overlay.useNamesAsLabels(false);
	Overlay.drawLabels(true); 	
	run("Show Overlay");
}

//Expects user to set global thresh for later
macro "Define Threshold [3]"{
	run("Threshold...");
	setThreshold(thresh, 255);
	waitForUser("Adjust theshold, then click OK");
	getThreshold(thresh, upper);
	showMessage("Threshold = "+ thresh +" will be used");
	resetThreshold;
}


//Creates TableA showing measurements for each circle in time
//timestamp: see comments in function makeTimeTable()
macro "Calc Wells[4]"{
	if(Overlay.size < 96)
		exit("96 circles expected");
	if(Overlay.size > 96){
		showMessageWithCancel("OK to delete Overlay > 96?");
		while(Overlay.size > 96)
			Overlay.removeSelection(96);
	}
	createTableA();
	tblPath = getDirectory("image") + Table.title;
	if(File.exists(tblPath)){
		showMessageWithCancel("Overwrite TableA on disk?");
	}
	saveAs("Results", tblPath);
}

//Of each circle, the mean of pixels above thresh 
//is measured in time and displayed in a table
function createTableA(){
	Overlay.hide;
	makeTimeTable();
	for(well = 0; well< 96; well++){//travel through wells
		Overlay.activateSelection(well);
		setThreshold(thresh, 255);
		for(frame = 1; frame <= nSlices; frame++){ //travel through frames
			setSlice(frame);
			mean = getValue("Mean limit");//only include pixels above thresh.
			columnTitle = "Mean_"+ (well + 1);
			Table.set(columnTitle, frame-1, mean, tabA);
			
			area = getValue("Area limit");//only include pixels above thresh.
			columnTitle = "Area_"+ (well + 1);
			if(includeArea)
				Table.set(columnTitle, frame-1, area, tabA);
		}
		run("Select None");
	}
	Table.update(tabA);
	setSlice(1);
	resetThreshold;
	Overlay.show;//09.04.24, 11:26
}

//Uses data of TableA to create a stack of 96 plots
//Each plot contins brightness vs time.
//If global show3Points=true, 3 statistical points (magenta) 
//are added with annotation
macro "TableA-to-Plots [5]" {
	close("Collective Plot*");
	close("Results");
	showTable(tabA);

	xx1 = Table.getColumn("Time[min]", tabA);
	Plot.create("Collective Plot", "Time [min]", "Brightness");
	Plot.setFrameSize(350, 300);
	cTitles = split(Table.headings);
	annotations = newArray(96);
	stats = newArray(96);
	well = 0;
	for(col=0; col <cTitles.length; col++){
		cTitle = cTitles[col];
		if(startsWith(cTitle, "Mean_")){
			well++;
     		Plot.setColor("black");
			Plot.setLineWidth(2);
  			yy1 = Table.getColumn(cTitle, tabA);
			Plot.add("connected", xx1, yy1);
			ranks = Array.rankPositions(yy1);
			Array.reverse(ranks);
			maxPos = ranks[0];
			maxFl = yy1[maxPos];
			
			leftPart = Array.trim(yy1, maxPos +1);// part left of peak
			ranks = Array.rankPositions(leftPart);
			minPos = ranks[0];
			minFl = yy1[minPos];
			amp = maxFl - minFl;

			yyPart = Array.slice(yy1, minPos, maxPos+1);// Part min to max
			xxPart = Array.slice(xx1, minPos, maxPos+1);// 
						
			t20 = 0; t50 = 0; t80 = 0;
			f20 = 0; f50 = 0; f80 = 0;
			if(maxPos > 0){
				hiResX = Array.resample(xxPart, 200);
				hiResY = Array.resample(yyPart, 200);
				for(jj=0; jj<200; jj++){
					if(hiResY[jj] >= (minFl + amp*0.2) && t20==0){
						t20 = hiResX[jj];
						f20 = hiResY[jj];
					}
					if(hiResY[jj] >= (minFl + amp*0.5) && t50==0){
						t50 = hiResX[jj];
						f50 = hiResY[jj];
					}
					if(hiResY[jj] >= (minFl + amp*0.8) && t80==0){
						t80 = hiResX[jj];
						f80 = hiResY[jj];
					}
				}
			}

			if(show3Points){
				Plot.setColor("magenta");
				xx = newArray(t20, t50, t80);
				yy = newArray(f20, f50, f80);
				Plot.add("circles", xx, yy);
			}
			txt = "";
			txt += "well=" + well + "_";
			txt += "t20=" + d2s(t20, 2) + "_";
			txt += "t50=" + d2s(t50, 2) + "_";
			txt += "t80=" + d2s(t80, 2) + "_";
			txt += "amp=" + d2s(amp, 1) + "_";
			annotations[well - 1] = txt;
			
			Plot.setLineWidth(1);
			Plot.setXYLabels("Time [min]   well #"  + well, "Brightness");
			Plot.setLimits(NaN,NaN, plotMinY, plotMaxY);			
			Plot.appendToStack;
		}
	}
	Plot.show();
	rename("dummy");
	close("Plots.tif");
	run("Duplicate...", "duplicate title=Plots.tif");//make non-virtual
	close("dummy");
	if(show3Points){
		setFont("SansSerif", 10);
		setColor("magenta");
		for(slc =1; slc < nSlices; slc++){
			setSlice(slc);
			lines = replace(annotations[slc-1], "_", "\n");
			drawString(lines, 100, 35);
		}
		setSlice(1);
	}
	path = basePath;
	plotTitle = "Plots_" + baseName + ".tif";
	plotPath = path + plotTitle;
	if(File.exists(plotPath)){
		showMessageWithCancel("Overwrite Plots on disk?");
	}
	close(plotTitle); //avoid older plots
	saveAs("Tiff", plotPath);

	flags = newArray(0);//rescue old flags
	tablePath = path + tabB;
	if(File.exists(tablePath)){
		open(tablePath);
		if(Table.columnExists("Flags", tabB)){
			flags = Table.getColumn("Flags", tabB);
			if(flags.length != Table.size(tabB))
				exit("update to daily build 1.54j15?");
		}
	}

	close(tabB);
	Table.create("tmp-B");
	for(row=0; row < annotations.length; row++){
		parts = split(annotations[row], "=_");
		for(col = 0; col < parts.length; col+=2){
			colTitle = parts[col];
			data = parts[col + 1];
			Table.set(colTitle, row, data, "tmp-B");
		}
	}
	//Table.update("tmp-B");
	Table.rename("tmp-B", tabB);	
	Table.setColumn("Flags", flags, tabB);
	Table.save(tablePath, tabB);	
}


//Reads Metadata from all stack frames and prepares a table:
//columns: for file, time, well_1 - well_96
//one row per frame 
//Example for timestamp in getMetadata("label") : "PFL 2022-08-30 12h52m56s(Alexa 488)"
function makeTimeTable(){
	if(nSlices<2)
		exit("Stack expected");
	title = getTitle;
	tabA = "TableA_" + replace(title, ".tif", ".csv");
	Table.create(tabA);
	withTimeStamp = true;
	for(frame = 1; frame <= nSlices; frame++){
		setSlice(frame);
		info = getMetadata("label");
		info = replace(info, ".jpg", "");
		ss =split(info, " -hms()");//separating hours, minutes, seconds etc
		if(ss.length < 7){
			withTimeStamp = false;
			break;
		}
		year = parseInt(ss[1]);
		if(!(year > 1990) || !(year < 2100)){
			withTimeStamp = false;
			break;
		}		
		day = parseInt(ss[3]);
		hours = parseInt(ss[4]);
		if(frame == 1)
			prevDay = day;
		if(day != prevDay){
			hours += 24;//survives midnights, intervals must be <24h
			prevDay = day;
		}
		minutes = parseInt(ss[5]);
		seconds = parseInt(ss[6]);
		gMinutes = hours * 60 + minutes + seconds/60;
		if(frame == 1 )
			startMinutes = gMinutes;
		Table.set("frame", frame-1, frame, tabA);
		Table.set("File", frame-1, info, tabA);
		Table.set("Thr", frame-1, thresh, tabA);
		Table.set("Back", frame-1, background, tabA);
		Table.set("Time[min]", frame-1, gMinutes - startMinutes, tabA);
		if(frame > 1){
			interval = gMinutes - prevMinutes;
			Table.set("Interval[min]", frame-1, interval, tabA);
			if(interval > 2)
				Table.set("Interaction", frame-1, "interval > 2", tabA);	
		}
		prevMinutes = gMinutes;
	}
	if(!withTimeStamp){
		interval = getNumber("No valid time stamps: enter frame interval [min]", 1);
		for(frame = 1; frame <= nSlices; frame++){
			Table.set("frame", frame-1, frame, tabA);
			Table.set("File", frame-1, "", tabA);
			Table.set("Thr", frame-1, thresh, tabA);
			Table.set("Back", frame-1, background, tabA);
			Table.set("Time[min]", frame-1, (frame - 1) *interval, tabA);
		}
	}
	Table.update(tabA);
	setSlice(1);
}

		
macro "Toggle Overlay [6]" {
	run("Toggle Overlay");
	Overlay.useNamesAsLabels(false);
	Overlay.drawLabels(true);
}

macro "Toggle Labels/Flags [7]" {
	checkBaseName();
	if(labelFlag)
		roiManager("Show None");
	else{
		if(!isOpen(tabB))
				showTable(tabB);
		if(roiManager("count") == 0)
			updateFlagRois();
		roiManager("Show All");
	}
	Overlay.drawLabels(labelFlag);
	labelFlag = !labelFlag;
}

//Duplicates stack and sets pixel outside circles and pixels<thresh
//to zero and applies greenLUT. This stack is intended only for visualisation
macro "Show Isolated Leaves"{
	checkBaseName();
	selectImage(baseName + ".tif");
	if(Overlay.size < 96)
		exit("Time-lapse stack with 96 overlay circles should be in front");
	Overlay.show;//09.04.24, 11:34
	if(isOpen(tabA))
		thr2 = Table.get("Thr", 0, tabA);
	else
		thr2 = getNumber("Threshold (see TableA): ", thresh);
	close("IsolatedLeaves");
	run("Duplicate...", "ignore title=IsolatedLeaves duplicate");
	roiManager("reset");
	for(jj = 0; jj < 96; jj++){
		Overlay.activateSelection(jj);
		roiManager("Add");
	}
	roiManager("combine");
	roiManager("reset");
	roiManager("Add");
	setBackgroundColor(0, 0, 0);
	run("Clear Outside", "stack");
	Overlay.remove;
	run("Select None");
	run("Green");
	for(slc = 1; slc <= nSlices; slc++){
		setSlice(slc);
		changeValues(0, thr2-1, 0);
	}
	setSlice(1);
	roiManager("reset");
	setOption("Changes", false);
}


//Creates from a stack of 96 plots eg. a stack of 12 combined plots
function projectPlots(){
	close("MIN_Plots*");
	wellsPerPlot = 8;//number of projections per output slice
	if(nSlices != 96)
		exit("Stack of 96 plots expected");
	srcID = getImageID;
	setBatchMode(true);//11.04.24, 15:52
	for(start = 1; start <= 96; start += wellsPerPlot){
		stop= start + wellsPerPlot -1;
		run("Z Project...", "start=&start stop=&stop projection=[Min Intensity]");
		if(start ==1){
			destID = getImageID;
		}
		else{
			run("Select All");
			run("Copy");
			close;
			selectImage(destID);
			run("Add Slice");
			run("Paste");
			run("Select None");
		}
		setMetadata("label", "wells "+ start + "-" + stop);
		selectImage(srcID);
	}
	setBatchMode("exit and display");
	selectImage(destID);
	setBatchMode("show");
	setSlice(1);
	setOption("Changes", false);
}

function dynamicPlots(){
	plotColors = split("blue red magenta darkGray orange green black ");
	showTable(tabA);
	qualified = getString("Leaf numbers to plot:", "3 6 14");
	numbers = split(qualified, " ,");
	plotTitle = "Plot Leaves #" + qualified;
	plotTitle = plotTitle.replace(" ", "_");
	close(plotTitle);
	Plot.create(plotTitle, "Time [min]", "Brightness");
	Plot.setFrameSize(350, 300);
	xx1 = Table.getColumn("Time[min]", tabA);
	cTitles = split(Table.headings);
	Plot.setLineWidth(2);
	legend = "";
	for(n = 0; n < numbers.length; n++){
		for(col = 0; col < cTitles.length; col++){
			if(cTitles[col] == "Mean_" + numbers[n]){
				Plot.setColor(plotColors[n], plotColors[n]);
				yy1 = Table.getColumn(cTitles[col], tabA);
				Plot.add("connected", xx1, yy1, "Leaf#" + numbers[n]);
				legend += "leaf " + numbers[n] + "\t";
			}
		}
	}
	Plot.setLegend(legend);
	Plot.setLimitsToFit();
	Plot.setLineWidth(1);
	Plot.show;
}


macro "Show Plots [p]"{
	if(baseName == ""){
		checkBaseName();
		exit("Try again");
	}
	checkBaseName();
	plotsName = "Plots_" + baseName + ".tif";
	if(!isOpen(plotsName))
		open(basePath + plotsName);
	selectWindow(plotsName);
	
	Dialog.create("Show Plots");
	Dialog.addMessage("Additionally show:");
	Dialog.addCheckbox("Montage 12x8", false);
	Dialog.addCheckbox("Projections", false);
	Dialog.addCheckbox("Dynamic Plots", false);
	Dialog.show;
	mont = Dialog.getCheckbox;
	proj = Dialog.getCheckbox;
	dyn = Dialog.getCheckbox;
	if(proj)
 		projectPlots();
 	if(mont){
		selectWindow(plotsName);
 		close("Montage");
			run("Make Montage...", "columns=12 rows=8 scale=0.3 label");
	}
 	if(dyn)
		dynamicPlots();
}

function showTable(theTab){
	if(baseName == ""){
		checkBaseName();
		exit("Try again");
	}
	checkBaseName();
	if(!isOpen(theTab))
		open(basePath + theTab);
	selectWindow(theTab);
}

macro "Set Flag A [8]"{
	setFlag("A");
}
macro "Set Flag B [9]"{
	setFlag("B");
}
macro "Set Flag C [0]"{
	setFlag("C");
}


function setFlag(flag){
	cmd = getInfo("command.name");
	if (!startsWith(cmd, "^"))
		exit("Command must be invoked via shortcut key, \nwhile cursor is at desired position");
	if(baseName == ""){
		checkBaseName();
		exit("Try again");
	}
	checkBaseName();
	if(isKeyDown("alt")){
		flag = getString("Set or Clear Flag", flag);
	}
	if(!isOpen(tabB))
		open(basePath + tabB);
	selectImage(baseName + ".tif");
	getCursorLoc(x0, y0, z, flags);
	showStatus("x="+x0);
	well = getClosestObject(x0, y0);//0-based
	if(well == -1)
		exit("Cursor not close to a well");

	if(!Table.columnExists("Flags", tabB))
		Table.setColumn("Flags", newArray(0), tabB); 
	oldFlag = Table.getString("Flags", well, tabB);

	if(isKeyDown("alt")){
		newFlag = getString("Change Flags:", oldFlag);
	}
	else {
		if(oldFlag.contains(flag))
			newFlag = oldFlag.replace(flag, "");
		else
			newFlag = oldFlag + " " + flag + " ";
	}
	newFlag = sortFlags(newFlag);
	Table.set("Flags", well, newFlag, tabB);
	Table.update;//Michael
 	Table.setSelection(well, well, tabB);
	updateFlagRois();
	Table.save(basePath + tabB, tabB);
	selectImage(getImageID);//
	//labelFlag = true;
}

 function updateFlagRois(){
 	selectImage(baseName + ".tif");//11.04.24, 17:04
  	setFont("", 18);
 	roiManager("reset");
	for(well = 0; well < 96; well++){//zero based
		flag = Table.getString("Flags", well, tabB);
		if(flag != ""){
			Overlay.getBounds(well, x, y, width, height);
			makeText(flag, x, y-20);
			Roi.setName("" + (well + 1));
			Roi.setStrokeColor(flagColor);
			roiManager("add");
		}
	}
	run("Select None");
	roiManager("Show All");
 }

//Returns index of closest circle (zero-based)
function getClosestObject(x0, y0){
	found = -1;
	deltaMin = 1e6;
	for(well= 0; well < 96; well++){
		Overlay.getBounds(well, x1, y1, w, h);
		delta = Math.sqrt(Math.sqr(x1 +w/2 - x0) + Math.sqr(y1 + h/2 - y0));
		if (delta < deltaMin){
			found = well;
			deltaMin = delta;
		}
	}
	if(deltaMin > 70)
		return -1;
	return found;
}

macro "Edit Flags in TableB [e]"{
	checkBaseName();
	showTable(tabB);
	waitForUser("Change Flags", "Optionally, select range of wells in TableB, \nthen click OK");
	wellsRange = "";
	start = Table.getSelectionStart(tabB);
	stop = Table.getSelectionEnd(tabB);
	if(start >=0)
		wellsRange =  "" + (start + 1) + "-" + (stop + 1);
	wait(50);
	Dialog.create("Change Flags");
	Dialog.addString("Wells range:", wellsRange);
	Dialog.addString("Flags to clear *", "");
	Dialog.addString("Flags to set", "");
	Dialog.addCheckbox("Auto-save", true);
	msg = "Range is preset by table selection";
	msg = "* Enter 'all' to clear existing flags";
	Dialog.addMessage(msg, 9, "darkgray"); 
	Dialog.show();
	range = Dialog.getString();
	toDelete = Dialog.getString();
	toAdd = Dialog.getString();
	if(toAdd.contains("all"))
		exit("'all' not allowed here");
	forgetOld = toDelete.contains("all");
	autoSave = Dialog.getCheckbox();

	num = newArray(1, 96);
	num = Array.resample(num, 96);
	mask = newArray(96);
	if(range != ""){
		ranges = split(range, ",,");
		for(kk = 0; kk < ranges.length; kk++){
			startStop = split(ranges[kk], "--");
			start = parseInt(startStop[0])-1;
			stop = start;
			if(startStop.length == 2)
				stop = parseInt(startStop[1])-1;
			for(index = start; index <= stop; index++)
				mask[index] = 1;
		}
	}
	for(index = 0; index < 96; index++){
		if(mask[index] ==1){
			flags = Table.getString("Flags", index, tabB);
			if(forgetOld)
				flags = "";
			else  for(jj = 0; jj < toDelete.length; jj++){
				chr = substring(toDelete, jj, jj+1);
				flags = replace(flags, chr, "");
			}
			for(jj = 0; jj < toAdd.length; jj++){
				chr = substring(toAdd, jj, jj+1);
				if(!flags.contains(chr))
					flags += chr;
			}
			flags = sortFlags(flags);
			Table.set("Flags", index, flags, tabB);	
		}
	}
	Table.update(tabB);
	Table.setSelection(start, stop, tabB);
	if(autoSave)
		Table.save(basePath + tabB, tabB);//11.04.24, 1:12
	updateFlagRois();
}

function sortFlags(flags){
	chars = newArray(0);
	flags = replace(flags, " ", "");
	chars = newArray(0);
	for(jj = 0; jj < flags.length; jj++)
		chars[jj] = substring(flags, jj, jj+1);
	Array.sort(chars);
	flags = "";
	for(jj = 0; jj < chars.length; jj++)
		flags = flags + " " + chars[jj];
	return flags;
}

function checkBaseName(){
	requires("1.54j");//avoid NaNs in table
	if(isOpen(baseName + ".tif"))
		return;
	if(nImages == 0)
		exit("First open Plate image");
	good =  (nSlices > 2) && (Overlay.size >= 96);
	title = getTitle;
	if(!endsWith(title, ".tif"))
		exit("Expecting Plate stack ending with '.tif'");
	if(substring(title, title.length -6, title.length -5) == "-")
		exit("Avoid minus sign at the end of file name:\n" + title);
	if(File.exists(basePath + title))
		exit("Not found on disk: \n" + basePath + title);
	if(good){
		showMessageWithCancel("Register Plate", "Register file:\n" + getTitle +"\nas current plate?");
		baseName = File.getNameWithoutExtension(getTitle);
		basePath = getDir("image");
		tabA = "TableA_" + baseName + ".csv";
		tabB = "TableB_" + baseName + ".csv";
	}
	else 
		exit("Expecting stack with 96 overlay circles in front"); 
}

macro "Show Related Files [f]"{
	checkBaseName();
	Dialog.create("Show file:");
	Dialog.addMessage("BaseName = " + baseName);
	
  	items = split("Plate,TableA,TableB,Show All,Show in folder", ",");
  	Dialog.addRadioButtonGroup("Show linked file:", items,  10, 1, "Plate");
	Dialog.show;
	choice = Dialog.getRadioButton;
	all = (choice == "Show All");
	if(choice == "Plate" || all)
		selectImage(baseName + ".tif");
	if(choice == "TableA"|| all)
		showTable(tabA);
	if(choice == "TableB"|| all)
		showTable(tabB);
	if(choice == "Show in folder"){
		selectImage(baseName + ".tif");
		run("Image"); 
	}
}