treediam.c

任意给定三维空间的点集

                     pair.get_right()->get_right() );
        return;
    }
}

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

    //printf( "compute_by_dir: %.4f\n", eps );
    // if a pair have 2r/l <= threshold then the error is small
    // enough, and one can use the direction between the centers to do
    // the computations.
    threshold = sqrt( 1.0 - (1.0 - eps) * (1.0 - eps) );
    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 );
    
    int  count =0;
    while  ( heap.size() > 0 ) {
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        count++;
        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_dir( pair, heap, eps, threshold );
        
        heap.pop();
    }
}


void  TreeDiamAlg::addPair( FSPairArray  & queue_out, 
                            FSTreeNode  * left,
                            FSTreeNode  * right )
{
    const Point3d  & p( tree.getPoint( left->index_left() ) );
    const Point3d  & q( tree.getPoint( right->index_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 );
    
    FSPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    queue_out.pushx( pair );
}


void  TreeDiamAlg::split_level( FSPairArray  & queue_in,
                                FSPairArray  & queue_out, 
                                double  eps )
{
    bool  f_is_left_splitable, f_is_right_splitable;
    
    for  ( int  ind = 0; ind < queue_in.size(); ind++ ) {
        FSPair  & 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  TreeDiamAlg::split_pair( FSPair  & pair,
                               heap_pairs  & heap, double  eps )
{
    bool  f_is_left_splitable, f_is_right_splitable;
    
    if  ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
        return;
    
    //printf( "threshold_brute: %d\n", threshold_brute );
    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  TreeDiamAlg::split_pair_wspd( FSPair  & pair,
                                    heap_pairs  & heap, 
                                    double  eps,
                                    int  & count_self_pairs )
{
    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 ) {
        addPairHeapWSPD( heap, 
                         pair.get_left()->get_left(),
                         pair.get_right()->get_left(),
                         eps,
                         count_self_pairs );
        addPairHeapWSPD( heap, 
                         pair.get_left()->get_left(),
                         pair.get_right()->get_right(),
                         eps,
                         count_self_pairs  );
        // to avoid exponential blowup
        if  ( pair.get_left() != pair.get_right() )
            addPairHeapWSPD( heap, 
                             pair.get_left()->get_right(),
                             pair.get_right()->get_left(),
                             eps,
                             count_self_pairs  );
        addPairHeapWSPD( heap, 
                         pair.get_left()->get_right(),
                         pair.get_right()->get_right(),
                         eps,
                         count_self_pairs  );
        return;
    }
    if  ( f_is_left_splitable ) {
        addPairHeapWSPD( heap,
                         pair.get_left()->get_left(),
                         pair.get_right(),
                         eps,
                         count_self_pairs  );
        addPairHeapWSPD( heap, 
                         pair.get_left()->get_right(),
                         pair.get_right(),
                         eps,
                         count_self_pairs  );
        return;
    }
    if  ( f_is_right_splitable ) {
        addPairHeapWSPD( heap, 
                         pair.get_left(),
                         pair.get_right()->get_left(),
                         eps,
                         count_self_pairs  );
        addPairHeapWSPD( heap, 
                         pair.get_left(),
                         pair.get_right()->get_right(),
                         eps,
                         count_self_pairs  );
        return;
    }
}



// Splits a pair only one side - avoiding a huge number of pairs...?
void  TreeDiamAlg::split_pair_wspd_one_side( FSPair  & pair,
                                             heap_pairs  & heap, 
                                             double  eps,
                                             int  & count_self_pairs )
{
    bool  f_is_left_splitable;

    if  ( pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance ))
        return;
    
    if  ( pair.isBelowThreshold( threshold_brute ) ) {
        pair.directDiameter( tree, pair_diam );
        return;
    }
    
    pair.setLeftBigger();

    f_is_left_splitable = pair.get_left()->isSplitable();
    assert( f_is_left_splitable );
    if ( f_is_left_splitable )
        tree.split_node( pair.get_left() );
    
    addPairHeapWSPD( heap, 
                     pair.get_left()->get_left(),
                     pair.get_right(),
                     eps,
                     count_self_pairs );
    addPairHeapWSPD( heap, 
                     pair.get_left()->get_right(),
                     pair.get_right(),
                     eps,
                     count_self_pairs  );
}


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

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

    return  diam.distance;
}

void  TreeDiamAlg::addPairHeap( heap_pairs  & heap,
                             FSTreeNode  * left,
                             FSTreeNode  * right )
{
    const Point3d  & p( tree.getPoint( left->index_left() ) );
    const Point3d  & 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 );
    
    FSPair  pair;
    
    pair.init( left, right );
    //printf( "__pair max diameter: %g\n", pair.maxDiam() );
    
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    //printf( "pair max diameter: %g\n", pair.maxDiam() );
    heap.push( pair );
}


void  TreeDiamAlg::addPairHeapWSPD( heap_pairs  & heap,
                                    FSTreeNode  * left,
                                    FSTreeNode  * right,
                                    double  eps,
                                    int  & count_self_pairs
                                    )
{
    const Point3d  & p( tree.getPoint( left->index_left() ) );
    const Point3d  & 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 );
    
    FSPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= (1.0 + eps) * pair_diam.distance )
        return;
    //printf( "iis: %d, pair diam: %10g  estimate: %g\n", 
    //        (int)( pair.isIdenticalSides() ),
    //        pair.maxDiam(),
    //        pair_diam.distance );
    heap.push( pair );
    if  ( pair.isIdenticalSides() )
        count_self_pairs++;
}


void  pnt_array_bound( ArrPoint3d  & arr, BBox3d  & bb )
{
    bb.init();

    for  ( int  ind = 0; ind < arr.size(); ind++ ) 
        bb.bound( arr[ ind ] );
}


void  FSTree::findExtremeInProjection( FSPair  & pair, 
                                       PointPair  & max_pair_diam )
{
    PointPair  candid;
    Point3d  dir;
    
    dir = pnt_subtract( pair.get_left()->getCenter(), 
                        pair.get_right()->getCenter() );
    
    candid.init( arr[ pair.get_left()->index_left() ] );
    computeExtremePairDir( arr, dir, candid, 
                           pair.get_left()->index_left(),
                           pair.get_left()->index_right() );
    computeExtremePairDir( arr, dir, candid, 
                           pair.get_right()->index_left(),
                           pair.get_right()->index_right() );
    max_pair_diam.update_diam( candid );
}


//--------------------------------------------------------------------------
// TreeDiamAlg2d
//--------------------------------------------------------------------------

void  computeExtremePair( const ArrPoint2dPtr  & arr,
                          int  dim, PointPair2d  & pp )
{
    pp.init( *(arr[ 0 ]) );

    for  ( int  ind = 1; ind < arr.size(); ind++ ) {
        const Point2d  & 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 );
}

class  heap_pairs_2d :
          public std::priority_queue<FS2dPair, 
          std::vector<FS2dPair>, lt_fspair_2d>
{