macro "Add Polar Constrictions"{
//N. Vischer 	28.8.14, 12:48 
	rad = 5; //roof radius
	roofAngle = 5; //negative angles (degrees) are concave
	setBatchMode(true);
	ojHideResults();
	ojInitColumn("PCon1");//Polar constriction
	ojInitColumn("PCon2");
	run("Close All");
	hits = 0;
	time = getTime;
	for(obj = 1; obj <= ojNObjects(); obj++){
		showProgress(obj/ojNObjects());
		ojShowObject(obj);
		close("\\Others");
		ojSetTarget("frontimage");
		ojDeleteItem("Constr", 1);
		ojDeleteItem("Constr", 1);
		ojSelectObject(obj);
		roof1X = newArray(3);
		roof1Y = newArray(3);
		roof2X = newArray(3);
		roof2Y = newArray(3);

		roof1X[1] = ojXPos(1);
		roof1Y[1] = ojYPos(1);
		roof2X[1] = ojXPos(ojNPoints());
		roof2Y[1] = ojYPos(ojNPoints());

		xstart = ojResult("startX", obj);
		ystart = ojResult("startY", obj);
		thr = ojResult("Thr", obj);
		setThreshold(0, thr);
		setThreshold(0, thr);//work-around
		doWand(xstart,ystart);
		run("Interpolate", "interval=0.3 smooth");
		resetThreshold;
		getSelectionCoordinates(xx, yy);
		run("Select None");
		len = xx.length;
		ptr1 = 0;
		ptr2 = 0;
		prevDist1 = NaN;
		prevDist2 = NaN;

		for(jj = 0; jj < (len * 1.2); jj++){
			ii = jj % len;
			dx = roof1X[1] - xx[ii];
			dy = roof1Y[1] - yy[ii];
			dist = sqrt(dx * dx + dy * dy);
			if(dist <= rad && prevDist1 > rad){
				roof1X[0] = xx[ii];
				roof1Y[0] = yy[ii];
				index1Left = ii;
			}
			if(dist > rad && prevDist1 <= rad){
				roof1X[2] = xx[ii];
				roof1Y[2] = yy[ii];
				index1Right = ii;
			}
			prevDist1 = dist;
		
			dx = roof2X[1] - xx[ii];
			dy = roof2Y[1] - yy[ii];
			dist = sqrt(dx * dx + dy * dy);
			if(dist <= rad && prevDist2 > rad){
				roof2X[0] = xx[ii];
				roof2Y[0] = yy[ii];
				index2Left = ii;
			}
			if(dist > rad && prevDist2 <= rad){
				roof2X[2] = xx[ii];
				roof2Y[2] = yy[ii];
				index2Right = ii; //look at cell 3527 !
			}
			prevDist2 = dist;
		}
		if (index1Left > index1Right)
			index1Right += len;
		index1Mid = round((index1Left + index1Right)/2)%len;

		roof1X[1] = xx[index1Mid];
		roof1Y[1] = yy[index1Mid];

		
		if (index2Left > index2Right)
			index2Right += len;
		index2Mid = round((index2Left + index2Right)/2)%len;
		
		roof2X[1] = xx[index2Mid];
		roof2Y[1] = yy[index2Mid];

		ojvInitStack("2d");
		for (kk = 0; kk <3; kk++)
			ojvPushVertex("x=" + roof1X[kk] + " y=" + roof1Y[kk]);
		deviation1 = ojvCalculate("deviation");
		if(deviation1 < roofAngle){
			ojOpenObject(obj);
			ojSwitchToItem("Constr");
			ojSetMarker(roof1X[0], roof1Y[0]);
			ojSetMarker(roof1X[1], roof1Y[1]);
			ojSetMarker(roof1X[2], roof1Y[2]);
			ojCloseObject();
			hits++;
		}

		ojvInitStack("2d");
		for (kk = 0; kk <3; kk++)
			ojvPushVertex("x=" + roof2X[kk] + " y=" + roof2Y[kk]);
		deviation2 = ojvCalculate("deviation");
		if(deviation2 < roofAngle){
			ojOpenObject(obj);
			ojSwitchToItem("Constr");
			ojSetMarker(roof2X[0], roof2Y[0]);
			ojSetMarker(roof2X[1], roof2Y[1]);
			ojSetMarker(roof2X[2], roof2Y[2]);
			ojCloseObject();
			hits++;
		}
		strong = minOf(deviation1, deviation2);
		weak = maxOf(deviation1, deviation2);
		if (strong < roofAngle)
			ojSetResult("PCon1", obj, strong);
		if (weak < roofAngle)
			ojSetResult("PCon2", obj, weak);
	}
	setBatchMode("exit and display");
	ojSetTarget("exit");
	ojShowResults();
	nPCons = ojGetStatistics("PCon1", "count")

	minutes = d2s((getTime-time)/60000, 1);
	showMessage("Cells with polar Constrictions: " + nPCons + " of "+ ojNObjects() +"\n(Time used: " + minutes + " min)");
}