treediam.c

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * TreeDiam.C -
 *     Implementation of algorithms for computing the diameter 
 * using tree heirarchies.
 * 
 * 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)
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/

#include  <iostream>
#include  <assert.h>
#include  <stdlib.h>
#include  <stdio.h>
#include  <time.h>
#include  <math.h>

#include  <queue>


//#include  "dmalloc.h"

#include  "globjwin.h"

#include  "_generic.h"
#include  "sarray.h"
#include  "real.h"
#include  "point.h"
#include  "bbox.h"
#include  "ftree.h"
#include  "ftree2d.h"
#include  "ftreeptr.h"
#include  "TreeDiam.h"
#include  "mat3d.h"
#include  "misc.h"
//#include  "terrain.h"


#define  HEAP_LIMIT  100000
#define  HEAP_DELTA  100000

//class  heap_pairs;
//class  FSPairArray;
struct lt_fspair
{
  bool operator()(const FSPair  & a , const FSPair  & b ) const
  {
      bool  a_iis, b_iis;

      if  ( ( ! a.isIdenticalSides() )  
            &&  ( ! b.isIdenticalSides() ) )
          return   a.maxDiam() < b.maxDiam();
      
      a_iis = a.isIdenticalSides();
      b_iis = b.isIdenticalSides();
      if  ( a_iis  &&  b_iis )
          return   a.maxDiam() < b.maxDiam();
      if  ( a_iis  &&  (! b_iis ) )
          return  false;
      else
          return  true;
      //return   a.minAprxDiam() < b.minAprxDiam();
  }
};
struct lt_fspair_wspd
{
  bool operator()(const FSPair  & a , const FSPair  & b ) const
  {
      return   a.maxError() < b.maxError();
      //return   a.minAprxDiam() < b.minAprxDiam();
  }
};
struct lt_fsppair
{
    bool operator()(const FSPPair  & a , const FSPPair  & b ) const
    {
        return   a.maxDiam() < b.maxDiam();
        //return   a.minAprxDiam() < b.minAprxDiam();
    }
};
struct lt_fspair_2d
{
  bool operator()(const FS2dPair  & a , const FS2dPair  & b ) const
  {
      return   a.maxDiam() < b.maxDiam();
      //return   a.minAprxDiam() < b.minAprxDiam();
  }
};

/*--- Constants ---*/
#define   THRESHOLD_LOCAL   400
#define  THRESHOLD   10

class  heap_pairs : public std::priority_queue<FSPair, 
                    std::vector<FSPair>, lt_fspair>
{
};

//   heap_pairs;
class  FSPairArray : public  Array<FSPair> {};

/*--- Start of Code ---*/

real   brute_diameter_range( const ArrPoint3d  & arr,
                             int  start,
                             int  end,
                             PointPair  & diam )
{
    for  ( int  ind = start; ind < end; ind++ )
        for  ( int  jnd = ind + 1; jnd <= end; jnd++ ) 
            diam.update_diam( arr[ ind ], arr[ jnd ] );

    return  diam.distance;
}


#define   THRESHOLD_LOCAL   400

real   brute_diameter_local_bi( const ArrPoint3d  & arr,
                                int  a_start,
                                int  a_end,
                                int  b_start,
                                int  b_end,
                                PointPair  & diam )
{
    int   a_mid, b_mid;

    if  ( ( a_end - a_start ) < THRESHOLD_LOCAL ) {
        for  ( int  ind = a_start; ind <= a_end; ind++ )
            for  ( int  jnd = b_start; jnd <= b_end; jnd++ ) 
                diam.update_diam( arr[ ind ], arr[ jnd ] );

        return  diam.distance;
    }

    a_mid = (a_start + a_end) / 2;
    b_mid = (b_start + b_end) / 2;
    brute_diameter_local_bi( arr, 
                             a_start, a_mid, 
                             b_start, b_mid, diam );
    brute_diameter_local_bi( arr,
                             a_start, a_mid, 
                             b_mid + 1, b_end, diam );
    brute_diameter_local_bi( arr,
                             a_mid + 1, a_end, 
                             b_start, b_mid, diam );
    brute_diameter_local_bi( arr,
                             a_mid + 1, a_end,
                             b_mid + 1, b_end, diam );

    return  diam.distance;
}



real   brute_diameter_local( const ArrPoint3d  & arr,
                             int  start,
                             int  end,
                             PointPair  & out )
{
    int  mid;

    if  ( ( end - start ) < THRESHOLD_LOCAL )
        return   brute_diameter_range( arr, start, end, out );
        
    mid = ( start + end ) / 2;

    brute_diameter_local_bi( arr, start, mid, mid + 1, end, out );
    brute_diameter_local( arr, start, mid, out );
    brute_diameter_local( arr, mid + 1, end, out );

    return  out.distance;
}


real   brute_diameter_local( const ArrPoint3d  & arr,
                             PointPair  & out )
{
    out.p = out.q = arr[ 0 ];
    out.distance = 0;

    return  brute_diameter_local( arr, 0, arr.size() - 1, out );
}



void   fill_array_cube( ArrPoint3d  & arr, int  num )
{
    arr.init( num + 10, 43897543 );

    for  ( int  ind = 0; ind < num; ind++ ) {
        Point3d  pnt;
        
        pnt.fillRandom();
        arr.pushx( pnt );
    }

}


void   fill_array_box( ArrPoint3d  & arr, int  num,
                       double  a, double  b, double  c )
{
    arr.init( num + 10, 43897543 );

    for  ( int  ind = 0; ind < num; ind++ ) {
        Point3d  pnt;
        
        pnt.fillRandom();
        Point3d  pntx( pnt[ 0 ] * a, pnt[ 1 ] * b, pnt[ 2 ] * c );

        arr.pushx( pntx );
    }

}


void   fill_array_sphere( ArrPoint3d  & arr, int  num )
{
    arr.init( num + 10, 43897543 );

    for  ( int  ind = 0; ind < num; ind++ ) {
        double  alpha, beta;
        Point3d  pnt;
        
        alpha = realRand() * 2.0 * 3.1415926535;
        beta = realRand() * 2.0 * 3.1415926535;

        pnt.setCoord( 0, cos( alpha ) * sin( beta ) );
        pnt.setCoord( 1, cos( alpha ) * cos( beta ) );
        pnt.setCoord( 2, sin( alpha ) );
        
        //pnt.fillRandom();
        arr.pushx( pnt );
    }

}


void   fill_array_sphere_caps( ArrPoint3d  & arr, int  num,
                            double  eps )
{
    arr.init( num + 10, 43897543 );

    for  ( int  ind = 0; ind < num; ind++ ) {
        double  alpha, beta;
        Point3d  pnt;
        
        alpha = (1.0 - eps * realRand()) * 3.1415926535 / 2.0;
        if  ( realRand() >= 0.5 ) 
            alpha = -alpha;

        beta = realRand() * 2.0 * 3.1415926535;

        pnt.setCoord( 0, cos( alpha ) * sin( beta ) );
        pnt.setCoord( 1, cos( alpha ) * cos( beta ) );
        pnt.setCoord( 2, sin( alpha ) );
        
        //pnt.fillRandom();
        arr.pushx( pnt );
    }

}


void   fill_array_two_arcs( ArrPoint3d  & arr, int  num )
{
    double  alpha_limit, alpha_increment, alpha;
    Point3d  pnt;

    arr.init( num + 10, 43897543 );

    num  = num / 2;

    alpha_limit = 1.0 /((double)num*num);
    alpha_increment = 2.0 * alpha_limit  / (double)num;
    
    Transform3d  trans_x, trans_z;
    trans_x.init_rot_around_x( 3.14150265358979323 / 4.0 );
    trans_z.init_rot_around_z( 3.14150265358979323 / 4.0 );

    alpha = -alpha_limit;
    for  ( int  ind = 0; ind < num; ind++ ) {
        pnt = Point3d( -cos( alpha ), 0, sin( alpha ) );
        trans_x.transform( pnt );
        trans_z.transform( pnt );
        arr.pushx( pnt );

        pnt = Point3d( cos( alpha ), sin( alpha ), 0 );
        trans_x.transform( pnt );
        trans_z.transform( pnt );
        arr.pushx( pnt );

        alpha += alpha_increment;
    }

}

void   fill_array_mutated_ellipsoid( ArrPoint3d  & arr, int  num,
                                     double  a, double  b, double  c )
{
    Point3d  new_pnt, pnt;

    arr.init( num + 10, 43897543 );

    Point3d  vec1( 1, 1, 1 );
    Point3d  vec2( 1, -1, 0 );
    Point3d  vec3;

    vec3 = cross_prod( vec1, vec2 );

    vec1.normalize();
    vec2.normalize();
    vec3.normalize();

    int  limit;

    limit = num / 11;
    for  ( int  ind = 0; ind < limit; ind++ ) {
        double  alpha, beta;
        Point3d  pnt;
        
        alpha = realRand() * 2.0 * 3.1415926535;
        beta = realRand() * 2.0 * 3.1415926535;

        pnt.setCoord( 0, a * cos( alpha ) * sin( beta ) );
        pnt.setCoord( 1, b * cos( alpha ) * cos( beta ) );
        pnt.setCoord( 2, c * sin( alpha ) );
        
        new_pnt = changeBase( vec1, vec2, vec3, pnt );
        
        //pnt.fillRandom();
        arr.pushx( new_pnt );
    }

    limit = 5 * num / 11;
    for  ( int  ind = 0; ind < limit; ind++ ) {
        pnt.setCoord( 0, 0 );
        pnt.setCoord( 1, 0 );
        pnt.setCoord( 2, 0.999 * c );

        new_pnt = changeBase( vec1, vec2, vec3, pnt );        
        arr.pushx( new_pnt );
    }
    for  ( int  ind = 0; ind < limit; ind++ ) {
        pnt.setCoord( 0, 0 );
        pnt.setCoord( 1, 0.999 * b );
        pnt.setCoord( 2, 0 );

        new_pnt = changeBase( vec1, vec2, vec3, pnt );        
        arr.pushx( new_pnt );
    }
}


static void   addPoint( ArrPoint3d  & arr, double  x, double  y, double  z )
{
    Point3d  pnt( x, y, z );
    Point3d  pnt1;

    pnt1.fillRandom();
    pnt.add_scale( pnt1, 0.00001 );
    arr.pushx( pnt );
}


// Cheats the regular bbox/pca algorithm. so they realize 1.72
// approximation in 3d
void   fill_array_strange_cube( ArrPoint3d  & arr, int  num )
{
    int limit;
    //Point3d  pnt, pnt1;
        
    arr.init( num + 10, 43897543 );

    addPoint( arr, -1, -1, -1 );
    addPoint( arr, -1, -1,  1 );
    addPoint( arr, -1,  1, -1 );
    addPoint( arr, -1,  1,  1 );
    addPoint( arr,  1, -1, -1 );
    addPoint( arr,  1, -1,  1 );
    addPoint( arr,  1,  1, -1 );
    addPoint( arr,  1,  1,  1 );

    addPoint( arr,  0,  0,  1.01000001 );
    addPoint( arr,  0,  0, -1.0100001 );
    addPoint( arr,  0,  1.01, 0 );
    addPoint( arr,  0, -1.01, 0 );
    addPoint( arr,  1.01, 0, 0 );
    addPoint( arr, -1.01, 0, 0 );

    limit = num / 12;
    for  ( int  ind = 0; ind < limit; ind++ ) {
        addPoint( arr,  0.001,  0,  1.011 );
        addPoint( arr,  0,  0.0001, -1.011 );
        addPoint( arr,  0.001,  0,  1.011 );
        addPoint( arr,  0,  0.0001, -1.011 );
        addPoint( arr,  0.001,  0,  1.011 );
        addPoint( arr,  0,  0.0001, -1.011 );

        addPoint( arr,  0,  1.011, 0 );
        addPoint( arr,  0, -1.011, 0 );
        addPoint( arr,  0,  1.011, 0 );
        addPoint( arr,  0, -1.011, 0 );

        addPoint( arr,   1.011, 0, 0 );
        addPoint( arr,  -1.011, 0, 0 );
    }
}



void   fill_array_strange_cube_diag( ArrPoint3d  & arr, int  num )
{
    int limit;
    //Point3d  pnt, pnt1;
        
    arr.init( num + 10, 43897543 );

    addPoint( arr, -1, -1, -1 );
    addPoint( arr, -1, -1,  1 );
    addPoint( arr, -1,  1, -1 );
    addPoint( arr, -1,  1,  1 );
    addPoint( arr,  1, -1, -1 );
    addPoint( arr,  1, -1,  1 );
    addPoint( arr,  1,  1, -1 );
    addPoint( arr,  1,  1,  1 );

    addPoint( arr,  0,  0,  1.1000001 );
    addPoint( arr,  0,  0, -1.100001 );
    addPoint( arr,  0,  1.1, 0 );
    addPoint( arr,  0, -1.1, 0 );
    addPoint( arr,  1.1, 0, 0 );
    addPoint( arr, -1.1, 0, 0 );