ftree.h

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * ftree.h -
 *     My implementation of fair split tree.
 * 
 * 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  __FTREE__H
#define  __FTREE__H


class FSPair;

void  computeExtremePairDir( const ArrPoint3d  & arr,
                             const Point3d  & dir,
                             PointPair  & pp,
                             int  start,
                             int  end
                             );

class   FSTreeNode {
private:
    BBox3d    bbox;
    int  left_ind, right_ind; // range in the array
    real  bbox_diam;
    FSTreeNode  * left, * right;
    Point3d  center;

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

    const Point3d  & 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;
    }

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

    real  maxDiam() const {
        return  bbox_diam;
    }

    FSTreeNode  * get_left() {
        return  left;
        }
    FSTreeNode  * 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  FSTree;

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


class  FSTree { 
private:
    ArrPoint3d  arr;
    FSTreeNode  * root;

public:
    void  init( const ArrPoint3d  & _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( FSTreeNode  * node ) {
        if  ( node == NULL )
            return;
        terminate( node->left );
        terminate( node->right );

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

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

    const  ArrPoint3d   & 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,
                           PointPair  & 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 Point3d  & getPoint( int  ind ) const {
        return  arr[ ind ];
    }

    FSTreeNode  * getRoot() {
        return  root;
    }

    FSTreeNode  * build_node( int   left, int   right ) {
        if  ( left > right ) {
            printf( "what!?\n" );
            fflush( stdout );
            assert( left <= right );
        }
        while  ( ( right > left ) 
                 &&  ( arr[ right ].isEqual( arr[ left ] ) ) )
            right--;
        FSTreeNode  * node = new  FSTreeNode( 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( FSTreeNode  * node ) 
    {
        int  dim, left_size;
        real  split_val;

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

        BBox3d  & 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 );
        if  ( left_size <= 0.0 ) {
            printf( "bb: %g   %g\n",
                    bb.min_coord( dim ), bb.max_coord( dim ) );
            printf( "left: %d, right: %d\n",
                    node->left_ind, node->right_ind );
            for  ( int  ind = node->left_ind; ind <= node->right_ind; ind++ ) {
                printf( "ind: %d: \n", ind );
                arr[ ind ].print();
                printf( "\n" );
            }
            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( FSPair  & pair, 
                                   PointPair  & max_pair_diam );
};

class  FSPair 
{
private:
    FSTreeNode  * left, * right;
    double  max_diam;

public:
    bool  isIdenticalSides() const {
        return  left == right;
    }

    void  setLeftBigger() {
        if  ( left->maxDiam() < right->maxDiam() )
            exchange( left, right );
    }

    void  init( FSTreeNode  * _left, FSTreeNode  * _right )
    {
        left = _left;
        right = _right;

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

    FSTreeNode  * get_left() {
        return  left;
        }
    FSTreeNode  * 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 FSTree  & tree, PointPair  & diam ) const {
        return  tree.brute_diameter( left->index_left(),
                                     left->index_right(),
                                     right->index_left(),
                                     right->index_right(),
                                     diam );
    }

    real  maxError() const {
        return  max( left->maxDiam(), right->maxDiam() );
    }

    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   /* __FTREE__H */
#error  Header file ftree.h included twice
#endif  /* __FTREE__H */

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