/* Copyright (C) 2001-2007 Peter Selinger and nitoyon.
   Original code(Potrace v1.8) by Peter Selinger.
   Ported to ActionScript 3.0 by nitoyon.
   This file is part of PotrAs. It is free software and it is covered
   by the GNU General Public License. See the file COPYING for details. */

package com.nitoyon.potras
{
	import flash.display.BitmapData;
	import flash.geom.*;

	public class ProcessPath
	{
		/**
		 *  it suffices that this is longer than any
		 *  path; it need not be really infinite
		 */
		private static const INFTY:int = 10000000;

		public static const POTRACE_CURVETO:int = 1;
		public static const POTRACE_CORNER:int = 2;

		/**
		 * On success, returns a Potrace state st with st->status ==
		 * POTRACE_STATUS_OK. On failure, returns NULL if no Potrace state
		 * could be created (with errno set), or returns an incomplete Potrace
		 * state (with st->status == POTRACE_STATUS_INCOMPLETE). Complete or
		 * incomplete Potrace state can be freed with potrace_state_free().
		 */
		public static function processPath(pathList:Array):ClosedPathList
		{
			var curveList:ClosedPathList = new ClosedPathList();

			for(var i:int = 0; i < pathList.length; i++)
			{
				var sums:Array       = ProcessPath.calcSums(pathList[i].priv as Array) as Array;
				var lons:Array       = ProcessPath.calcLon(pathList[i].priv as Array);
				var po:Array         = ProcessPath.bestPolygon(pathList[i].priv as Array, lons, sums);
				var vertex:Array     = ProcessPath.adjustVertices(pathList[i].priv as Array, sums, po);
				var curve:ClosedPath = ProcessPath.smooth(vertex, pathList[i].sign, .9);

				curveList.$a.push(curve);
			}

			return curveList;
		}

		/**
		 * Preparation: fill in the sum* fields of a path (used for later
		 * rapid summing). Return 0 on success, 1 with errno set on
		 * failure.
		 */
		public static function calcSums(pt:Array):Array
		{
			var n:int = pt.length;

			var sums:Array = new Array(n + 1);

			// origin
			var x0:int = pt[0].x;
			var y0:int = pt[0].y;

			// preparatory computation for later fast summing
			sums[0] = new Sums();
			sums[0].x2 = sums[0].xy = sums[0].y2 = sums[0].x = sums[0].y = 0;
			for(var i:int = 0; i < n; i++)
			{
				var x:int = pt[i].x - x0;
				var y:int = pt[i].y - y0;
				sums[i + 1] = new Sums();
				sums[i + 1].x  = sums[i].x  + x;
				sums[i + 1].y  = sums[i].y  + y;
				sums[i + 1].x2 = sums[i].x2 + x*x;
				sums[i + 1].xy = sums[i].xy + x*y;
				sums[i + 1].y2 = sums[i].y2 + y*y;
			}


			return sums;
		}

		/**
		 * Stage 1: determine the straight subpaths (Sec. 2.2.1). Fill in the
		 * "lon" component of a path object (based on pt/len).        For each i,
		 * lon[i] is the furthest index such that a straight line can be drawn
		 * from i to lon[i]. Return 1 on error with errno set, else 0.
		 */
		public static function calcLon(pt:Array):Array
		{
			var lon:Array = [];
			var n:int = pt.length;
			var j:Number;

			// initialize the nc data structure. Point from each point to the
			// furthest future point to which it is connected by a vertical or
			// horizontal segment. We take advantage of the fact that there is
			// always a direction change at 0 (due to the path decomposition
			// algorithm). But even if this were not so, there is no harm, as
			// in practice, correctness does not depend on the word "furthest"
			// above.
			var nc:Array = [];
			var k:int = 0;
			for(var i:int = n - 1; i >= 0; i--)
			{
				if(pt[i].x != pt[k].x && pt[i].y != pt[k].y)
				{
					k = i + 1; // necessarily i<n-1 in this case
				}
				nc[i] = k;
			}

			// determine pivot points: for each i, let pivk[i] be the furthest k
			// such that all j with i<j<k lie on a line connecting i,k. */
			var pivk:Array = [];
			var ct:Array;
			var dir:int;
			var k1:int;
			var constraint0:Point = new Point();
			var constraint1:Point = new Point();
			var cur:Point = new Point();
			var off:Point = new Point();
			var dk:Point = new Point();
			for(i = n - 1; i >= 0; i--)
			{
				ct = [0, 0, 0, 0];

				// keep track of "directions" that have occurred
				dir = (3 + 3 * (pt[mod(i + 1, n)].x - pt[i].x) + (pt[mod(i + 1, n)].y - pt[i].y)) / 2;
				ct[dir]++;
				constraint0.x = constraint0.y = constraint1.x = constraint1.y = 0;

				// find the next k such that no straight line from i to k
				k = nc[i];
				k1 = i;
				var foundk:Boolean = false;
				while(true)
				{
					dir = (3 + 3 * sign(pt[k].x - pt[k1].x) + sign(pt[k].y - pt[k1].y)) / 2;
					ct[dir]++;

					// if all four "directions" have occurred, cut this path
					if(ct[0] && ct[1] && ct[2] && ct[3])
					{
						pivk[i] = k1;
						foundk = true;
						break;
					}

					cur.x = pt[k].x - pt[i].x;
					cur.y = pt[k].y - pt[i].y;

					// see if current constraint is violated
					if(xprod(constraint0, cur) > 0 || xprod(constraint1, cur) < 0)
					{
						break;
					}

					// else, update constraint
					if (Math.abs(cur.x) <= 1 && Math.abs(cur.y) <= 1)
					{
						// no constraint
					}
					else
					{
						off.x = cur.x + ((-cur.y >= 0 && (-cur.y > 0 || cur.x < 0)) ? 1 : -1);
						off.y = cur.y + ((-cur.x <= 0 && (-cur.x < 0 || cur.y < 0)) ? 1 : -1);
						if(xprod(constraint0, off) <= 0)
						{
							constraint0.x = off.x;
							constraint0.y = off.y;
						}

						off.x = cur.x + ((-cur.y <= 0 && (-cur.y < 0 || cur.x < 0)) ? 1 : -1);
						off.y = cur.y + ((-cur.x >= 0 && (-cur.x > 0 || cur.y < 0)) ? 1 : -1);
						if(xprod(constraint1, off) >= 0)
						{
							constraint1.x = off.x;
							constraint1.y = off.y;
						}
					}

					k1 = k;
					k = nc[k1];
					if(!cyclic(k, i, k1))
					{
						break;
					}
				}

				// constraint
				if(!foundk)
				{
					dk.x = sign(pt[k].x-pt[k1].x);
					dk.y = sign(pt[k].y-pt[k1].y);
					cur.x = pt[k1].x - pt[i].x;
					cur.y = pt[k1].y - pt[i].y;

					// find largest integer j such that xprod(constraint0, cur+j*dk) >= 0 
					// and xprod(constraint1, cur+j*dk) <= 0. Use bilinearity of xprod.
					var a:int = xprod(constraint0, cur);
					var b:int = xprod(constraint0, dk);
					var c:int = xprod(constraint1, cur);
					var d:int = xprod(constraint1, dk);

					// find largest integer j such that a+j*b>=0 and c+j*d<=0.
					// This can be solved with integer arithmetic.
					j = INFTY;
					if(b > 0)
					{
						j = floordiv(-a, b);
					}
					if(d < 0)
					{
						j = Math.min(j, floordiv(c, -d));
					}
					pivk[i] = mod(k1 + j, n);
				}

				// foundk
			} // for i

			// clean up: for each i, let lon[i] be the largest k such that for
			// all i' with i<=i'<k, i'<k<=pivk[i'].
			j = pivk[n - 1];
			lon[n - 1] = j;
			for(i = n - 2; i >= 0; i--)
			{
				if(cyclic(i + 1, pivk[i], j))
				{
					j = pivk[i];
				}
				lon[i] = j;
			}

			for(i = n - 1; cyclic(mod(i + 1, n), j, lon[i]); i--)
			{
				lon[i] = j;
			}

			return lon;
		}

		/**
		 * Stage 2: calculate the optimal polygon (Sec. 2.2.2-2.2.4).
		 *
		 * find the optimal polygon. Fill in the m and po components. Return 1
		 * on failure with errno set, else 0. Non-cyclic version: assumes i=0
		 * is in the polygon. Fixme: ### implement cyclic version.
		 */
		public static function bestPolygon(pt:Array, lon:Array, sums:Array):Array
		{
			var i:int, j:int, k:int;
			var n:int = pt.length;

			var clip0:Array = new Array(n);		// clip0[n]: longest segment pointer, non-cyclic
			var clip1:Array = new Array(n + 1);	// clip1[n+1]: backwards segment pointer, non-cyclic

			// calculate clipped paths
			for(i = 0; i < n; i++)
			{
				var c:int = mod(lon[mod(i - 1, n)] - 1, n);
				if(c == i)
				{
					c = mod(i + 1, n);
				}
				if(c < i)
				{
					clip0[i] = n;
				}
				else
				{
					clip0[i] = c;
				}
			}

			// calculate backwards path clipping, non-cyclic. j <= clip0[i] iff clip1[j] <= i, for i,j=0..n.
			j = 1;
			for(i = 0; i < n; i++)
			{
				while(j <= clip0[i])
				{
					clip1[j] = i;
					j++;
				}
			}

			var seg0:Array = new Array(n + 1);		// seg0[m+1]: forward segment bounds, m<=n
			var seg1:Array = new Array(n + 1);	 	// seg1[m+1]: backward segment bounds, m<=n

			// calculate seg0[j] = longest path from 0 with j segments
			i = 0;
			for(j = 0; i < n; j++)
			{
				seg0[j] = i;
				i = clip0[i];
			}
			seg0[j] = n;
			var m:int = j;

			// calculate seg1[j] = longest path to n with m-j segments
			i = n;
			for(j = m; j > 0; j--)
			{
				seg1[j] = i;
				i = clip1[i];
			}
			seg1[0] = 0;

			// now find the shortest path with m segments, based on penalty3
			// note: the outer 2 loops jointly have at most n interations, thus
			// the worst-case behavior here is quadratic. In practice, it is
			// close to linear since the inner loop tends to be short.
			var pen:Array   = new Array(n + 1); // pen[n+1]: penalty vector
			var prev:Array  = new Array(n + 1);	// prev[n+1]: best path pointer vector
			var thispen:Number;
			pen[0] = 0;
			for(j = 1; j <= m; j++)
			{
				for(i = seg1[j]; i <= seg0[j]; i++)
				{
					var best:Number = -1;
					for(k = seg0[j - 1]; k >= clip1[i]; k--)
					{
						thispen = penalty3(pt, k, i, sums) + pen[k];
						if(best < 0 || thispen < best)
						{
							prev[i] = k;
							best = thispen;
						}
					}
					pen[i] = best;
				}
			}

			// read off shortest path
			var po:Array = new Array(m);
			for(i = n, j = m - 1; i > 0; j--)
			{
				po[j] = i = prev[i];
			}

			return po;
		}

		/**
		 * Auxiliary function: calculate the penalty of an edge from i to j in
		 * the given path. This needs the "lon" and "sum*" data.
		 */
		private static function penalty3(pt:Array, i:int, j:int, sums:Array):Number
		{
			var n:int = pt.length;

			// assume 0 <= i < j <= n
			var r:int = 0; // rotations from i to j

			if(j >= n)
			{
				j -= n;
				r += 1;
			}
			
			var x:Number  = sums[j + 1].x  - sums[i].x  + r * sums[n].x;
			var y:Number  = sums[j + 1].y  - sums[i].y  + r * sums[n].y;
			var x2:Number = sums[j + 1].x2 - sums[i].x2 + r * sums[n].x2;
			var xy:Number = sums[j + 1].xy - sums[i].xy + r * sums[n].xy;
			var y2:Number = sums[j + 1].y2 - sums[i].y2 + r * sums[n].y2;
			var k:Number  = j + 1 - i + r * n;

			var px:Number = (pt[i].x + pt[j].x) / 2.0 - pt[0].x;
			var py:Number = (pt[i].y + pt[j].y) / 2.0 - pt[0].y;
			var ey:Number = (pt[j].x - pt[i].x);
			var ex:Number =  -(pt[j].y - pt[i].y);

			var a:Number = ((x2 - 2 * x * px) / k + px * px);
			var b:Number = ((xy - x * py - y * px) / k + px * py);
			var c:Number = ((y2 - 2 * y * py) / k + py * py);
			
			var s:Number = ex * ex * a  +  2 * ex * ey * b  +  ey * ey * c;

			return Math.sqrt(s);
		}

		/**
		/* Stage 3: vertex adjustment (Sec. 2.3.1).
		 *
		 * Adjust vertices of optimal polygon: calculate the intersection of
		 * the two "optimal" line segments, then move it into the unit square
		 * if it lies outside. Return 1 with errno set on error; 0 on
		 * success.
		 */
		public static function adjustVertices(pt:Array, sums:Array, po:Array):Array
		{
			var m:int = po.length;
			var n:int = pt.length;

			var i:int, j:int, k:int, l:int;

			// represent each line segment as a singular quadratic form; the
			// distance of a point (x,y) from the line segment will be
			// (x,y,1)Q(x,y,1)^t, where Q=q[i].
			var q:Array = [];
			var v:Point3D = new Point3D();
			var ctr:Point = new Point();
			var dir:Point = new Point();
			for(i = 0; i < m; i++)
			{
				// calculate "optimal" point-slope representation for each line segment
				j = po[(i + 1) % m];
				j = mod(j - po[i], n) + po[i];
				ctr.x = ctr.y = dir.x = dir.y = 0;
				pointslope(pt, sums, po[i], j, ctr, dir);

				q[i] = new QuadraticForm();
				var d1:Number = dir.x * dir.x + dir.y * dir.y;
				if(d1 != 0.0)
				{
					v[0] = dir.y;
					v[1] = -dir.x;
					v[2] = -v[1] * ctr.y - v[0] * ctr.x
					q[i].fromVectorMultiply(v).scalar(1 / d1);
				}
			}

			// now calculate the "intersections" of consecutive segments.
			// Instead of using the actual intersection, we find the point
			// within a given unit square which minimizes the square distance to
			// the two lines.
			var vertex:Array = [];
			var p0:Point = pt[0].clone();
			var s:Point = new Point();
			var w:Point = new Point();
			var _q:QuadraticForm = new QuadraticForm();
			var minCoord:Point = new Point;
			for(i = 0; i < m; i++)
			{
				var z:int;
				vertex[i] = new Point();

				// let s be the vertex, in coordinates relative to x0/y0
				s.x = pt[po[i]].x - p0.x;
				s.y = pt[po[i]].y - p0.y;

				// intersect segments i-1 and i
				j = mod(i - 1, m);

				// add quadratic forms
				var Q:QuadraticForm = q[j].clone().add(q[i]);

				while(1)
				{
					// minimize the quadratic form Q on the unit square
					// find intersection
					var det:Number = Q[0][0] * Q[1][1] - Q[0][1] * Q[1][0];
					if(det != 0.0)
					{
						w.x = (-Q[0][2] * Q[1][1] + Q[1][2] * Q[0][1]) / det;
						w.y = ( Q[0][2] * Q[1][0] - Q[1][2] * Q[0][0]) / det;
						break;
					}

					// matrix is singular - lines are parallel. Add another,
					// orthogonal axis, through the center of the unit square
					if(Q[0][0] > Q[1][1])
					{
						v[0] = -Q[0][1];
						v[1] = Q[0][0];
					}
					else if(Q[1][1])
					{
						v[0] = -Q[1][1];
						v[1] = Q[1][0];
					}
					else
					{
						v[0] = 1;
						v[1] = 0;
					}

					var d:Number = v[0] * v[0] + v[1] * v[1];
					v[2] = - v[1] * s.y - v[0] * s.x;
					Q.add(_q.fromVectorMultiply(v)).scalar(1 / d);
				}

				var dx:Number = Math.abs(w.x - s.x);
				var dy:Number = Math.abs(w.y - s.y);
				if(dx <= .5 && dy <= .5)
				{
					vertex[i].x = w.x + p0.x;
					vertex[i].y = w.y + p0.y;
					continue;
				}

				// the minimum was not in the unit square; now minimize quadratic
				// on boundary of square
				var min:Number, cand:Number; // minimum and candidate for minimum of quad. form
				minCoord.x = s.x; minCoord.y = s.y; // coordinates of minimum
				min = Q.apply(s);

				if(Q[0][0] != 0.0)
				{
					for(z = 0; z < 2; z++) // value of the y-coordinate
					{
						w.y = s.y - 0.5 + z;
						w.x = - (Q[0][1] * w.y + Q[0][2]) / Q[0][0];
						dx = (w.x - s.x > 0 ? w.x - s.x : s.x - w.x);
						cand = Q.apply(w);
						if(dx <= .5 && cand < min)
						{
							min = cand;
							minCoord.x = w.x; minCoord.y = w.y;
						}
					}
				}

				if(Q[1][1] != 0.0)
				{
					for(z = 0; z < 2; z++)	 // value of the x-coordinate
					{
						w.x = s.x - 0.5 + z;
						w.y = - (Q[1][0] * w.x + Q[1][2]) / Q[1][1];
						dy = (w.y - s.y > 0 ? w.y - s.y : s.y - w.y);
						cand = Q.apply(w);
						if(dy <= .5 && cand < min)
						{
							min = cand;
							minCoord.x = w.x; minCoord.y = w.y;
						}
					}
				}

				// check four corners
				for(l = 0; l < 2; l++)
				{
					for(k = 0; k < 2; k++)
					{
						w.x = s.x - 0.5 + l;
						w.y = s.y - 0.5 + k;
						cand = Q.apply(w);
						if(cand < min)
						{
							min = cand;
							minCoord.x = w.x; minCoord.y = w.y;
						}
					}
				}

				vertex[i].x = minCoord.x + p0.x;
				vertex[i].y = minCoord.y + p0.y;
				continue;
			}

			return vertex;
		}

		/**
		 *  Stage 4: smoothing and corner analysis (Sec. 2.3.3)
		 *
		 *  Always succeeds and returns 0
		 */
		public static function smooth(vertex:Array, sign:String, alphamax:Number = 1.0):ClosedPath
		{
			var m:int = vertex.length;
			var closedPath:ClosedPath = new ClosedPath();

			var i:int;

			if(sign == '-')
			{
				// reverse orientation of negative paths
				var tmp:Point;
				for(i = 0, j = m - 1; i < j; i++, j--)
				{
					tmp = vertex[i];
					vertex[i] = vertex[j];
					vertex[j] = tmp;
				}
			}

			// examine each vertex and find its best fit
			var p2:Point = new Point();
			var p3:Point = new Point();
			var p4:Point = new Point();
			for(i = 0; i < m; i++)
			{
				var j:int = (i + 1) % m;
				var k:int = (i + 2) % m;
				interval(1 / 2.0, vertex[k], vertex[j], p4);

				var curve:Curve = new Curve();
				curve.vertex.x =vertex[j].x;
				curve.vertex.y =vertex[j].y;

				var denom:Number = ddenom(vertex[i], vertex[k]);
				var alpha:Number;
				if(denom != 0.0)
				{
					var dd:Number = dpara(vertex[i], vertex[j], vertex[k]) / denom;
					dd = dd < 0 ? -dd : dd;
					alpha = dd > 1 ? (1 - 1.0 / dd) / 0.75 : 0;
				}
				else
				{
					alpha = 4 / 3.0;
				}
				curve.alpha0 = alpha;   // remember "original" value of alpha

				if(alpha > alphamax)  // pointed corner
				{
					curve.tag  = POTRACE_CORNER;
					curve.c[0].x = 0; curve.c[0].y = 0;
					curve.c[1].x = vertex[j].x; curve.c[1].y = vertex[j].y;
					curve.c[2].x = p4.x; curve.c[2].y = p4.y;
				}
				else
				{
					if(alpha < 0.55)
					{
						alpha = 0.55;
					}
					else if(alpha > 1)
					{
						alpha = 1;
					}
					interval(.5 + .5 * alpha, vertex[i], vertex[j], p2);
					interval(.5 + .5 * alpha, vertex[k], vertex[j], p3);
					curve.tag = POTRACE_CURVETO;
					curve.c[0].x = p2.x; curve.c[0].y = p2.y;
					curve.c[1].x = p3.x; curve.c[1].y = p3.y;
					curve.c[2].x = p4.x; curve.c[2].y = p4.y;
				}
				curve.alpha = alpha; // store the "cropped" value of alpha
				curve.beta  = 0.5;

				closedPath.$a[j] = curve;
			}

			return closedPath;
		}

		/**
		 *  determine the center and slope of the line i..j. Assume i<j. Needs
		 *	"sum" components of p to be set.
		 */
		public static function pointslope(pt:Array, sums:Array, i:int, j:int, ctr:Point, dir:Point):void
		{
			// assume i<j

			var n:int = pt.length;

			var x:Number, y:Number, x2:Number, xy:Number, y2:Number;
			var k:Number;
			var a:Number, b:Number, c:Number, lambda2:Number, l:Number;
			var r:int = 0; // rotations from i to j

			while (j >= n)
			{
				j -= n;
				r += 1;
			}
			while(i >= n)
			{
				i -= n;
				r -= 1;
			}
			while(j < 0)
			{
				j += n;
				r -= 1;
			}
			while(i < 0)
			{
				i += n;
				r += 1;
			}
			
			x  = sums[j + 1].x  - sums[i].x  + r * sums[n].x;
			y  = sums[j + 1].y  - sums[i].y  + r * sums[n].y;
			x2 = sums[j + 1].x2 - sums[i].x2 + r * sums[n].x2;
			xy = sums[j + 1].xy - sums[i].xy + r * sums[n].xy;
			y2 = sums[j + 1].y2 - sums[i].y2 + r * sums[n].y2;
			k = j + 1 - i + r * n;
			
			ctr.x = x / k;
			ctr.y = y / k;

			a = (x2 - x * x / k) / k;
			b = (xy - x * y / k) / k;
			c = (y2 - y * y / k) / k;
			
			lambda2 = (a + c + Math.sqrt((a - c)*(a - c) + 4 * b * b)) / 2; // larger e.value 

			// now find e.vector for lambda2
			a -= lambda2;
			c -= lambda2;

			if(Math.abs(a) >= Math.abs(c))
			{
				l = Math.sqrt(a * a + b * b);
				if(l != 0)
				{
					dir.x = -b / l;
					dir.y =  a / l;
				}
			}
			else
			{
				l = Math.sqrt(c * c + b * b);
				if (l != 0)
				{
					dir.x = -c / l;
					dir.y =  b / l;
				}
			}
			if(l == 0)
			{
				dir.x = dir.y = 0;  // sometimes this can happen when k=4:
				                      // the two eigenvalues coincide
			}
		}

		/**
		 *  range over the straight line segment [a,b] when lambda ranges over [0,1]
		 */
		public static function interval(lambda:Number, a:Point, b:Point, ret:Point):void
		{
			ret.x = a.x + lambda * (b.x - a.x);
			ret.y = a.y + lambda * (b.y - a.y);
		}

		/**
		 *  some useful macros. Note: the "mod" macro works correctly for
		 *  negative a. Also note that the test for a>=n, while redundant,
		 *  speeds up the mod function by 70% in the average case (significant
		 *  since the program spends about 16% of its time here - or 40%
		 *  without the test). The "floordiv" macro returns the largest integer
		 *  <= a/n, and again this works correctly for negative a, as long as
		 *  a,n are integers and n>0.
		 */
		private static function mod(a:int, n:int):int
		{
			return a >= n ? a % n : a >= 0 ? a : n - 1 - (-1 - a) % n;
		}

		private static function floordiv(a:int, n:int):int
		{
			return a >= 0 ? a / n : -1 - (-1 - a) / n;
		}

		/**
		 *  calculate p1 x p2
		 */
		private static function xprod(p1:Point, p2:Point):int
		{
			return p1.x * p2.y - p1.y * p2.x;
		}

		/**
		 *  sign
		 */
		private static function sign(x:Number):int
		{
			return (x > 0 ? 1 : x < 0 ? -1 : 0);
		}

		/**
		 *  return a direction that is 90 degrees counterclockwise from p2-p0,
		 *  but then restricted to one of the major wind directions (n, nw, w, etc)
		 */
		public static function dorth_infty(p0:Point, p2:Point):Point
		{
			return new Point(
				sign(p2.x  - p0.x), 
				-sign(p2.y - p0.y)
				);
		}

		/**
		 *  return (p1-p0)x(p2-p0), the area of the parallelogram
		 */
		static public function dpara(p0:Point, p1:Point, p2:Point):Number
		{
			var x1:Number = p1.x - p0.x;
			var y1:Number = p1.y - p0.y;
			var x2:Number = p2.x - p0.x;
			var y2:Number = p2.y - p0.y;

			return x1 * y2 - x2 * y1;
		}

		/**
		 *  ddenom/dpara have the property that the square of radius 1 centered
		 *  at p1 intersects the line p0p2 iff |dpara(p0,p1,p2)| <= ddenom(p0,p2)
		 */
		public static function ddenom(p0:Point, p2:Point):Number
		{
			var r:Point = dorth_infty(p0, p2);
			return r.y * (p2.x - p0.x) - r.x * (p2.y - p0.y);
		}

		/**
		 *  return 1 if a <= b < c < a, in a cyclic sense (mod n)
		 */
		private static function cyclic(a:int, b:int, c:int):Boolean
		{
			if (a <= c)
			{
				return (a <= b) && (b < c);
			}
			else
			{
				return (a <= b) || (b < c);
			}
		}
	}
}


internal class Sums
{
	public var x:Number;
	public var y:Number;
	public var x2:Number;
	public var xy:Number;
	public var y2:Number;
}

import flash.geom.Point;
import flash.utils.Proxy;
import flash.utils.flash_proxy;
import flash.errors.IllegalOperationError;

internal class Point3D extends Proxy
{
	private var $v:Array;

	public function Point3D(x:Number = 0.0, y:Number = 0.0, z:Number = 0.0):void
	{
		$v = [x, y, z];
	}

	flash_proxy override function getProperty(name:*):*
	{
		if(name == 0 || name == "x") return $v[0];
		if(name == 1 || name == "y") return $v[1];
		if(name == 2 || name == "z") return $v[2];
		return undefined;
	}

	flash_proxy override function hasProperty(name:*):Boolean
	{
		return flash_proxy.getProperty(name) != undefined;
	}

	flash_proxy override function setProperty(name:*, value:*):void
	{
		if(name == 0 || name == "x"){$v[0] = Number(value); return;}
		if(name == 1 || name == "y"){$v[1] = Number(value); return;}
		if(name == 2 || name == "z"){$v[2] = Number(value); return;}
		throw new IllegalOperationError();
	}

	public function toString():String
	{
		return "(" + $v[0] + "," + $v[1] + "," + $v[2] + ")";
	}
}

/**
 *  the type of (affine) quadratic forms, represented as symmetric 3x3
 *  matrices.  The value of the quadratic form at a vector (x,y) is v^t
 *  Q v, where v = (x,y,1)^t.
 */
internal class QuadraticForm extends Proxy
{
	private var $m:Array;

	public function QuadraticForm():void
	{
		$m = [];
		$m[0] = new Point3D();
		$m[1] = new Point3D();
		$m[2] = new Point3D();
	}

	flash_proxy override function getProperty(name:*):*
	{
		if(name == 0) return $m[0];
		if(name == 1) return $m[1];
		if(name == 2) return $m[2];
		return undefined;
	}

	flash_proxy override function hasProperty(name:*):Boolean
	{
		return flash_proxy.getProperty(name) != undefined;
	}

	flash_proxy override function setProperty(name:*, value:*):void
	{
		throw new IllegalOperationError();
	}

	/**
	 *  Apply quadratic form Q to vector w = (w.x,w.y)
	 */
	public function apply(w:Point):Number
	{
		var v:Array = [w.x, w.y, 1];
		var sum:Number = 0.0;

		for(var i:int = 0; i < 3; i++)
		{
			for(var j:int = 0; j < 3; j++)
			{
				sum += v[i] * $m[i][j] * v[j];
			}
		}
		return sum;
	}

	public function clone():QuadraticForm
	{
		var ret:QuadraticForm = new QuadraticForm();
		for(var i:int = 0; i < 3; i++)
		{
			for(var j:int = 0; j < 3; j++)
			{
				ret[i][j] = $m[i][j];
			}
		}
		return ret;
	}

	public function add(m2:QuadraticForm):QuadraticForm
	{
		for(var i:int = 0; i < 3; i++)
		{
			for(var j:int = 0; j < 3; j++)
			{
				$m[i][j] += m2[i][j];
			}
		}
		return this;
	}

	public function scalar(s:Number):QuadraticForm
	{
		for(var i:int = 0; i < 3; i++)
		{
			for(var j:int = 0; j < 3; j++)
			{
				$m[i][j] *= s;
			}
		}
		return this;
	}

	public function fromVectorMultiply(v:Point3D):QuadraticForm
	{
		for(var i:int = 0; i < 3; i++)
		{
			for(var j:int = 0; j < 3; j++)
			{
				$m[i][j] = v[i] * v[j];
			}
		}
		return this;
	}

	public function toString():String
	{
		return "[" + 
			$m[0][0] + "," + $m[0][1] + "," + $m[0][2] + "," + 
			$m[1][0] + "," + $m[1][1] + "," + $m[1][2] + "," + 
			$m[2][0] + "," + $m[2][1] + "," + $m[2][2] + "]";
	}
}