#ifndef  __JBXL_CPP_VECTOR_H_
#define  __JBXL_CPP_VECTOR_H_


/**
@brief    ベクトルライブラリ for C++
@file     Vector.h
@author   Fumi.Iseki (C)
*/



/*
	template <typename T=double> class DllExport Vector		// ３次元ベクトルの定義
	template <typename T=double> class DllExport Vector2D	// ２次元ベクトルの定義, パラメトリック曲面用
	template <typename T=int> class DllExport  PCoordinate	// 多角形の座標の値を格納するクラス
	template <typename T=int> class DllExport RBound		// 境界構造体

	template <typename T> inline bool disJunctBounds(RBound<T> b1, RBound<T> b2);
*/



#include "tools++.h"
#include <math.h>
#include "Tolerance.h"



//
namespace  jbxl {



#define  VECTOR    Vector
#define  VECTOR2D  Vector2D


#define  BOUNDARY_BLANK 5		///< 境界 RBound<T> の余白





/**
template <typename T=double> class Vector

	３次元ベクトルの定義

*/
template <typename T=double> class DllExport Vector
{
public:
	T x;
	T y;
	T z;

	double n;		///< ノルム
	double c;		///< 信頼度

public:
	Vector(T X=0, T Y=0, T Z=0, double N=0.0, double C=1.0) { set(X, Y, Z, N, C);}
	~Vector(void) {}

	T      norm2(void) { return (x*x + y*y + z*z);}
	double norm(void)  { n = sqrt((double)(x*x + y*y + z*z)); return n;}
	void   normalize(void);

	void   set(T X, T Y=0, T Z=0, double N=0.0, double C=1.0);
	void   init(double C=1.0) { x = y = z = (T)0; n = 0.0; c = C;}
};




template <typename T> void Vector<T>::normalize(void)
{
	double nrm = sqrt((double)x*x+y*y+z*z);

	if (nrm>=Zero_Eps) {
		x  = (T)(x/nrm);
		y  = (T)(y/nrm);
		z  = (T)(z/nrm);
		n  = 1.0;
	}
	else {
		init();
	}
}



/**
template <typename T> void Vector<T>::set(T X, T Y, T Z, double N) 

3次元ベクトルに値をセット．

@attention
ノルムの計算無し（高速化のため）
*/
template <typename T> void Vector<T>::set(T X, T Y, T Z, double N, double C) 
{
	x  = X;
	y  = Y;
	z  = Z;
	n  = N;
	c  = C;

	//if (n<=0.0) {
	//	n = sqrt((double)(x*x + y*y + z*z));
	//}
}





/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// オペレータ

template <typename T> inline Vector<T> operator - (const Vector<T> a)
{ return Vector<T>(-a.x, -a.y, -a.z, a.n, a.c); }


template <typename T> inline Vector<T> operator + (const Vector<T> a, const Vector<T> b)
{ return Vector<T>(a.x+b.x, a.y+b.y, a.z+b.z, 0.0, Min(a.c, b.c)); }


template <typename T> inline Vector<T> operator - (const Vector<T> a, const Vector<T> b)
{ return Vector<T>(a.x-b.x, a.y-b.y, a.z-b.z, 0.0, Min(a.c, b.c)); }


template <typename T, typename R> inline Vector<T> operator * (const R d, const Vector<T> a)
{ return Vector<T>((T)d*a.x, (T)d*a.y, (T)d*a.z, d*a.n, a.c); }


template <typename T, typename R> inline Vector<T> operator * (const Vector<T> a, const R d)
{ return Vector<T>(a.x*(T)d, a.y*(T)d, a.z*(T)d, a.n*d, a.c); }


template <typename T, typename R> inline Vector<T> operator / (const Vector<T> a, const R d)
{ return Vector<T>(a.x/(T)d, a.y/(T)d, a.z/(T)d, a.n/d, a.c); }


/// Cross product  外積
template <typename T> inline Vector<T> operator ^ (const Vector<T> a, const Vector<T> b)
{ return Vector<T>(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x, 0.0, Min(a.c, b.c)); }


/// Dot product    内積
template <typename T> inline T operator * (const Vector<T> a, const Vector<T> b)
{ return (a.x*b.x + a.y*b.y + a.z*b.z); }


template <typename T> inline bool operator == (const Vector<T> v1, const Vector<T> v2)
{ return (v1.x==v2.x && v1.y==v2.y && v1.z==v2.z); }


//
template <typename T> inline Vector<T> MidPoint(const Vector<T> a, const Vector<T> b)
{ return 0.5*(a+b); }


//
/// 点a と b の距離　（a,b は位置ベクトル）
template <typename T> inline double VectorDist(const Vector<T> a, const Vector<T> b)
{
	//return (b-a).n;
	return (b-a).norm();
}



template <typename T> inline bool operator < (const Vector<T> v1, const Vector<T> v2)
{
	if (v1.x != v2.x) return v1.x < v2.x;
	if (v1.y != v2.y) return v1.y < v2.y;
	if (v1.z != v2.z) return v1.z < v2.z;
	return false;
}





/////////////////////////////////////////////////////////////////////////////////////////////////////////////

/**
template <typename T=double> class Vector2D

２次元ベクトルの定義@n
パラメトリック曲面用
*/
template <typename T=double> class DllExport Vector2D
{
public:
	T  u;
	T  v;

public:
	Vector2D(const T U=0, const T V=0) : u(U), v(V) {}
};





/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 3D多角形座標

/**
template <typename T=int> class DllExport  PCoordinate

多角形の座標の値を格納するクラス

*/
template <typename T=int> class DllExport  PCoordinate
{
public:
	int  dim;
	Vector<T>** point;

public:
	
	PCoordinate(void)  { dim = 0; point = NULL;}
	PCoordinate(int n) { init(n);}

	~PCoordinate(void) {}


	//
	void  init(int n) {
		dim   = 0;
		point = NULL;

		if (n>0) {
			point = (Vector<T>**)malloc(sizeof(Vector<T>*)*n);
			if (point!=NULL) {
				dim = n;
				for (int i=0; i<dim; i++) point[i] = new Vector<T>();
			}
		}
	}


	void  set(int m, T x, T y=0, T z=0, double n=0.0) {
		if (m>=0 && m<dim && point!=NULL) point[m]->set(x, y, z, n); 
	}		
	

	void  clear(void) {
		if (point!=NULL) {
			for (int i=0; i<dim; i++) point[i]->set((T)0);
		}
	}


	void  free(void) {
		if (point!=NULL) {
			for (int i=0; i<dim; i++) delete(point[i]);
			::free(point);
			point = NULL;
		}
		dim = 0;
	}

};






/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 境界構造体

/**
template <typename T=int> class RBound 

3Dデータの境界情報を保存する型． 
*/
template <typename T=int> class DllExport RBound 
{
public:
    T   xmin;		///< x軸境界の最小値．
    T   xmax;		///< x軸境界の最大値．
    T   ymin;		///< y軸境界の最小値．
    T   ymax;		///< y軸境界の最大値．
    T   zmin;		///< z軸境界の最小値．
    T   zmax;		///< z軸境界の最大値．
	T   tmin;		///< 汎用．
	T   tmax;		///< 汎用．

public:

	RBound(T XMin=(T)0, T XMax=(T)0, T YMin=(T)0, T YMax=(T)0, T ZMin=(T)0, T ZMax=(T)0, T TMin=(T)0, T TMax=(T)0) {
		set(XMin, XMax, YMin, YMax, ZMin, ZMax, TMin, TMax);
	}

	~RBound() {}


	//
	void  set(T XMin=(T)0, T XMax=(T)0, T YMin=(T)0, T YMax=(T)0, T ZMin=(T)0, T ZMax=(T)0, T TMin=(T)0, T TMax=(T)0){
		xmin = XMin;
		ymin = YMin;
		zmin = ZMin;
		xmax = XMax;
		ymax = YMax;
		zmax = ZMax;
		tmin = TMin;
		tmax = TMax;
	}		

	void init() { xmin = xmax = ymin = ymax = zmin = zmax = tmin = tmax = (T)0;}




	///////////////////////////////////////////
	// 拡大

	/// 境界構造体 bound を広げる
	void  enlarge(T f) {
		xmin -= f;
		ymin -= f;
		zmin -= f;
		xmax += f;
		ymax += f;
		zmax += f;
	}



	///////////////////////////////////////////
	// 融合

	/// 境界構造体 bound と融合させる
	void  fusion(RBound<T> bound) {
		xmin = Min(xmin, bound.xmin);
		ymin = Min(ymin, bound.ymin);
		zmin = Min(zmin, bound.zmin);
		xmax = Max(xmax, bound.xmax);
		ymax = Max(ymax, bound.ymax);
		zmax = Max(zmax, bound.zmax);
	}


	/// ベクトル vect と融合させる  
	void  fusion(Vector<T> vect){
		xmin = Min(xmin, vect.x);
		ymin = Min(ymin, vect.y);
		zmin = Min(zmin, vect.z);
		xmax = Max(xmax, vect.x);
		ymax = Max(ymax, vect.y);
		zmax = Max(zmax, vect.z);
	}


	/// ポイントと融合させる  
	void  fusion(T x, T y, T z){
		xmin = Min(xmin, x);
		ymin = Min(ymin, y);
		zmin = Min(zmin, z);
		xmax = Max(xmax, x);
		ymax = Max(ymax, y);
		zmax = Max(zmax, z);
	}




	///////////////////////////////////////////
	// 共通領域

	/// 境界構造体 bound との共通領域  共通領域がない場合は，min>max になる
	void  commonarea(RBound<T> bound) {
		xmin = Max(xmin, bound.xmin);
		ymin = Max(ymin, bound.ymin);
		zmin = Max(zmin, bound.zmin);
		xmax = Min(xmax, bound.xmax);
		ymax = Min(ymax, bound.ymax);
		zmax = Min(zmax, bound.zmax);
	}



	///////////////////////////////////////////
	/// 切り出した場合の境界
	void  cutdown(RBound<T> bound) {
		xmin += bound.xmin;
		ymin += bound.ymin;
		zmin += bound.zmin;
		xmax += bound.xmin;
		ymax += bound.ymin;
		zmax += bound.zmin;
	}


	//   
	void  cutdown(Vector<T> vect){
		xmin += vect.x;
		ymin += vect.y;
		zmin += vect.z;
		xmax += vect.x;
		ymax += vect.y;
		zmax += vect.z;
	}

	
	//  
	void  cutdown(T x, T y, T z){
		xmin += x;
		ymin += y;
		zmin += z;
		xmax += x;
		ymax += y;
		zmax += z;
	}



	/// ベクトルが位置ベクトルの場合，その点は境界外か？ 境界外：true，境界内：false
	bool  outofBounds(Vector<T> vect) {
		return	vect.x < xmin || vect.x > xmax ||
				vect.y < ymin || vect.y > ymax ||
				vect.z < zmin || vect.z > zmax;
	}


	/// その点は境界外か？ 境界外：true，境界内：false
	bool  outofBounds(T x, T y, T z) {
		return	x < xmin || x > xmax ||
				y < ymin || y > ymax ||
				z < zmin || z > zmax;
	}

};




/**
template <typename T> inline bool disJunctBounds(RBound<T> b1, RBound<T> b2)
	
境界領域b1とb2が重なっているかどうかをチェック
   
@retval true   重なっていない
@retval false  重なっている 
*/
template <typename T> inline bool disJunctBounds(RBound<T> b1, RBound<T> b2)
{
	return	(b1.xmin >= b2.xmax) || (b2.xmin >= b1.xmax) || 
			(b1.ymin >= b2.ymax) || (b2.ymin >= b1.ymax) || 
			(b1.zmin >= b2.zmax) || (b2.zmin >= b1.zmax);
}




//
/// Normal Vector with Newell Mothod
template <typename T> inline Vector<T> NewellMethod(Vector<T> v1, Vector<T> v2, Vector<T> v3)
{
	Vector<T> vect;

	vect.x = (v1.y-v2.y)*(v1.z+v2.z) + (v2.y-v3.y)*(v2.z+v3.z) + (v3.y-v1.y)*(v3.z+v1.z);
	vect.y = (v1.z-v2.z)*(v1.x+v2.x) + (v2.z-v3.z)*(v2.x+v3.x) + (v3.z-v1.z)*(v3.x+v1.x);
	vect.z = (v1.x-v2.x)*(v1.y+v2.y) + (v2.x-v3.x)*(v2.y+v3.y) + (v3.x-v1.x)*(v3.y+v1.y);
	vect.n = 0.0;
	vect.c = Min(v1.c, v2.c);
	vect.c = Min(v3.c, vect.c);
	return vect;
}



template <typename T> inline Vector<T> NewellMethod3(Vector<T> v1, Vector<T> v2, Vector<T> v3)
{
	return NewellMethod(v1, v2, v3);
}



template <typename T> inline Vector<T> NewellMethod4(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4)
{
	Vector<T> vect;

	vect.x = (v1.y-v2.y)*(v1.z+v2.z) + (v2.y-v3.y)*(v2.z+v3.z) + (v3.y-v4.y)*(v3.z+v4.z) + (v4.y-v1.y)*(v4.z+v1.z);
	vect.y = (v1.z-v2.z)*(v1.x+v2.x) + (v2.z-v3.z)*(v2.x+v3.x) + (v3.z-v4.z)*(v3.x+v4.x) + (v4.z-v1.z)*(v4.x+v1.x);
	vect.z = (v1.x-v2.x)*(v1.y+v2.y) + (v2.x-v3.x)*(v2.y+v3.y) + (v3.x-v4.x)*(v3.y+v4.y) + (v4.x-v1.x)*(v4.y+v1.y);
	vect.n = 0.0;
	vect.c = Min(v1.c, v2.c);
	vect.c = Min(v3.c, vect.c);
	vect.c = Min(v4.c, vect.c);
	return vect;
}




}		// namespace


#endif