Logo Search packages:      
Sourcecode: madlib version File versions  Download package

kd_util.cpp

//----------------------------------------------------------------------
// File:                kd_util.cpp
// Programmer:          Sunil Arya and David Mount
// Description:         Common utilities for kd-trees
// Last modified: 01/04/05 (Version 1.0)
//----------------------------------------------------------------------
// Copyright (c) 1997-2005 University of Maryland and Sunil Arya and
// David Mount.  All Rights Reserved.
// 
// This software and related documentation is part of the Approximate
// Nearest Neighbor Library (ANN).  This software is provided under
// the provisions of the Lesser GNU Public License (LGPL).  See the
// file ../ReadMe.txt for further information.
// 
// The University of Maryland (U.M.) and the authors make no
// representations about the suitability or fitness of this software for
// any purpose.  It is provided "as is" without express or implied
// warranty.
//----------------------------------------------------------------------
// History:
//    Revision 0.1  03/04/98
//          Initial release
//----------------------------------------------------------------------

#include "kd_util.h"                            // kd-utility declarations

#include <ANN/ANNperf.h>                        // performance evaluation

//----------------------------------------------------------------------
// The following routines are utility functions for manipulating
// points sets, used in determining splitting planes for kd-tree
// construction.
//----------------------------------------------------------------------

//----------------------------------------------------------------------
//    NOTE: Virtually all point indexing is done through an index (i.e.
//    permutation) array pidx.  Consequently, a reference to the d-th
//    coordinate of the i-th point is pa[pidx[i]][d].  The macro PA(i,d)
//    is a shorthand for this.
//----------------------------------------------------------------------
                                                            // standard 2-d indirect indexing
#define PA(i,d)               (pa[pidx[(i)]][(d)])
                                                            // accessing a single point
#define PP(i)                 (pa[pidx[(i)]])

//----------------------------------------------------------------------
//    annAspectRatio
//          Compute the aspect ratio (ratio of longest to shortest side)
//          of a rectangle.
//----------------------------------------------------------------------

double annAspectRatio(
      int                           dim,              // dimension
      const ANNorthRect &bnd_box)         // bounding cube
{
      ANNcoord length = bnd_box.hi[0] - bnd_box.lo[0];
      ANNcoord min_length = length;                   // min side length
      ANNcoord max_length = length;                   // max side length
      for (int d = 0; d < dim; d++) {
            length = bnd_box.hi[d] - bnd_box.lo[d];
            if (length < min_length) min_length = length;
            if (length > max_length) max_length = length;
      }
      return max_length/min_length;
}

//----------------------------------------------------------------------
//    annEnclRect, annEnclCube
//          These utilities compute the smallest rectangle and cube enclosing
//          a set of points, respectively.
//----------------------------------------------------------------------

void annEnclRect(
      ANNpointArray           pa,                     // point array
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           dim,              // dimension
      ANNorthRect             &bnds)                  // bounding cube (returned)
{
      for (int d = 0; d < dim; d++) {           // find smallest enclosing rectangle
            ANNcoord lo_bnd = PA(0,d);          // lower bound on dimension d
            ANNcoord hi_bnd = PA(0,d);          // upper bound on dimension d
            for (int i = 0; i < n; i++) {
                  if (PA(i,d) < lo_bnd) lo_bnd = PA(i,d);
                  else if (PA(i,d) > hi_bnd) hi_bnd = PA(i,d);
            }
            bnds.lo[d] = lo_bnd;
            bnds.hi[d] = hi_bnd;
      }
}

void annEnclCube(                               // compute smallest enclosing cube
      ANNpointArray           pa,                     // point array
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           dim,              // dimension
      ANNorthRect             &bnds)                  // bounding cube (returned)
{
      int d;
                                                            // compute smallest enclosing rect
      annEnclRect(pa, pidx, n, dim, bnds);

      ANNcoord max_len = 0;                     // max length of any side
      for (d = 0; d < dim; d++) {               // determine max side length
            ANNcoord len = bnds.hi[d] - bnds.lo[d];
            if (len > max_len) {                // update max_len if longest
                  max_len = len;
            }
      }
      for (d = 0; d < dim; d++) {               // grow sides to match max
            ANNcoord len = bnds.hi[d] - bnds.lo[d];
            ANNcoord half_diff = (max_len - len) / 2;
            bnds.lo[d] -= half_diff;
            bnds.hi[d] += half_diff;
      }
}

//----------------------------------------------------------------------
//    annBoxDistance - utility routine which computes distance from point to
//          box (Note: most distances to boxes are computed using incremental
//          distance updates, not this function.)
//----------------------------------------------------------------------

ANNdist annBoxDistance(             // compute distance from point to box
      const ANNpoint          q,                      // the point
      const ANNpoint          lo,                     // low point of box
      const ANNpoint          hi,                     // high point of box
      int                           dim)              // dimension of space
{
      register ANNdist dist = 0.0;        // sum of squared distances
      register ANNdist t;

      for (register int d = 0; d < dim; d++) {
            if (q[d] < lo[d]) {                       // q is left of box
                  t = ANNdist(lo[d]) - ANNdist(q[d]);
                  dist = ANN_SUM(dist, ANN_POW(t));
            }
            else if (q[d] > hi[d]) {            // q is right of box
                  t = ANNdist(q[d]) - ANNdist(hi[d]);
                  dist = ANN_SUM(dist, ANN_POW(t));
            }
      }
      ANN_FLOP(4*dim)                                 // increment floating op count

      return dist;
}

//----------------------------------------------------------------------
//    annSpread - find spread along given dimension
//    annMinMax - find min and max coordinates along given dimension
//    annMaxSpread - find dimension of max spread
//----------------------------------------------------------------------

ANNcoord annSpread(                       // compute point spread along dimension
      ANNpointArray           pa,                     // point array
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           d)                      // dimension to check
{
      ANNcoord min = PA(0,d);                   // compute max and min coords
      ANNcoord max = PA(0,d);
      for (int i = 1; i < n; i++) {
            ANNcoord c = PA(i,d);
            if (c < min) min = c;
            else if (c > max) max = c;
      }
      return (max - min);                             // total spread is difference
}

void annMinMax(                           // compute min and max coordinates along dim
      ANNpointArray           pa,                     // point array
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           d,                      // dimension to check
      ANNcoord                &min,             // minimum value (returned)
      ANNcoord                &max)             // maximum value (returned)
{
      min = PA(0,d);                                  // compute max and min coords
      max = PA(0,d);
      for (int i = 1; i < n; i++) {
            ANNcoord c = PA(i,d);
            if (c < min) min = c;
            else if (c > max) max = c;
      }
}

int annMaxSpread(                               // compute dimension of max spread
      ANNpointArray           pa,                     // point array
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           dim)              // dimension of space
{
      int max_dim = 0;                          // dimension of max spread
      ANNcoord max_spr = 0;                     // amount of max spread

      if (n == 0) return max_dim;               // no points, who cares?

      for (int d = 0; d < dim; d++) {           // compute spread along each dim
            ANNcoord spr = annSpread(pa, pidx, n, d);
            if (spr > max_spr) {                // bigger than current max
                  max_spr = spr;
                  max_dim = d;
            }
      }
      return max_dim;
}

//----------------------------------------------------------------------
//    annMedianSplit - split point array about its median
//          Splits a subarray of points pa[0..n] about an element of given
//          rank (median: n_lo = n/2) with respect to dimension d.  It places
//          the element of rank n_lo-1 correctly (because our splitting rule
//          takes the mean of these two).  On exit, the array is permuted so
//          that:
//
//          pa[0..n_lo-2][d] <= pa[n_lo-1][d] <= pa[n_lo][d] <= pa[n_lo+1..n-1][d].
//
//          The mean of pa[n_lo-1][d] and pa[n_lo][d] is returned as the
//          splitting value.
//
//          All indexing is done indirectly through the index array pidx.
//
//          This function uses the well known selection algorithm due to
//          C.A.R. Hoare.
//----------------------------------------------------------------------

                                                            // swap two points in pa array
#define PASWAP(a,b) { int tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }

void annMedianSplit(
      ANNpointArray           pa,                     // points to split
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           d,                      // dimension along which to split
      ANNcoord                &cv,              // cutting value
      int                           n_lo)             // split into n_lo and n-n_lo
{
      int l = 0;                                      // left end of current subarray
      int r = n-1;                                    // right end of current subarray
      while (l < r) {
            register int i = (r+l)/2;           // select middle as pivot
            register int k;

            if (PA(i,d) > PA(r,d))              // make sure last > pivot
                  PASWAP(i,r)
            PASWAP(l,i);                              // move pivot to first position

            ANNcoord c = PA(l,d);               // pivot value
            i = l;
            k = r;
            for(;;) {                                 // pivot about c
                  while (PA(++i,d) < c) ;
                  while (PA(--k,d) > c) ;
                  if (i < k) PASWAP(i,k) else break;
            }
            PASWAP(l,k);                              // pivot winds up in location k

            if (k > n_lo)        r = k-1;       // recurse on proper subarray
            else if (k < n_lo) l = k+1;
            else break;                               // got the median exactly
      }
      if (n_lo > 0) {                                 // search for next smaller item
            ANNcoord c = PA(0,d);               // candidate for max
            int k = 0;                                // candidate's index
            for (int i = 1; i < n_lo; i++) {
                  if (PA(i,d) > c) {
                        c = PA(i,d);
                        k = i;
                  }
            }
            PASWAP(n_lo-1, k);                        // max among pa[0..n_lo-1] to pa[n_lo-1]
      }
                                                            // cut value is midpoint value
      cv = (PA(n_lo-1,d) + PA(n_lo,d))/2.0;
}

//----------------------------------------------------------------------
//    annPlaneSplit - split point array about a cutting plane
//          Split the points in an array about a given plane along a
//          given cutting dimension.  On exit, br1 and br2 are set so
//          that:
//          
//                      pa[ 0 ..br1-1] <  cv
//                      pa[br1..br2-1] == cv
//                      pa[br2.. n -1] >  cv
//
//          All indexing is done indirectly through the index array pidx.
//
//----------------------------------------------------------------------

void annPlaneSplit(                       // split points by a plane
      ANNpointArray           pa,                     // points to split
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           d,                      // dimension along which to split
      ANNcoord                cv,                     // cutting value
      int                           &br1,             // first break (values < cv)
      int                           &br2)             // second break (values == cv)
{
      int l = 0;
      int r = n-1;
      for(;;) {                                       // partition pa[0..n-1] about cv
            while (l < n && PA(l,d) < cv) l++;
            while (r >= 0 && PA(r,d) >= cv) r--;
            if (l > r) break;
            PASWAP(l,r);
            l++; r--;
      }
      br1 = l;                            // now: pa[0..br1-1] < cv <= pa[br1..n-1]
      r = n-1;
      for(;;) {                                       // partition pa[br1..n-1] about cv
            while (l < n && PA(l,d) <= cv) l++;
            while (r >= br1 && PA(r,d) > cv) r--;
            if (l > r) break;
            PASWAP(l,r);
            l++; r--;
      }
      br2 = l;                            // now: pa[br1..br2-1] == cv < pa[br2..n-1]
}


//----------------------------------------------------------------------
//    annBoxSplit - split point array about a orthogonal rectangle
//          Split the points in an array about a given orthogonal
//          rectangle.  On exit, n_in is set to the number of points
//          that are inside (or on the boundary of) the rectangle.
//
//          All indexing is done indirectly through the index array pidx.
//
//----------------------------------------------------------------------

void annBoxSplit(                   // split points by a box
      ANNpointArray           pa,                     // points to split
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           dim,              // dimension of space
      ANNorthRect             &box,             // the box
      int                           &n_in)                  // number of points inside (returned)
{
      int l = 0;
      int r = n-1;
      for(;;) {                                       // partition pa[0..n-1] about box
            while (l < n && box.inside(dim, PP(l))) l++;
            while (r >= 0 && !box.inside(dim, PP(r))) r--;
            if (l > r) break;
            PASWAP(l,r);
            l++; r--;
      }
      n_in = l;                           // now: pa[0..n_in-1] inside and rest outside
}

//----------------------------------------------------------------------
//    annSplitBalance - compute balance factor for a given plane split
//          Balance factor is defined as the number of points lying
//          below the splitting value minus n/2 (median).  Thus, a
//          median split has balance 0, left of this is negative and
//          right of this is positive.  (The points are unchanged.)
//----------------------------------------------------------------------

int annSplitBalance(                // determine balance factor of a split
      ANNpointArray           pa,                     // points to split
      ANNidxArray             pidx,             // point indices
      int                           n,                      // number of points
      int                           d,                      // dimension along which to split
      ANNcoord                cv)                     // cutting value
{
      int n_lo = 0;
      for(int i = 0; i < n; i++) {        // count number less than cv
            if (PA(i,d) < cv) n_lo++;
      }
      return n_lo - n/2;
}

//----------------------------------------------------------------------
//    annBox2Bnds - convert bounding box to list of bounds
//          Given two boxes, an inner box enclosed within a bounding
//          box, this routine determines all the sides for which the
//          inner box is strictly contained with the bounding box,
//          and adds an appropriate entry to a list of bounds.  Then
//          we allocate storage for the final list of bounds, and return
//          the resulting list and its size.
//----------------------------------------------------------------------

void annBox2Bnds(                               // convert inner box to bounds
      const ANNorthRect &inner_box,       // inner box
      const ANNorthRect &bnd_box,         // enclosing box
      int                           dim,              // dimension of space
      int                           &n_bnds,          // number of bounds (returned)
      ANNorthHSArray          &bnds)                  // bounds array (returned)
{
      int i;
      n_bnds = 0;                                                 // count number of bounds
      for (i = 0; i < dim; i++) {
            if (inner_box.lo[i] > bnd_box.lo[i])      // low bound is inside
                        n_bnds++;
            if (inner_box.hi[i] < bnd_box.hi[i])      // high bound is inside
                        n_bnds++;
      }

      bnds = new ANNorthHalfSpace[n_bnds];            // allocate appropriate size

      int j = 0;
      for (i = 0; i < dim; i++) {                           // fill the array
            if (inner_box.lo[i] > bnd_box.lo[i]) {
                        bnds[j].cd = i;
                        bnds[j].cv = inner_box.lo[i];
                        bnds[j].sd = +1;
                        j++;
            }
            if (inner_box.hi[i] < bnd_box.hi[i]) {
                        bnds[j].cd = i;
                        bnds[j].cv = inner_box.hi[i];
                        bnds[j].sd = -1;
                        j++;
            }
      }
}

//----------------------------------------------------------------------
//    annBnds2Box - convert list of bounds to bounding box
//          Given an enclosing box and a list of bounds, this routine
//          computes the corresponding inner box.  It is assumed that
//          the box points have been allocated already.
//----------------------------------------------------------------------

void annBnds2Box(
      const ANNorthRect &bnd_box,         // enclosing box
      int                           dim,              // dimension of space
      int                           n_bnds,                 // number of bounds
      ANNorthHSArray          bnds,             // bounds array
      ANNorthRect             &inner_box)       // inner box (returned)
{
      annAssignRect(dim, inner_box, bnd_box);         // copy bounding box to inner

      for (int i = 0; i < n_bnds; i++) {
            bnds[i].project(inner_box.lo);                  // project each endpoint
            bnds[i].project(inner_box.hi);
      }
}

Generated by  Doxygen 1.6.0   Back to index