//13.12.16 15:00
//version_1m

var 
wing0 = 35,//max angle for invag detector
bendThr= -10,
xx=newArray(0), yy, angles360, 
origAngle,
showDeepInvag = true,
showDeepSensor = false,
running = 0,
lastTime = 0,

eightPos,
invagPos, bumpPos,
xc, yc,
sensitivity = 0.025, //for sensor lines, e.g. 2.5% of min-max range
srcID, clickX, clickY,

macro "ThioBac Tool - C04fT1e12TCf00T8e12T" {

	if(startsWith(getTitle, "Plot")){
		 markPerimeter();
		 exit;
	}
	if(indexOf(toLowerCase(getTitle), ".tif") < 0) 
		exit("click in *.tif Window");
	getCursorLoc(clickX, clickY, z, flags);
	analyzeShape(clickX, clickY);
}


macro "ThioBac Tool Options [1]" {
	Dialog.create("Options");
	Dialog.addNumber("Bending threshold:",bendThr);
	Dialog.addCheckbox("Show Deep Invagination", showDeepInvag);
	Dialog.addCheckbox("Show Deep Sensors", showDeepSensor);
	Dialog.show();	
	bendThr = Dialog.getNumber();
	showDeepInvag = Dialog.getCheckbox();
	showDeepSensor = Dialog.getCheckbox();
	if(isOpen(srcID)){
		selectImage(srcID);
		analyzeShape(clickX, clickY);
	}
		
}

var view = 0;
var viewOptions = split("Grays,Fire,Invert LUT", ",");
macro "Change LUT [2]"{ 
	selectImage("Box");
	run(viewOptions[++view%3]); 
}

macro "List Elements [3]"{ 
	selectImage("Box");
	run("List Elements");
}



function analyzeShape(x, y){
	srcID = getImageID;
	view = 0;
	run("Line Width...", "line=1");
	setBatchMode(true);
	close("Box*");//removes threshold!
	close("Plot*");
	run("Select None");
	setAutoThreshold("Default");
	getThreshold(lower, upper);
	getVoxelSize(pxWidth, dummy, dummy, dummy);
	if(pxWidth > 0.125) exit("Scale should be >8 pixels/um");
	
	doWand(x, y);
	resetThreshold;//01.12.16 0:12
	getSelectionBounds(x, y, w, h);
	if(w==0) exit("Failed to create Box");
	box = round(4/ pxWidth);//ca 4 um
	factor = round(100 * pxWidth);
	makeRectangle(x + w/2 - box/2, y+h/2- box/2, box, box);
	run("Duplicate...", "title=Tmp");
	run("Grays");
	run("Scale...", "x=&factor y=&factor interpolation=Bilinear create title=Box");
	close("Tmp");
	run("Enhance Contrast", "saturated=0.35");
	getMinAndMax(dispMin, dispMax);
	run("Scale Bar...", "width=1 height=4 font=14 color=White background=None location=[Lower Right] bold");

	run("Set Scale...", "distance=0");

	setThreshold(lower, upper);
	doWand(box/2 * factor, box/2 * factor);
	resetThreshold;
	if(selectionType < 0) exit("Could not outline Cell \n(close to border?)");
	run("Interpolate", "interval=5 smooth adjust");
	run("Properties... ", "name=Contour  stroke=magenta");
	roiManager("reset");
	roiManager("add");
	List.setMeasurements;
	
	origAngle = List.getValue("Angle");
	ar2 = List.getValue("AR")/2;//half aspect ratio
	xc = List.getValue("X");
	yc = List.getValue("Y");
	wing1 = wing0/ar2;//angular window for invag finder
	makePoint(xc, yc);
	run("Properties... ", "name=Center  stroke=blue point=Dot size=Medium");
	roiManager("add");



	roiManager("select", 0);//Evaluate contour and radius along perimeter
	getSelectionCoordinates(xx, yy);
	len = xx.length;
	gamma = newArray(len);
	radii = newArray(len);
	relAngles = newArray(len);
	angles360 = newArray(len);
	rollArrays(xx, yy, xc, yc, origAngle, relAngles, angles360);//+++++ roll Arrays
	arm = 4;//default: 4 px;  angle of 2 connected arms is calculated 
	curvature = getVertexAngles(xx, yy, arm);
	for(jj = 0; jj < len; jj++)
		radii[jj] = sqrt(pow((xx[jj] - xc), 2) + pow((yy[jj] - yc), 2));
	
	startIndex = 0;
	for(jj = 0; jj < len; jj++){
		if(relAngles[jj] < 0){
			startIndex = jj;
			break;
		}
	}
	//eight positions along perimeter; 0-90-180-270 deg and "left-right" evaluation angles
	eightPos = newArray(8);//pos0, left1, pos90, right1, pos180, left2, pos270, right2
	for(jj = startIndex; jj < len; jj++){
		psi = angles360[jj];
		eightPos[0] = 0;//starts at ellipse short axis
		if(psi >  90 - wing1  && eightPos[1] == 0)
			eightPos[1] = jj-1;
		if(psi >  90  && eightPos[2] == 0)
			eightPos[2] = jj-1;
		if(psi >  90 + wing1  && eightPos[3] == 0)
			eightPos[3] = jj-1;

		if(psi > 180  && eightPos[4] == 0)
			eightPos[4] = jj-1;

		if(psi >  270 - wing1  && eightPos[5] == 0)
			eightPos[5] = jj-1;
		if(psi >  270  && eightPos[6] == 0)
			eightPos[6] = jj-1;
		if(psi >  270 + wing1  && eightPos[7] == 0)
			eightPos[7] = jj-1;
	}

	
	invagPos = newArray(2);
	invagVal= newArray(999, 999);
	bumpPos = newArray(4);
	bumpVal = newArray(4);


	for(jj = startIndex; jj < len; jj++){
		inSector1 = jj > eightPos[1] && jj < eightPos[3];
		if(inSector1 && curvature[jj] < invagVal[0] && curvature[jj] < bendThr ){
			invagVal[0] = curvature[jj];
			invagPos[0] = jj;
		}
		inSector2 = jj > eightPos[5] && jj < eightPos[7];
		if(inSector2 && curvature[jj] < invagVal[1] && curvature[jj] < bendThr ){
			invagVal[1] = curvature[jj];
			invagPos[1] = jj;
		}
	}
	bumpID = 0;
//debug;
	for(jj = startIndex; jj < len; jj++){
		phi = relAngles[jj];
		rad = radii[jj];
		if(rad > bumpVal[bumpID]){
			bumpVal[bumpID] = rad;
			bumpPos[bumpID] = jj;
		}
		if(jj > 0 && (jj == invagPos[0] || jj == invagPos[1]))
			bumpID = bumpID | 1;//odd ID
		if(phi > 0 && bumpID <2)
			bumpID = 2;
	}
	if(invagPos[0] == 0){
		bumpPos[0] = eightPos[2];
		bumpPos[1] = 0;
	}
	if(invagPos[1] == 0){
		bumpPos[2] = eightPos[6];
		bumpPos[3] = 0;
	}	
	
	if(eightPos[1]== eightPos[3] ||eightPos[5]== eightPos[7]) 
		exit("Failed to proceed");
	pole1xx = Array.slice(xx, eightPos[1], eightPos[3]);
	pole1yy = Array.slice(yy, eightPos[1], eightPos[3]);
	makeSelection("line", pole1xx, pole1yy);
	run("Properties... ", "name=PoleU  stroke=green");
	roiManager("add");
	pole2xx = Array.slice(xx, eightPos[5], eightPos[7]);
	pole2yy = Array.slice(yy, eightPos[5], eightPos[7]);
	makeSelection("line", pole2xx, pole2yy);
	run("Properties... ", "name=PoleL  stroke=cyan");
	roiManager("add");
	
	

	makeLine(xx[eightPos[0]], yy[eightPos[0]], xx[eightPos[4]], yy[eightPos[4]]);
	run("Properties... ", "name=CrossAxis  stroke=blue");
	roiManager("add");

	if(bumpPos[0] > 0){
		makePoint(xx[bumpPos[0]], yy[bumpPos[0]]);
		run("Properties... ", "name=BumpU1  stroke=orange point=Dot size=Medium");
		roiManager("add");
	}
	
	if(bumpPos[1] > 0){
		makePoint(xx[bumpPos[1]], yy[bumpPos[1]]);
		run("Properties... ", "name=BumpU2  stroke=orange point=Dot size=Medium");
		roiManager("add");
	}
	
	if(bumpPos[2] > 0){
		makePoint(xx[bumpPos[2]], yy[bumpPos[2]]);
		run("Properties... ", "name=BumpL1  stroke=yellow point=Dot size=Medium");
		roiManager("add");
	}
	
	if(bumpPos[3] > 0){
		makePoint(xx[bumpPos[3]], yy[bumpPos[3]]);
		run("Properties... ", "name=BumpL2  stroke=yellow point=Dot size=Medium");
		roiManager("add");
	}
	
	if(invagPos[0] > 0){
		makePoint(xx[invagPos[0]], yy[invagPos[0]]);
		run("Properties... ", "name=InvagU  stroke=green point=Dot size=Medium");
		roiManager("add");
	}
	
	if(invagPos[1] > 0){
		makePoint(xx[invagPos[1]], yy[invagPos[1]]);
		run("Properties... ", "name=InvagL  stroke=cyan point=Dot size=Medium");
		roiManager("add");
	}
	
	
	doDeepSensor(invagPos[0], 1);
	doDeepSensor(invagPos[1], 2);

	Plot.create("Plot1", "Perimeter", "Curvature");
	Plot.setLineWidth(2);
	Plot.setColor("blue");

	Plot.add("line", curvature);
	Plot.setColor("red");
	Plot.add("line", radii);

	Plot.setColor("orange");
	bendX = newArray(0, len);
	bendY = newArray(bendThr, bendThr);
	
	yRange = 180;
	Plot.add("line", bendX, bendY);
	Plot.setLegend("Curvature\nRadius\nBend.Thr");

	Plot.setLimits(0, len, -yRange, yRange);

	Plot.setFrameSize(400, 250);
	Plot.show;
	plotShadeY = newArray(yRange, -yRange, yRange, -yRange);
	plotShadeX = newArray(4);//borders of shaded areas in plot
	plotShadeX[0] = eightPos[1];
	plotShadeX[1] = eightPos[3];
	plotShadeX[2] = eightPos[5];
	plotShadeX[3] = eightPos[7];
	toUnscaled(plotShadeX, plotShadeY);//shaded areas in plot
	makeRectangle(plotShadeX[0],plotShadeY[0], plotShadeX[1] - plotShadeX[0], plotShadeY[1] - plotShadeY[0]);
	changeValues(0xffffff, 0xffffff, 0xddffdd); 
	makeRectangle(plotShadeX[2],plotShadeY[2], plotShadeX[3] - plotShadeX[2], plotShadeY[3] - plotShadeY[2]);
	changeValues(0xffffff, 0xffffff, 0xddddff); 
	run("Select None");
	setBatchMode("show");//Plot
	setOption("Changes", false);
	setLocation(20, 80);
	selectImage("Box");
	setBatchMode("show");//Box
	setLocation(20, 420);
	Overlay.remove;
	run("From ROI Manager");
	roiManager("Show All without labels");
	setOption("Changes", false);
	selectImage(srcID);
	run("Enlarge...", "enlarge=0");//force update
}


//Make sure perimeter starts at minor axis
function rollArrays(xx, yy, xc, yc, origAngle, relAngles, angles360){
	alpha = origAngle + 90;
	if(alpha >= 180)
		alpha -= 360;
	len = xx.length;
	minDPhi = 360;
	for(jj = 0; jj <len; jj++){
		dx = xx[jj] - xc;
		dy = -(yy[jj] - yc);
		phi = atan2(dy, dx) * 180/PI;
		dPhi = (phi - alpha);
		if(dPhi > 180)
			dPhi -= 360;
		if(dPhi < -180)
			dPhi += 360;
		relAngles[jj] = dPhi;
		if(minDPhi > abs(dPhi)){
			minDPhi = abs(dPhi);
			startIndex = jj;
		}
	}
	tmpXx = Array.copy(xx);
	tmpYy = Array.copy(yy);
	tmpAngles = Array.copy(relAngles);

	for(jj = 0; jj <len; jj++){
		kk = (jj + startIndex + len)%len;
		xx[jj] = tmpXx[kk];
		yy[jj] = tmpYy[kk];
		relAngles[jj] = tmpAngles[kk];
		angles360[jj] = -relAngles[jj];
		if(angles360[jj]< 0)
			angles360[jj]+= 360;
	}
}

//return curvature in x an y arrays, using arm length
function  getVertexAngles(x, y, arm) {
	len = x.length;
	vAngles = newArray(len);
	for ( mid = 0; mid < len; mid++) {
		left = (mid + len - arm) % len;
		right = (mid + arm) % len;
		dotprod = (x[right] - x[mid]) * (x[left] - x[mid]) + (y[right] - y[mid]) * (y[left] - y[mid]);
		crossprod = (x[right] - x[mid]) * (y[left] - y[mid]) - (y[right] - y[mid]) * (x[left] - x[mid]);
		phi = 180 - 180/PI *  atan2(crossprod, dotprod);
		while(phi>180)
			phi -= 360;
		vAngles[mid] = phi;
	}
	return vAngles;
}

//marks point on perimeter  when user clicks in Plot window
function markPerimeter() {
	getCursorLoc(x1, y1, z, flags);

	makeRectangle(x1 -1, 15,  2, 250);
	run("Properties... ", "fill=red");
	toScaled(x1, y1);
	outside = x1 < 0 || x1 >= xx.length;
	if(outside) 
		run("Select None");
	else{
		phiLocal = round(angles360[x1]);
		phiGlobal = round(origAngle-phiLocal +90);
		if(phiGlobal >= 180)
			phiGlobal -= 360;
		setMetadata("label", "phiLocal=" + phiLocal +" phiGlobal=" + phiGlobal + " deg   ");
	}
	selectWindow("Box");
	if(outside){ 
		run("Select None");
		return;
	}
	makeOval(xx[x1] -5, yy[x1] - 5, 10, 10);
	run("Properties... ", "fill=magenta  ");
}





//send sensor waves from invagination to center to evaluate cross profiles
function doDeepSensor(minPos, pole){//pole is 1 or 2
	if(showDeepInvag +  showDeepSensor == 0) return;
	if(minPos == 0) return;
	alpha = relAngles[minPos] + origAngle + 90;
	rad = radii[minPos];
	stdDev = 0;
	str = "";
	n = 1;

	midBumpsX = 0.5 * (xx[bumpPos[pole * 2 -2]] + xx[bumpPos[pole * 2 -1]]);
	midBumpsY = 0.5 * (yy[bumpPos[pole * 2 -2]] + yy[bumpPos[pole * 2 -1]]);
	
	xOpp = xc + (xc - midBumpsX);//opposite side
	yOpp = yc + (yc - midBumpsY);

	k10 = 15;
	maximaX = newArray(20); 
	maximaY = newArray(20); 
	lastLine = 0;
	for(line = 0; line < k10; line++){
		a = 1 - line / k10/2; //1 .. 0.5
		b = 1 - a; //0 .. 0.5
		x1 = a * xx[bumpPos[pole * 2 -2]] +b * xOpp;
		y1 = a  * yy[bumpPos[pole * 2 -2]] + b * yOpp;
		
		x2 = a * xx[bumpPos[pole * 2 -1]] + b * xOpp;
		y2 = a * yy[bumpPos[pole * 2 -1]] + b * yOpp;
		makeLine(x1, y1, x2, y2);//left
		profile = getProfile;
		if (stdDev == 0)
			Array.getStatistics(profile,min, max, mean, stdDev);
		maxima = Array.findMaxima(profile, (dispMax - dispMin) * sensitivity, true);
		found = maxima.length >= 1;
		if(!found) break;//07.12.16 13:10
		if (found){
			lastLine = line;
			maximaX[line] = x1 + maxima[0]/profile.length *(x2-x1);
			maximaY[line]  = y1 + maxima[0]/profile.length *(y2-y1);

			if( angles360[minPos] < 180) label = "sensorU" + n;
			else label = "sensorL" + n;
			if(	showDeepSensor){
				run("Properties... ", "name=&label  stroke=#88ffff00 ");
				roiManager("add");
				
				makePoint(maximaX[line], maximaY[line]);
				run("Properties... ", "name=&label  stroke=#88ff00 point=Dot size=Small ");
				roiManager("add");
			}
		}
	}
	if(lastLine > 0  && 	showDeepInvag){	
		makePoint(maximaX[lastLine], maximaY[lastLine]);
		run("Properties... ", "name=&label  stroke=white point=Dot size=Large ");
		roiManager("add");
	}
}