treediam.c

任意给定三维空间的点集

    limit = num / 2;
    for  ( int  ind = 0; ind < limit; ind++ ) {
        addPoint( arr, -1, -1, -1 );
        addPoint( arr,  1,  1,  1 );
    }
}


void  computeExtremePairDir( const ArrPoint3d  & arr,
                             const Point3d  & dir,
                             PointPair  & pp
                             )
{
    real  min_dot, max_dot, prod;
    int  ind_min_dot, ind_max_dot;

    pp.init( arr[ 0 ] );
    min_dot = max_dot = dir.dot_prod( arr[ 0 ] );
    ind_min_dot = ind_max_dot = 0;

    for  ( int  ind = 1; ind < arr.size(); ind++ ) {
        prod = dir.dot_prod( arr[ ind ] );
        if  ( prod < min_dot ) {
            min_dot = prod;
            ind_min_dot = ind;
        }
        if  ( prod > max_dot ) {
            max_dot = prod;
            ind_max_dot = ind;
        }
    }

    pp.init( arr[ ind_min_dot ], arr[ ind_max_dot ] );
}


void  computeExtremePairDirs( const ArrPoint3d  & arr,
                              const Point3d  & dir_0,
                              const Point3d  & dir_1,
                              const Point3d  & dir_2,
                              PointPair  & pp_0,
                              PointPair  & pp_1,
                              PointPair  & pp_2
                             )
{
    real  min_dot_0, max_dot_0, prod_0;
    int  ind_min_dot_0, ind_max_dot_0;
    real  min_dot_1, max_dot_1, prod_1;
    int  ind_min_dot_1, ind_max_dot_1;
    real  min_dot_2, max_dot_2, prod_2;
    int  ind_min_dot_2, ind_max_dot_2;

    pp_0.init( arr[ 0 ] );
    pp_1.init( arr[ 0 ] );
    pp_2.init( arr[ 0 ] );
    min_dot_0 = max_dot_0 = dir_0.dot_prod( arr[ 0 ] );
    ind_min_dot_0 = ind_max_dot_0 = 0;
 
    min_dot_1 = max_dot_1 = dir_1.dot_prod( arr[ 0 ] );
    ind_min_dot_1 = ind_max_dot_1 = 0;

    min_dot_2 = max_dot_2 = dir_2.dot_prod( arr[ 0 ] );
    ind_min_dot_2 = ind_max_dot_2 = 0;

    for  ( int  ind = 1; ind < arr.size(); ind++ ) {
        const Point3d  & pnt( arr[ ind ] );

        prod_0 = dir_0.dot_prod( pnt );
        if  ( prod_0 < min_dot_0 ) {
            min_dot_0 = prod_0;
            ind_min_dot_0 = ind;
        }
        if  ( prod_0 > max_dot_0 ) {
            max_dot_0 = prod_0;
            ind_max_dot_0 = ind;
        }

        prod_1 = dir_1.dot_prod( pnt );
        if  ( prod_1 < min_dot_1 ) {
            min_dot_1 = prod_1;
            ind_min_dot_1 = ind;
        }
        if  ( prod_1 > max_dot_1 ) {
            max_dot_1 = prod_1;
            ind_max_dot_1 = ind;
        }

        prod_2 = dir_2.dot_prod( pnt );
        if  ( prod_2 < min_dot_2 ) {
            min_dot_2 = prod_2;
            ind_min_dot_2 = ind;
        }
        if  ( prod_2 > max_dot_2 ) {
            max_dot_2 = prod_2;
            ind_max_dot_2 = ind;
        }
    }

    pp_0.init( arr[ ind_min_dot_0 ], arr[ ind_max_dot_0 ] );
    pp_1.init( arr[ ind_min_dot_1 ], arr[ ind_max_dot_1 ] );
    pp_2.init( arr[ ind_min_dot_2 ], arr[ ind_max_dot_2 ] );
}


void  computeExtremePairDir( const ArrPoint3d  & arr,
                             const Point3d  & dir,
                             PointPair  & pp,
                             int  start,
                             int  end
                             )
{
    real  min_dot, max_dot, prod;
    const Point3d  *p_pnt_min, *p_pnt_max;

    p_pnt_min = &pp.p;
    p_pnt_max = &pp.q;

    if  ( dir.dot_prod( *p_pnt_min ) > dir.dot_prod( *p_pnt_max ) )
        exchange( p_pnt_min, p_pnt_max );

    min_dot = dir.dot_prod( *p_pnt_min );
    max_dot = dir.dot_prod( *p_pnt_max );
 
    for  ( int  ind = start; ind < arr.size(); ind++ ) {
        prod = dir.dot_prod( arr[ ind ] );
        if  ( prod < min_dot ) {
            min_dot = prod;
            p_pnt_min = &( arr[ ind ] );
        }
        if  ( prod > max_dot ) {
            max_dot = prod;
            p_pnt_max = &( arr[ ind ] );
        }
    }
    Point3d  p, q;

    p = *p_pnt_max;
    q = *p_pnt_min;
    pp.init( p, q );
}



void  computeExtremePair( const ArrPoint3d  & arr,
                          int  dim, PointPair  & pp )
{
    pp.init( arr[ 0 ] );

    for  ( int  ind = 1; ind < arr.size(); ind++ ) {
        const Point3d  & pnt( arr[ ind ] );

        if  ( pnt[ dim ] < pp.p[ dim ] )
            pp.p = pnt;
        if  ( pnt[ dim ] > pp.q[ dim ] )
            pp.q = pnt;
    }

    pp.distance = pnt_distance( pp.p, pp.q );
}


void  computeExtremePairAllDim( const ArrPoint3d  & arr,
                                PointPair  & pp0,
                                PointPair  & pp1,
                                PointPair  & pp2
                                )
{
    pp0.init( arr[ 0 ] );
    pp1.init( arr[ 0 ] );
    pp2.init( arr[ 0 ] );

    for  ( int  ind = 1; ind < arr.size(); ind++ ) {
        const Point3d  & pnt( arr[ ind ] );
        
        if  ( pnt[ 0 ] < pp0.p[ 0 ] )
            pp0.p = pnt;
        if  ( pnt[ 0 ] > pp0.q[ 0 ] )
            pp0.q = pnt;

        if  ( pnt[ 1 ] < pp1.p[ 1 ] )
            pp1.p = pnt;
        if  ( pnt[ 1 ] > pp1.q[ 1 ] )
            pp1.q = pnt;

        if  ( pnt[ 2 ] < pp2.p[ 2 ] )
            pp2.p = pnt;
        if  ( pnt[ 2 ] > pp2.q[ 2 ] )
            pp2.q = pnt;
    }

    pp0.distance = pnt_distance( pp0.p, pp0.q );
    pp1.distance = pnt_distance( pp1.p, pp1.q );
    pp2.distance = pnt_distance( pp2.p, pp2.q );
}


void  computeExtremePair( CPoint3dPtr  arr,
                          int  size,
                          int  dim, PointPair  & pp )
{
    //printf( "arr: %p size: %d\n", arr, size );

    pp.init( arr[ 0 ] );

    for  ( int  ind = 1; ind < size; ind++ ) {
        const Point3d  & pnt( arr[ ind ] );

        if  ( pnt[ dim ] < pp.p[ dim ] )
            pp.p = pnt;
        if  ( pnt[ dim ] > pp.q[ dim ] )
            pp.q = pnt;
    }

    pp.distance = pnt_distance( pp.p, pp.q );
}


void  TreeDiamAlg::compute_by_level( double  eps ) 
{
    FSPairArray  queue_a, queue_b;
    FSPairArray  * p_queue_in, * p_queue_out;
    int  heap_limit, heap_delta;
    
    queue_a.init( 100, 43432 );
    queue_b.init( 100, 43431 );
    
    FSTreeNode  * root = tree.getRoot();

    computeExtremePair( tree.getPoints(), root->getBBox().getLongestDim(),
                        pair_diam );

    FSPair  root_pair;
    
    root_pair.init( root, root );
    
    queue_a.pushx( root_pair );
    
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
    p_queue_in = &queue_a;
    p_queue_out = &queue_b;
    for  ( int  ind = 0; ind < 1000; ind++ ) {
        //printf( "ind: %d\n", ind );
        if   ( p_queue_in->isEmpty() ) {
            //printf( "breaking!\n" );
            break;
        }
        p_queue_out->reset();
        if  ( p_queue_in->size() > heap_limit ) {
                threshold_brute *= 4;
                printf( "threshold_brute: %d\n", threshold_brute );
                heap_limit += heap_delta;
        }

        split_level( *p_queue_in, *p_queue_out, eps );
        
        exchange( p_queue_in, p_queue_out );
    }
    
}


void  TreeDiamAlg::compute_by_heap( double  eps ) 
{
    heap_pairs  heap;
    int  heap_limit, heap_delta;
    
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;

    printf( "lll\n" );
    FSTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), root->getBBox().getLongestDim(),
                        pair_diam );
    //diam = root->getBBox().getLongestEdge();
    
    FSPair  root_pair;
    
    root_pair.init( root, root );
    
    printf( "root diam: %g\n", root_pair.maxDiam() );
    heap.push( root_pair );
    //queue_a.pushx( root_pair );
    
    int  count =0;
    while  ( heap.size() > 0 ) {
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        FSPair  pair = heap.top();
        split_pair( pair, heap, eps );
    
        if  ( ( count & 0x3ff ) == 0 ) {
            if  ( ((int)heap.size()) > heap_limit ) {
                threshold_brute *= 2;
                printf( "threshold_brute: %d\n", threshold_brute );
                heap_limit += heap_delta;
            }
            //printf( "size: %d  [%d]\n", (int)heap.size(),
            //        (int)sizeof( FSPair ) );
            //mem_report();
            //fflush( stdout );
        }

        count++;
        heap.pop();
    }
    printf( "\tOverall count:%d\n", count );
}


void  TreeDiamAlg::compute_by_heap_wspd( double  eps ) 
{
    heap_pairs  heap;
    int  count_self_pairs;
    int  heap_limit, heap_delta;
    
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;

    FSTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), root->getBBox().getLongestDim(),
                        pair_diam );
    //diam = root->getBBox().getLongestEdge();
    
    FSPair  root_pair;
    
    root_pair.init( root, root );
    
    heap.push( root_pair );
    //queue_a.pushx( root_pair );
    
    //t  count =0;
    int  count = 0;

    count_self_pairs = 1;
    while  ( count_self_pairs > 0 ) {
        //intf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        FSPair  pair = heap.top();

        if  ( ( count & 0x3ff ) == 0 ) {
            if  ( ((int)heap.size()) > heap_limit ) {
                threshold_brute *= 2;
                printf( "threshold_brute: %d\n", threshold_brute );
                heap_limit += heap_delta;
            } 
       }

        if  ( pair.isIdenticalSides() )
            --count_self_pairs;
        split_pair_wspd( pair, heap, eps, count_self_pairs );
        count++;
        
        heap.pop();
    }

    //printf( "no self pairs in heap!\ncount:%d\n", count );

    while  ( heap.size() > 0 ) {
        //intf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        FSPair  pair = heap.top();

        if  ( ( count & 0x3ff ) == 0 ) {
            if  ( ((int)heap.size()) > heap_limit ) {
                threshold_brute *= 2;
                printf( "threshold_brute: %d\n", threshold_brute );
                heap_limit += heap_delta;
            } 
       }

        ///split_pair_wspd( pair, heap, eps, count_self_pairs );
        split_pair_wspd_one_side( pair, heap, eps, count_self_pairs );
        count++;

        heap.pop();
    }
    printf( "\tOverall count:%d\n", count );
}


void  TreeDiamAlg::split_pair_dir( FSPair  & pair,
                                   heap_pairs  & heap, 
                                   double  eps,
                                   double  threshold )
{
    bool  f_is_left_splitable, f_is_right_splitable;
    
    //printf( "split_pair_dir: (%g), ?, %g, %g\n", 
    //        pair_diam.distance, eps, threshold );
    if  ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
        return;
    
    if  ( pair.isBelowThreshold( threshold_brute ) ) {
        pair.directDiameter( tree, pair_diam );
        return;
    }
    if  ( pair.getProjectionError() <= threshold ) {
        tree.findExtremeInProjection( pair, pair_diam );
        return;
    }

    f_is_left_splitable = pair.get_left()->isSplitable();
    f_is_right_splitable = pair.get_right()->isSplitable();
    assert( f_is_left_splitable  ||  f_is_right_splitable );
    if ( f_is_left_splitable )
        tree.split_node( pair.get_left() );
    if ( f_is_right_splitable )
        tree.split_node( pair.get_right() );
    
    if  ( f_is_left_splitable 
          &&  f_is_right_splitable ) {
        addPairHeap( heap, 
                     pair.get_left()->get_left(),
                     pair.get_right()->get_left() );
        addPairHeap( heap, 
                     pair.get_left()->get_left(),
                     pair.get_right()->get_right() );
        // to avoid exponential blowup
        if  ( pair.get_left() != pair.get_right() )
            addPairHeap( heap, 
                         pair.get_left()->get_right(),
                         pair.get_right()->get_left() );
        addPairHeap( heap, 
                     pair.get_left()->get_right(),
                     pair.get_right()->get_right() );
        return;
    }
    if  ( f_is_left_splitable ) {
        addPairHeap( heap,
                     pair.get_left()->get_left(),
                     pair.get_right() );
        addPairHeap( heap, 
                     pair.get_left()->get_right(),
                     pair.get_right() );
        return;
    }
    if  ( f_is_right_splitable ) {
        addPairHeap( heap, 
                     pair.get_left(),
                     pair.get_right()->get_left() );
        addPairHeap( heap, 
                     pair.get_left(),