Skip to content
Snippets Groups Projects
Commit e2669093 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent 509e02db
No related branches found
No related tags found
No related merge requests found
Showing with 0 additions and 4606 deletions
//----------------------------------------------------------------------
// File: kd_split.cpp
// Programmer: Sunil Arya and David Mount
// Description: Methods for splitting 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
// Revision 1.0 04/01/05
//----------------------------------------------------------------------
#include "kd_tree.h" // kd-tree definitions
#include "kd_util.h" // kd-tree utilities
#include "kd_split.h" // splitting functions
//----------------------------------------------------------------------
// Constants
//----------------------------------------------------------------------
const double ERR = 0.001; // a small value
const double FS_ASPECT_RATIO = 3.0; // maximum allowed aspect ratio
// in fair split. Must be >= 2.
//----------------------------------------------------------------------
// kd_split - Bentley's standard splitting routine for kd-trees
// Find the dimension of the greatest spread, and split
// just before the median point along this dimension.
//----------------------------------------------------------------------
void kd_split(
ANNpointArray pa, // point array (permuted on return)
ANNidxArray pidx, // point indices
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo) // num of points on low side (returned)
{
// find dimension of maximum spread
cut_dim = annMaxSpread(pa, pidx, n, dim);
n_lo = n/2; // median rank
// split about median
annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
}
//----------------------------------------------------------------------
// midpt_split - midpoint splitting rule for box-decomposition trees
//
// This is the simplest splitting rule that guarantees boxes
// of bounded aspect ratio. It simply cuts the box with the
// longest side through its midpoint. If there are ties, it
// selects the dimension with the maximum point spread.
//
// WARNING: This routine (while simple) doesn't seem to work
// well in practice in high dimensions, because it tends to
// generate a large number of trivial and/or unbalanced splits.
// Either kd_split(), sl_midpt_split(), or fair_split() are
// recommended, instead.
//----------------------------------------------------------------------
void midpt_split(
ANNpointArray pa, // point array
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo) // num of points on low side (returned)
{
int d;
ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
for (d = 1; d < dim; d++) { // find length of longest box side
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (length > max_length) {
max_length = length;
}
}
ANNcoord max_spread = -1; // find long side with most spread
for (d = 0; d < dim; d++) {
// is it among longest?
if (double(bnds.hi[d] - bnds.lo[d]) >= (1-ERR)*max_length) {
// compute its spread
ANNcoord spr = annSpread(pa, pidx, n, d);
if (spr > max_spread) { // is it max so far?
max_spread = spr;
cut_dim = d;
}
}
}
// split along cut_dim at midpoint
cut_val = (bnds.lo[cut_dim] + bnds.hi[cut_dim]) / 2;
// permute points accordingly
int br1, br2;
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
//------------------------------------------------------------------
// On return: pa[0..br1-1] < cut_val
// pa[br1..br2-1] == cut_val
// pa[br2..n-1] > cut_val
//
// We can set n_lo to any value in the range [br1..br2].
// We choose split so that points are most evenly divided.
//------------------------------------------------------------------
if (br1 > n/2) n_lo = br1;
else if (br2 < n/2) n_lo = br2;
else n_lo = n/2;
}
//----------------------------------------------------------------------
// sl_midpt_split - sliding midpoint splitting rule
//
// This is a modification of midpt_split, which has the nonsensical
// name "sliding midpoint". The idea is that we try to use the
// midpoint rule, by bisecting the longest side. If there are
// ties, the dimension with the maximum spread is selected. If,
// however, the midpoint split produces a trivial split (no points
// on one side of the splitting plane) then we slide the splitting
// (maintaining its orientation) until it produces a nontrivial
// split. For example, if the splitting plane is along the x-axis,
// and all the data points have x-coordinate less than the x-bisector,
// then the split is taken along the maximum x-coordinate of the
// data points.
//
// Intuitively, this rule cannot generate trivial splits, and
// hence avoids midpt_split's tendency to produce trees with
// a very large number of nodes.
//
//----------------------------------------------------------------------
void sl_midpt_split(
ANNpointArray pa, // point array
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo) // num of points on low side (returned)
{
int d;
ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
for (d = 1; d < dim; d++) { // find length of longest box side
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (length > max_length) {
max_length = length;
}
}
ANNcoord max_spread = -1; // find long side with most spread
for (d = 0; d < dim; d++) {
// is it among longest?
if ((bnds.hi[d] - bnds.lo[d]) >= (1-ERR)*max_length) {
// compute its spread
ANNcoord spr = annSpread(pa, pidx, n, d);
if (spr > max_spread) { // is it max so far?
max_spread = spr;
cut_dim = d;
}
}
}
// ideal split at midpoint
ANNcoord ideal_cut_val = (bnds.lo[cut_dim] + bnds.hi[cut_dim])/2;
ANNcoord min, max;
annMinMax(pa, pidx, n, cut_dim, min, max); // find min/max coordinates
if (ideal_cut_val < min) // slide to min or max as needed
cut_val = min;
else if (ideal_cut_val > max)
cut_val = max;
else
cut_val = ideal_cut_val;
// permute points accordingly
int br1, br2;
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
//------------------------------------------------------------------
// On return: pa[0..br1-1] < cut_val
// pa[br1..br2-1] == cut_val
// pa[br2..n-1] > cut_val
//
// We can set n_lo to any value in the range [br1..br2] to satisfy
// the exit conditions of the procedure.
//
// if ideal_cut_val < min (implying br2 >= 1),
// then we select n_lo = 1 (so there is one point on left) and
// if ideal_cut_val > max (implying br1 <= n-1),
// then we select n_lo = n-1 (so there is one point on right).
// Otherwise, we select n_lo as close to n/2 as possible within
// [br1..br2].
//------------------------------------------------------------------
if (ideal_cut_val < min) n_lo = 1;
else if (ideal_cut_val > max) n_lo = n-1;
else if (br1 > n/2) n_lo = br1;
else if (br2 < n/2) n_lo = br2;
else n_lo = n/2;
}
//----------------------------------------------------------------------
// fair_split - fair-split splitting rule
//
// This is a compromise between the kd-tree splitting rule (which
// always splits data points at their median) and the midpoint
// splitting rule (which always splits a box through its center.
// The goal of this procedure is to achieve both nicely balanced
// splits, and boxes of bounded aspect ratio.
//
// A constant FS_ASPECT_RATIO is defined. Given a box, those sides
// which can be split so that the ratio of the longest to shortest
// side does not exceed ASPECT_RATIO are identified. Among these
// sides, we select the one in which the points have the largest
// spread. We then split the points in a manner which most evenly
// distributes the points on either side of the splitting plane,
// subject to maintaining the bound on the ratio of long to short
// sides. To determine that the aspect ratio will be preserved,
// we determine the longest side (other than this side), and
// determine how narrowly we can cut this side, without causing the
// aspect ratio bound to be exceeded (small_piece).
//
// This procedure is more robust than either kd_split or midpt_split,
// but is more complicated as well. When point distribution is
// extremely skewed, this degenerates to midpt_split (actually
// 1/3 point split), and when the points are most evenly distributed,
// this degenerates to kd-split.
//----------------------------------------------------------------------
void fair_split(
ANNpointArray pa, // point array
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo) // num of points on low side (returned)
{
int d;
ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
cut_dim = 0;
for (d = 1; d < dim; d++) { // find length of longest box side
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (length > max_length) {
max_length = length;
cut_dim = d;
}
}
ANNcoord max_spread = 0; // find legal cut with max spread
cut_dim = 0;
for (d = 0; d < dim; d++) {
ANNcoord length = bnds.hi[d] - bnds.lo[d];
// is this side midpoint splitable
// without violating aspect ratio?
if (((double) max_length)*2.0/((double) length) <= FS_ASPECT_RATIO) {
// compute spread along this dim
ANNcoord spr = annSpread(pa, pidx, n, d);
if (spr > max_spread) { // best spread so far
max_spread = spr;
cut_dim = d; // this is dimension to cut
}
}
}
max_length = 0; // find longest side other than cut_dim
for (d = 0; d < dim; d++) {
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (d != cut_dim && length > max_length)
max_length = length;
}
// consider most extreme splits
ANNcoord small_piece = max_length / FS_ASPECT_RATIO;
ANNcoord lo_cut = bnds.lo[cut_dim] + small_piece;// lowest legal cut
ANNcoord hi_cut = bnds.hi[cut_dim] - small_piece;// highest legal cut
int br1, br2;
// is median below lo_cut ?
if (annSplitBalance(pa, pidx, n, cut_dim, lo_cut) >= 0) {
cut_val = lo_cut; // cut at lo_cut
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = br1;
}
// is median above hi_cut?
else if (annSplitBalance(pa, pidx, n, cut_dim, hi_cut) <= 0) {
cut_val = hi_cut; // cut at hi_cut
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = br2;
}
else { // median cut preserves asp ratio
n_lo = n/2; // split about median
annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
}
}
//----------------------------------------------------------------------
// sl_fair_split - sliding fair split splitting rule
//
// Sliding fair split is a splitting rule that combines the
// strengths of both fair split with sliding midpoint split.
// Fair split tends to produce balanced splits when the points
// are roughly uniformly distributed, but it can produce many
// trivial splits when points are highly clustered. Sliding
// midpoint never produces trivial splits, and shrinks boxes
// nicely if points are highly clustered, but it may produce
// rather unbalanced splits when points are unclustered but not
// quite uniform.
//
// Sliding fair split is based on the theory that there are two
// types of splits that are "good": balanced splits that produce
// fat boxes, and unbalanced splits provided the cell with fewer
// points is fat.
//
// This splitting rule operates by first computing the longest
// side of the current bounding box. Then it asks which sides
// could be split (at the midpoint) and still satisfy the aspect
// ratio bound with respect to this side. Among these, it selects
// the side with the largest spread (as fair split would). It
// then considers the most extreme cuts that would be allowed by
// the aspect ratio bound. This is done by dividing the longest
// side of the box by the aspect ratio bound. If the median cut
// lies between these extreme cuts, then we use the median cut.
// If not, then consider the extreme cut that is closer to the
// median. If all the points lie to one side of this cut, then
// we slide the cut until it hits the first point. This may
// violate the aspect ratio bound, but will never generate empty
// cells. However the sibling of every such skinny cell is fat,
// and hence packing arguments still apply.
//
//----------------------------------------------------------------------
void sl_fair_split(
ANNpointArray pa, // point array
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo) // num of points on low side (returned)
{
int d;
ANNcoord min, max; // min/max coordinates
int br1, br2; // split break points
ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
cut_dim = 0;
for (d = 1; d < dim; d++) { // find length of longest box side
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (length > max_length) {
max_length = length;
cut_dim = d;
}
}
ANNcoord max_spread = 0; // find legal cut with max spread
cut_dim = 0;
for (d = 0; d < dim; d++) {
ANNcoord length = bnds.hi[d] - bnds.lo[d];
// is this side midpoint splitable
// without violating aspect ratio?
if (((double) max_length)*2.0/((double) length) <= FS_ASPECT_RATIO) {
// compute spread along this dim
ANNcoord spr = annSpread(pa, pidx, n, d);
if (spr > max_spread) { // best spread so far
max_spread = spr;
cut_dim = d; // this is dimension to cut
}
}
}
max_length = 0; // find longest side other than cut_dim
for (d = 0; d < dim; d++) {
ANNcoord length = bnds.hi[d] - bnds.lo[d];
if (d != cut_dim && length > max_length)
max_length = length;
}
// consider most extreme splits
ANNcoord small_piece = max_length / FS_ASPECT_RATIO;
ANNcoord lo_cut = bnds.lo[cut_dim] + small_piece;// lowest legal cut
ANNcoord hi_cut = bnds.hi[cut_dim] - small_piece;// highest legal cut
// find min and max along cut_dim
annMinMax(pa, pidx, n, cut_dim, min, max);
// is median below lo_cut?
if (annSplitBalance(pa, pidx, n, cut_dim, lo_cut) >= 0) {
if (max > lo_cut) { // are any points above lo_cut?
cut_val = lo_cut; // cut at lo_cut
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = br1; // balance if there are ties
}
else { // all points below lo_cut
cut_val = max; // cut at max value
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = n-1;
}
}
// is median above hi_cut?
else if (annSplitBalance(pa, pidx, n, cut_dim, hi_cut) <= 0) {
if (min < hi_cut) { // are any points below hi_cut?
cut_val = hi_cut; // cut at hi_cut
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = br2; // balance if there are ties
}
else { // all points above hi_cut
cut_val = min; // cut at min value
annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
n_lo = 1;
}
}
else { // median cut is good enough
n_lo = n/2; // split about median
annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
}
}
//----------------------------------------------------------------------
// File: kd_split.h
// Programmer: Sunil Arya and David Mount
// Description: Methods for splitting 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
//----------------------------------------------------------------------
#ifndef ANN_KD_SPLIT_H
#define ANN_KD_SPLIT_H
#include "kd_tree.h" // kd-tree definitions
//----------------------------------------------------------------------
// External entry points
// These are all splitting procedures for kd-trees.
//----------------------------------------------------------------------
void kd_split( // standard (optimized) kd-splitter
ANNpointArray pa, // point array (unaltered)
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo); // num of points on low side (returned)
void midpt_split( // midpoint kd-splitter
ANNpointArray pa, // point array (unaltered)
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo); // num of points on low side (returned)
void sl_midpt_split( // sliding midpoint kd-splitter
ANNpointArray pa, // point array (unaltered)
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo); // num of points on low side (returned)
void fair_split( // fair-split kd-splitter
ANNpointArray pa, // point array (unaltered)
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo); // num of points on low side (returned)
void sl_fair_split( // sliding fair-split kd-splitter
ANNpointArray pa, // point array (unaltered)
ANNidxArray pidx, // point indices (permuted on return)
const ANNorthRect &bnds, // bounding rectangle for cell
int n, // number of points
int dim, // dimension of space
int &cut_dim, // cutting dimension (returned)
ANNcoord &cut_val, // cutting value (returned)
int &n_lo); // num of points on low side (returned)
#endif
//----------------------------------------------------------------------
// File: kd_tree.cpp
// Programmer: Sunil Arya and David Mount
// Description: Basic methods 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
// Revision 1.0 04/01/05
// Increased aspect ratio bound (ANN_AR_TOOBIG) from 100 to 1000.
// Fixed leaf counts to count trivial leaves.
// Added optional pa, pi arguments to Skeleton kd_tree constructor
// for use in load constructor.
// Added annClose() to eliminate KD_TRIVIAL memory leak.
//----------------------------------------------------------------------
#include "kd_tree.h" // kd-tree declarations
#include "kd_split.h" // kd-tree splitting rules
#include "kd_util.h" // kd-tree utilities
#include <ANN/ANNperf.h> // performance evaluation
//----------------------------------------------------------------------
// Global data
//
// For some splitting rules, especially with small bucket sizes,
// it is possible to generate a large number of empty leaf nodes.
// To save storage we allocate a single trivial leaf node which
// contains no points. For messy coding reasons it is convenient
// to have it reference a trivial point index.
//
// KD_TRIVIAL is allocated when the first kd-tree is created. It
// must *never* deallocated (since it may be shared by more than
// one tree).
//----------------------------------------------------------------------
static int IDX_TRIVIAL[] = {0}; // trivial point index
ANNkd_leaf *KD_TRIVIAL = NULL; // trivial leaf node
//----------------------------------------------------------------------
// Printing the kd-tree
// These routines print a kd-tree in reverse inorder (high then
// root then low). (This is so that if you look at the output
// from the right side it appear from left to right in standard
// inorder.) When outputting leaves we output only the point
// indices rather than the point coordinates. There is an option
// to print the point coordinates separately.
//
// The tree printing routine calls the printing routines on the
// individual nodes of the tree, passing in the level or depth
// in the tree. The level in the tree is used to print indentation
// for readability.
//----------------------------------------------------------------------
void ANNkd_split::print( // print splitting node
int level, // depth of node in tree
ostream &out) // output stream
{
child[ANN_HI]->print(level+1, out); // print high child
out << " ";
for (int i = 0; i < level; i++) // print indentation
out << "..";
out << "Split cd=" << cut_dim << " cv=" << cut_val;
out << " lbnd=" << cd_bnds[ANN_LO];
out << " hbnd=" << cd_bnds[ANN_HI];
out << "\n";
child[ANN_LO]->print(level+1, out); // print low child
}
void ANNkd_leaf::print( // print leaf node
int level, // depth of node in tree
ostream &out) // output stream
{
out << " ";
for (int i = 0; i < level; i++) // print indentation
out << "..";
if (this == KD_TRIVIAL) { // canonical trivial leaf node
out << "Leaf (trivial)\n";
}
else{
out << "Leaf n=" << n_pts << " <";
for (int j = 0; j < n_pts; j++) {
out << bkt[j];
if (j < n_pts-1) out << ",";
}
out << ">\n";
}
}
void ANNkd_tree::Print( // print entire tree
ANNbool with_pts, // print points as well?
ostream &out) // output stream
{
out << "ANN Version " << ANNversion << "\n";
if (with_pts) { // print point coordinates
out << " Points:\n";
for (int i = 0; i < n_pts; i++) {
out << "\t" << i << ": ";
annPrintPt(pts[i], dim, out);
out << "\n";
}
}
if (root == NULL) // empty tree?
out << " Null tree.\n";
else {
root->print(0, out); // invoke printing at root
}
}
//----------------------------------------------------------------------
// kd_tree statistics (for performance evaluation)
// This routine compute various statistics information for
// a kd-tree. It is used by the implementors for performance
// evaluation of the data structure.
//----------------------------------------------------------------------
#define MAX(a,b) ((a) > (b) ? (a) : (b))
void ANNkdStats::merge(const ANNkdStats &st) // merge stats from child
{
n_lf += st.n_lf; n_tl += st.n_tl;
n_spl += st.n_spl; n_shr += st.n_shr;
depth = MAX(depth, st.depth);
sum_ar += st.sum_ar;
}
//----------------------------------------------------------------------
// Update statistics for nodes
//----------------------------------------------------------------------
const double ANN_AR_TOOBIG = 1000; // too big an aspect ratio
void ANNkd_leaf::getStats( // get subtree statistics
int dim, // dimension of space
ANNkdStats &st, // stats (modified)
ANNorthRect &bnd_box) // bounding box
{
st.reset();
st.n_lf = 1; // count this leaf
if (this == KD_TRIVIAL) st.n_tl = 1; // count trivial leaf
double ar = annAspectRatio(dim, bnd_box); // aspect ratio of leaf
// incr sum (ignore outliers)
st.sum_ar += float(ar < ANN_AR_TOOBIG ? ar : ANN_AR_TOOBIG);
}
void ANNkd_split::getStats( // get subtree statistics
int dim, // dimension of space
ANNkdStats &st, // stats (modified)
ANNorthRect &bnd_box) // bounding box
{
ANNkdStats ch_stats; // stats for children
// get stats for low child
ANNcoord hv = bnd_box.hi[cut_dim]; // save box bounds
bnd_box.hi[cut_dim] = cut_val; // upper bound for low child
ch_stats.reset(); // reset
child[ANN_LO]->getStats(dim, ch_stats, bnd_box);
st.merge(ch_stats); // merge them
bnd_box.hi[cut_dim] = hv; // restore bound
// get stats for high child
ANNcoord lv = bnd_box.lo[cut_dim]; // save box bounds
bnd_box.lo[cut_dim] = cut_val; // lower bound for high child
ch_stats.reset(); // reset
child[ANN_HI]->getStats(dim, ch_stats, bnd_box);
st.merge(ch_stats); // merge them
bnd_box.lo[cut_dim] = lv; // restore bound
st.depth++; // increment depth
st.n_spl++; // increment number of splits
}
//----------------------------------------------------------------------
// getStats
// Collects a number of statistics related to kd_tree or
// bd_tree.
//----------------------------------------------------------------------
void ANNkd_tree::getStats( // get tree statistics
ANNkdStats &st) // stats (modified)
{
st.reset(dim, n_pts, bkt_size); // reset stats
// create bounding box
ANNorthRect bnd_box(dim, bnd_box_lo, bnd_box_hi);
if (root != NULL) { // if nonempty tree
root->getStats(dim, st, bnd_box); // get statistics
st.avg_ar = st.sum_ar / st.n_lf; // average leaf asp ratio
}
}
//----------------------------------------------------------------------
// kd_tree destructor
// The destructor just frees the various elements that were
// allocated in the construction process.
//----------------------------------------------------------------------
ANNkd_tree::~ANNkd_tree() // tree destructor
{
if (root != NULL) delete root;
if (pidx != NULL) delete [] pidx;
if (bnd_box_lo != NULL) annDeallocPt(bnd_box_lo);
if (bnd_box_hi != NULL) annDeallocPt(bnd_box_hi);
}
//----------------------------------------------------------------------
// This is called with all use of ANN is finished. It eliminates the
// minor memory leak caused by the allocation of KD_TRIVIAL.
//----------------------------------------------------------------------
void annClose() // close use of ANN
{
if (KD_TRIVIAL != NULL) {
delete KD_TRIVIAL;
KD_TRIVIAL = NULL;
}
}
//----------------------------------------------------------------------
// kd_tree constructors
// There is a skeleton kd-tree constructor which sets up a
// trivial empty tree. The last optional argument allows
// the routine to be passed a point index array which is
// assumed to be of the proper size (n). Otherwise, one is
// allocated and initialized to the identity. Warning: In
// either case the destructor will deallocate this array.
//
// As a kludge, we need to allocate KD_TRIVIAL if one has not
// already been allocated. (This is because I'm too dumb to
// figure out how to cause a pointer to be allocated at load
// time.)
//----------------------------------------------------------------------
void ANNkd_tree::SkeletonTree( // construct skeleton tree
int n, // number of points
int dd, // dimension
int bs, // bucket size
ANNpointArray pa, // point array
ANNidxArray pi) // point indices
{
dim = dd; // initialize basic elements
n_pts = n;
bkt_size = bs;
pts = pa; // initialize points array
root = NULL; // no associated tree yet
if (pi == NULL) { // point indices provided?
pidx = new ANNidx[n]; // no, allocate space for point indices
for (int i = 0; i < n; i++) {
pidx[i] = i; // initially identity
}
}
else {
pidx = pi; // yes, use them
}
bnd_box_lo = bnd_box_hi = NULL; // bounding box is nonexistent
if (KD_TRIVIAL == NULL) // no trivial leaf node yet?
KD_TRIVIAL = new ANNkd_leaf(0, IDX_TRIVIAL); // allocate it
}
ANNkd_tree::ANNkd_tree( // basic constructor
int n, // number of points
int dd, // dimension
int bs) // bucket size
{ SkeletonTree(n, dd, bs); } // construct skeleton tree
//----------------------------------------------------------------------
// rkd_tree - recursive procedure to build a kd-tree
//
// Builds a kd-tree for points in pa as indexed through the
// array pidx[0..n-1] (typically a subarray of the array used in
// the top-level call). This routine permutes the array pidx,
// but does not alter pa[].
//
// The construction is based on a standard algorithm for constructing
// the kd-tree (see Friedman, Bentley, and Finkel, ``An algorithm for
// finding best matches in logarithmic expected time,'' ACM Transactions
// on Mathematical Software, 3(3):209-226, 1977). The procedure
// operates by a simple divide-and-conquer strategy, which determines
// an appropriate orthogonal cutting plane (see below), and splits
// the points. When the number of points falls below the bucket size,
// we simply store the points in a leaf node's bucket.
//
// One of the arguments is a pointer to a splitting routine,
// whose prototype is:
//
// void split(
// ANNpointArray pa, // complete point array
// ANNidxArray pidx, // point array (permuted on return)
// ANNorthRect &bnds, // bounds of current cell
// int n, // number of points
// int dim, // dimension of space
// int &cut_dim, // cutting dimension
// ANNcoord &cut_val, // cutting value
// int &n_lo) // no. of points on low side of cut
//
// This procedure selects a cutting dimension and cutting value,
// partitions pa about these values, and returns the number of
// points on the low side of the cut.
//----------------------------------------------------------------------
ANNkd_ptr rkd_tree( // recursive construction of kd-tree
ANNpointArray pa, // point array
ANNidxArray pidx, // point indices to store in subtree
int n, // number of points
int dim, // dimension of space
int bsp, // bucket space
ANNorthRect &bnd_box, // bounding box for current node
ANNkd_splitter splitter) // splitting routine
{
if (n <= bsp) { // n small, make a leaf node
if (n == 0) // empty leaf node
return KD_TRIVIAL; // return (canonical) empty leaf
else // construct the node and return
return new ANNkd_leaf(n, pidx);
}
else { // n large, make a splitting node
int cd; // cutting dimension
ANNcoord cv; // cutting value
int n_lo; // number on low side of cut
ANNkd_node *lo, *hi; // low and high children
// invoke splitting procedure
(*splitter)(pa, pidx, bnd_box, n, dim, cd, cv, n_lo);
ANNcoord lv = bnd_box.lo[cd]; // save bounds for cutting dimension
ANNcoord hv = bnd_box.hi[cd];
bnd_box.hi[cd] = cv; // modify bounds for left subtree
lo = rkd_tree( // build left subtree
pa, pidx, n_lo, // ...from pidx[0..n_lo-1]
dim, bsp, bnd_box, splitter);
bnd_box.hi[cd] = hv; // restore bounds
bnd_box.lo[cd] = cv; // modify bounds for right subtree
hi = rkd_tree( // build right subtree
pa, pidx + n_lo, n-n_lo,// ...from pidx[n_lo..n-1]
dim, bsp, bnd_box, splitter);
bnd_box.lo[cd] = lv; // restore bounds
// create the splitting node
ANNkd_split *ptr = new ANNkd_split(cd, cv, lv, hv, lo, hi);
return ptr; // return pointer to this node
}
}
//----------------------------------------------------------------------
// kd-tree constructor
// This is the main constructor for kd-trees given a set of points.
// It first builds a skeleton tree, then computes the bounding box
// of the data points, and then invokes rkd_tree() to actually
// build the tree, passing it the appropriate splitting routine.
//----------------------------------------------------------------------
ANNkd_tree::ANNkd_tree( // construct from point array
ANNpointArray pa, // point array (with at least n pts)
int n, // number of points
int dd, // dimension
int bs, // bucket size
ANNsplitRule split) // splitting method
{
SkeletonTree(n, dd, bs); // set up the basic stuff
pts = pa; // where the points are
if (n == 0) return; // no points--no sweat
ANNorthRect bnd_box(dd); // bounding box for points
annEnclRect(pa, pidx, n, dd, bnd_box);// construct bounding rectangle
// copy to tree structure
bnd_box_lo = annCopyPt(dd, bnd_box.lo);
bnd_box_hi = annCopyPt(dd, bnd_box.hi);
switch (split) { // build by rule
case ANN_KD_STD: // standard kd-splitting rule
root = rkd_tree(pa, pidx, n, dd, bs, bnd_box, kd_split);
break;
case ANN_KD_MIDPT: // midpoint split
root = rkd_tree(pa, pidx, n, dd, bs, bnd_box, midpt_split);
break;
case ANN_KD_FAIR: // fair split
root = rkd_tree(pa, pidx, n, dd, bs, bnd_box, fair_split);
break;
case ANN_KD_SUGGEST: // best (in our opinion)
case ANN_KD_SL_MIDPT: // sliding midpoint split
root = rkd_tree(pa, pidx, n, dd, bs, bnd_box, sl_midpt_split);
break;
case ANN_KD_SL_FAIR: // sliding fair split
root = rkd_tree(pa, pidx, n, dd, bs, bnd_box, sl_fair_split);
break;
default:
annError("Illegal splitting method", ANNabort);
}
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This directory contains a modified version of GNU libmatheval.
Copyright (C) 1999, 2002, 2003 Aleksandar B. Samardzic
Copying and distribution of this file, with or without modification, are
permitted in any medium without royalty provided the copyright notice
and this notice are preserved.
WHAT IS IT?
GNU libmatheval is a library which contains several procedures that make
it possible to create an in-memory tree from the string representation
of a mathematical function over single or multiple variables. This tree
can be used later to evaluate a function for specified variable values,
to create a corresponding tree for the function derivative over a
specified variable, or to write a textual tree representation to a
specified string.
BUGS
Please report bugs and eventually send patches to
<bug-libmatheval@gnu.org> mailing list.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
typedef union {
Node *node;
Record *record;
} YYSTYPE;
#define NUMBER 257
#define VARIABLE 258
#define FUNCTION 259
#define NEG 260
extern YYSTYPE melval;
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment