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


macro "ThioBac Tool - C04fT1e12TCf00T8e12T" {
	getCursorLoc(x, y, z, flags);

	if(startsWith(getTitle, "Plot")){
		 markPerimeter(x, y);
		 exit;
	}
	if(!endsWith(toLowerCase(getTitle), ".tif")) exit("click in *.tif Window");
	analyzeShape(x, y);
}
	
function analyzeShape(x, y){	
	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;
		}
	}

	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;
	
	for(jj = startIndex; jj < len; jj++){
		phi = relAngles[jj];
		rad = radii[jj];
		if(rad > bumpVal[bumpID]){
			bumpVal[bumpID] = rad;
			bumpPos[bumpID] = jj;
		}
		if(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");
	}
	
	
	deepSensor(invagPos[0]);
	deepSensor(invagPos[1]);

	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");
	setOption("Changes", false);
	setLocation(20, 80);
	selectImage("Box");
	setBatchMode("show");
	setLocation(20, 420);
	Overlay.remove;
	run("From ROI Manager");
	roiManager("Show All without labels");
	setOption("Changes", false);
}

macro "ThioBac Tool Options" {
	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();
}


//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(x1, y1) {
	//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
		setMetadata("label", "Phi=" + round(angles360[x1]) + " deg");
	selectWindow("Box");
	if(outside){ 
		run("Select None");
		return;
	}
	makeOval(xx[x1] -4, yy[x1] - 4, 8, 8);
	run("Properties... ", "fill=magenta  size=Large");
}

//send sensor waves from invagination to center to evaluate cross profiles
function deepSensor(minPos){
	if(showDeepInvag +  showDeepSensor == 0) return;
	if(minPos == 0) return;
	wing2 = wing1 * 0.3;
	alpha = relAngles[minPos] + origAngle + 90;
	rad = radii[minPos];
	stdDev = 0;
	deepRad = 0;
	str = "";
	n = 1;
	do{
		x1 = xc + cos( (alpha -wing2)/180 * PI) * rad;
		x2 = xc + cos( (alpha +wing2)/180 * PI) * rad;
		y1 = yc - sin( (alpha -wing2)/180 * PI) * rad;
		y2 = yc - sin( (alpha +wing2)/180 * PI) * rad;
 		makeLine(x1, y1, x2, y2);
		profile = getProfile;
		if (stdDev == 0)
			Array.getStatistics(profile,min, max, mean, stdDev);
		maxima = Array.findMaxima(profile, (dispMax - dispMin)/20, true);
		found = maxima.length >= 1;
		if (found){
			xm = x1 + maxima[0]/profile.length *(x2-x1);
			ym = 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(xm, ym);
				run("Properties... ", "name=&label  stroke=#88ff00 point=Dot size=Small ");
				roiManager("add");
				
			}
			deepRad = rad;
		}
		rad -=5;//
		n++;
	}while(rad >5 && found);
	if(deepRad > 0  && 	showDeepInvag){
		x = xc + cos(alpha/180 * PI) * deepRad;
		y = yc - sin(alpha/180 * PI) * deepRad;
		makeLine(xx[minPos], yy[minPos], x, y) ;
		run("Properties... ", "name=deepL  stroke=#aaff0000");
		roiManager("add");
		makePoint(x, y) ;
 		run("Properties... ", "name=deepP  width=1 stroke=white  point=Dot size=Medium");
		roiManager("add");
	}
}


macro "Analyze Previous Object [1]"{
	showZoomedAndAnalyze(ojSelectedObject() - 1);
}

macro "Analyze Next Object [2]"{
	showZoomedAndAnalyze(ojSelectedObject() + 1);
}

macro "Grays [3]"{ run("Grays"); }

macro "Fire [4]"{run("Fire");}

macro "List Elements [5]"{ run("List Elements");}

macro "Invert LUT [6]"{ run("Invert LUT");}



	



//shows object zoomed and temporarily highlighted by circle roi
function showZoomedAndAnalyze(obj){
	if(running == 1 &&  getTime - lastTime < 2500)
		{beep; exit;}
	lastTime = getTime;
	running=1;
	
	if(obj<1 || obj>ojNObjects()){
		beep;
		exit;
	}

	rrOuter = 100;
	ojShowImage(ojOwnerIndex(obj));
	ojSelectObject(obj);
	z = ojZPos(1);
	setSlice(z);


	ojShowObject(obj);
	getLocationAndSize(winX, winY, winW, winH);
	setLocation(600, 80, 600, 600);

	x1 = ojXPos(1);
	y1 = ojYPos(1);
	x2 = ojXPos(ojNPoints());
	y2 = ojYPos(ojNPoints());
	xc = (x1 + x2)/ 2;
	yc = (y1 + y2)/ 2;
	rrInner = 10 + 0.5*sqrt(pow((x2 - x1), 2) + pow((y2 - y1), 2));
	zoomFactor = 3;
	ojZoom(zoomFactor, xc, yc);//zoom at factor 3
	makeRectangle(xc-10, yc - 10, 20, 20);
	analyzeShape(xc, yc);
	ojShowObject(obj);
	running = 0;

}