treediam.c

任意给定三维空间的点集

};

void  TreeDiamAlg2d::addPairHeap( heap_pairs_2d  & heap,
                                  FS2dTreeNode  * left,
                                  FS2dTreeNode  * right )
{
    const Point2d  & p( tree.getPoint( left->index_left() ) );
    const Point2d  & q( tree.getPoint( right->index_left() ) );

    pair_diam.update_diam( p, q );
    //printf( "old_diam; %g\n", diam );
    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
    
    FS2dPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    heap.push( pair );
}


void  TreeDiamAlg2d::split_pair( FS2dPair  & pair,
                                 heap_pairs_2d  & heap, double  eps )
{
    bool  f_is_left_splitable, f_is_right_splitable;

    //printf( "xx\n" );
    if  ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
        return;
    
    if  ( pair.isBelowThreshold( threshold_brute ) ) {
        pair.directDiameter( tree, 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(),
                     pair.get_right()->get_right() );
        return;
    }
}


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

    FS2dTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), root->getBBox().getLongestDim(),
                        pair_diam );
    
    FS2dPair  root_pair;
    
    root_pair.init( root, root );
    
    heap.push( root_pair );
    //queue_a.pushx( root_pair );
    
    int  count =0;
    while  ( heap.size() > 0 ) {
        //intf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        FS2dPair  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( pair, heap, eps );
        count++;
        
        heap.pop();
    }
}

void  computeExtremePairDir2d( const ArrPoint2dPtr  & arr,
                               const Point2d  & dir,
                               PointPair2d  & 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 ]) );
}


//----------------------------------------------------------------------
// TreeDiamPAlg -
// fair split tree that uses pointers
//----------------------------------------------------------------------

class  heap_pairs_p : public std::priority_queue<FSPPair, 
                      std::vector<FSPPair>, 
                      lt_fsppair>
{
};
//   heap_pairs;
class  FSPPairArray : public  Array<FSPPair> {};

void  TreeDiamPAlg::compute_by_level( double  eps ) 
{
    FSPPairArray  queue_a, queue_b;
    FSPPairArray  * p_queue_in, * p_queue_out;
    int  heap_limit, heap_delta;
 
    //printf( "xa\n" );
    queue_a.init( 100, 43432 );
    queue_b.init( 100, 43431 );
    
    FSPTreeNode  * root = tree.getRoot();

    //printf( "xb\n" );
    computeExtremePair( tree.getPoints(), 
                        points_num,
                        root->getBBox().getLongestDim(),
                        pair_diam );

    FSPPair  root_pair;
    
    root_pair.init( root, root );
    
    //printf( "xc\n" );
    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();
        //printf( "xd\n" );
        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  TreeDiamPAlg::compute_by_heap( double  eps ) 
{
    heap_pairs_p  heap;
    int  heap_limit, heap_delta;
    
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
    
    FSPTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), points_num,
                        root->getBBox().getLongestDim(),
                        pair_diam );
    //diam = root->getBBox().getLongestEdge();
    
    FSPPair  root_pair;
    
    root_pair.init( root, root );
    
    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 );
        FSPPair  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;
            } 
        }
        count++;

        heap.pop();
    }
}


void  TreeDiamPAlg::addPair( FSPPairArray  & queue_out, 
                            FSPTreeNode  * left,
                            FSPTreeNode  * right )
{
    const Point3d  & p( **(left->ptr_pnt_left()) );
    const Point3d  & q( **(right->ptr_pnt_left()) );
    
    //printf( "old_diam; %g\n", diam );
    pair_diam.update_diam( p, q );

    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
    
    FSPPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    queue_out.pushx( pair );
}


void  TreeDiamPAlg::split_level( FSPPairArray  & queue_in,
                                FSPPairArray  & queue_out, 
                                double  eps )
{
    bool  f_is_left_splitable, f_is_right_splitable;
    
    for  ( int  ind = 0; ind < queue_in.size(); ind++ ) {
        FSPPair  & in( queue_in[ ind ] );
        if  ( in.maxDiam() <= (1.0 + eps) * pair_diam.distance )
            continue;
        
        if  ( in.isBelowThreshold( threshold_brute ) ) {
            in.directDiameter( tree, pair_diam );
            continue;
        }
        
        //printf( "\n---curr_diam: %g\n", diam );
        //in.dump();
        //printf( "-------------\n\n" );
        // Oh well, we need to split, if possible...
        f_is_left_splitable = in.get_left()->isSplitable();
        f_is_right_splitable = in.get_right()->isSplitable();
        assert( f_is_left_splitable ||  f_is_right_splitable );
        if ( f_is_left_splitable )
            tree.split_node( in.get_left() );
        if ( f_is_right_splitable )
            tree.split_node( in.get_right() );
        
        if  ( f_is_left_splitable 
              &&  f_is_right_splitable ) {
            addPair( queue_out, 
                     in.get_left()->get_left(),
                     in.get_right()->get_left() );
            addPair( queue_out, 
                     in.get_left()->get_left(),
                     in.get_right()->get_right() );
            // to avoid exponential blowup
            if  ( in.get_left() != in.get_right() )
                addPair( queue_out, 
                         in.get_left()->get_right(),
                         in.get_right()->get_left() );
            addPair( queue_out, 
                     in.get_left()->get_right(),
                     in.get_right()->get_right() );
            continue;
        }
        if  ( f_is_left_splitable ) {
            addPair( queue_out, 
                     in.get_left()->get_left(),
                     in.get_right() );
            addPair( queue_out, 
                     in.get_left()->get_right(),
                     in.get_right() );
            continue;
        }
        if  ( f_is_right_splitable ) {
            addPair( queue_out, 
                     in.get_left(),
                     in.get_right()->get_left() );
            addPair( queue_out, 
                     in.get_left(),
                     in.get_right()->get_right() );
            continue;
        }
    }
}


void  TreeDiamPAlg::split_pair( FSPPair  & pair,
                               heap_pairs_p  & heap, double  eps )
{
    bool  f_is_left_splitable, f_is_right_splitable;
    
    if  ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
        return;
    
    if  ( pair.isBelowThreshold( threshold_brute ) ) {
        pair.directDiameter( tree, 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(),
                     pair.get_right()->get_right() );
        return;
    }
}

/*
real   brute_diameter( const ArrPoint3d  & arr,
                       PointPair  & diam )
{
    diam.init( arr[ 0 ] );

    for  ( int  ind = 0; ind < ( arr.size() - 1 ); ind++ )
        for  ( int  jnd = ind + 1; jnd < arr.size(); jnd++ ) 
            diam.update_diam( arr[ ind ], arr[ jnd ] );

    return  diam.distance;
}*/

void  TreeDiamPAlg::addPairHeap( heap_pairs_p  & heap,
                             FSPTreeNode  * left,
                             FSPTreeNode  * right )
{
    const Point3d  & p( **(left->ptr_pnt_left()) );
    const Point3d  & q( **(right->ptr_pnt_left()) );

    pair_diam.update_diam( p, q );
    //printf( "old_diam; %g\n", diam );
    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
    
    FSPPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    heap.push( pair );
}

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