gdiam.h

任意给定三维空间的点集

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

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

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

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

    gdiam_real   get_diam() const
    {
        gdiam_real  sum, val;

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

        return  sqrt( sum );
    }

    // 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....
    gdiam_real   get_diam_proj( gdiam_point  dir ) const
    {
        gdiam_real  sum, coord;
        gdiam_real  prod, val;
        
        //printf( "get_diam_proj: " );
        sum = 0;
        for  ( int  ind = 0; ind < GDIAM_DIM; 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();
    }

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


// gdiam_bbox is the famous - arbitrary orieted bounding box
class  gdiam_bbox
{
private:
    gdiam_point_t  dir_1, dir_2, dir_3;
    
    gdiam_real  low_1, high_1;
    gdiam_real  low_2, high_2;
    gdiam_real  low_3, high_3;
    bool  f_init;

public:
    gdiam_bbox() { 
        f_init = false;
    }

    gdiam_real  volume() const {
        gdiam_real  len1, len2, len3;

        len1 = ( high_1 - low_1 );
        len2 = ( high_2 - low_2 );
        len3 = ( high_3 - low_3 );

        return  len1 * len2 * len3;
    }


    void  set_third_dim_longest() {
        gdiam_point_t  pnt_tmp;

        if  ( ( high_1 - low_1 ) > ( high_3 - low_3 ) ) {
            gdiam_exchange( high_1, high_3 );
            gdiam_exchange( low_1, low_3 );

            memcpy( pnt_tmp, dir_1, sizeof( dir_1 ) );
            //pnt_tmp = dir_1;
            memcpy( dir_1, dir_3, sizeof( dir_3) );
            //dir_1 = dir_3;
            memcpy( dir_3, pnt_tmp, sizeof( dir_3 ) );
            //dir_3 = pnt_tmp;
        }
        if  ( ( high_2 - low_2 ) > ( high_3 - low_3 ) ) {
            gdiam_exchange( high_2, high_3 );
            gdiam_exchange( low_2, low_3 );

            pnt_copy( pnt_tmp, dir_2 );
            pnt_copy( dir_2, dir_3 );
            pnt_copy( dir_3, pnt_tmp );
        }
    }

    gdiam_point  get_dir( int  ind ) {
        if  ( ind == 0 )
            return  dir_1;
        if  ( ind == 1 )
            return  dir_2;
        if  ( ind == 2 )
            return  dir_3;
     
        assert( false );

        return  NULL;
    }
    
    void  combine( gdiam_point  out,
                   double  a_coef, double  b_coef,
                   double  c_coef ) {
        for  ( int  ind = 0; ind < GDIAM_DIM; ind++ ) 
            out[ ind ] = a_coef * dir_1[ ind ] 
                + b_coef * dir_2[ ind ] 
                + c_coef * dir_3[ ind ];
    }



    void  dump_vertex( double  sel1, double  sel2, double  sel3 ) const {
        gdiam_real  coef1, coef2, coef3;

        //printf( "selection: (%g, %g, %g)\n", sel1, sel2, sel3 );
        coef1 = low_1 + sel1 * ( high_1 - low_1 );
        coef2 = low_2 + sel2 * ( high_2 - low_2 );
        coef3 = low_3 + sel3 * ( high_3 - low_3 );

        //printf( "coeficients: (%g, %g,  %g)\n", 
        //        coef1, coef2, coef3 );

        gdiam_point_t  pnt;
        pnt_zero( pnt );
        
        //printf( "starting...\n" );
        //pnt_dump( pnt );
        pnt_scale_and_add( pnt, coef1, dir_1 );
        //pnt_dump( pnt );
        pnt_scale_and_add( pnt, coef2, dir_2 );
        //pnt_dump( pnt );
        pnt_scale_and_add( pnt, coef3, dir_3 );
        //pnt_dump( pnt );

        pnt_dump( pnt );
    }

    void  dump() const {
        printf( "-----------------------------------------------\n" );
        dump_vertex( 0, 0, 0 );
        dump_vertex( 0, 0, 1 );
        dump_vertex( 0, 1, 0 );
        dump_vertex( 0, 1, 1 );
        dump_vertex( 1, 0, 0 );
        dump_vertex( 1, 0, 1 );
        dump_vertex( 1, 1, 0 );
        dump_vertex( 1, 1, 1 );
        printf( "volume: %g\n", volume() );
        printf( "Directions:\n" );
        pnt_dump( dir_1 );
        pnt_dump( dir_2 );
        pnt_dump( dir_3 );
        printf( "prods: %g, %g, %g\n",
                pnt_dot_prod( dir_1, dir_2 ),
                pnt_dot_prod( dir_1, dir_3 ),
                pnt_dot_prod( dir_2, dir_3 ) );

        printf( "--------------------------------------------------\n" );
        //printf( "range_1: %g... %g\n", low_1, high_1 );
        //printf( "range_2: %g... %g\n", low_2, high_2 );
        //printf( "range_3: %g... %g\n", low_3, high_3 );
        //printf( "prd: %g\n", pnt_dot_prod( dir_1, dir_2 ) );
        //printf( "prd: %g\n", pnt_dot_prod( dir_2, dir_3 ) );
        //printf( "prd: %g\n", pnt_dot_prod( dir_1, dir_3 ) );
    }


    void  init( const GBBox  & bb ) {
        pnt_init( dir_1, 1, 0, 0 );
        pnt_init( dir_2, 0, 1, 0 );
        pnt_init( dir_3, 0, 0, 1 );
    
        low_1 = bb.min_coord( 0 );
        high_1 = bb.max_coord( 0 );

        low_2 = bb.min_coord( 1 );
        high_2 = bb.max_coord( 1 );

        low_3 = bb.min_coord( 2 );
        high_3 = bb.max_coord( 2 );
        f_init = true;
    }

    void  init( const gdiam_point  _dir_1,
                const gdiam_point  _dir_2,
                const gdiam_point  _dir_3 ) {
        memset( this, 0, sizeof( gdiam_bbox ) );

        pnt_copy( dir_1, _dir_1 );
        pnt_copy( dir_2, _dir_2 );
        pnt_copy( dir_3, _dir_3 );
        pnt_normalize( dir_1 );
        pnt_normalize( dir_2 );
        pnt_normalize( dir_3 );

        if  ( ( ! (fabs( pnt_dot_prod( dir_1, dir_2 ) ) < 1e-6 ) )
              ||  ( ! (fabs( pnt_dot_prod( dir_1, dir_3 ) ) < 1e-6 ) )
              ||  ( ! (fabs( pnt_dot_prod( dir_2, dir_3 ) ) < 1e-6 ) ) ) {
            printf( "should be all close to zero: %g, %g, %g\n",
                    pnt_dot_prod( dir_1, dir_2 ),
                    pnt_dot_prod( dir_1, dir_3 ),
                    pnt_dot_prod( dir_2, dir_3 ) );
            pnt_dump( _dir_1 );
            pnt_dump( _dir_2 );
            pnt_dump( _dir_3 );
            fflush( stdout );
            fflush( stderr );
            assert( fabs( pnt_dot_prod( dir_1, dir_2 ) ) < 1e-6 );
            assert( fabs( pnt_dot_prod( dir_1, dir_3 ) ) < 1e-6 );
            assert( fabs( pnt_dot_prod( dir_2, dir_3 ) ) < 1e-6 );
        }

        // The following reduce the error by slightly improve the
        // orthoginality of the third vector. Doing to the second
        // vector is not efficient...
        pnt_scale_and_add( dir_3, -pnt_dot_prod( dir_3, dir_1 ),
                           dir_1 );
        pnt_scale_and_add( dir_3, -pnt_dot_prod( dir_3, dir_2 ),
                           dir_2 );
        pnt_normalize( dir_3 );
        //printf( "__should be all close to zero: %g, %g, %g\n",
        //         pnt_dot_prod( dir_1, dir_2 ),
        //        pnt_dot_prod( dir_1, dir_3 ),
        //        pnt_dot_prod( dir_2, dir_3 ) );        
    }

    void  bound( const gdiam_point  pnt ) {
        gdiam_real  prod_1, prod_2, prod_3;

        prod_1 = pnt_dot_prod( dir_1, pnt );
        prod_2 = pnt_dot_prod( dir_2, pnt );
        prod_3 = pnt_dot_prod( dir_3, pnt );

        if  ( ! f_init ) {
            f_init = true;
            low_1 = high_1 = prod_1;
            low_2 = high_2 = prod_2;
            low_3 = high_3 = prod_3;
            return;
        }
        if  ( prod_1 < low_1 )
            low_1 = prod_1;
        if  ( prod_2 < low_2 )
            low_2 = prod_2;
        if  ( prod_3 < low_3 )
            low_3 = prod_3;

        if  ( prod_1 > high_1 )
            high_1 = prod_1;
        if  ( prod_2 > high_2 )
            high_2 = prod_2;
        if  ( prod_3 > high_3 )
            high_3 = prod_3;
    }

    // compute the coordinates as if the bounding box is unit cube of
    // the grid
    void  get_normalized_coordinates( gdiam_point  in,
                                      gdiam_point  out ) {
        out[ 0 ] = ( pnt_dot_prod( dir_1, in ) - low_1 ) / ( high_1 - low_1 );
        out[ 1 ] = ( pnt_dot_prod( dir_2, in ) - low_2 ) / ( high_2 - low_2 );
        out[ 2 ] = ( pnt_dot_prod( dir_3, in ) - low_3 ) / ( high_3 - low_3 );
    }

};

GPointPair   gdiam_approx_diam( gdiam_point  * start, int  size,
                                    gdiam_real  eps );
GPointPair   gdiam_approx_diam( gdiam_point  * start, int  size,
                                gdiam_real  eps );
gdiam_real   gdiam_approx_diam( gdiam_real  * start, int  size,
                                gdiam_real  eps );
GPointPair   gdiam_approx_diam_pair( gdiam_real  * start, int  size,
                                     gdiam_real  eps );
GPointPair   gdiam_approx_diam_pair_UDM( gdiam_real  * start, int  size,
                                         gdiam_real  eps );
gdiam_bbox   gdiam_approx_const_mvbb( gdiam_point  * start, int  size,
                                      gdiam_real  eps, 
                                      GBBox  * p_ap_bbox );    
gdiam_point  * gdiam_convert( gdiam_real  * start, int  size );
gdiam_bbox   gdiam_approx_mvbb( gdiam_point  * start, int  size,
                                gdiam_real  eps ) ;
gdiam_bbox   gdiam_approx_mvbb_grid( gdiam_point  * start, int  size,
                                     int  grid_size );
gdiam_bbox   gdiam_approx_mvbb_grid_sample( gdiam_point  * start, int  size,
                                            int  grid_size, int  sample_size );
gdiam_bbox   gdiam_approx_mvbb_grid_sample( gdiam_real  * start, int  size,
                                            int  grid_size, int  sample_size );


void  gdiam_generate_orthonormal_base( gdiam_point  in,
                                       gdiam_point  out1,
                                       gdiam_point  out2 );

#else   /* __GDIAM__H */
#error  Header file gdiam.h included twice
#endif  /* __GDIAM__H */

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