// © Copyright 2008 Владимир Фролов
#define MGML_MATH_INTERSECTIONS_GUARDIAN

#ifndef MGML_MATH_GUARDIAN
	#include "MGML_MATH.h"
#endif

#ifndef __CUDACC__
	#ifndef MGML_MATH_SIMD_GUARDIAN
		#include "MGML_MATH_SIMD.h"
	#endif
#endif


#ifndef EXT_TRIBOXOVERLAP_GUARDIAN
	#include "EXT_TriBoxOverlap.h"
#endif


namespace MGML_MATH
{


// интерсекторы, то есть классы, отвечающие за 
// результат пересечений.

struct ray_box_face_intersector
{
	int				face_number[2];
	VECTOR<2,float>	flat_coord[2];	
};

///////////////////////////////////////////////////////////////////////////////////
////
template<class A,class B>
struct Intersect
{ 
 static  bool exec(const A& obj1,const B& obj2)
 { ERROR_INTERSECTION_UNDEFINED(obj1,obj2); return false; }

};


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

template<>
struct Intersect<RAY<3,float>, AABB<4,float> >
{
	typedef AABB<4,float>   Box;
	typedef RAY<3,float>    Ray;
	typedef VECTOR<2,float> vec2f;
	typedef VECTOR<4,float> vec4f;

	static universal_call bool exec(const Ray& r,const Box& box,float* p_tmin,float* p_tmax)
	{
		register int sign1, sign2, sign3;
		// fucking errors of nvcc with constructors
		// i know - thats bad
		const VECTOR<4,float>* V = box.CVPointer();
		// \\ now we can use V

		sign1 = (r.dir[0] < 0);
		sign2 = (r.dir[1] < 0);

		register float tmin = ( V[sign1][0] - r.pos[0]) / r.dir[0];
		register float tmax = ( V[1 - sign1][0] - r.pos[0]) / r.dir[0];
		register float tymin = ( V[sign2][1] - r.pos[1]) / r.dir[1];
		register float tymax = ( V[1 - sign2][1] - r.pos[1]) / r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		sign3 = ( r.dir[2] < 0);

		register float tzmin = ( V[sign3][2] - r.pos[2]) / r.dir[2];
		register float tzmax = ( V[1 - sign3][2] - r.pos[2]) / r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		*p_tmin = tmin;
		*p_tmax = tmax;
		return (tmax > 0) && ( tmin < tmax);
	}


  static universal_call bool exec(const Ray& r,const Box& box)
	{
		register int sign1, sign2, sign3;
		// fucking errors of nvcc with constructors
		// i know - thats bad
		const VECTOR<4,float>* V = box.CVPointer();
		// \\ now we can use V

		sign1 = (r.dir[0] < 0);
		sign2 = (r.dir[1] < 0);

		register float tmin = ( V[sign1][0] - r.pos[0]) / r.dir[0];
		register float tmax = ( V[1 - sign1][0] - r.pos[0]) / r.dir[0];
		register float tymin = ( V[sign2][1] - r.pos[1]) / r.dir[1];
		register float tymax = ( V[1 - sign2][1] - r.pos[1]) / r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		sign3 = ( r.dir[2] < 0);

		register float tzmin = ( V[sign3][2] - r.pos[2]) / r.dir[2];
		register float tzmax = ( V[1 - sign3][2] - r.pos[2]) / r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		return (tmax > 0) && ( tmin < tmax);
	}
};


///////////////////////////////////////////////////////////////////////////////////
////
template<int n,class T>
struct Intersect< SEGMENT<n,T>, SPHERE<n,T> >
{
	typedef VECTOR<n,T>  vec;

	inline static universal_call bool exec(const SEGMENT<n,T> seg, const SPHERE<n,T>& sph)
	{
		if( (seg.begin() - sph.pos).lenSquare() < sph.rSquare() || 
			(seg.end() - sph.pos).lenSquare() < sph.rSquare() )
			return true;

		const vec& rpos = seg.begin();
		vec rdir = (seg.end() - seg.begin());
		float seg_len = rdir.len();
		rdir /= seg_len;
        const vec& spos = sph.pos;

		register float b,c,d;
		register vec k = rpos - spos;

		//a == 1; // because rdir must be normalized
		b = dot(k,rdir);
		c = dot(k,k) - sph.rSquare();
		d = b*b - c;
		
		if( d < 0 || b > 0 || c < 0) 
			return false;

		float sqrtfd = sqrtf(d);
		// нам нужен минимальный t, a == 1
		float t1 =  -sqrtfd - b;
		float t2 =  -sqrtfd + b;

		if(t1 < 0 || t1 > seg_len)
			return false;

		return true;
	}

};








///////////////////////////////////////////////////////////////////////////////////
////
template<int n,class T>
struct Intersect< RAY<n,T>, SPHERE<n,T> >
{
    typedef SPHERE<n,T> _sphere;
	typedef RAY<n,T>    _ray;
	typedef VECTOR<n,T> vec;

// ray-sphere
	inline static universal_call bool exec(const _ray& ray,const _sphere& sph)
	{
		const vec& rpos = ray.pos;
		const vec& rdir = ray.dir;
        const vec& spos = sph.pos;

		register float b,c,d;
		register vec k = rpos - spos;

		//a == 1; // because rdir must be normalized
		b = dot(k,rdir);
		c = dot(k,k) - sph.rSquare();
		d = b*b - c;

		return (d < 0 || b > 0 || c < 0); 
	}

	inline static universal_call bool exec(const _ray& ray,
										   const vec& spos,
										   float sphere_rSquare,
										   vec* pPos,
										   vec* pNorm,
										   float* t_near)
	{

		float b,c,d;
		vec k = ray.pos - spos;

		//a == 1; // because rdir must be normalized
		b = dot(k,ray.dir);
		c = dot(k,k) - sphere_rSquare;
		d = b*b - c;

		if( d < 0 || b > 0 || c < 0) 
			return false;
	
		float sqrtfd = sqrtf(d);
		// нам нужен минимальный t, a == 1
		float t =  -(sqrtfd + b);

		*t_near = t;
		*pPos  = ray.pos + ray.dir*t;  

		*pNorm = (*pPos - spos);
		pNorm->Normalize();

		return true;
	}

	inline static universal_call bool exec(const _ray& ray,
										   const vec& spos,
										   float sphere_rSquare,
										   float* t_near,
									       float* t_far)
	{
		const vec& rpos = ray.pos;
		const vec& rdir = ray.dir;

		float b,c,d;
		vec k = rpos - spos;

		//a == 1; // because rdir must be normalized
		b = dot(k,rdir);
		c = dot(k,k) - sphere_rSquare;
		d = b*b - c;

		if( d < 0 ) 
			return false;
	
		float sqrtfd = sqrtf(d);
		// нам нужен минимальный t, a == 1
    float tfar = -sqrtfd + b;
    float tnear = -sqrtfd - b;

    *t_far = -sqrtfd - b;
		*t_near= sqrtfd + b;

		return true;
	}

  inline static universal_call bool exec(const _ray& ray,
										   const vec& spos,
										   float sphere_rSquare,
										   float* t_near)
	{
		const vec& rpos = ray.pos;
		const vec& rdir = ray.dir;

		float b,c,d;
		vec k = rpos - spos;

		//a == 1; // because rdir must be normalized
		b = dot(k,rdir);
		c = dot(k,k) - sphere_rSquare;
		d = b*b - c;

		if( d < 0 ) 
			return false;
	
		float sqrtfd = sqrtf(d);
		// нам нужен минимальный t, a == 1
		//*t_near= -sqrtfd - b;

    //*t_far = -sqrtfd - b;
		*t_near= sqrtfd + b;

		return true;
	}

};

///////////////////////////////////////////////////////////////////////////////////
////
template<int n, class T>
struct Intersect<RAY<n,T>, AABB<n,T> >
{
	typedef AABB<n,T>   Box;
	typedef RAY<n,T>    Ray;
	typedef VECTOR<2,T> vec2f;
	typedef VECTOR<n,T> vec;

	static inline universal_call bool exec(const Ray& r, const Box& box, float* t_min, float* t_max)
	{
		register int sign1, sign2;
		// fucking errors of nvcc with constructors
		// i know - thats bad
		const VECTOR<n,T>* V = box.CVPointer();
		// \\ now we can use V

		sign1 = (r.dir[0] < 0);
		sign2 = (r.dir[1] < 0);

		register float tmin = ( V[ sign1][0] - r.pos[0]) / r.dir[0];
		register float tmax = ( V[ 1 - sign1][0] - r.pos[0]) / r.dir[0];
		register float tymin = ( V[ sign2][1] - r.pos[1]) / r.dir[1];
		register float tymax = ( V[ 1 - sign2][1] - r.pos[1]) / r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		register int sign3 = ( r.dir[2] < 0);

		register float tzmin = ( V[sign3][2] - r.pos[2]) / r.dir[2];
		register float tzmax = ( V[1 - sign3][2] - r.pos[2]) / r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		*t_min = tmin;
		*t_max = tmax;
		return (tmax > 0) && ( tmin <= tmax);
	}

  static universal_call bool execRayDirInverse(const Ray& r, const Box& box, T& tmin, T& tmax)
	{
	  	// i know - thats bad
		const VECTOR<n,T>* V = box.CVPointer();
		// \\ now we can use V

		tmin = ( V[(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		tmax = ( V[1-(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		float tymin = ( V[(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];
		float tymax = ( V[1-(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		int sign3 = (r.dir[2]<0);

		float tzmin = (V[sign3][2] - r.pos[2]) * r.dir[2];
		float tzmax = (V[1-sign3][2] - r.pos[2]) * r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		return (tmax > 0) && ( tmin <= tmax);
	}

  static inline universal_call bool exec(const Ray& r, const Box& box)
	{
		// i know - thats bad
		const VECTOR<n,T>* V = box.CVPointer();
		// \\ now we can use V

		register int sign1 = (r.dir[0] < 0);
		register int sign2 = (r.dir[1] < 0);

		register float tmin = ( V[sign1][0] - r.pos[0]) / r.dir[0];
		register float tmax = ( V[1-sign1][0] - r.pos[0]) / r.dir[0];
		register float tymin = ( V[sign2][1] - r.pos[1]) / r.dir[1];
		register float tymax = ( V[1-sign2][1] - r.pos[1]) / r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		register int sign3 = ( r.dir[2] < 0);

		register float tzmin = ( V[sign3][2] - r.pos[2]) / r.dir[2];
		register float tzmax = ( V[1-sign3][2] - r.pos[2]) / r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		return (tmax > 0) && ( tmin <= tmax);
	}

  static inline universal_call bool execRayDirInverse(const Ray& r, const Box& box)
	{
	 	// i know - thats bad
		const VECTOR<n,T>* V = box.CVPointer();
		// \\ now we can use V

		float tmin = ( V[(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		float tmax = ( V[1-(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		float tymin = ( V[(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];
		float tymax = ( V[1-(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		int sign3 = (r.dir[2]<0);

		float tzmin = (V[sign3][2] - r.pos[2]) * r.dir[2];
		float tzmax = (V[1-sign3][2] - r.pos[2]) * r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		return (tmax > 0) && ( tmin <= tmax);
	}

  static inline universal_call bool execRayDirInverse(const Ray& r, const Box& box, T& tmin)
	{
	 	// i know - thats bad
		const VECTOR<n,T>* V = box.CVPointer();
		// \\ now we can use V

		tmin = ( V[(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		float tmax = ( V[1-(r.dir[0] < 0)][0] - r.pos[0]) * r.dir[0];
		float tymin = ( V[(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];
		float tymax = ( V[1-(r.dir[1] < 0)][1] - r.pos[1]) * r.dir[1];

		if ( (tmin > tymax) || (tymin > tmax) )
			return false;

		if (tymin > tmin)
			tmin = tymin;

		if (tymax < tmax)
			tmax = tymax;

		int sign3 = (r.dir[2]<0);

		float tzmin = (V[sign3][2] - r.pos[2]) * r.dir[2];
		float tzmax = (V[1-sign3][2] - r.pos[2]) * r.dir[2];

		if ( (tmin > tzmax) || (tzmin > tmax) )
			return false;

		if (tzmin > tmin)
			tmin = tzmin;
		if (tzmax < tmax)
			tmax = tzmax;

		return (tmax > 0) && ( tmin <= tmax);
	}


};

///////////////////////////////////////////////////////////////////////////////////
//// Plane-AABB
template<int n, class T>
struct Intersect< PLANE<n,T>, AABB<n,T> >
{
	static universal_call bool exec(const PLANE<n,T>& plane, const AABB<n,T>& box)
	{
		return false;
	}
};



#ifndef __CUDACC__
///////////////////////////////////////////////////////////////////////////////////
//// ray-AABB
template<>
struct Intersect<RAY<4,float>,AABB<4,float> >
{
	typedef AABB<4,float>   _box;
	typedef RAY<4,float>    _ray;

	static const float flt_plus_inf;
    static __declspec(align(16)) const float ps_cst_plus_inf[4];
	static __declspec(align(16)) const float ps_cst_minus_inf[4];
	static __declspec(align(16)) const float ps_cst_one[4];


	static bool exec(const _ray& ray,const _box& box,float* t_near,float* t_far)
	{
		// you may already have those values hanging around somewhere
		const __m128
			plus_inf	= _mm_load_ps(ps_cst_plus_inf),
			minus_inf	= _mm_load_ps(ps_cst_minus_inf),
			one         = _mm_load_ps(ps_cst_one);

	   VECTOR<4,float> box_min = box.vmin;
	   VECTOR<4,float> box_max = box.vmax;
	   VECTOR<4,float> pos     = ray.pos;
	   VECTOR<4,float> dir	   = ray.dir;
		//const __m128 box_min2 = temp.ss;
		// use whatever's apropriate to load.
		//const __m128
		//	box_min	= _mm_load_ps(box.vmin.M),
		//	box_max	= _mm_load_ps(box.vmax.M),
		//	pos	= _mm_load_ps(ray.pos.M),
		//	inv_dir = _mm_load_ps(ray.dir.M); 

		// use a div if inverted directions aren't available
		const __m128 l1 = _mm_mul_ps(_mm_sub_ps(box_min, pos), _mm_div_ps(one,dir) );
		const __m128 l2 = _mm_mul_ps(_mm_sub_ps(box_max, pos), _mm_div_ps(one,dir) );

		// the order we use for those min/max is vital to filter out
		// NaNs that happens when an inv_dir is +/- inf and
		// (box_min - pos) is 0. inf * 0 = NaN
		const __m128 filtered_l1a = _mm_min_ps(l1, plus_inf);
		const __m128 filtered_l2a = _mm_min_ps(l2, plus_inf);

		const __m128 filtered_l1b = _mm_max_ps(l1, minus_inf);
		const __m128 filtered_l2b = _mm_max_ps(l2, minus_inf);

		// now that we're back on our feet, test those slabs.
		__m128 lmax = _mm_max_ps(filtered_l1a, filtered_l2a);
		__m128 lmin = _mm_min_ps(filtered_l1b, filtered_l2b);


		// unfold back. try to hide the latency of the shufps & co.
		const __m128 lmax0 = rotatelps(lmax);
		const __m128 lmin0 = rotatelps(lmin);
		lmax = _mm_min_ss(lmax, lmax0);
		lmin = _mm_max_ss(lmin, lmin0);

		const __m128 lmax1 = _mm_movehl_ps(lmax,lmax);
		const __m128 lmin1 = _mm_movehl_ps(lmin,lmin);
		lmax = _mm_min_ss(lmax, lmax1);
		lmin = _mm_max_ss(lmin, lmin1);

#pragma warning(disable:4800)
		const int ret = _mm_comige_ss(lmax, _mm_setzero_ps()) & _mm_comige_ss(lmax,lmin);

		_mm_store_ss(t_near,lmin);
		_mm_store_ss(t_far,lmax);

		return  ret;

	}
};

#endif

/*

// ray-box
template<>
struct Intersect<RAY<4,float>,BOX<4,float>,ray_box_face_intersector>
{
	typedef BOX<4,float>    _box;
	typedef RAY<4,float>    _ray;
	typedef VECTOR<4,float> vec4f;

	static bool exec(const _ray& ray,const _box& box,ray_box_face_intersector* res) // подлежит оптимизации
	{
	
		//_ray ray2 = ray*box.matrix_inverse();
		_ray ray2 = ray;

		VECTOR<2,float> t;

		bool ret =  Intersect
					< RAY<4,float>,AABB<4,float>,VECTOR<2,float> 
					>::exec(ray2,box,&t);

		if(!ret) return false;
		
		VECTOR<4,float> rpos = ray2.pos;
		VECTOR<4,float> rdir = ray2.dir;

		VECTOR<4,float> p1 = _mm_add_ps(rpos,rdir*t.M[0]);
		VECTOR<4,float> p2 = _mm_add_ps(rpos,rdir*t.M[1]);


		//  НУЖНА ОПТИМИЗАЦИЯ!!!!!!!!!!
		if( fabsf(p1.z - box.vmin.z) < EPSILON_E5 ) 
			res->face_number[0] = _box::FRONT;
		else if(fabsf(p1.z - box.vmax.z) < EPSILON_E5 )
			res->face_number[0] = _box::BACK;

		else if(fabsf(p1.x - box.vmin.x) < EPSILON_E5 ) 
			res->face_number[0] = _box::LEFT;
		else if(fabsf(p1.x - box.vmax.x) < EPSILON_E5 ) 
			res->face_number[0] = _box::RIGHT;

		else if(fabsf(p1.y - box.vmin.y) < EPSILON_E5 ) 
			res->face_number[0] = _box::TOP;
		else 
			res->face_number[0] = _box::BOTTOM;


		if( fabsf(p2.z - box.vmin.z) < EPSILON_E5 ) 
			res->face_number[1] = _box::FRONT;
		else if(fabsf(p2.z - box.vmax.z) < EPSILON_E5 ) 
			res->face_number[1] = _box::BACK;

		else if(fabsf(p2.x - box.vmin.x) < EPSILON_E5 ) 
			res->face_number[1] = _box::LEFT;
		else if(fabsf(p2.x - box.vmax.x) < EPSILON_E5 ) 
			res->face_number[1] = _box::RIGHT;

		else if(fabsf(p2.y - box.vmin.y) < EPSILON_E5 ) 
			res->face_number[1] = _box::TOP;
		else 
			res->face_number[1] = _box::BOTTOM;

		res->flat_coord[0] = box.get_flat_coord(res->face_number[0],p1);
		res->flat_coord[1] = box.get_flat_coord(res->face_number[1],p2);
	

		return ret;
	};
};*/

///////////////////////////////////////////////////////////////////////////////////
//// ray - triangle
template<int n,class T>
struct Intersect< RAY<n,T>,
			      TRIANGLE<n,T> >
{

 typedef VECTOR<n,T>	   vec;
 typedef const VECTOR<n,T> cvec;

 static universal_call bool exec(const RAY<n,T>& ray, 
								 cvec& A, cvec& B, cvec& C,
								 vec* bar)
 {
  register vec   edge1, edge2,pvec,tvec,qvec;
  register float det,inv_det;

  edge1 = B - A;
  edge2 = C - A;

  pvec  = cross3(ray.dir,edge2);
 
  det  =  dot(edge1,pvec);
  //culling is disabled 
  // То enable cull next [if] must be 
   // if(det < EPSILON_E6M)
   if (det > -EPSILON_E6M && det < EPSILON_E6M)
     return false;
  inv_det = 1.0f / det;
 
  tvec = ray.pos - A;

  bar->M[1] = dot(tvec,pvec)*inv_det;
  if (bar->M[1] < -0.001f || bar->M[1] > 1.001f)
   return false;
 
  qvec = cross3(tvec,edge1);
  bar->M[2] = dot(ray.dir,qvec)*inv_det;
  if (bar->M[2] < -0.001f || bar->M[1] + bar->M[2] > 1.001f)
     return false;

  bar->M[0] = 1.0f - bar->M[2] - bar->M[1];
  return true;
 }

// static bool exec(const RAY<4,float>& ray, const  VECTOR<4,float> v[3], VECTOR<2,float>* bar)
// {
//  register __m128 edge1, edge2,pvec,tvec,qvec;
//  register float det,inv_det;
//
//  edge1 = _mm_sub_ps(v[1],v[0]);
//  edge2 = _mm_sub_ps(v[2],v[0]);

//  pvec  = _mm_cross3(ray.dir,edge2);
 
//  det  =  _mm_dot3(edge1,pvec);
//  //culling is disabled 
//  // То enable cull next [if] must be if(det < MGML_MATH::EPSILON_E6)
//  if (det > -MGML_MATH::EPSILON_E6 && det < MGML_MATH::EPSILON_E6)
//  if(det < MGML_MATH::EPSILON_E6)
//     return false;
//  inv_det = 1.0f / det;
 
//  tvec = _mm_sub_ps(ray.pos,v[0]);

//  bar->M[1] = _mm_dot3(tvec,pvec)*inv_det;
//  if (bar->M[1] < -0.001 || bar->M[1] > 1.001f)
//   return false;
 
//  qvec = _mm_cross3(tvec,edge1);
//  bar->M[2] = _mm_dot3(ray.dir,qvec)*inv_det;
//  if (bar->M[2] < -0.001 || bar->M[1] + bar->M[2] > 1.001f)
//     return false;

//  bar->M[0] = 1.0f - bar->M[2] - bar->M[1];
//  return true;
// }


};

///////////////////////////////////////////////////////////////////////////////////
//// box-box
template<int n,class T>
struct Intersect<AABB<n,T>,AABB<n,T> >
{
 typedef AABB<n,T>      _box;

 static universal_call bool exec(const _box& box,const _box& box2)
 {
  const int k = MGML_MATH::VECTOR<n,T>::k;
  
  for(int i=0;i<k;i++)
  {
		if(box.vmax.M[i] > box2.vmax.M[i])
		{
			if(box2.vmax.M[i] < box.vmin.M[i])
				return false;
		}
		else
		{
			if(box.vmax.M[i] < box2.vmin.M[i])
				return false;
		}
  }
  return true;
 }

};


///////////////////////////////////////////////////////////////////////////////////
//// box-triangle
template<>
struct Intersect< AABB<4,float>,TRIANGLE<4,float> >
{
 typedef AABB<4,float>     _box;
 typedef TRIANGLE<4,float> _triangle; 
 typedef VECTOR<4,float>   vec4f;

 static bool exec(const _box& box,const _triangle& tri)
 {
	 vec4f	center    = box.center();
	 vec4f half_size = (box.vmax - box.vmin)*0.5;
	 float triverts[3][3];
	 for(int i=0;i<3;i++)
	 {
		triverts[0][i] = tri.A.M[i];
		triverts[1][i] = tri.B.M[i];
		triverts[2][i] = tri.C.M[i];
	 }
#pragma warning(disable:4800)
	 return EXTERNAL_TOOLS::triBoxOverlap(center.M,half_size.M,triverts);
 }


};






///////////////////////////////////////////////////////////////////////////////////
//// box-sphere
template<int n, class T>
struct Intersect< AABB<n,T>,SPHERE<n,T> >
{
 typedef AABB<n,T>		   _box;
 typedef SPHERE<n,T>       _sphere;
 typedef VECTOR<n,T>	   vec;

 static universal_call bool exec(const _box& box,const _sphere& sph)
 {
	 const int k = MGML_MATH::VECTOR<n,T>::k;

	 bool res = true;
	 vec diff = box.center() - sph.pos;
	 vec half_size = (box.vmax - box.vmin)*0.5;
	
	 for(int i=0;i<k;i++)
	 {
		res = res & ( abs(diff.M[i]) < sph.r + half_size.M[i]);
	 }

	 return res;
 }

};


///////////////////////////////////////////////////////////////////////////////////
//// segment-plane
template<int n, class T>
struct Intersect< SEGMENT<n,T>, PLANE<n,T> >
{
	typedef SEGMENT<n,T> _line;
	typedef PLANE<n,T> _plane;
	typedef VECTOR<n,T>	_vec;

	static universal_call bool exec(const _line& line, const _plane& plane, VECTOR<n,T>* pRes)
	{
		_vec norm = line.end() - line.begin();
		norm.Normalize();

		T t = (-1)*(plane.substitute(line.begin()))/(plane.norm()*norm);

		*pRes = line.begin() + t*norm;
		return true;
	}


};


}

