kd_split.cpp

c++实现的KNN库:建立高维度的K-d tree,实现K邻域搜索

//----------------------------------------------------------------------
// 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;
}

//----------------------------------------------------------------------