source: trunk/VUT/GtpVisibilityPreprocessor/src/AxisAlignedBox3.cpp @ 448

// GOLEM library
#include <assert.h>
#include <iostream>
using namespace std;
#include "AxisAlignedBox3.h"
#include "Ray.h"
#include "Polygon3.h"

#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 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
{
  register float minx, maxx;
  ray.ComputeInvertedDir();
  
  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.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) < 0.001) {
    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) < 0.001) {
    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

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.0 * (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.0 + 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.0 - 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.0)) - atan22(b*c/sqrt(b*b + c*c + 1.0))
      - atan22(d*a/sqrt(a*a + d*d + 1.0)) + atan22(a*c/sqrt(a*a + c*c + 1.0));
  }

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

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 intersectioins 
        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 = 1e6;
        minT = 0;
        
        // test with  bounding box
        if (!GetMinMaxT(ray, &minT, &maxT))
                return false;

        if (minT < 0) // start ray from origin
                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;
}