#ifndef GEOM_H
#define GEOM_H

#include <math.h>

class Vector;
class Line;
class Plane;
class Matrix;
class Transform;

static float PI = 3.14159265359f;		//180 degrees
static float TWOPI = PI*2.0f;			//360 degrees
static float HALFPI = PI*.5f;			//90  degrees
static float QUARTERPI = PI*.25f;		//45  degrees

class Vector {
public:
	float x, y, z;

	Vector() :x(0), y(0), z(0) {
	}
	Vector(float x, float y, float z) :x(x), y(y), z(z) {
	}
	operator float*() {
		return &x;
	}
	operator const float *() {
		return &x;
	}
	float &operator[](int n) {
		return (&x)[n];
	}
	float operator[](int n)const {
		return (&x)[n];
	}
	Vector operator-()const {
		return Vector(-x, -y, -z);
	}
	Vector operator*(float scale)const {
		return Vector(x*scale, y*scale, z*scale);
	}
	Vector operator*(const Vector &q)const {
		return Vector(x*q.x, y*q.y, z*q.z);
	}
	Vector operator/(float scale)const {
		return Vector(x / scale, y / scale, z / scale);
	}
	Vector operator/(const Vector &q)const {
		return Vector(x / q.x, y / q.y, z / q.z);
	}
	Vector operator+(const Vector &q)const {
		return Vector(x + q.x, y + q.y, z + q.z);
	}
	Vector operator-(const Vector &q)const {
		return Vector(x - q.x, y - q.y, z - q.z);
	}
	Vector &operator*=(float scale) {
		x *= scale; y *= scale; z *= scale; return *this;
	}
	Vector &operator*=(const Vector &q) {
		x *= q.x; y *= q.y; z *= q.z; return *this;
	}
	Vector &operator/=(float scale) {
		x /= scale; y /= scale; z /= scale; return *this;
	}
	Vector &operator/=(const Vector &q) {
		x /= q.x; y /= q.y; z /= q.z; return *this;
	}
	Vector &operator+=(const Vector &q) {
		x += q.x; y += q.y; z += q.z; return *this;
	}
	Vector &operator-=(const Vector &q) {
		x -= q.x; y -= q.y; z -= q.z; return *this;
	}
	bool operator<(const Vector &q)const {
		if (fabs(x - q.x) > FLT_EPSILON) return x < q.x ? true : false;
		if (fabs(y - q.y) > FLT_EPSILON) return y < q.y ? true : false;
		return fabs(z - q.z) > FLT_EPSILON && z < q.z;
	}
	bool operator==(const Vector &q)const {
		return fabs(x - q.x) <= FLT_EPSILON && fabs(y - q.y) <= FLT_EPSILON && fabs(z - q.z) <= FLT_EPSILON;
	}
	bool operator!=(const Vector &q)const {
		return fabs(x - q.x) > FLT_EPSILON || fabs(y - q.y) > FLT_EPSILON || fabs(z - q.z) > FLT_EPSILON;
	}
	float dot(const Vector &q)const {
		return x*q.x + y*q.y + z*q.z;
	}
	Vector cross(const Vector &q)const {
		return Vector(y*q.z - z*q.y, z*q.x - x*q.z, x*q.y - y*q.x);
	}
	float length()const {
		return sqrtf(x*x + y*y + z*z);
	}
	float distance(const Vector &q)const {
		float dx = x - q.x, dy = y - q.y, dz = z - q.z; return sqrtf(dx*dx + dy*dy + dz*dz);
	}
	Vector normalized()const {
		float l = length(); return Vector(x / l, y / l, z / l);
	}
	void normalize() {
		float l = length(); x /= l; y /= l; z /= l;
	}
	float yaw()const {
		return -atan2f(x, z);
	}
	float pitch()const {
		return -atan2f(y, sqrtf(x*x + z*z));
	}
	void clear() {
		x = y = z = 0;
	}
};

class Line {
public:
	Vector o, d;
	Line() {
	}
	Line(const Vector &o, const Vector &d) :o(o), d(d) {
	}
	Line operator+(const Vector &q)const {
		return Line(o + q, d);
	}
	Line operator-(const Vector &q)const {
		return Line(o - q, d);
	}
	Vector operator*(float q)const {
		return o + d*q;
	}
	Vector nearest(const Vector &q)const {
		return o + d*(d.dot(q - o) / d.dot(d));
	}
};

class Plane {
public:
	Vector n;
	float d;

	Plane() :d(0) {
	}
	//normal/offset form
	Plane(const Vector &n, float d) :n(n), d(d) {
	}
	//point/normal form
	Plane(const Vector &p, const Vector &n) :n(n), d(-n.dot(p)) {
	}
	//create plane from tri
	Plane(const Vector &v0, const Vector &v1, const Vector &v2) {
		n = (v1 - v0).cross(v2 - v0).normalized(); d = -n.dot(v0);
	}
	Plane operator-()const {
		return Plane(-n, -d);
	}
	float t_intersect(const Line &q)const {
		return -distance(q.o) / n.dot(q.d);
	}
	Vector intersect(const Line &q)const {
		return q*t_intersect(q);
	}
	Line intersect(const Plane &q)const {
		Vector lv = n.cross(q.n).normalized();
		return Line(q.intersect(Line(nearest(n*-d), n.cross(lv))), lv);
	}
	Vector nearest(const Vector &q)const {
		return q - n*distance(q);
	}
	void negate() {
		n = -n; d = -d;
	}
	float distance(const Vector &q)const {
		return n.dot(q) + d;
	}
};

struct Quat {
	float w;
	Vector v;
	Quat() :w(1) {
	}
	Quat(float w, const Vector &v) :w(w), v(v) {
	}
	Quat operator-()const {
		return Quat(w, -v);
	}
	Quat operator+(const Quat &q)const {
		return Quat(w + q.w, v + q.v);
	}
	Quat operator-(const Quat &q)const {
		return Quat(w - q.w, v - q.v);
	}
	Quat operator*(const Quat &q)const {
		return Quat(w*q.w - v.dot(q.v), q.v.cross(v) + q.v*w + v*q.w);
	}
	Vector operator*(const Vector &q)const {
		return (*this * Quat(0, q) * -*this).v;
	}
	Quat operator*(float q)const {
		return Quat(w*q, v*q);
	}
	Quat operator/(float q)const {
		return Quat(w / q, v / q);
	}
	float dot(const Quat &q)const {
		return v.x*q.v.x + v.y*q.v.y + v.z*q.v.z + w*q.w;
	}
	float length()const {
		return sqrtf(w*w + v.x*v.x + v.y*v.y + v.z*v.z);
	}
	void normalize() {
		*this = *this / length();
	}
	Quat normalized()const {
		return *this / length();
	}
	Quat slerpTo(const Quat &q, float a)const {
		Quat t = q;
		float d = dot(q), b = 1 - a;
		if (d < 0) { t.w = -t.w; t.v = -t.v; d = -d; }
		if (d < 1 - FLT_EPSILON) {
			float om = acosf(d);
			float si = sinf(om);
			a = sinf(a*om) / si;
			b = sinf(b*om) / si;
		}
		return *this*b + t*a;
	}
	Vector i()const {
		float xz = v.x*v.z, wy = w*v.y;
		float xy = v.x*v.y, wz = w*v.z;
		float yy = v.y*v.y, zz = v.z*v.z;
		return Vector(1 - 2 * (yy + zz), 2 * (xy - wz), 2 * (xz + wy));
	}
	Vector j()const {
		float yz = v.y*v.z, wx = w*v.x;
		float xy = v.x*v.y, wz = w*v.z;
		float xx = v.x*v.x, zz = v.z*v.z;
		return Vector(2 * (xy + wz), 1 - 2 * (xx + zz), 2 * (yz - wx));
	}
	Vector k()const {
		float xz = v.x*v.z, wy = w*v.y;
		float yz = v.y*v.z, wx = w*v.x;
		float xx = v.x*v.x, yy = v.y*v.y;
		return Vector(2 * (xz - wy), 2 * (yz + wx), 1 - 2 * (xx + yy));
	}
};

class Matrix {
	static Matrix tmps[64];
	static Matrix &alloc_tmp() { static int tmp = 0; return tmps[tmp++ & 63]; }
	friend class Transform;
public:
	Vector i, j, k;

	Matrix() :i(Vector(1, 0, 0)), j(Vector(0, 1, 0)), k(Vector(0, 0, 1)) {
	}
	Matrix(const Vector &i, const Vector &j, const Vector &k) :i(i), j(j), k(k) {
	}
	Matrix(const Quat &q) {
		float xx = q.v.x*q.v.x, yy = q.v.y*q.v.y, zz = q.v.z*q.v.z;
		float xy = q.v.x*q.v.y, xz = q.v.x*q.v.z, yz = q.v.y*q.v.z;
		float wx = q.w*q.v.x, wy = q.w*q.v.y, wz = q.w*q.v.z;
		i = Vector(1 - 2 * (yy + zz), 2 * (xy - wz), 2 * (xz + wy)),
			j = Vector(2 * (xy + wz), 1 - 2 * (xx + zz), 2 * (yz - wx)),
			k = Vector(2 * (xz - wy), 2 * (yz + wx), 1 - 2 * (xx + yy));
	}
	Matrix(float angle, const Vector &axis) {
		const Vector &u = axis;
		float c = cosf(angle), s = sinf(angle);
		float x2 = axis.x*axis.x, y2 = axis.y*axis.y, z2 = axis.z*axis.z;
		i = Vector(x2 + c*(1 - x2), u.x*u.y*(1 - c) - u.z*s, u.z*u.x*(1 - c) + u.y*s);
		j = Vector(u.x*u.y*(1 - c) + u.z*s, y2 + c*(1 - y2), u.y*u.z*(1 - c) - u.x*s);
		k = Vector(u.z*u.x*(1 - c) - u.y*s, u.y*u.z*(1 - c) + u.x*s, z2 + c*(1 - z2));
	}
	Vector &operator[](int n) {
		return (&i)[n];
	}
	const Vector &operator[](int n)const {
		return (&i)[n];
	}
	Matrix &operator~()const {
		Matrix &m = alloc_tmp();
		m.i.x = i.x; m.i.y = j.x; m.i.z = k.x;
		m.j.x = i.y; m.j.y = j.y; m.j.z = k.y;
		m.k.x = i.z; m.k.y = j.z; m.k.z = k.z;
		return m;
	}
	float determinant()const {
		return i.x*(j.y*k.z - j.z*k.y) - i.y*(j.x*k.z - j.z*k.x) + i.z*(j.x*k.y - j.y*k.x);
	}
	Matrix &operator-()const {
		Matrix &m = alloc_tmp();
		float t = 1.0f / determinant();
		m.i.x = t*(j.y*k.z - j.z*k.y); m.i.y = -t*(i.y*k.z - i.z*k.y); m.i.z = t*(i.y*j.z - i.z*j.y);
		m.j.x = -t*(j.x*k.z - j.z*k.x); m.j.y = t*(i.x*k.z - i.z*k.x); m.j.z = -t*(i.x*j.z - i.z*j.x);
		m.k.x = t*(j.x*k.y - j.y*k.x); m.k.y = -t*(i.x*k.y - i.y*k.x); m.k.z = t*(i.x*j.y - i.y*j.x);
		return m;
	}
	Matrix &cofactor()const {
		Matrix &m = alloc_tmp();
		m.i.x = (j.y*k.z - j.z*k.y); m.i.y = -(j.x*k.z - j.z*k.x); m.i.z = (j.x*k.y - j.y*k.x);
		m.j.x = -(i.y*k.z - i.z*k.y); m.j.y = (i.x*k.z - i.z*k.x); m.j.z = -(i.x*k.y - i.y*k.x);
		m.k.x = (i.y*j.z - i.z*j.y); m.k.y = -(i.x*j.z - i.z*j.x); m.k.z = (i.x*j.y - i.y*j.x);
		return m;
	}
	bool operator==(const Matrix &q)const {
		return i == q.i && j == q.j && k == q.k;
	}
	bool operator!=(const Matrix &q)const {
		return i != q.i || j != q.j || k != q.k;
	}
	Vector operator*(const Vector &q)const {
		return Vector(i.x*q.x + j.x*q.y + k.x*q.z, i.y*q.x + j.y*q.y + k.y*q.z, i.z*q.x + j.z*q.y + k.z*q.z);
	}
	Matrix &operator*(const Matrix &q)const {
		Matrix &m = alloc_tmp();
		m.i.x = i.x*q.i.x + j.x*q.i.y + k.x*q.i.z; m.i.y = i.y*q.i.x + j.y*q.i.y + k.y*q.i.z; m.i.z = i.z*q.i.x + j.z*q.i.y + k.z*q.i.z;
		m.j.x = i.x*q.j.x + j.x*q.j.y + k.x*q.j.z; m.j.y = i.y*q.j.x + j.y*q.j.y + k.y*q.j.z; m.j.z = i.z*q.j.x + j.z*q.j.y + k.z*q.j.z;
		m.k.x = i.x*q.k.x + j.x*q.k.y + k.x*q.k.z; m.k.y = i.y*q.k.x + j.y*q.k.y + k.y*q.k.z; m.k.z = i.z*q.k.x + j.z*q.k.y + k.z*q.k.z;
		return m;
	}
	void orthogonalize() {
		k.normalize();
		i = j.cross(k).normalized();
		j = k.cross(i);
	}
	Matrix &orthogonalized()const {
		Matrix &m = alloc_tmp();
		m = *this; m.orthogonalize();
		return m;
	}
};

class Box {
public:
	Vector a, b;
	Box() :a(Vector(INFINITY, INFINITY, INFINITY)), b(Vector(-INFINITY, -INFINITY, -INFINITY)) {
	}
	Box(const Vector &q) :a(q), b(q) {
	}
	Box(const Vector &a, const Vector &b) :a(a), b(b) {
	}
	Box(const Line &l) :a(l.o), b(l.o) {
		update(l.o + l.d);
	}
	void clear() {
		a.x = a.y = a.z = INFINITY;
		b.x = b.y = b.z = -INFINITY;
	}
	bool empty()const {
		return b.x < a.x || b.y < a.y || b.z < a.z;
	}
	Vector centre()const {
		return Vector((a.x + b.x)*.5f, (a.y + b.y)*.5f, (a.z + b.z)*.5f);
	}
	Vector corner(int n)const {
		return Vector(((n & 1) ? b : a).x, ((n & 2) ? b : a).y, ((n & 4) ? b : a).z);
	}
	void update(const Vector &q) {
		if (q.x < a.x) a.x = q.x; if (q.y < a.y) a.y = q.y; if (q.z < a.z) a.z = q.z;
		if (q.x > b.x) b.x = q.x; if (q.y > b.y) b.y = q.y; if (q.z > b.z) b.z = q.z;
	}
	void update(const Box &q) {
		if (q.a.x < a.x) a.x = q.a.x; if (q.a.y < a.y) a.y = q.a.y; if (q.a.z < a.z) a.z = q.a.z;
		if (q.b.x > b.x) b.x = q.b.x; if (q.b.y > b.y) b.y = q.b.y; if (q.b.z > b.z) b.z = q.b.z;
	}
	bool overlaps(const Box &q)const {
		return
			(b.x < q.b.x ? b.x : q.b.x) >= (a.x > q.a.x ? a.x : q.a.x) &&
			(b.y < q.b.y ? b.y : q.b.y) >= (a.y > q.a.y ? a.y : q.a.y) &&
			(b.z < q.b.z ? b.z : q.b.z) >= (a.z > q.a.z ? a.z : q.a.z);
	}
	void expand(float n) {
		a.x -= n; a.y -= n; a.z -= n; b.x += n; b.y += n; b.z += n;
	}
	float width()const {
		return b.x - a.x;
	}
	float height()const {
		return b.y - a.y;
	}
	float depth()const {
		return b.z - a.z;
	}
	bool contains(const Vector &q) {
		return q.x >= a.x && q.x <= b.x && q.y >= a.y && q.y <= b.y && q.z >= a.z && q.z <= b.z;
	}
};

class Transform {
	static Transform tmps[64];
	static Transform &alloc_tmp() { static int tmp = 0; return tmps[tmp++ & 63]; }
public:
	Matrix m;
	Vector v;

	Transform() {
	}
	Transform(const Matrix &m) :m(m) {
	}
	Transform(const Vector &v) :v(v) {
	}
	Transform(const Matrix &m, const Vector &v) :m(m), v(v) {
	}
	Transform &operator-()const {
		Transform &t = alloc_tmp();
		t.m = -m; t.v = t.m*-v;
		return t;
	}
	Transform &operator~()const {
		Transform &t = alloc_tmp();
		t.m = ~m; t.v = t.m*-v;
		return t;
	}
	Vector operator*(const Vector &q)const {
		return m*q + v;
	}
	Line operator*(const Line &q)const {
		Vector t = (*this)*q.o;
		return Line(t, (*this)*(q.o + q.d) - t);
	}
	Box operator*(const Box &q)const {
		Box t((*this*q.corner(0)));
		for (int k = 1; k < 8; ++k) t.update(*this*q.corner(k));
		return t;
	}
	Transform &operator*(const Transform &q)const {
		Transform &t = alloc_tmp();
		t.m = m*q.m; t.v = m*q.v + v;
		return t;
	}
	bool operator==(const Transform &q)const {
		return m == q.m && v == q.v;
	}
	bool operator!=(const Transform &q)const {
		return !operator==(q);
	}
};

inline float transformRadius(float r, const Matrix &t) {
	static const float sq_3 = sqrtf(1.0f / 3.0f);
	return (t * Vector(sq_3, sq_3, sq_3)).length()*r;
}

inline Matrix pitchMatrix(float q) {
	return Matrix(Vector(1, 0, 0), Vector(0, cosf(q), sinf(q)), Vector(0, -sinf(q), cosf(q)));
}

inline Matrix yawMatrix(float q) {
	return Matrix(Vector(cosf(q), 0, sinf(q)), Vector(0, 1, 0), Vector(-sinf(q), 0, cosf(q)));
}

inline Matrix rollMatrix(float q) {
	return Matrix(Vector(cosf(q), sinf(q), 0), Vector(-sinf(q), cosf(q), 0), Vector(0, 0, 1));
}

inline float matrixPitch(const Matrix &m) {
	return m.k.pitch();
	//	return asinf( -m.k.y );
}

inline float matrixYaw(const Matrix &m) {
	return m.k.yaw();
	//return atan2f( -m.k.x,m.k.z );
}

inline float matrixRoll(const Matrix &m) {
	return atan2f(m.i.y, m.j.y);
	//Matrix t=pitchMatrix( -matrixPitch(m) )*yawMatrix( -matrixYaw(m) )*m;
	//return atan2f( t.i.y,t.i.x );
}

inline Matrix scaleMatrix(float x, float y, float z) {
	return Matrix(Vector(x, 0, 0), Vector(0, y, 0), Vector(0, 0, z));
}

inline Matrix scaleMatrix(const Vector &scale) {
	return Matrix(Vector(scale.x, 0, 0), Vector(0, scale.y, 0), Vector(0, 0, scale.z));
}

inline Quat pitchQuat(float p) {
	return Quat(cosf(p / -2), Vector(sinf(p / -2), 0, 0));
}

inline Quat yawQuat(float y) {
	return Quat(cosf(y / 2), Vector(0, sinf(y / 2), 0));
}

inline Quat rollQuat(float r) {
	return Quat(cosf(r / -2), Vector(0, 0, sinf(r / -2)));
}

//inline Quat rotationQuat( float p,float y,float r ){
//	return yawQuat(y)*pitchQuat(p)*rollQuat(r);
//}

Quat rotationQuat(float p, float y, float r);

inline Matrix rotationMatrix(float p, float y, float r) {
	return yawMatrix(y)*pitchMatrix(p)*rollMatrix(r);
}

inline Matrix rotationMatrix(const Vector &rot) {
	return yawMatrix(rot.y)*pitchMatrix(rot.x)*rollMatrix(rot.z);
}

inline float quatPitch(const Quat &q) {
	return q.k().pitch();
}

inline float quatYaw(const Quat &q) {
	return q.k().yaw();
}

inline float quatRoll(const Quat &q) {
	//	Vector i=q.i(),j=q.j();
	//	return atan2f( i.y,j.y );
	return matrixRoll(q);
}

inline Quat matrixQuat(const Matrix &p) {
	Matrix m = p;
	m.orthogonalize();
	float t = m.i.x + m.j.y + m.k.z, w, x, y, z;
	if (t > FLT_EPSILON) {
		t = sqrtf(t + 1) * 2;
		x = (m.k.y - m.j.z) / t;
		y = (m.i.z - m.k.x) / t;
		z = (m.j.x - m.i.y) / t;
		w = t / 4;
	} else if (m.i.x > m.j.y && m.i.x > m.k.z) {
		t = sqrtf(m.i.x - m.j.y - m.k.z + 1) * 2;
		x = t / 4;
		y = (m.j.x + m.i.y) / t;
		z = (m.i.z + m.k.x) / t;
		w = (m.k.y - m.j.z) / t;
	} else if (m.j.y > m.k.z) {
		t = sqrtf(m.j.y - m.k.z - m.i.x + 1) * 2;
		x = (m.j.x + m.i.y) / t;
		y = t / 4;
		z = (m.k.y + m.j.z) / t;
		w = (m.i.z - m.k.x) / t;
	} else {
		t = sqrtf(m.k.z - m.j.y - m.i.x + 1) * 2;
		x = (m.i.z + m.k.x) / t;
		y = (m.k.y + m.j.z) / t;
		z = t / 4;
		w = (m.j.x - m.i.y) / t;
	}
	return Quat(w, Vector(x, y, z));
}

#endif