/*******************************************************************************
 * Point<T> (Point3D : Point<double>)
 *
 * 極々ありがちなポイントクラス。外積、内積、角度、長さ、正規化、スカラーとの演
 * 算、行列との演算など一通りサポート。
 *******************************************************************************/

#if !defined ATS_POINT_H
#define ATS_POINT_H

#include <iostream>
#include <limits>
#include <cmath>

// inlineについて
// VC++ではvirtualでもできるところはinlineにしてくれるらしい？

// virtualについて
// 演算子ってvirtualにして良かったっけ？

// 0による除算対策
// numeric_limit<T>::epsilon()
// 例外や無視ではなく最小単位の数とみなし、計算する。
// 整数型を入れたときなどはepsilon()に依存する。

// 注意
// Point型同士だけでなく、Point型 + テンプレート型もあります。
// うっかり足してしまわないように。

// ==演算子は誤差を考慮する。
// 誤差の範囲はGet or Set RangeOfError()関数で定義する。
// ただし、スタティック変数・関数であるため他のインスタンスも
// それに従う。範囲は
// top > value < underで、topとunderはvalueに誤差を加差したもので
// top = value + epsilon(); under = value - epsilon()となる。 

// 誤差の初期値
// numeric_limits<T>::epsilon() * 5。

// 掛け算はデフォルトで外積になります。

namespace ats
{
	// 循環参照解決
	template <class T_FLOAT>
	class Matrix3D;

	inline double clamp(double a, double l, double h)
	{
		return ((a) < (l) ? (l) : (a) > (h) ? (h) : (a));
	}

	template <class T>
	class Point
	{
	//friend Point<T> operator*(const Point<T>& a, const Matrix3D<T>& b);
	private:
		T m_x, m_y, m_z;
		// 誤差の範囲
		static T m_range_of_error;

	public:
		inline explicit Point(const T x=0, const T y=0, const T z=0)
		: m_x(x), m_y(y), m_z(z) {}
		inline explicit Point(const Point<T>& a)
		: m_x(a.m_x), m_y(a.m_y), m_z(a.m_z) {}
		inline ~Point() {}

		inline T X() const { return m_x; }
		inline T Y() const { return m_y; }
		inline T Z() const { return m_z; }
		inline void X(const T a) { m_x = a; }
		inline void Y(const T a) { m_y = a; }
		inline void Z(const T a) { m_z = a; }
		inline void SetXYZ(const T x, const T y, const T z){ m_x = x; m_y = y; m_z = z; }

		inline T Length() const               { return sqrt(m_x*m_x + m_y*m_y + m_z*m_z); }
		inline Point<T>& Normalize()          { return *this /= Length(); }
		inline Point<T> GetNormal() const     { return *this / Length(); }
		inline T Dot(const Point<T>& a) const { return (m_x * a.m_x + m_y * a.m_y + m_z * a.m_z); }
		T Angle(const Point<T>& a) const;

		static T GetRangeOfError();
		static void SetRangeOfError(const T a);

		bool operator==(const Point<T>& a) const;

		/* 意味ないから消そうかな */
		bool operator> (const Point<T>& a) const;
		bool operator>=(const Point<T>& a) const;
		bool operator< (const Point<T>& a) const;
		bool operator<=(const Point<T>& a) const;

		Point<T>& operator= (const Point<T>&);
		Point<T>& operator+=(const Point<T>&);
		Point<T>& operator-=(const Point<T>&);
		Point<T>& operator*=(const Point<T>&);
		Point<T>& operator/=(const Point<T>&);

		inline Point<T> operator+(const Point<T>& a) const { return Point<T>(m_x + a.m_x, m_y + a.m_y, m_z + a.m_z); }
		inline Point<T> operator-(const Point<T>& a) const { return Point<T>(m_x - a.m_x, m_y - a.m_y, m_z - a.m_z); }
		Point<T> operator*(const Point<T>&) const;
		Point<T> operator/(const Point<T>&) const;

		Point<T>& operator+=(const T a);
		Point<T>& operator-=(const T a);
		Point<T>& operator*=(const T a);
		Point<T>& operator/=(const T a);

		Point<T> operator+(const T a) const;
		Point<T> operator-(const T a) const;
		Point<T> operator*(const T a) const;
		Point<T> operator/(const T a) const;

		inline Point<T> operator*(const Matrix3D<T>& a)
		{
			return Point3D(
				m_x * a(0,0) + m_y * a(0,1) + m_z * a(0,2) + a(0,3),
				m_x * a(1,0) + m_y * a(1,1) + m_z * a(1,2) + a(1,3),
				m_x * a(2,0) + m_y * a(2,1) + m_z * a(2,2) + a(2,3)
			);
		}

		inline Point<T>& operator*=(const Matrix3D<T>& a) { return (*this = *this * a); }

		inline Point<T> operator*(const T a[4][4])
		{
			return Point3D(
				m_x * a[0][0] + m_y * a[0][1] + m_z * a[0][2] + a[0][3],
				m_x * a[1][0] + m_y * a[1][1] + m_z * a[1][2] + a[1][3],
				m_x * a[2][0] + m_y * a[2][1] + m_z * a[2][2] + a[2][3]
			);
		}

		inline Point<T>& operator*=(const T a[4][4])
		{
			return (*this = *this * a);
		}
	};

	//#include "ats_Point.cpp"

	// 静的変数 --------------------------------------------------------
	template<class T>
	T Point<T>::m_range_of_error = std::numeric_limits<T>::epsilon() * 5;

	template<class T>
	T Point<T>::Angle(const Point<T>& a) const
	{
		// a･b = |a||b|cosθ
		// (a･b) / |a||b| = cosθ
		// acos((a･b) / |a||b|) = θ
		//
		// a = 60(* pi / 180)
		// cos(a) = b;
		// acos(b) = c;
		// (a == c) is true
		//return acos(clamp(Dot(a) / (Length()*a.Length()), -1.0, 1.0));
		return acos(Dot(a) / (Length()*a.Length()));
		// clampはよくわからん
	}

	// 比較演算子 ------------------------------------------------------
	template<class T>
	bool Point<T>::operator==(const Point<T>& a) const
	{
		Point<T> range_top(*this + m_range_of_error);
		Point<T> range_under(*this - m_range_of_error);
		return (range_top > a && range_under > a);
	}

	template<class T>
	bool Point<T>::operator>(const Point<T>& a) const
	{
		return (m_x > a.m_x && m_y > a.m_y && m_z > a.m_z);
	}

	template<class T>
	bool Point<T>::operator>=(const Point<T>& a) const
	{
		return (m_x >= a.m_x && m_y >= a.m_y && m_z >= a.m_z);
	}

	template<class T>
	bool Point<T>::operator<(const Point<T>& a) const
	{
		return (m_x < a.m_x && m_y < a.m_y && m_z < a.m_z);
	}

	template<class T>
	bool Point<T>::operator<=(const Point<T>& a) const
	{
		return (m_x <= a.m_x && m_y <= a.m_y && m_z <= a.m_z);
	}

	// 誤差 ------------------------------------------------------------
	template<class T>
	T Point<T>::GetRangeOfError()
	{
		return m_range_of_error;
	}

	template<class T>
	void Point<T>::SetRangeOfError(const T a)
	{
		m_range_of_error = a;
	}

	// 代入演算子 --------------------------------------------------------
	template<class T>
	Point<T>& Point<T>::operator=(const Point<T>& a) 
	{
		m_x = a.m_x;
		m_y = a.m_y;
		m_z = a.m_z;
		return *this;
	}

	template<class T>
	Point<T>& Point<T>::operator+=(const Point<T>& a) 
	{
		m_x += a.m_x;
		m_y += a.m_y;
		m_z += a.m_z;
		return *this;
	}

	template<class T>
	Point<T>& Point<T>::operator-=(const Point<T>& a) 
	{
		m_x -= a.m_x;
		m_y -= a.m_y;
		m_z -= a.m_z;
		return *this;
	}

	template<class T>
	Point<T>& Point<T>::operator*=(const Point<T>& a) 
	{
		*this = Point<T>(
			m_y * a.m_z - m_z * a.m_y,
			m_z * a.m_x - m_x * a.m_z,
			m_x * a.m_y - m_y * a.m_x
		);
		return *this;
	}

	template<class T>
	Point<T>& Point<T>::operator/=(const Point<T>& a)
	{
		T x, y, z;
		// 0による除算対策 : 0の場合は最も小さな数とする
		if(0 == a.m_x){ x = std::numeric_limits<T>::epsilon(); }
		else          { x = a.m_x; }
		if(0 == a.m_y){ y = std::numeric_limits<T>::epsilon(); }
		else          { y = a.m_y; }
		if(0 == a.m_z){ z = std::numeric_limits<T>::epsilon(); }
		else          { z = a.m_z; }
		m_x /= x;
		m_y /= y;
		m_z /= z;
		return *this;
	}

	// 代入演算子 対 テンプレート型 ---------------------------------------------------
	template<class T>
	Point<T>& Point<T>::operator+=(const T a) 
	{
		m_x += a;
		m_y += a;
		m_z += a;
		return *this;
	}
	template<class T>
	Point<T>& Point<T>::operator-=(const T a) 
	{
		m_x -= a;
		m_y -= a;
		m_z -= a;
		return *this;
	}
	template<class T>
	Point<T>& Point<T>::operator*=(const T a) 
	{
		m_x *= a;
		m_y *= a;
		m_z *= a;
		return *this;
	}
	template<class T>
	Point<T>& Point<T>::operator/=(const T a) 
	{
		T tmp;
		if(0 == a){ tmp = std::numeric_limits<T>::epsilon(); }
		else { tmp = a; }
		m_x /= tmp;
		m_y /= tmp;
		m_z /= tmp;
		return *this;
	}

	// 四則演算子 ---------------------------------------------------------

	template<class T>
	Point<T> Point<T>::operator*(const Point<T>& a) const 
	{
		return Point<T>(m_y * a.m_z - m_z * a.m_y,
		                m_z * a.m_x - m_x * a.m_z,
		                m_x * a.m_y - m_y * a.m_x);
	}

	template<class T>
	Point<T> Point<T>::operator/(const Point<T>& a) const
	{
		T x, y, z;
		// 0による除算対策 : 0の場合は最も小さな数とする
		if(0 == a.m_x){ x = std::numeric_limits<T>::epsilon(); }
		else          { x = a.m_x; }
		if(0 == a.m_y){ y = std::numeric_limits<T>::epsilon(); }
		else          { y = a.m_y; }
		if(0 == a.m_z){ z = std::numeric_limits<T>::epsilon(); }
		else          { z = a.m_z; }
		return Point<T>(m_x/x, m_y/y, m_z/z);
	}

	// 四則演算子 対 テンプレート型 ----------------------------------------
	template<class T>
	Point<T> Point<T>::operator+(const T a) const 
	{
		return Point<T>(m_x + a, m_y + a, m_z + a);
	}

	template<class T>
	Point<T> Point<T>::operator-(const T a) const 
	{
		return Point<T>(m_x - a, m_y - a, m_z - a);
	}

	template<class T>
	Point<T> Point<T>::operator*(const T a) const 
	{
		return Point<T>(m_x * a, m_y * a, m_z * a);
	}

	template<class T>
	Point<T> Point<T>::operator/(const T a) const
	{
		T tmp;
		// 0による除算対策 : 0の場合は最も小さな数とする
		if(0 == a){ tmp = std::numeric_limits<T>::epsilon(); }
		else      { tmp = a; }
		return Point<T>(m_x/tmp, m_y/tmp, m_z/tmp);
	}

	template <class T>
	std::ostream& operator<<(std::ostream& out, const Point<T>& a)
	{
		out << a.X() << " : " << a.Y() << " : " << a.Z() << endl;
		return out;
	}

	typedef Point<double> Point3D;

	template <class T>
	Point<T> operator-(const Point<T>& a)
	{
		return a * -1;
	}
}

#endif // #if !defined ATS_POINT_H