algorithm.c

任意给定三维空间的点集

    //printf( "\n" );
 
    time = tm.seconds();    
}


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


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

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

    return  pDiam;  //(DiameterAlgorithm *)
}


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

    return  name;
}


void   DiameterMassPCA::compute( const ArrPoint3d  & arr )
{
    Timer   tm;
    //Point3d  base_pnt, side_0, side_1, side_2;
    FILE  * fl;
    Point3d  base_pnt, side_0, side_1, side_2;
    int  args;
    double  x, y, z;
    double  mCoVars[3][3];

    printf( "void   DiameterMassPCA::compute( const ArrPoint3d  & arr )\n" );
    tm.start();

    fl = fopen( "points_out.txt", "wt" );
    assert( fl != NULL );

    fprintf( fl, "3\n%d", arr.size() );
    for  ( int  ind = 0; ind < arr.size(); ind++ ) {
        fprintf( fl, "\n%.20f %.20f %.20f ",
                 arr[ ind ][ 0 ],
                 arr[ ind ][ 1 ],
                 arr[ ind ][ 2 ] );
    }
    fprintf( fl, "\n" );
    fclose( fl );
    system( "./qhuller points_out.txt" );

    // Read the Mass PCA results
    fl = fopen( "out2.txt", "rt" );
    assert( fl != NULL );


    args = fscanf( fl, "%lg %lg %lg\n", &x, &y, &z );
    assert( args == 3 );

    //Point3d  mass_center( x, y, z );

    args = fscanf( fl, "%lg %lg %lg\n", &x, &y, &z );
    assert( args == 3 );
                         
    mCoVars[ 0 ][ 0 ] = x;
    mCoVars[ 0 ][ 1 ] = y;
    mCoVars[ 0 ][ 2 ] = z;

    args = fscanf( fl, "%lg %lg %lg\n", &x, &y, &z );
    assert( args == 3 );
    mCoVars[ 1 ][ 0 ] = x;
    mCoVars[ 1 ][ 1 ] = y;
    mCoVars[ 1 ][ 2 ] = z;

    args = fscanf( fl, "%lg %lg %lg\n", &x, &y, &z );
    assert( args == 3 );
    mCoVars[ 2 ][ 0 ] = x;
    mCoVars[ 2 ][ 1 ] = y;
    mCoVars[ 2 ][ 2 ] = z;

    fclose( fl );

    double  max_entry;

    max_entry = 0;
    for  ( int  ind = 0; ind < 3; ind++ ) 
        for  ( int  jnd = 0; jnd < 3; jnd++ )
            max_entry = max( max_entry, fabs( mCoVars[ ind ][ jnd ] ) );

    if  ( max_entry != 0.0 ) {
        for  ( int  ind = 0; ind < 3; ind++ ) 
            for  ( int  jnd = 0; jnd < 3; jnd++ ) {
                mCoVars[ ind ][ jnd ] = mCoVars[ ind ][ jnd ] / max_entry;
                if  ( fabs( mCoVars[ ind ][ jnd ] ) < 1e-8 ) 
                    mCoVars[ ind ][ jnd ] = 0;
            }
    }


    //printf( "\nCenter of mass: " );
    //mass_center.print();
    //printf( "\n" );

    /*
    for  ( int  ind = 0; ind < 3; ind++ ) {
        printf( "--- [ " );
        for  ( int  jnd = 0; jnd < 3; jnd++ )
            printf( "%g ", mCoVars[ ind ][ jnd ] );
        printf( "]\n" );
    }
    */

    Point3d  vec1, vec2, vec3;
    EigenSymetric3vec( mCoVars, vec1, vec2, vec3 );

    vec1.normalize();
    vec2.normalize();
    vec3.normalize();
    //extractOrthonormalBase( vec1, vec2, vec3 );

    /*
    printf( "\n\n----------------- directions \n" );
    vec1.print();
    printf( "\n" );
    vec2.print();
    printf( "\n" );
    vec3.print();
    printf( "\n" );
    */
    computeParallelogram( arr, vec1, vec2, vec3,
                          base_pnt, side_0, side_1, side_2 );

    diam.init( arr[ 0 ] );
    /*
    printf( "sides...\n" );
    side_0.print();
    printf( "\n" );
    side_1.print();
    printf( "\n" );
    side_2.print();
    printf( "\n" );
    printf( "\n" );
    */
    PointPair  pp;

    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();

    time = tm.seconds();    
}

#define  LINE_SIZE  1024

static void   get_line( FILE  * fl, char  * line )
{
    char  * out;

    out = fgets( line, 1024, fl );
    assert( out != NULL );
}

static  bool  prefix_eq( char  * line, char  * prefix )
{
    while  ( *prefix ) {
        if  ( *line != *prefix )
            return  false;
        line++;
        prefix++;
    }
    return  true;
}


static  bool  suffix_eq( const char  * line, const char  * suffix )
{
    int  len_line, len_suffix;

    len_line = strlen( line );
    len_suffix = strlen( suffix );
    if  ( len_suffix == 0 )
        return  true;
    if  ( len_suffix > len_line )
        return  false;

    return ( strcmp( line + len_line - len_suffix, suffix ) == 0 );
}


void  read_poly_file( ArrPoint3d  & arr,
                      const char  * file_name )
{
    FILE  * fl;
    char  line[ LINE_SIZE ];
    int  vertices, args;

    fl = fopen( file_name, "rt" );
    if  ( fl == NULL ) {
        fprintf( stderr, "Unable to open file: [%s]\n",
                 file_name );
        exit( -1 );
    }

    while  ( ( ! feof( fl ) )  && 
             ( ! prefix_eq( line, "element vertex" ) ) ) 
        get_line( fl, line );
    if  ( feof( fl ) ) {
        fprintf( stderr, "Vertex num info not found!\n" );
        exit( -1 );
    }
    args = sscanf( line, "element vertex %d", &vertices );
    if  ( args < 1 ) {
        fprintf( stderr, "Vertex num info non decipherable!\n" );
        exit( -1 );
    }
    assert( vertices > 0 );
    
    // skip till end of header
    while  ( ( ! feof( fl ) )  && 
             ( ! prefix_eq( line, "end_header" ) ) ) 
        get_line( fl, line );
    if  ( feof( fl ) ) {
        fprintf( stderr, "header termination not found!\n" );
        exit( -1 );
    }

    arr.init( vertices + 10, 3454 );

    // read the vertices
    for  ( int  ind = 0; ind < vertices; ind++ ) {
        double  x, y, z;

        get_line( fl, line );
        if  ( feof( fl ) ) {
            fprintf( stderr, "unexpected end of file!\n" );
            exit( -1 );
        }
        args = sscanf( line, "%lg %lg %lg", &x, &y, &z );
        Point3d  pnt;
        pnt.setCoord( 0, x  );
        pnt.setCoord( 1, y );
        pnt.setCoord( 2, z );

        arr.pushx( pnt );
    }
}

PointPair  * tester_inner( ArrPoint3d   * p_arr,
                           const char  * description )
{
    DiameterAlgorithm  * alg;

    if  ( p_arr == NULL )
        return  NULL;

    TesterDiamAlg  tester;

    printf( "\n\n" );
    printf( "Performing extended tests. Thist might take a while...\n" );
    printf( "Points: %d\n", p_arr->size() );
    tester.init();
    tester.points_num = p_arr->size();

    //-compute_mvbb( *p_arr );

    //printf( "XXXXQuad by heap gdiam ptr ...\n" );
    //tester.push( alg = QuadDiameterByGPtrHeap::CreateAndCompute( 
    //                                               *p_arr, 0.0 ) );
    //-----------------------------------------------------------------
    // 0.0
    //-----------------------------------------------------------------
    printf( "Quad by heap...\n" );
    tester.push( alg = QuadDiameterByHeap::CreateAndCompute( *p_arr, 0.0 ) );

    printf( "Quad by heap gdiam ptr ...\n" );
    tester.push( alg = QuadDiameterByGPtrHeap::CreateAndCompute( 
                                                   *p_arr, 0.0 ) );
    //printf( "Quad by heap gdiam ptr UDM ...\n" );
    //tester.push( alg = QuadDiameterByGPtrHeapUDMSample::CreateAndCompute( 
    //                                               *p_arr, 0.0 ) );

    printf( "Quad by heap WSPD...\n" );
    tester.push( alg = QuadDiameterByHeapWSPD::CreateAndCompute( *p_arr,
                                                                 0.0 ) );
    printf( "Quad by level_ptr 0.0...\n" );
    tester.push( QuadDiameterPtrByLevel::CreateAndCompute( *p_arr, 0.0 ) );

    printf( "Quad by level...\n" );
    tester.push( QuadDiameterByLevel::CreateAndCompute( *p_arr, 0.0 ) );

    if   ( p_arr->size() < 4000 ) {
        printf( "Brute force...\n" );
        tester.push( BruteDiameter::CreateAndCompute( *p_arr ) );

        printf( "Brute force local...\n" );
        tester.push( BruteDiameterLocal::CreateAndCompute( *p_arr ) );
    }

    //-----------------------------------------------------------------
    // 0.01
    //-----------------------------------------------------------------
    printf( "Quad by heap (0.01)...\n" );
    tester.push( alg = QuadDiameterByHeap::CreateAndCompute( *p_arr, 0.01 ) );

    printf( "Quad by heap gdiam ptr ...\n" );
    tester.push( alg = QuadDiameterByGPtrHeap::CreateAndCompute( 
                                                   *p_arr, 0.01 ) );

    printf( "Quad by heap WSPD (0.01)...\n" );
    tester.push( alg = QuadDiameterByHeapWSPD::CreateAndCompute( *p_arr,
                                                                 0.01 ) );
    printf( "Quad by level_ptr 0.01...\n" );
    tester.push( QuadDiameterPtrByLevel::CreateAndCompute( *p_arr, 0.01 ) );

    printf( "Quad by level (0.01)...\n" );
    tester.push( QuadDiameterByLevel::CreateAndCompute( *p_arr, 0.01 ) );
    fflush( stdout );

    printf( "Chan algorithm (0.01)...\n" );
    tester.push( DiameterChan::CreateAndCompute( *p_arr, 0.01 ) );
    fflush( stdout );

    printf( "Chan modified algorithm (0.01)...\n" );
    tester.push( DiameterModChan::CreateAndCompute( *p_arr, 0.01 ) );
    fflush( stdout );

    /*
    if  ( p_arr->size() <= 3000 ) {
        printf( "Quad by directions (0.01)...\n" );
        tester.push( DiameterDirections::CreateAndCompute( *p_arr, 0.01 ) );
        fflush( stdout );
        }*/
    printf( "Quad by grid (0.01)...\n" );
    tester.push( DiameterGrid::CreateAndCompute( *p_arr, 0.01 ) );

    //printf( "Quad by grid + fs dir (0.01)...\n" );
    //tester.push( DiameterGridQuadDir::CreateAndCompute( *p_arr, 0.01 ) );

    printf( "Quad by dir (0.01)...\n" );
    tester.push( QuadDiameterByDir::CreateAndCompute( *p_arr, 0.01 ) );
    fflush( stdout );

    //-----------------------------------------------------------------
    // 0.1
    //-----------------------------------------------------------------
    printf( "Quad by heap (0.1)...\n" );
    tester.push( alg = QuadDiameterByHeap::CreateAndCompute( *p_arr, 0.1 ) );

    printf( "Quad by heap gdiam ptr ...\n" );
    tester.push( alg = QuadDiameterByGPtrHeap::CreateAndCompute( 
                                                   *p_arr, 0.1 ) );

    printf( "Quad by heap WSPD (0.1)...\n" );
    tester.push( alg = QuadDiameterByHeapWSPD::CreateAndCompute( *p_arr, 
                                                                 0.1 ) );

    printf( "Quad by level_ptr 0.1...\n" );
    tester.push( QuadDiameterPtrByLevel::CreateAndCompute( *p_arr, 0.1 ) );

    printf( "Quad by level (0.1)...\n" );
    tester.push( QuadDiameterByLevel::CreateAndCompute( *p_arr, 0.1 ) );
    printf( "Chan algorithm (0.1)...\n" );
    tester.push( DiameterChan::CreateAndCompute( *p_arr, 0.1 ) );
    fflush( stdout );

    printf( "Chan modified algorithm (0.1)...\n" );
    tester.push( DiameterModChan::CreateAndCompute( *p_arr, 0.1 ) );
    fflush( stdout );

    /*
    if  ( p_arr->size() <= 3  0 ) {
        printf( "Quad by directions (0.1)...\n" );
        tester.push( DiameterDirections::CreateAndCompute( *p_arr, 0.1 ) );
        fflush( stdout );
        }*/

    printf( "Quad by grid (0.1)...\n" );
    tester.push( DiameterGrid::CreateAndCompute( *p_arr, 0.1 ) );
    printf( "Quad by grid + fs dir (0.1)...\n" );
    tester.push( DiameterGridQuadDir::CreateAndCompute( *p_arr, 0.1 ) );
    printf( "Quad by dir (0.1)...\n" );
    fflush( stdout );
    tester.push( QuadDiameterByDir::CreateAndCompute( *p_arr, 0.1 ) );

    //-----------------------------------------------------------------
    // 1.72
    //-----------------------------------------------------------------
    //if  ( p_arr->size() <= 30000 ) {
        //printf( "PCA Mass algorithm ...\n" );
        //tester.push( DiameterMassPCA::CreateAndCompute( *p_arr ) );
        //fflush( stdout );
    //}
    printf( "Diameter by axis parallel bounding box\n" );
    fflush( stdout );
    tester.push( DiameterAPBBox::CreateAndCompute( *p_arr ) );

    printf( "PCA algorithm ...\n" );
    tester.push( DiameterPCA::CreateAndCompute( *p_arr ) );
    fflush( stdout );

    printf( "Points: %d\n", p_arr->size() );
    tester.dump_to_file( "out.txt", description );
    tester.dump();

    fflush( stdout );
    printf( "\n\nTest done!\n\n" );

    return   new  PointPair( alg->diam );
}

void  tester( const char  * file_name,
              const char  * description )
{
    ArrPoint3d   * p_arr;

    p_arr = new ArrPoint3d;
    p_arr->init( 10, 43283 );

    //p_arr->reset();

    if  ( ( suffix_eq( file_name, ".3ds" ) )
          ||  ( suffix_eq( file_name, ".3DS" ) ) )
        read_3ds( file_name, *p_arr );
    else
        read_poly_file( *p_arr, file_name );

    pnt_array_normalize( *p_arr );
    
    tester_inner( p_arr, description );

    delete p_arr;
}



void   compute_mvbb( const ArrPoint3d  & arr )
{
    Timer  tm;
    gdiam_real  * p_arr, * p_start;
    GBBox  ap_bbox;
    //double  num;

    tm.start();

    p_arr = (gdiam_real *)malloc( sizeof( gdiam_real ) * arr.size() * 3 );
    p_start = p_arr;
    for  ( int  ind = 0; ind < arr.size(); ind++ ) {
        p_arr[ 0 ] = arr[ ind ][ 0 ];
        p_arr[ 1 ] = arr[ ind ][ 1 ];
        p_arr[ 2 ] = arr[ ind ][ 2 ];
        p_arr += 3;
    }

    
    gdiam_bbox  bb, bb1;
    gdiam_point  * ptr_arr;

    Timer  tm_reg;

    tm_reg.start();

    ptr_arr = gdiam_convert( p_start, arr.size() );
    bb = gdiam_approx_const_mvbb( ptr_arr, arr.size(), 0, &ap_bbox );

    tm_reg.end();

    Timer  tm_reg2;

    tm_reg2.start();

    //bb1 = gdiam_approx_mvbb_grid( ptr_arr, arr.size(), 5 );
    bb1 = gdiam_approx_mvbb_grid_sample( ptr_arr, arr.size(), 
                                         5, 400 );
    
    tm_reg2.end();

    free( ptr_arr );
    free( p_start );
    
    tm.end();

    printf( "Overall time: %g\n", tm.seconds() );
    printf( "Time for const approx   : %g\n", tm_reg.seconds() );

    printf( "Volume on axis para bbox  : %g\n", ap_bbox.volume() );
            
    printf( "Volume of const approx    : %g\n", bb.volume() );
    printf( "Volume of grid            : %g\n", bb1.volume() );
    printf( "\n\n\n\n" );
    fflush( stdout );
}



/* algorithm.C - End of File ------------------------------------------*/