//17.11.16 11:11
//version 01
var radius = 10, 
wing0 = 30,
bendThr=0,
xx=newArray(0), yy,

macro "ThioBac Tool - C04fV117fV717f" {
	if(startsWith(getTitle, "Plot")){
		 markPerimeter();
		 exit;
	}
	if(!endsWith(getTitle, ".tif")) exit("click in *.tif Window");
	
	setBatchMode(true);
	close("Box*");//removes threshold!
	close("Plot*");

	getThreshold(lower, upper);
	if(lower < 0){ 
		run("Select None");
		setAutoThreshold("Default");
		getThreshold(lower, upper);
	}
	getVoxelSize(pxWidth, dummy, dummy, dummy);
	if(pxWidth > 0.125) exit("Scale should be >8 pixels/um");
	
	getCursorLoc(x, y, z, flags);
	doWand(x, y);
	
	getSelectionBounds(x, y, w, h);
	box = round(6/ pxWidth);// 6um
	factor = round(60 * pxWidth);//15.11.16 12:21 
	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("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 auto-outline Cell");
	run("Interpolate", "interval=5 smooth adjust");
	run("Properties... ", "  stroke=magenta");
	roiManager("reset");
	roiManager("add");
	List.setMeasurements;
	
	alpha = List.getValue("Angle");
	ar2 = List.getValue("AR")/2;//half aspect ratio
	xc = List.getValue("X");
	yc = List.getValue("Y");
	
	
	wing1 = wing0/ar2;//Create evaluation triangles
	colors = split("#00dd00 #00aaff");
	angles = newArray(-wing1, wing1);
	for (tt = 0; tt < 2; tt++){
		triangleX = newArray(3);
		triangleY = newArray(3);

		triangleX[2] = xc;
		triangleY[2] = yc;

		for (pp = 0; pp < 2; pp++){
			beta = (angles[pp ]+ alpha) / 180 * PI;
			dx = cos(beta) * factor * box/3;
			dy = -sin(beta) * factor * box/3;
			triangleX[pp] = xc + dx;
			triangleY[pp] = yc + dy;
		}
		makeSelection("polygon", triangleX, triangleY); 
		run("Properties... ", "  stroke=" +colors[tt]);
		roiManager("add");
		angles[0] += 180;
		angles[1] += 180;
	}

	roiManager("select", 0);//Evaluate contour and radius along perimeter
	getSelectionCoordinates(xx, yy);
	len = xx.length;
	gamma = newArray(len);
	distance = newArray(len);
	relAngles = newArray(len);
	rollArrays(xx, yy, xc, yc, alpha, relAngles);
	arm = 3;
	curvature = getVertexAngles(xx, yy, arm);
	for(jj = 0; jj < len; jj++)
		distance[jj] = sqrt(pow((xx[jj] - xc), 2) + pow((yy[jj] - yc), 2));

	poleTags = newArray(len);  //Mark poles in arrays
	startStop = newArray(len, 0, len, 0);
	for(jj = 0; jj < len; jj++){
		psi =relAngles[jj];
		if(psi >-90 - wing1 && psi <-90 + wing1 ){
			poleTags[jj] = 1;
			startStop[0] = minOf(startStop[0], jj);
			startStop[1] = maxOf(startStop[1], jj);
		}
		if(psi >90 - wing1 && psi <90 + wing1 ){
			poleTags[jj] = 2;
			startStop[2] = minOf(startStop[2], jj);
			startStop[3] = maxOf(startStop[3], jj);
		}
	}
	
	minima = Array.findMinima(curvature, 2); //Find invaginations
	found= "";
	for(kk = 0; kk<minima.length; kk++){
		minPos = minima[kk];
		
		if(poleTags[minPos] > 0 && curvature[minPos]< bendThr ){
			pole = ""+ poleTags[minPos];
			if((indexOf(found, pole)) == -1)//only strongest on either pole
				makePoint(xx[minPos], yy[minPos]);
			found += pole;
			roiManager("add");
		}
	}
	
	Plot.create("Plot1", "Perimeter", "Curvature");
	Plot.setColor("blue");
	Plot.add("line", curvature);
	Plot.setColor("red");
	Plot.add("line", distance);
	Plot.setColor("orange");
	bendX = newArray(0, len);
	bendY = newArray(bendThr, bendThr);
	yRects = newArray(120, -120, 120, -120);
	Plot.add("line", bendX, bendY);
	Plot.setLimits(0, len, -120, 120);
	Plot.setLegend("Curvature\nRadius\nBend.Thr");
	Plot.setFrameSize(400, 250);
	Plot.show;

	toUnscaled(startStop, yRects);//shaded areas in plot
	makeRectangle(startStop[0],yRects[0], startStop[1] - startStop[0], yRects[1] - yRects[0]);
	changeValues(0xffffff, 0xffffff, 0xddffdd); 
	makeRectangle(startStop[2],yRects[2], startStop[3] - startStop[2], yRects[3] - yRects[2]);
	changeValues(0xffffff, 0xffffff, 0xddddff); 
	run("Select None");
	setBatchMode("show");
	setLocation(20, 80);
	selectImage("Box");
	setBatchMode("show");
	setLocation(20, 420);

	roiManager("Show All without labels");
}

macro "ThioBac Tool Options" {
   bendThr = getNumber("Bending Threshold: ", bendThr);
}


function rollArrays(xx, yy, xc, yc, alpha, relAngles){//poles should be symmetrical in plot
	alpha += 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];
	}
}


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;
}


function markPerimeter() {
	getCursorLoc(x1, y1, z, flags);
	makeRectangle(x1 -1, 15,  2, 250);
	run("Properties... ", "fill=red");

	toScaled(x1, y1);
	selectWindow("Box");
	if(x1 >= 0 && x1 < xx.length){
		makeOval(xx[x1] -4, yy[x1] - 4, 8, 8);
		run("Properties... ", "fill=red");
	}
}