glwindow_geometry.h

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * geometry.h -
 *     Defines basic geometry prmitives.
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/

#ifndef  __GEOMETRY__H
#define  __GEOMETRY__H

//#ifndef  ___XGENERIC__H
//#include  "_generic.h"
//#endif

#define  X  0
#define  Y  1
#define  Z  2
#define  W  3

typedef  double greal;
#define  GDIM 3
typedef  greal    gpoint3_t[ 3 ];
typedef  greal  * gpoint3;
typedef  greal    gpoint4_t[ 4 ];
typedef  greal  * gpoint4;
typedef  const greal  * gpoint3_cnt;

class  GPoint3
{
public:
    gpoint3_t  pnt;

    GPoint3() {}
    GPoint3( greal  a, greal  b, greal  c ) {
        pnt[ 0 ] = a;
        pnt[ 1 ] = b;
        pnt[ 2 ] = c;
    }
};


class  GPoint4
{
public:
    gpoint4_t  pnt;

    GPoint4( const GPoint3  & p, double  _w ) {
        pnt[ 0 ] = p.pnt[ 0 ];
        pnt[ 1 ] = p.pnt[ 1 ];
        pnt[ 2 ] = p.pnt[ 2 ];
        pnt[ 3 ] = _w;
    }

    GPoint4( greal  a, greal  b, greal  c, greal  d ) {
        pnt[ 0 ] = a;
        pnt[ 1 ] = b;
        pnt[ 2 ] = c;
        pnt[ 3 ] = d;
    }
    GPoint4() {}
    GPoint4( const gpoint4_t  & p ) {
        pnt[ 0 ] = p[ 0 ];
        pnt[ 1 ] = p[ 1 ];
        pnt[ 2 ] = p[ 2 ];
        pnt[ 3 ] = p[ 3 ];
    }

    greal  operator[]( int pos ) const {
        return  pnt[ pos ];
    }
    const GPoint4  & operator+=( const GPoint4  & src ) {
        pnt[ 0 ] += src.pnt[ 0 ];
        pnt[ 1 ] += src.pnt[ 1 ];
        pnt[ 2 ] += src.pnt[ 2 ];
        pnt[ 3 ] += src.pnt[ 3 ];        

        return  *this;
    }
    const GPoint4  & operator-=( const GPoint4  & src ) {
        pnt[ 0 ] -= src.pnt[ 0 ];
        pnt[ 1 ] -= src.pnt[ 1 ];
        pnt[ 2 ] -= src.pnt[ 2 ];
        pnt[ 3 ] -= src.pnt[ 3 ];        

        return  *this;
    }
    const GPoint4  & operator*=( const greal  num ) {
        pnt[ 0 ] *= num;
        pnt[ 1 ] *= num;
        pnt[ 2 ] *= num;
        pnt[ 3 ] *= num;        

        return  *this;
    }
    const GPoint4  & operator/=( const greal  num ) {
        pnt[ 0 ] /= num;
        pnt[ 1 ] /= num;
        pnt[ 2 ] /= num;
        pnt[ 3 ] /= num;        

        return  *this;
    }

};

inline GPoint4  operator+( const GPoint4  & a,
                          const GPoint4  & b ) {
    GPoint4  tmp = a;
    
    tmp += b;

    return  tmp;
}

inline greal  operator*( const GPoint4  & a,
                         const GPoint4  & b ) {
    greal  sum;

    sum = a.pnt[ 0 ] * b.pnt[ 0 ] 
        + a.pnt[ 1 ] * b.pnt[ 1 ] 
        + a.pnt[ 2 ] * b.pnt[ 2 ] 
        + a.pnt[ 3 ] * b.pnt[ 3 ];

    return  sum;
}


inline GPoint4  operator-( const GPoint4  & a ) {
    GPoint4  tmp( 0, 0, 0, 0 );
    
    tmp -= a;

    return  tmp;
}


inline GPoint4  operator-( const GPoint4  & a, const GPoint4  & b ) {
    GPoint4  tmp( a );
    
    tmp -= b;

    return  tmp;
}


inline GPoint4  operator*( const GPoint4  & a, greal  num ) {
    GPoint4  tmp = a;
    
    tmp *= num;

    return  tmp;
}

inline GPoint4  operator/( const GPoint4  & a, greal  num ) {
    GPoint4  tmp = a;
    
    tmp /= num;

    return  tmp;
}

inline greal   norm2( const GPoint4  & a )
{
    return  a * a;
}


// GBBox -  Axis parallel bounding box.
class   GBBox {
private:
    greal  min_coords[ GDIM ];
    greal  max_coords[ GDIM ];

public:
    void  init( const GBBox   & a,
                const GBBox   & b ) {
        for  (int  ind = 0; ind < GDIM; ind++ ) { 
            min_coords[ ind ] = min( a.min_coords[ ind ],
                                     b.min_coords[ ind ] );
            max_coords[ ind ] = max( a.max_coords[ ind ],
                                     b.max_coords[ ind ] );
        }
    }
    void  dump() const {
        greal  prod, diff;

        prod = 1.0;
        printf( "__________________________________________\n" );
        for  (int  ind = 0; ind < GDIM; ind++ ) { 
            printf( "%d: [%g...%g]\n", ind, min_coords[ ind ],
                    max_coords[ ind ] );
            diff = max_coords[ ind ] - min_coords[ ind ];
            prod *= diff;
        }
        printf( "volume = %g\n", prod );
        printf( "\\__________________________________________\n" );
    }

    void  center( gpoint3  out ) const {
        for  (int  ind = 0; ind < GDIM; ind++ ) {
            out[ ind ] = ( min_coords[ ind ] + max_coords[ ind ] ) / 2.0;
        }
    }

    void  center( GPoint3  & out ) const {
        for  (int  ind = 0; ind < GDIM; ind++ ) {
            out.pnt[ ind ] = ( min_coords[ ind ] + max_coords[ ind ] ) / 2.0;
        }
    }

    greal  volume() const {
        greal  prod, val;

        prod = 1;        
        for  ( int  ind = 0; ind < GDIM; ind++ ) {
            val = length_dim( ind );
            prod *= val;
        }

        return  prod;
    }
    void  init() { 
        for  (int  ind = 0; ind < GDIM; ind++ ) {
            min_coords[ ind ] = 1e20;
            max_coords[ ind ] = -1e20;
        }
    }
    /*
    void  dump() const {
        printf( "---(" );
        for  (int  ind = 0; ind < GDIM; ind++ ) {
            printf( "[%g, %g] ",
                    min_coords[ ind ],
                    max_coords[ ind ] );
        }
        printf( ")\n" );
        }*/

    const greal  & min_coord( int  coord ) const {
        return  min_coords[ coord ];
    }
    const greal  & max_coord( int  coord ) const {
        return  max_coords[ coord ];
    }
    
    void  bound( const double  x, 
                 const double  y, 
                 const double  z ) {
        gpoint3_t  pnt;

        pnt[ 0 ] = x;
        pnt[ 1 ] = y;
        pnt[ 2 ] = z;

        bound( pnt );        
    }


    void  bound( const gpoint3   pnt ) {
        //cout << "bounding: " << pnt << "\n";
        for  ( int  ind = 0; ind < GDIM; ind++ ) {
            if  ( pnt[ ind ] < min_coords[ ind ] )
                min_coords[ ind ] = pnt[ ind ];
            if  ( pnt[ ind ] > max_coords[ ind ] )
                max_coords[ ind ] = pnt[ ind ];
        }
    }

    greal  length_dim( int  dim )  const {
        return  max_coords[ dim ] - min_coords[ dim ];
    }
    int  getLongestDim() const {
        int  dim = 0;
        greal  len = length_dim( 0 );

        for  ( int  ind = 1; ind < GDIM; ind++ ) 
            if  ( length_dim( ind ) > len ) {
                len = length_dim( ind );
                dim = ind;
            }

        return  dim;
    }
    greal  getLongestEdge() const {
        return  length_dim( getLongestDim() );
    }

    greal   get_diam() const
    {
        greal  sum, val;

        sum = 0;        
        for  ( int  ind = 0; ind < GDIM; ind++ ) {
            val = length_dim( ind );
            sum += val * val;
        }

        return  sqrt( sum );
    }

    greal  get_radius() const 
    {
        return  get_diam() / 2.0;
    }


    // in the following we assume that the length of dir is ONE!!!
    // Note that the following is an overestaime - the diameter of
    // projection of a cube, is bounded by the length of projections
    // of its edges....
    greal   get_diam_proj( gpoint3  dir ) const
    {
        greal  sum, coord;
        greal  prod, val;
        
        //printf( "get_diam_proj: " );
        sum = 0;
        for  ( int  ind = 0; ind < GDIM; ind++ ) {
            coord = length_dim( ind );
            //printf( "coord[%d]: %g\n",ind, coord );
            prod = coord * dir [ ind ];
            val = coord * coord - prod * prod;
            assert( val >= 0 );
            sum += sqrt( val );
        }
        //printf( "sum: %g, %g\n", sum, get_diam() );
        
        // sum = squard diameter of the bounding box
        // prod = length of projection of the diameter of cube on the
        //      direction.
        
        return  get_diam();
    }

    greal  get_min_coord( int  dim ) const {
        return  min_coords[ dim ];
    }
    greal  get_max_coord( int  dim ) const {
        return  max_coords[ dim ];
    }
};

inline void  pnt_scale_and_add( gpoint3  dest,
                         greal  coef,
                         gpoint3_cnt   vec ) {
    dest[ 0 ] += coef * vec[ 0 ];
    dest[ 1 ] += coef * vec[ 1 ];
    dest[ 2 ] += coef * vec[ 2 ];
}

inline void  pnt_scale( gpoint3  dest,
                        greal  coef ) {
    dest[ 0 ] = coef * dest[ 0 ];
    dest[ 1 ] = coef * dest[ 1 ];
    dest[ 2 ] = coef * dest[ 2 ];
}

#else   /* __GEOMETRY__H */
#error  Header file geometry.h included twice
#endif  /* __GEOMETRY__H */

/* geometry.h - End of File ------------------------------------------*/