/*
SporeTrackerB06f
14.11.19, 22:36
*/

var	radP = 5,//roi radius for phase-bright detection, default = 5 px
	sigmaForPeak = 1.5,// gaussian blur for peak detection, default = 1.5 px
	sigmaForView = 0.8,// gaussian blur for viewing only, default = 0.8 px
	intFlRad = 5, //radius for integrated fluorescence of spots, default = 5 pixels
	medianKernelForBackgnd = 7, //* radius for median kernel for cell background calculation, default=10 px
	medianKernelForView = medianKernelForBackgnd, //for viewing only
	thrFraction = 0.50, //set threshold to fraction of min-max range for particle analyzer (default=0.50)
	sizeLimit = "size=[0.2 -0.7]", //[um] particle Analyzer, default = "size=[0.15 -0.7]"
	minNeighborDistance = 0.8, //contour to contour, default = 1 um
	tolerance = 20, //noise tolerance for maximum finder, default = 20
	binWidth = 10, //default = 10 minutes (histo germination)
	expandContour = 0,//if > 0,  enlarge phase-bright contour and set outside area to img-backgnd, default = 10 px
	boxW = 44, // * boxWidth for montage
	plotLimitTop = 35,//y plotrange Fluorescence at leader position vs time
	minPeakForMap = 7,//weaker leaders are not included in Map (default=10)
	stepConfidence = 1, //density of error bars, 0=omit them
	traceRad = 0,//neighbor searching for maxima, default = 0 px =disabled

	//-- end user selectable parameters

	boxR = boxW/2, //boxRadius
	viewWhat = "raw",//when browsing through fluor-montage
	k4 = 4, //matrix rows
	k9 = 9, //matrix columns (k9*k4=36 frames)
	k36 = k4*k9,
	nTime = 0,//avoid cascading shortcuts
	auxPath = "",
	acronym = "",
	withBatch = true, //set false for testing
	rawGermX, germY,
	currentTitle = "",
	dummy = doQueue(),
	;

//Analyse particles, mark object centers in slice 1, and 
//store value pairs : phC & time in column "Phase" througout the movie 
macro "Mark Spores"{ 
	if (ojNObjects() > 0)
		showMessageWithCancel("Delete " + ojNObjects() + " objects?");
	markSpores();
	showStatus("Finished marking spores, n= " + ojNObjects());
	ojShowImage(1);
}
function markSpores(){
	requires("1.52p");
	close("Montage*");
	close("Map*");
	close("boxStack*");
	close("*_Plot*");
	close("Fluor_at_LeaderPos.tif");
	
	ojHideResults();
	ojSetVirtualFlag(true);//we don't change pixels
	ojDeleteColumn("*");
	ojDeleteAllObjects();
	ojInitTextColumn("Phase");//for arrays as text
	ojInitTextColumn("Spots");//for arrays as text

	ojInitColumn("Thr");
	ojInitColumn("PrepTime");
	ojInitColumn("RawGTime");//germination time excluding prepTime
	ojInitColumn("GermTime");//germination time including prepTime
	ojInitColumn("RawMaxFlTime");//time of max spot excluding prepTime
	ojInitColumn("MaxFlTime");//time of max spot including prepTime
	ojInitColumn("PeakToGermTime");//time of max spot including prepTime
	ojInitColumn("Neighbor");
	ojInitColumn("MaxFlA");//Integrated Fluoescence of leader spot
	ojInitColumn("FlT0");//Fluorescence at leader pos at t = 0 minutes
	ojSetColumnProperty("*", "width", 100);
	ojSetColumnProperty("MaxFlA", "color", 0x0000ff);
	ojSetColumnProperty("Phase FLuor Spots Neighbor Thr PrepTime", "visible", 0);
	run("Set Measurements...", "area mean centroid fit shape decimal=2");

	if(ojNImages() == 0)
		exit("No images are linked");
	for(img = 1; img <= ojNImages(); img++){
		if (ojGetImageValue(img, "frameinterval") > 100)
			exit("Img " + img + " frameInterval=" + ojGetImageValue(img, "frameinterval"));
		if (ojGetVoxelSize(img, "x") > 0.5)
			exit("Img " + img + " pxWidth> 0.5 ");//microns
	}
	for(img = 1; img <= ojNImages(); img++){
		showProgress(-img/ojNImages());
		frameInterval = ojGetImageValue(img, "frameinterval");//must be in minutes
		setBatchMode(false);
		ojShowImage(img);
		run("Enhance Contrast", "saturated=0.35");
		run("Set... ", "zoom=100");
		prepTime = 0;
		key= "preparationtime=";
		info = toLowerCase(getMetadata("Info") + "\n");
		a = indexOf(info, key);
		if(a>=0){
			b = indexOf(info, "\n", a);
			prepTime = parseFloat(substring(info, a+lengthOf(key), b));
		}		
		
		setSlice(1);
		run("Remove Overlay");
				
		run("Select All");
		run("Median...", "radius=8 slice");
		getStatistics(area, mean, min, max);
		run("Revert");
		run("Median...", "radius=1.5 slice");

		lowerThr = min + thrFraction *(max - min); 
		upperThr = 1e9;
		setThreshold(lowerThr, upperThr);
		if(is("Caps Lock Down")){
			run("Threshold...");
			waitForUser("Caps Lock is Down: \n \nOptionally adjust threshold, \nthen click Ok");;
			getThreshold(lowerThr, upperThr);
			close("Threshold");
		}
		setThreshold(lowerThr, upperThr);
		createVoronoi();
		run("Select All");
		setThreshold(lowerThr, upperThr);
		run("Analyze Particles...", sizeLimit + "  exclude add clear include slice");
 		pxSize = ojGetVoxelSize(img, "x");
		dNeighbor = newArray(nResults);
		selectImage("Voronoi");
		for(jj = 0; jj <nResults; jj++){
			roiManager("select", jj);
			getStatistics(area, mean, min, max);
			dNeighbor[jj] = min;//half neighbor distance
		}
		close("Voronoi");
		ojShowImage(img);		
		for(jj = 0; jj <nResults; jj++){
			x = round(getResult("X", jj)/pxSize);
			y = round(getResult("Y", jj)/pxSize);
			inside = x > boxR && y >boxR && x < getWidth - boxR && y < getHeight - boxR;
			solidity = getResult("Solidity", jj);
			if(inside && solidity > 0.9){ 
				ojSetMarker(x, y);
				ojCloseObject();
				ojSetResult("Thr", ojNewestObject(), lowerThr);
				ojSetResult("Neighbor", ojNewestObject(), 2 * dNeighbor[jj] * pxSize);
				ojSetResult("PrepTime", ojNewestObject(), prepTime);
			}
		}
		ojSetTool(0);
		roiManager("reset");
		resetThreshold;
		if(withBatch) setBatchMode("hide");
		for (obj = ojFirstObject(img); obj <= ojLastObject(img); obj++){
			showProgress(obj/(ojNObjects()+50));
			ojSelectObject(obj);
			makeOval(ojXPos(1) - radP, ojYPos(1)-radP, radP*2, radP*2);
			Stack.getDimensions(ww, hh, channels, slices, frames);
			phaseStr = "phase ";
			timeStr = "time ";
			for(frame = 1; frame <= frames; frame++){
				time = (frame-1) * frameInterval;
				Stack.setPosition(1, 1, frame);
				run("Median...", "radius=3 slice");//tmp median filter before max calculation
				List.setMeasurements;
				max = List.getValue("Max");
				timeStr += d2s(time, 2) + " ";
				phaseStr += d2s(max, -4) + " ";
			}
			
			ojSetResult("Phase", obj, timeStr + ";" + phaseStr + ";");
		}
		close;
	}
	setBatchMode(false);
	ojShowImage(1);
	ojSetTool(0);
	close("Threshold");
}

macro "Delete Neighbors"{deleteNeighbors();}
function deleteNeighbors(){
	n1 = ojNObjects();
	for(obj = n1; obj >=1; obj--){
		if(ojResult("Neighbor", obj) < minNeighborDistance) // default 0.8 um
			ojDeleteObject(obj);
	}
	n2 = ojNObjects();
	selectWindow("ImageJ");
	showStatus("Deleted Neighbors: " + (n1 - n2) +  ", Remaining=" + round(n2/n1 * 100) + "%");
}


//Calculate germination point per object
//create phC plot family with marked germination point
//put arrays as strings into column "Phase"
macro "Germination Plots"{germinationPlots(true);}
function germinationPlots(ask){
	if(ojNObjects() == 0)
		exit("No objects marked yet");
	sync = true;
	if(ask)
		sync = getBoolean("Synchronize germination times?");
	ojHideResults();
	
	close("Phase_Plots.tif");
	if (sync) xStr = "Time relative to germiantion [min]";
	else xStr = "Time incl. preparation [min]";
	Plot.create("Phase_Plots.tif", xStr, "PhC_Brightness");
	
	Plot.setFrameSize(500, 200);
	labels = "";
	for(obj = 1; obj <= ojNObjects(); obj++){
		showStatus(-obj/ojNObjects()/2);
		xValues = extractArray("Phase", obj, "time");
		yValues = extractArray("Phase", obj, "phase");
		calcGerminationPoint(xValues, yValues);
		ojSetResult("RawGTime", obj, rawGermX);
		pTime = ojResult("PrepTime", obj);
		gTime = rawGermX + pTime;
		ojSetResult("GermTime", obj, gTime);

		if(obj == 1){
			Array.getStatistics(yValues, min, max, mean, stdDev);
			limitTop = max * 1.5; 
			Plot.setLimits(0, 180, 0, limitTop );
			if(sync)
				Plot.setLimits(-180, 180, 0, limitTop );
		}
		if(sync){
			for(jj = 0; jj < xValues.length; jj++)
				xValues[jj] -= rawGermX;
		}
		
		else{
			for(jj = 0; jj < xValues.length; jj++)
				xValues[jj] += pTime;		
			if(rawGermX > 0){//draw red bar first
				Plot.setColor("red");
				Plot.setLineWidth(3);
				Plot.drawLine(gTime-2, germY, gTime + 2, germY);
			}
		}
		Plot.setColor("blue");
		
		Plot.setLineWidth(1);
		if(!sync)
			Plot.add("connected", xValues, yValues);
		if(sync && rawGermX > 0)
			Plot.add("line", xValues, yValues);
		Plot.setColor("black");
		
		Plot.appendToStack;
		labels = labels + " " + obj;
		
	}
	Plot.show;
	labels = split(labels);
	for(slc = 1; slc <= ojNObjects(); slc++){
		setSlice(slc);
		setMetadata("label", "obj=" + labels[slc - 1]);
	}
	addProjection("Phase_Plots.tif");
	setLocation(720, 100);
	selectWindow("ImageJ");
	showStatus("--- Finished with Germination Plots");
}

//add superimposition of all plots to last slice
function addProjection(newTitle){
	rename("tmp");
	run("Duplicate...", "duplicate");
	rename(newTitle);
	close("tmp");

	run("Z Project...", "projection=[Min Intensity]");
	run("Select All");
	run("Copy");
	close;
	setSlice(nSlices);
	run("Add Slice");
	
	run("Paste");
	run("Select None");
	setOption("changes", false);
	selectImage(newTitle);
}
 

//Example for a string with two named arrays: "fluor 22 33 44;time 0 5 10;"
function extractArray(colTitle, obj, arrTitle){
	ss = ojResultString(colTitle, obj);
	lines = split(ss, ";;");
	for(jj = 0; jj<lines.length; jj++){
		if(startsWith(lines[jj], arrTitle + " ")){
			parts = split(lines[jj]);
			arr = newArray(parts.length-1);
			for(kk = 1; kk < parts.length; kk++)
				arr[kk-1] = parseFloat(parts[kk]);
			return arr;
		}
	}
	return newArray(0);
}

//For detecting closest distance between contours
function createVoronoi(){
	run("Select None");
	id = getImageID;
	run("Duplicate...", "title=Voronoi");
	setThreshold(lowerThr, upperThr);

	setOption("BlackBackground", true);
	run("Make Binary");
	run("Erode");
	run("Dilate");
	run("Voronoi");
	run("Multiply...", "value=255");
	run("Invert");
	if(ojNObjects() == 0){
		run("Restore Selection");
		wait(100);
	}
	run("Distance Map");
	if(ojNObjects() == 0){
		run("Enhance Contrast", "saturated=0.35");
		run("Restore Selection");	
		wait(100);
	}
	selectImage(id);
	wait(10);
}

function calcGerminationPoint(xValues, yValues){
	len = xValues.length * 10;
	xx = Array.resample(xValues, len);
	yy = Array.resample(yValues, len);

	Array.getStatistics(yValues, min, max, mean, stdDev);
	if(min > 0.5 * max){//no prominent drop
		rawGermX = NaN;
		return;
	}
	for(jj = 0; jj < len; jj++){
		if(yy[jj] < (min + max)/2){
			rawGermX = xx[jj];
			germY = yy[jj];
			return;
		}
	}
}



macro "Make Montages"{makeMontages();}
//Create phase and fluor montage and add triangles overlay to phase montage
function makeMontages(){
	if(ojNObjects() == 0)
		exit("No objects marked yet");
 	getAuxPath();
	dummy = File.makeDirectory(auxPath);//unless already exists
	fluorPath = auxPath +"/Montage_Fluor.tif";

	boxTitles = newArray("", "boxStackPhase", "boxStackFluor");
	montageTitles  = newArray("", "Montage_Phase.tif", "Montage_Fluor.tif");
	ojHideResults();
	close("Montage*");
	roiManager("reset");
	roiManager("Associate", "true");
	if(withBatch)  setBatchMode(true);
	for(img = 1; img <= ojNImages(); img++){
		ojShowImage(img);
		srcID = getImageID;
		setSlice(2);
		run("Select None");
		run("Smooth", "slice");
		medianFr1 = getValue("Median");//median of fluor frame 1 used as background
		run("Undo");
		
		if(withBatch) setBatchMode("hide"); 
		frameInterval = ojGetImageValue(img, "frameinterval");
		for(obj = ojFirstObject(img); obj <= ojLastObject(img); obj++){
			showProgress(-obj/ojNObjects());
			ojShowObject(obj);
			x = ojXPos(1);
			y = ojYPos(1);
			close("boxStack*");
			makeRectangle(x-boxR, y- boxR, boxR*2, boxR*2);
			run("Duplicate...", "title=boxStackFluor duplicate channels=2");
			selectImage(srcID);
			run("Duplicate...", "title=boxStackPhase duplicate channels=1");
			lowerThr = ojResult("Thr", obj);

			if(expandContour > 0){//if > 0,  enlarge phase-bright contour and set outside area to img-backgnd
				setThreshold(lowerThr, 1e9);
				doWand(boxR,boxR);		
				selectImage("boxStackFluor");
				run("Restore Selection");
				run("Enlarge...", "enlarge=&expandContour pixel");//default 10 px
				run("Make Inverse");
				run("Set...", "value=&medianFr1 stack");//set surrounding to image background
			}
			run("Select None");

			for(chn = 1; chn <=2; chn++){
				selectImage(boxTitles[chn]);
				run("Make Montage...", "columns=&k9 rows=&k4 scale=1 border=0");
				if(obj == 1)
						rename(montageTitles[chn]);
				else {
					run("Select All");
					run("Copy");
					selectImage(montageTitles[chn]);
					run("Add Slice");
					run("Paste");
					close("Montage");
				}
				setMetadata("label", "obj=" + obj);
			}
		}	
	}
	close("boxStack*");
	getAuxPath();
	for(chn = 1; chn <=2; chn++){
		selectImage(montageTitles[chn]);
		for(obj = 1; obj <= ojNObjects(); obj++){
			rawGermT = ojResult("RawGTime", obj);
			setSlice(obj);
			addGermTriangle(rawGermT);
		}
		if(chn==1)
			run("Grays");
		if(chn==2)
			run("Green");
		run("Select None");
		fullPath = auxPath + "/" + montageTitles[chn];
		run("Enhance Contrast", "saturated=0.35");
		//run("Labels...", "color=white font=12");//remove labels
		saveAs("Tiff", fullPath);
		setBatchMode("show");			
	}
	selectWindow("ImageJ");
	showStatus("--- Finished Make Montages");

}	



//adds germination triangle to overlay
function addGermTriangle(rawGermT){
	if(!isNaN(rawGermT)){
		len = rawGermT/frameInterval * boxR*2 + boxR;
		x = len % getWidth;
		y = floor( len/getWidth) * boxR * 2;
		setColor(255,0,0);
		makePolygon(x-3, y, x+3, y, x, y+3);
		run("Properties... ", "  fill=magenta");
		run("Add Selection...");
		run("Labels...", "color=white font=12");//remove labels
	}
}

function doQueue() {
	setOption("QueueMacros", true);
	return true;
}


function getAuxPath() {
	prName = ojGetProjectName();
	dir = ojGetProjectPath();
	auxPath =  dir+prName + "Aux";
	if(ojNImages() > 0){
		if(startsWith(ojGetImageName(1), "mcc"))
			acronym = "mcc";
		if(startsWith(ojGetImageName(1), "rcc"))
			acronym = "rcc";
	}
  }

macro "Load Montages"{loadMontages();}

function loadMontages(){
	titlesTif = split("Montage_Phase.tif Montage_Fluor.tif");
 	getAuxPath();
	msg = "loaded: ";
	for(img = 0; img <titlesTif.length; img++){
		fullPath = auxPath + "/" +titlesTif[img];
		if(File.exists(fullPath)){
			close(titlesTif[img]);
			run("TIFF Virtual Stack...", "open=[" + fullPath +"]");
			msg += titlesTif[img] + ",  ";
			showStatus(msg);
		}
	}
	arrangeWindows();
	navigate(0);
}

macro "Mark Fluor. Peaks"{markPeaks();}

function markPeaks(){
	if(ojNObjects() == 0)
		exit("No objects marked yet");
	getAuxPath();
	ojHideResults();
	ojInitColumn("MaxFlA");

	fluorPath = auxPath + "/Montage_Fluor.tif";
	phasePath = auxPath + "/Montage_Phase.tif";
	close("Montage_*");
	run("TIFF Virtual Stack...", "open=[" + phasePath + "]");
	run("TIFF Virtual Stack...", "open=[" + fluorPath + "]");
	if(nSlices != ojNObjects())
		exit("Error: objects ("+ ojNObjects() + ") do not match  montage slices ("+ nSlices + ")");
	run("Select None");
	Overlay.clear;
	roiManager("reset");//will later accumulate all circles in all slices
	if(withBatch) 
		setBatchMode("hide");
	run("Enhance Contrast", "saturated=0.35");
	for(obj = 1; obj <= ojNObjects(); obj++){
		showProgress(-obj/ojNObjects());
		ojSelectObject(obj);
		setSlice(obj);
		markObjectSpots(obj);	
	}
	close("Montage_Fluor.tif");
	open(fluorPath);//open normally
	
	Overlay.clear;
	run("From ROI Manager");//circles
	selectImage("Montage_Phase.tif");
	setBatchMode("show");
	run("Labels...", "color=white font=12");//remove labels
	if (Overlay.size > 0)
		run("To ROI Manager");//triangles ___
	selectImage("Montage_Fluor.tif");
	run("From ROI Manager");//add triangles to fluor
	//run("Labels...", "color=white font=12");//remove labels
	run("Select None");
	run("Enhance Contrast", "saturated=0.35");
	saveAs(fluorPath);
	
	close("Montage*");	//close and reopen virtually
	setBatchMode(false);
	run("TIFF Virtual Stack...", "open=[" + fluorPath + "]");
	run("TIFF Virtual Stack...", "open=[" + phasePath + "]");
	ojSelectObject(1);
	navigate(0);
	arrangeWindows();
	selectWindow("ImageJ");
	showStatus("--- Finished Marking Peaks");

}

/*we perform gauss smoothing (default 1.5 px) 
- create a "Back" image (cell background after heavy median filter)
- subtract Back from smoothed image
- find "leader" = stongest peak position (Find Maxima)
- restore raw image
- sample mean of circular roi at leader position in all 36 frames => totalFl
- use same rois in "Back" image to record cell backgrounds
- 36 times, 36 totalFl and 36 cellBack are squeezed into single string
*/

function markObjectSpots(obj){
	ojHideResults();
	selectImage("Montage_Fluor.tif");
	run("Revert");
	frameInterval = ojGetImageValue(1, "frameinterval");
	prepTime = ojResult("PrepTime", obj);
	rawGermT = ojResult("rawGTime", obj);
  	intFls = "intFl ";
 	times = "time ";
	
 	run("Gaussian Blur...", "sigma=&sigmaForPeak slice");

	close("Back");
	run("Select None");
	run("Duplicate...", "title=Back");
	
	run("Median...", "radius=&medianKernelForBackgnd");
	selectImage("Montage_Fluor.tif");
	imageCalculator("Subtract", "Montage_Fluor.tif","Back");	
	
	run("Find Maxima...", "noise=&tolerance exclude output=[Point Selection]");
	if( selectionType  != 10)
		return;
	getSelectionCoordinates(xx, yy);
	
 	xL = xx[0] % boxW - boxR;//local from center
 	yL = yy[0] % boxW - boxR; 

	netFlArr = newArray(k36);
	roisStart = roiManager("count");
	for(frame = 0; frame < k36; frame++){
		col = frame%k9;
		row = floor(frame/k9);
		hookX = col * boxW + boxR + xL - intFlRad;
		hookY = row * boxW + boxR + yL - intFlRad;
		
		foundDx = 0;
		foundDy = 0;
		if (traceRad >0){//neighbor searching for highest fluorescence
			foundMean = 0;
			for(dy = -traceRad; dy <traceRad; dy++){
				for(dx = -traceRad; dx <traceRad; dx++){
					if(sqrt(dx * dx + dy *dy) < traceRad){
						makeOval(hookX + dx, hookY + dy,  intFlRad*2,  intFlRad*2);
						getRawStatistics(area, mean, min, max);
						if(mean > foundMean){
							foundMean = mean;
							foundDx = dx;
							foundDy = dy;
						}
					}
				}
			}
		}
		makeOval(hookX + foundDx, hookY + foundDy,  intFlRad*2,  intFlRad*2);
		getRawStatistics(area, mean, min, max);


		
		
		
		//makeOval(hookX, hookY,  intFlRad*2,  intFlRad*2);
		run("Properties... ", "  size=Small stroke=magenta");
		roiManager("add");
		getRawStatistics(area, mean, min, max);
		factor = intFlRad * intFlRad / 36;// =1 in our case
		mean*= factor;
		netFlArr[frame] = mean;
		time = frame * frameInterval;
		intFls = intFls + d2s(mean, 1) + " ";	
		times = times + time + " ";
	}
	ranks = Array.rankPositions(netFlArr);
	leaderFrame = ranks[ranks.length - 1];
	rLeaderTime = leaderFrame * frameInterval;//raw
	ojSetResult("FlT0",obj, netFlArr[0]);
	ojSetResult("MaxFlA",obj, netFlArr[leaderFrame]);
	ojSetResult("RawMaxFlTime",obj, rLeaderTime);
	ojSetResult("MaxFlTime",obj, rLeaderTime + prepTime);
	ojSetResult("PeakToGermTime",obj, rawGermT - rLeaderTime);
	roiManager("select", roisStart + leaderFrame);
	run("Properties... ", "  size=Small stroke=red  width=1.5");
	
	
//may be different from xx[0] after background subtraction
	hookX = (leaderFrame% k9) * boxW;
	hookY = floor(leaderFrame / k9) * boxW;
	xG =  hookX + boxR + xL;//global in matrix
	yG = hookY + boxR + yL;//global in matrix

	xGStr = "xG " + xG + ";";
	yGStr = "yG " + yG + ";";
	xLStr = "xL " + xL + ";";
	yLStr = "yL " + yL + ";";

	objArrays = times + ";"  + intFls +";" + xGStr + ";" + yGStr   + ";" + xLStr + ";" + yLStr   + ";";
	ojSetResult("Spots", obj, objArrays);
	run("Select None");
	selectImage("Montage_Fluor.tif");
	run("Select None");
	close("Back");
	roiManager("Show All");
}	

function capsdebug(){
	if(is("Caps Lock Set")){
		withBatch = false;
		debug;
	}
	
}
function arrangeWindows(){
	titles = split("Montage_Phase.tif Montage_Fluor.tif");
	locations = newArray(10, 100, 10, 500);
	for(img = 0; img < titles.length; img++){
		if(isOpen(titles[img])){
			selectImage(titles[img]);
			setLocation(locations[img*2], locations[img*2+1]);
			run("Set... ", "zoom=199.9");//force refresh
			run("Set... ", "zoom=200");
		}
	}	
}	

macro "-"{}
macro "Batch (all above)"{
	msg =  "These tasks will be executed:";
	msg =  "\n---";
	msg =  "\nmarkSpores";
	msg =  "markSpores";
	msg += "\ndeleteNeighbors";
	msg += "\ngerminationPlots";
	msg += "\nmakeMontages";
	msg += "\nmarkPeaks";
	if (ojNObjects() > 0)
		msg += "\n---\n " + ojNObjects() + "  objects will be deleted.";
	showMessageWithCancel("Process sequence:", msg);
	close("*_Peak*");
	close("*Histo_*");
	time2 = getTime;
	markSpores();
	deleteNeighbors();
	germinationPlots(false);
	makeMontages();
	markPeaks();		
	setBatchMode("exit and display");

	print("\\Clear");
	print("Finished:\n"+ ojNObjects() + " objects in "+ d2s((getTime - time2)/60000, 1) + " min");
	arrangeWindows();
	selectWindow("Log");
}
	
macro "-"{}



macro "Plot Fluorescence at LeaderPos vs Time "{
	if(ojNObjects() == 0)
		exit("No objects marked yet");
	germAtCenter = true;
	close("Fluor_at_LeaderPos.tif");
	Dialog.create("Plot Options");
	Dialog.addCheckbox("Germination at Center", true);
	Dialog.addNumber("Y-Plot Limit", plotLimitTop);
	Dialog.show;
	plotLimitTop = Dialog.getNumber;
	//withPrepTime =  Dialog.getCheckbox;
	germAtCenter =  Dialog.getCheckbox;
	withPrepTime = !germAtCenter; 
	if(!germAtCenter)
		xLabel = "Time[min] incl. Preparation";
	else
		xLabel = "Time[min] excl. Preparation";
	if(germAtCenter)
		xLabel = "Time[min] relative to germination";
		
 	Plot.create("Fluor_at_LeaderPos.tif", xLabel, "Fluorescence leader pos");
	
	if(germAtCenter){
		Plot.setFrameSize(700, 200);
		Plot.setLimits(-180, 180, -5, plotLimitTop );
	}
	else{
		Plot.setFrameSize(500, 200);
		Plot.setLimits(0, 190, -5, plotLimitTop );
	}
	first = 1; 
	last = ojNObjects();
	slc = 1;
	displayObj = 0;
	for(obj = first; obj <= last; obj++){
		prepTime = ojResult("prepTime", obj);
		intFls = extractArray("Spots", obj, "intFl");
		Array.getStatistics(intFls, min, max, mean, stdDev);
 		len = intFls.length;
 		rawGermT = ojResult("RawGTime", obj);
 		msg = "obj= " + obj;
		if(len > 0 && rawGermT>=0){
			if(displayObj ==0)
				displayObj = obj;
			times = extractArray("Spots", obj, "time");
			for(jj = 0; jj < len; jj++){
 				if(withPrepTime)
 					times[jj] += prepTime;
 				//else 
				if(germAtCenter)
					times[jj] -= rawGermT;
			}
			Plot.setColor("blue");
			Plot.add("line",  times, intFls);
			gTime = ojResult("RawGTime", obj);
			maxFlTime = ojResult("RawMaxFlTime", obj);
			peakToGerm = ojResult("PeakToGermTime", obj);
			msg = "obj = " + obj + "\nPeakToGerm = " + d2s(peakToGerm, 1) + " min";
			msg += "\nPeak Fluor = green, \nGermTime = red";
			if(withPrepTime){
				gTime += prepTime;
				maxFlTime += prepTime;
			}
			Plot.setLineWidth(3);
			if(germAtCenter){
				Plot.setColor("red");
				Plot.drawLine(0, -1e6, 0, 1e6);
				Plot.setColor("green");
				Plot.drawLine(maxFlTime-rawGermT, -1e6, maxFlTime-rawGermT, 1e6);
			}
			else{
				Plot.setColor("red");
				Plot.drawLine(gTime, -1e6, gTime, 1e6);
				Plot.setColor("green");
				Plot.drawLine(maxFlTime, -1e6, maxFlTime, 1e6);
			}
			
			Plot.setLineWidth(1);
			Plot.setColor("black");
		}
		Plot.addText(msg, 0.7, 0.12);
		if(first < last)
			Plot.appendToStack;

	}
	Plot.show;
	rename("Fluor_at_LeaderPos.tif");
	setSlice(displayObj);
}


*/
//Map of objects, sorted for germination time,
//showing each leader profile vs time as one pixel row.
//Non-germinating spores, weak leaders (<minPeakForMap) and unqualified objects are ignored
//Left-most pixel shows negative object number
//Fixed to 36 frames and 5 min/frame
macro "Map of Leaders"{

	getAuxPath();
	mapTitle = "Map_of_Leaders.tif";
	profileTitle = "Collective Leader Profile_" + acronym + ".tif";
	close(mapTitle);
	close(profileTitle);
	newImage(mapTitle, "32-bit black", 72, 300, 1);
	mapID = getImageID;
	run("Set...", "value=NaN");
	run("Properties...", "pixel_width=5.0000 pixel_height=1.0000 unit=min voxel_depth=1 origin=36,0");
	gTimes = ojGetColumn("RawGTime", true);
	rankPosArr = Array.rankPositions(gTimes);
	y = 0;
	for(rr = 0; rr < rankPosArr.length; rr++){
		obj = rankPosArr[rr] + 1;
		gTime = ojResult("RawGTime", obj);
		maxFlA = ojResult("maxFlA", obj);
		if(gTime >=0 && maxFlA > minPeakForMap  && ojQualified(obj)){//29.09.19, 18:57
			intFls = extractArray("Spots", obj, "intFl");
			times = extractArray("Spots", obj, "time");
			for(jj = 0; jj < times.length; jj++){
				fluor = intFls[jj];
				time = (180 + times[jj] - gTime)/5;
				setPixel(time, y, fluor);
			}
			setPixel(0, y, -obj);
			y++;
		}
	}
	height = y;
	makeRectangle(0, 0, 72, height);
	run("Crop");
	makeLine(36, 0, 36, height);
	Overlay.addSelection("cyan");
	run("Select None");
	run("Set... ", "zoom=400") ;
	setMinAndMax(-3, 32);
	run("Fire");
	
	getAuxPath();
	fullPath = auxPath + "/" + mapTitle;
	setMetadata("label", "" + height + " spores");
	saveAs("Tiff", fullPath);
	makeRectangle(1, 0, getWidth-1, getHeight);//omit first column
	peakProfile = getProfile;
	run("Select None");
	Plot.create(profileTitle, "Time relative to germination[min]", "Leader intensity");
	Plot.setLimits(-120, 120, 0, NaN);
	Plot.setFrameSize(240, 240);//30.09.19, 19:53
	time3 = newArray(-175, 180);
	time3 = Array.resample(time3, 72);
	plotID = getImageID;
	profilePath = auxPath + "/" + profileTitle;
	
	if(stepConfidence > 0){
		selectImage(mapID);
		size = 72;
		timeArr = newArray(size);
		timeArr = Array.fill(timeArr, NaN);
		meanArr = Array.copy(timeArr);
		stdArr = Array.copy(timeArr);
		nPixelsArr = Array.copy(timeArr);
		confidArr = Array.copy(timeArr);
		index = 0;
		for(x=stepConfidence; x< size-stepConfidence; x+=stepConfidence){
			makeRectangle(x, 0, stepConfidence, getHeight);
			x2 = x+stepConfidence/2;
			dummy = 0;
			toScaled(x2, dummy);
			timeArr[index] = x2;	
			mean = getValue("Mean");
			meanArr[index] = mean;
			stDev = getValue("StdDev");
			stdArr[index] = stDev;
			index++;
		}	
		run("Select None");
		run("Properties...", "pixel_width=1 pixel_height=1");
		changeValues(-9999, 1e6, 1);//leave NaN intact
		run("Enhance Contrast", "saturated=0.35");
		index = 0;
		for(x=stepConfidence; x< size-stepConfidence; x+=stepConfidence){
			makeRectangle(x, 0, stepConfidence, getHeight);
			mean= getValue("Mean");
			area= getValue("Area");
			nCount = area* mean;
			nPixelsArr[index] = nCount;
			confid = 1.96 * stdArr[index]/(sqrt(nCount-1));//95% confidence
			confidArr[index] = confid;
			index++;
		}
		run("Revert");
		run("Select None");
		timeArr = Array.trim(timeArr, index);
		meanArr = Array.trim(meanArr, index);
		stdArr = Array.trim(stdArr, index);
		nPixelsArr = Array.trim(nPixelsArr, index);
		confidArr = Array.trim(confidArr, index);
		//Array.show("Confidence Table", stdArr, nPixelsArr, confidArr, timeArr);
		selectImage(plotID);
		Plot.setColor("#bbbbbb", "#bbbbbb");
		Plot.add("dots", timeArr, meanArr);
		Plot.add("error bars",  confidArr);
	}
	Plot.setLineWidth(3);
	Plot.setColor("red", "red");
	Plot.add("line", time3, peakProfile);
	Plot.setLineWidth(1);
	Plot.show;

	setMetadata("label", "" + height + " spores");
	saveAs("Tiff", profilePath);
	close("Plot Values");
}

macro "Plot All Peaks vs GTime"{
	title = "All_Peaks-" + substring(ojGetImageName(1), 0, 3);
	allTimes = newArray(0);
	allPeaks = newArray(0);
	count = 0;
	for(obj = 1; obj <= ojNObjects(); obj++){
		rawGTime = ojResult("RawGTime", obj);
		if(!isNaN(rawGTime) && ojQualified(obj)){
			count++;
			xArr = extractArray("Spots", obj, "time");
			yArr = extractArray("Spots", obj, "intFl");
			if(xArr.length > 0){
				for(jj = 0; jj < xArr.length; jj++){
					xArr[jj] -= rawGTime;
				}
				allTimes = Array.concat(allTimes, xArr);
				allPeaks = Array.concat(allPeaks, yArr);
			}
		}
	}
	close(title);
	Plot.create(title, "time relative to germination [min]", "Intensity");
	Plot.setLimits(-180, 180, 0, NaN );
	Plot.setFrameSize(800, 500 );
	Plot.setColor("red", "red");
	Plot.add("circles", allTimes, allPeaks);
	msg = "Showing " + count + " of " + ojNObjects() + " spores";
	Plot.setColor("black");
	Plot.addText(msg, 0.7, 0.12);
	Plot.show;
	close("Plot Values");
}

macro "Plot Leader Peaks (Three Plots)"{
	getAuxPath();
	close("Leader_*");
	for(loop = 1; loop <= 3; loop++){//three different plots
		theTimes = newArray(ojNObjects());
		intensities = newArray(ojNObjects());
		objects = newArray(ojNObjects());
		hits = 0;
		for(obj = 1; obj <= ojNObjects(); obj++){
			if(!ojQualified(obj))
				continue;
			peaks = extractArray("Spots", obj, "intFl");
			nPeaks = peaks.length;//zero if not germinating
			if(nPeaks > 0){
				times = extractArray("Spots", obj, "time");
				rankPositions = Array.rankPositions(peaks);
				maxPos = rankPositions[nPeaks-1];//the largest peak
				relPeakTime = times[maxPos]- ojResult("RawGTime", obj);	//typically negative
				peak = peaks[maxPos];
				netPeak = peak ;
				objects[hits] = obj;

				if(loop == 1  || loop == 3)
					theTimes[hits] = ojResult("GermTime", obj);
				if(loop == 2)
					theTimes[hits] = relPeakTime;	
				if(loop != 3){		
					intensities[hits] = peaks[maxPos];
				}
				if(loop == 3){
					intensities[hits] = peaks[0];
				}
				hits++;
			}
		}
		theTimes = Array.trim(theTimes, hits);
		intensities = Array.trim(intensities, hits);
		objects = Array.trim(objects, hits);
		if(loop == 1){
			Plot.create("Leader_Peak_A_" + acronym , "Intensity Strongest Peak", "Time to Germination [min]");
			Plot.setColor("red", "red");
			Plot.add("circles", intensities, theTimes);
			Plot.setFrameSize(300, 250);
			Plot.setLimits(NaN, NaN, 0, NaN);
			Plot.show;
			setLocation(50, 100);
		}
		if(loop == 2){
			Plot.create("Leader_Peak_B_" + acronym, "Peak time relative to Germination  [min]", "Intensity Strongest Peak");
			Plot.setColor("blue", "blue");
			Plot.add("circles",  theTimes, intensities);
			Plot.setFrameSize(300, 250);
			Plot.show;
			setLocation(450, 100);
		}
		if(loop == 3){
			Plot.create("Leader_Peak_C_" + acronym, "Time to Germination [min]", "Intensity on Leader Position at t=0");
			Plot.setColor("#00aa00", "#00aa00");
			Plot.add("circles",  theTimes, intensities);
			Plot.setFrameSize(300, 250);
			Plot.show;
			setLocation(850, 100);
		}
	}
 	close("Plot Values");
}


macro "Histogram Germinations"{
	getAuxPath();
	title = "Histo_Mid-GermTime-" + acronym;
	close(title);
	Dialog.create("Histogram Germination Time");
	Dialog.addNumber("BinWidth [min]", binWidth);
	Dialog.addCheckbox("Relative Frequency [%]", 1);
	Dialog.show;
	binWidth = Dialog.getNumber;
	relative = Dialog.getCheckbox;
	gTimes = ojGetColumn("GermTime", false);//Set unqualified to NaN
	gTimes = Array.deleteValue(gTimes, NaN);
	gCount = ojGetStatistics("GermTime", "count");
	mean = ojGetStatistics("GermTime", "mean");
	stDev = ojGetStatistics("GermTime", "stdev");
	yLabel = "Frequency";

	if(relative)
		yLabel += " [%]";
	for(loop = 1; loop <= 2; loop++){
		Plot.create(title, "Time until mid-germination (incl. PrepTime) [min]", yLabel);
		Plot.setFrameSize(500, 300 );
		Plot.setColor("blue", "#ccccff");
		if(loop == 1){
			Plot.addHistogram(gTimes, binWidth, binWidth/2);
			Plot.show;
			Plot.getValues(xpoints, ypoints);
			if(relative){
				close;
				continue;
			}
			else
				break;
		}
		if(loop == 2 && relative){
			len = xpoints.length;
			for(jj= 0; jj < len; jj++)
				ypoints[jj] = ypoints[jj]/gCount * 100;
			Plot.setColor("blue", "#bbffdd");
			Plot.add("bar", xpoints, ypoints);

		}
	}
	all = ojNObjects() ;
	txt = "";
	txt += "\nn= " + gCount + "  (" + d2s(100*gCount/all, 1) + "% of " + all + ")";
	txt += "\nmean= " + d2s(mean, 1) + " min";
	txt += "\nstDev= " + d2s(stDev, 1) + " min";
	Plot.addText(txt, 0.7, 0.1);
	Plot.setLimits(0, 200, 0, NaN );
	if(relative)
		Plot.show;
	close("Plot Values");
}
	
macro "<Navigate To:>Previous Object [1]"{
	navigate(-1);
}

macro "<Navigate To:>Next Object [2]"{
	if(getTime - nTime < 3000)
		return;
	navigate(1);
}

macro "<Navigate To:>Current Object [3]"{
	navigate(0);
}

macro "<Navigate To:>Object # ...[4]"{
	obj = getNumber("Navigate ot Object #", ojSelectedObject());
	selectImage("Montage_Fluor.tif");
	setSlice(obj);
	navigate(0);
}


function navigate(delta){
	if(getTime - nTime < 3000)
		return;
	
	nTime = getTime;
	if(nImages == 0)
		exit;
	frontTitle = getTitle;
	goodSrc = (startsWith(frontTitle, "Montage_") || endsWith(frontTitle, "_Plots.tif") || frontTitle == "Fluor_at_LeaderPos.tif");
	
	if(!goodSrc){
		nTime = 0;
		exit("Navigating works with Montage or Time Plots in front");;
	}

	if (ojSelectedObject() > 0)
		current = ojSelectedObject();
	
	else 
		current = 0;
	obj = current + delta;
	if (obj > ojNObjects() || obj <1){
		beep;
		nTime = 0;
		exit;
	}
	setBatchMode(true);
	ojSelectObject(obj);
	titles = split("Montage_Phase.tif Montage_Fluor.tif Phase_Plots.tif Fluor_at_LeaderPos.tif");
	for (jj = 0; jj < titles.length; jj++){
		title = titles[jj];
		if(isOpen(title)){
			selectImage(title);
			setSlice(obj);
			run("Revert");
 			run("Select None");
			if(title == "Montage_Fluor.tif" && viewWhat == "smooth")
 				run("Gaussian Blur...", "sigma=&sigmaForView slice");
 			if(title == "Montage_Fluor.tif" && viewWhat == "smooth & median"){
 				run("Gaussian Blur...", "sigma=&sigmaForPeak slice");
 				run("Median...", "radius=" + medianKernelForView + " slice");
			}
		}
	}
	msg = "Obj=" + obj ;
	showStatus(msg);
	selectWindow("ImageJ");
	selectImage(frontTitle);
	nTime = 0;
	setBatchMode(false);
}


macro "<Utils>Selected Object as Table [T]"{
	close(currentTitle);
	close("Arrays");
	obj = ojSelectedObject();
	if(obj == 0) 
		exit("No object is selected");
	ojShowObject(obj);
	s = ojResultString("Spots", obj);
	print("\\Clear");
	print("---Spots string of obj = " + obj);
	s2 = replace(s, ";", ";\n");
	print(s2);
	if(s == "")
		exit("Cell " + obj + " is not germinating");
 	xG = extractArray("Spots", obj, "xG");//global (in montage)
	yG = extractArray("Spots", obj, "yG");
	xL = extractArray("Spots", obj, "xL");//local (from box center)
 	yL = extractArray("Spots", obj, "yL");
 	intFls = extractArray("Spots", obj, "intFl");
 	times = extractArray("Spots", obj, "time");
 	currentTitle = "Spots Table obj "+ obj;
  	Array.show(currentTitle, times, intFls, xG, yG, xL, yL);
 	if(intFls.length >0){
 		ranks = Array.rankPositions(intFls);
 		Array.reverse(ranks);
 		Table.setSelection(ranks[0], ranks[0], currentTitle);//select strongest
 	}
}

macro "<Utils>Arrange [F2]"{
	arrangeWindows();
}

macro "<Utils>Convolve Map"{
	selectImage("Map_of_Leaders.tif")	
	close("Convolved*");
	run("Select None");
	run("Duplicate...", "title=Convolved_1");	
	kernel = "text1=[" + getString("Kernel:", "1 1 1 1 1") + "\n] normalize";
	run("Convolve...", kernel);
	run("Select None");
	run("Set... ", "zoom=400") ;
	run("Duplicate...", "title=Convolved_2");
	run("Grays");
	run("Set... ", "zoom=400") ;
	setThreshold(10, 1e9);
	run("Set... ", "zoom=400") ;
}

macro "<Utils>Test S-Curve"{
	obj = ojSelectedObject(); 
	obj = getNumber("Object No:", obj);
	if(obj == 0) obj = 1;

 	xx = extractArray("Phase", obj, "time");
 	yy = extractArray("Phase", obj, "phase");

	len = xx.length;

	Array.getStatistics(yy, low, hi, mean, stdDev);
	sumLow=0; countLow=0; sumHi=0; countHi=0; gt = NaN;
	for(jj = 0; jj < len; jj++){
		y = yy[jj];
		if(y < low + (hi-low) * 0.1){ 
			sumLow += y;
			countLow++;
		}
		if(y > low + (hi-low) * 0.9){ 
			sumHi += y;
			countHi++;
		}
		if(y<(low+hi)/2 && isNaN(gt))
			gt = xx[jj];
	}
	low = sumLow/countLow;
	hi = sumHi/countHi;

	leftPart = low + (hi-low)/2;
	factor =  (hi-low)/2;
	simpleSEq = "y=" + leftPart + " + " + factor + "* (a-x)/(sqrt(b + (a-x) * (a-x)));";
	Fit.doFit(simpleSEq, xx, yy);
	td = Fit.p(0);
	slope =  Fit.p(1);

	Fit.plot();
	rename("S-test Obj=" + obj);
	waitForUser;
	close;
}


//prints the content of an ObjectJ Result column
macro "<Utils>Column to Log"{
	colTitle = getString("Column title", "");
	close("Log");
	print(colTitle);
	print("---");
	for(obj = 1; obj <= ojNObjects(); obj++){
		str = ojResultString(colTitle, obj);
		print(obj, "\t", str);
	}
} 

macro "Choose View Type [5]"{
	Dialog.create("Viewing Options");
	Dialog.addNumber("Gaussian Smooth   Sigma (px)", sigmaForView);
	Dialog.addNumber("Median Radius (px):", medianKernelForView);
	items = newArray("raw", "smooth", "smooth & median");
  	Dialog.addRadioButtonGroup("View:", items, 3, 1, viewWhat);
  	Dialog.addCheckbox("Show Surface Plot", false);
  	
  	msg = "---\nThese settings are only for viewing."
  	msg += "\nThey do not affect the results."
  	msg += "\nImage analysis uses: 
  	msg += "\n  Sigma for peak detection = " + sigmaForPeak +  " px";"
  	msg += "\n  Median kernel rad for cell background = "+ medianKernelForView +  " px";
  	
  	Dialog.addMessage(msg);
  	Dialog.show();
	sigmaForView = Dialog.getNumber();;
	medianKernelForView = Dialog.getNumber();
	viewWhat = Dialog.getRadioButton;
	obj = ojSelectedObject();
	if(Dialog.getCheckbox && obj> 0){//make surface plot
		close("Surface*");
		ojShowObject(obj);
		x = ojXPos(1);
		y = ojYPos(1);
		makeRectangle(x - boxR, y - boxR, boxW, boxW);
		run("Duplicate...", "duplicate channels=2 title=Tmp_obj#" +  obj);
		run("Gaussian Blur...", "sigma=&sigmaForView stack");
		if(indexOf(viewWhat, "median") >=0)
			run("Median...", "radius=&medianKernelForView stack");
		theMax = 0;
		theMin = 9e9;
		for(jj= 1; jj<= nSlices; jj++){
			setSlice(jj);
			getStatistics(area, mean, min, max);
			theMax = maxOf(max, theMax);
			theMin = minOf(min, theMin);
		}	
		run("Subtract...", "value=theMin stack");
		setMinAndMax(0, theMax-theMin);
		
		run("Surface Plot...", "polygon=66 draw_wireframe shade draw_axis stack");
		rename("Surface Plot obj=" +obj);
		setMinAndMax(-20, 160);
		run("Apply LUT", "stack");
		surfaceID = getImageID;
		setBackgroundColor(0xffffff);
		run("Canvas Size...", "height=660 position=Bottom-Left");
		close("Tmp*");

		//--- paint matrix at top of surface plot
		selectImage("Montage_Fluor.tif");
		ww = getWidth;
		hh = getHeight;
		run("View 100%");
		run("Flatten", "slice");
		run("Select All");
		run("Copy");
		close;
		selectWindow("Montage_Fluor.tif");
		run("Set... ", "zoom=200 x=180 y=80");
		selectImage(surfaceID);
		setColor("red");
		for(slc = nSlices; slc >=1; slc--){
			setSlice(slc);
			makeRectangle(100, 0, ww, hh);
			run("Paste");
			x = 100 + (slc - 1) % k9 * boxW;
			y = floor( (slc - 1)/k9) * boxW;
			drawRect(x, y, boxW, boxW);
		}
		run("Select None");
		setLocation(720,100);
		setOption("Changes", false);
		
		ojShowObject(obj);
		close;
	}
	navigate(0);
}


//Combines rcc and mcc profile plots 
//and saves result as "Combined_Profiles.tif"
// in the same folder as the rcc plot

macro "<Utils>Combine Collective Leader Profiles"{	
	waitForUser("1. Select rcc profile plot, then click OK");
	if(getInfo("window.type") != "Plot")
		exit("Plot must be in front");
	rccTitle = getTitle;
	dir1 = getDirectory("image");
	if(dir1 == "")
		exit("No path associated to " + rccTitle);
	waitForUser("2. Select mcc profile plot, then click OK");
	if(getInfo("window.type") != "Plot")
		exit("plot must be in front");
	mccTitle = getTitle;
	
	selectImage(rccTitle);
	close("Combined_Profiles.tif");
	rename("Combined_Profiles.tif");
	saveAs( dir1 + "Combined_Profiles.tif");
	
	Plot.addFromPlot(mccTitle, 0);
	Plot.addFromPlot(mccTitle, 1);	
	Plot.setStyle(1, "#eeaa22,eeaa22,3.0,Line");
	Plot.setStyle(3, "blue,blue,3.0,Line");
	Plot.setLegend(".\trcc\t.\tmcc");
	saveAs( dir1 + "Combined_Profiles.tif");
	open(dir1 + rccTitle);
}



