ftree2d.h

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * ftree2d.h -
 * 
 * Copyright 2000 Sariel Har-Peled ([e-mail])
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation; either version 2, or (at your option) any
 * later version.
 * 
 * Code is based on the paper:
 *   A Practical Approach for Computing the Diameter of a Point-Set.
 *   Sariel Har-Peled (http://www.uiuc.edu/~sariel)
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/

#ifndef  __FTREE2D__H
#define  __FTREE2D__H

class  FS2dPair;
void  computeExtremePairDir( const ArrPoint2d  & arr,
                             const Point2d  & dir,
                             PointPair2d  & pp,
                             int  start,
                             int  end
                             );

class   FS2dTreeNode {
private:
    BBox2d    bbox;
    int  left_ind, right_ind; // range in the array
    real  bbox_diam;
    FS2dTreeNode  * left, * right;
    Point2d  center;

    friend class  FS2dTree;
public:
    void  dump() const {
        printf( "\tNode: Range: [%d, %d]\n", left_ind, right_ind );
    }

    const Point2d  & getCenter() const {
        return  center;
    }
    int  nodes_number() const {
        int  sum;

        if  ( ( left == NULL )  &&  ( right == NULL ) )
            return  1;
        sum = 1;

        if  ( left != NULL )
            sum += left->nodes_number();
        if  ( right != NULL )
            sum += right->nodes_number();

        return  sum;
    }

        
    FS2dTreeNode( int _left, int  _right ) {
        bbox.init();
        left = right = NULL;
        left_ind = _left;
        right_ind = _right;
        bbox_diam = -1;
    }

    real  maxDiam() const {
        return  bbox_diam;
    }

    FS2dTreeNode  * get_left() {
        return  left;
        }
    FS2dTreeNode  * get_right() {
        return   right;
    }

    int  size() const {
        return  right_ind - left_ind + 1;
    }

    bool  isSplitable() const {
        return  size() > 1;
    }

    int  index_left() {
        return  left_ind;
    }
    int  index_right() {
        return  right_ind;
    }


    //friend class  FS2dTree;

    const BBox2d  & getBBox() const {
        return  bbox;
    }
};


class  FS2dTree { 
private:
    ArrPoint2dPtr  arr;
    FS2dTreeNode  * root;

public:
    void  init( const ArrPoint2dPtr  & _arr ) {
        arr.init( _arr.size() + 10, 4353443 );

        for  ( int  ind = 0; ind < _arr.size(); ind++ )
            arr.pushx( _arr[ ind ] );

        root = build_node( 0, arr.size() - 1 );
    }

    static void  terminate( FS2dTreeNode  * node ) {
        if  ( node == NULL )
            return;
        terminate( node->left );
        terminate( node->right );

        node->left = node->right = NULL;

        delete node;
    }
    void  term() {
        arr.term();
        FS2dTree::terminate( root );
        root = NULL;
    }


    const  ArrPoint2dPtr   & getPoints() const {
        return  arr;
    }

    int  nodes_number() const {
        return  root->nodes_number();
    }

    real   brute_diameter( int  a_lo, int  a_hi,
                           int  b_lo, int  b_hi,
                           PointPair2d  & diam ) const { 
        double  max_dist;

        max_dist = 0;
        for  ( int  ind = a_lo; ind <= a_hi; ind++ )
            for  ( int  jnd = b_lo; jnd <= b_hi; jnd++ )
                diam.update_diam( *(arr[ ind ]), *(arr[ jnd ]) );
        
        return  diam.distance;
    }

    const Point2d  & getPoint( int  ind ) const {
        return  *(arr[ ind ]);
    }

    FS2dTreeNode  * getRoot() {
        return  root;
    }

    FS2dTreeNode  * build_node( int   left, int   right ) {
        if  ( left > right ) {
            printf( "what!?\n" );
            fflush( stdout );
            assert( left <= right );
        }
        FS2dTreeNode  * node = new  FS2dTreeNode( left, right );

        node->bbox.init();
        for  ( int  ind = left; ind <= right; ind++ ) 
            node->bbox.bound( *(arr[ ind ]) );
        node->bbox_diam = node->bbox.get_diam();

        node->center = node->bbox.center();

        return  node;
    }
    
    // rturn how many elements on left side.
    int    split_array( int  left, int  right, int  dim, const real  & val ) {
        int  start_left = left;
        while  ( left < right ) {
            if  ( (*(arr[ left ]))[ dim ] < val ) 
                left++;
            else
                if  ( ( *(arr[ right ]) )[ dim ] >= val ) 
                    right--;
                else {
                    exchange( arr[ right ], arr[ left ] );
                }
        }

        return  left - start_left;
    }


    void   split_node( FS2dTreeNode  * node ) 
    {
        int  dim, left_size;
        real  split_val;

        if  ( node->left != NULL )
            return;

        BBox2d  & bb( node->bbox );

        dim = bb.getLongestDim();
        split_val = ( bb.min_coord( dim ) + bb.max_coord( dim ) ) / 2.0;

        left_size = split_array( node->left_ind, node->right_ind, 
                              dim, split_val );
        assert( left_size > 0 );
        if  ( left_size >= (node->right_ind - node->left_ind + 1 ) ) {
            printf( "left size too large?\n" );
            fflush( stdout );
            assert( left_size < (node->right_ind - node->left_ind + 1 ) );
        }

        node->left = build_node( node->left_ind, 
                                 node->left_ind + left_size - 1 );
        node->right = build_node( node->left_ind + left_size, 
                                  node->right_ind );
    }

    void  findExtremeInProjection( FS2dPair  & pair, 
                                   PointPair2d  & max_pair_diam );
};

class  FS2dPair 
{
private:
    FS2dTreeNode  * left, * right;
    double  max_diam;

public:
    void  init( FS2dTreeNode  * _left, FS2dTreeNode  * _right )
    {
        left = _left;
        right = _right;

        BBox2d  bb;

        bb.init( left->getBBox(), right->getBBox() );
        //printf( "\t\tbbox: \n" );
        //bb.dump();
        max_diam = bb.get_diam();
        //printf( "\t\tmax_diam: %g\n", max_diam );
    }

    // error 2r/l - approximate cos of the max possible angle in the pair
    double  getProjectionError() const 
    {
        double  l, two_r;

        l = pnt_distance( left->getCenter(), right->getCenter() );
        two_r = max( left->maxDiam(), right->maxDiam() );

        if  ( l == 0.0 )
            return  10;

        return  two_r / l;
    }

    FS2dTreeNode  * get_left() {
        return  left;
        }
    FS2dTreeNode  * get_right() {
        return   right;
    }

    void  dump() const {
        printf( "[\n" );
        left->dump();
        right->dump();
        printf( "\tmax_diam; %g\n", max_diam );
    }

    bool  isBelowThreshold( int  threshold ) const {
        if  ( left->size() > threshold )
            return  false;
        if  ( right->size() > threshold )
            return  false;

        return  true;
    }

    double   directDiameter( const FS2dTree  & tree, 
                             PointPair2d  & diam ) const {
        return  tree.brute_diameter( left->index_left(),
                                     left->index_right(),
                                     right->index_left(),
                                     right->index_right(),
                                     diam );
    }

    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;
    }
};

#else   /* __FTREE2D__H */
#error  Header file ftree2d.h included twice
#endif  /* __FTREE2D__H */

/* ftree2d.h - End of File ------------------------------------------*/