source: GTP/trunk/App/Demos/Vis/CHC_revisited/src/AxisAlignedBox3.cpp @ 2782

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


using namespace std;


namespace CHCDemoEngine
{

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


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


AxisAlignedBox3::AxisAlignedBox3() 
{ }


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



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


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


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


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


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

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

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

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


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

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


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


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

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

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


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


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

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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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


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

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

        return true;
}

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

        return false;
}

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

        return false; // has non-zero volume
}


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

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

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

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

        minDistance = sumMin;
        maxDistance = sumMax;
}


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


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


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


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


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


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


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


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


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


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

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


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


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


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

        assert(dir >= 0.f);

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

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


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


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


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


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


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


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


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

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


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

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

                        const float di = *ei - *si;

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

                        const float di = *ei - *si;

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

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

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

        //*time = fst;
        return true;
}


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

        return v;
}


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

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

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


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

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

        /*int index;

        if (vecPlaneNormal.x <= 0.0f)
                index = (vecPlaneNormal.y <= 0.0f) ? 0 : 3;
        else
                index = (vecPlaneNormal.y <= 0.0f) ? 1 : 2;

        return (vecPlaneNormal.z <= 0.0f) ? index : 7 - index;*/
}


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


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

}