//////////////////////////////////////////////////////////////////////
//
//	Crytek Common Source code
//	
//	File:Cry_Intersection.h
//	Description: Common intersection-tests
//
//	History:
//	-March 15,2003: Created by Ivo Herzeg
//
//////////////////////////////////////////////////////////////////////

#ifndef CRYINTERSECTION_H
#define CRYINTERSECTION_H

#if _MSC_VER > 1000
# pragma once
#endif


///////////////////////////////////////////////////////////////////////////////
// Forward declarations                                                      //
///////////////////////////////////////////////////////////////////////////////

struct Line;
struct Ray;
struct Lineseg;
struct Plane;

struct AABB;
struct OBB;

struct Sphere;
struct AAEllipsoid;
struct Ellipsoid;

///////////////////////////////////////////////////////////////////////////////
// Definitions                                                               //
///////////////////////////////////////////////////////////////////////////////

#define FINDMINMAX(x0,x1,x2,min,max) \
	min = max = x0;   \
	if(x1<min) min=x1;\
	if(x1>max) max=x1;\
	if(x2<min) min=x2;\
	if(x2>max) max=x2;



///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct Line
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
struct Line {

	Vec3 pointonline;
	Vec3 direction; //caution: the direction is important for any intersection test

	//default Line constructor (without initialisation)
	inline Line( void ) {}
	inline Line( const Vec3 &o, const Vec3 &d ) {  pointonline=o; direction=d; }
	inline void operator () (  const Vec3 &o, const Vec3 &d  ) {  pointonline=o; direction=d; }

	~Line( void ) {};
};



///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct Ray
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
struct Ray {

	Vec3 origin;
	Vec3 direction;

	//default Ray constructor (without initialisation)
	inline Ray( void ) {}
	inline Ray( const Vec3 &o, const Vec3 &d ) {  origin=o; direction=d; }
	inline void operator () (  const Vec3 &o, const Vec3 &d  ) {  origin=o; direction=d; }

	~Ray( void ) {};
};


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct Lineseg
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
struct Lineseg {

	Vec3 start;
	Vec3 end;

	//default Lineseg constructor (without initialisation)
	inline Lineseg( void ) {}
	inline Lineseg( const Vec3 &s, const Vec3 &e ) {  start=s; end=e; }
	inline void operator () (  const Vec3 &s, const Vec3 &e  ) {  start=s; end=e; }

	~Lineseg( void ) {};
};

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct AABB
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

struct AABB {

	Vec3 min;
	Vec3 max;

	//default AABB constructor (without initialisation)
	inline AABB( void ) {}
	inline AABB( const Vec3 &vmin, const Vec3 &vmax ) {  min=vmin; max=vmax; }
	inline void operator () (  const Vec3 &vmin, const Vec3 &vmax  ) {  min=vmin; max=vmax; }

	~AABB( void ) {};

	//! Reset Bounding box before calculating bounds.
	void Reset()
	{
		min = Vec3( 100000,100000,100000 );
		max = Vec3( -100000,-100000,-100000 );
	}

	void Add( const Vec3 &v )
	{
		min.x = __min( min.x,v.x );
		min.y = __min( min.y,v.y );
		min.z = __min( min.z,v.z );

		max.x = __max( max.x,v.x );
		max.y = __max( max.y,v.y );
		max.z = __max( max.z,v.z );
	}

	//! Check if this bounding box overlap with bounding box of sphere.
	bool IsOverlapSphereBounds( const Vec3 &pos,float radius ) const
	{
		if (pos.x > min.x && pos.x < max.x &&
			pos.y > min.y && pos.y < max.y &&
			pos.z > min.z && pos.z < max.z) 
			return true;

		if (pos.x+radius < min.x) return false;
		if (pos.y+radius < min.y) return false;
		if (pos.z+radius < min.z) return false;
		if (pos.x-radius > max.x) return false;
		if (pos.y-radius > max.y) return false;
		if (pos.z-radius > max.z) return false;
		return true;
	}

	//! Check if this bounding box contain sphere within itself.
	bool IsContainSphere( const Vec3 &pos,float radius ) const
	{
		if (pos.x-radius < min.x) return false;
		if (pos.y-radius < min.y) return false;
		if (pos.z-radius < min.z) return false;
		if (pos.x+radius > max.x) return false;
		if (pos.y+radius > max.y) return false;
		if (pos.z+radius > max.z) return false;
		return true;
	}

	//! Check if bounding box is empty (Zero volume).
	bool IsEmpty() const { return min == max; }

	// Check two bounding boxes for intersection.
	inline bool	IsIntersectBox( const AABB &b ) const
	{
		// Check for intersection on X axis.
		if ((min.x > b.max.x)||(b.min.x > max.x)) return false;
		// Check for intersection on Y axis.
		if ((min.y > b.max.y)||(b.min.y > max.y)) return false;
		// Check for intersection on Z axis.
		if ((min.z > b.max.z)||(b.min.z > max.z)) return false;

		// Boxes overlap in all 3 axises.
		return true;
	}

	//! Transforms AABB with specified matrix.
	void Transform( const Matrix44 &tm )
	{
		Vec3 m = tm.TransformPointOLD( min );
		Vec3 vx = Vec3(tm(0,0),tm(0,1),tm(0,2))*(max.x-min.x);
		Vec3 vy = Vec3(tm(1,0),tm(1,1),tm(1,2))*(max.y-min.y);
		Vec3 vz = Vec3(tm(2,0),tm(2,1),tm(2,2))*(max.z-min.z);
		min = m;
		max = m;
		if (vx.x < 0) min.x += vx.x; else max.x += vx.x;
		if (vx.y < 0) min.y += vx.y; else max.y += vx.y;
		if (vx.z < 0) min.z += vx.z; else max.z += vx.z;

		if (vy.x < 0) min.x += vy.x; else max.x += vy.x;
		if (vy.y < 0) min.y += vy.y; else max.y += vy.y;
		if (vy.z < 0) min.z += vy.z; else max.z += vy.z;

		if (vz.x < 0) min.x += vz.x; else max.x += vz.x;
		if (vz.y < 0) min.y += vz.y; else max.y += vz.y;
		if (vz.z < 0) min.z += vz.z; else max.z += vz.z;
	}
};


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct OBB
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

struct OBB {

	Matrix33 m33; //orientation vectors
	Vec3 h;				//half-length-vector
	Vec3 c;				//center of obb 

	//default OBB constructor (without initialisation)
	inline OBB( void ) {}
	inline OBB( const Matrix33& mat33, const Vec3& hlv, const Vec3& center  ) {  m33=mat33; h=hlv; c=center;  }
	inline void operator () (  const Matrix33& mat33, const Vec3& hlv, const Vec3& center  ) {  m33=mat33; h=hlv; c=center; }

	~OBB( void ) {};
};




///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct Sphere
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
struct Sphere {

	Vec3 center;
	float radius;

	//default Sphere constructor (without initialisation)
	inline Sphere( void ) {}
	inline Sphere( const Vec3 &c, const float &r ) {  center=c; radius=r; }
	inline void operator () (  const Vec3 &c, const float &r  ) {  center=c; radius=r; }

	~Sphere( void ) {};
};


///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct AAEllipsoid
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

struct AAEllipsoid {

	Vec3 center;
	Vec3 radius_vec;

	//default AAEllipsoid constructor (without initialisation)
	inline AAEllipsoid( void ) {}
	inline AAEllipsoid( const Vec3 &c, const Vec3 &rv  ) {  radius_vec=rv; center=c; }
	inline void operator () ( const Vec3 &c, const Vec3 &rv ) {  radius_vec=rv; center=c; }

	~AAEllipsoid( void ) {};
};



///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
// struct Ellipsoid
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

struct Ellipsoid {

	Matrix34 ExtensionPos;

	//default Ellipsoid constructor (without initialisation)
	inline Ellipsoid( void ) {}
	inline Ellipsoid( const Matrix34 &ep ) {  ExtensionPos=ep; }
	inline void operator () (  const Matrix34 &ep  ) {  ExtensionPos=ep; }

	~Ellipsoid( void ) {};
};






namespace Overlap {

	////////////////////////////////////////////////////////////////		
	//! check if the point is inside an AABB
	inline bool	Point_AABB(const Vec3 &p, const Vec3 &mins,const Vec3 &maxs)
	{
		if ((p.x>=mins.x && p.x<=maxs.x) && (p.y>=mins.y && p.y<=maxs.y) && (p.z>=mins.z && p.z<=maxs.z))	return (true);
		return (false);
	}


//-----------------------------------------------------------------------------------------


	//! check if a Lineseg and a Sphere overlap
	inline bool	Lineseg_Sphere(const Lineseg& ls,const Sphere& s)
	{

		float radius2=s.radius*s.radius;

		Vec3 diff;
		Vec3 AC, AB;
		Vec3 point;

		//check if one of the two edpoints of the line is inside the sphere  
		diff.x = ls.end.x-s.center.x;
		diff.y = ls.end.y-s.center.y;
		diff.z = ls.end.z-s.center.z;
		if (diff.x*diff.x+diff.y*diff.y+diff.z*diff.z <= radius2)	return true;

		AC.x = s.center.x-ls.start.x;
		AC.y = s.center.y-ls.start.y;
		AC.z = s.center.z-ls.start.z;
		if (AC.x*AC.x+AC.y*AC.y+AC.z*AC.z <= radius2)	return true;

		//check distance from the sphere to the line
		AB.x = ls.end.x-ls.start.x;
		AB.y = ls.end.y-ls.start.y;
		AB.z = ls.end.z-ls.start.z;

		float r = (AC.x*AB.x+AC.y*AB.y+AC.z*AB.z) / (AB.x*AB.x+AB.y*AB.y+AB.z*AB.z);

		//projection falls outside the line
		if (r<0 || r>1)	return false;

		//check if the distance from the line to the center of the sphere is less than radius
		point.x = ls.start.x + r*AB.x;
		point.y = ls.start.y + r*AB.y;
		point.z = ls.start.z + r*AB.z;

		if ((point.x-s.center.x)*(point.x-s.center.x) + (point.y-s.center.y)*(point.y-s.center.y) + (point.z-s.center.z)*(point.z-s.center.z) > radius2)	return false;

		return true;
	}


	/*!
	* we use the SEPARATING AXIS TEST to check if a Linesegment overlap an AABB.
	*
	* Example:
	*  bool result=Overlap::Lineseg_AABB( ls, pos,aabb );
	*
	*/
	inline bool Lineseg_AABB ( const Lineseg &ls, const Vec3 &pos, const AABB &aabb ) {

		//calculate the half-length-vectors of the AABB
		Vec3 h	=	(aabb.max-aabb.min)*0.5f;

		//"t" is the transfer-vector from one center to the other
		Vec3 t		= ((ls.start+ls.end)*0.5f - ((aabb.max+aabb.min)*0.5f+pos));

		//calculate line-direction
		Vec3 ld  =  (ls.end-ls.start)*0.5f;

		if( fabsf(t.x) > (h.x + fabsf(ld.x)) ) return 0;
		if( fabsf(t.y) > (h.y + fabsf(ld.y)) ) return 0;
		if( fabsf(t.z) > (h.z + fabsf(ld.z)) ) return 0;
		if( fabsf(t.z*ld.y-t.y*ld.z) > (fabsf(h.y*ld.z) + fabsf(h.z*ld.y)) )	return 0;
		if( fabsf(t.x*ld.z-t.z*ld.x) > (fabsf(h.x*ld.z) + fabsf(h.z*ld.x)) )	return 0;
		if( fabsf(t.y*ld.x-t.x*ld.y) > (fabsf(h.x*ld.y) + fabsf(h.y*ld.x)) )	return 0;

		//no separating axis found, we have an intersection
		return 1;
	}



	/*!
	* we use the SEPARATING AXIS TEST to check if two OBB's overlap.
	*
	* Example:
	*  bool result=Overlap::Lineseg_OBB( lineseg, pos,obb );
	*
	*/
	inline bool Lineseg_OBB ( const Lineseg &ls, const Vec3 &pos, const OBB &obb ) {

		//------ convert OBB into half-length OBB -----------

		//calculate the half-length-vectors of the OBB
//		Vec3 h	=	(obb.aabb.max-obb.aabb.min)*0.5f;

		//the center of OBB is relative to the PIVOT
//		Vec3 c	=	(obb.aabb.max+obb.aabb.min)*0.5f;

		//the new center-position of Lineseg and OBB in world-space
		Vec3 wposobb	=	obb.m33*obb.c + pos;
		Vec3 wposls	=	(ls.start+ls.end)*0.5f;

		//"t" is the transfer-vector from one center to the other
		Vec3 t		= Matrix33::GetTransposed33(obb.m33)*(wposls - wposobb);

		//calculate line-direction in local obb-space
		Vec3 ld  =  (Matrix33::GetTransposed33(obb.m33)*(ls.end-ls.start))*0.5f;
		if( fabsf(t.x) > (obb.h.x + fabsf(ld.x)) ) return 0;
		if( fabsf(t.y) > (obb.h.y + fabsf(ld.y)) ) return 0;
		if( fabsf(t.z) > (obb.h.z + fabsf(ld.z)) ) return 0;
		if( fabsf(t.z*ld.y-t.y*ld.z) > (fabsf(obb.h.y*ld.z) + fabsf(obb.h.z*ld.y)) ) return 0;
		if( fabsf(t.x*ld.z-t.z*ld.x) > (fabsf(obb.h.x*ld.z) + fabsf(obb.h.z*ld.x)) ) return 0;
		if( fabsf(t.y*ld.x-t.x*ld.y) > (fabsf(obb.h.x*ld.y) + fabsf(obb.h.y*ld.x)) ) return 0;

		//no separating axis found, we have an intersection
		return 1;
	}





	/*----------------------------------------------------------------------------------
	* Sphere_AABB
	*	Sphere and AABB are assumed to be in the same space
	*
	* Example:
	*  bool result=Overlap::Sphere_AABB_Inside( sphere, aabb );
	*
	* 0 = no intersection             
	* 1 = intersection                
	*----------------------------------------------------------------------------------*/
	ILINE bool Sphere_AABB( const Sphere &s, const AABB &aabb ) {

		//we are using Arvo's method, to check if the objects are intersecting 
		float quatradius = s.radius * s.radius;
		Vec3 quat(0,0,0);

		if(s.center.x < aabb.min.x) { quat.x = s.center.x - aabb.min.x;	} 
		else if(s.center.x > aabb.max.x) {	quat.x = s.center.x - aabb.max.x;}
		if(s.center.y < aabb.min.y) { quat.y = s.center.y - aabb.min.y;	}	
		else if(s.center.y > aabb.max.y) {	quat.y = s.center.y - aabb.max.y;	}
		if(s.center.z < aabb.min.z) { quat.z=s.center.z-aabb.min.z;	} 
		else if(s.center.z > aabb.max.z) {	quat.z = s.center.z - aabb.max.z; }

		//return((quat.x * quat.x + quat.y * quat.y + quat.z * quat.z) < quatradius);
		return( (quat|quat) < quatradius);
	}

	/*!
	*
	* conventional method to check if a Sphere and an AABB overlap, 
	* or if the Sphere is completely inside the AABB.
	*	Sphere and AABB are assumed to be in the same space
	*
	* Example:
	*  bool result=Overlap::Sphere_AABB_Inside( sphere, aabb );
	*
	* return values:
	* 0x00 = no intersection         
	* 0x01 = Sphere and AABB overlap
	* 0x02 = Sphere in inside AABB 
	*/
	ILINE char Sphere_AABB_Inside(  const Sphere &s, const AABB& aabb ) {
		if ( Sphere_AABB(s,aabb) ) {
			Vec3 amin=aabb.min-s.center;
			Vec3 amax=aabb.max-s.center;
			if  (amin.x>=(-s.radius)) return 1;
			if  (amin.y>=(-s.radius)) return 1; 
			if 	(amin.z>=(-s.radius)) return 1; 
			if  (amax.x<=(+s.radius)) return 1;
			if  (amax.y<=(+s.radius)) return 1;
			if  (amax.z<=(+s.radius)) return 1;
			//yes, its inside 
			return 2;
		}
		return 0;
	}

	//----------------------------------------------------------------------------------
	//  Sphere_OBB
	//	VERY IMPORTANT: Sphere is assumed to be in the space of the OBB, otherwise it won't work 
	//
	//--- 0 = no intersection                                ---------------------------
	//--- 1 = intersection                                             -----------------
	//----------------------------------------------------------------------------------
	//----------------------------------------------------------------------------------
	//  Sphere_OBB
	//	VERY IMPORTANT: Sphere is assumed to be in the space of the OBB, otherwise it won't work 
	//
	//--- 0 = no intersection                                ---------------------------
	//--- 1 = intersection                                             -----------------
	//----------------------------------------------------------------------------------
	inline bool Sphere_OBB( const Sphere &s, const OBB &obb ) {

		//first we transform the sphere-center into the AABB-space of the OBB
		Vec3 SphereInOBBSpace	=	Matrix33::GetTransposed33(obb.m33) * s.center;

		//the rest ist the same as the "SphereAABBIntersection" calculation
		float quatradius = s.radius * s.radius;
		Vec3 quat(0,0,0);
		AABB aabb=AABB(obb.c-obb.h,obb.c+obb.h);
		if(SphereInOBBSpace.x < aabb.min.x) { quat.x = SphereInOBBSpace.x - aabb.min.x;	} 
		else if(SphereInOBBSpace.x > aabb.max.x) {	quat.x = SphereInOBBSpace.x - aabb.max.x;}
		if(SphereInOBBSpace.y < aabb.min.y) { quat.y = SphereInOBBSpace.y - aabb.min.y;	}	
		else if(SphereInOBBSpace.y > aabb.max.y) {	quat.y = SphereInOBBSpace.y - aabb.max.y;	}
		if(SphereInOBBSpace.z < aabb.min.z) { quat.z=SphereInOBBSpace.z-aabb.min.z;	} 
		else if(SphereInOBBSpace.z > aabb.max.z) {	quat.z = SphereInOBBSpace.z - aabb.max.z; }

		return((quat|quat) < quatradius);
	}


	/*!
	*
	* we use the SEPARATING-AXIS-TEST for OBB/Plane intersection.
	*
	* Example:
	*  bool result=Overlap::OBB_Plane( pos,obb, plane );
	*
	*/
	inline bool OBB_Plane( const Vec3 &pos, const OBB &obb, const Plane &plane ) {

		//calculate the half-length-vectors
		//Vec3 h	=	(obb.aabb.max-obb.aabb.min)*0.5f;

		//the new center-position in world-space
		Vec3 p	=	obb.m33*obb.c + pos;

		//extract the orientation-vectors from the columns of the 3x3 matrix
		//and scale them by the half-lengths
		Vec3 ax	= Vec3(obb.m33(0,0), obb.m33(1,0), obb.m33(2,0))*obb.h.x;
		Vec3 ay	= Vec3(obb.m33(0,1), obb.m33(1,1), obb.m33(2,1))*obb.h.y;
		Vec3 az	= Vec3(obb.m33(0,2), obb.m33(1,2), obb.m33(2,2))*obb.h.z;

		//check OBB against Plane, using the plane-normal as separating axis
		return	fabsf(plane|p) < (fabsf(plane.n|ax) + fabsf(plane.n|ay) + fabsf(plane.n|az));
	}


	/*!
	*
	* we use the SEPARATING AXIS TEST to check if a triangle intersects an AABB.
	*
	* Example:
	*  bool result=Overlap::AABB_Triangle( pos,aabb, tv0,tv1,tv2 );
	*
	*/
	inline bool AABB_Triangle ( const Vec3 &pos, const AABB &aabb,  const Vec3 &tv0, const Vec3 &tv1, const Vec3 &tv2 ) {

		//------ convert AABB into half-length AABB -----------
		Vec3 h	=	(aabb.max-aabb.min)*0.5f;	//calculate the half-length-vectors
		Vec3 c	=	(aabb.max+aabb.min)*0.5f;	//the center is relative to the PIVOT
		Vec3 p	=	c + pos;									//the new center-position in world-space

		//move everything so that the boxcenter is in (0,0,0)
		Vec3 v0	=	tv0-p;
		Vec3 v1	=	tv1-p;
		Vec3 v2	=	tv2-p;

		//compute triangle edges
		Vec3 e0	=	v1-v0;
		Vec3 e1	=	v2-v1;
		Vec3 e2	=	v0-v2;

		//--------------------------------------------------------------------------------------
		//    use SEPARATING AXIS THEOREM to test intersection between AABB and triangle
		//    cross-product(edge from triangle, {x,y,z}-direction),  this are 3x3=9 tests
		//--------------------------------------------------------------------------------------
		float min,max,p0,p1,p2,rad,fex,fey,fez;  
		fex = fabsf(e0.x);
		fey = fabsf(e0.y);
		fez = fabsf(e0.z);

		//AXISTEST_X01(e0.z, e0.y, fez, fey);
		p0 = e0.z*v0.y - e0.y*v0.z;  
		p2 = e0.z*v2.y - e0.y*v2.z;   
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} 
		rad = fez * h.y + fey * h.z; 
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y02(e0.z, e0.x, fez, fex);
		p0 = -e0.z*v0.x + e0.x*v0.z;
		p2 = -e0.z*v2.x + e0.x*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * h.x + fex * h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z12(e0.y, e0.x, fey, fex);
		p1 = e0.y*v1.x - e0.x*v1.y;
		p2 = e0.y*v2.x - e0.x*v2.y;
		if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
		rad = fey * h.x + fex * h.y;
		if(min>rad || max<-rad) return 0;

		//-----------------------------------------------

		fex = fabsf(e1.x);
		fey = fabsf(e1.y);
		fez = fabsf(e1.z);
		//AXISTEST_X01(e1.z, e1.y, fez, fey);
		p0 = e1.z*v0.y - e1.y*v0.z;
		p2 = e1.z*v2.y - e1.y*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * h.y + fey * h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y02(e1.z, e1.x, fez, fex);
		p0 = -e1.z*v0.x + e1.x*v0.z;
		p2 = -e1.z*v2.x + e1.x*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * h.x + fex * h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z0(e1.y, e1.x, fey, fex);
		p0 = e1.y*v0.x - e1.x*v0.y;
		p1 = e1.y*v1.x - e1.x*v1.y;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fey * h.x + fex * h.y;
		if(min>rad || max<-rad) return 0;

		//-----------------------------------------------

		fex = fabsf(e2.x);
		fey = fabsf(e2.y);
		fez = fabsf(e2.z);
		//AXISTEST_X2(e2.z, e2.y, fez, fey);
		p0 = e2.z*v0.y - e2.y*v0.z;
		p1 = e2.z*v1.y - e2.y*v1.z;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fez * h.y + fey * h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y1(e2.z, e2.x, fez, fex);
		p0 = -e2.z*v0.x + e2.x*v0.z;
		p1 = -e2.z*v1.x + e2.x*v1.z;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fez * h.x + fex * h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z12(e2.y, e2.x, fey, fex);
		p1 = e2.y*v1.x - e2.x*v1.y;
		p2 = e2.y*v2.x - e2.x*v2.y;
		if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
		rad = fey * h.x + fex * h.y;
		if(min>rad || max<-rad) return 0;

		//the {x,y,z}-directions (actually, since we use the AABB of the triangle we don't even need to test these) 
		//first test overlap in the {x,y,z}-directions
		//find min, max of the triangle each direction, and test for overlap in that direction -- 
		//this is equivalent to testing a minimal AABB around the triangle against the AABB
		AABB taabb;
		FINDMINMAX(v0.x, v1.x, v2.x, taabb.min.x,taabb.max.x);
		FINDMINMAX(v0.y, v1.y, v2.y, taabb.min.y,taabb.max.y);
		FINDMINMAX(v0.z, v1.z, v2.z, taabb.min.z,taabb.max.z);

		//test in X-direction
		FINDMINMAX(v0.x, v1.x, v2.x, taabb.min.x,taabb.max.x);
		if(taabb.min.x>h.x || taabb.max.x<-h.x) return 0;

		//test in Y-direction
		FINDMINMAX(v0.y, v1.y, v2.y, taabb.min.y,taabb.max.y);
		if(taabb.min.y>h.y || taabb.max.y<-h.y) return 0;

		//test in Z-direction
		FINDMINMAX(v0.z, v1.z, v2.z, taabb.min.z,taabb.max.z);
		if(taabb.min.z>h.z || taabb.max.z<-h.z) return 0;

		//test if the box intersects the plane of the triangle
		//compute plane equation of triangle: normal*x+d=0
		Plane plane=GetPlane( (e0%e1), v0);

		Vec3 vmin,vmax;
		if(plane.n.x>0.0f) {  vmin.x=-h.x; vmax.x=+h.x; } 
		else { vmin.x=+h.x; vmax.x=-h.x; }
		if(plane.n.y>0.0f) {  vmin.y=-h.y; vmax.y=+h.y; } 
		else { vmin.y=+h.y; vmax.y=-h.y; }
		if(plane.n.z>0.0f) {	vmin.z=-h.z; vmax.z=+h.z; } 
		else {	vmin.z=+h.z;	vmax.z=-h.z;	}
		if( (plane|vmin) > 0.0f) return 0;
		if( (plane|vmax) < 0.0f) return 0;
		return 1;
	}



	/*!
	*
	* we use the SEPARATING AXIS TEST to check if a triangle intersects an OBB.
	*
	* Example:
	*  bool result=Overlap::OBB_Trinagle( pos1,obb1, tv0,tv1,tv2 );
	*
	*/
	inline bool OBB_Triangle( const Vec3 &pos, const OBB &obb,  const Vec3 &tv0, const Vec3 &tv1, const Vec3 &tv2 ) {

		//------ convert AABB into half-length AABB -----------
		//Vec3 h	=	(obb.aabb.max-obb.aabb.min)*0.5f;	//calculate the half-length-vectors
		//Vec3 c	=	(obb.aabb.max+obb.aabb.min)*0.5f;	//the center is relative to the PIVOT
		Vec3 p	=	obb.m33*obb.c + pos;									//the new center-position in world-space

		Matrix33 transpose=Matrix33::GetTransposed33(obb.m33);

		//move everything so that the boxcenter is in (0,0,0)
		Vec3 v0	=	transpose*(tv0-p);
		Vec3 v1	=	transpose*(tv1-p);
		Vec3 v2	=	transpose*(tv2-p);

		//compute triangle edges
		Vec3 e0	=	v1-v0;
		Vec3 e1	=	v2-v1;
		Vec3 e2	=	v0-v2;

		//--------------------------------------------------------------------------------------
		//    use SEPARATING AXIS THEOREM to test intersection between AABB and triangle
		//    cross-product(edge from triangle, {x,y,z}-direction),  this are 3x3=9 tests
		//--------------------------------------------------------------------------------------
		float min,max,p0,p1,p2,rad,fex,fey,fez;  
		fex = fabsf(e0.x);
		fey = fabsf(e0.y);
		fez = fabsf(e0.z);

		//AXISTEST_X01(e0.z, e0.y, fez, fey);
		p0 = e0.z*v0.y - e0.y*v0.z;  
		p2 = e0.z*v2.y - e0.y*v2.z;   
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;} 
		rad = fez * obb.h.y + fey * obb.h.z; 
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y02(e0.z, e0.x, fez, fex);
		p0 = -e0.z*v0.x + e0.x*v0.z;
		p2 = -e0.z*v2.x + e0.x*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * obb.h.x + fex * obb.h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z12(e0.y, e0.x, fey, fex);
		p1 = e0.y*v1.x - e0.x*v1.y;
		p2 = e0.y*v2.x - e0.x*v2.y;
		if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
		rad = fey * obb.h.x + fex * obb.h.y;
		if(min>rad || max<-rad) return 0;

		//-----------------------------------------------

		fex = fabsf(e1.x);
		fey = fabsf(e1.y);
		fez = fabsf(e1.z);
		//AXISTEST_X01(e1.z, e1.y, fez, fey);
		p0 = e1.z*v0.y - e1.y*v0.z;
		p2 = e1.z*v2.y - e1.y*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * obb.h.y + fey * obb.h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y02(e1.z, e1.x, fez, fex);
		p0 = -e1.z*v0.x + e1.x*v0.z;
		p2 = -e1.z*v2.x + e1.x*v2.z;
		if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
		rad = fez * obb.h.x + fex * obb.h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z0(e1.y, e1.x, fey, fex);
		p0 = e1.y*v0.x - e1.x*v0.y;
		p1 = e1.y*v1.x - e1.x*v1.y;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fey * obb.h.x + fex * obb.h.y;
		if(min>rad || max<-rad) return 0;

		//-----------------------------------------------

		fex = fabsf(e2.x);
		fey = fabsf(e2.y);
		fez = fabsf(e2.z);
		//AXISTEST_X2(e2.z, e2.y, fez, fey);
		p0 = e2.z*v0.y - e2.y*v0.z;
		p1 = e2.z*v1.y - e2.y*v1.z;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fez * obb.h.y + fey * obb.h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Y1(e2.z, e2.x, fez, fex);
		p0 = -e2.z*v0.x + e2.x*v0.z;
		p1 = -e2.z*v1.x + e2.x*v1.z;
		if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
		rad = fez * obb.h.x + fex * obb.h.z;
		if(min>rad || max<-rad) return 0;

		//AXISTEST_Z12(e2.y, e2.x, fey, fex);
		p1 = e2.y*v1.x - e2.x*v1.y;
		p2 = e2.y*v2.x - e2.x*v2.y;
		if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
		rad = fey * obb.h.x + fex * obb.h.y;
		if(min>rad || max<-rad) return 0;

		//the {x,y,z}-directions (actually, since we use the AABB of the triangle we don't even need to test these) 
		//first test overlap in the {x,y,z}-directions
		//find min, max of the triangle each direction, and test for overlap in that direction -- 
		//this is equivalent to testing a minimal AABB around the triangle against the AABB
		AABB taabb;
		FINDMINMAX(v0.x, v1.x, v2.x, taabb.min.x,taabb.max.x);
		FINDMINMAX(v0.y, v1.y, v2.y, taabb.min.y,taabb.max.y);
		FINDMINMAX(v0.z, v1.z, v2.z, taabb.min.z,taabb.max.z);

		// test in X-direction
		FINDMINMAX(v0.x, v1.x, v2.x, taabb.min.x,taabb.max.x);
		if(taabb.min.x>obb.h.x || taabb.max.x<-obb.h.x) return 0;

		// test in Y-direction
		FINDMINMAX(v0.y, v1.y, v2.y, taabb.min.y,taabb.max.y);
		if(taabb.min.y>obb.h.y || taabb.max.y<-obb.h.y) return 0;

		// test in Z-direction
		FINDMINMAX(v0.z, v1.z, v2.z, taabb.min.z,taabb.max.z);
		if(taabb.min.z>obb.h.z || taabb.max.z<-obb.h.z) return 0;

		//test if the box intersects the plane of the triangle
		//compute plane equation of triangle: normal*x+d=0
		Plane plane=GetPlane( (e0%e1), v0);

		Vec3 vmin,vmax;
		if(plane.n.x>0.0f) {  vmin.x=-obb.h.x; vmax.x=+obb.h.x; } 
		else { vmin.x=+obb.h.x; vmax.x=-obb.h.x; }
		if(plane.n.y>0.0f) {  vmin.y=-obb.h.y; vmax.y=+obb.h.y; } 
		else { vmin.y=+obb.h.y; vmax.y=-obb.h.y; }
		if(plane.n.z>0.0f) {	vmin.z=-obb.h.z; vmax.z=+obb.h.z; } 
		else {	vmin.z=+obb.h.z;	vmax.z=-obb.h.z;	}
		if( (plane|vmin) > 0.0f) return 0;
		if( (plane|vmax) < 0.0f) return 0;
		return 1;
	}









	/*!
	*
	* conventional method to check if two AABB's overlap.
	* both AABBs are assumed to be in the same space
	*
	* Example:
	*  bool result=Overlap::AABB_AABB( aabb1, aabb2 );
	*
	*/
	ILINE bool AABB_AABB( const AABB &aabb1, const AABB &aabb2 ) {
		if  (aabb1.min.x>=aabb2.max.x) return 0;
		if  (aabb1.min.y>=aabb2.max.y) return 0; 
		if 	(aabb1.min.z>=aabb2.max.z) return 0; 
		if  (aabb1.max.x<=aabb2.min.x) return 0;
		if  (aabb1.max.y<=aabb2.min.y) return 0;
		if  (aabb1.max.z<=aabb2.min.z) return 0; 
		return 1; //the aabb's overlap
	}

	/*!
	*
	* Conventional method to check if two AABB's overlap.
	* Both AABBs are in local object space. Used the position-vector 
	* to translate them into world-space
	*
	* Example:
	*  bool result=Overlap::AABB_AABB( pos1,aabb1, pos2,aabb2 );
	*
	*/
	ILINE bool AABB_AABB( const Vec3 &pos1,const AABB &aabb1,  const Vec3 &pos2,const AABB &aabb2 ) {
		AABB waabb1(aabb1.min+pos1,aabb1.max+pos1);
		AABB waabb2(aabb2.min+pos2,aabb2.max+pos2);
		return AABB_AABB( waabb1, waabb2 );
	}


	/*!
	*
	* conventional method to check if two AABB's overlap 
	* or if AABB1 is comletely inside AABB2.
	* both AABBs are assumed to be in the same space
	*
	* Example:
	*  bool result=Overlap::AABB_AABB_Inside( aabb1, aabb2 );
	*
	* return values:
	* 0x00 = no intersection         
	* 0x01 = both AABBs one overlap
	* 0x02 = AABB1 in inside AABB2 
	*/
	ILINE char AABB_AABB_Inside( const AABB& aabb1, const AABB& aabb2 ) {
		if ( AABB_AABB(aabb1,aabb2) ) {
			if  (aabb1.min.x<=aabb2.min.x) return 1;
			if  (aabb1.min.y<=aabb2.min.y) return 1; 
			if 	(aabb1.min.z<=aabb2.min.z) return 1; 
			if  (aabb1.max.x>=aabb2.max.x) return 1;
			if  (aabb1.max.y>=aabb2.max.y) return 1;
			if  (aabb1.max.z>=aabb2.max.z) return 1;
			//yes, its inside 
			return 2;
		}
		return 0;
	}


	/*!
	*
	* we use the SEPARATING AXIS TEST to check if two OBB's overlap.
	*
	* Example:
	*  bool result=Overlap::OBB_OBB( pos1,obb1, pos2,obb2 );
	*
	*/
	inline bool OBB_OBB ( const Vec3 &pos1,const OBB &obb1,  const Vec3 &pos2,const OBB &obb2 ) {

		//tranform obb2 in local space of obb1
		Matrix33 M=obb1.m33.T()*obb2.m33;

		//the new center-position in world-space
		Vec3 p1	=	obb1.m33*obb1.c + pos1;
		Vec3 p2	=	obb2.m33*obb2.c + pos2;

		//"t" is the transfer-vector from one center to the other
		Vec3 t		= (p2-p1)*obb1.m33;

		float ra,rb;

		//--------------------------------------------------------------------------
		//--  we use the vectors "1,0,0","0,1,0" and "0,0,1" as separating axis
		//--------------------------------------------------------------------------
		rb	=	fabsf(M(0,0)*obb2.h.x) + fabsf(M(0,1)*obb2.h.y) + fabsf(M(0,2)*obb2.h.z);
		if(   fabsf(t.x) > (fabsf(obb1.h.x)+rb) ) return 0;
		rb	=	fabsf(M(1,0)*obb2.h.x) + fabsf(M(1,1)*obb2.h.y) + fabsf(M(1,2)*obb2.h.z);
		if(   fabsf(t.y) > (fabsf(obb1.h.y)+rb) ) return 0;
		rb	=	fabsf(M(2,0)*obb2.h.x) + fabsf(M(2,1)*obb2.h.y) + fabsf(M(2,2)*obb2.h.z);
		if(   fabsf(t.z) > (fabsf(obb1.h.z)+rb) ) return 0;

		//--------------------------------------------------------------------------
		//--  we use the orientation-vectors "Mx","My" and "Mz" as separating axis
		//--------------------------------------------------------------------------
		ra	= fabsf(M(0,0)*obb1.h.x) + fabsf(M(1,0)*obb1.h.y) + fabsf(M(2,0)*obb1.h.z);
		if( fabsf(t|Vec3(M(0,0),M(1,0),M(2,0))) > (ra+obb2.h.x) ) return 0;
		ra	= fabsf(M(0,1)*obb1.h.x) + fabsf(M(1,1)*obb1.h.y) + fabsf(M(2,1)*obb1.h.z);
		if( fabsf(t|Vec3(M(0,1),M(1,1),M(2,1))) > (ra+obb2.h.y) ) return 0;
		ra	= fabsf(M(0,2)*obb1.h.x) + fabsf(M(1,2)*obb1.h.y) + fabsf(M(2,2)*obb1.h.z);
		if( fabsf(t|Vec3(M(0,2),M(1,2),M(2,2))) > (ra+obb2.h.z) ) return 0;

		//---------------------------------------------------------------------
		//----  using 9 cross products we generate new separating axis
		//---------------------------------------------------------------------
		ra = obb1.h.y*fabsf(M(2,0)) + obb1.h.z*fabsf(M(1,0));
		rb = obb2.h.y*fabsf(M(0,2)) + obb2.h.z*fabsf(M(0,1));
		if( fabsf(t.z*M(1,0)-t.y*M(2,0)) > (ra+rb) ) return 0;
		ra = obb1.h.y*fabsf(M(2,1)) + obb1.h.z*fabsf(M(1,1));
		rb = obb2.h.x*fabsf(M(0,2)) + obb2.h.z*fabsf(M(0,0));
		if( fabsf(t.z*M(1,1)-t.y*M(2,1)) > (ra+rb) ) return 0;
		ra = obb1.h.y*fabsf(M(2,2)) + obb1.h.z*fabsf(M(1,2));
		rb = obb2.h.x*fabsf(M(0,1)) + obb2.h.y*fabsf(M(0,0));
		if( fabsf(t.z*M(1,2)-t.y*M(2,2)) > (ra+rb) ) return 0;


		ra = obb1.h.x*fabsf(M(2,0)) + obb1.h.z*fabsf(M(0,0));
		rb = obb2.h.y*fabsf(M(1,2)) + obb2.h.z*fabsf(M(1,1));
		if( fabsf(t.x*M(2,0)-t.z*M(0,0)) > (ra+rb) ) return 0;
		ra = obb1.h.x*fabsf(M(2,1)) + obb1.h.z*fabsf(M(0,1));
		rb = obb2.h.x*fabsf(M(1,2)) + obb2.h.z*fabsf(M(1,0));
		if( fabsf(t.x*M(2,1)-t.z*M(0,1)) > (ra+rb) ) return 0;
		ra = obb1.h.x*fabsf(M(2,2)) + obb1.h.z*fabsf(M(0,2));
		rb = obb2.h.x*fabsf(M(1,1)) + obb2.h.y*fabsf(M(1,0));
		if( fabsf(t.x*M(2,2)-t.z*M(0,2)) > (ra+rb) ) return 0;


		ra = obb1.h.x*fabsf(M(1,0)) + obb1.h.y*fabsf(M(0,0));
		rb = obb2.h.y*fabsf(M(2,2)) + obb2.h.z*fabsf(M(2,1));
		if( fabsf(t.y*M(0,0)-t.x*M(1,0)) > (ra+rb) ) return 0;
		ra = obb1.h.x*fabsf(M(1,1)) + obb1.h.y*fabsf(M(0,1));
		rb = obb2.h.x*fabsf(M(2,2)) + obb2.h.z*fabsf(M(2,0));
		if( fabsf(t.y*M(0,1)-t.x*M(1,1)) > (ra+rb) ) return 0;
		ra = obb1.h.x*fabsf(M(1,2)) + obb1.h.y*fabsf(M(0,2));
		rb = obb2.h.x*fabsf(M(2,1)) + obb2.h.y*fabsf(M(2,0));
		if( fabsf(t.y*M(0,2)-t.x*M(1,2)) > (ra+rb) ) return 0;

		//no separating axis found, we have an intersection
		return 1;
	}






#define	PLANE_X			0
#define	PLANE_Y			1
#define	PLANE_Z			2
#define	PLANE_NON_AXIAL 3

	//! check if the point is inside a triangle
	inline bool	PointInTriangle(const Vec3& point, const Vec3& v0,const Vec3& v1,const Vec3& v2,const Vec3& normal)
	{

		float xt,yt;
		Vec3 nn;

		int p1,p2;

		nn = normal;
		nn.x = (float)fabs(nn.x);
		nn.y = (float)fabs(nn.y);
		nn.z = (float)fabs(nn.z);

		if ((nn.x>=nn.y) && (nn.x>=nn.z)) 
		{
			xt=point.y; yt=point.z;
			p1=PLANE_Y;p2=PLANE_Z;
		}
		else
			if ((nn.y>=nn.x) && (nn.y>=nn.z))
			{
				xt=point.x;yt=point.z;
				p1=PLANE_X;p2=PLANE_Z;
			}
			else
			{
				xt=point.x;yt=point.y;
				p1=PLANE_X;p2=PLANE_Y;
			}

			float Ax,Ay,Bx,By;
			float s;

			bool front=false;
			bool back=false;


			Ax=(v0)[p1];Bx=(v1)[p1];
			Ay=(v0)[p2];By=(v1)[p2];

			s=((Ay-yt)*(Bx-Ax)-(Ax-xt)*(By-Ay));

			if (s>=0) 
			{ 
				if (back) 
					return (false); 
				front=true;
			}
			else 
			{ 
				if (front) 
					return (false); 
				back=true; 
			}

			Ax=(v1)[p1];Bx=(v2)[p1];
			Ay=(v1)[p2];By=(v2)[p2];

			s=((Ay-yt)*(Bx-Ax)-(Ax-xt)*(By-Ay));

			if (s>=0) 
			{ 
				if (back) 
					return (false); 
				front=true;
			}
			else 
			{ 
				if (front) return (false); 
				back=true; 
			}

			Ax=(v2)[p1];Bx=(v0)[p1];
			Ay=(v2)[p2];By=(v0)[p2];

			s=((Ay-yt)*(Bx-Ax)-(Ax-xt)*(By-Ay));

			if (s>=0) 
			{ 
				if (back) 
					return (false); 
				front=true;
			}
			else 
			{ 
				if (front) 
					return (false); 
				back=true; 
			}

			return (true);
	}





}




namespace Intersect {



inline bool Ray_Plane(const Ray &ray, const Plane &plane, Vec3 &output ) {
	float numer		= plane|ray.origin;
	float cosine	=	plane.n|ray.direction;
	//REJECTION 1: if "line-direction" is perpendicular to "plane-normal", an intersection is not possible!
	//REJECTION 2: we deal with single-sided planes. 
	//             if "line-direction" is pointing in the same direction as "the plane-normal", 
	//             an intersection is not possible!
	if (cosine > 0) return 0;					//normal is orthogonal to vector, cant intersect
	output				=	ray.origin+(ray.direction*(-numer/cosine));
	//skip, if cutting-point is "behind" ray.origin
	if (((output-ray.origin)|ray.direction)<0) return 0;
	return 1;														//intersection occured	
}

inline bool Line_Plane(const Line &line, const Plane &plane, Vec3 &output ) {
	double perpdist	= plane|line.pointonline;
	double cosine		=	plane.n|line.direction;
	//REJECTION 1: if "line-direction" is perpendicular to "plane-normal", an intersection is not possible!
	//REJECTION 2: we deal with single-sided planes. 
	//             if "line-direction" is pointing in the same direction as "the plane-normal", 
	//             an intersection is not possible!
	if (cosine > 0) return 0;
	//an intersection is possible: calculate the exact point!
	float pd_c			=	(float)(-perpdist/cosine);
	output					=	line.pointonline+(line.direction*pd_c);
	return 1;														//intersection occured	
}





/*!
*
* overlap-test between a line and a triangle.
* IMPORTANT: this is a single-sided test. That means its not enough 
* that the triangle and line overlap, its also important that the triangle 
* is "visible" when you are looking along the line-direction.  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* return "true" if line and triangle overlap.
* 
*/
inline bool Line_Triangle( const Line &line, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2 ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  line.direction % edge_1; 
	//if determinat is near zero, line lies in plane of triangel 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to line origin
	Vec3 tvec=line.pointonline-v0;
	float u=tvec | pvec;
	if ( (u<0.0f) || (u>det)) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	float v= (line.direction | qvec);
	if ( (v<0.0f) || ((u+v)>det)) return 0;
	return 1;
}

/*
* calculates intersection between a line and a triangle.
* IMPORTANT: this is a single-sided intersection test. That means its not enough 
* that the triangle and line overlap, its also important that the triangle 
* is "visible" when you are looking along the line-direction.  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* if there is an intertection the functions return "true" and stores the  
* 3d-intersection point in "output". if the function returns "false" the value in 
* "output" is undefined  
* 
*/
inline bool Line_Triangle( const Line &line, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2, Vec3 &output ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  line.direction % edge_1; 
	//if determinat is near zero, ray lies in plane of triangel 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to ray origin
	Vec3 tvec=line.pointonline-v0;
	//calculate U parameter and test bounds
	float u=tvec | pvec;
	if (u<0.0f || u>det) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	//calculate V parameter and test bounds
	float v= (line.direction | qvec);
	if ( v<0.0f || (u+v)>det) return 0;
	//------------------------------------------------------
	//we have an intersection and now we can calculate t
	float t = (edge_1 | qvec) / det;
	//we use t als a scale parameter, to get the 3D-intersection point
	output = (line.direction*t)+line.pointonline;
	return 1;
}






/*!
*
* overlap-test between a ray and a triangle.
* IMPORTANT: this is a single-sided test. That means its not sufficient 
* that the triangle and ray overlap, its also important that the triangle 
* is "visible" when you are looking from the origin along the ray-direction.  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* return "true" if ray and triangle overlap.
*/
inline bool Ray_Triangle( const Ray &ray, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2 ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  ray.direction % edge_1; 
	//if determinat is near zero, ray lies in plane of triangle 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to ray origin
	Vec3 tvec=ray.origin-v0;
	//calculate U parameter and test bounds
	float u=tvec | pvec;
	if ( u<0.0f || u>det) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	//calculate V parameter and test bounds
	float v= (ray.direction | qvec);
	if ( v<0.0f || (u+v)>det) return 0;
	//------------------------------------------------------
	//We have an intersection and now we can calculate t
	float t = (edge_1 | qvec) / det;
	//we use t als a scale parameter, to get the 3D-intersection point
	Vec3 output = (ray.direction*t)+ray.origin;
	//skip, if cutting-point is "behind" ray.origin
	if (((output-ray.origin)|ray.direction)<0) return 0;
	return 1;
}

/*
* calculates intersection between a ray and a triangle.
* IMPORTANT: this is a single-sided intersection test. That means its not sufficient 
* that the triangle and rayt overlap, its also important that the triangle 
* is "visible" when you from the origin along the ray-direction.  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* if there is an intertection the functions return "true" and stores the  
* 3d-intersection point in "output". if the function returns "false" the value in 
* "output" is undefined  
*/
inline bool Ray_Triangle( const Ray &ray, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2, Vec3 &output ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  ray.direction % edge_1; 
	//if determinat is near zero, ray lies in plane of triangle 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to ray origin
	Vec3 tvec=ray.origin-v0;
	//calculate U parameter and test bounds
	float u=tvec | pvec;
	if (u<0.0f || u>det) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	//calculate V parameter and test bounds
	float v= (ray.direction | qvec);
	if ( v<0.0f || (u+v)>det) return 0;
	//------------------------------------------------------
	//We have an intersection and we can calculate t
	float t = (edge_1 | qvec) / det;
	//we use t als a scale parameter, to get the 3D-intersection point
	output = (ray.direction*t)+ray.origin;
	//skip, if cutting-point is "behind" ray.origin
	if (((output-ray.origin)|ray.direction)<0) return 0;
	return 1;
}







/*!
*
* overlap-test between line-segment and a triangle.
* IMPORTANT: this is a single-sided test. That means its not sufficient 
* that the triangle and line-segment overlap, its also important that the triangle 
* is "visible" when you are looking along the linesegment from "start" to "end".  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* return "true" if linesegment and triangle overlap.
*/
inline bool Lineseg_Triangle( const Lineseg &lineseg, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2 ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//direction of lineseg. normalizing is not necessary
	Vec3 direction=(lineseg.end-lineseg.start);
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  direction % edge_1; 
	//if determinat is near zero, ray lies in plane of triangle 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to ray origin
	Vec3 tvec=lineseg.start-v0;
	//calculate U parameter and test bounds
	float u=tvec | pvec;
	if (u<0.0f || u>det) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	//calculate V parameter and test bounds
	float v= (direction | qvec);
	if ( v<0.0f || (u+v)>det) return 0;
	//------------------------------------------------------
	//If we ignore the limitations of linesegment, then we have an intersection and we can calculate t
	float t = (edge_1 | qvec) / det;
	//skip, if lineseg and triangle are not overlapping
	Vec3 output = (direction*t)+lineseg.start;
	if (((output-lineseg.start)|direction)<0) return 0;
	if (((output-lineseg.end)|direction)>0) return 0;
	return 1;
}



/*!
*
* calculates intersection between a line-segment and a triangle.
* IMPORTANT: this is a single-sided intersection test. That means its not sufficient 
* that the triangle and line-segment overlap, its also important that the triangle 
* is "visible" when you are looking along the linesegment from "start" to "end".  
* 
* If you need a double-sided test, you'll have to call this function twice with 
* reversed order of triangle vertices.   
* 
* return values 
* if there is an intertection the the functions return "true" and stores the  
* 3d-intersection point in "output". if the function returns "false" the value in 
* "output" is undefined  
* 
*/
inline bool Lineseg_Triangle( const Lineseg &lineseg, const Vec3 &v0, const Vec3 &v1, const Vec3 &v2, Vec3 &output ) {
	//find vectors for two edges sharing v0
	Vec3 edge_1 = v1-v0;
	Vec3 edge_2 = v2-v0;
	//direction of lineseg. normalizing is not necessary
	Vec3 direction=(lineseg.end-lineseg.start);
	//begin calculating determinant - also used to calculate U parameter
	Vec3 pvec  =  direction % edge_1; 
	//if determinat is near zero, ray lies in plane of triangle 
	float det = edge_2 | pvec;
	if (det<=0) return 0;
	//calculate distance from v0 to ray origin
	Vec3 tvec=lineseg.start-v0;
	//calculate U parameter and test bounds
	float u=tvec | pvec;
	if (u<0.0f || u>det) return 0;
	//prepare to test V parameter
	Vec3 qvec=tvec % edge_2;
	//calculate V parameter and test bounds
	float v= (direction | qvec);
	if ( v<0.0f || (u+v)>det) return 0;
	//------------------------------------------------------
	//If we ignore the limitations of linesegment, then we have an intersection and we calcultate t
	float t = (edge_1 | qvec) / det;
	//we use t als a scale parameter, to get the 3D-intersection point
	output = (direction*t)+lineseg.start;
	//skip, if lineseg and triangle are not overlapping
	if (((output-lineseg.start)|direction)<0) return 0;
	if (((output-lineseg.end)|direction)>0) return 0;
	return 1;
}





//----------------------------------------------------------------------------------
//--- 0x00 = no intersection                               --------------------------
//--- 0x01 = not possible   --
//--- 0x02 = not possible  --
//--- 0x03 = two intersection, lineseg has ENTRY and EXIT point  --
//----------------------------------------------------------------------------------

inline unsigned char Line_Sphere( const Line &line, const Sphere &s, Vec3 &i0, Vec3 &i1 ) {

	Vec3 end=line.pointonline+line.direction;

	float a = line.direction|line.direction;
	float b = (line.direction|(line.pointonline-s.center))*2.0f;
	float c = ((line.pointonline-s.center)|(line.pointonline-s.center)) - (s.radius*s.radius); 

	float desc = (b*b) - (4 * a *c); 

	unsigned char intersection=0;
	if (desc >= 0.0f)
	{
		float lamba0 = (-b - cry_sqrtf(desc)) / (2.0f * a); 
		//_stprintf(d3dApp.token,"lamba0: %20.12f",lamba0);
		//d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
		i0 = line.pointonline + ((end-line.pointonline)*lamba0); 
		intersection=1;

		float lamba1 = (-b + cry_sqrtf(desc))/(2.0f*a); 
		//_stprintf(d3dApp.token,"lamba1: %20.12f",lamba1);
		//d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
		i1 = line.pointonline + ((end-line.pointonline)*lamba1); 
		intersection|=2;
	}

	return intersection; 
}



//----------------------------------------------------------------------------------
//--- 0x00 = no intersection                               --------------------------
//--- 0x01 = not possible   --
//--- 0x02 = one intersection, lineseg has just an EXIT point but no ENTRY point (ls.start is inside the sphere)  --
//--- 0x03 = two intersection, lineseg has ENTRY and EXIT point  --
//----------------------------------------------------------------------------------

inline unsigned char Ray_Sphere( const Ray &ray, const Sphere &s, Vec3 &i0, Vec3 &i1 ) {
		Vec3 end=ray.origin+ray.direction;
		float a = ray.direction|ray.direction;
		float b = (ray.direction|(ray.origin-s.center))*2.0f;
		float c = ((ray.origin-s.center)|(ray.origin-s.center)) - (s.radius*s.radius); 

		float desc = (b*b) - (4 * a *c); 
				
		unsigned char intersection=0;
		if (desc >= 0.0f)
		{
			float lamba0 = (-b - cry_sqrtf(desc)) / (2.0f * a); 
			//	_stprintf(d3dApp.token,"lamba0: %20.12f",lamba0);
			//	d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
			if (lamba0>0.0f)	{
				i0 = ray.origin + ((end-ray.origin)*lamba0); 
				intersection=1;
			}

			float lamba1 = (-b + cry_sqrtf(desc))/(2.0f*a); 
			//	_stprintf(d3dApp.token,"lamba1: %20.12f",lamba1);
			//	d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
			if (lamba1>0.0f)	{
				i1 = ray.origin + ((end-ray.origin)*lamba1); 
				intersection|=2;
			}
	}
	return intersection; 
}

inline bool Ray_SphereFirst( const Ray &ray, const Sphere &s, Vec3 &intPoint )
{
	Vec3 p2;
	unsigned char res = Ray_Sphere( ray,s,intPoint,p2 );
	if (res == 2) {	intPoint = p2; }
	if (res > 1)	return true;
	return false;
}



//----------------------------------------------------------------------------------
//--- 0x00 = no intersection                               --------------------------
//--- 0x01 = one intersection, lineseg has just an ENTRY point but no EXIT point (ls.end is inside the sphere)  --
//--- 0x02 = one intersection, lineseg has just an EXIT point but no ENTRY point (ls.start is inside the sphere)  --
//--- 0x03 = two intersection, lineseg has ENTRY and EXIT point  --
//----------------------------------------------------------------------------------
inline unsigned char Lineseg_Sphere( const Lineseg &ls, const Sphere &s, Vec3 &i0, Vec3 &i1 ) {

		Vec3 dir = (ls.end - ls.start);

		float a = dir|dir;
		float b = (dir|(ls.start-s.center))*2.0f;
		float c = ((ls.start-s.center)|(ls.start-s.center)) - (s.radius*s.radius); 
		float desc = (b*b) - (4 * a *c); 
				
		unsigned char intersection=0;
		if (desc >= 0.0f)	{

		float lamba0 = (-b - cry_sqrtf(desc)) / (2.0f * a); 
		//_stprintf(d3dApp.token,"lamba0: %20.12f",lamba0);
		//d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
		if (lamba0 >0.0f)	{
			i0 = ls.start + ((ls.end-ls.start)*lamba0); 
			//skip, if 1st cutting-point is "in front" of ls.end
			if (((i0-ls.end)|dir)>0) return 0;
			intersection=0x01;
		}
					
		float lamba1 = (-b + cry_sqrtf(desc)) / (2.0f * a); 
		//_stprintf(d3dApp.token,"lamba1: %20.12f",lamba1);
		//d3dApp.m_pFont->DrawText( 2, d3dApp.PrintY, D3DCOLOR_ARGB(255,255,255,0), d3dApp.token );	d3dApp.PrintY+=20;
		if ( lamba1 > 0.0f)	{
			i1 = ls.start + ((ls.end-ls.start)*lamba1); 
			//skip, if 2nd cutting-point is "in front" of ls.end (=ls.end is inside sphere)
			if (((i1-ls.end)|dir)>0) return intersection;
			intersection|=0x02;
		}
	}
	return intersection; 
}



// line/sphere-intersection
// p1: 1st point of line
// p2: 2nd point of line
// p3: center of sphere
// r: radius of sphere
// i1: 1st intersection point
// i2: 2nd intersection point
// return number of intersections
inline int Lineseg_Sphere(Vec3 p1, Vec3 p2, Vec3 p3, float r, Vec3 &i1, Vec3 i2)
{
	float dx=p2.x-p1.x;
	float dy=p2.y-p1.y;
	float dz=p2.z-p1.z;
	float a=dx*dx+dy*dy+dz*dz;
	float b=2.0f*(dx*(p1.x-p3.x)+dy*(p1.y-p3.y)+dz*(p1.z-p3.z));
	float c=p3.Dot(p3)+p2.Dot(p2)-2.0f*(p3.x*p1.x+p3.y*p1.y+p3.z-p1.z)-r*r;
	float d=b*b-4.0f*a*c;
	if (d<0.0f)
		return 0;
	float u;
	u=(-b+(float)sqrt((double)d))/a*a;
	i1=p1+((p2-p1)*u);
	if (d==0.0)
		return 1;
	u=(-b-(float)sqrt((double)d))/a*a;
	i2=p1+((p2-p1)*u);
	return 2;
}


}; //CIntersect


















/////////////////////////////////////////////////////////////////////////
//this is some special engine stuff, should be moved to a better location
/////////////////////////////////////////////////////////////////////////

// for bbox's checks and calculations
#define MAX_BB	+99999.0f
#define MIN_BB	-99999.0f

//! checks if this has been set to minBB
inline bool IsMinBB( const Vec3& v ) {
	if (v.x<=MIN_BB) return (true);
	if (v.y<=MIN_BB) return (true);
	if (v.z<=MIN_BB) return (true);
	return (false);
}

//! checks if this has been set to maxBB
inline bool IsMaxBB( const Vec3& v ) {
	if (v.x>=MAX_BB) return (true);
	if (v.y>=MAX_BB) return (true);
	if (v.z>=MAX_BB) return (true);
	return (false);
}

inline Vec3	SetMaxBB( void ) { return Vec3(MAX_BB,MAX_BB,MAX_BB); }
inline Vec3	SetMinBB( void ) { return Vec3(MIN_BB,MIN_BB,MIN_BB); }





inline void AddToBounds (const Vec3& v, Vec3& mins, Vec3& maxs) {
	if (v.x < mins.x)	mins.x = v.x;
	if (v.x > maxs.x)	maxs.x = v.x;
	if (v.y < mins.y)	mins.y = v.y;
	if (v.y > maxs.y)	maxs.y = v.y;
	if (v.z < mins.z)	mins.z = v.z;
	if (v.z > maxs.z)	maxs.z = v.z;
}

////////////////////////////////////////////////////////////////		
//! calc the area of a polygon giving a list of vertices and normal
inline float CalcArea(const Vec3 *vertices,int numvertices,const Vec3 &normal)
{	
	Vec3 csum(0,0,0);

	int n=numvertices;
	for (int i = 0, j = 1; i <= n-2; i++, j++)
	{
		csum.x += vertices[i].y*vertices[j].z-
			vertices[i].z*vertices[j].y;
		csum.y += vertices[i].z*vertices[j].x-
			vertices[i].x*vertices[j].z;
		csum.z += vertices[i].x*vertices[j].y-
			vertices[i].y*vertices[j].x;
	}

	csum.x += vertices[n-1].y*vertices[0].z-
		vertices[n-1].z*vertices[0].y;
	csum.y += vertices[n-1].z*vertices[0].x-
		vertices[n-1].x*vertices[0].z;
	csum.z += vertices[n-1].x*vertices[0].y-
		vertices[n-1].y*vertices[0].x;

	float area=0.5f*(float)fabs(normal*csum);
	return (area);
}



#endif //vector
