#include <cassert>
#include <iostream>
#include "AxisAlignedBox3.h"
#include "Plane3.h"


using namespace std;


namespace CHCDemoEngine
{

// Overload << operator for C++-style output
ostream& operator<< (ostream &s, const AxisAlignedBox3 &A)
{
	return s << '[' << A.mMin.x << ", " << A.mMin.y << ", " << A.mMin.z << "]["
		<< A.mMax.x << ", " << A.mMax.y << ", " << A.mMax.z << ']';
}


// Overload >> operator for C++-style input
istream& operator>> (istream &s, AxisAlignedBox3 &A)
{
	char a;
	// read "[min.x, min.y, min.z][mMax.x, mMax.y, mMax.z]"
	return s >> a >> A.mMin.x >> a >> A.mMin.y >> a >> A.mMin.z >> a >> a
		>> A.mMax.x >> a >> A.mMax.y >> a >> A.mMax.z >> a;
}


AxisAlignedBox3::AxisAlignedBox3() 
{ }


AxisAlignedBox3::AxisAlignedBox3(const Vector3 &nMin, const Vector3 &nMax)
: mMin(nMin), mMax(nMax)
{}



bool AxisAlignedBox3::Unbounded() const
{
	return (mMin == Vector3(-MAXFLOAT)) || 
		(mMax == Vector3(-MAXFLOAT));
}


void AxisAlignedBox3::Include(const Vector3 &newpt)
{
	Minimize(mMin, newpt);
	Maximize(mMax, newpt);
}


void AxisAlignedBox3::Include(const AxisAlignedBox3 &bbox)
{
	Minimize(mMin, bbox.mMin);
	Maximize(mMax, bbox.mMax);
}


bool AxisAlignedBox3::IsCorrect()
{
	if ( (mMin.x > mMax.x) ||
		(mMin.y > mMax.y) ||
		(mMin.z > mMax.z) )
		return false; // box is not formed
	return true;
}


void AxisAlignedBox3::GetEdge(const int edge, Vector3 *a, Vector3 *b) const
{
	switch(edge) {
  case 0:
	  a->SetValue(mMin.x, mMin.y, mMin.z);
	  b->SetValue(mMin.x, mMin.y, mMax.z);
	  break;
  case 1:
	  a->SetValue(mMin.x, mMin.y, mMin.z);
	  b->SetValue(mMin.x, mMax.y, mMin.z);
	  break;
  case 2:
	  a->SetValue(mMin.x, mMin.y, mMin.z);
	  b->SetValue(mMax.x, mMin.y, mMin.z);
	  break;
  case 3:
	  a->SetValue(mMax.x, mMax.y, mMax.z);
	  b->SetValue(mMax.x, mMax.y, mMin.z);
	  break;
  case 4:
	  a->SetValue(mMax.x, mMax.y, mMax.z);
	  b->SetValue(mMax.x, mMin.y, mMax.z);
	  break;
  case 5:
	  a->SetValue(mMax.x, mMax.y, mMax.z);
	  b->SetValue(mMin.x, mMax.y, mMax.z);
	  break;

  case 6:
	  a->SetValue(mMin.x, mMin.y, mMax.z);
	  b->SetValue(mMin.x, mMax.y, mMax.z);
	  break;
  case 7:
	  a->SetValue(mMin.x, mMin.y, mMax.z);
	  b->SetValue(mMax.x, mMin.y, mMax.z);
	  break;
  case 8:
	  a->SetValue(mMin.x, mMax.y, mMin.z);
	  b->SetValue(mMin.x, mMax.y, mMax.z);
	  break;
  case 9:
	  a->SetValue(mMin.x, mMax.y, mMin.z);
	  b->SetValue(mMax.x, mMax.y, mMin.z);
	  break;
  case 10:
	  a->SetValue(mMax.x, mMin.y, mMin.z);
	  b->SetValue(mMax.x, mMax.y, mMin.z);
	  break;
  case 11:
	  a->SetValue(mMax.x, mMin.y, mMin.z);
	  b->SetValue(mMax.x, mMin.y, mMax.z);
	  break;
	}
}

// returns the vertex indices in the range <0..7>, v = 4.x + 2.y + z, where
// x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
void
AxisAlignedBox3::GetEdge(const int edge, int  &aIdx, int &bIdx) const
{
	switch(edge) {
  case 0: aIdx = 0; bIdx = 1; break;
  case 1: aIdx = 0; bIdx = 2; break;
  case 2: aIdx = 0; bIdx = 4; break;
  case 3: aIdx = 7; bIdx = 6; break;
  case 4: aIdx = 7; bIdx = 5; break;
  case 5: aIdx = 7; bIdx = 3; break;
  case 6: aIdx = 1; bIdx = 3; break;
  case 7: aIdx = 1; bIdx = 5; break;
  case 8: aIdx = 2; bIdx = 3; break;
  case 9: aIdx = 2; bIdx = 6; break;
  case 10: aIdx = 4; bIdx = 6; break;
  case 11: aIdx = 4; bIdx = 5; break;
	}
}

void
AxisAlignedBox3::Include(const int &axis, const float &newBound)
{
	switch (axis) {
	case 0: { // x-axis
		if (mMin.x > newBound)
			mMin.x = newBound;
		if (mMax.x < newBound)
			mMax.x = newBound;
		break;
			}
	case 1: { // y-axis
		if (mMin.y > newBound)
			mMin.y = newBound;
		if (mMax.y < newBound)
			mMax.y = newBound;
		break;
			}
	case 2: { // z-axis
		if (mMin.z > newBound)
			mMin.z = newBound;
		if (mMax.z < newBound)
			mMax.z = newBound;
		break;
			}
	}
}


void AxisAlignedBox3::Describe(ostream &app, int ind) const
{
	indent(app, ind);
	app << "AxisAlignedBox3: min at(" << mMin << "), max at(" << mMax << ")\n";
}

int
AxisAlignedBox3::IsInside(const Vector3 &v) const
{
	return ! ((v.x < mMin.x) ||
		(v.x > mMax.x) ||
		(v.y < mMin.y) ||
		(v.y > mMax.y) ||
		(v.z < mMin.z) ||
		(v.z > mMax.z));
}


bool AxisAlignedBox3::Includes(const AxisAlignedBox3 &b) const
{
	return (b.mMin.x >= mMin.x &&
		b.mMin.y >= mMin.y &&
		b.mMin.z >= mMin.z &&
		b.mMax.x <= mMax.x &&
		b.mMax.y <= mMax.y &&
		b.mMax.z <= mMax.z);
}


// compute the coordinates of one vertex of the box
Vector3 AxisAlignedBox3::GetVertex(int xAxis, int yAxis, int zAxis) const
{
	Vector3 p;
	if (xAxis)
		p.x = mMax.x;
	else
		p.x = mMin.x;

	if (yAxis)
		p.y = mMax.y;
	else
		p.y = mMin.y;

	if (zAxis)
		p.z = mMax.z;
	else
		p.z = mMin.z;
	return p;
}


// compute the vertex for number N = <0..7>, N = 4.x + 2.y + z, where
// x,y,z are either 0 or 1; (0 .. min coordinate, 1 .. max coordinate)
void AxisAlignedBox3::GetVertex(const int N, Vector3 &vertex) const
{
	switch (N) 
	{
	case 0: vertex = mMin; break;
	case 1: vertex.SetValue(mMin.x, mMin.y, mMax.z); break;
	case 2: vertex.SetValue(mMin.x, mMax.y, mMin.z); break;
	case 3: vertex.SetValue(mMin.x, mMax.y, mMax.z); break;
	case 4: vertex.SetValue(mMax.x, mMin.y, mMin.z); break;
	case 5: vertex.SetValue(mMax.x, mMin.y, mMax.z); break;
	case 6: vertex.SetValue(mMax.x, mMax.y, mMin.z); break;
	case 7: vertex = mMax; break;
	default: 
		{
			cerr << "ERROR in AxisAlignedBox3::GetVertex N=" << N <<  "\n";
		}
	}
}


float AxisAlignedBox3::SurfaceArea() const
{
	Vector3 ext = mMax - mMin;

	return 2.0f * (ext.x * ext.y +
		ext.x * ext.z +
		ext.y * ext.z);
}


float AxisAlignedBox3::GetVolume() const 
{
	return (mMax.x - mMin.x) * (mMax.y - mMin.y) * (mMax.z - mMin.z);
}


const int AxisAlignedBox3::bvertices[27][9] =
{                   // region number.. position
	{5,1,3,2,6,4,-1,-1,-1},      // 0 .. x=0 y=0 z=0
	{4,5,1,3,2,0,-1,-1,-1},      // 1 .. x=0 y=0 z=1
	{4,5,7,3,2,0,-1,-1,-1},      // 2 .. x=0 y=0 z=2

	{0,1,3,2,6,4,-1,-1,-1},      // 3 .. x=0 y=1 z=0
	{0,1,3,2,-1,-1,-1,-1,-1},    // 4 .. x=0 y=1 z=1
	{1,5,7,3,2,0,-1,-1,-1},      // 5 .. x=0 y=1 z=2

	{0,1,3,7,6,4,-1,-1,-1},      // 6 .. x=0 y=2 z=0
	{0,1,3,7,6,2,-1,-1,-1},      // 7 .. x=0 y=2 z=1
	{1,5,7,6,2,0,-1,-1,-1},      // 8 .. x=0 y=2 z=2

	// the regions number <9,17>
	{5,1,0,2,6,4,-1,-1,-1},      // 9  .. x=1 y=0 z=0
	{5,1,0,4,-1,-1,-1,-1,-1},    // 10 .. x=1 y=0 z=1
	{7,3,1,0,4,5,-1,-1,-1},      // 11 .. x=1 y=0 z=2

	{4,0,2,6,-1,-1,-1,-1,-1},    // 12 .. x=1 y=1 z=0
	{0,2,3,1,5,4,6,7,-1},        // 13 .. x=1 y=1 z=1 .. inside the box
	{1,5,7,3,-1,-1,-1,-1,-1},    // 14 .. x=1 y=1 z=2

	{4,0,2,3,7,6,-1,-1,-1},      // 15 .. x=1 y=2 z=0
	{6,2,3,7,-1,-1,-1,-1,-1},    // 16 .. x=1 y=2 z=1
	{1,5,7,6,2,3,-1,-1,-1},      // 17 .. x=1 y=2 z=2

	// the regions number <18,26>
	{1,0,2,6,7,5,-1,-1,-1},      // 18 .. x=2 y=0 z=0
	{1,0,4,6,7,5,-1,-1,-1},      // 19 .. x=2 y=0 z=1
	{0,4,6,7,3,1,-1,-1,-1},      // 20 .. x=2 y=0 z=2

	{4,0,2,6,7,5,-1,-1,-1},      // 21 .. x=2 y=1 z=0
	{5,4,6,7,-1,-1,-1,-1,-1},    // 22 .. x=2 y=1 z=1
	{5,4,6,7,3,1,-1,-1,-1},      // 23 .. x=2 y=1 z=2

	{4,0,2,3,7,5,-1,-1,-1},      // 24 .. x=2 y=2 z=0
	{5,4,6,2,3,7,-1,-1,-1},      // 25 .. x=2 y=2 z=1
	{5,4,6,2,3,1,-1,-1,-1},      // 26 .. x=2 y=2 z=2
};

// the visibility of boundary faces from a given region
// one to three triples: (axis, min_vertex, max_vertex), axis==-1(terminator)
const int AxisAlignedBox3::bfaces[27][10] =
{                              // region number .. position
	{0,0,3,1,0,5,2,0,6,-1},      // 0 .. x=0 y=0 z=0
	{0,0,3,1,0,5,-1,-1,-1,-1},   // 1 .. x=0 y=0 z=1
	{0,0,3,1,0,5,2,1,7,-1},      // 2 .. x=0 y=0 z=2

	{0,0,3,2,0,6,-1,-1,-1,-1},   // 3 .. x=0 y=1 z=0
	{0,0,3,-1,-1,-1,-1,-1,-1,-1},// 4 .. x=0 y=1 z=1
	{0,0,3,2,1,7,-1,-1,-1,-1},   // 5 .. x=0 y=1 z=2

	{0,0,3,1,2,7,2,0,6,-1},      // 6 .. x=0 y=2 z=0
	{0,0,3,1,2,7,-1,-1,-1,-1},   // 7 .. x=0 y=2 z=1
	{0,0,3,1,2,7,2,1,7,-1},      // 8 .. x=0 y=2 z=2

	// the regions number <9,17>
	{1,0,5,2,0,6,-1,-1,-1,-1},   // 9  .. x=1 y=0 z=0
	{1,0,5,-1,-1,-1,-1,-1,-1,-1},// 10 .. x=1 y=0 z=1
	{1,0,5,2,1,7,-1,-1,-1,-1},   // 11 .. x=1 y=0 z=2

	{2,0,6,-1,-1,-1,-1,-1,-1,-1},// 12 .. x=1 y=1 z=0
	{-1,-1,-1,-1,-1,-1,-1,-1,-1,-1},// 13 .. x=1 y=1 z=1 .. inside the box
	{2,1,7,-1,-1,-1,-1,-1,-1,-1},// 14 .. x=1 y=1 z=2

	{1,2,7,2,0,6,-1,-1,-1,-1},   // 15 .. x=1 y=2 z=0
	{1,2,7,-1,-1,-1,-1,-1,-1,-1},// 16 .. x=1 y=2 z=1
	{1,2,7,2,1,7,-1,-1,-1,-1},   // 17 .. x=1 y=2 z=2

	// the region number <18,26>
	{0,4,7,1,0,5,2,0,6,-1},      // 18 .. x=2 y=0 z=0
	{0,4,7,1,0,5,-1,-1,-1,-1},   // 19 .. x=2 y=0 z=1
	{0,4,7,1,0,5,2,1,7,-1},      // 20 .. x=2 y=0 z=2

	{0,4,7,2,0,6,-1,-1,-1,-1},   // 21 .. x=2 y=1 z=0
	{0,4,7,-1,-1,-1,-1,-1,-1,-1},// 22 .. x=2 y=1 z=1
	{0,4,7,2,1,7,-1,-1,-1,-1},   // 23 .. x=2 y=1 z=2

	{0,4,7,1,2,7,2,0,6,-1},      // 24 .. x=2 y=2 z=0
	{0,4,7,1,2,7,-1,-1,-1,-1},   // 25 .. x=2 y=2 z=1
	{0,4,7,1,2,7,2,1,7,-1},      // 26 .. x=2 y=2 z=2
};

// the correct corners indexing from entry face to exit face
// first index determines entry face, second index exit face, and
// the two numbers (indx, inc) determines: ind = the index on the exit
// face, when starting from the vertex 0 on entry face, 'inc' is
// the increment when we go on entry face in order 0,1,2,3 to create
// convex shaft with the rectangle on exit face. That is, inc = -1 or 1.
const int AxisAlignedBox3::pairFaceRects[6][6][2] = {
	{ // entry face = 0
		{-1,0}, // exit face 0 .. no meaning
		{0,-1}, // 1
		{0,-1}, // 2
		{0,1}, // 3 .. opposite face
		{3,1}, // 4
		{1,1} // 5
	}, 
	{ // entry face = 1
		{0,-1}, // exit face 0
		{-1,0}, // 1 .. no meaning
		{0,-1}, // 2
		{1,1}, // 3
		{0,1}, // 4 .. opposite face
		{3,1} // 5
	},
	{ // entry face = 2
		{0,-1}, // 0
		{0,-1}, // 1 
		{-1,0}, // 2 .. no meaning
		{3,1}, // 3
		{1,1}, // 4
		{0,1} // 5 .. opposite face
	},
	{ // entry face = 3
		{0,1}, // 0 .. opposite face
		{3,-1}, // 1 
		{1,1}, // 2
		{-1,0}, // 3  .. no meaning
		{0,-1}, // 4
		{0,-1} // 5
	},
	{ // entry face = 4
		{1,1}, // 0
		{0,1}, // 1 .. opposite face
		{3,1}, // 2
		{0,-1}, // 3
		{-1,0}, // 4 .. no meaning
		{0,-1} // 5
	},
	{ // entry face = 5
		{3,-1}, // 0
		{1,1}, // 1
		{0,1}, // 2  .. opposite face
		{0,-1}, // 3
		{0,-1}, // 4
		{-1,0} // 5 .. no meaning
	}
};


// ------------------------------------------------------------
// The vertices that form CLOSEST points with respect to the region
// for all the regions possible, number of regions is 3^3 = 27,
// since two parallel sides of bbox forms three disjoint spaces.
// The vertices are given in anti-clockwise order, stopped by value 15,
// at most 8 points, at least 1 point.
// The table includes the closest 1/2/4/8 points, followed possibly
// by the set of coordinates that should be used for testing for
// the proximity queries. The coordinates to be tested are described by
// the pair (a,b), when a=0, we want to test min vector of the box,
//                 when a=1, we want to test max vector of the box
//                 b=0,1,2 corresponds to the axis (0=x,1=y,2=z)
// The sequence is ended by 15, number -1 is used as the separator
// between the vertices and coordinates.
const int
AxisAlignedBox3::cvertices[27][9] = 
{                   // region number.. position
	{0,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D one vertex
	{0,1,-1,0,0,0,1,15,15},      // 1 .. x=0 y=0 z=1 D two vertices foll. by 2
	{1,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D one vertex

	{0,2,-1,0,0,0,2,15,15},      // 3 .. x=0 y=1 z=0 D
	{0,1,3,2,-1,0,0,15,15},      // 4 .. x=0 y=1 z=1 D
	{1,3,-1,0,0,1,2,15,15},      // 5 .. x=0 y=1 z=2 D

	{2,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D
	{2,3,-1,0,0,1,1,15,15},      // 7 .. x=0 y=2 z=1 D
	{3,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D

	// the regions number <9,17>
	{0,4,-1,0,1,0,2,15,15},      // 9  .. x=1 y=0 z=0 D
	{5,1,0,4,-1,0,1,15,15},      // 10 .. x=1 y=0 z=1 D
	{1,5,-1,0,1,1,2,15,15},      // 11 .. x=1 y=0 z=2 D

	{4,0,2,6,-1,0,2,15,15},      // 12 .. x=1 y=1 z=0 D
	{0,2,3,1,5,4,6,7,15},        // 13 .. x=1 y=1 z=1 .. inside the box
	{1,5,7,3,-1,1,2,15,15},      // 14 .. x=1 y=1 z=2 D

	{6,2,-1,0,2,1,1,15,15},      // 15 .. x=1 y=2 z=0 D
	{6,2,3,7,-1,1,1,15,15},      // 16 .. x=1 y=2 z=1 D
	{3,7,-1,1,1,1,2,15,15},      // 17 .. x=1 y=2 z=2 D

	// the regions number <18,26>
	{4,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D
	{4,5,-1,0,1,1,0,15,15},      // 19 .. x=2 y=0 z=1 D
	{5,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D

	{4,6,-1,0,2,1,0,15,15},      // 21 .. x=2 y=1 z=0 D
	{5,4,6,7,-1,1,0,15,15},      // 22 .. x=2 y=1 z=1 D
	{7,5,-1,1,0,1,2,15,15},      // 23 .. x=2 y=1 z=2 D

	{6,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D
	{6,7,-1,1,0,1,1,15,15},      // 25 .. x=2 y=2 z=1 D
	{7,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D
};

// Table for Sphere-AABB intersection based on the region knowledge
// Similar array to previous cvertices, but we omit the surfaces
// which are not necessary for testing. First are vertices,
// they are finished with -1. Second, there are indexes in
// the pair (a,b), when a=0, we want to test min vector of the box,
//                 when a=1, we want to test max vector of the box
//                 b=0,1,2 corresponds to the axis (0=x,1=y,2=z)
//
// So either we check the vertices or only the distance in specified
// dimensions. There are at all four possible cases:
//
//   1) we check one vertex - then sequence start with non-negative index
//      and is finished with 15
//   2) we check two coordinates of min/max vector describe by the pair
//         (a,b) .. a=min/max(0/1) b=x/y/z (0/1/2), sequence starts with 8
//      and finishes with 15
//   3) we check only one coordinate of min/max, as for 2), sequence start
//      with 9 and ends with 15
//   4) Position 13 - sphere is inside the box, intersection always exist
//        the sequence start with 15 .. no further testing is necessary
//        in this case
const int
AxisAlignedBox3::csvertices[27][6] = 
{                   // region number.. position
	{0,15,15,15,15,15},  // 0 .. x=0 y=0 z=0 D vertex only
	{8,0,0,0,1,15},      // 1 .. x=0 y=0 z=1 D two coords.
	{1,15,15,15,15,15},  // 2 .. x=0 y=0 z=2 D vertex only

	{8,0,0,0,2,15},      // 3 .. x=0 y=1 z=0 D two coords
	{9,0,0,15,15,15},    // 4 .. x=0 y=1 z=1 D one coord
	{8,0,0,1,2,15},      // 5 .. x=0 y=1 z=2 D two coords.

	{2,15,15,15,15,15},  // 6 .. x=0 y=2 z=0 D one vertex
	{8,0,0,1,1,15},      // 7 .. x=0 y=2 z=1 D two coords
	{3,15,15,15,15,15},  // 8 .. x=0 y=2 z=2 D one vertex

	// the regions number <9,17>
	{8,0,1,0,2,15},      // 9  .. x=1 y=0 z=0 D two coords
	{9,0,1,15,15,15},    // 10 .. x=1 y=0 z=1 D one coord
	{8,0,1,1,2,15},      // 11 .. x=1 y=0 z=2 D two coords

	{9,0,2,15,15,15},    // 12 .. x=1 y=1 z=0 D one coord
	{15,15,15,15,15,15}, // 13 .. x=1 y=1 z=1 inside the box, special case/value
	{9,1,2,15,15,15},    // 14 .. x=1 y=1 z=2 D one corrd

	{8,0,2,1,1,15},      // 15 .. x=1 y=2 z=0 D two coords
	{9,1,1,15,15},       // 16 .. x=1 y=2 z=1 D one coord
	{8,1,1,1,2,15},      // 17 .. x=1 y=2 z=2 D two coords

	// the regions number <18,26>
	{4,15,15,15,15,15},  // 18 .. x=2 y=0 z=0 D one vertex
	{8,0,1,1,0,15},      // 19 .. x=2 y=0 z=1 D two coords
	{5,15,15,15,15,15},  // 20 .. x=2 y=0 z=2 D one vertex

	{8,0,2,1,0,15},      // 21 .. x=2 y=1 z=0 D two coords
	{9,1,0,15,15,15},    // 22 .. x=2 y=1 z=1 D one coord
	{8,1,0,1,2,15},      // 23 .. x=2 y=1 z=2 D two coords

	{6,15,15,15,15,15},  // 24 .. x=2 y=2 z=0 D one vertex
	{8,1,0,1,1,15},      // 25 .. x=2 y=2 z=1 D two coords
	{7,15,15,15,15,15},  // 26 .. x=2 y=2 z=2 D one vertex
};


// The vertices that form all FARTHEST points with respect to the region
// for all the regions possible, number of regions is 3^3 = 27,
// since two parallel sides of bbox forms three disjoint spaces.
// The vertices are given in anti-clockwise order, stopped by value 15,
// at most 8 points, at least 1 point.
// For testing, if the AABB is whole in the sphere, it is enough
// to test only vertices, either 1,2,4, or 8.
const int
AxisAlignedBox3::fvertices[27][9] = 
{                   // region number.. position
	{7,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D
	{6,7,15,15,15,15,15,15,15},  // 1 .. x=0 y=0 z=1 D
	{6,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D

	{5,7,15,15,15,15,15,15,15},  // 3 .. x=0 y=1 z=0 D
	{4,5,7,6,15,15,15,15,15},    // 4 .. x=0 y=1 z=1 D
	{4,6,15,15,15,15,15,15,15},  // 5 .. x=0 y=1 z=2 D

	{5,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D
	{4,5,15,15,15,15,15,15,15},  // 7 .. x=0 y=2 z=1 D
	{4,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D

	// the regions number <9,17>
	{3,7,15,15,15,15,15,15,15},  // 9  .. x=1 y=0 z=0 D
	{7,3,2,6,15,15,15,15,15},    // 10 .. x=1 y=0 z=1 D
	{2,6,15,15,15,15,15,15,15},  // 11 .. x=1 y=0 z=2 D

	{5,1,3,7,15,15,15,15,15},    // 12 .. x=1 y=1 z=0 D
	{0,7,1,6,3,4,5,2,15},        // 13 .. x=1 y=1 z=1 .. inside the box
	{0,4,6,2,15,15,15,15,15},    // 14 .. x=1 y=1 z=2 D

	{5,1,15,15,15,15,15,15,15},  // 15 .. x=1 y=2 z=0 D
	{4,0,1,5,15,15,15,15,15},    // 16 .. x=1 y=2 z=1 D
	{4,0,15,15,15,15,15,15,15},  // 17 .. x=1 y=2 z=2 D

	// the regions number <18,26>
	{3,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D
	{2,3,15,15,15,15,15,15,15},  // 19 .. x=2 y=0 z=1 D
	{2,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D

	{1,3,15,15,15,15,15,15,15},  // 21 .. x=2 y=1 z=0 D
	{1,0,2,3,15,15,15,15,15},    // 22 .. x=2 y=1 z=1 D
	{2,0,15,15,15,15,15,15,15},  // 23 .. x=2 y=1 z=2 D

	{1,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D
	{0,1,15,15,15,15,15,15,15},  // 25 .. x=2 y=2 z=1 D
	{0,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D
};

// Similar table as above, farthest points, but only the ones
// necessary for testing the intersection problem. If we do
// not consider the case 13, center of the sphere is inside the
// box, then we can always test at most 2 box vertices to say whether
// the whole box is inside the sphere.
// The number of vertices is minimized using some assumptions
// about the ortogonality of vertices and sphere properties.
const int
AxisAlignedBox3::fsvertices[27][9] = 
{                   // region number.. position
	{7,15,15,15,15,15,15,15,15}, // 0 .. x=0 y=0 z=0 D 1 vertex
	{6,7,15,15,15,15,15,15,15},  // 1 .. x=0 y=0 z=1 D 2 vertices
	{6,15,15,15,15,15,15,15,15}, // 2 .. x=0 y=0 z=2 D 1 vertex

	{5,7,15,15,15,15,15,15,15},  // 3 .. x=0 y=1 z=0 D 2 vertices
	{4,7,15,5,6,15,15,15,15},    // 4 .. x=0 y=1 z=1 D 4/2 vertices
	{4,6,15,15,15,15,15,15,15},  // 5 .. x=0 y=1 z=2 D 2 vertices

	{5,15,15,15,15,15,15,15,15}, // 6 .. x=0 y=2 z=0 D 1 vertex
	{4,5,15,15,15,15,15,15,15},  // 7 .. x=0 y=2 z=1 D 2 vertices
	{4,15,15,15,15,15,15,15,15}, // 8 .. x=0 y=2 z=2 D 1 vertex

	// the regions number <9,17>
	{3,7,15,15,15,15,15,15,15},  // 9  .. x=1 y=0 z=0 D 2 vertices
	{7,2,15,3,6,15,15,15,15},    // 10 .. x=1 y=0 z=1 D 4/2 vertices
	{2,6,15,15,15,15,15,15,15},  // 11 .. x=1 y=0 z=2 D 2 vertices

	{5,3,15,1,7,15,15,15,15},    // 12 .. x=1 y=1 z=0 D 4/2 vertices
	{0,7,1,6,3,4,5,2,15},        // 13 .. x=1 y=1 z=1 .. inside the box
	{0,6,15,4,2,15,15,15,15},    // 14 .. x=1 y=1 z=2 D 4/2 vertices

	{5,1,15,15,15,15,15,15,15},  // 15 .. x=1 y=2 z=0 D 2 vertices
	{4,1,15,0,5,15,15,15,15},    // 16 .. x=1 y=2 z=1 D 4/2 vertices
	{4,0,15,15,15,15,15,15,15},  // 17 .. x=1 y=2 z=2 D 2 vertices

	// the regions number <18,26>
	{3,15,15,15,15,15,15,15,15}, // 18 .. x=2 y=0 z=0 D 1 vertex
	{2,3,15,15,15,15,15,15,15},  // 19 .. x=2 y=0 z=1 D 2 vertices
	{2,15,15,15,15,15,15,15,15}, // 20 .. x=2 y=0 z=2 D 1 vertex

	{1,3,15,15,15,15,15,15,15},  // 21 .. x=2 y=1 z=0 D 2 vertices
	{1,2,15,0,3,15,15,15,15},    // 22 .. x=2 y=1 z=1 D 4/2 vertices
	{2,0,15,15,15,15,15,15,15},  // 23 .. x=2 y=1 z=2 D 2 vertices

	{1,15,15,15,15,15,15,15,15}, // 24 .. x=2 y=2 z=0 D 1 vertex
	{0,1,15,15,15,15,15,15,15},  // 25 .. x=2 y=2 z=1 D 2 vertices
	{0,15,15,15,15,15,15,15,15}, // 26 .. x=2 y=2 z=2 D 1 vertex
};


// The fast computation of arctangent .. the maximal error is less
// than 4.1 degrees, according to Graphics GEMSII, 1991, pages 389--391
// Ron Capelli: "Fast approximation to the arctangent"
float atan22(const float& y)
{
	const float x = 1.0;
	const float c = (float)(M_PI * 0.25);

	if (y < 0.0) {
		if (y < -1.0)
			return c * (-2.0f + x / y); // for angle in <-PI/2, -PI/4)
		else
			return c * (y / x); // for angle in <-PI/4 , 0>
	}
	else {
		if (y > 1.0)
			return c * (2.0f - x / y); // for angle in <PI/4, PI/2>
		else
			return c * (y / x); // for angle in <0, PI/2>
	}
}


// This definition allows to be a point when answering true
bool
AxisAlignedBox3::IsCorrectAndNotPoint() const
{
	if ( (mMin.x > mMax.x) ||
		(mMin.y > mMax.y) ||
		(mMin.z > mMax.z) )
		return false; // box is not formed

	if ( (mMin.x == mMax.x) &&
		(mMin.y == mMax.y) &&
		(mMin.z == mMax.z) )
		return false; // degenerates to a point

	return true;
}

// This definition allows to be a point when answering true
bool AxisAlignedBox3::IsPoint() const
{
	if ( (mMin.x == mMax.x) &&
		(mMin.y == mMax.y) &&
		(mMin.z == mMax.z) )
		return true; // degenerates to a point

	return false;
}

// This definition requires shape of non-zero volume
bool AxisAlignedBox3::IsSingularOrIncorrect() const
{
	if ( (mMin.x >= mMax.x) ||
		(mMin.y >= mMax.y) ||
		(mMin.z >= mMax.z) )
		return true; // box is not formed

	return false; // has non-zero volume
}


void AxisAlignedBox3::GetSqrDistances(const Vector3 &point,
									  float &minDistance,
									  float &maxDistance) const
{
	// some overlap is possible, check the distance of vertices
	float sumMin = 0;
	float sumMax = 0;

	// Try to minimize the function of a distance
	// from the sphere center

	const float *minp = &(mMin.x);
	const float *maxp = &(mMax.x);
	const float *pcenter = &(point.x);

	// for x-axis
	for (int i = 0; i < 3; i++, minp++, maxp++, pcenter++) {
		float minsqr = sqr(*minp - *pcenter);
		float maxsqr = sqr(*maxp - *pcenter);
		if (*pcenter < *minp)
			sumMin += minsqr;
		else
			if (*pcenter > *maxp)
				sumMin += maxsqr;
		sumMax += (minsqr > maxsqr) ? minsqr : maxsqr;
	}

	minDistance = sumMin;
	maxDistance = sumMax;
}


void AxisAlignedBox3::Scale(const float scale) 
{
	Vector3 newSize = Size()*(scale*0.5f);
	Vector3 center = Center();
	mMin = center - newSize;
	mMax = center + newSize;
}


void AxisAlignedBox3::Scale(const Vector3 &scale) 
{
	Vector3 newSize = Size()*(scale*0.5f);
	Vector3 center = Center();
	mMin = center - newSize;
	mMax = center + newSize;
}


void AxisAlignedBox3::Translate(const Vector3 &shift) 
{
	mMin += shift;
	mMax += shift;
}


void AxisAlignedBox3::Initialize() 
{
	mMin = Vector3(MAXFLOAT);
	mMax = Vector3(-MAXFLOAT);
}


Vector3 AxisAlignedBox3::Center() const 
{ 
	return 0.5 * (mMin + mMax); 
}


Vector3 AxisAlignedBox3::Diagonal() const 
{ 
	return (mMax - mMin); 
}


float AxisAlignedBox3::Center(const int axis) const 
{
	return  0.5f * (mMin[axis] + mMax[axis]);
}


float AxisAlignedBox3::Min(const int axis) const 
{
	return mMin[axis];
}


float AxisAlignedBox3::Max(const int axis) const 
{
	return  mMax[axis];
}


float AxisAlignedBox3::Size(const int axis) const 
{
	return  Max(axis) - Min(axis);
}

// Read-only const access tomMin and max vectors using references
const Vector3& AxisAlignedBox3::Min() const 
{ 
	return mMin;
}


const Vector3& AxisAlignedBox3::Max() const 
{ 
	return mMax;
}


void AxisAlignedBox3::Enlarge (const Vector3 &v) 
{
	mMax += v;
	mMin -= v;
}


void AxisAlignedBox3::EnlargeToMinSize()
{
	const float epsMin = 1e-7f;
	const float epsAdd = 1e-6f;
	const float epsMul = 1.000005f;
	float dir = mMax.x - mMin.x;

	assert(dir >= 0.f);

	if (dir < epsMin) 
	{
		if (mMax.x > 0) 
		{
			mMax.x *= epsMul;
			mMax.x += epsAdd;
		}
		else 
		{
			mMin.x *= epsMul;
			mMin.x -= epsAdd;
		}
		
		assert(mMin.x < mMax.x);
	}

	dir = mMax.y - mMin.y;
	assert(dir >= 0.f);
	
	if (dir < epsMin) 
	{
		if (mMax.y > 0) 
		{
			mMax.y *= epsMul;
			mMax.y += epsAdd;
		}
		else 
		{
			mMin.y *= epsMul;
			mMin.y -= epsAdd;
		}
		assert(mMin.y < mMax.y);
	}
	
	dir = mMax.z - mMin.z;
	assert(dir >= 0.f);
	
	if (dir < epsMin) 
	{
		if (mMax.z > 0) 
		{
			mMax.z *= epsMul;
			mMax.z += epsAdd;
		}
		else 
		{
			mMin.z *= epsMul;
			mMin.z -= epsAdd;
		}
		
		assert(mMin.z < mMax.z);
	}
}


void AxisAlignedBox3::SetMin(const Vector3 &v) 
{
	mMin = v;
}


void AxisAlignedBox3::SetMax(const Vector3 &v) 
{
	mMax = v;
}


void AxisAlignedBox3::SetMin(int axis, const float value) 
{
	mMin[axis] = value;
}


void AxisAlignedBox3::SetMax(int axis, const float value) 
{
	mMax[axis] = value;
}


// Decrease box by given splitting plane
void AxisAlignedBox3::Reduce(int axis, int right, float value) 
{
	if ((value >=mMin[axis]) && (value <= mMax[axis]))
	{
		if (right)
			mMin[axis] = value;
		else
			mMax[axis] = value;
	}
}


Vector3 AxisAlignedBox3::Size() const 
{ 
	return mMax - mMin; 
}


Vector3 AxisAlignedBox3::GetRandomPoint() const
{
	Vector3 size = Size();

	return mMin + Vector3(RandomValue(0.0f, 1.0f),
		                  RandomValue(0.0f, 1.0f),
						  RandomValue(0.0f, 1.0f)) * Size();
}


bool AxisAlignedBox3::Intersects(const Vector3 &lStart, 
								 const Vector3 &lEnd) const
{
	float st, et, fst = 0, fet = 1;
	float const *bmin = &mMin.x;
	float const *bmax = &mMax.x;
	float const *si = &lStart.x;
	float const *ei = &lEnd.x;

	for (int i = 0; i < 3; ++ i) 
	{
		if (*si < *ei) 
		{
			if (*si > *bmax || *ei < *bmin)
				return false;

			const float di = *ei - *si;

			st = (*si < *bmin)? (*bmin - *si) / di : 0;
			et = (*ei > *bmax)? (*bmax - *si) / di : 1;
		}
		else 
		{
			if (*ei > *bmax || *si < *bmin)
				return false;

			const float di = *ei - *si;

			st = (*si > *bmax)? (*bmax - *si) / di : 0;
			et = (*ei < *bmin)? (*bmin - *si) / di : 1;
		}

		if (st > fst) fst = st;
		if (et < fet) fet = et;
		if (fet < fst)
			return false;

		++ bmin; ++ bmax;
		++ si; ++ ei;
	}

	//*time = fst;
	return true;
}


Vector3 AxisAlignedBox3::GetVertex(const int N) const 
{
	Vector3 v;
	GetVertex(N, v);

	return v;
}


float AxisAlignedBox3::GetExtent(int face) const 
{
	if (face < 3)
		return mMin[face];
	else
		return mMax[face - 3];
}


int AxisAlignedBox3::GetIndexNearestVertex(const Vector3 &vecPlaneNormal)
{
	int x, y, z;

	x = (vecPlaneNormal.x > 0.0f) ? 0 : 1;
	y = (vecPlaneNormal.y > 0.0f) ? 0 : 1;
	z = (vecPlaneNormal.z > 0.0f) ? 0 : 1;
	
	return 4 * x + 2 * y + z;
}


int AxisAlignedBox3::GetIndexFarthestVertex(const Vector3 &vecPlaneNormal)
{
	int x, y, z;

	x = (vecPlaneNormal.x <= 0.0f) ? 0 : 1;
	y = (vecPlaneNormal.y <= 0.0f) ? 0 : 1;
	z = (vecPlaneNormal.z <= 0.0f) ? 0 : 1;
	
	return 4 * x + 2 * y + z;
}


float AxisAlignedBox3::GetDistance(int index, const Plane3 &plane) const
{
	return plane.Distance(GetVertex(index));
}


float AxisAlignedBox3::GetMinVisibleDistance(const Plane3 &near) const
{
	return	near.mNormal.x * ((near.mNormal.x > 0.0f) ? mMin.x : mMax.x) +
			near.mNormal.y * ((near.mNormal.y > 0.0f) ? mMin.y : mMax.y) +
			near.mNormal.z * ((near.mNormal.z > 0.0f) ? mMin.z : mMax.z) +
			near.mD;
}

}