gdiam.c

任意给定三维空间的点集

    }

    double   directDiameter( const GFSPTree  & tree,
                             GPointPair  & diam,
                             const gdiam_point  dir ) const {
        return  tree.brute_diameter( left->ptr_pnt_left(),
                                     left->ptr_pnt_right(), 
                                     right->ptr_pnt_left(),
                                     right->ptr_pnt_right(),
                                     diam, dir );
    }

    double  maxDiam()  const {
        return  max_diam;
    }
    double  minAprxDiam()  const {
        double   sub_diam;
        
        sub_diam = left->maxDiam() + right->maxDiam();
        if  ( sub_diam > (max_diam / 10.0) )
            return  max_diam;

        return  max_diam - sub_diam / 2.0;
    }
};


class g_heap_pairs_p;

#define  UDM_SIMPLE      1
#define  UDM_BIG_SAMPLE  2 

class  GTreeDiamAlg
{
private:
    GFSPTree  tree;
    //double  diam;
    GPointPair   pair_diam;
    int  points_num;
    int  threshold_brute;
    int  update_diam_mode;
public:
    GTreeDiamAlg( gdiam_point  * arr, int  size, int  mode ) {
        tree.init( arr, size );
        //diam = 0;
        pair_diam.init( arr[ 0 ] );
        points_num = size;
        threshold_brute = 40;
        update_diam_mode = mode;
    }

    void  term() {
        tree.term();
    }
    
    int  size() const {
        return  points_num;
    }
    const GPointPair  & getDiameter() const {
        return  pair_diam;
    }

    int  nodes_number() const {
        return  tree.nodes_number(); 
    }

    void  addPairHeap( g_heap_pairs_p  & heap,
                       GFSPTreeNode  * left,
                       GFSPTreeNode  * right );
    void  addPairHeap( g_heap_pairs_p  & heap,
                       GFSPTreeNode  * left,
                       GFSPTreeNode  * right,
                       gdiam_point  proj,
                       GFSPPair  & father );
    void  split_pair( GFSPPair  & pair,
                      g_heap_pairs_p  & heap, double  eps );
    void  split_pair_proj( GFSPPair  & pair,
                           g_heap_pairs_p  & heap, double  eps,
                           gdiam_point  proj );
    void  compute_by_heap( double  eps );
    void  compute_by_heap_proj( double  eps, gdiam_point  proj );
 
    double  diameter() {
        return  pair_diam.distance;
    }
};


/*--- Heap implementation ---*/

typedef  int   ( *ptrCompareFunc )( void  * aPtr, void  * bPtr );
typedef  void  * voidPtr_t;

struct  heap_t
{
    voidPtr_t  * pArr;
    int  curr_size, max_size;   
    ptrCompareFunc  pCompFunc;
};

void  heap_init( heap_t  * pHeap, ptrCompareFunc  _pCompFunc )
{
    assert( pHeap != NULL );
    assert( _pCompFunc != NULL );
    
    memset( pHeap, 0, sizeof( heap_t ) );

    pHeap->pCompFunc = _pCompFunc;
    pHeap->max_size = 100;
    pHeap->pArr = (voidPtr_t *)malloc( sizeof( void * ) 
                                       * pHeap->max_size );
    assert( pHeap->pArr != NULL );
    pHeap->curr_size = 0;
}


static void  resize( heap_t  * pHeap, int  size ) 
{
    int  max_sz;
    voidPtr_t  * pTmp;

    if   ( size <= pHeap->max_size ) 
        return;

    max_sz = size * 2;
    pTmp = (voidPtr_t *)malloc( max_sz * sizeof( void * ) );
    assert( pTmp != NULL );
    memset( pTmp, 0, max_sz * sizeof( void * ) );
    memcpy( pTmp, pHeap->pArr, pHeap->curr_size * sizeof( void * ) );
    free( pHeap->pArr );
    pHeap->pArr = pTmp;
}


void  heap_term( heap_t  * pHeap )
{
    assert( pHeap != NULL );

    if  ( pHeap->pArr != NULL ) 
        free( pHeap->pArr );
    memset( pHeap, 0, sizeof( heap_t ) );
}


void  heap_insert( heap_t  * pHeap, void  * pData )
{
    int  ind, father;
    voidPtr_t  pTmp;

    resize( pHeap, pHeap->curr_size + 1 );
    
    ind = pHeap->curr_size;
    pHeap->curr_size++;

    pHeap->pArr[ ind ] = pData;

    while  ( ind > 0 ) {
        father = ( ind - 1 ) / 2;
        
        if  ( (*pHeap->pCompFunc)( pHeap->pArr[ ind ], 
                                   pHeap->pArr[ father ] ) <= 0 )
            break;

        pTmp = pHeap->pArr[ ind ];
        pHeap->pArr[ ind ] = pHeap->pArr[ father ];
        pHeap->pArr[ father ] = pTmp;

        ind = father;
    }    
}


void  * heap_top( heap_t  * pHeap ) 
{
    return  pHeap->pArr[ 0 ];
}


void  * heap_delete_max( heap_t  * pHeap )
{
    void  * res;
    void  * pTmp;
    int  ind, left, right, max_ind;

    if  ( pHeap->curr_size <= 0 )
        return  NULL;

    res = pHeap->pArr[ 0 ];

    //printf( "pHeap->curr_size: %d\n", pHeap->curr_size );
    //printf( "res= %p\n", res );

    pHeap->curr_size--;
    pHeap->pArr[ 0 ] = pHeap->pArr[ pHeap->curr_size ]; 
    pHeap->pArr[ pHeap->curr_size ] = NULL;
    ind = 0;

    while  ( ind < pHeap->curr_size ) {
        left = 2 * ind + 1;
        right = 2 * ind + 2;
        if  ( left >= pHeap->curr_size ) 
            break;
        if  ( right >= pHeap->curr_size ) 
            right = left;
        
        if  ( (*pHeap->pCompFunc)( pHeap->pArr[ left ], 
                                   pHeap->pArr[ right ] ) <= 0 )
            max_ind = right;
        else
            max_ind = left;
        if  ( (*pHeap->pCompFunc)( pHeap->pArr[ ind  ], 
                                   pHeap->pArr[ max_ind ] ) >= 0 )
            break;
        
        pTmp = pHeap->pArr[ ind ];
        pHeap->pArr[ ind ] = pHeap->pArr[ max_ind ];
        pHeap->pArr[ max_ind ] = pTmp;

        ind = max_ind;
    }

    return  res;
}


bool  heap_is_empty( heap_t  * pHeap )
{
    assert( pHeap != NULL );

    return( pHeap->curr_size <= 0 );
}


int   compare_pairs( void  * _a, void  * _b )
{
    const GFSPPair  & a( *(GFSPPair *)_a );
    const GFSPPair  & b( *(GFSPPair *)_b );

    if  ( a.maxDiam() < b.maxDiam() )
        return  -1;
    if  ( a.maxDiam() > b.maxDiam() )
        return  1;

    return  0;
}



class  g_heap_pairs_p
{
private:
    heap_t  heap;

public:
    g_heap_pairs_p() {
        heap_init( &heap, compare_pairs );
    }

    ~g_heap_pairs_p() {
        while  ( size() > 0 ) {
            pop();
        }
        heap_term( &heap );
    }

    void  push( GFSPPair  & pair ) 
    {
        //printf( "pushing: (%p, %p)\n", pair.get_left(),
        //        pair.get_right() );
        GFSPPair  * p_pair = new GFSPPair( pair );
        heap_insert( &heap, (void *)p_pair );
    }

    int  size() const { 
        return  heap.curr_size;
    }
    GFSPPair  top() {
        return  *(GFSPPair *)heap_top( &heap );
    }
    void  pop() {
        GFSPPair  * ptr = (GFSPPair *)heap_delete_max( &heap );

        delete  ptr;
    }
};


/* heap.implementation end */


void  computeExtremePair( const gdiam_point  * arr, int  size,
                          int  dim, GPointPair  & pp )
{
    pp.init( arr[ 0 ] );

    for  ( int  ind = 1; ind < size; ind++ ) {
        const gdiam_point  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  GTreeDiamAlg::compute_by_heap( double  eps ) 
{
    g_heap_pairs_p  heap;
    int  heap_limit, heap_delta;
    
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
    
    GFSPTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), points_num,
                        root->getBBox().getLongestDim(),
                        pair_diam );
    
    GFSPPair  root_pair;
    
    root_pair.init( root, root );
    
    heap.push( root_pair );
    //queue_a.pushx( root_pair );
    
    int  count =0;
    while  ( heap.size() > 0 ) {
        //printf( "heap.size: %d\n", heap.size() );
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        GFSPPair  pair = heap.top();
        heap.pop();
        split_pair( pair, heap, eps );
        //if  ( ( count & 0x3ff ) == 0 ) {
        //    printf( "heap.size: %d\n", heap.size() );
        //}
        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++;
    }
}


void  GTreeDiamAlg::compute_by_heap_proj( double  eps,
                                          gdiam_point  proj ) 
{
    g_heap_pairs_p  heap;
    int  heap_limit, heap_delta;
    GPointPair   pair_diam_x;
    

    //pnt_dump( proj );
    //printf( "length = %g\n", pnt_length( proj ) );

    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
    
    GFSPTreeNode  * root = tree.getRoot();
    
    computeExtremePair( tree.getPoints(), points_num,
                        root->getBBox().getLongestDim(),
                        pair_diam_x );
    pair_diam.init( pair_diam_x.p, pair_diam_x.q, proj );

    GFSPPair  root_pair;
    
    root_pair.init( root, root, proj, 0 );
    
    //printf( "root pair distance: %g\n", root_pair.maxDiam() );
    heap.push( root_pair );
    //queue_a.pushx( root_pair );
    
    int  count =0;
    while  ( heap.size() > 0 ) {
        //printf( "heap.size: %d\n", heap.size() );
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        GFSPPair  pair = heap.top();
        heap.pop();
        //printf( "pair distance: %g\n", pair.maxDiam() );
        split_pair_proj( pair, heap, eps, proj );
        //if  ( ( count & 0x3ff ) == 0 ) {
        //    printf( "heap.size: %d\n", heap.size() );
        //}
        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++;
        //printf( "Poping!\n" );
    }
}


void  GTreeDiamAlg::addPairHeap( g_heap_pairs_p  & heap,
                                 GFSPTreeNode  * left,
                                 GFSPTreeNode  * right )
{
    if  ( update_diam_mode == UDM_SIMPLE ) {
        const gdiam_point  p( *(left->ptr_pnt_left()) );
        const gdiam_point  q( *(right->ptr_pnt_left()) );
        
        pair_diam.update_diam( p, q );
    } else
        if  ( update_diam_mode == UDM_BIG_SAMPLE ) {
            if  ( ( left->size() < 100 )  ||  ( right->size() < 100 ) ) {
                const gdiam_point  p( *(left->ptr_pnt_left()) );
                const gdiam_point  q( *(right->ptr_pnt_left()) );
                
                pair_diam.update_diam( p, q );
                } else
            {                
#define  UMD_SIZE        5
                gdiam_point  arr_left[ UMD_SIZE ], 
                    arr_right[ UMD_SIZE ];
                for  ( int  ind = 0; ind < UMD_SIZE; ind++ ) {                
                    const gdiam_point  p( *(left->ptr_pnt_rand()) );
                    const gdiam_point  q( *(right->ptr_pnt_rand()) );
                    arr_left[ ind ] = p;
                    arr_right[ ind ] = q;

                    //pair_diam.update_diam( p, q );
                }
                for  ( int  ind = 0; ind < UMD_SIZE; ind++ )
                    for  ( int  jnd = 0; jnd < UMD_SIZE; jnd++ )
                        pair_diam.update_diam( arr_left[ ind ], 
                                               arr_right[ jnd ] );
            }
        } else {
            assert( false );
        }
    //printf( "old_diam; %g\n", diam );
    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
    
    GFSPPair  pair;
    
    pair.init( left, right );
    if  ( pair.maxDiam() <= pair_diam.distance )
        return;
    heap.push( pair );
}


void  GTreeDiamAlg::addPairHeap( g_heap_pairs_p  & heap,
                                 GFSPTreeNode  * left,
                                 GFSPTreeNode  * right,
                                 gdiam_point  proj,
                                 GFSPPair  & father )
{
    const gdiam_point  p( *(left->ptr_pnt_left()) );
    const gdiam_point  q( *(right->ptr_pnt_left()) );