using System;

namespace EllipseTruncate {
	public struct AngleRange {
		public double angle, size;
		
		public AngleRange(double angle, double size)
			{ this.angle = angle; this.size = size; }
		public bool isEmpty()
			{ return size <= Geometry.precision; }
		public bool isFull()
			{ return size >= 2.0*Math.PI - Geometry.precision; }
		public double end()
			{ return wrap(angle + size); }
		public bool contains(double a)
			{ return diff(a, angle) < size; }
		public bool intersects(AngleRange r)
			{ return contains(r.angle) || r.contains(angle); }
		public AngleRange union(AngleRange r) {
			double da = diff(r.angle, angle);
			if (da <= size + Geometry.precision)
				return new AngleRange(angle, Math.Max(size, da + r.size));
			da = 2.0*Math.PI - da;
			if (da <= r.size + Geometry.precision)
				return new AngleRange(r.angle, Math.Max(r.size, da + size));
			da = Math.Min(angle, r.angle);
			return new AngleRange(da, Math.Max(angle + size, r.angle + r.size) - da);
		}

		public static AngleRange byAngles(double a0, double a1)
			{ return new AngleRange(a0, diff(a1, a0)); }
		public static double diff(double a1, double a0)
			{ return a0 < a1 ? a1 - a0 : a1 - a0 + 2.0*Math.PI; }
		public static double wrap(double angle) {
			return angle >  Math.PI ? angle - 2.0*Math.PI
			     : angle < -Math.PI ? angle + 2.0*Math.PI
			     : angle;
		}
	}

	public class Ellipse {
		public Point center;
		public Point p1;
		public Point p2;
		public double r1;
		public double r2;
		public double angle;
		
		public Matrix matrix;
		public Matrix matrixInv;
	
		public Ellipse(Point p0, Point p1, Point p2) {
			center = p0;
			this.p1 = p1;
			Point d = p1 - p0;
			r1 = d.len();
			angle = d.atan();
			Point dp = d.rotate90()/r1;
			r2 = Math.Abs(dp*(p2 - p0));
			this.p2 = p0 + dp*r2;
			
			matrix = Matrix.identity()
					.translate(center)
			        .rotate(angle)
			        .scale(r1, r2);
			matrixInv = matrix.invert();
		}
		
		public void drawFull(Cairo.Context context) {
    		int segments = 100;
    		double da = 2.0*Math.PI/segments;
    		double s = Math.Sin(da);
    		double c = Math.Cos(da);
    		Point r = new Point(1.0, 0.0);

			context.Save();
			context.Matrix = (new Matrix(context.Matrix)*matrix).toCairo();
			context.SetSourceRGBA(1.0, 0.0, 0.0, 0.1);
			context.Arc(0.0, 0.0, 1.0, 0.0, 2.0*Math.PI);
			context.ClosePath();
			context.Fill();
			context.Restore();

    		context.Save();
			context.LineWidth = 0.5;
			context.SetSourceRGBA(1.0, 0.0, 0.0, 0.5);
			Point p = matrix*r;
			for(int i = 0; i < segments; ++i) {
				r = new Point(r.x*c - r.y*s, r.y*c + r.x*s);
  				p = matrix*r;
				context.LineTo(p.x, p.y);
			}
			context.ClosePath();
			context.Stroke();
    		context.Restore();
		}
		
		private static void swap(ref double a, ref double b)
			{ double c = a; a = b; b = c; }
		private static void swap(ref int a, ref int b)
			{ int c = a; a = b; b = c; }

		private bool cutRange(ref int index, AngleRange[] ranges, double da, double h) {
			if (h <= Geometry.precision - 1.0) {
				return true;
			} else
			if (h >= 1.0 - Geometry.precision) {
				return false;
			} else {
				double a = Math.Asin(h);
				ranges[index].angle = AngleRange.wrap( a < 0.0
				                                     ? da - a - Math.PI
				                                     : da - a + Math.PI );
				ranges[index].size = Math.PI + a + a;
				int j = index;
				for(int i = 0, r = 0; i < index;) {
					if (r > 0) ranges[i] = ranges[i + r];
					if (ranges[i].intersects(ranges[j])) {
						if (j < i) { ranges[j] = ranges[i].union(ranges[j]); ++r; --index; }
						      else { ranges[i] = ranges[i].union(ranges[j]); j = i; ++i; }
						if (ranges[j].isFull())
							return false;
					} else ++i;
				}
				if (j == index) ++index;
			}
			return true;
		}
		
		public void drawTruncated(Cairo.Context context, Point b0, Point b1, Point b2) {
			Point dx = matrixInv.turn(b1 - b0);
			Point dy = matrixInv.turn(b2 - b0);
			Point nx = dx.rotate90().normalize();
			Point ny = dy.rotate90().normalize();
			Point o = matrixInv*b0;
			
			double ax = dx.atan();
			double ay = dy.atan();
			double aax = AngleRange.wrap(ax + Math.PI);
			double aay = AngleRange.wrap(ay + Math.PI);

			double sign = nx*dy;
			if (Math.Abs(sign) <= Geometry.precision) return;
			if (sign < 0.0) {
				double a;
				nx *= -1.0; ny *= -1.0;
				a = ax; ax = aax; aax = a;
				a = ay; ay = aay; aay = a;
			}

			context.Save();
			context.Matrix = (new Matrix(context.Matrix)*matrix).toCairo();
			context.SetSourceRGBA(1.0, 0.0, 0.0, 0.1);
			context.MoveTo(o.x, o.y);
			context.RelLineTo(dx.x, dx.y);
			context.RelLineTo(dy.x, dy.y);
			context.RelLineTo(-dx.x, -dx.y);
			context.ClosePath();
			context.Fill();
			context.Restore();
			
			// gather invisible ranges
			AngleRange[] cutRanges = new AngleRange[4];
			int count = 0;
			if ( !cutRange(ref count, cutRanges, ax , (o*nx))          ) return;
			if ( !cutRange(ref count, cutRanges, aax, -((o+dx+dy)*nx)) ) return;
			if ( !cutRange(ref count, cutRanges, ay , (o+dx)*ny)       ) return;
			if ( !cutRange(ref count, cutRanges, aay, -((o+dy)*ny))    ) return;
			
			// sort bounds
			for(int i = 0; i < count; ++i)
				for(int j = i+1; j < count; ++j)
					if (cutRanges[j].angle < cutRanges[i].angle)
						{ AngleRange r = cutRanges[i]; cutRanges[i] = cutRanges[j]; cutRanges[j] = r; }
						
			// invert bounds
			AngleRange[] ranges = new AngleRange[4];
			for(int i = 0; i < count; ++i)
				ranges[i] = AngleRange.byAngles(cutRanges[(i > 0 ? i : count) - 1].end(), cutRanges[i].angle);
			
			// dummy 
			if (count == 0)
				ranges[count++] = new AngleRange(0.0, 2.0*Math.PI);
			
			// draw
    		int segments = 100;
    		double da = 2.0*Math.PI/segments;
    		double s = Math.Sin(da);
    		double c = Math.Cos(da);

    		context.Save();
			context.LineWidth = 2.0;
			context.SetSourceRGBA(0.0, 0.0, 1.0, 1.0);
			for(int i = 0; i < count; ++i) {
				double angle = ranges[i].angle;
				double size = ranges[i].size;
    			int cnt = (int)Math.Floor(size/da);
    			Point r = new Point(1.0, 0.0).rotate(angle);
    			Point p = matrix*r;
    			context.MoveTo(p.x, p.y);
				for(int j = 0; j < cnt; ++j) {
					r = new Point(r.x*c - r.y*s, r.y*c + r.x*s);
	  				p = matrix*r;
					context.LineTo(p.x, p.y);
				}
    			r = new Point(1.0, 0.0).rotate(angle + size);
    			p = matrix*r;
    			context.LineTo(p.x, p.y);
			}
			context.Stroke();
			context.Restore();
		}
	}
}