algorithm.c

任意给定三维空间的点集

//------------------------------------------------------------------
// DiameterDirections -
//     Approx diameter by projecting the points into alot of
//   directions
//------------------------------------------------------------------
DiameterAlgorithm  * DiameterDirections::CreateAndCompute( ArrPoint3d  & arr,
                         double  _eps )
{
    DiameterDirections   * pDiam = new DiameterDirections;

    pDiam->eps = _eps;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}

char   * DiameterDirections::getID()
{
    static char  * name = "Diametr from directions";

    return  name;
}


void   DiameterDirections::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    int  alpha_limit, beta_limit;
    real  delta, alpha, beta;

    diam.init( arr[ 0 ] );
    tm.start();
    
    delta = sqrt( eps );
    alpha_limit = (int)(2.0 * 3.1415926535 / delta);
    beta_limit = (int)(3.1415926535 / delta);
    
    alpha = 0.0;
    for  ( int  ind_alpha = 0; ind_alpha < alpha_limit; 
           ind_alpha++, alpha += delta ) {
        beta = -3.1415926535 / 2.0;
        for  ( int  ind_beta = 0; ind_beta < beta_limit; 
               ind_beta++, beta += delta ) {
            Point3d  pnt;

            pnt.setCoord( 0, cos( beta ) * sin( alpha ) );
            pnt.setCoord( 1, cos( beta ) * cos( alpha ) );
            pnt.setCoord( 2, sin( beta ) );

            PointPair  pp;
            computeExtremePairDir( arr, pnt, pp );
            diam.update_diam( pp );
        }
    }
    
    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// DiameterGrid -
//     We filter out points by using a grid,
//------------------------------------------------------------------

class  GridCleaner
{
public:
    ArrPoint3d   arr, out_arr;
    Point3d  left_bottom;
    int  side_size, overall_size;
    real  x_delta, y_delta, z_delta;
    Point3d  ** shafts_bottom, ** shafts_top;

    ~GridCleaner() {
        arr.term();
        out_arr.term();
        free( shafts_bottom );
        free( shafts_top );
        shafts_bottom = shafts_top = NULL;
    }

    void  init( const ArrPoint3d  & _arr, 
                double   eps );
    void  map_xy( Point3d  pnt, int  & x, int  & y ) {
        x = (int)( ( pnt[ 0 ] - left_bottom[ 0 ] ) / x_delta );
        y = (int)( ( pnt[ 1 ] - left_bottom[ 1 ] ) / y_delta );
    }
    int  cell( Point3d  pnt ) {
        int  x, y;

        map_xy( pnt, x, y );
        //printf( "x: %d, y: %d  size_side: %d\n", x, y,
        //        side_size );
        if  ( y == side_size )
            y--;
        return  y * side_size + x;
    }
    void   register_point( Point3d  & pnt ) {
        int  ind = cell( pnt );

        if  ( ( ind < 0 )  ||  ( ind >= overall_size ) ) {
            printf( "ind: %d, overall size: %d\n",
                    ind, overall_size );
            assert( (0 <= ind )  &&  ( ind < overall_size ) );
        }

        if   ( ( shafts_bottom[ ind ] == NULL )
               ||  ( (*(shafts_bottom[ ind ]))[ 2 ] < pnt[ 2 ] ) )
            shafts_bottom[ ind ] = &pnt;

        if   ( ( shafts_top[ ind ] == NULL )
               ||  ( (*(shafts_top[ ind ]))[ 2 ] > pnt[ 2 ] ) )
            shafts_top[ ind ] = &pnt;
    }
};


void  GridCleaner::init( const ArrPoint3d  & _arr, 
                         double   eps )
{ 
    arr.init( _arr.size(), 43543 );

    for  ( int  ind = 0; ind <  _arr.size(); ind++ ) 
        arr.pushx( _arr[ ind ] );

    BBox3d  bb;

    pnt_array_bound( arr, bb );
    
    left_bottom.setCoord( 0, bb.min_coord( 0 ) );
    left_bottom.setCoord( 1, bb.min_coord( 1 ) );
    left_bottom.setCoord( 2, bb.min_coord( 2 ) );

    side_size = (int)( 1.0 + sqrt( 6.0 ) / eps );
    x_delta = (1.00000001 * bb.length_dim( 0 ) ) / (double)side_size;
    y_delta = (1.00000001 * bb.length_dim( 1 ) ) / (double)side_size;
    z_delta = (1.00000001 * bb.length_dim( 2 ) ) / (double)side_size;
 
    overall_size = side_size * side_size;

    shafts_bottom = (Point3d **)malloc( sizeof( Point3d * ) * overall_size );
    assert( shafts_bottom != NULL );
    memset( shafts_bottom, 0, sizeof( Point3d * ) * overall_size  );

    shafts_top = (Point3d **)malloc( sizeof( Point3d * ) * overall_size );
    assert( shafts_top != NULL );
    memset( shafts_top, 0, sizeof( Point3d * ) * overall_size  );
    
    for  ( int  ind = 0; ind <  arr.size(); ind++ ) 
        register_point( arr[ ind ] );

    out_arr.init( 100, 45344354 );
    for  ( int  ind = 0; ind < overall_size; ind++ ) { 
        if  ( ( shafts_bottom[ ind ] != NULL )
              &&  ( shafts_top[ ind ] != NULL ) ) {
            out_arr.pushx( *( shafts_bottom[ ind ] ) );
            if  ( shafts_bottom[ ind ] != shafts_top[ ind ] )
                out_arr.pushx( *( shafts_top[ ind ] ) );
        } else {
            if  ( shafts_bottom[ ind ] != NULL )
                out_arr.pushx( *( shafts_bottom[ ind ] ) );
            if  ( shafts_top[ ind ] != NULL )
                out_arr.pushx( *( shafts_top[ ind ] ) );
            
        }
    }
}


DiameterAlgorithm  * DiameterGrid::CreateAndCompute( ArrPoint3d  & arr,
                                                     double  _eps )
{
    DiameterGrid  * pDiam = new DiameterGrid;

    pDiam->eps = _eps;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * DiameterGrid::getID() 
{
    static char  * name = "Diameter grid";

    return  name;
}


void   DiameterGrid::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    GridCleaner  grid;
    
    diam.init( arr[ 0 ] );

    tm.start();

    grid.init( arr, eps / 2.0 );
    DiameterAlgorithm  * pDiamHeap;

    printf( "grid cleaning\n"
            "\tin_size : %d\n"
            "\tout_size: %d\n", 
            arr.size(), grid.out_arr.size() );
    pDiamHeap = QuadDiameterByLevel::CreateAndCompute( grid.out_arr, 
                              
                         eps / 2.0 );
    
    diam = pDiamHeap->diam;

    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// DiameterGridQuadDir -
//     We filter out points by using a grid, and then use fiar 
// split with direcitons
//------------------------------------------------------------------


DiameterAlgorithm  * DiameterGridQuadDir::CreateAndCompute( ArrPoint3d  & arr,
                                                            double  _eps )
{
    DiameterGridQuadDir  * pDiam = new DiameterGridQuadDir;

    pDiam->eps = _eps;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * DiameterGridQuadDir::getID() 
{
    static char  * name = "Diameter grid->FSTree dir";

    return  name;
}


void   DiameterGridQuadDir::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    GridCleaner  grid;
    
    diam.init( arr[ 0 ] );

    tm.start();

    grid.init( arr, eps / 2.0 );
    DiameterAlgorithm  * pDiamHeap;

    //printf( "\tgrid cleanring in_size: %d\nout_size: %d\n", 
    //        arr.size(), grid.out_arr.size() );
    pDiamHeap = QuadDiameterByDir::CreateAndCompute( grid.out_arr, eps/2 );
    
    diam = pDiamHeap->diam;

    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// ModifiedDiameterChan - 
//     Describes the algorithm due to Timothy Chan Socg2000
//     On the 2d level we use cleaup & fiar split tree to get
//   better resuls.
//------------------------------------------------------------------

DiameterAlgorithm  * DiameterModChan::CreateAndCompute(
                                              ArrPoint3d  & arr,
                                              double  _eps )
{
    DiameterModChan  * pDiam = new DiameterModChan;

    pDiam->eps = _eps;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * DiameterModChan::getID()
{
    static char  * name = "Chan's mod algorithm";

    return  name;
}


void   DiameterModChan::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    GridCleaner  grid;
    double  delta, alpha, alpha_limit;

    diam.init( arr[ 0 ] );

    tm.start();

    //for  ( int  ind = 0; ind < arr.size(); ind++ ) 
    //    printf( "___id: %d\n", arr[ ind ].getID() );
        
    grid.init( arr, eps / 3.0 );


    //for  ( int  ind = 0; ind < grid.out_arr.size(); ind++ ) 
    //    printf( "x___id: %d\n", grid.out_arr[ ind ].getID() );

    //printf( "in_size: %d\nout_size: %d\n", 
    //        arr.size(), grid.out_arr.size() );
    
    delta = sqrt( eps / 3.0 );
    alpha_limit = (int)(3.1415926535 / delta);

    alpha = 0.0;
    for  ( int  ind_alpha = 0; ind_alpha < alpha_limit; 
           ind_alpha++, alpha += delta ) {
        ProjPoints  prj;
        Point3d  dir;
        ArrPoint2dPtr  out_arr;

        dir[ 0 ] = sin( alpha );
        dir[ 1 ] = cos( alpha );
        dir[ 3 ] = 0;
        
        prj.init( grid.out_arr, dir );

        // Figure out the exact constants...
        prj.grid_clean( 1 + (int)(3.0 * sqrt(2) / eps),
                        out_arr );
        
        TreeDiamAlg2d  * p_alg;

        p_alg = new TreeDiamAlg2d( out_arr );
        p_alg->compute_by_heap( eps / 2.0 );
        
        PointPair2d  diam2d = p_alg->getDiameter();

        PointPair   new_diam;
        new_diam.init( arr[ diam2d.p.getID() ],
                       arr[ diam2d.q.getID() ] );
        //printf( "id: %d, %d\n", diam2d.p.getID(), diam2d.q.getID() );

        diam.update_diam( new_diam );
        p_alg->term();
        delete  p_alg;
    }

    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// DiameterChan - 
//     Describes the algorithm due to Timothy Chan Socg2000
//------------------------------------------------------------------

DiameterAlgorithm  * DiameterChan::CreateAndCompute(
                                              ArrPoint3d  & arr,
                                              double  _eps )
{
    DiameterChan  * pDiam = new DiameterChan;

    pDiam->eps = _eps;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * DiameterChan::getID()
{
    static char  * name = "Chan's algorithm";

    return  name;
}


void   DiameterChan::compute( const ArrPoint3d  & arr )
{
    Timer   tm;
    GridCleaner   grid;
    double   delta, alpha, beta;
    int   alpha_limit, beta_limit;

    diam.init( arr[ 0 ] );

    tm.start();

    //for  ( int  ind = 0; ind < arr.size(); ind++ ) 
    //    printf( "___id: %d\n", arr[ ind ].getID() );
        
    grid.init( arr, eps / 3.0 );


    //for  ( int  ind = 0; ind < grid.out_arr.size(); ind++ ) 
    //    printf( "x___id: %d\n", grid.out_arr[ ind ].getID() );

    //printf( "in_size: %d\nout_size: %d\n", 
    //        arr.size(), grid.out_arr.size() );
    
    delta = sqrt( eps / 3.0 );
    alpha_limit = (int)(3.1415926535 / delta);
    beta_limit = (int)(3.1415926535 / delta);

    alpha = 0.0;
    for  ( int  ind_alpha = 0; ind_alpha < alpha_limit; 
           ind_alpha++, alpha += delta ) {
        ProjPoints  prj;
        Point3d  dir;
        ArrPoint2dPtr  out_arr;

        dir[ 0 ] = sin( alpha );
        dir[ 1 ] = cos( alpha );
        dir[ 3 ] = 0;
        
        prj.init( grid.out_arr, dir );

        // Figure out the exact constants...
        prj.grid_clean( 1 + (int)(3.0 * sqrt(2) / eps),
                        out_arr );        
        
        beta = 0;
        for  ( int  ind_beta = 0; ind_beta < beta_limit; 
               ind_beta++, beta += delta ) {
            Point2d  pnt2d;

            pnt2d[ 0 ] = sin( beta );
            pnt2d[ 1 ] = cos( beta );
            
            PointPair2d  diam_tmp;
            //find_extreme_diameter_in_dir( out_arr, pnt2d, diam_tmp );
            computeExtremePairDir2d( out_arr, pnt2d, diam_tmp );

            PointPair   new_diam;
            new_diam.init( arr[ diam_tmp.p.getID() ],
                           arr[ diam_tmp.q.getID() ] );
            diam.update_diam( new_diam );
        }
    }

    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// DiameterPCA - 
//       Diameter using Principal component analysis.
//------------------------------------------------------------------


DiameterAlgorithm  * DiameterPCA::CreateAndCompute( ArrPoint3d  & arr )
{
    DiameterPCA  * pDiam = new DiameterPCA;

    pDiam->eps = sqrt(3);
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * DiameterPCA::getID()
{
    static char  * name = "PCA";

    return  name;
}


void   DiameterPCA::compute( const ArrPoint3d  & arr )
{
    Timer   tm;
    Point3d  base_pnt, side_0, side_1, side_2;

    diam.init( arr[ 0 ] );

    tm.start();
    
    MakeCovarOBB( arr, base_pnt, side_0, side_1, side_2 );
    
    PointPair  pp_0, pp_1, pp_2;
    
    diam.init( arr[ 0 ] );
    tm.start();

    computeExtremePairDirs( arr, 
                            side_0, side_1, side_2, 
                            pp_0, pp_1, pp_2 );
    diam.update_diam( pp_0 );
    diam.update_diam( pp_1 );
    diam.update_diam( pp_2 );
        /*
    computeExtremePairDir( arr, side_0, pp );
    diam.update_diam( pp );

    computeExtremePairDir( arr, side_1, pp );
    diam.update_diam( pp );

    computeExtremePairDir( arr, side_2, pp );
    diam.update_diam( pp );
        */

    tm.end();
    //printf( "PCA: \n" );
    //side_0.print();
    //printf( "\n" );
    //side_1.print();
    ////printf( "\n" );
    //side_2.print();