gdiam.cpp

Go to the documentation of this file.

/*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
 * gdiam.C -
 *   Computes diameter, and a tight-fitting bounding box of a
 *   point-set in 3d.
 *
 * Copyright 2000 Sariel Har-Peled (ssaarriieell@cs.uiuc.edu)
 *
 * 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 (in ACM Sym. on Computation Geometry
 *                   SoCG 2001)
 *   Sariel Har-Peled (http://www.uiuc.edu/~sariel)
 *--------------------------------------------------------------
 * History
 * 3/28/01 -
 *     This is a more robust version of the code. It should
 *     handlereally abnoxious inputs well (i.e., points with equal
 *     coordinates, etc.
\*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*/
 
#include  <stdlib.h>
#include  <stdio.h>
#include  <assert.h>
#include  <memory.h>
#include  <math.h>
 
#include  <vector>
#include  <algorithm>
 
#include  "gdiam.h"
 
/*--- Constants ---*/
 
#define  HEAP_LIMIT  10000
#define  HEAP_DELTA  10000
#define  PI 3.14159265958979323846
 
/*--- Start of Code ---*/
 
typedef long double  ldouble;
 
class GFSPPair;
 
/*void  computeExtremePairDir( const ArrPoint3d  & arr,
                             const Point3d  & dir,
                             GPointPair  & pp,
                             int  start,
                             int  end
                             );
*/
 
class   GFSPTreeNode
{
private:
    GBBox    bbox;
    //int  left_ind, p_pnt_right; // range in the array
    gdiam_point*    p_pnt_left, * p_pnt_right;
    gdiam_real  bbox_diam, bbox_diam_proj;
    GFSPTreeNode*   left, * right;
    gdiam_point_t  center;
 
public:
    void  dump() const
    {
        printf("dump...\n");
        //printf( "\tNode: Range: [%d, %d]\n", left_ind, p_pnt_right );
    }
    int  size() const
    {
        return (int)(p_pnt_right - p_pnt_left);
    }
    gdiam_point  getCenter() const
    {
        return (gdiam_point)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;
    }
 
 
    GFSPTreeNode(gdiam_point*    _p_pnt_left,
                 gdiam_point*    _p_pnt_right)
    {
        bbox.init();
        left = right = NULL;
        p_pnt_left = _p_pnt_left;
        p_pnt_right = _p_pnt_right;
        bbox_diam = bbox_diam_proj = -1;
    }
 
    gdiam_real  maxDiam() const
    {
        return  bbox_diam;
    }
    gdiam_real  maxDiamProj() const
    {
        return  bbox_diam_proj;
    }
 
    GFSPTreeNode*   get_left()
    {
        return  left;
    }
    GFSPTreeNode*   get_right()
    {
        return   right;
    }
 
    bool  isSplitable() const
    {
        return (size() > 1);
    }
 
    const gdiam_point*    ptr_pnt_left()
    {
        return  p_pnt_left;
    }
    const gdiam_point*    ptr_pnt_rand()
    {
        return  p_pnt_left + (rand() % (p_pnt_right - p_pnt_left + 1));
    }
    const gdiam_point*    ptr_pnt_right()
    {
        return  p_pnt_right;
    }
 
 
    friend class  GFSPTree;
 
    const GBBox&   getBBox() const
    {
        return  bbox;
    }
};
 
 
class  GFSPTree
{
private:
    gdiam_point*   arr;
    GFSPTreeNode*   root;
 
public:
    void  init(const gdiam_point*   _arr,
               int  size)
    {
        arr = (gdiam_point*)malloc(sizeof(gdiam_point) * size);
        assert(arr != NULL);
 
        for (int  ind = 0; ind < size; ind++)
        {
            arr[ ind ] = _arr[ ind ];
        }
 
        root = build_node(arr, arr + size - 1);
    }
 
    static void  terminate(GFSPTreeNode*   node)
    {
        if (node == NULL)
        {
            return;
        }
 
        terminate(node->left);
        terminate(node->right);
 
        node->left = node->right = NULL;
 
        delete node;
    }
    void  term()
    {
        free(arr);
        arr = NULL;
 
        GFSPTree::terminate(root);
        root = NULL;
    }
 
    const gdiam_point*    getPoints() const
    {
        return  arr;
    }
 
    int  nodes_number() const
    {
        return  root->nodes_number();
    }
 
    gdiam_real   brute_diameter(int  a_lo, int  a_hi,
                                int  b_lo, int  b_hi,
                                GPointPair&   diam) const
    {
        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;
    }
    gdiam_real   brute_diameter(int  a_lo, int  a_hi,
                                int  b_lo, int  b_hi,
                                GPointPair&   diam,
                                const gdiam_point  dir) const
    {
        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 ], dir);
            }
 
        return  diam.distance;
    }
    gdiam_real   brute_diameter(const gdiam_point*   a_lo,
                                const gdiam_point*   a_hi,
                                const gdiam_point*   b_lo,
                                const gdiam_point*   b_hi,
                                GPointPair&   diam) const
    {
        for (const gdiam_point*   ind = a_lo; ind <= a_hi; ind++)
            for (const gdiam_point*    jnd = b_lo; jnd <= b_hi; jnd++)
            {
                diam.update_diam(*ind, *jnd);
            }
 
        return  diam.distance;
    }
    gdiam_real   brute_diameter(const gdiam_point*   a_lo,
                                const gdiam_point*   a_hi,
                                const gdiam_point*   b_lo,
                                const gdiam_point*   b_hi,
                                GPointPair&   diam,
                                const gdiam_point  dir) const
    {
        for (const gdiam_point*   ind = a_lo; ind <= a_hi; ind++)
            for (const gdiam_point*    jnd = b_lo; jnd <= b_hi; jnd++)
            {
                diam.update_diam(*ind, *jnd, dir);
            }
 
        return  diam.distance;
    }
 
    gdiam_point  getPoint(int  ind) const
    {
        return arr[ ind ];
    }
 
    GFSPTreeNode*   getRoot()
    {
        return  root;
    }
 
    GFSPTreeNode*   build_node(gdiam_point*   left,
                               gdiam_point*   right)
    {
        if (left > right)
        {
            printf("what!?\n");
            fflush(stdout);
            assert(left <= right);
        }
 
        while ((right > left)
               && (pnt_isEqual(*right, *left)))
        {
            right--;
        }
 
        GFSPTreeNode*   node = new  GFSPTreeNode(left, right);
 
        node->bbox.init();
 
        for (const gdiam_point*   ind = left; ind <= right; ind++)
        {
            node->bbox.bound(*ind);
        }
 
        node->bbox_diam = node->bbox.get_diam();
 
        node->bbox.center(node->center);
 
        return  node;
    }
 
    // return how many elements on left side.
    int    split_array(gdiam_point*    left,
                       gdiam_point*    right,
                       int  dim, const gdiam_real&   val)
    {
        const gdiam_point*    start_left = left;
 
        while (left < right)
        {
            if ((*left)[ dim ] < val)
            {
                left++;
            }
            else if ((*right)[ dim ] >= val)
            {
                right--;
            }
            else
            {
                gdiam_exchange(*right, *left);
            }
        }
 
        return  left - start_left;
    }
 
 
    void   split_node(GFSPTreeNode*   node)
    {
        int  dim, left_size;
        gdiam_real  split_val;
 
        if (node->left != NULL)
        {
            return;
        }
 
        GBBox&   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(GFSPPair&   pair,
                                  GPointPair&   max_pair_diam);
};
 
void  bbox_vertex(const GBBox&   bb, gdiam_point_t&   ver,
                  int  vert)
{
    ver[ 0 ] = ((vert & 0x1) != 0) ?
               bb.min_coord(0) : bb.max_coord(0);
    ver[ 1 ] = ((vert & 0x2) != 0) ?
               bb.min_coord(1) : bb.max_coord(1);
    ver[ 2 ] = ((vert & 0x4) != 0) ?
               bb.min_coord(2) : bb.max_coord(2);
}
 
 
gdiam_real  bbox_proj_dist(const GBBox&   bb1,
                           const GBBox&   bb2,
                           gdiam_point_cnt  proj)
{
    gdiam_point_t  a, b;
    gdiam_real  rl;
 
    rl = 0;
 
    for (int  ind = 0; ind < 8; ind++)
    {
        bbox_vertex(bb1, a, ind);
 
        //printf( "\n\n" );
        for (int  jnd = 0; jnd < 8; jnd++)
        {
            bbox_vertex(bb2, b, jnd);
            //pnt_dump( b );
            rl = max(rl, pnt_distance(a, b, proj));
        }
    }
 
    return  rl;
}
 
 
 
class  GFSPPair
{
private:
    GFSPTreeNode*   left, * right;
    double  max_diam;
 
public:
    void  init(GFSPTreeNode*   _left, GFSPTreeNode*   _right)
    {
        left = _left;
        right = _right;
 
        GBBox  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 );
    }
 
    void  init(GFSPTreeNode*   _left, GFSPTreeNode*   _right,
               gdiam_point  proj_dir, gdiam_real dist)
    {
        left = _left;
        right = _right;
 
        GBBox  bb;
 
        bb.init(left->getBBox(), right->getBBox());
        //printf( "\t\tbbox: \n" );
        //bb.dump();
 
        //max_diam = dist + _left->getBBox().get_diam()
        //    + _right->getBBox().get_diam();
        max_diam = max(max(bbox_proj_dist(_left->getBBox(),
                                          _right->getBBox(),
                                          proj_dir),
                           bbox_proj_dist(_left->getBBox(),
                                          _left->getBBox(),
                                          proj_dir)),
                       bbox_proj_dist(_right->getBBox(),
                                      _right->getBBox(),
                                      proj_dir));
 
        //printf( "diam: %g, smart diam: %g\n",
        //        max_diam,
        //        bbox_proj_dist( _left->getBBox(),
        //                        _right->getBBox(),
        //                        proj_dir ) );
 
        //printf( "dist: %g, %g, %g:\n",
        //        dist, _left->getBBox().get_diam(),
        //         _right->getBBox().get_diam() );
 
        //bb.get_diam_proj( proj_dir );
        //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;
    }
 
    GFSPTreeNode*   get_left()
    {
        return  left;
    }
    GFSPTreeNode*   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 GFSPTree&   tree,
                            GPointPair&   diam) const
    {
        return  tree.brute_diameter(left->ptr_pnt_left(),
                                    left->ptr_pnt_right(),
                                    right->ptr_pnt_left(),
                                    right->ptr_pnt_right(),
                                    diam);
    }
 
    double   directDiameter(const GFSPTree&   tree,
                            GPointPair&   diam,
                            const gdiam_point  dir) const
    {
        return  tree.brute_diameter(left->ptr_pnt_left(),
                                    left->ptr_pnt_right(),
                                    right->ptr_pnt_left(),
                                    right->ptr_pnt_right(),
                                    diam, dir);
    }
 
    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;
    }
};
 
 
class g_heap_pairs_p;
 
#define  UDM_SIMPLE      1
#define  UDM_BIG_SAMPLE  2
 
class  GTreeDiamAlg
{
private:
    GFSPTree  tree;
    //double  diam;
    GPointPair   pair_diam;
    int  points_num;
    int  threshold_brute;
    int  update_diam_mode;
public:
    GTreeDiamAlg(gdiam_point*   arr, int  size, int  mode)
    {
        tree.init(arr, size);
        //diam = 0;
        pair_diam.init(arr[ 0 ]);
        points_num = size;
        threshold_brute = 40;
        update_diam_mode = mode;
    }
 
    void  term()
    {
        tree.term();
    }
 
    int  size() const
    {
        return  points_num;
    }
    const GPointPair&   getDiameter() const
    {
        return  pair_diam;
    }
 
    int  nodes_number() const
    {
        return  tree.nodes_number();
    }
 
    void  addPairHeap(g_heap_pairs_p&   heap,
                      GFSPTreeNode*   left,
                      GFSPTreeNode*   right);
    void  addPairHeap(g_heap_pairs_p&   heap,
                      GFSPTreeNode*   left,
                      GFSPTreeNode*   right,
                      gdiam_point  proj,
                      GFSPPair&   father);
    void  split_pair(GFSPPair&   pair,
                     g_heap_pairs_p&   heap, double  eps);
    void  split_pair_proj(GFSPPair&   pair,
                          g_heap_pairs_p&   heap, double  eps,
                          gdiam_point  proj);
    void  compute_by_heap(double  eps);
    void  compute_by_heap_proj(double  eps, gdiam_point  proj);
 
    double  diameter()
    {
        return  pair_diam.distance;
    }
};
 
 
/*--- Heap implementation ---*/
 
typedef  int (*ptrCompareFunc)(void*   aPtr, void*   bPtr);
typedef  void*   voidPtr_t;
 
struct  heap_t
{
    voidPtr_t*   pArr;
    int  curr_size, max_size;
    ptrCompareFunc  pCompFunc;
};
 
void  heap_init(heap_t*   pHeap, ptrCompareFunc  _pCompFunc)
{
    assert(pHeap != NULL);
    assert(_pCompFunc != NULL);
 
    memset(pHeap, 0, sizeof(heap_t));
 
    pHeap->pCompFunc = _pCompFunc;
    pHeap->max_size = 100;
    pHeap->pArr = (voidPtr_t*)malloc(sizeof(void*)
                                     * pHeap->max_size);
    assert(pHeap->pArr != NULL);
    pHeap->curr_size = 0;
}
 
 
static void  resize(heap_t*   pHeap, int  size)
{
    int  max_sz;
    voidPtr_t*   pTmp;
 
    if (size <= pHeap->max_size)
    {
        return;
    }
 
    max_sz = size * 2;
    pTmp = (voidPtr_t*)malloc(max_sz * sizeof(void*));
    assert(pTmp != NULL);
    memset(pTmp, 0, max_sz * sizeof(void*));
    memcpy(pTmp, pHeap->pArr, pHeap->curr_size * sizeof(void*));
    free(pHeap->pArr);
    pHeap->pArr = pTmp;
}
 
 
void  heap_term(heap_t*   pHeap)
{
    assert(pHeap != NULL);
 
    if (pHeap->pArr != NULL)
    {
        free(pHeap->pArr);
    }
 
    memset(pHeap, 0, sizeof(heap_t));
}
 
 
void  heap_insert(heap_t*   pHeap, void*   pData)
{
    int  ind, father;
    voidPtr_t  pTmp;
 
    resize(pHeap, pHeap->curr_size + 1);
 
    ind = pHeap->curr_size;
    pHeap->curr_size++;
 
    pHeap->pArr[ ind ] = pData;
 
    while (ind > 0)
    {
        father = (ind - 1) / 2;
 
        if ((*pHeap->pCompFunc)(pHeap->pArr[ ind ],
                                pHeap->pArr[ father ]) <= 0)
        {
            break;
        }
 
        pTmp = pHeap->pArr[ ind ];
        pHeap->pArr[ ind ] = pHeap->pArr[ father ];
        pHeap->pArr[ father ] = pTmp;
 
        ind = father;
    }
}
 
 
void*   heap_top(heap_t*   pHeap)
{
    return  pHeap->pArr[ 0 ];
}
 
 
void*   heap_delete_max(heap_t*   pHeap)
{
    void*   res;
    void*   pTmp;
    int  ind, left, right, max_ind;
 
    if (pHeap->curr_size <= 0)
    {
        return  NULL;
    }
 
    res = pHeap->pArr[ 0 ];
 
    //printf( "pHeap->curr_size: %d\n", pHeap->curr_size );
    //printf( "res= %p\n", res );
 
    pHeap->curr_size--;
    pHeap->pArr[ 0 ] = pHeap->pArr[ pHeap->curr_size ];
    pHeap->pArr[ pHeap->curr_size ] = NULL;
    ind = 0;
 
    while (ind < pHeap->curr_size)
    {
        left = 2 * ind + 1;
        right = 2 * ind + 2;
 
        if (left >= pHeap->curr_size)
        {
            break;
        }
 
        if (right >= pHeap->curr_size)
        {
            right = left;
        }
 
        if ((*pHeap->pCompFunc)(pHeap->pArr[ left ],
                                pHeap->pArr[ right ]) <= 0)
        {
            max_ind = right;
        }
        else
        {
            max_ind = left;
        }
 
        if ((*pHeap->pCompFunc)(pHeap->pArr[ ind  ],
                                pHeap->pArr[ max_ind ]) >= 0)
        {
            break;
        }
 
        pTmp = pHeap->pArr[ ind ];
        pHeap->pArr[ ind ] = pHeap->pArr[ max_ind ];
        pHeap->pArr[ max_ind ] = pTmp;
 
        ind = max_ind;
    }
 
    return  res;
}
 
 
bool  heap_is_empty(heap_t*   pHeap)
{
    assert(pHeap != NULL);
 
    return (pHeap->curr_size <= 0);
}
 
 
int   compare_pairs(void*   _a, void*   _b)
{
    const GFSPPair&   a(*(GFSPPair*)_a);
    const GFSPPair&   b(*(GFSPPair*)_b);
 
    if (a.maxDiam() < b.maxDiam())
    {
        return  -1;
    }
 
    if (a.maxDiam() > b.maxDiam())
    {
        return  1;
    }
 
    return  0;
}
 
 
 
class  g_heap_pairs_p
{
private:
    heap_t  heap;
 
public:
    g_heap_pairs_p()
    {
        heap_init(&heap, compare_pairs);
    }
 
    ~g_heap_pairs_p()
    {
        while (size() > 0)
        {
            pop();
        }
 
        heap_term(&heap);
    }
 
    void  push(GFSPPair&   pair)
    {
        //printf( "pushing: (%p, %p)\n", pair.get_left(),
        //        pair.get_right() );
        GFSPPair*   p_pair = new GFSPPair(pair);
        heap_insert(&heap, (void*)p_pair);
    }
 
    int  size() const
    {
        return  heap.curr_size;
    }
    GFSPPair  top()
    {
        return  *(GFSPPair*)heap_top(&heap);
    }
    void  pop()
    {
        GFSPPair*   ptr = (GFSPPair*)heap_delete_max(&heap);
 
        delete  ptr;
    }
};
 
 
/* heap.implementation end */
 
 
void  computeExtremePair(const gdiam_point*   arr, int  size,
                         int  dim, GPointPair&   pp)
{
    pp.init(arr[ 0 ]);
 
    for (int  ind = 1; ind < size; ind++)
    {
        const gdiam_point  pnt(arr[ ind ]);
 
        if (pnt[ dim ] < pp.p[ dim ])
        {
            pp.p = pnt;
        }
 
        if (pnt[ dim ] > pp.q[ dim ])
        {
            pp.q = pnt;
        }
    }
 
    pp.distance = pnt_distance(pp.p, pp.q);
}
 
 
 
void  GTreeDiamAlg::compute_by_heap(double  eps)
{
    g_heap_pairs_p  heap;
    int  heap_limit, heap_delta;
 
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
 
    GFSPTreeNode*   root = tree.getRoot();
 
    computeExtremePair(tree.getPoints(), points_num,
                       root->getBBox().getLongestDim(),
                       pair_diam);
 
    GFSPPair  root_pair;
 
    root_pair.init(root, root);
 
    heap.push(root_pair);
    //queue_a.pushx( root_pair );
 
    int  count = 0;
 
    while (heap.size() > 0)
    {
        //printf( "heap.size: %d\n", heap.size() );
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        GFSPPair  pair = heap.top();
        heap.pop();
        split_pair(pair, heap, eps);
 
        //if  ( ( count & 0x3ff ) == 0 ) {
        //    printf( "heap.size: %d\n", heap.size() );
        //}
        if ((count & 0x3ff) == 0)
        {
            if (((int)heap.size()) > heap_limit)
            {
                threshold_brute *= 2;
                printf("threshold_brute: %d\n", threshold_brute);
                heap_limit += heap_delta;
            }
        }
 
        count++;
    }
}
 
 
void  GTreeDiamAlg::compute_by_heap_proj(double  eps,
        gdiam_point  proj)
{
    g_heap_pairs_p  heap;
    int  heap_limit, heap_delta;
    GPointPair   pair_diam_x;
 
 
    //pnt_dump( proj );
    //printf( "length = %g\n", pnt_length( proj ) );
 
    heap_limit = HEAP_LIMIT;
    heap_delta = HEAP_DELTA;
 
    GFSPTreeNode*   root = tree.getRoot();
 
    computeExtremePair(tree.getPoints(), points_num,
                       root->getBBox().getLongestDim(),
                       pair_diam_x);
    pair_diam.init(pair_diam_x.p, pair_diam_x.q, proj);
 
    GFSPPair  root_pair;
 
    root_pair.init(root, root, proj, 0);
 
    //printf( "root pair distance: %g\n", root_pair.maxDiam() );
    heap.push(root_pair);
    //queue_a.pushx( root_pair );
 
    int  count = 0;
 
    while (heap.size() > 0)
    {
        //printf( "heap.size: %d\n", heap.size() );
        //printf( "%d: curr_diam: %g\n", count++, pair_diam.distance );
        GFSPPair  pair = heap.top();
        heap.pop();
        //printf( "pair distance: %g\n", pair.maxDiam() );
        split_pair_proj(pair, heap, eps, proj);
 
        //if  ( ( count & 0x3ff ) == 0 ) {
        //    printf( "heap.size: %d\n", heap.size() );
        //}
        if ((count & 0x3ff) == 0)
        {
            if (((int)heap.size()) > heap_limit)
            {
                threshold_brute *= 2;
                printf("threshold_brute: %d\n", threshold_brute);
                heap_limit += heap_delta;
            }
        }
 
        count++;
        //printf( "Poping!\n" );
    }
}
 
 
void  GTreeDiamAlg::addPairHeap(g_heap_pairs_p&   heap,
                                GFSPTreeNode*   left,
                                GFSPTreeNode*   right)
{
    if (update_diam_mode == UDM_SIMPLE)
    {
        const gdiam_point  p(*(left->ptr_pnt_left()));
        const gdiam_point  q(*(right->ptr_pnt_left()));
 
        pair_diam.update_diam(p, q);
    }
    else if (update_diam_mode == UDM_BIG_SAMPLE)
    {
        if ((left->size() < 100)  || (right->size() < 100))
        {
            const gdiam_point  p(*(left->ptr_pnt_left()));
            const gdiam_point  q(*(right->ptr_pnt_left()));
 
            pair_diam.update_diam(p, q);
        }
        else
        {
#define  UMD_SIZE        5
            gdiam_point  arr_left[ UMD_SIZE ],
                         arr_right[ UMD_SIZE ];
 
            for (int  ind = 0; ind < UMD_SIZE; ind++)
            {
                const gdiam_point  p(*(left->ptr_pnt_rand()));
                const gdiam_point  q(*(right->ptr_pnt_rand()));
                arr_left[ ind ] = p;
                arr_right[ ind ] = q;
 
                //pair_diam.update_diam( p, q );
            }
 
            for (int  ind = 0; ind < UMD_SIZE; ind++)
                for (int  jnd = 0; jnd < UMD_SIZE; jnd++)
                    pair_diam.update_diam(arr_left[ ind ],
                                          arr_right[ jnd ]);
        }
    }
    else
    {
        assert(false);
    }
 
    //printf( "old_diam; %g\n", diam );
    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
 
    GFSPPair  pair;
 
    pair.init(left, right);
 
    if (pair.maxDiam() <= pair_diam.distance)
    {
        return;
    }
 
    heap.push(pair);
}
 
 
void  GTreeDiamAlg::addPairHeap(g_heap_pairs_p&   heap,
                                GFSPTreeNode*   left,
                                GFSPTreeNode*   right,
                                gdiam_point  proj,
                                GFSPPair&   father)
{
    const gdiam_point  p(*(left->ptr_pnt_left()));
    const gdiam_point  q(*(right->ptr_pnt_left()));
 
    //printf( "addPairHeap( ?, %p, %p, ? )\n", left, right );
    pair_diam.update_diam(p, q, proj);
    //printf( "old_diam; %g\n", diam );
    //diam = max( diam, pnt_distance( p, q ) );
    //printf( "new_diam; %g\n", diam );
 
    GFSPPair  pair;
 
    pair.init(left, right, proj,
              pnt_distance(p, q, proj));
 
    if (pair.maxDiam() <= pair_diam.distance)
    {
        return;
    }
 
    //printf( "pair.maxDiam: %g, father.maxDiam: %g\n",
    //        pair.maxDiam(),  father.maxDiam() );
    //printf( "pair_diam.distance: %g\n", pair_diam.distance );
    //pnt_dump( proj );
    heap.push(pair);
}
 
 
void  GTreeDiamAlg::split_pair_proj(GFSPPair&   pair,
                                    g_heap_pairs_p&   heap, double  eps,
                                    gdiam_point  proj)
{
    bool  f_is_left_splitable, f_is_right_splitable;
 
    //printf( "pair.maxDiam: %g, limit: %g\n",
    //        pair.maxDiam(), ((1.0 + eps) * pair_diam.distance ));
 
    if (pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance))
    {
        return;
    }
 
    if (pair.isBelowThreshold(threshold_brute))
    {
        pair.directDiameter(tree, pair_diam, proj);
        return;
    }
 
    //printf( "Curr pair: (%p, %p)\n", pair.get_left(), pair.get_right() );
    f_is_left_splitable = pair.get_left()->isSplitable();
    f_is_right_splitable = pair.get_right()->isSplitable();
    assert(f_is_left_splitable  ||  f_is_right_splitable);
 
    if (f_is_left_splitable)
    {
        tree.split_node(pair.get_left());
    }
 
    if (f_is_right_splitable)
    {
        tree.split_node(pair.get_right());
    }
 
    if (f_is_left_splitable
        &&  f_is_right_splitable)
    {
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right()->get_left(),
                    proj, pair);
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right()->get_right(),
                    proj, pair);
 
        // to avoid exponential blowup
        if (pair.get_left() != pair.get_right())
            addPairHeap(heap,
                        pair.get_left()->get_right(),
                        pair.get_right()->get_left(),
                        proj, pair);
 
        addPairHeap(heap,
                    pair.get_left()->get_right(),
                    pair.get_right()->get_right(),
                    proj, pair);
        return;
    }
 
    if (f_is_left_splitable)
    {
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right(), proj, pair);
        addPairHeap(heap,
                    pair.get_left()->get_right(),
                    pair.get_right(), proj, pair);
        return;
    }
 
    if (f_is_right_splitable)
    {
        addPairHeap(heap,
                    pair.get_left(),
                    pair.get_right()->get_left(), proj, pair);
        addPairHeap(heap,
                    pair.get_left(),
                    pair.get_right()->get_right(), proj, pair);
        return;
    }
}
 
 
void  GTreeDiamAlg::split_pair(GFSPPair&   pair,
                               g_heap_pairs_p&   heap,
                               double  eps)
{
    bool  f_is_left_splitable, f_is_right_splitable;
 
    if (pair.maxDiam() <= ((1.0 + eps) * pair_diam.distance))
    {
        return;
    }
 
    if (pair.isBelowThreshold(threshold_brute))
    {
        pair.directDiameter(tree, pair_diam);
        return;
    }
 
    f_is_left_splitable = pair.get_left()->isSplitable();
    f_is_right_splitable = pair.get_right()->isSplitable();
    assert(f_is_left_splitable  ||  f_is_right_splitable);
 
    if (f_is_left_splitable)
    {
        tree.split_node(pair.get_left());
    }
 
    if (f_is_right_splitable)
    {
        tree.split_node(pair.get_right());
    }
 
    if (f_is_left_splitable
        &&  f_is_right_splitable)
    {
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right()->get_left());
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right()->get_right());
 
        // to avoid exponential blowup
        if (pair.get_left() != pair.get_right())
            addPairHeap(heap,
                        pair.get_left()->get_right(),
                        pair.get_right()->get_left());
 
        addPairHeap(heap,
                    pair.get_left()->get_right(),
                    pair.get_right()->get_right());
        return;
    }
 
    if (f_is_left_splitable)
    {
        addPairHeap(heap,
                    pair.get_left()->get_left(),
                    pair.get_right());
        addPairHeap(heap,
                    pair.get_left()->get_right(),
                    pair.get_right());
        return;
    }
 
    if (f_is_right_splitable)
    {
        addPairHeap(heap,
                    pair.get_left(),
                    pair.get_right()->get_left());
        addPairHeap(heap,
                    pair.get_left(),
                    pair.get_right()->get_right());
        return;
    }
}
 
 
GPointPair   gdiam_approx_diam(gdiam_point*   start, int  size,
                               gdiam_real  eps)
{
    //gdiam_real  diam;
    GPointPair  pair;
 
    GTreeDiamAlg*   pAlg = new  GTreeDiamAlg(start, size, UDM_SIMPLE);
    pAlg->compute_by_heap(eps);
 
    pair = pAlg->getDiameter();
 
    delete  pAlg;
 
    return  pair;
}
 
 
GPointPair   gdiam_approx_diam_UDM(gdiam_point*   start, int  size,
                                   gdiam_real  eps)
{
    //gdiam_real  diam;
    GPointPair  pair;
 
    GTreeDiamAlg*   pAlg = new  GTreeDiamAlg(start, size, UDM_BIG_SAMPLE);
    pAlg->compute_by_heap(eps);
 
    pair = pAlg->getDiameter();
 
    delete  pAlg;
 
    return  pair;
}
 
 
gdiam_point*   gdiam_convert(gdiam_real*   start, int  size)
{
    gdiam_point*   p_arr;
 
    assert(start != NULL);
    assert(size > 0);
 
    p_arr = (gdiam_point*)malloc(sizeof(gdiam_point) * size);
    assert(p_arr != NULL);
 
    for (int  ind = 0; ind < size; ind++)
    {
        p_arr[ ind ] = (gdiam_point) & (start[ ind * 3 ]);
    }
 
    return  p_arr;
}
 
 
GPointPair   gdiam_approx_diam_pair(gdiam_real*   start, int  size,
                                    gdiam_real  eps)
{
    gdiam_point*   p_arr;
    GPointPair  pair;
 
    p_arr = gdiam_convert(start, size);
    pair = gdiam_approx_diam(p_arr, size, eps);
    free(p_arr);
 
    return  pair;
}
 
 
GPointPair   gdiam_approx_diam_pair_UDM(gdiam_real*   start, int  size,
                                        gdiam_real  eps)
{
    gdiam_point*   p_arr;
    GPointPair  pair;
 
    p_arr = gdiam_convert(start, size);
    pair = gdiam_approx_diam_UDM(p_arr, size, eps);
    free(p_arr);
 
    return  pair;
}
 
 
gdiam_bbox   gdiam_approx_const_mvbb(gdiam_point*   start, int  size,
                                     gdiam_real  eps,
                                     GBBox*   p_ap_bbox)
{
    //gdiam_real  diam;
    GPointPair  pair, pair_2;
 
    GTreeDiamAlg*   pAlg = new  GTreeDiamAlg(start, size, UDM_SIMPLE);
    pAlg->compute_by_heap(eps);
 
    pair = pAlg->getDiameter();
 
    /* Comput ethe resulting diameter */
    gdiam_point_t  dir, dir_2, dir_3;
 
    dir[ 0 ] = pair.p[ 0 ] - pair.q[ 0 ];
    dir[ 1 ] = pair.p[ 1 ] - pair.q[ 1 ];
    dir[ 2 ] = pair.p[ 2 ] - pair.q[ 2 ];
    pnt_normalize(dir);
 
    //    printf( "Before computing second direction!\n" );
    //fflush( stdout );
    pAlg->compute_by_heap_proj(eps, dir);
    //printf( "second direction computed!\n" );
    //fflush( stdout );
 
    pair_2 = pAlg->getDiameter();
    dir_2[ 0 ] = pair_2.p[ 0 ] - pair_2.q[ 0 ];
    dir_2[ 1 ] = pair_2.p[ 1 ] - pair_2.q[ 1 ];
    dir_2[ 2 ] = pair_2.p[ 2 ] - pair_2.q[ 2 ];
 
    gdiam_real prd;
 
    prd = pnt_dot_prod(dir, dir_2);
    dir_2[ 0 ] = dir_2[ 0 ] - prd * dir[ 0 ];
    dir_2[ 1 ] = dir_2[ 1 ] - prd * dir[ 1 ];
    dir_2[ 2 ] = dir_2[ 2 ] - prd * dir[ 2 ];
 
    pnt_normalize(dir);
    pnt_normalize(dir_2);
 
    pnt_cross_prod(dir, dir_2, dir_3);
    pnt_normalize(dir_3);
 
    gdiam_bbox  bbox;
    GBBox  ap_bbox;
 
    if ((pnt_length(dir_2) < 1e-6)
        && (pnt_length(dir_3) < 1e-6))
    {
        gdiam_generate_orthonormal_base(dir, dir_2, dir_3);
    }
 
    if ((pnt_length(dir) == 0.0)
        || (pnt_length(dir_2) < 1e-6)
        || (pnt_length(dir_3) < 1e-6))
    {
        gdiam_generate_orthonormal_base(dir, dir_2, dir_3);
        pnt_dump(dir);
        pnt_dump(dir_2);
        pnt_dump(dir_3);
 
        fflush(stdout);
        fflush(stderr);
        assert(false);
    }
 
    bbox.init(dir, dir_2, dir_3);
    ap_bbox.init();
 
    for (int  ind = 0; ind < size; ind++)
    {
        bbox.bound(start[ ind ]);
        ap_bbox.bound(start[ ind ]);
    }
 
    //printf( "Axis parallel bounding box:\n" );
    //ap_bbox.dump();
 
    //printf( "Arbitrarily oriented bounding box:\n" );
    //bbox.dump();
 
    delete  pAlg;
 
    if (ap_bbox.volume() < bbox.volume())
    {
        bbox.init(ap_bbox);
    }
 
    if (p_ap_bbox != NULL)
    {
        *p_ap_bbox = ap_bbox;
    }
 
    return  bbox;
}
 
 
gdiam_bbox   gdiam_approx_const_mvbb(gdiam_real*   start, int  size,
                                     gdiam_real  eps)
{
    gdiam_point*   p_arr;
    gdiam_bbox  gbbox;
 
    p_arr = gdiam_convert(start, size);
    gbbox = gdiam_approx_const_mvbb(p_arr, size, eps, NULL);
    free(p_arr);
 
    return   gbbox;
}
 
/*-----------------------------------------------------------------------
 * Code for computing best bounding box in a given drection
\*-----------------------------------------------------------------------*/
 
void  gdiam_generate_orthonormal_base(gdiam_point  in,
                                      gdiam_point  out1,
                                      gdiam_point  out2)
{
    assert(pnt_length(in) > 0.0);
 
    pnt_normalize(in);
 
    // stupid cases..
    if (in[ 0 ] == 0.0)
    {
        if (in[ 1 ] == 0.0)
        {
            pnt_init_normalize(out1, 1, 0, 0);
            pnt_init_normalize(out2, 0, 1, 0);
            return;
        }
 
        if (in[ 2 ] == 0.0)
        {
            pnt_init_normalize(out1, 1, 0, 0);
            pnt_init_normalize(out2, 0, 0, 1);
            return;
        }
 
        pnt_init_normalize(out1, 0, -in[ 2 ], in[ 1 ]);
        pnt_init_normalize(out2, 1, 0, 0);
        return;
    }
 
    if ((in[ 1 ] == 0.0)  && (in[ 2 ] == 0.0))
    {
        pnt_init_normalize(out1, 0, 1, 0);
        pnt_init_normalize(out2, 0, 0, 1);
        return;
    }
 
    if (in[ 1 ] == 0.0)
    {
        pnt_init_normalize(out1, -in[ 2 ], 0, in[ 0 ]);
        pnt_init_normalize(out2, 0, 1, 0);
        return;
    }
 
    if (in[ 2 ] == 0.0)
    {
        pnt_init_normalize(out1, -in[ 1 ], in[ 0 ], 0);
        pnt_init_normalize(out2, 0, 0, 1);
        return;
    }
 
    // all entries in the vector are not zero.
    pnt_init_normalize(out1, -in[ 1 ], in[ 0 ], 0);
    pnt_cross_prod(in, out1, out2);
 
    pnt_normalize(out2);
}
 
 
class  point2d
{
public:
    gdiam_real  x, y;
    gdiam_point  src;
 
    void  init(gdiam_point  pnt,
               gdiam_point  base_x,
               gdiam_point  base_y)
    {
        src = pnt;
        x = pnt_dot_prod(pnt, base_x);
        y = pnt_dot_prod(pnt, base_y);
    }
 
    gdiam_real  dist(const point2d&   pnt)
    {
        return  sqrt((x - pnt.x) * (x - pnt.x)
                     + (y - pnt.y) * (y - pnt.y));
    }
 
    void  dump() const
    {
        printf("(%g, %g)\n", x, y);
    }
    bool  equal(const point2d&   pnt) const
    {
        return ((x == pnt.x)  && (y == pnt.y));
    }
    bool  equal_real(const point2d&   pnt) const
    {
        return ((fabs(x - pnt.x) < 1e-8)
                && (fabs(y - pnt.y) < 1e-8));
    }
};
 
typedef  point2d*   point2d_ptr;
class  vec_point_2d : public std::vector<point2d_ptr> {};
 
 
 
inline ldouble     Area(const point2d&   a,
                        const point2d&   b,
                        const point2d&   c)
{
    double  x1, y1, x2, y2, len;
 
    x1 = b.x - a.x;
    y1 = b.y - a.y;
 
    x2 = c.x - a.x;
    y2 = c.y - a.y;
 
    if ((x1 != 0.0)  || (y1 != 0.0))
    {
        len = sqrt(x1 * x1 + y1 * y1);
        x1 /= len;
        y1 /= len;
    }
 
    if ((x2 != 0.0)  || (y2 != 0.0))
    {
        len = sqrt(x2 * x2 + y2 * y2);
        x2 /= len;
        y2 /= len;
    }
 
    ldouble  area;
 
    area = x1 * y2 - x2 * y1;
    //area = (ldouble)a.x * (ldouble)b.y - (ldouble)a.y * (ldouble)b.x +
    //    (ldouble)a.y * (ldouble)c.x - (ldouble)a.x * (ldouble)c.y +
    //    (ldouble)b.x * (ldouble)c.y - (ldouble)c.x * (ldouble)b.y;
 
    //printf( "area: %g\n", area );
    return  area;
}
 
 
inline int     AreaSign(const point2d&   a,
                        const point2d&   b,
                        const point2d&   c)
{
    ldouble  area;
 
    area = a.x * b.y - a.y * b.x +
           a.y * c.x - a.x * c.y +
           b.x * c.y - c.x * b.y;
 
    //printf( "area: %g\n", area );
    return ((area < (ldouble)0.0) ? -1 :
            ((area > (ldouble)0.0) ? 1 : 0));
}
 
inline  bool    Left(const point2d&   a,
                     const point2d&   b,
                     const point2d&   c)
{
    return  AreaSign(a, b, c) > 0;
}
 
inline  bool    Left_colinear(const point2d&   a,
                              const point2d&   b,
                              const point2d&   c)
{
    return  AreaSign(a, b, c) >= 0;
}
 
 
struct CompareByAngle
{
public:
    point2d  base;
 
    bool operator()(const point2d_ptr&   a, const point2d_ptr&    b)
    {
        int  sgn;
        gdiam_real  len1, len2;
 
        if (a->equal(*b))
        {
            return  false;
        }
 
        assert(a != NULL);
        assert(b != NULL);
 
        if (a->equal(base))
        {
            //assert( false );
            return  true;
        }
 
        if (b->equal(base))
        {
            //   assert( false );
            return  false;
        }
 
        //printf( "before: %p %p...\n", a, b );
        //fflush( stdout );
        /*
        if  ( a->equal( base ) )
            return  true;
        if  ( b->equal( base ) )
            return  false;
        */
        //printf( "a: (%g, %g)\n", a->x, a->y );
        //printf( "b: (%g, %g)\n", b->x, b->y );
        //fflush( stdout );
 
        sgn = AreaSign(base, *a, *b);
 
        if (sgn != 0)
        {
            return (sgn > 0);
        }
 
        len1 = base.dist(*a);
        len2 = base.dist(*b);
 
        //printf( "bogi!\n" );
        //fflush( stdout );
 
        return (len1 > len2);
    }
};
 
 
point2d_ptr  get_min_point(vec_point_2d&   in,
                           int&   index)
{
    point2d_ptr  min_pnt = in[ 0 ];
 
    index = 0;
 
    for (int  ind = 0; ind < (int)in.size(); ind++)
    {
        if (in[ ind ]->y < min_pnt->y)
        {
            min_pnt = in[ ind ];
            index = ind;
        }
        else if ((in[ ind ]->y == min_pnt->y)
                 && (in[ ind ]->x < min_pnt->x))
        {
            min_pnt = in[ ind ];
            index = ind;
        }
    }
 
    return  min_pnt;
}
 
void  dump(vec_point_2d&    vec)
{
    for (int  ind = 0; ind < (int)vec.size(); ind++)
    {
        printf("-- %11d (%-11g, %-11g)\n",
               ind, vec[ ind ]->x,
               vec[ ind ]->y);
    }
 
    printf("\n\n");
}
 
 
static void  print_pnt(point2d_ptr  pnt)
{
    printf("(%g, %g)\n", pnt->x, pnt->y);
}
 
 
void   verify_convex_hull(vec_point_2d&   in, vec_point_2d&   ch)
{
    ldouble  area;
 
    //dump( ch );
    //fflush( stdout );
    //fflush( stderr );
    for (int  ind = 0; ind < (int)(ch.size() - 2); ind++)
        assert(Left(*(ch[ ind ]), *(ch[ ind + 1 ]),
                    *(ch[ ind + 2 ])));
 
    assert(Left(*(ch[ ch.size() - 2 ]), *(ch[ ch.size() - 1 ]),
                *(ch[ 0 ])));
    assert(Left(*(ch[ ch.size() - 1 ]), *(ch[ 0 ]),
                *(ch[ 1 ])));
 
    for (int  ind = 0; ind < (int)(in.size() - 2); ind++)
    {
        for (int  jnd = 0; jnd < (int)(ch.size() - 1); jnd++)
        {
            if (ch[ jnd ]->equal_real(*(in[ ind ])))
            {
                continue;
            }
 
            if (ch[ jnd + 1 ]->equal(*(in[ ind ])))
            {
                continue;
            }
 
            if (! Left_colinear(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
                                *(in[ ind ])))
            {
                area = Area(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
                            *(in[ ind ]));
 
                if (fabs(area) < 1e-12)
                {
                    continue;
                }
 
                printf("Failure in progress!\n\n");
                print_pnt(ch[ jnd ]);
                print_pnt(ch[ jnd + 1 ]);
                print_pnt(in[ ind ]);
 
                //dump( ch );
                printf("ch[ jnd ]: (%g, %g)\n",
                       ch[ jnd ]->x, ch[ jnd ]->y);
                printf("ch[ jnd + 1 ]: (%g, %g)\n",
                       ch[ jnd + 1 ]->x, ch[ jnd + 1 ]->y);
                printf("ch[ ind ]: (%g, %g)\n",
                       in[ ind ]->x, in[ ind ]->y);
                printf("Area: %g\n", (double)Area(*(ch[ jnd ]),
                                                  *(ch[ jnd + 1 ]),
                                                  *(in[ ind ])));
                printf("jnd: %d, jnd + 1: %d, ind: %d\n",
                       jnd, jnd + 1, ind);
                fflush(stdout);
                fflush(stderr);
 
                assert(Left(*(ch[ jnd ]), *(ch[ jnd + 1 ]),
                            *(in[ ind ])));
            }
        }
 
        if (ch[ 0 ]->equal(*(in[ ind ])))
        {
            continue;
        }
 
        if (ch[ ch.size() - 1 ]->equal(*(in[ ind ])))
        {
            continue;
        }
 
        assert(Left(*(ch[ ch.size() - 1 ]),
                    *(ch[ 0 ]),
                    *(in[ ind ])));
    }
 
    printf("Convex hull seems to be ok!\n");
}
 
 
 
static  void   remove_duplicate(vec_point_2d&   in,
                                point2d_ptr  ptr,
                                int  start)
{
    int  dest;
 
    dest = start;
 
    while (start < (int)in.size())
    {
        if (in[ start ]->equal_real(*ptr))
        {
            start++;
            continue;
        }
 
        in[ dest++ ] = in[ start++ ];
    }
 
    while ((int)in.size() > dest)
    {
        in.pop_back();
    }
}
 
 
static  void   remove_consecutive_dup(vec_point_2d&   in)
{
    int  dest, pos;
 
    if (in.size() < 1)
    {
        return;
    }
 
    dest = pos = 0;
 
    // copy first entry
    in[ dest++ ] = in[ pos++ ];
 
    while (pos < ((int)in.size() - 1))
    {
        if (in[ pos - 1 ]->equal_real(*(in[ pos ])))
        {
            pos++;
            continue;
        }
 
        in[ dest++ ] = in[ pos++ ];
    }
 
    in[ dest++ ] = in[ pos++ ];
 
    while ((int)in.size() > dest)
    {
        in.pop_back();
    }
}
 
 
void  convex_hull(vec_point_2d&   in, vec_point_2d&   out)
{
 
    CompareByAngle  comp;
    int  ind, position;
 
    assert(in.size() > 1);
 
    comp.base = *(get_min_point(in, position));
    std::swap(in[ 0 ], in[ position ]);
 
    remove_duplicate(in, in[ 0 ], 1);
 
    /*
    printf( "comp.base: (%g, %g)\n",
            comp.base.x,
            comp.base.y );
    */
    for (int  ind = 0; ind < (int)in.size(); ind++)
    {
        assert(in[ ind ] != NULL);
    }
 
    /*
    if  ( in.size() == 24 ) {
        dump( in );
        fflush( stdout );
        fflush( stderr );
        }*/
 
    //printf( "sort( %d, %d, comp )\n", 1, in.size() );
    sort(in.begin() + 1, in.end(), comp);
    remove_consecutive_dup(in);
 
    //dump( in );
    /*
    for  ( int  ind = 0; ind < (int)in.size(); ind++ ) {
        double  x_delta, y_delta;
 
        x_delta = in[ ind ]->x - comp.base.x;
        y_delta = in[ ind ]->y - comp.base.y;
 
        printf( "-- %11d (%-11g, %-11g)  atan2: %-11g\n",
                ind, in[ ind ]->x,
                in[ ind ]->y,
                atan2( y_delta, x_delta ) );
                }*/
 
    //fflush( stdout );
    //printf( "pushing...\n" );
    //fflush( stdout );
    out.push_back(in[ in.size() - 1 ]);
    out.push_back(in[ 0 ]);
 
    // perform the graham scan
    ind = 1;
    int  last;
 
    while (ind < (int)in.size())
    {
        if (out[ out.size() - 1 ]->equal_real(*(in[ ind ])))
        {
            ind++;
            continue;
        }
 
        last = out.size();
        assert(last > 1);
 
        if (Left(*(out[ last - 2 ]),
                 *(out[ last - 1 ]),
                 *(in[ ind ])))
        {
            if (! out[ last - 1 ]->equal(*(in[ ind ])))
            {
                out.push_back(in[ ind ]);
            }
 
            ind++;
        }
        else
        {
            if (out.size() < 3)
            {
                dump(out);
                printf("Input:\n");
                dump(in);
                printf("ind: %d\n", ind);
                fflush(stdout);
                assert(out.size() > 2);
            }
 
            out.pop_back();
        }
    }
 
    // we pushed in[ in.size() -1 ] twice in the output
    out.pop_back();
 
    //verify_convex_hull( in, out );
}
 
 
struct bbox_2d_info
{
    gdiam_point_2d_t  vertices[ 4 ];
    gdiam_real  area;
 
    void  get_dir(int  ind, gdiam_point_2d  out)
    {
        out[ 0 ] = vertices[ ind ][ 0 ]
                   - vertices[(ind + 1) % 4 ][ 0 ];
        out[ 1 ] = vertices[ ind ][ 1 ]
                   - vertices[(ind + 1) % 4 ][ 1 ];
    }
 
    void  dump()
    {
        printf("--- bbox 2d ----------------------------\n");
 
        for (int  ind = 0; ind < 4; ind++)
            printf("%d: (%g, %g)\n", ind, vertices[ ind ][0],
                   vertices[ ind ][1]);
 
        for (int  ind = 0; ind < 4; ind++)
        {
            gdiam_point_2d_t  dir;
 
            get_dir(ind, dir);
 
            printf("dir %d: (%g, %g)\n", ind, dir[0],
                   dir[1]);
        }
 
        printf("\\----------------------------------\n");
    }
 
};
 
#define  EPS  1e-6
inline  bool   eq_real(gdiam_real  a, gdiam_real  b)
{
    return (((b - EPS) < a)
            && (a < (b + EPS)));
}
 
class  MinAreaRectangle
{
private:
    vec_point_2d   ch;
    gdiam_real*   angles;
 
public:
    void  dump_ch()
    {
        for (int  ind = 0; ind < (int)ch.size(); ind++)
        {
            printf("%d: (%g, %g)\n", ind, ch[ ind ]->x,
                   ch[ ind ]->y);
        }
    }
 
    void  compute_edge_dir(int  u)
    {
        int  u2;
        double  ang;
        double  x1, y1, x2, y2;
 
        u2 = (u + 1) % ch.size();
 
        x1 = ch[ u ]->x;
        y1 = ch[ u ]->y;
        x2 = ch[ u2 ]->x;
        y2 = ch[ u2 ]->y;
 
        ang = atan2(y2 - y1, x2 - x1);
 
        if (ang < 0)
        {
            ang += 2.0 * PI;
        }
 
        angles[ u ] = ang;
 
        // floating point nonsence. A huck to solve this...
        if ((u > 0)
            && (angles[ u ] < angles[ u - 1 ])
            &&   eq_real(angles[ u ], angles[ u - 1 ]))
        {
            angles[ u ] = angles[ u - 1 ];
        }
 
        /*
        printf( "angles[ %4d ]: %.20f ", u, ang );
        if  ( u > 0 )
            printf( "%2d", (angles[ u ] >= angles[ u - 1 ] ) );
        if  ( u > 1 ) {
            printf( " area: %.20f ", (double)Area( *( ch[ u - 2 ] ),
                                                  *( ch[ u - 1 ] ),
                                                  *( ch[ u ] ) ) );
        }
        printf( "\n" );
        */
    }
 
    /* Finding the first vertex whose edge is over a given angle */
    int  find_vertex(int  start_vertex,
                     double  trg_angle)
    {
        double  prev_angle, curr_angle;
        bool  f_found;
        int  ver, prev_vertex;
 
        f_found = false;
 
        prev_vertex = start_vertex - 1;
 
        if (prev_vertex < 0)
        {
            prev_vertex = ch.size() - 1;
        }
 
        prev_angle = angles[ prev_vertex ];
        ver = start_vertex;
 
        while (! f_found)
        {
            curr_angle = angles[ ver ];
 
            if (prev_angle <= curr_angle)
                f_found = ((prev_angle < trg_angle)
                           && (trg_angle <= curr_angle));
            else
                f_found = ((prev_angle < trg_angle)
                           || (trg_angle <= curr_angle));
 
            if (f_found)
            {
                break;
            }
            else
            {
                prev_angle = curr_angle;
            }
 
            ver = (ver + 1) % ch.size();
        }
 
        return  ver;
    }
 
 
    /* Computing the intersection point of two lines, each given by a
    point and slope */
    void  compute_crossing(int  a_ind, double  a_angle,
                           int  b_ind, double  b_angle,
                           gdiam_real&   x,
                           gdiam_real&   y)
    {
        double  a_slope, b_slope, x_a, x_b, y_a, y_b;
 
        if (eq_real(a_angle, 0.0)  ||  eq_real(a_angle, PI))
        {
            x = ch[ b_ind ]->x;
            y = ch[ a_ind ]->y;
            return;
        }
 
        if (eq_real(a_angle, PI / 2.0)  ||  eq_real(a_angle, PI * 1.5))
        {
            x = ch[ a_ind ]->x;
            y = ch[ b_ind ]->y;
            return;
        }
 
        a_slope = tan(a_angle);
        b_slope = tan(b_angle);
        x_a = ch[ a_ind ]->x;
        y_a = ch[ a_ind ]->y;
        x_b = ch[ b_ind ]->x;
        y_b = ch[ b_ind ]->y;
 
        double  a_const, b_const, x_out, y_out;
 
        // Let see:
        //  l_a = y = a_slope * x + a_const ->
        //        y_a = a_slope * x_a + a_const  (since (x_a, y_a) in l_a
        //        a_const = y_a - a_slope * x_a
        a_const = y_a - a_slope * x_a;
        //printf( "a_const: %g, y_a: %g, a_slope: %g, x_a: %g\n",
        //        a_const, y_a, a_slope, x_a );
 
        b_const = y_b - b_slope * x_b;
        //printf( "b_const: %g, y_b: %g, b_slope: %g, x_b: %g\n",
        //        b_const, y_b, b_slope, x_b );
 
        x_out = - (a_const - b_const) / (a_slope - b_slope);
        y_out = a_slope * x_out + a_const;
        /*
        x = (y_b - y_a + x_a * a_slope - x_b * a_slope)
            / (a_slope - b_slope);
        y = ((y_a - x_a * a_slope) / a_slope - (y_b - x_b * b_slope)
             / b_slope) /
            (1.0 / a_slope - 1.0 / b_slope );
        printf( "gill: (%g, %g)\nMine: (%g, %g)\n",
                x, y, x_out, y_out );
        */
        //printf( "a_angle: %g, b_angle: %g\n",
        //        a_angle / PI, b_angle/ PI );
        //printf( "a_slope: %g, b_slope: %g\n",
        //        a_slope, b_slope );
        //printf( "on line_b: %g\n",
        //        b_slope * x_out + b_const - y_out );
        x = x_out;
        y = y_out;
        //printf( "line_a : y = %g * x + %g\n",
        //        a_slope, a_const );
        //printf( "line_b : y = %g * x + %g\n",
        //        b_slope, b_const );
        //printf( "a: (%g, %g)\n"
        //        "b: (%g, %g)\n",
        //        x_a, y_a, x_b, y_b );
 
        //printf( "out: (%g, %g)\n", x, y );
    }
 
 
    void  get_bbox(int   a_ind, gdiam_real   a_angle,
                   int   b_ind, gdiam_real   b_angle,
                   int   c_ind, gdiam_real   c_angle,
                   int   d_ind, gdiam_real   d_angle,
                   bbox_2d_info&    bbox,
                   gdiam_real&   area)
    {
        /*
        printf( "get_bbox: (%d, %g,\n"
                "%d, %g\n"
                "%d, %g\n"
                "%d, %g )\n",
                a_ind, a_angle / PI,
                b_ind, b_angle / PI,
                c_ind, c_angle / PI,
                d_ind, d_angle / PI );*/
 
        compute_crossing(a_ind, a_angle, b_ind, b_angle,
                         bbox.vertices[ 0 ][ 0 ],
                         bbox.vertices[ 0 ][ 1 ]);
        compute_crossing(b_ind, b_angle, c_ind, c_angle,
                         bbox.vertices[ 1 ][ 0 ],
                         bbox.vertices[ 1 ][ 1 ]);
        compute_crossing(c_ind, c_angle, d_ind, d_angle,
                         bbox.vertices[ 2 ][ 0 ],
                         bbox.vertices[ 2 ][ 1 ]);
        compute_crossing(d_ind, d_angle, a_ind, a_angle,
                         bbox.vertices[ 3 ][ 0 ],
                         bbox.vertices[ 3 ][ 1 ]);
 
        area = pnt_distance_2d(bbox.vertices[ 0 ],
                               bbox.vertices[ 1 ])
               * pnt_distance_2d(bbox.vertices[ 0 ],
                                 bbox.vertices[ 3 ]);
    }
 
    /* Given one angle (bbx direction), find the other three */
    void  get_angles(int  ind,
                     gdiam_real&    angle_1,
                     gdiam_real&    angle_2,
                     gdiam_real&    angle_3)
    {
        double   angle;
 
        angle = angles[ ind ];
        angle_1 = angle + PI * 0.5;
        //printf( "angle = %g\n", angle / PI );
        //printf( "angle1 = %g\n", angle_1 / PI );
 
        if (angle_1 >= (2.0 * PI))
        {
            angle_1 -= 2.0 * PI;
        }
 
        //printf( "angle1 = %g\n", angle_1 / PI );
 
 
        angle_2 = angle + PI;
 
        if (angle_2 >= (2.0 * PI))
        {
            angle_2 -= 2.0 * PI;
        }
 
        angle_3 = angle + 1.5 * PI;
 
        if (angle_3 >= (2.0 * PI))
        {
            angle_3 -= 2.0 * PI;
        }
    }
 
    void  compute_min_bbox_inner(bbox_2d_info&    min_bbox,
                                 gdiam_real&   min_area)
    {
        int        u, v, s, t;
        gdiam_real     ang1, ang2, ang3, tmp_area;
        bbox_2d_info  tmp_bbox;
 
        angles = (gdiam_real*)malloc(sizeof(gdiam_real)
                                     * (int)ch.size());
        assert(angles != NULL);
 
        /* Pre-computing all edge directions */
        for (u = 0; u < (int)ch.size(); u ++)
        {
            compute_edge_dir(u);
        }
 
        /* Initializing the search */
        //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
        get_angles(0, ang1, ang2, ang3);
 
        //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
        s = find_vertex(0, ang1);
        //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
        v = find_vertex(s, ang2);
        t = find_vertex(v, ang3);
 
        //printf( "_angles[ 0 ]: %g\n", angles[ 0 ] );
        get_bbox(0, angles[ 0 ],
                 s, ang1,
                 v, ang2,
                 t, ang3,
                 min_bbox, min_area);
        //min_bbox.dump();
 
        //printf( "min_area: %g\n", min_area );
        /* Performing a double rotating calipers */
        for (u = 1; u < (int)ch.size(); u ++)
        {
            get_angles(u, ang1, ang2, ang3);
            s = find_vertex(s, ang1);
            v = find_vertex(v, ang2);
            t = find_vertex(t, ang3);
            get_bbox(u, angles[ u ],
                     s, ang1,
                     v, ang2,
                     t, ang3,
                     tmp_bbox, tmp_area);
 
            //tmp_bbox.dump();
            //printf( "tmp_area: %g\n", tmp_area );
            if (min_area > tmp_area)
            {
                min_area = tmp_area;
                min_bbox = tmp_bbox;
            }
        }
 
        free(angles);
        angles = NULL;
    }
 
    void  compute_min_bbox(vec_point_2d&    points,
                           bbox_2d_info&    min_bbox,
                           gdiam_real&   min_area)
    {
        //printf( "before cunvex hull\n" );
        convex_hull(points, ch);
        //printf( "cunvex hull done!\n" );
        //dump_ch();
        compute_min_bbox_inner(min_bbox, min_area);
        //min_bbox.dump();
    }
 
};
 
void   dump_points(gdiam_point*   in_arr,
                   int  size)
{
    for (int  ind = 0; ind < size; ind++)
    {
        printf("%d: (%g, %g, %g)\n", ind,
               in_arr[ ind ][ 0 ],
               in_arr[ ind ][ 1 ],
               in_arr[ ind ][ 2 ]);
    }
}
 
#define  ZERO_EPS  (1e-6)
 
static void   construct_base_inner(gdiam_point  pnt1,
                                   gdiam_point  pnt2,
                                   gdiam_point  pnt3)
{
    pnt_normalize(pnt1);
    pnt_normalize(pnt2);
    pnt_normalize(pnt3);
 
    if ((pnt_length(pnt1) < ZERO_EPS)
        && (pnt_length(pnt2) < ZERO_EPS)
        && (pnt_length(pnt3) < ZERO_EPS))
    {
        assert(false);
    }
 
    if ((pnt_length(pnt1) < ZERO_EPS)
        && (pnt_length(pnt2) < ZERO_EPS))
    {
        gdiam_generate_orthonormal_base(pnt3, pnt1, pnt2);
        return;
    }
 
    if ((pnt_length(pnt1) < ZERO_EPS)
        && (pnt_length(pnt3) < ZERO_EPS))
    {
        gdiam_generate_orthonormal_base(pnt2, pnt1, pnt3);
        return;
    }
 
    if ((pnt_length(pnt2) < ZERO_EPS)
        && (pnt_length(pnt3) < ZERO_EPS))
    {
        gdiam_generate_orthonormal_base(pnt1, pnt2, pnt3);
        return;
    }
 
    if (pnt_length(pnt1) < ZERO_EPS)
    {
        pnt_cross_prod(pnt2, pnt3, pnt1);
        return;
    }
 
    if (pnt_length(pnt2) < ZERO_EPS)
    {
        pnt_cross_prod(pnt1, pnt3, pnt2);
        return;
    }
 
    if (pnt_length(pnt3) < ZERO_EPS)
    {
        pnt_cross_prod(pnt1, pnt2, pnt3);
        return;
    }
}
 
void   construct_base(gdiam_point  pnt1,
                      gdiam_point  pnt2,
                      gdiam_point  pnt3)
{
    construct_base_inner(pnt1, pnt2, pnt3);
    pnt_scale_and_add(pnt2, -pnt_dot_prod(pnt1, pnt2),
                      pnt1);
    pnt_normalize(pnt2);
 
    pnt_scale_and_add(pnt3, -pnt_dot_prod(pnt1, pnt3),
                      pnt1);
    pnt_scale_and_add(pnt3, -pnt_dot_prod(pnt2, pnt3),
                      pnt2);
    pnt_normalize(pnt3);
}
 
 
class  ProjPointSet
{
private:
    gdiam_point_t  base_x, base_y, base_proj;
    point2d*   arr;
    vec_point_2d   points, ch;
    gdiam_bbox   bbox;
    gdiam_point*   in_arr;
    int  size;
 
public:
    ~ProjPointSet()
    {
        term();
    }
    void  init(gdiam_point  dir,
               gdiam_point*   _in_arr,
               int  _size)
    {
        arr = NULL;
 
        if (pnt_length(dir) == 0.0)
        {
            dump_points(_in_arr, _size);
            pnt_dump(dir);
            fflush(stdout);
            fflush(stderr);
            assert(pnt_length(dir) > 0.0);
        }
 
        size = _size;
        in_arr = _in_arr;
        assert(size > 0);
 
        pnt_copy(base_proj, dir);
        gdiam_generate_orthonormal_base(base_proj, base_x, base_y);
 
        arr = (point2d*)malloc(sizeof(point2d) * size);
        assert(arr != 0);
 
        for (int  ind = 0; ind < size; ind++)
        {
            arr[ ind ].init(in_arr[ ind ], base_x, base_y);
            points.push_back(&(arr[ ind ]));
        }
    }
 
    void  compute()
    {
        MinAreaRectangle  mar;
        bbox_2d_info   min_bbox;
        gdiam_real   min_area;
        double  x1, y1, x2, y2;
        gdiam_point_t  out_1, out_2;
 
        mar.compute_min_bbox(points, min_bbox, min_area);
 
        // We take the resulting min_area rectangle and lift it back to 3d...
        x1 = min_bbox.vertices[ 0 ][ 0 ] -
             min_bbox.vertices[ 1 ][ 0 ];
        y1 = min_bbox.vertices[ 0 ][ 1 ] -
             min_bbox.vertices[ 1 ][ 1 ];
        x2 = min_bbox.vertices[ 0 ][ 0 ] -
             min_bbox.vertices[ 3 ][ 0 ];
        y2 = min_bbox.vertices[ 0 ][ 1 ] -
             min_bbox.vertices[ 3 ][ 1 ];
 
        //printf( "prod: %g\n", x1 * x2 + y1 * y2 );
 
        double  len;
 
        len = sqrt(x1 * x1 + y1 * y1);
 
        if (len > 0.0)
        {
            x1 /= len;
            y1 /= len;
        }
 
        len = sqrt(x2 * x2 + y2 * y2);
 
        if (len > 0.0)
        {
            x2 /= len;
            y2 /= len;
        }
 
        pnt_zero(out_1);
        pnt_scale_and_add(out_1, x1, base_x);
        pnt_scale_and_add(out_1, y1, base_y);
        pnt_normalize(out_1);
 
        pnt_zero(out_2);
        pnt_scale_and_add(out_2, x2, base_x);
        pnt_scale_and_add(out_2, y2, base_y);
        pnt_normalize(out_2);
 
        construct_base(base_proj, out_1, out_2);
 
        if ((!(fabs(pnt_dot_prod(base_proj, out_1)) < 1e-6))
            || (!(fabs(pnt_dot_prod(base_proj, out_2)) < 1e-6))
            || (!(fabs(pnt_dot_prod(out_1, out_2)) < 1e-6)))
        {
            printf("should be all close to zero: %g, %g, %g\n",
                   pnt_dot_prod(base_proj, out_1),
                   pnt_dot_prod(base_proj, out_2),
                   pnt_dot_prod(out_1, out_2));
            pnt_dump(base_proj);
            pnt_dump(out_1);
            pnt_dump(out_2);
 
 
            printf("real points:\n");
            dump_points(in_arr, size);
 
            printf("points on CH:\n");
            dump(points);
 
            printf("BBox:\n");
            min_bbox.dump();
 
            fflush(stdout);
            fflush(stderr);
            assert(fabs(pnt_dot_prod(base_proj, out_1)) < 1e-6);
            assert(fabs(pnt_dot_prod(base_proj, out_2)) < 1e-6);
            assert(fabs(pnt_dot_prod(out_1, out_2)) < 1e-6);
        }
 
        bbox.init(base_proj, out_1, out_2);
 
        for (int  ind = 0; ind < size; ind++)
        {
            bbox.bound(in_arr[ ind ]);
        }
    }
 
    gdiam_bbox&    get_bbox()
    {
        return  bbox;
    }
 
    void  term()
    {
        if (arr != NULL)
        {
            free(arr);
            arr = NULL;
        }
    }
    void  init()
    {
    }
};
 
 
gdiam_bbox   gdiam_mvbb_optimize(gdiam_point*   start, int  size,
                                 gdiam_bbox  bb_out, int  times)
{
    gdiam_bbox  bb_tmp;
 
    //printf( "gdiam_mvbb_optimize called\n" );
 
    for (int  ind = 0; ind < times; ind++)
    {
        ProjPointSet  pps;
 
        if (pnt_length(bb_out.get_dir(ind % 3)) == 0.0)
        {
            printf("Dumping!\n");
            bb_out.dump();
            fflush(stdout);
            continue;
        }
 
        pps.init(bb_out.get_dir(ind % 3), start, size);
 
        pps.compute();
 
        //printf( "Old volume: %g\n", bb_out.volume() );
        bb_tmp = pps.get_bbox();
        //printf( "New volume: %g\n", bb_tmp.volume() );
        //printf( "Old bounding box:\n" );
        //bb2.dump();
 
        if (bb_tmp.volume() < bb_out.volume())
        {
            bb_out = bb_tmp;
        }
    }
 
    return  bb_out;
}
 
 
gdiam_bbox   gdiam_approx_mvbb(gdiam_point*   start, int  size,
                               gdiam_real  eps)
{
    gdiam_bbox  bb, bb2;
 
    bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
    bb2 = gdiam_mvbb_optimize(start, size,  bb, 10);
 
    //bb2.dump();
 
    return  bb2;
}
 
 
static int  gcd2(int  a, int  b)
{
    if (a == 0  ||  a == b)
    {
        return  b;
    }
 
    if (b == 0)
    {
        return  a;
    }
 
    if (a > b)
    {
        return  gcd2(a % b, b);
    }
    else
    {
        return  gcd2(a, b % a);
    }
}
 
 
static int  gcd3(int  a, int  b, int  c)
{
    a = abs(a);
    b = abs(b);
    c = abs(c);
 
    if (a == 0)
    {
        return  gcd2(b, c);
    }
 
    if (b == 0)
    {
        return  gcd2(a, c);
    }
 
    if (c == 0)
    {
        return  gcd2(a, b);
    }
 
    return  gcd2(a, gcd2(b, c));
}
 
 
static void  try_direction(gdiam_bbox&   bb,
                           gdiam_bbox&   bb_out,
                           gdiam_point*   start,
                           int  size,
                           int  x_coef,
                           int  y_coef,
                           int  z_coef)
{
    gdiam_bbox  bb_new;
 
    if (gcd3(x_coef, y_coef, z_coef) > 1)
    {
        return;
    }
 
    if ((x_coef == 0)  && (y_coef == 0)
        && (z_coef == 0))
    {
        return;
    }
 
    //printf( "trying: (%d, %d, %d)\n",
    //        x_coef, y_coef, z_coef );
    //fflush( stdout );
    gdiam_point_t  new_dir;
    bb.combine(new_dir, x_coef, y_coef, z_coef);
 
    ProjPointSet  pps;
 
    pps.init(new_dir, start, size);
 
    pps.compute();
 
    bb_new = pps.get_bbox();
    bb_new = gdiam_mvbb_optimize(start, size, bb_new, 10);
 
    if (bb_new.volume() < bb_out.volume())
    {
        bb_out = bb_new;
    }
}
 
 
gdiam_bbox   gdiam_approx_mvbb_grid(gdiam_point*   start, int  size,
                                    int  grid_size)
{
    gdiam_bbox  bb, bb_out;
 
    bb_out = bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
 
    if (bb.volume() == 0)
    {
        dump_points(start, size);
        printf("1zero volume???\n");
        bb.dump();
    }
 
    bb_out = bb = gdiam_mvbb_optimize(start, size, bb_out, 10);
 
    if (bb.volume() == 0)
    {
        printf("2zero volume???\n");
        bb.dump();
    }
 
    //trying: (1, -4, 4)
    //try_direction( bb, bb_out, start, size,
    //               1, -4, 4 );
 
    for (int  x_coef = -grid_size; x_coef <= grid_size; x_coef++)
        for (int  y_coef = -grid_size; y_coef <= grid_size; y_coef++)
            for (int  z_coef = 0; z_coef <= grid_size; z_coef++)
                try_direction(bb, bb_out, start, size,
                              x_coef, y_coef, z_coef);
 
    bb_out = gdiam_mvbb_optimize(start, size, bb_out, 10);
 
    return  bb_out;
}
 
 
static void   register_point(gdiam_point  pnt,
                             gdiam_point*   tops,
                             gdiam_point*   bottoms,
                             int  grid_size,
                             gdiam_bbox&   bb)
{
    gdiam_point_t  pnt_bb, pnt_bottom, pnt_top;
    int  x_ind, y_ind, position;
 
    bb.get_normalized_coordinates(pnt, pnt_bb);
 
    // x_ind
    x_ind = (int)(pnt_bb[ 0 ] * grid_size);
    assert((-1 <= x_ind)  && (x_ind <= grid_size));
 
    if (x_ind < 0)
    {
        x_ind = 0;
    }
 
    if (x_ind >= grid_size)
    {
        x_ind = grid_size - 1;
    }
 
    // y_ind
    y_ind = (int)(pnt_bb[ 1 ] * grid_size);
    assert((-1 <= y_ind)  && (y_ind <= grid_size));
 
    if (y_ind < 0)
    {
        y_ind = 0;
    }
 
    if (y_ind >= grid_size)
    {
        y_ind = grid_size - 1;
    }
 
    position = x_ind + y_ind * grid_size;
 
    if (tops[ position ] == NULL)
    {
        tops[ position ] = bottoms[ position ] = pnt;
        return;
    }
 
    bb.get_normalized_coordinates(tops[ position ], pnt_top);
 
    if (pnt_top[ 2 ] < pnt_bb[ 2 ])
    {
        tops[ position ] = pnt;
    }
 
    bb.get_normalized_coordinates(bottoms[ position ], pnt_bottom);
 
    if (pnt_bottom[ 2 ] > pnt_bb[ 2 ])
    {
        bottoms[ position ] = pnt;
    }
}
 
 
 
// gdiam_convex_sample
 
//   We are given a point set, and (hopefully) a tight fitting
// bounding box. We compute a smplae of the given size that represents
// the point-set. The only guarenteed is that if we use ample of size m,
// we get an approximation of quality about 1/\sqrt{m}. Note that we pad
// the sample if necessary to get the desired size.
gdiam_point*   gdiam_convex_sample(gdiam_point*   start, int  size,
                                   gdiam_bbox&   bb,
                                   int  sample_size)
{
    int  grid_size, mem_size, grid_entries;
    gdiam_point*   bottoms, * tops, * out_arr;
    int  out_count;
 
    assert(sample_size > 1);
 
    // we want the grid to be on the two minor dimensions.
    bb.set_third_dim_longest();
 
    grid_size = (int)sqrt(sample_size / 2);
 
    grid_entries = grid_size * grid_size;
    mem_size = (int)(sizeof(gdiam_point) * grid_size * grid_size);
 
    bottoms = (gdiam_point*)malloc(mem_size);
    tops = (gdiam_point*)malloc(mem_size);
    out_arr = (gdiam_point*)malloc(sizeof(gdiam_point) * sample_size);
 
    assert(bottoms != NULL);
    assert(tops != NULL);
    assert(out_arr != NULL);
 
    for (int  ind = 0; ind < grid_entries; ind++)
    {
        tops[ ind ] = bottoms[ ind ] = NULL;
    }
 
    // No we stream the points registering them with the relevant
    // shaft in the grid.
    for (int  ind = 0; ind < size; ind++)
        register_point(start[ ind ], tops, bottoms,
                       grid_size, bb);
 
    out_count = 0;
 
    for (int  ind = 0; ind < grid_entries; ind++)
    {
        if (tops[ ind ] == NULL)
        {
            continue;
        }
 
        out_arr[ out_count++ ] = tops[ ind ];
 
        if (tops[ ind ] != bottoms[ ind ])
        {
            out_arr[ out_count++ ] = bottoms[ ind ];
        }
    }
 
    // We dont have neough entries in our sample - lets randomly pick
    // points.
    while (out_count < sample_size)
    {
        out_arr[ out_count++ ] = start[ rand() % size ];
    }
 
    free(tops);
    free(bottoms);
 
    return  out_arr;
}
 
 
gdiam_bbox   gdiam_approx_mvbb_grid_sample(gdiam_point*   start, int  size,
        int  grid_size, int  sample_size)
{
    gdiam_bbox  bb, bb_new;
    gdiam_point*   sample;
 
    if (sample_size >= size)
    {
        return   gdiam_approx_mvbb_grid(start, size, grid_size);
    }
 
    //printf( "gdiam_approx_mvbb_grid_sample(...) called\n" );
    //fflush( stdout );
 
    bb = gdiam_approx_const_mvbb(start, size, 0.0, NULL);
    //printf( "guess from const approximation:\n" );
 
    sample = gdiam_convex_sample(start, size, bb, sample_size);
 
    bb_new = gdiam_approx_mvbb_grid(sample, sample_size, grid_size);
 
    //printf( "new volume on sample: %g\n",
    //        bb_new.volume() );
    //fflush( stdout );
 
    for (int  ind = 0; ind < size; ind++)
    {
        bb_new.bound(start[ ind ]);
    }
 
    //printf( "new volume on resized input: %g\n",
    //        bb_new.volume() );
 
    //bb_new = gdiam_mvbb_optimize( start, size, bb_new, 10 );
    //printf(s "optimized new volume: %g\n",
    //        bb_new.volume() );
 
    return  bb_new;
}
 
 
 
gdiam_bbox   gdiam_approx_mvbb_grid_sample(gdiam_real*   start, int  size,
        int  grid_size, int  sample_size)
{
    gdiam_point*   p_arr;
    gdiam_bbox  bb;
 
    //printf( "gdiam_approx_mvbb_grid_sample( ?, %d, ...)\n", size );
    //fflush( stdout );
 
    p_arr = gdiam_convert(start, size);
    bb = gdiam_approx_mvbb_grid_sample(p_arr, size, grid_size, sample_size);
    free(p_arr);
 
    return  bb;
}
 
 
/* gdiam.C - End of File ------------------------------------------*/