#ifndef _POINT_H_
#define _POINT_H_

/* class Pt2, Pt3, and Pt<N>
 *
 * These classes represent a POINT in R^2, R^3, and R^N space,
 * respectively. There is NO way to cast a Pt into a Vec. Needing to do
 * so means you're thinking about it wrong. The way to do it would be
 * (ptN - PtN::origin) if you really have to.
 *
 * Things to note:
 *  -member e[] is the array of coordinates. Named components can be
 *   gotten from the component part of the union.
 *  -The components of a vector CANNOT be extracted by the [] operator.
 *   Doing so will invoke the OpenGL cast operator. Use the member e[]
 *   instead.
 *  -You can pass a PtN class straight to glVertexNfv().
 *
 * Affine geometry notes:
 *   Point +/- Vector = Point
 *   Vector + Vector  = Vector
 *   Point - Point    = Vector
 *   Point + Point    = Point (so that we can do affine combinations of points)
 *
 * Last modified: 2008-02-10
 */

#include <vDEC/include/vdec_types.h>
#include <vDEC/include/Vector.h>
#include <iostream>

struct Pt2{
	// The actual coordinate data
	union{
		real e[2];
		struct{
			real x, y;
		} component;
	};
	
	/// Returns the origin.
	inline static const Pt2 Origin();
	/// Returns a point whose components are very large positive numbers.
	inline static const Pt2 PositiveInfinity();
	/// Returns a point whose components are very large negative numbers.
	inline static const Pt2 NegativeInfinity();
	
	inline Pt2();
	inline Pt2(const real &a, const real &b);
	inline Pt2(const Pt2 &P);

	inline Pt2& operator +=(const Vec2 &V);
	inline Pt2& operator +=(const Pt2 &Q);
	inline Pt2& operator *=(const real &C);
	inline Pt2& operator -=(const Vec2 &V);
	inline Pt2& operator /=(const real &C);

	/// Cast into a raw array. Suitable for use in glVertex3().
	inline operator real*( );
	/// Cast into a raw array. Suitable for use in glVertex3().
	inline operator const real*( ) const;

	/// Sorts the two points in lexicographical order such that P0 < P1 returns -1, equality 0, else 1.
	inline static int Order(const Pt2 &P0, const Pt2 &P1);
};

inline bool operator<(const Pt2 &A, const Pt2 &B);
inline bool operator>(const Pt2 &A, const Pt2 &B);
inline Vec2 operator-(const Pt2 &A, const Pt2 &B);
inline Pt2 operator-(const Pt2 &P, const Vec2 &V);
inline Pt2 operator-(const Pt2 &P);
inline Pt2 operator+(const Pt2 &P, const Vec2 &D);
inline Pt2 operator+(const Pt2 &P, const Pt2 &Q);
inline Pt2 operator*(const real &C, const Pt2 &P);

struct Pt3{
	union{
		real e[3];
		struct{
			real x, y, z;
		} components;
	};

	/// Returns the origin.
	inline static const Pt3 Origin();
	/// Returns a point whose components are very large positive numbers.
	inline static const Pt3 PositiveInfinity();
	/// Returns a point whose components are very large negative numbers.
	inline static const Pt3 NegativeInfinity();

	inline Pt3();
	inline Pt3(real a, real b, real c);
	inline Pt3(const Pt3 &P);
	inline Pt3& operator +=(const Vec3 &V);
	inline Pt3& operator +=(const Pt3 &Q);
	inline Pt3& operator *=(const real &C);
	inline Pt3& operator -=(const Vec3 &V);
	inline Pt3& operator /=(const real &C);

	/// Cast into a raw array. Suitable for use in glVertex3().
	inline operator real*( );
	/// Cast into a raw array. Suitable for use in glVertex3().
	inline operator const real*( ) const;
};

inline Vec3 operator-(const Pt3 &A, const Pt3 &B);
inline Pt3 operator-(const Pt3 &P, const Vec3 &V);
inline Pt3 operator-(const Pt3 &P);
inline Pt3 operator+(const Pt3 &P, const Vec3 &D);
inline Pt3 operator+(const Pt3 &P, const Pt3 &Q);
inline Pt3 operator*(const real &C, const Pt3 &P);

inline std::ostream& operator<< (std::ostream& os, const Pt3 &P);

// Generic multidimensional point class

template <unsigned int N>
struct Pt{
	real e[N];
	
	Pt(){ for(unsigned int i = 0; i < N; ++i){ e[i] = 0; } }
	Pt(const Pt<N> &P){ for(unsigned int i = 0; i < N; ++i){ e[i] = P.e[i]; } }
	
	Pt<N>& operator +=(const Vec<N> &V){ for(unsigned int i = 0; i < N; ++i){ e[i] += V.e[i]; } return *this; }
	Pt<N>& operator +=(const Pt<N> &Q){ for(unsigned int i = 0; i < N; ++i){ e[i] += Q.e[i]; } return *this; }
	Pt<N>& operator *=(const real &C){ for(unsigned int i = 0; i < N; ++i){ e[i] *= C; } return *this; }
	Pt<N>& operator -=(const Vec<N> &V){ for(unsigned int i = 0; i < N; ++i){ e[i] -= V.e[i]; } return *this; }
	Pt<N>& operator /=(const real &C){ for(unsigned int i = 0; i < N; ++i){ e[i] /= C; } return *this; }
};

template <unsigned int N>
Vec<N> operator-(const Pt<N> &A, const Pt<N> &B){
	Vec<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = A.e[i] - B.e[i];
	}
	return Result;
}

template <unsigned int N>
Pt<N> operator-(const Pt<N> &P, const Vec<N> &V){
	Pt<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = P.e[i] - V.e[i];
	}
	return Result;
}

template <unsigned int N>
Pt<N> operator-(const Pt<N> &P){
	Pt<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = -P.e[i];
	}
	return Result;
}

template <unsigned int N>
Pt<N> operator+(const Pt<N> &P, const Vec<N> &D){
	Pt<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = P.e[i] + D.e[i];
	}
	return Result;
}

template <unsigned int N>
Pt<N> operator+(const Pt<N> &P, const Pt<N> &Q){
	Pt<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = P.e[i] + Q.e[i];
	}
	return Result;
}

template <unsigned int N>
Pt<N> operator*(const real &C, const Pt<N> &P){
	Pt<N> Result;
	for(unsigned int i = 0; i < N; ++i){
		Result.e[i] = C * P.e[i];
	}
	return Result;
}


/********************** Implementation for Pt2 ***********************/

const Pt2 Pt2::Origin(){ static const Pt2 O(0,0); return O; }
const Pt2 Pt2::PositiveInfinity(){ static const Pt2 inf(REAL_MAX,REAL_MAX); return inf; }
const Pt2 Pt2::NegativeInfinity(){ static const Pt2 inf(-REAL_MAX,-REAL_MAX); return inf; }

Pt2::Pt2(){
	e[0] = 0;
	e[1] = 0;
}
Pt2::Pt2(const real &a, const real &b){
	e[0] = a;
	e[1] = b;
}
Pt2::Pt2(const Pt2 &P){
	e[0] = P.e[0];
	e[1] = P.e[1];
}
Pt2& Pt2::operator +=(const Vec2 &V){ e[0] += V.x; e[1] += V.y; return *this; }
Pt2& Pt2::operator +=(const Pt2 &Q){ e[0] += Q.e[0]; e[1] += Q.e[1]; return *this; }
Pt2& Pt2::operator *=(const real &C){ e[0] *= C; e[1] *= C; return *this; }
Pt2& Pt2::operator -=(const Vec2 &V){ e[0] -= V.x; e[1] -= V.y; return *this; }
Pt2& Pt2::operator /=(const real &C){ e[0] /= C; e[1] /= C; return *this; }

Pt2::operator real*( ){ return (real*)e; }
Pt2::operator const real*( ) const{ return (real*)(const real *)e; }

int Pt2::Order(const Pt2 &P0, const Pt2 &P1){
	// test the x-coord first
	if (P0.e[0] > P1.e[0]) return 1;
	if (P0.e[0] < P1.e[0]) return (-1);
	// and test the y-coord second
	if (P0.e[1] > P1.e[1]) return 1;
	if (P0.e[1] < P1.e[1]) return (-1);
	// when you exclude all other possibilities, what remains is...
	return 0;  // they are the same point
}

bool operator<(const Pt2 &A, const Pt2 &B){
	return (A.e[0] < B.e[0] && A.e[1] < B.e[1]);
}

bool operator>(const Pt2 &A, const Pt2 &B){
	return (A.e[0] > B.e[0] && A.e[1] > B.e[1]);
}

Vec2 operator-(const Pt2 &A, const Pt2 &B){
	return Vec2(A.e[0]-B.e[0], A.e[1]-B.e[1]);
}

Pt2 operator-(const Pt2 &P, const Vec2 &V){
	return Pt2(P.e[0]-V.x, P.e[1]-V.y);
}

Pt2 operator-(const Pt2 &P){
	return Pt2(-P.e[0], -P.e[1]);
}

Pt2 operator+(const Pt2 &P, const Vec2 &D){
	return Pt2(P.e[0] + D.x, P.e[1] + D.y);
}

Pt2 operator+(const Pt2 &P, const Pt2 &Q){
	return Pt2(P.e[0] + Q.e[0], P.e[1] + Q.e[1]);
}

Pt2 operator*(const real &C, const Pt2 &P){
	return Pt2(C*P.e[0], C*P.e[1]);
}

/********************** Implementation for Pt3 ***********************/

const Pt3 Pt3::Origin(){ static const Pt3 O(0,0,0); return O; }
const Pt3 Pt3::PositiveInfinity(){ static const Pt3 inf(REAL_MAX,REAL_MAX,REAL_MAX); return inf; }
const Pt3 Pt3::NegativeInfinity(){ static const Pt3 inf(-REAL_MAX,-REAL_MAX,-REAL_MAX); return inf; }

Pt3::Pt3(){
	e[0] = 0;
	e[1] = 0;
	e[2] = 0;
}

Pt3::Pt3(real a, real b, real c){
	e[0] = a;
	e[1] = b;
	e[2] = c;
}

Pt3::Pt3(const Pt3 &P){
	e[0] = P.e[0];
	e[1] = P.e[1];
	e[2] = P.e[2];
}

Pt3& Pt3::operator +=(const Vec3 &V){ e[0] += V.x; e[1] += V.y; e[2] += V.z; return *this; }
Pt3& Pt3::operator +=(const Pt3 &Q){ e[0] += Q.e[0]; e[1] += Q.e[1]; e[2] += Q.e[2]; return *this; }
Pt3& Pt3::operator *=(const real &C){ e[0] *= C; e[1] *= C; e[2] *= C; return *this; }
Pt3& Pt3::operator -=(const Vec3 &V){ e[0] -= V.x; e[1] -= V.y; e[2] -= V.z; return *this; }
Pt3& Pt3::operator /=(const real &C){ e[0] /= C; e[1] /= C; e[2] /= C; return *this; }

Pt3::operator real*( ){ return (real*)e; }
Pt3::operator const real*( ) const{ return (const real*)(const real *)e; }

Vec3 operator-(const Pt3 &A, const Pt3 &B){
	return Vec3(A.e[0]-B.e[0], A.e[1]-B.e[1], A.e[2]-B.e[2]);
}

Pt3 operator-(const Pt3 &P, const Vec3 &V){
	return Pt3(P.e[0]-V.x, P.e[1]-V.y, P.e[2]-V.z);
}

Pt3 operator-(const Pt3 &P){
	return Pt3(-P.e[0], -P.e[1], -P.e[2]);
}

Pt3 operator+(const Pt3 &P, const Vec3 &D){
	return Pt3(P.e[0] + D.x, P.e[1] + D.y, P.e[2] + D.z);
}

Pt3 operator+(const Pt3 &P, const Pt3 &Q){
	return Pt3(P.e[0] + Q.e[0], P.e[1] + Q.e[1], P.e[2] + Q.e[2]);
}

Pt3 operator*(const real &C, const Pt3 &P){
	return Pt3(C*P.e[0], C*P.e[1], C*P.e[2]);
}

std::ostream& operator<< (std::ostream& os, const Pt3 &P){
	os << '(' << P.e[0] << ' ' << P.e[1] << ' ' << P.e[2] << ')';
	return os;
}

#endif