
#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);

template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b)
template <typename T> inline Vector<T> NewellMethod(Vector<T> v1, Vector<T> v2, Vector<T> v3)
template <typename T> inline Vector<T> BSplineInterp4(Vector<T> p0, Vector<T> p1, double t)
*/



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



//
namespace  jbxl {



#define  VECTOR		Vector
#define  UVMAP	  	UVMap


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





/////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// Vector
//

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

	３次元ベクトルの定義

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

	double n;		///< ノルム
	double c;		///< 信頼度
	int    d;		///< 汎用

public:
	Vector(T X=0, T Y=0, T Z=0, double N=0.0, double C=1.0, int D=0) { set(X, Y, Z, N, C, D);}
	virtual ~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;}
	Vector<T> normalize(void);

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

	// 別名
	T	   element1(void) { return x;}
	T	   element2(void) { return y;}
	T	   element3(void) { return z;}
	T	   element(int i) { if(i==1) return x; else if(i==2) return y; else if(i==3) return z; else return (T)0;} 
};



template <typename T> Vector<T> 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();
	}

	return *this;
}



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

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

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

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



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

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


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), a.d); }

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

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


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), a.d); }

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

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


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, a.d); }

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, a.d); }


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, a.d); }

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


/// 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 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;
}



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

/**
v1, v2 が同じ点かチェックする．信頼度は考慮しない．

*/
template <typename T> inline bool same_vector(Vector<T> v1, Vector<T> v2)
{
	T dx = v1.x - v2.x;
	T dy = v1.y - v2.y;
	T dz = v1.z - v2.z;
	T d2 = dx*dx + dy*dy + dz*dz;
	double t2 = Vector_Tolerance*Vector_Tolerance;

	if ((double)d2>t2) return false;
	return true;
}



template <typename T> inline Vector<T>* dupVector(Vector<T>* a, int n)
{
	Vector<T>* v = (Vector<T>*)malloc(sizeof(Vector<T>)*n);
	if (v==NULL) return NULL;

	for (int i=0; i<n; i++) v[i] = a[i];

	return v;
}


// angle between a and b
template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b)
{
	a.normalize();
	b.normalize();
	if (a.n<Zero_Eps || b.n<Zero_Eps) return 0.0;

	double cs = a*b;
	if      (cs>=1.0)  return 0.0;
	else if (cs<=-1.0) return PI;
	
	return acos(a*b);
}


template <typename T> inline double VectorAngle(Vector<T> a, Vector<T> b,  Vector<T> c)
{
	return VectorAngle(b-a, c-b);
}


//
/// Normal Vector of 3 Vectors 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;
}



/**
4階（3次）Bスプライン曲線補間．
2次の導関数が両端で 0の条件有り．

P(t) = 1/6*(1-t)^3*Q0 + (1/2*t^3 - t^2 + 2/3)*Q1 + (-1/2*t^3 + 1/2*t^2 + 1/2*t + 1/6)*Q2 + 1/6*t^3*Q3
	Q0 = 2*P(0) - P(1);
	Q1 = P(0)
	Q2 = P(1)
	Q3 = 2*P(1) - P(0)
*/
template <typename T> inline Vector<T> BSplineInterp4(Vector<T> p0, Vector<T> p1, double t)
{
	Vector<T> q0 = 2*p0 - p1;
//	Vector<T> q1 = p0;
//	Vector<T> q2 = p1;
	Vector<T> q3 = 2*p1 - p0;

	double t2 = t*t;
	double t3 = t2*t;
	
	double c0 = (t2 - t)*0.5 + (1.0 - t3)*0.1666666666666667;	// 1/6*(1-t)^3
	double c1 = t3*0.5 - t2 + 0.6666666666666667;				// 1/2*t3 - t2 + 2/3
	double c2 = (t + t2 - t3)*0.5 + 0.1666666666666667;			// -1/2*t3 + 1/2*t2 + 1/2*t + 1/6
	double c3 = t3*0.1666666666666667;							// 1/6*t3

	Vector<T> vect = c0*q0 + c1*p0 + c2*p1 + c3*q3;

	return vect;
}





/////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// 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);}

	virtual ~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);
	}

	virtual ~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;
	}

	void  multiple(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);
}





/////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
// UVMap
//

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

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

public:
	UVMap(T U=0, T V=0, int d=0) { set(U, V, d);}
	virtual ~UVMap(void) {}

	void   init(void) { u = v = (T)0; d = 0;}
	void   set(T U, T V=0, int D=0) { u = U, v = V; d = D;}

	UVMap flip(void) { u = 1.0 - u; v = 1.0 - v; return *this;}
};




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

template <typename T> inline UVMap<T> operator - (const UVMap<T> a)
{ return UVMap<T>(-a.u, -a.v, a.d); }


template <typename T> inline UVMap<T> operator + (const UVMap<T> a, const UVMap<T> b)
{ return UVMap<T>(a.u+b.u, a.v+b.v, a.d); }

template <typename T, typename R> inline UVMap<T> operator + (const UVMap<T> a, R c)
{ return UVMap<T>(a.u+(T)c, a.v+(T)c, a.d); }

template <typename T, typename R> inline UVMap<T> operator + (const R c, UVMap<T> a)
{ return UVMap<T>((T)c+a.u, (T)c+a.v, a.d); }


template <typename T> inline UVMap<T> operator - (const UVMap<T> a, const UVMap<T> b)
{ return UVMap<T>(a.u-b.u, a.v-b.v, a.d); }

template <typename T, typename R> inline UVMap<T> operator - (const UVMap<T> a, R c)
{ return UVMap<T>(a.u-(T)c, a.v-(T)c, a.d); }

template <typename T, typename R> inline UVMap<T> operator - (const R c, UVMap<T>a)
{ return UVMap<T>((T)c-a.u, (T)c-a.v, a.d); }


template <typename T, typename R> inline UVMap<T> operator * (const R d, const UVMap<T> a)
{ return UVMap<T>((T)d*a.u, (T)d*a.v, a.d); }

template <typename T, typename R> inline UVMap<T> operator * (const UVMap<T> a, const R d)
{ return UVMap<T>(a.u*(T)d, a.v*(T)d, a.d); }


template <typename T, typename R> inline UVMap<T> operator / (const UVMap<T> a, const R d)
{ return UVMap<T>(a.u/(T)d, a.v/(T)d, a.d); }

template <typename T, typename R> inline UVMap<T> operator / (const R d, const UVMap<T> a)
{ return UVMap<T>((T)d/a.u, (T)d/a.v, a.d); }


template <typename T> inline bool operator == (const UVMap<T> a, const UVMap<T> b)
{ return (a.u==b.u && a.v==b.v); }

template <typename T> inline bool operator != (const UVMap<T> a, const UVMap<T> b)
{ return (a.u!=b.u || a.v!=b.v); }



}		// namespace

#endif 


