algorithm.c

任意给定三维空间的点集

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * algorithm.C -
 *     Algorithms for computing the diameter
 * 
 * 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  <qapplication.h>
#include  <qgl.h>

//#include  "dmalloc.h"

//#include  "scape/vector3d.h"
//#include  "line.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  "proj_points.h"
#include  "misc.h"
#include  "obb.h"

#include  "algorithm.h"
#include  "gdiam.h"


/*--- Constants ---*/

//------------------------------------------------------------------
// DiameterAlgorithm - Generic class for algorithms computing diam
//------------------------------------------------------------------
void   DiameterAlgorithm::dump()
{
    printf( "%-25s  %10.3f   %10.3g  eps: %10.4g\n",                   
            getID(),
            time,
            diam.distance,
            eps );
}


void   DiameterAlgorithm::dump_to_file( FILE  * fl )
{
    fprintf( fl, "%10.2f\n",                   
             time );
}


DiameterAlgorithm  * BruteDiameter::CreateAndCompute( ArrPoint3d  & arr )
{
    BruteDiameter  * pDiam = new BruteDiameter;

    pDiam->eps = 0;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


char   * BruteDiameter::getID()
{
    static char  * name = "Brute force diameter";

    return  name;
}

void   BruteDiameter::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    diam.init( arr[ 0 ] );

    tm.start();

    brute_diameter( arr, diam );

    tm.end();
    time = tm.seconds();
}

//------------------------------------------------------------------
// BruteDiameterLocal
//    Compute diameter in a brute force manner - using cache locality
//    to improve speed.
//------------------------------------------------------------------

char   * BruteDiameterLocal::getID()
{
    static char  * name = "Brute force diam local";

    return  name;
}


DiameterAlgorithm  * BruteDiameterLocal::CreateAndCompute( 
                                                 ArrPoint3d  & arr )
{
    BruteDiameterLocal  * pDiam = new BruteDiameterLocal;

    pDiam->eps = 0;
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}

    
void   BruteDiameterLocal::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    diam.init( arr[ 0 ] );
    tm.start();

    brute_diameter_local( arr, diam );

    tm.end();
    time = tm.seconds();    
}


//------------------------------------------------------------------
// QuadDiameterByLevel
//------------------------------------------------------------------
char   * QuadDiameterByLevel::getID()
{
    static char  * name = "FairSplit using levels";

    return  name;
}


void   QuadDiameterByLevel::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    tm.start();

    p_alg = new TreeDiamAlg( arr );
    p_alg->compute_by_level( eps );
    tm.end();

    diam = p_alg->getDiameter();

    time = tm.seconds();   
    p_alg->term();
    delete  p_alg;
    p_alg = NULL;
}


DiameterAlgorithm  * QuadDiameterByLevel::CreateAndCompute( 
                         ArrPoint3d  & arr,
                         double  _eps )
{
    QuadDiameterByLevel   * pDiam = new QuadDiameterByLevel;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}

//------------------------------------------------------------------
// QuadDiameterPtrByLevel
//------------------------------------------------------------------
char   * QuadDiameterPtrByLevel::getID()
{
    static char  * name = "FairSplitPtr using levels";

    return  name;
}


void   QuadDiameterPtrByLevel::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    tm.start();

    //printf( "a\n" );
    p_alg = new TreeDiamPAlg( arr );
    //printf( "b\n" );-
    p_alg->compute_by_level( eps );
    //printf( "c\n" );
    tm.end();

    diam = p_alg->getDiameter();

    time = tm.seconds();    
    p_alg->term();
    delete  p_alg;
    p_alg = NULL;
}

DiameterAlgorithm  * QuadDiameterPtrByLevel::CreateAndCompute( 
                         ArrPoint3d  & arr,
                         double  _eps )
{
    QuadDiameterPtrByLevel   * pDiam = new QuadDiameterPtrByLevel;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


//------------------------------------------------------------------
// QuadDiameterByDir
//------------------------------------------------------------------
char   * QuadDiameterByDir::getID()
{
    static char  * name = "FairSplit by direction";

    return  name;
}


void   QuadDiameterByDir::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    tm.start();

    p_alg = new TreeDiamAlg( arr );
    p_alg->compute_by_dir( eps );
    tm.end();

    diam = p_alg->getDiameter();

    time = tm.seconds();    
    p_alg->term();
    delete  p_alg;
    p_alg = NULL;
}

DiameterAlgorithm  * QuadDiameterByDir::CreateAndCompute( 
                         ArrPoint3d  & arr,
                         double  _eps )
{
    QuadDiameterByDir   * pDiam = new QuadDiameterByDir;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}



//------------------------------------------------------------------
// QuadDiameterByHeap
//------------------------------------------------------------------
char   * QuadDiameterByGPtrHeap::getID()
{
    static char  * name = "FS gdiam ptrheap";

    return  name;
}


void   QuadDiameterByGPtrHeap::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    gdiam_real  * p_arr, * p_start;
    //double  num;

    tm.start();

    p_arr = (gdiam_real *)malloc( sizeof( gdiam_real ) * arr.size() * 3 );
    p_start = p_arr;
    for  ( int  ind = 0; ind < arr.size(); ind++ ) {
        p_arr[ 0 ] = arr[ ind ][ 0 ];
        p_arr[ 1 ] = arr[ ind ][ 1 ];
        p_arr[ 2 ] = arr[ ind ][ 2 ];
        p_arr += 3;
    }

    GPointPair  pair;

    pair = gdiam_approx_diam_pair( p_start, arr.size(), eps );
    //gdiam_approx_const_mvbb( p_start, arr.size(), eps );
    
    diam.init( arr[ (pair.p - p_start) / 3 ],
               arr[ (pair.q - p_start) / 3 ] );

    free( p_start );
    
    tm.end();

    time = tm.seconds();    
}


DiameterAlgorithm  * QuadDiameterByGPtrHeap::CreateAndCompute( 
                             ArrPoint3d  & arr, double  _eps )
{
    QuadDiameterByGPtrHeap  * pDiam;

    pDiam = new QuadDiameterByGPtrHeap;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


//------------------------------------------------------------------
// QuadDiameterByHeapUDMSample
//------------------------------------------------------------------
char   * QuadDiameterByGPtrHeapUDMSample::getID()
{
    static char  * name = "FS(UDM) gdiam ptrheap";

    return  name;
}



void   QuadDiameterByGPtrHeapUDMSample::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    gdiam_real  * p_arr, * p_start;
    //double  num;

    tm.start();

    p_arr = (gdiam_real *)malloc( sizeof( gdiam_real ) * arr.size() * 3 );
    p_start = p_arr;
    for  ( int  ind = 0; ind < arr.size(); ind++ ) {
        p_arr[ 0 ] = arr[ ind ][ 0 ];
        p_arr[ 1 ] = arr[ ind ][ 1 ];
        p_arr[ 2 ] = arr[ ind ][ 2 ];
        p_arr += 3;
    }

    GPointPair  pair;

    pair = gdiam_approx_diam_pair_UDM( p_start, arr.size(), eps );
    //gdiam_approx_const_mvbb( p_start, arr.size(), eps );
    
    diam.init( arr[ (pair.p - p_start) / 3 ],
               arr[ (pair.q - p_start) / 3 ] );

    free( p_start );
    
    tm.end();

    time = tm.seconds();    
}

DiameterAlgorithm  * QuadDiameterByGPtrHeapUDMSample::CreateAndCompute( 
                             ArrPoint3d  & arr, double  _eps )
{
    QuadDiameterByGPtrHeapUDMSample  * pDiam;

    pDiam = new QuadDiameterByGPtrHeapUDMSample;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


//------------------------------------------------------------------
// QuadDiameterByHeap
//------------------------------------------------------------------
char   * QuadDiameterByHeap::getID()
{
    static char  * name = "FairSplit using heap";

    return  name;
}


void   QuadDiameterByHeap::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    tm.start();

    p_alg = new TreeDiamAlg( arr );
    p_alg->compute_by_heap( eps );
    tm.end();

    diam = p_alg->getDiameter();

    time = tm.seconds();    
    p_alg->term();
    delete  p_alg;
    p_alg = NULL;
}

DiameterAlgorithm  * QuadDiameterByHeap::CreateAndCompute( 
                             ArrPoint3d  & arr, double  _eps )
{
    QuadDiameterByHeap  * pDiam = new QuadDiameterByHeap;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}

//------------------------------------------------------------------
// QuadDiameterByHeapWSPD
//------------------------------------------------------------------
char   * QuadDiameterByHeapWSPD::getID()
{
    static char  * name = "FairSplit using heap WSPD";

    return  name;
}


void   QuadDiameterByHeapWSPD::compute( const ArrPoint3d  & arr )
{
    Timer  tm;

    tm.start();


    p_alg = new TreeDiamAlg( arr );
    p_alg->compute_by_heap_wspd( eps );
    tm.end();

    diam = p_alg->getDiameter();

    time = tm.seconds();    
    p_alg->term();
    delete  p_alg;
    p_alg = NULL;
}

DiameterAlgorithm  * QuadDiameterByHeapWSPD::CreateAndCompute( 
                             ArrPoint3d  & arr, double  _eps )
{
    QuadDiameterByHeapWSPD  * pDiam = new QuadDiameterByHeapWSPD;

    pDiam->eps = _eps;

    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}


//------------------------------------------------------------------
// DiameterAPBBox -
//    Compute the largest pair from the bounding box...
//------------------------------------------------------------------
DiameterAlgorithm  * DiameterAPBBox::CreateAndCompute( ArrPoint3d  & arr )
{
    DiameterAPBBox  * pDiam = new DiameterAPBBox;

    pDiam->eps = sqrt(3);
    pDiam->compute( arr );

    return  pDiam;  //(DiameterAlgorithm *)
}

char   * DiameterAPBBox::getID()
{
    static char  * name = "Diameter from BBox";

    return  name;
}


void   DiameterAPBBox::compute( const ArrPoint3d  & arr )
{
    Timer  tm;
    PointPair  pp0, pp1, pp2;
    
    diam.init( arr[ 0 ] );
    tm.start();
    
    computeExtremePairAllDim( arr, pp0, pp1, pp2 );

    diam.update_diam( pp0 );
    diam.update_diam( pp1 );
    diam.update_diam( pp2 );

    tm.end();
 
    time = tm.seconds();    
}