#include <iostream>
#include "AxisAlignedBox3.h"
#include "Ray.h"
#include "Polygon3.h"
#include "Mesh.h"
#include "Halton.h"
#include "Triangle3.h"
#include "VssRay.h"



#include <cassert>


using namespace std;

namespace GtpVisibilityPreprocessor {

#define FATAL Debug
#define FATAL_ABORT exit(1)


// AxisAlignedBox3 implementations

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

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 Polygon3 &newpoly)
{
	VertexContainer::const_iterator it, it_end = newpoly.mVertices.end();

	for (it = newpoly.mVertices.begin(); it != it_end; ++ it)
		Include(*it);
}

void
AxisAlignedBox3::Include(const PolygonContainer &polys)
{
	PolygonContainer::const_iterator it, it_end = polys.end();

	for (it = polys.begin(); it != it_end; ++ it)
		Include(*(*it));
}

void
AxisAlignedBox3::Include(Mesh *mesh)
{
	VertexContainer::const_iterator it, it_end = mesh->mVertices.end();

	for (it = mesh->mVertices.begin(); it != it_end; ++ it)
		Include(*it);
}


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

#if 0
// ComputeMinMaxT computes the minimum and maximum signed distances
// of intersection with the ray; it returns 1 if the ray hits
// the bounding box and 0 if it does not.
int
AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
				float *tmin,
				float *tmax) const
{
  float minx, maxx, miny, maxy, minz, maxz;

  if (fabs(ray.dir.x) < 0.001) {
    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.x - ray.loc.x) / ray.dir.x;
    float t2 = (mMax.x - ray.loc.x) / ray.dir.x;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
    if (maxx < 0)
      return 0;
  }

  if (fabs(ray.dir.y) < 0.001) {
    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
      miny = -MAXFLOAT;
      maxy = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.y - ray.loc.y) / ray.dir.y;
    float t2 = (mMax.y - ray.loc.y) / ray.dir.y;
    if (t1 < t2) {
      miny = t1;
      maxy = t2;
    }
    else {
      miny = t2;
      maxy = t1;
    }
    if (maxy < 0.0)
      return 0;
  }

  if (fabs(ray.dir.z) < 0.001) {
    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
      minz = -MAXFLOAT;
      maxz = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.z - ray.loc.z) / ray.dir.z;
    float t2 = (mMax.z - ray.loc.z) / ray.dir.z;
    if (t1 < t2) {
      minz = t1;
      maxz = t2;
    }
    else {
      minz = t2;
      maxz = t1;
    }
    if (maxz < 0.0)
      return 0;
  }

  *tmin = minx;
  if (miny > *tmin)
    *tmin = miny;
  if (minz > *tmin)
    *tmin = minz;

  *tmax = maxx;
  if (maxy < *tmax)
    *tmax = maxy;
  if (maxz < *tmax)
    *tmax = maxz;

  return 1; // yes, intersection was found
}
#else
// another variant of the same, with less variables
int
AxisAlignedBox3::ComputeMinMaxT(const Ray &ray,
				float *tmin,
				float *tmax) const
{
  const float dirEps = 1e-8f;
  register float minx, maxx;
  
  if (fabs(ray.dir.x) < dirEps) {
    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.x - ray.loc.x) * ray.invDir.x;
    float t2 = (mMax.x - ray.loc.x) * ray.invDir.x;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  *tmin = minx;
  *tmax = maxx;
  
  if (fabs(ray.dir.y) < dirEps) {
    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.y - ray.loc.y) * ray.invDir.y;
    float t2 = (mMax.y - ray.loc.y) * ray.invDir.y;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  if (minx > *tmin)
    *tmin = minx;
  if (maxx < *tmax)
    *tmax = maxx;
  
  if (fabs(ray.dir.z) < dirEps) {
    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.z - ray.loc.z) * ray.invDir.z;
    float t2 = (mMax.z - ray.loc.z) * ray.invDir.z;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  if (minx > *tmin)
    *tmin = minx;
  if (maxx < *tmax)
    *tmax = maxx;

  return 1; // yes, intersection was found
}
#endif

// ComputeMinMaxT computes the minimum and maximum parameters
// of intersection with the ray; it returns 1 if the ray hits
// the bounding box and 0 if it does not.
int AxisAlignedBox3::ComputeMinMaxT(const Ray &ray, float *tmin, float *tmax,
				    EFaces &entryFace, EFaces &exitFace) const
{
  float minx, maxx, miny, maxy, minz, maxz;
  int swapped[3];
  
  if (fabs(ray.dir.x) < 0.001) {
    if (mMin.x < ray.loc.x && mMax.x > ray.loc.x) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.x - ray.loc.x) / ray.dir.x;
    float t2 = (mMax.x - ray.loc.x) / ray.dir.x;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
      swapped[0] = 0;
    }
    else {
      minx = t2;
      maxx = t1;
      swapped[0] = 1;
    }
    if (maxx < 0)
      return 0;
  }

  if (fabs(ray.dir.y) < 0.001) {
    if (mMin.y < ray.loc.y && mMax.y > ray.loc.y) {
      miny = -MAXFLOAT;
      maxy = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.y - ray.loc.y) / ray.dir.y;
    float t2 = (mMax.y - ray.loc.y) / ray.dir.y;
    if (t1 < t2) {
      miny = t1;
      maxy = t2;
      swapped[1] = 0;
    }
    else {
      miny = t2;
      maxy = t1;
      swapped[1] = 1;
    }
    if (maxy < 0.0)
      return 0;
  }

  if (fabs(ray.dir.z) < 0.001) {
    if (mMin.z < ray.loc.z && mMax.z > ray.loc.z) {
      minz = -MAXFLOAT;
      maxz = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.z - ray.loc.z) / ray.dir.z;
    float t2 = (mMax.z - ray.loc.z) / ray.dir.z;
    if (t1 < t2) {
      minz = t1;
      maxz = t2;
      swapped[2] = 0;
    }
    else {
      minz = t2;
      maxz = t1;
      swapped[2] = 1;
    }
    if (maxz < 0.0)
      return 0;
  }

  *tmin = minx;
  entryFace = ID_Back;
  if (miny > *tmin) {
    *tmin = miny;
    entryFace = ID_Left;
  }
  if (minz > *tmin) {
    *tmin = minz;
    entryFace = ID_Bottom;
  }
  
  *tmax = maxx;
  exitFace = ID_Back;
  if (maxy < *tmax) {
    *tmax = maxy;
    exitFace = ID_Left;
  }
  if (maxz < *tmax) {
    *tmax = maxz;
    exitFace = ID_Bottom;
  }

  if (swapped[entryFace])
    entryFace = (EFaces)(entryFace + 3);
  
  if (!swapped[exitFace])
    exitFace = (EFaces)(exitFace + 3);
  
  return 1; // yes, intersection was found
}

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

// computes the passing through parameters for case tmin<tmax and tmax>0
int
AxisAlignedBox3::GetMinMaxT(const Ray &ray, float *tmin, float *tmax) const
{
  if (!ComputeMinMaxT(ray, tmin, tmax))
    return 0;
  if ( *tmax < *tmin)
    return 0; // the ray passes outside the box

  if ( *tmax < 0.0)
    return 0; // the intersection is not on the positive halfline

  return 1; // ray hits the box .. origin can be outside or inside the box
}

// computes the signed distances for case tmin<tmax and tmax>0
int AxisAlignedBox3::GetMinMaxT(const Ray &ray, 
								float *tmin, 
								float *tmax,
								EFaces &entryFace,
								EFaces &exitFace) const
{
  if (!ComputeMinMaxT(ray, tmin, tmax, entryFace, exitFace))
    return 0;
  if ( *tmax < *tmin)
    return 0; // the ray passes outside the box
  
  if ( *tmax < 0.0)
    return 0; // the intersection is not on the positive halfline
  
  return 1; // ray hits the box .. origin can be outside or inside the box
}

#if 0
int
AxisAlignedBox3::IsInside(const Vector3 &point) const
{
  return (point.x >= mMin.x) && (point.x <= mMax.x) &&
         (point.y >= mMin.y) && (point.y <= mMax.y) &&
         (point.z >= mMin.z) && (point.z <= mMax.z);
}
#else
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));
}
#endif

int AxisAlignedBox3::IsInside(const VssRay &v) const
{
	return IsInside(v.mTermination) && IsInside(v.mOrigin);
}


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: {
      FATAL << "ERROR in AxisAlignedBox3::GetVertex N=" << N <<  "\n";
      FATAL_ABORT;
    }
  }
}

// Returns region 0 .. 26 ; R = 9*x + 3*y + z ; (x,y,z) \in {0,1,2}
int
AxisAlignedBox3::GetRegionID(const Vector3 &point) const
{
  int ID = 0;

  if (point.z >= mMin.z) {
    if (point.z <= mMax.z)
      ID += 1; // inside the two boundary planes
    else
      ID += 2; // outside
  }

  if (point.y >= mMin.y) {
    if (point.y <= mMax.y)
      ID += 3; // inside the two boundary planes
    else
      ID += 6; // outside
  }

  if (point.x >= mMin.x) {
    if (point.x <= mMax.x)
      ID += 9; // inside the two boundary planes
    else
      ID += 18; // outside
  }
  return ID;
}



// computes if a given box (smaller) lies at least in one
// projection whole in the box (larger) = this encompasses given
// ;-----;
// | ;-; |
// | '-' |
// '-----'
int
AxisAlignedBox3::IsPiercedByBox(const AxisAlignedBox3 &box, int &axis) const
{
  // test on x-axis
  if ( (mMax.x < box.mMin.x) ||
       (mMin.x > box.mMax.x) )
    return 0; // the boxes do not overlap at all at x-axis
  if ( (box.mMin.y > mMin.y) &&
       (box.mMax.y < mMax.y) &&
       (box.mMin.z > mMin.z) &&
       (box.mMax.z < mMax.z) ) {
    axis = 0;
    return 1; // the boxes overlap in x-axis
  }
  // test on y-axis
  if ( (mMax.y < box.mMin.y) ||
       (mMin.y > box.mMax.y) )
    return 0; // the boxes do not overlap at all at y-axis
  if ( (box.mMin.x > mMin.x) &&
       (box.mMax.x < mMax.x) &&
       (box.mMin.z > mMin.z) &&
       (box.mMax.z < mMax.z) ) {
    axis = 1;
    return 1; // the boxes overlap in y-axis
  }
  // test on z-axis
  if ( (mMax.z < box.mMin.z) ||
       (mMin.z > box.mMax.z) )
    return 0; // the boxes do not overlap at all at y-axis
  if ( (box.mMin.x > mMin.x) &&
       (box.mMax.x < mMax.x) &&
       (box.mMin.y > mMin.y) &&
       (box.mMax.y < mMax.y) ) {
    axis = 2;
    return 1; // the boxes overlap in z-axis
  }
  return 0;
}

float
AxisAlignedBox3::SurfaceArea() const
{
  Vector3 ext = mMax - mMin;
 
  return 2.0f * (ext.x * ext.y +
	        ext.x * ext.z +
	        ext.y * ext.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>
  }
}


float
AxisAlignedBox3::ProjectToSphereSA(const Vector3 &viewpoint, int *tcase) const
{
  int id = GetRegionID(viewpoint);
  *tcase = id;

  // spherical projection .. SA represents solid angle
  if (id == 13) // .. inside the box
    return (float)(4.0*M_PI); // the whole sphere
  float SA = 0.0; // inital value
    
  int i = 0; // the pointer in the array of vertices
  while (bfaces[id][i] >= 0) {
    int axisO = bfaces[id][i++];
    int minvIdx = bfaces[id][i++];
    int maxvIdx = bfaces[id][i++];
    Vector3 vmin, vmax;
    GetVertex(minvIdx, vmin);
    GetVertex(maxvIdx, vmax);
    float h = fabs(vmin[axisO] - viewpoint[axisO]);
    int axis = (axisO + 1) % 3; // next axis 
    float a = (vmin[axis] - viewpoint[axis]) / h; // minimum for v-range
    float b = (vmax[axis] - viewpoint[axis]) / h; // maximum for v-range
    //if (a > b) {
    //  FATAL << "ProjectToSphereSA::Error a > b\n";
    //  FATAL_ABORT;
    //}
    //if (vmin[axisO] != vmax[axisO]) {
    //  FATAL << "ProjectToSphereSA::Error a-axis != b-axis\n";
    //  FATAL_ABORT;     
    //}
    axis = (axisO + 2) % 3; // next second axis       
    float c = (vmin[axis] - viewpoint[axis]) / h; // minimum for u-range
    float d = (vmax[axis] - viewpoint[axis]) / h; // maximum for u-range
    //if (c > d) {
    //  FATAL << "ProjectToSphereSA::Error c > d\n";
    //  FATAL_ABORT;
    //}
    SA +=atan22(d*b/sqrt(b*b + d*d + 1.0f)) - atan22(b*c/sqrt(b*b + c*c + 1.0f))
      - atan22(d*a/sqrt(a*a + d*d + 1.0f)) + atan22(a*c/sqrt(a*a + c*c + 1.0f));
  }

#if 0
  if ((SA > 2.0*M_PI) ||
      (SA < 0.0)) {
    FATAL << "The solid angle has strange value: ";
    FATAL << "SA = "<< SA << endl;
    FATAL_ABORT;
  }
#endif
  
  return SA;
}

// Projects the box to a plane given a normal vector only and
// computes the surface area of the projected silhouette
// no clipping of the box is performed.
float
AxisAlignedBox3::ProjectToPlaneSA(const Vector3 &normal) const
{
  Vector3 size = Size();

  // the surface area of the box to a yz-plane - perpendicular to x-axis
  float sax = size.y * size.z;

  // the surface area of the box to a zx-plane - perpendicular to y-axis
  float say = size.z * size.x;

  // the surface area of the box to a xy-plane - perpendicular to z-axis
  float saz = size.x * size.y;
  
  return sax * fabs(normal.x) + say * fabs(normal.y) + saz * fabs(normal.z);
}



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

// returns true, when the sphere specified by the origin and radius
// fully contains the box
bool
AxisAlignedBox3::IsFullyContainedInSphere(const Vector3 &center, float rad) const
{
  int region = GetRegionID(center);
  float rad2 = rad*rad;

  // vertex of the box
  Vector3 vertex;

  int i = 0;
  for (i = 0 ; ; i++) {
    int a = fsvertices[region][i];
    if (a == 15)
      return true; // if was not false untill now, it must be contained

    assert( (a>=0) && (a<8) );
    
    // normal vertex
    GetVertex(a, vertex);

    if (SqrMagnitude(vertex - center) > rad2)
      return false;    
  } // for
  
}

// returns true, when the volume of the sphere and volume of the
// axis aligned box has no intersection
bool
AxisAlignedBox3::HasNoIntersectionWithSphere(const Vector3 &center, float rad) const
{
  int region = GetRegionID(center);
  float rad2 = rad*rad;

  // vertex of the box
  Vector3 vertex;

  switch (csvertices[region][0]) {
  case 8: {
    // test two coordinates described within the field
    int face = 3*csvertices[region][1] + csvertices[region][2];
    float dist = GetExtent(face) - center[csvertices[region][2]];
    dist *= dist;
    face = 3 * (csvertices[region][3]) + csvertices[region][4];
    float dist2 = GetExtent(face) - center[csvertices[region][4]];
    dist += (dist2 * dist2);
    if (dist > rad2)
      return true; // no intersection is possible
  }
  case 9: {
    // test one coordinate described within the field
    int face = 3*csvertices[region][1] + csvertices[region][2];
    float dist = fabs(GetExtent(face) - center[csvertices[region][2]]);
    if (dist > rad)
      return true; // no intersection is possible
  }
  case 15:
    return false; // box and sphere surely has intersection
  default: {
    // test using normal vertices 
    assert( (csvertices[region][0]>=0) && (csvertices[region][0]<8) );

    // normal vertex
    GetVertex(csvertices[region][0], vertex);

    if (SqrMagnitude(vertex - center) > rad2)
      return true; // no intersectino is possible    
    }
  } // switch

  return false; // partial or full containtment
}

#if 0
// Given the sphere, determine the mutual position between the
// sphere and box

// SOME BUG IS INSIDE !!!! V.H. 25/4/2001
int
AxisAlignedBox3::MutualPositionWithSphere(const Vector3 &center, float rad) const
{
  int region = GetRegionID(center);
  float rad2 = rad*rad;

  // vertex of the box
  Vector3 vertex;

  // first testing for  full containtment - whether sphere fully
  // contains the box
  int countInside = 0; // how many points were found inside
  
  int i = 0;
  for (i = 0 ; ; i++) {
    int a = fsvertices[region][i];
    if (a == 15)
      return 1; // the sphere fully contain the box

    assert( (a>=0) && (a<8) );
    
    // normal vertex
    GetVertex(a, vertex);

    if (SqrMagnitude(vertex - center) <= rad2)
      countInside++; // the number of vertices inside the sphere
    else {
      if (countInside)
	return 0; // partiall overlap has been found
      // the sphere does not fully contain the box .. only way to break
      // this loop and go for other testing
      break;
    }
  } // for

  // now only box and sphere can partially overlap or no intersection
  switch (csvertices[region][0]) {
  case 8: {
    // test two coordinates described within the field
    int face = 3*csvertices[region][1] + csvertices[region][2];
    float dist = GetExtent(face) - center[csvertices[region][2]];
    dist *= dist;
    face = 3 * (csvertices[region][3]) + csvertices[region][4];
    float dist2 = GetExtent(face) - center[csvertices[region][4]];
    dist += (dist2 * dist2);
    if (dist > rad2 )
      return -1; // no intersection is possible
  }
  case 9: {
    // test one coordinate described within the field
    int face = 3*csvertices[region][1] + csvertices[region][2];
    float dist = fabs(GetExtent(face) - center[csvertices[region][2]]);
    if (dist > rad)
      return -1; // no intersection is possible
  }
  case 15:
    return 0 ; // partial overlap is now guaranteed
  default: {
    // test using normal vertices 
    assert( (csvertices[region][0]>=0) && (csvertices[region][0]<8) );
    
    // normal vertex
    GetVertex(csvertices[region][0], vertex);

    if (SqrMagnitude(vertex - center) > rad2)
      return -1; // no intersection is possible    
    }
  } // switch

  return 0; // partial intersection is guaranteed
}
#else

// Some maybe smarter version, extendible easily to d-dimensional
// space!
// Given a sphere described by the center and radius,
// the fullowing function returns:
//   -1 ... the sphere and the box are completely separate
//    0 ... the sphere and the box only partially overlap
//    1 ... the sphere contains fully the box
//  Note: the case when box fully contains the sphere is not reported
//        since it was not required.
int
AxisAlignedBox3::MutualPositionWithSphere(const Vector3 &center, float rad) const
{
//#define SPEED_UP

#ifndef SPEED_UP
  // slow version, instructively written  
#if 0
  // does it make sense to test
  // checking the sides of the box for possible non-intersection
  if ( ((center.x + rad) < mMin.x) ||
       ((center.x - rad) > mMax.x) ||
       ((center.y + rad) < mMin.y) ||
       ((center.y - rad) > mMax.y) ||
       ((center.z + rad) < mMin.z) ||
       ((center.z - rad) > mMax.z) ) {
    // cout << "r ";
    return -1;  // no overlap is possible
  }
#endif
  
  // someoverlap is possible, check the distance of vertices
  rad = rad*rad;
  float sumMin = 0;
  // Try to minimize the function of a distance
  // from the sphere center

  // for x-axis
  float minSqrX = sqr(mMin.x - center.x);
  float maxSqrX = sqr(mMax.x - center.x);
  if (center.x < mMin.x)
    sumMin = minSqrX; 
  else
    if (center.x > mMax.x)
      sumMin = maxSqrX;

  // for y-axis
  float minSqrY = sqr(mMin.y - center.y);
  float maxSqrY = sqr(mMax.y - center.y);
  if (center.y < mMin.y)
    sumMin += minSqrY; 
  else
    if (center.y > mMax.y)
      sumMin += maxSqrY;

  // for z-axis
  float minSqrZ = sqr(mMin.z - center.z);
  float maxSqrZ = sqr(mMax.z - center.z);
  if (center.z < mMin.z)
    sumMin += minSqrZ; 
  else
    if (center.z > mMax.z)
      sumMin += maxSqrZ;

  if (sumMin > rad)
    return -1; // no intersection between sphere and box

  // try to find out the maximum distance between the
  // sphere center and vertices
  float sumMax = 0;
  
  if (minSqrX > maxSqrX)
    sumMax = minSqrX;
  else
    sumMax = maxSqrX;
  
  if (minSqrY > maxSqrY)
    sumMax += minSqrY;
  else
    sumMax += maxSqrY;
  
  if (minSqrZ > maxSqrZ)
    sumMax += minSqrZ;
  else
    sumMax += maxSqrZ;
  
  // sumMin < rad
  if (sumMax < rad)
    return 1; // the sphere contains the box completely

  // partial intersection, part of the box is outside the sphere
  return 0;
#else

  // Optimized version of the test
  
#ifndef __VECTOR_HACK
#error "__VECTOR_HACK for Vector3 was not defined"
#endif

  // some overlap is possible, check the distance of vertices
  rad = rad*rad;
  float sumMin = 0;
  float sumMax = 0;
  // Try to minimize the function of a distance
  // from the sphere center

  const float *minp = &(min[0]);
  const float *maxp = &(max[0]);
  const float *pcenter = &(center[0]);
  
  // 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;
  }

  if (sumMin > rad)
    return -1; // no intersection between sphere and box

  // sumMin < rad
  if (sumMax < rad)
    return 1; // the sphere contains the box completely

  // partial intersection, part of the box is outside the sphere
  return 0;
#endif
}
#endif

// Given the cube describe by the center and the half-sie,
// determine the mutual position between the cube and the box
int
AxisAlignedBox3::MutualPositionWithCube(const Vector3 &center, float radius) const
{
  // the cube is described by the center and the distance to the any face
  // along the axes

  // Note on efficiency!
  // Can be quite optimized using tables, but I do not have time
  // V.H. 18/11/2001

  AxisAlignedBox3 a = 
    AxisAlignedBox3(Vector3(center.x - radius, center.y - radius, center.z - radius),
	   Vector3(center.x + radius, center.y + radius, center.z + radius));

  if (a.Includes(*this))
    return 1; // cube contains the box

  if (OverlapS(a,*this))
    return 0; // cube partially overlap the box

  return -1; // completely separate  
}

void
AxisAlignedBox3::GetSqrDistances(const Vector3 &point,
				 float &minDistance,
				 float &maxDistance
				 ) const
{

  
#ifndef __VECTOR_HACK
#error "__VECTOR_HACK for Vector3 was not defined"
#endif

  // 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[0]);
  const float *maxp = &(mMax[0]);
  const float *pcenter = &(point[0]);
  
  // 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;
}


int
AxisAlignedBox3::Side(const Plane3 &plane) const
{
  Vector3 v;
  int i, m=3, M=-3, s;
  
  for (i=0;i<8;i++) {
    GetVertex(i, v);
    if((s = plane.Side(v)) < m)
      m=s;
    if(s > M)
      M=s;
    if (m && m==-M)
      return 0;
  }
  
  return (m == M) ? m : m + M;
}


Plane3 AxisAlignedBox3::GetPlane(const int face) const
{
	switch (face) 
	{
        
	case 0:
		return Plane3(Vector3(-1, 0, 0), mMin);
	case 1: 
		return Plane3(Vector3(1, 0, 0), mMax);
	case 2:
		return Plane3(Vector3(0, -1, 0), mMin);
	case 3:
		return Plane3(Vector3(0, 1, 0), mMax);
	case 4:
		return Plane3(Vector3(0, 0, -1), mMin);
	case 5:
		return Plane3(Vector3(0, 0, 1), mMax);
	}

	// should not come here
	return Plane3();
}


int
AxisAlignedBox3::GetFaceVisibilityMask(const Rectangle3 &rectangle) const
{
  int mask = 0;
  for (int i=0; i < 4; i++)
    mask |= GetFaceVisibilityMask(rectangle.mVertices[i]);
  return mask;
}


int
AxisAlignedBox3::GetFaceVisibilityMask(const Vector3 &position) const {
  
  // assume that we are not inside the box
  int c=0;
  
  if (position.x<(mMin.x-Limits::Small))
    c|=1;
  else
    if (position.x>(mMax.x+Limits::Small))
      c|=2;
  
  if (position.y<(mMin.y-Limits::Small))
    c|=4;
  else
    if (position.y>(mMax.y+Limits::Small))
      c|=8;
  
  if (position.z<(mMin.z-Limits::Small))
    c|=16;
  else
    if (position.z>(mMax.z+Limits::Small))
      c|=32;
  
  return c;
}


Rectangle3
AxisAlignedBox3::GetFace(const int face) const
{
  Vector3 v[4];
  switch (face) {
    
  case 0:
    v[3].SetValue(mMin.x,mMin.y,mMin.z);
    v[2].SetValue(mMin.x,mMax.y,mMin.z);
    v[1].SetValue(mMin.x,mMax.y,mMax.z);
    v[0].SetValue(mMin.x,mMin.y,mMax.z);
    break;
    
  case 1: 
    v[0].SetValue(mMax.x,mMin.y,mMin.z);
    v[1].SetValue(mMax.x,mMax.y,mMin.z);
    v[2].SetValue(mMax.x,mMax.y,mMax.z);
    v[3].SetValue(mMax.x,mMin.y,mMax.z);
    break;
    
  case 2:
    v[3].SetValue(mMin.x,mMin.y,mMin.z);
    v[2].SetValue(mMin.x,mMin.y,mMax.z);
    v[1].SetValue(mMax.x,mMin.y,mMax.z);
    v[0].SetValue(mMax.x,mMin.y,mMin.z);
    break;
  
  case 3:
    v[0].SetValue(mMin.x,mMax.y,mMin.z);
    v[1].SetValue(mMin.x,mMax.y,mMax.z);
    v[2].SetValue(mMax.x,mMax.y,mMax.z);
    v[3].SetValue(mMax.x,mMax.y,mMin.z);
    break;

  case 4:
    v[3].SetValue(mMin.x,mMin.y,mMin.z);
    v[2].SetValue(mMax.x,mMin.y,mMin.z);
    v[1].SetValue(mMax.x,mMax.y,mMin.z);
    v[0].SetValue(mMin.x,mMax.y,mMin.z);
    break;
  
  case 5:
    v[0].SetValue(mMin.x,mMin.y,mMax.z);
    v[1].SetValue(mMax.x,mMin.y,mMax.z);
    v[2].SetValue(mMax.x,mMax.y,mMax.z);
    v[3].SetValue(mMin.x,mMax.y,mMax.z);
    break;
  }
  
  return Rectangle3(v[0], v[1], v[2], v[3]);
}

struct VertexData
{
	Vector3 mVertex;
	float mAngle;
	
	VertexData(Vector3 vtx, float angle): mVertex(vtx), mAngle(angle)
	{}

	bool operator<(const VertexData &b) const 
	{
		return mAngle > b.mAngle;
	}
};

// TODO: use a table to avoid normal and distance computations
Polygon3 *AxisAlignedBox3::CrossSection(const Plane3 &plane) const
{
	Polygon3 *planePoly = new Polygon3();
	
	int side[8];
	bool onFrontSide = false, onBackSide = false;
	
	Vector3 vtx;
	
	//////////////
	//-- compute classification of vertices

	for (int i = 0; i < 8; ++i)
	{
		GetVertex(i, vtx);
		side[i] = plane.Side(vtx);
		if (side[i] > 0)
			onFrontSide = true;
		else if (side[i] < 0)
			onBackSide = true;
		else // vertex coincident => push_back
			planePoly->mVertices.push_back(vtx);
	}

	///////////
	//-- find intersections

	if (onFrontSide && onBackSide)
	{
		Vector3 ptA, ptB;
		for (int i = 0; i < 12; ++ i)
		{
			int aIdx, bIdx;
			GetEdge(i, aIdx, bIdx);

			ptA = GetVertex(aIdx);
			ptB = GetVertex(bIdx);

			int sideA = side[aIdx];
			int sideB = side[bIdx];

			if (((sideA > 0) && (sideB < 0)) || (sideA < 0) && (sideB > 0))
				planePoly->mVertices.push_back(plane.FindIntersection(ptA, ptB));	
		}
	}

	// order intersections	
	if (planePoly->mVertices.size() > 3)
	{
		Vector3 centerOfMass(0);

		int i;
		// compute center of mass
		for (i = 0; i < (int)planePoly->mVertices.size(); ++ i)
			centerOfMass += planePoly->mVertices[i];
		
		centerOfMass /= (float)planePoly->mVertices.size();

		vector<VertexData> vertexData;
		Vector3 refVec = Normalize(centerOfMass - planePoly->mVertices[0]);

		// compute angle to reference point
		for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
		{
		    float angle = 
		      Angle(refVec, centerOfMass - planePoly->mVertices[i], plane.mNormal);
		    
		    vertexData.push_back(VertexData(planePoly->mVertices[i], angle));
		}
		
		std::stable_sort(vertexData.begin(), vertexData.end());

		// update vertices
		for (i = 1; i < (int)planePoly->mVertices.size(); ++ i)
			planePoly->mVertices[i] = vertexData[i - 1].mVertex;
	}
	else if (planePoly->mVertices.size() == 3)
	{
		// fix orientation if needed
		if (DotProd(planePoly->GetNormal(), plane.mNormal) < 0)
		{
			Vector3 v = planePoly->mVertices[1];
			planePoly->mVertices[1] = planePoly->mVertices[2];
			planePoly->mVertices[2] = v;
		}
	}
	
	return planePoly;
}


bool AxisAlignedBox3::GetRaySegment(const Ray &ray, 
									float &minT, 
									float &maxT) const
{
	maxT = 1e15f;
	minT = 0;
	
	// test with bounding box
	if (!GetMinMaxT(ray, &minT, &maxT))
		return false;

	if (minT < 0)
		minT = 0;

	// bound ray or line segment
	if (//(ray.GetType() == Ray::LOCAL_RAY) && 
	    !ray.intersections.empty() && 
		(ray.intersections[0].mT <= maxT))
	{
		maxT = ray.intersections[0].mT;
	}

	return true;
}


  // Compute tmin and tmax for a ray, whenever required .. need not pierce box
int
AxisAlignedBox3::ComputeMinMaxT(const Vector3 &origin,
				const Vector3 &direction,
				float *tmin,
				float *tmax) const
{
  register float minx, maxx;

  
  Vector3 invDirection;
  const float eps = 1e-6f;
  const float invEps = 1e6f;
  
  // it does change the ray direction very slightly,
  // but the size direction vector is not practically changed
  
  if (fabs(direction.x) < eps) {
    if (direction.x < 0.0f)
      invDirection.x = -invEps;
    else
      invDirection.x = invEps;
  }
  else
    invDirection.x = 1.0f / direction.x;
  
  if (fabs(direction.y) < eps) {
    if (direction.y < 0.0)
      invDirection.y = -invEps;
    else
      invDirection.y = invEps;
  }
  else
    invDirection.y = 1.0f / direction.y;
  
  if (fabs(direction.z) < eps) {
    if (direction.z < 0.0f)
      invDirection.z = -invEps;
    else
      invDirection.z = invEps;
  }
  else
    invDirection.z = 1.0f / direction.z;


  
  if (fabs(direction.x) < 0.001f) {
    if (mMin.x < origin.x && mMax.x > origin.x) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.x - origin.x) * invDirection.x;
    float t2 = (mMax.x - origin.x) * invDirection.x;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  *tmin = minx;
  *tmax = maxx;
  
  if (fabs(direction.y) < 0.001) {
    if (mMin.y < origin.y && mMax.y > origin.y) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.y - origin.y) * invDirection.y;
    float t2 = (mMax.y - origin.y) * invDirection.y;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  if (minx > *tmin)
    *tmin = minx;
  if (maxx < *tmax)
    *tmax = maxx;
  
  if (fabs(direction.z) < 0.001) {
    if (mMin.z < origin.z && mMax.z > origin.z) {
      minx = -MAXFLOAT;
      maxx = MAXFLOAT;
    }
    else
      return 0;
  }
  else {
    float t1 = (mMin.z - origin.z) * invDirection.z;
    float t2 = (mMax.z - origin.z) * invDirection.z;
    if (t1 < t2) {
      minx = t1;
      maxx = t2;
    }
    else {
      minx = t2;
      maxx = t1;
    }
	//    if (maxx < 0.0)
	//      return 0;
  }

  if (minx > *tmin)
    *tmin = minx;
  if (maxx < *tmax)
    *tmax = maxx;

  return 1; // yes, intersection was found
}


bool AxisAlignedBox3::GetIntersectionFace(Rectangle3 &face, 
										  const AxisAlignedBox3 &neighbour) const
										   
{ 
	if (EpsilonEqual(mMin[0], neighbour.Max(0)))
	{
		float maxy = min(mMax.y, neighbour.mMax.y);
		float maxz = min(mMax.z, neighbour.mMax.z);
		float miny = max(mMin.y, neighbour.mMin.y);
		float minz = max(mMin.z, neighbour.mMin.z);

		face.mVertices[3].SetValue(mMin.x, miny, minz);
		face.mVertices[2].SetValue(mMin.x, maxy, minz);
		face.mVertices[1].SetValue(mMin.x, maxy, maxz);
		face.mVertices[0].SetValue(mMin.x, miny, maxz);

        return true;
	}
    if (EpsilonEqual(mMax[0], neighbour.Min(0)))
	{
		float maxy = min(mMax.y, neighbour.mMax.y);
		float maxz = min(mMax.z, neighbour.mMax.z);
		float miny = max(mMin.y, neighbour.mMin.y);
		float minz = max(mMin.z, neighbour.mMin.z);

		face.mVertices[0].SetValue(mMax.x, miny, minz);
		face.mVertices[1].SetValue(mMax.x, maxy, minz);
		face.mVertices[2].SetValue(mMax.x, maxy, maxz);
		face.mVertices[3].SetValue(mMax.x, miny, maxz);

		return true;
	}
    if (EpsilonEqual(mMin[1], neighbour.Max(1)))
	{
		float maxx = min(mMax.x, neighbour.mMax.x);
		float maxz = min(mMax.z, neighbour.mMax.z);
		float minx = max(mMin.x, neighbour.mMin.x);
		float minz = max(mMin.z, neighbour.mMin.z);

		face.mVertices[3].SetValue(minx, mMin.y, minz);
		face.mVertices[2].SetValue(minx, mMin.y, maxz);
		face.mVertices[1].SetValue(maxx, mMin.y, maxz);
		face.mVertices[0].SetValue(maxx, mMin.y, minz);
		
		return true;
	}
	if (EpsilonEqual(mMax[1], neighbour.Min(1)))
	{		
		float maxx = min(mMax.x, neighbour.mMax.x);
		float maxz = min(mMax.z, neighbour.mMax.z);
		float minx = max(mMin.x, neighbour.mMin.x);
		float minz = max(mMin.z, neighbour.mMin.z);

		face.mVertices[0].SetValue(minx, mMax.y, minz);
		face.mVertices[1].SetValue(minx, mMax.y, maxz);
		face.mVertices[2].SetValue(maxx, mMax.y, maxz);
		face.mVertices[3].SetValue(maxx, mMax.y, minz);

		return true;
	}
	if (EpsilonEqual(mMin[2], neighbour.Max(2)))
	{
		float maxx = min(mMax.x, neighbour.mMax.x);
		float maxy = min(mMax.y, neighbour.mMax.y);
		float minx = max(mMin.x, neighbour.mMin.x);
		float miny = max(mMin.y, neighbour.mMin.y);

		face.mVertices[3].SetValue(minx, miny, mMin.z);
		face.mVertices[2].SetValue(maxx, miny, mMin.z);
		face.mVertices[1].SetValue(maxx, maxy, mMin.z);
		face.mVertices[0].SetValue(minx, maxy, mMin.z);
	
		return true;
	}
	if (EpsilonEqual(mMax[2], neighbour.Min(2)))
	{
		float maxx = min(mMax.x, neighbour.mMax.x);
		float maxy = min(mMax.y, neighbour.mMax.y);
		float minx = max(mMin.x, neighbour.mMin.x);
		float miny = max(mMin.y, neighbour.mMin.y);

		face.mVertices[0].SetValue(minx, miny, mMax.z);
		face.mVertices[1].SetValue(maxx, miny, mMax.z);
		face.mVertices[2].SetValue(maxx, maxy, mMax.z);
		face.mVertices[3].SetValue(minx, maxy, mMax.z);

		return true;
	}

	return false;
}


void IncludeBoxInMesh(const AxisAlignedBox3 &box, Mesh &mesh)
{
	// add 6 vertices of the box
	int index = (int)mesh.mVertices.size();
	
	for (int i=0; i < 8; ++ i) 
	{
		Vector3 v;
		box.GetVertex(i, v);
		mesh.mVertices.push_back(v);
	}
	
	mesh.AddFace(new Face(index + 0, index + 1, index + 3, index + 2) );
	mesh.AddFace(new Face(index + 0, index + 2, index + 6, index + 4) );
	mesh.AddFace(new Face(index + 4, index + 6, index + 7, index + 5) );
	
	mesh.AddFace(new Face(index + 3, index + 1, index + 5, index + 7) );
	mesh.AddFace(new Face(index + 0, index + 4, index + 5, index + 1) );
	mesh.AddFace(new Face(index + 2, index + 3, index + 7, index + 6) );
}


void AxisAlignedBox3::ExtractPolys(PolygonContainer &polys) const
{
	Polygon3 *face1 = new Polygon3();
	polys.push_back(face1);   

    face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
	face1->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
	face1->mVertices.push_back(Vector3(mMin.x, mMax.y ,mMin.z));
	face1->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));

	Polygon3 *face2 = new Polygon3();  
	polys.push_back(face2);
    
    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
    face2->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
    face2->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
  
	Polygon3 *face3 = new Polygon3();  
	polys.push_back(face3);

    face3->mVertices.push_back(Vector3(mMax.x, mMin.y ,mMin.z));
	face3->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
	face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
	face3->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));

	Polygon3 *face4 = new Polygon3();  
	polys.push_back(face4);

	face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
	face4->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));
	face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
	face4->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
    
	Polygon3 *face5 = new Polygon3();  
	polys.push_back(face5);

	face5->mVertices.push_back(Vector3(mMin.x, mMax.y, mMin.z));
    face5->mVertices.push_back(Vector3(mMax.x, mMax.y, mMin.z));
    face5->mVertices.push_back(Vector3(mMax.x, mMin.y, mMin.z));
	face5->mVertices.push_back(Vector3(mMin.x, mMin.y, mMin.z));

	Polygon3 *face6 = new Polygon3();  
	polys.push_back(face6);  
  
    face6->mVertices.push_back(Vector3(mMin.x, mMin.y, mMax.z));
    face6->mVertices.push_back(Vector3(mMax.x, mMin.y, mMax.z));
    face6->mVertices.push_back(Vector3(mMax.x, mMax.y, mMax.z));
    face6->mVertices.push_back(Vector3(mMin.x, mMax.y, mMax.z));

}


void AxisAlignedBox3::Split(const int axis, 
							const float value, 
							AxisAlignedBox3 &front, 
							AxisAlignedBox3 &back) const
{
	if ((value >= mMin[axis]) && (value <= mMax[axis]))
	{
		front.mMin = mMin; front.mMax = mMax;
		back.mMin = mMin; back.mMax = mMax;

		back.mMax[axis] = value;
		front.mMin[axis] = value;
	}
}

 
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-7;
  const float epsAdd = 1e-6;
  const float epsMul = 1.000005;
  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();

#if 0
  static HaltonSequence halton;
  
  halton.GenerateNext();
  
  return GetRandomPoint(Vector3(halton.GetNumber(1),
								halton.GetNumber(2),
								halton.GetNumber(3)));
  
#else
  return GetRandomPoint(Vector3(RandomValue(0.0f, 1.0f),
								RandomValue(0.0f, 1.0f),
								RandomValue(0.0f, 1.0f)));
#endif
}

Vector3
AxisAlignedBox3::GetRandomPoint(const Vector3 &r) const
{
  return mMin + Size()*r;
}

Vector3 AxisAlignedBox3::GetRandomSurfacePoint() const
{
	const int idx = Random(6);

	const Rectangle3 face = GetFace(idx);

	Vector3 point = Vector3(0,0,0);
	float sum = 0.0f;

	for (int i = 0; i < 4; ++ i) 
	{
		const float r = RandomValue(0, 1);
		sum += r;
		point += face.mVertices[i] * r;
	}

	point *= 1.0f / sum;

	//normal = face.GetNormal();

	return point;
}

Vector3
AxisAlignedBox3::GetUniformRandomSurfacePoint() const
{
	// get face based on its surface area
	float pdf[7];
	pdf[0] = 0.0f;
	int i;
	for (i=0; i < 6; i++) {
		pdf[i+1] = pdf[i] + GetFace(i).GetArea();
	}
	
	float x = RandomValue(0, pdf[6]);
	for (i=0; i < 6; i++) {
		if (x < pdf[i])
			break;
	}

	const int idx = i-1;
	
	const Rectangle3 face = GetFace(idx);
	
	Vector3 point = Vector3(0,0,0);
	float sum = 0.0f;

	for (i = 0; i < 4; ++ i) 
	{
		const float r = RandomValue(0, 1);
		sum += r;
		point += face.mVertices[i] * r;
	}

	point *= 1.0f / sum;
	//normal = face.GetNormal();

	return point;
}


 bool AxisAlignedBox3::Intersects(const Mesh &mesh) const
 {
	 VertexContainer::const_iterator vit, vit_end = mesh.mVertices.end();

	 for (vit = mesh.mVertices.begin(); vit != vit_end; ++ vit)
	 {
		if (IsInside(*vit))
			return true;
	 }
	 
	 return false;
 }


 bool AxisAlignedBox3::Intersects(const Triangle3 &tri) const
 {
	 if (IsInside(tri.mVertices[0]) ||
		 IsInside(tri.mVertices[1]) ||
		 IsInside(tri.mVertices[2]))
	 {
			return true;
	 }

	 return false;
 }


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

}