#ifndef CURVE_H
#define CURVE_H


#include <cassert>

#include "vector.h"


class Bezier;

class Hermite {
public:
	Vector p0, p1, t0, t1;

	explicit Hermite(
		const Vector &p0 = Vector(),
		const Vector &p1 = Vector(),
		const Vector &t0 = Vector(),
		const Vector &t1 = Vector()
	):
		p0(p0), p1(p1), t0(t0), t1(t1) { }

	bool operator< (const Hermite &a) const {
		return p0 < a.p0 ? true : a.p0 < p0 ? false
			 : p1 < a.p1 ? true : a.p1 < p1 ? false
			 : t0 < a.t0 ? true : a.t0 < t0 ? false
			 : t1 < a.t1;
	}
	bool operator>(const Hermite &a) const { return a < *this; }
	bool operator==(const Hermite &a) const { return p0 == a.p0 && p1 == a.p1 && t0 == a.t0 && t1 == a.t1; }
	bool operator!=(const Hermite &a) const { return !(*this == a); }

	Vector p(const Real &l) const
		{ return p(l, p0, p1, t0, t1); }
	Vector t(const Real &l) const
		{ return t(l, p0, p1, t0, t1); }

	Vector operator() (const Real &l) const  { return p(l); }
	
	Vector d(Real l) const {
		Vector t = this->t(l);
		Real lensqr = t.lensqr();
		if (precision*precision < lensqr) return t/sqrt(lensqr);
		return (p(l + differential) - p(l - differential)).norm();
	}
	
	Vector pp0() const { return p0 + t0/3; }
	Vector pp1() const { return p1 - t1/3; }
	
	Hermite flip() const
		{ return Hermite(p1, p0, -t1, -t0); }
	
	Hermite sub(const Real &a, const Real &b) const {
		Real k = b - a;
		return Hermite(p(a), p(b), t(a)*k, t(b)*k);
	}
	
	void split(const Real &l, Hermite *a, Hermite *b) const {
		if (!a && !b) return;
		assert(a != b);
		Vector point = p(l);
		Vector tangent = t(l);
		Hermite h = *this;
		if (a) {
			a->p0 = h.p0;
			a->t0 = h.t0*l;
			a->p1 = point;
			a->t1 = tangent*l;
		}
		if (b && b != a) {
			Real ll = 1 - l;
			b->p0 = point;
			b->t0 = tangent*ll;
			b->p1 = h.p1;
			b->t1 = h.t1*ll;
		}
	}

	Rect bounds() const
		{ return Rect(p0).expand(p1).expand(pp0()).expand(pp1()); }

	int bends_x(Real *l) const
		{ return bends(l, p0.x, p1.x, t0.x, t1.x); }
	int bends_y(Real *l) const
		{ return bends(l, p0.y, p1.y, t0.y, t1.y); }

	Range bounds_accurate_x() const
		{ return bounds_accurate(p0.x, p1.x, t0.x, t1.x); }
	Range bounds_accurate_y() const
		{ return bounds_accurate(p0.y, p1.y, t0.y, t1.y); }
	Rect bounds_accurate() const
		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }

	int inflection_x(Real *l) const
		{ return inflection(l, p0.x, p1.x, t0.x, t1.x); }
	int inflection_y(Real *l) const
		{ return inflection(l, p0.y, p1.y, t0.y, t1.y); }

	int intersections_x(Real *l, const Real &x) const
		{ return intersections(l, x, p0.x, p1.x, t0.x, t1.x); }
	int intersections_y(Real *l, const Real &y) const
		{ return intersections(l, y, p0.y, p1.y, t0.y, t1.y); }
	
	inline Bezier get_bezier() const;
	
	static Hermite linear(const Vector &p, const Vector &t)
		{ return Hermite(p, p + t, t, t);  }
	static Hermite point(const Vector &p)
		{ return Hermite(p, p);  }

	static Real p(const Real &l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
		return p0*((( 2*l - 3)*l    )*l + 1)
			 + p1*(((-2*l + 3)*l    )*l    )
			 + t0*(((   l - 2)*l + 1)*l    )
			 + t1*(((   l - 1)*l    )*l    );
	}

	static Real t(const Real &l, const Real &p0, const Real &p1, const Real &t0, const Real &t1) {
		return p0*(( 6*l - 6)*l    )
			 + p1*((-6*l + 6)*l    )
			 + t0*(( 3*l - 4)*l + 1)
			 + t1*(( 3*l - 2)*l    );
	}

	static Vector p(const Real &l, const Vector &p0, const Vector &p1, const Vector &t0, const Vector &t1)
		{ return Vector( p(l, p0.x, p1.x, t0.x, t1.x), p(l, p0.y, p1.y, t0.y, t1.y) ); }
	static Vector t(const Real &l, const Vector &p0, const Vector &p1, const Vector &t0, const Vector &t1)
		{ return Vector( t(l, p0.x, p1.x, t0.x, t1.x), t(l, p0.y, p1.y, t0.y, t1.y) ); }
	
	static int inflection(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
	static int bends(Real *l, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
	static int intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &t0, const Real &t1);
	static Range bounds_accurate(const Real &p0, const Real &p1, const Real &t0, const Real &t1);
};


class Bezier {
public:
	Vector p0, pp0, pp1, p1;

	Bezier(const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1):
		p0(p0), pp0(pp0), pp1(pp1), p1(p1) { }
	Bezier(const Vector &p0, const Vector &p1, const Vector &pp):
		Bezier(p0, p1, pp, pp) { }
	Bezier(const Vector &p0, const Vector &p1):
		Bezier(p0, p1, (p1 - p0)/3 + p0, (p0 - p1)/3 + p1) { }
	explicit Bezier(const Vector &p = Vector()):
		Bezier(p, p, p, p) { }

	bool operator< (const Bezier &a) const {
		return p0  < a.p0  ? true : a.p0  < p0  ? false
			 : p1  < a.p1  ? true : a.p1  < p1  ? false
			 : pp0 < a.pp0 ? true : a.pp0 < pp0 ? false
			 : pp1 < a.pp1;
	}
	bool operator>(const Bezier &a) const { return a < *this; }
	bool operator==(const Bezier &a) const { return p0 == a.p0 && p1 == a.p1 && pp0 == a.pp0 && pp1 == a.pp1; }
	bool operator!=(const Bezier &a) const { return !(*this == a); }

	Vector p(const Real &l) const
		{ return p(l, p0, p1, pp0, pp1); }
	Vector t(const Real &l) const
		{ return t(l, p0, p1, pp0, pp1); }

	Vector operator() (const Real &l) const  { return p(l); }
	
	Vector d(Real l) const {
		Vector t = this->t(l);
		Real lensqr = t.lensqr();
		if (precision*precision < lensqr) return t/sqrt(lensqr);
		return (p(l + differential) - p(l - differential)).norm();
	}
	
	Vector t0() const { return 3*(pp0 - p0); }
	Vector t1() const { return 3*(p1 - pp1); }
	
	void split(const Real &l, Bezier *a, Bezier *b) const {
		if (!a && !b) return;
		assert(a != b);
		Real ll = 1 - l;
		Vector p0 = this->p0;
		Vector p1 = this->p1;
		Vector a0 =  p0*ll + pp0*l;
		Vector a1 = pp0*ll + pp1*l;
		Vector a2 = pp1*ll +  p1*l;
		Vector b0 =  a0*ll +  a1*l;
		Vector b1 =  a1*ll +  a2*l;
		Vector c  =  b0*ll +  b1*l;
		if (a) {
			a->p0  = p0;
			a->pp0 = a0;
			a->pp1 = b0;
			a->p1  = c;
		}
		if (b && b != a) {
			b->p0  = c;
			b->pp0 = b1;
			b->pp1 = a2;
			b->p1  = p1;
		}
	}
	
	Bezier flip() const
		{ return Bezier(p1, p0, pp1, pp0); }

	Bezier sub(const Real &a, const Real &b) const {
		if (equal(a, b)) return Bezier( p(0.5*(a + b)) );

		Bezier s;
		if (equal(a, 1)) {
			split(b, nullptr, &s);
			return s.flip();
		}
		
		split(a, nullptr, &s);
		s.split((b - a)/(1 - a), &s, nullptr);
		return s;
	}
	
	Rect bounds() const
		{ return Rect(p0).expand(p1).expand(pp0).expand(pp1); }

	int bends_x(Real *l) const
		{ return bends(l, p0.x, p1.x, pp0.x, pp1.x); }
	int bends_y(Real *l) const
		{ return bends(l, p0.y, p1.y, pp0.y, pp1.y); }

	Range bounds_accurate_x() const
		{ return bounds_accurate(p0.x, p1.x, pp0.x, pp1.x); }
	Range bounds_accurate_y() const
		{ return bounds_accurate(p0.y, p1.y, pp0.y, pp1.y); }
	Rect bounds_accurate() const
		{ return Rect(bounds_accurate_x(), bounds_accurate_y()); }

	int inflection_x(Real *l) const
		{ return inflection(l, p0.x, p1.x, pp0.x, pp1.x); }
	int inflection_y(Real *l) const
		{ return inflection(l, p0.y, p1.y, pp0.y, pp1.y); }

	int intersections_x(Real *l, const Real &x) const
		{ return intersections(l, x, p0.x, p1.x, pp0.x, pp1.x); }
	int intersections_y(Real *l, const Real &y) const
		{ return intersections(l, y, p0.y, p1.y, pp0.y, pp1.y); }

	inline Hermite get_hermite() const;

	static Bezier linear(const Vector &p0, const Vector &p1)
		{ return Bezier(p0, p1, p0 + (p1 - p0)/3, p1 + (p0 - p1)/3);  }
	static Bezier point(const Vector &p)
		{ return Bezier(p, p);  }

	static Real p(const Real &l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
		Real ll = 1 - l;
		const Real a0 =  p0*ll + pp0*l;
		const Real a1 = pp0*ll + pp1*l;
		const Real a2 = pp1*ll +  p1*l;
		const Real b0 =  a0*ll +  a1*l;
		const Real b1 =  a1*ll +  a2*l;
		return b0*ll +  b1*l;
	}

	static Real t(const Real &l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1) {
		Real ll = 1 - l;
		const Real a0 =  p0*ll + pp0*l;
		const Real a1 = pp0*ll + pp1*l;
		const Real a2 = pp1*ll +  p1*l;
		const Real b0 =  a0*ll +  a1*l;
		const Real b1 =  a1*ll +  a2*l;
		return 3*(b1 - b0)*l;
	}

	static Vector p(const Real &l, const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1)
		{ return Vector( p(l, p0.x, p1.x, pp0.x, pp1.x), p(l, p0.y, p1.y, pp0.y, pp1.y) ); }
	static Vector t(const Real &l, const Vector &p0, const Vector &p1, const Vector &pp0, const Vector &pp1)
		{ return Vector( t(l, p0.x, p1.x, pp0.x, pp1.x), t(l, p0.y, p1.y, pp0.y, pp1.y) ); }
	
	static int inflection(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
	static int bends(Real *l, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
	static int intersections(Real *l, const Real &p, const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
	static Range bounds_accurate(const Real &p0, const Real &p1, const Real &pp0, const Real &pp1);
};


inline Bezier Hermite::get_bezier() const { return Bezier(p0, p1, pp0(), pp1()); }
inline Hermite Bezier::get_hermite() const { return Hermite(p0, p1, t0(), t1()); }

#endif