text: //----------------------------------------------------------------------
  // 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);
  	}
  }
  