//RotatingCalipers.txt
//using method descibed here:
//http://cgm.cs.mcgill.ca/~orm/diam.html
//measures diameter (maximum) and width (minimum)
//distances between parallel lines of support (~tangents)
//This program is in the public domain
//Michael Doube
//Imperial College London 2008

pi = 3.14159265359;
var dmax, dmin;
run("Convex Hull");
getSelectionCoordinates(perimx, perimy);

//Compute the polygon's extreme points in the y direction
ymin = 999999;
ymax = 0;
npoints = lengthOf(perimx);
for (n = 0; n < npoints; n++){
	if (perimy[n] < ymin){
		ymin = perimy[n];
		podal = n;    		//n is always 0 in ImageJ, start point of selection is minimum x value within set with minimum y
	}
	if (perimy[n] > ymax){
		ymax = perimy[n];
		antipodal = n;
	}
}

//create list of angles beween points
angles = newArray(npoints);
for (n = 0; n < npoints - 1; n++){
	angles[n] = atan2(perimy[n]-perimy[n+1] , perimx[n]-perimx[n+1]);
}

//angle between last point and first point of selection
angles[npoints-1] = atan2(perimy[npoints-1]-perimy[0] , perimx[npoints-1]-perimx[0]);

//correct any angles < -pi
for (n=0; n<lengthOf(angles); n++) {
	if (angles[n] < -pi) angles[n] = angles[n] + 2*pi;
}

//Construct two horizontal lines of support through ymin and ymax.
//Since this is already an anti-podal pair, compute the distance, and keep as maximum.

dmax = ymax - ymin;
dmin = dmax;

end = antipodal;
while (antipodal < npoints  && podal <= end){
	podincr = 0;
	antipodincr = 0;
//find the start condition for theta
	if (podal > 0) {
		starttheta = endtheta;
	} else {
		if (angles[antipodal - 1] > 0) { 
			testangle = angles[antipodal - 1] - pi;
		} else {
			testangle = angles[antipodal - 1] + pi;
		}
		if (angles[npoints -1] <= testangle){
			starttheta = angles[npoints - 1];
		} else {
			starttheta = testangle;
		}
	}

//Rotate the lines until one is flush with an edge of the polygon.
//find the end condition for theta
//increment the podal or antipodal point depending on
//which caliper blade touches the next point first

	if (angles[antipodal] > 0) { 
		testangle = angles[antipodal] - pi;
	} else {
		testangle = angles[antipodal] + pi;
	}
	if (antipodal <= npoints - 1) {
		if (angles[podal] > testangle && angles[podal] * testangle > 0) {
			endtheta = angles[podal];
			podincr = 1;
		} else if (angles[podal] < testangle) {
			endtheta = testangle;
			antipodincr = 1;
		} else if (angles[podal] > testangle && angles[podal] * testangle < 0) {
			endtheta = testangle;
			antipodincr = 1;
		} else {
			endtheta = angles[podal];
			podincr = 1;
			antipodincr = 1;
		}
	} else {
		endtheta = 0;
	}

//A new anti-podal pair is determined. Compute the new distance,
//compare to old maximum, and update if necessary.

	thetah = atan2(perimy[antipodal] - perimy[podal] , perimx[antipodal] - perimx[podal]);
	dp = sqrt(sqr((perimx[antipodal]-perimx[podal])) + sqr((perimy[antipodal]-perimy[podal])));
	if (dp > dmax) dmax = dp;

	if (starttheta >= endtheta){
		ds = dp*abs(cos(-thetah + starttheta - pi/2));
		de = dp*abs(cos(-thetah + endtheta - pi/2));
		if (ds < dmin) dmin = ds;
		else if (de < dmin) dmin = de;
	}	
	else if (starttheta < endtheta && starttheta * endtheta < 0){
		ds = dp*abs(cos(-thetah + starttheta - pi/2));
		de = dp*abs(cos(-thetah + - 3*pi/2));
		if (ds < dmin) dmin = ds;
		else if (de < dmin) dmin = de;

		ds = dp*abs(cos(-thetah + pi/2));
		de = dp*abs(cos(-thetah + endtheta - pi/2));
		if (ds < dmin) dmin = ds;
		else if (de < dmin) dmin = de;
	}

//increment either or both podal or antipodal
podal = podal + podincr;
antipodal = antipodal + antipodincr;

}
setResult("Label", 0, getTitle());
setResult("RCmax", 0, dmax);
setResult("RCmin", 0, dmin);
updateResults();

function sqr(n) {
	return n*n;
}