Cubic_Bounding_Box.h

Go to the documentation of this file.

#ifndef _H_CUBIC_BOUNDING_BOX
#define _H_CUBIC_BOUNDING_BOX

#include "Bounding_Box.h"
#include "randomUtils.h"
#include "orderings.h"

#include <iterator>
#include <vector>
#include <sstream>
#include <iostream>
#include <fstream>

//#define _ASSERTS_ON

//#define _HP

#define _SHELL_CODE

#define _SHELL_CODE_MOD

//#define _DEBUG

class Cubic_Bounding_Box : public Bounding_Box< Cubic_Bounding_Box > {

friend class CubicLowCoorOrdering< 0, Cubic_Bounding_Box >;

friend class CubicLowCoorOrdering< 1, Cubic_Bounding_Box >;

friend class CubicLowCoorOrdering< 2, Cubic_Bounding_Box >;

friend class DistOrdering< Cubic_Bounding_Box >;

//the below should actually not be necessary
friend class Bisection_Tree;
friend class Cart;
friend class Basic_R_tree;

private:

        bool contains( const vector< double > * x ) const {

                if ( isEmpty() )
                        return false;

#ifndef _HP
                //dimension check
                if ( _mins.size() != x->size() ) {
                        cerr << "Warning: cannot check if a point with dimension not equal to bounding box dimension is contained within current bounding box!" << endl;
                        return false;
                }
#endif

                //per dimension check
                for ( unsigned int i=0; i<_mins.size(); i++ ) {
                        double current = (*x)[i];
                        if ( _maxs[i] < current )
                                return false;
                        if ( _mins[i] > current )
                                return false;
                }

                //no exit means containment:
                return true;
        }

        vector< double > axpy( const double a, const vector< double > * x, const vector< double > * y ) const {
                vector< double > ret;
                ret.resize( x->size() );
                for( unsigned int i=0; i<x->size(); i++ )
                        ret[ i ] = a * (*x)[ i ] + (*y)[ i ];
                return ret;
        }

protected:

        static vector<double> temp1;

        static vector<double> temp2;

        vector<double> _mins;

        vector<double> _maxs;

#ifdef _SHELL_CODE

        static double UtilDivideSafe( const double num, const double denom, const double defaultRes ) {
                return ( denom != 0.0 ? num / denom : defaultRes );
        }

        static bool BBOverlap( const double *bb1Min, const double *bb1Max, const double *bb2Min, const double *bb2Max,
                                                const int dim ) {
#ifndef _SHELL_CODE_MOD
                bool result = true;
#endif
                int d;

                if ( dim == 3 )
                {
#ifndef _SHELL_CODE_MOD
                        result = ( bb1Min[0] <= bb2Max[0] ) && ( bb1Max[0] >= bb2Min[0] ) && ( bb1Min[1] <= bb2Max[1] ) && ( bb1Max[1] >= bb2Min[1] ) && ( bb1Min[2] <= bb2Max[2] ) && ( bb1Max[2] >= bb2Min[2] );
#else
                        return ( bb1Min[0] <= bb2Max[0] ) && ( bb1Max[0] >= bb2Min[0] ) && ( bb1Min[1] <= bb2Max[1] ) && ( bb1Max[1] >= bb2Min[1] ) && ( bb1Min[2] <= bb2Max[2] ) && ( bb1Max[2] >= bb2Min[2] );
#endif
                }
                else
                {
#ifndef _SHELL_CODE_MOD
                        for ( d = 0; ( d < dim ) && result; d++ )
#else
                        for ( d = 0; ( d < dim ); d++ )
#endif
                        {
#ifndef _SHELL_CODE_MOD
                                result = ( bb1Min[d] <= bb2Max[d] ) && ( bb1Max[d] >= bb2Min[d] );
#else
                                if ( ( bb1Min[d] > bb2Max[d] ) || ( bb1Max[d] < bb2Min[d] ) )
                                        return false;
#endif
                        }
                }
#ifndef _SHELL_CODE_MOD
                return ( result );
#else
                return true;
#endif
        }

        static bool BBEdgeIntersect( const int dim, const double *point1, const double *point2, const double *minBB,
                                                                        const double *maxBB, const int dir, const double fixed_val )
        {
#ifndef _SHELL_CODE_MOD
                bool   intersect = false;
#endif
                const double cor = fixed_val;
                const double cor1 = point1[dir];
                const double cor2 = point2[dir];

                //assert( cor == maxBB[dir] );

#ifdef _SHELL_CODE_MOD
                if ( ( cor1 > cor ) && ( cor2 > cor ) )
                        return false;

                if ( ( cor1 < cor ) && (cor2 < cor ) )
                        return false;
#else
                if ( ( cor1 <= cor ) || ( cor2 <= cor ) )
                {
                        if ( ( cor1 >= cor ) || ( cor2 >= cor ) )
                        {
#endif
                                const double frac = UtilDivideSafe( cor - cor1, cor2 - cor1, 0.5 );
                                int          d;

#ifdef _ASSERTS_ON
                                assert( ( frac >= 0. ) && ( frac <= 1. ) );
#endif

                                /* assume that there is intersection unless proven otherwise */
#ifdef _SHELL_CODE_MOD
                                for ( d = 0; ( d < dim ); d++ )
#else
                                intersect = true;
                                for ( d = 0; ( d < dim ) && intersect; d++ )
#endif
                                {
                                        if ( d != dir )
                                        {
                                                const double c = ( 1 - frac ) * point1[d] + ( frac ) * point2[d];
#ifdef _SHELL_CODE_MOD
                                                if ( ( c < minBB[d] ) || ( c > maxBB[d] ) )
                                                        return false;
#else
                                                intersect = ( ( c >= minBB[d] ) && ( c <= maxBB[d] ) );
#endif
                                        }
#ifndef _SHELL_CODE_MOD
                                }
                        }
#endif
                }

#ifdef _SHELL_CODE_MOD
                return true;
#else
                return ( intersect );
#endif

        }

#else

        bool edgeLineIntersect( const vector< double > &start, const vector< double > &end, const vector< double > &direction, const int &dim_fixed, const double &fixed_coordinate ) const {

                if ( start[dim_fixed] > fixed_coordinate && end[dim_fixed] > fixed_coordinate )
                        return false;
                if ( start[dim_fixed] < fixed_coordinate && end[dim_fixed] < fixed_coordinate )
                        return false;

                const double scale = ( fixed_coordinate - start[dim_fixed] ) / direction[dim_fixed];

#ifdef _DEBUG
                cout << "Fixing coordinate " << dim_fixed << endl;
                cout << "Scalemin is " << scale << endl;
                cout << "Scalemax is " << scale << endl;
#endif

                //build yo, yt and do containment checking on the fly:
                for ( int d=0; d<dim_fixed; d++ ) {
                        const double yo = scale * direction[d] + start[d];
                        //if out range, we have no edge intersection on this dimension
                        if ( yo < _mins[d] || yo > _maxs[d] )
                                return false;
                }

                //for-loop is split to avoid checking for and avoiding (d==i)
                for ( unsigned int d=dim_fixed+1; d<start.size(); d++ ) {
                        const double yo = scale * direction[d] + start[d];
                        //if out range, we have no edge intersection on this dimension
                        if ( yo < _mins[d] || yo > _maxs[d] )
                                return false;
                }

                return true;

        }

        bool infiniteLineIntersect( const vector< double > const * start, const vector< double > const * end ) const {

                //old method is inefficient. It is more readable though...
                //vector< double > yo = axpy( ( _mins[i] - start[i] ) / direction[i], direction, start );
                //vector< double > yt = axpy( ( _maxs[i] - start[i] ) / direction[i], direction, start );
                //if( contains( yo ) || contains( yt ) )
                //      return true;

#ifdef _DEBUG
                cout << "MBR " << toString() << endl;
                cout << "Intersecting with line < ";
                for( unsigned int i=0; i<start->size(); i++ )
                        cout << (*start)[ i ] << " ";
                cout << "to < ";
                for( unsigned int i=0; i<start->size(); i++ )
                        cout << (*end)[ i ] << " ";
                cout << endl;
                cout << "Direction vector: < ";
#endif

                //build direction vector
                vector< double > direction;
                for( unsigned int i=0; i<start->size(); i++ ) {
                        direction.push_back( (*end)[ i ] - (*start)[ i ] );

#ifdef _DEBUG
                        cout << direction[ i ] << " ";
#endif

                        //fast contradiction detection
                        if ( direction[ i ] == 0 )
                                if ( (*start)[ i ] < _mins[ i ] || (*start)[ i ] > _maxs[ i ] )
                                        return false;
                }

#ifdef _DEBUG
                cout << endl;
                cout << "No contradiction found." << endl;
#endif

                //fix the i-th edge and check if the line intersects that edge
                for( unsigned int i=0; i<start->size(); i++ ) {

                        if ( edgeLineIntersect( *start, *end, direction, i, _mins[ i ] ) )
                                return true;
                        if ( edgeLineIntersect( *start, *end, direction, i, _maxs[ i ] ) )
                                return true;

                }

                return false;
        }
#endif

public:

        static CubicLowCoorOrdering< 0, Cubic_Bounding_Box > LCX_ORDERING;
        
        static CubicLowCoorOrdering< 1, Cubic_Bounding_Box > LCY_ORDERING;
        
        static CubicLowCoorOrdering< 2, Cubic_Bounding_Box > LCZ_ORDERING;
        
        static DistOrdering< Cubic_Bounding_Box > DIST_ORDERING;

        Cubic_Bounding_Box() {
                clear();
        }

        virtual void writeToFile( ofstream &ofs ) const {

                ofs << "-1 ";
                ofs << (_mins[0]) << " " << (_mins[1]) << " " << (_mins[2]) << " ";
                ofs << (_maxs[0]) << " " << (_maxs[1]) << " " << (_maxs[2]) << endl;

        }
        
        virtual bool lineIntersect( const vector< double > &start, const vector< double > &end ) const {

#ifdef _SHELL_CODE
                        bool intersect = false;

                        const unsigned int dim = start.size();

                        unsigned int d;

                        /* Any point inside the bounding box? */
                        if ( !intersect )
                        {
                                //intersect = BBOverlap( point1, point1, minBB, maxBB, dim );
                                intersect = BBOverlap( &(start)[ 0 ], &(start)[ 0 ], &_mins[ 0 ], &_maxs[ 0 ], dim );
                        }

                        if ( !intersect )
                        {
                                intersect = BBOverlap( &(end)[ 0 ], &(end)[ 0 ], &_mins[ 0 ], &_maxs[ 0 ], dim );
                        }

                        /* test left edge of bounding box */
                        for ( d = 0; ( d < dim ) && !intersect; d++ )
                        {
                                //const double corStore = _mins[d];

                                //_mins[d] = _maxs[d];
                                intersect = BBEdgeIntersect( dim, &(start)[ 0 ], &(end)[ 0 ], &_mins[ 0 ], &_maxs[ 0 ], d, _maxs[ d ] );
                                //minBB[d] = corStore;
                        }

                        /* test right edge of bounding box */
                        for ( d = 0; ( d < dim ) && !intersect; d++ )
                        {
                                //const double corStore = _maxs[d];

                                //_maxs[d] = _mins[d];
                                intersect = BBEdgeIntersect( dim, &(start)[ 0 ], &(end)[ 0 ], &_mins[ 0 ], &_maxs[ 0 ], d, _mins[ d ] );
                                //_maxs[d] = corStore;
                        }

                return ( intersect );
#else

                if ( contains( Point( start ) ) || contains( Point( end ) ) )
                        return true;

                for( unsigned int i=0; i<start->size(); i++ ) {

                        //if no interval overlap, return false
                        if ( (*start)[ i ] < (*end)[ i ] ) {
                                if ( (*end)[ i ] < _mins[i] || (*start)[ i ] > _maxs[i] )
                                        return false;
                        } else {
                                if ( (*start)[ i ] < _mins[i] || (*end)[ i ] > _maxs[i] )
                                        return false;
                        }

                }

                //check if the line truly intersects and does not have a 'fake' overlap
                return infiniteLineIntersect( start, end );
#endif

        }

        Cubic_Bounding_Box( vector<double> &mins, vector<double> &maxs ) {

#ifndef _HP
                //dimension check
                if ( mins.size() != maxs.size() ) {
                        cerr << "Dimension mis-match in Cubic_Bounding_Box constructor!" << endl;
                        exit(1);
                }
#endif

                //set variables
                _mins = mins; _maxs = maxs;
        }

        Cubic_Bounding_Box( const vector<double> &mins, const vector<double> &maxs ) {

#ifndef _HP
                //dimension check
                if ( mins.size() != maxs.size() ) {
                        cerr << "Dimension mis-match in Cubic_Bounding_Box constructor!" << endl;
                        exit(1);
                }
#endif

                //set variables
                _mins = mins; _maxs = maxs;
        }

        ~Cubic_Bounding_Box() {}


        Cubic_Bounding_Box( const Cubic_Bounding_Box *toCopy ) {
                become( toCopy );
        }

        virtual void clear() {
                _mins.resize( 1 ); _maxs.resize( 1 );
                _mins[ 0 ]=0; _maxs[ 0 ]=0;
        }

        virtual void become( const Cubic_Bounding_Box *toCopy ) {
#ifdef _DEBUG
                cout << "Becoming: "; fflush( stdout );
                cout << toCopy->toString() << endl; fflush( stdout );
                cout << "Dimension check... " << endl; fflush( stdout );
#endif

#ifndef _HP
                //dimension check
                if ( toCopy->_mins.size() != toCopy->_maxs.size() ) {
                        cerr << "Dimension mismatch in Cubic_Bounding_Box constructor!" << endl;
                        exit(1);
                }
#endif

#ifdef _DEBUG
                cout << "Setting variables..." << endl; fflush( stdout );
#endif

                //set variables
                _mins = toCopy->_mins; _maxs = toCopy->_maxs;

#ifdef _DEBUG
                cout << "Become successful" << endl; fflush( stdout );
#endif

        }

        double getMinCoordinateOnDimension( unsigned int i ) const { return _mins[ i ]; }

        double getMaxCoordinateOnDimension( unsigned int i ) const { return _maxs[ i ]; }

        virtual bool isEmpty() const {
                return ( _mins.size()==1 && _mins[0]==0 && _maxs[0]==0 );
        }

        virtual void getCenterCoordinate( vector< double > * into ) const {
                //vector< double > ret;

                into->resize( _mins.size() );
                for( unsigned int i=0; i<_mins.size(); i++ )
                        (*into)[ i ] = ( _mins[i] + _maxs[i] ) / 2;
                        //ret.push_back( ( _mins[i] + _maxs[i] ) / 2 );

                //return &ret;
        }

        virtual bool intersects( const Cubic_Bounding_Box &other ) const {
                if ( isEmpty() || other.isEmpty() )
                        return false;

#ifndef _SHELL_CODE

#ifndef _HP
                //dimension check
                if ( _mins.size() != other._mins.size() ) {
                        cerr << "Warning: cannot compare two bounding boxes of different dimensions!" << endl;
                        return false;
                }
#endif

                //per dimension check
                for ( int i=0; i<_mins.size(); i++ ) {
                        if ( _maxs[i] < other._mins[i] )
                                return false;
                        if ( other._maxs[i] < _mins[i] )
                                return false;
                }

                //did not exit during check; we intersect:
                return true;
#else
                return BBOverlap( &_mins[0], &_maxs[0], &(other._mins[0]), &(other._maxs[0]), _mins.size() );
#endif
        }


        virtual bool contains( const Point &x ) const {

                if ( isEmpty() )
                        return false;

                //dimension check
                if ( _mins.size() != x.getDimension() ) {
                        cerr << "Warning: cannot check if a point with dimension not equal to bounding box dimension is contained within current bounding box!" << endl;
                        return false;
                }

                //per dimension check
                for ( unsigned int i=0; i<_mins.size(); i++ ) {
                        double current = x.getCoordinate(i);
                        if ( _maxs[i] < current )
                                return false;
                        if ( _mins[i] > current )
                                return false;
                }

                //no exit means containment:
                return true;
        }

        virtual bool contains( const Cubic_Bounding_Box &other ) const {
                return ( contains( &( other._mins ) ) && contains( &( other._maxs ) ) );
        }

        virtual void unite( const Cubic_Bounding_Box* other ) {
                if ( other->isEmpty() )
                        return;

                if ( isEmpty() )
                        become( other );
                        
#ifndef _HP
                //dimension check
                if ( _mins.size() != other->_mins.size() ) {
                        cerr << "Error: Dimension mismatch in Cubic_Bounding_Box::unite() !" << endl;
                        exit(1);
                }
#endif

                //per dimension builder
                for ( unsigned int i=0; i<_mins.size(); i++ ) {
                        if ( _mins[ i ] > other->_mins[ i ] )
                                _mins[ i ] = other->_mins[ i ];
                        if ( _maxs[ i ] < other->_maxs[ i ] )
                                _maxs[ i ] = other->_maxs[ i ];
                }
        }

        virtual Cubic_Bounding_Box unionWith( const Cubic_Bounding_Box &other ) const {
#ifdef _DEBUG
                cout << "Unite called" << endl; fflush( stdout );
#endif

                if ( isEmpty() && other.isEmpty() )
                        return Cubic_Bounding_Box();

                if ( isEmpty() )
                        return Cubic_Bounding_Box( other._mins, other._maxs );

                if ( other.isEmpty() )
                        return Cubic_Bounding_Box( _mins, _maxs );

#ifndef _HP
                //dimension check
                if ( _mins.size() != other._mins.size() ) {
                        cerr << "Error: Dimension mismatch in Cubic_Bounding_Box::unite() !" << endl;
                        exit(1);
                }
#endif
                
                temp1.resize( _mins.size() ); temp2.resize( _mins.size() );

                //per dimension builder
                for ( unsigned int i=0; i<_mins.size(); i++ ) {
                        if ( _mins[i] < other._mins[i] )
                                temp1[i] = _mins[i];
                        else
                                temp1[i] = other._mins[i];
                        if ( _maxs[i] > other._maxs[i] )
                                temp2[i] = _maxs[i];
                        else
                                temp2[i] = other._maxs[i];
                }

                return Cubic_Bounding_Box( temp1, temp2 );
        }

        virtual Cubic_Bounding_Box intersect( const Cubic_Bounding_Box &other ) const {

#ifdef _DEBUG
                cout << "THIS: " << endl << toString();
                cout << "THAT: " << endl << other.toString();
#endif

                if ( isEmpty() || other.isEmpty() ) {
                        return Cubic_Bounding_Box();
                }

#ifndef _HP
                //dimension check
                if ( _mins.size() != other._mins.size() ) {
                        cerr << "Error: Dimension mismatch in Cubic_Bounding_Box::unite() !" << endl;
                        exit(1);
                }
#endif
                
                temp1.resize( _mins.size() ); temp2.resize( _mins.size() );

                //per dimension builder
                for ( unsigned int i=0; i<_mins.size(); i++ ) {
                        if ( ( other._maxs[i] < _mins[i] ) || ( other._mins[i] > _maxs[i] ) ) {
                                temp1[i] = 1; temp2[i] = 0;
                        } else {

                                if ( _mins[i] > other._mins[i] )
                                        temp1[i] = _mins[i];
                                else
                                        temp1[i] = other._mins[i];

                                if ( _maxs[i] < other._maxs[i] )
                                        temp2[i] = _maxs[i];
                                else
                                        temp2[i] = other._maxs[i];
                        }
                }

#ifdef _DEBUG
                Cubic_Bounding_Box ret( Cubic_Bounding_Box( temp1, temp2 ) );
                cout << "RET : " << endl << ret.toString();
                return ret;
#else
                return Cubic_Bounding_Box( temp1, temp2 );
#endif

        }

        static Ordering< Cubic_Bounding_Box > * getDefaultOrdering();
        
        static vector< Ordering< Cubic_Bounding_Box > * > getOrderings();

        virtual string toString() const {

                ostringstream os;
                os << _mins.size() << ": ";
                for ( unsigned int i=0; i<_mins.size(); i++ )
                        os << "< " << _mins[i] << " , " << _maxs[i] << " >" << endl;
                return os.str();
        }

        virtual unsigned int getDimension() const {
                return _mins.size();
        }

        virtual bool isDefining( Cubic_Bounding_Box * other ) const {
                bool ret = _mins[ 0 ] < other->_mins[ 0 ];
                ret = ret && ( _maxs[ 0 ] > other->_maxs[ 0 ] );
                for ( unsigned int i=1; ( i<_mins.size() ) && ret; i++ ) {
                        ret = ret && ( _mins[ i ] < other->_mins[ i ] );
                        ret = ret && ( _maxs[ i ] > other->_maxs[ i ] );
                }
                return !ret;
        }

};



/***********************************************
 *                   PUBLIC                    *
 ***********************************************/

ostream& operator<<( ostream &os, const Cubic_Bounding_Box box );

Cubic_Bounding_Box operator&&( const Cubic_Bounding_Box left,  const Cubic_Bounding_Box right );

Cubic_Bounding_Box operator||( const Cubic_Bounding_Box left,  const Cubic_Bounding_Box right );

#endif

---