ftreeptr.h

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * ftreeptr.h -
 *     My implementation of fair split tree using pointers.
 * 
 * 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  __FTREEPTR__H
#define  __FTREEPTR__H


class FSPPair;
typedef  const Point3d  * CPoint3dPtr;

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

class   FSPTreeNode {
private:
    BBox3d    bbox;
    //int  left_ind, p_pnt_right; // range in the array
    const Point3d   ** p_pnt_left, ** p_pnt_right;
    real  bbox_diam;
    FSPTreeNode  * left, * right;
    Point3d  center;

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

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

        
    FSPTreeNode( CPoint3dPtr   * _p_pnt_left, 
                 CPoint3dPtr   * _p_pnt_right ) {
        bbox.init();
        left = right = NULL;
        p_pnt_left = _p_pnt_left;
        p_pnt_right = _p_pnt_right;
        bbox_diam = -1;
    }

    real  maxDiam() const {
        return  bbox_diam;
    }

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

    int  size() const {
        return  p_pnt_right - p_pnt_left + 1;
    }

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

    const Point3d   ** ptr_pnt_left() {
        return  p_pnt_left;
    }
    const Point3d   ** ptr_pnt_right() {
        return  p_pnt_right;
    }


    friend class  FSPTree;

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


class  FSPTree { 
private:
    Point3d  * pnt_arr;
    CPoint3dPtr  * arr;
    FSPTreeNode  * root;

public:
    void  init( const ArrPoint3d  & _arr ) {
        pnt_arr = (Point3d *)malloc( sizeof( Point3d ) * _arr.size() );
        assert( pnt_arr != NULL );

        arr = (CPoint3dPtr *)malloc( sizeof( CPoint3dPtr ) * _arr.size() );
        assert( arr != NULL );
        //arr.init( _arr.size() + 10, 4353443 );

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

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

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

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

        delete node;
    }
    void  term() {
        free( pnt_arr);
        free( arr );
        pnt_arr = NULL;
        arr = NULL;

        FSPTree::terminate( root );
        root = NULL;
    }

    const Point3d   * getPoints() const {
        return  pnt_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;
    }
    real   brute_diameter( CPoint3dPtr  * a_lo, 
                           CPoint3dPtr  * a_hi,
                           CPoint3dPtr  * b_lo, 
                           CPoint3dPtr  * b_hi,
                           PointPair  & diam ) const { 
        double  max_dist;

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

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

    FSPTreeNode  * getRoot() {
        return  root;
    }

    FSPTreeNode  * build_node( CPoint3dPtr  * left,
                               CPoint3dPtr  * right ) {
        if  ( left > right ) {
            printf( "what!?\n" );
            fflush( stdout );
            assert( left <= right );
        }
        while  ( ( right > left ) 
                 &&  ( (*right)->isEqual( **left ) ) )
            right--;
        FSPTreeNode  * node = new  FSPTreeNode( left, right );

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

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

        return  node;
    }
    
    // return how many elements on left side.
    int    split_array( CPoint3dPtr  * left, 
                        CPoint3dPtr  * right, 
                        int  dim, const real  & val ) {
        CPoint3dPtr  * start_left = left;

        while  ( left < right ) {
            if  ( (**left)[ dim ] < val )
                left++;
            else
                if  ( (**right)[ dim ] >= val )
                    right--;
                else {
                    exchange( *right, *left );
                }
        }

        return  left - start_left;
    }


    void   split_node( FSPTreeNode  * 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->p_pnt_left, node->p_pnt_right, 
                              dim, split_val );
        if  ( left_size <= 0.0 ) {
            printf( "bb: %g   %g\n",
                    bb.min_coord( dim ), bb.max_coord( dim ) );
            printf( "left: %p, right: %p\n",
                    node->p_pnt_left, node->p_pnt_right );
            assert( left_size > 0 );
        }
        if  ( left_size >= (node->p_pnt_right - node->p_pnt_left + 1 ) ) {
            printf( "left size too large?\n" );
            fflush( stdout );
            assert( left_size < (node->p_pnt_right - node->p_pnt_left + 1 ) );
        }

        node->left = build_node( node->p_pnt_left, 
                                 node->p_pnt_left + left_size - 1 );
        node->right = build_node( node->p_pnt_left + left_size, 
                                  node->p_pnt_right );
    }

    void  findExtremeInProjection( FSPPair  & pair, 
                                   PointPair  & max_pair_diam );
};

class  FSPPair 
{
private:
    FSPTreeNode  * left, * right;
    double  max_diam;

public:
    void  init( FSPTreeNode  * _left, FSPTreeNode  * _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;
    }

    FSPTreeNode  * get_left() {
        return  left;
        }
    FSPTreeNode  * 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 FSPTree  & tree,
                             PointPair  & diam ) const {
        return  tree.brute_diameter( left->ptr_pnt_left(),
                                     left->ptr_pnt_right(), 
                                     right->ptr_pnt_left(),
                                     right->ptr_pnt_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   /* __FTREEPTR__H */
#error  Header file ftreeptr.h included twice
#endif  /* __FTREEPTR__H */

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