// "buckytree" (c)2007 Bryan McNett. All Rights Reserved.
// bmcnett at gmail dot com

// sort triangles in O(NlogN) time, 
// then search them in O(logN)? time,
// without requiring a spatial index

// the database contains "minimum" values for each convex hull plane,
// and the search query contains "maximum" values.

#define DIMENSIONS 3
#define AABB

#ifdef AABB
#define DIRECTIONS (DIMENSIONS*2) // AABB is easier to understand
#else
#define DIRECTIONS (DIMENSIONS+1) // simplex ( triangle, tetrahedron ) is more efficient
#endif

struct vector { float v[DIMENSIONS]; };
struct triangle { vector v[3]; };
struct hull { float plane[DIRECTIONS]; }; // for 3D AABB, this is [ max(x), max(-x), max(y), max(-y), max(z), max(-z) ]

vector basis[DIRECTIONS] 
#if DIMENSIONS == 3
= { 
#if defined(AABB)
    {{  1,  0,  0 }}, {{ -1,  0,  0 }}, 
    {{  0,  1,  0 }}, {{  0, -1,  0 }}, 
    {{  0,  0,  1 }}, {{  0,  0, -1 }} 
#else
    {{  0.000,  0.000,  1.000 }},
    {{  0.943,  0.000, -0.333 }},
    {{ -0.471,  0.816, -0.333 }},
    {{ -0.471, -0.816, -0.333 }}
#endif // defined(AABB)
  }
#endif // DIMENSIONS == 3
;

float dot( const vector& a, const vector& b ) 
{
  float result = a.v[0] * b.v[0];
  for( int i=1; i<DIMENSIONS; ++i )
    result += a.v[i] * b.v[i];
  return result;
}

struct comparer
{
  int direction;
  comparer( int d ) : direction(d) {}
  float minimum( const triangle& a ) const
  {
    vector b = basis[direction];
    return min( dot( a.v[0], b ), dot( a.v[1], b ), dot( a.v[2], b ) );
  }
  bool operator<( const triangle& a, const triangle& b ) const
  {
    return minimum(a) < minimum(b);
  }
};

// sort an array of triangles in O(NlogN) time
//
void buckysort( triangle* begin, triangle* end, int direction )
{
  int triangles = end - begin;
  if( triangles < 2 ) 
    return;
  triangle* median = begin + triangles / 2;

  comparer predicate( direction );
  std::nth_element( begin, median, end, predicate );

  // rotate vertices of triangle until "smallest" vector is in index 0

  float minimum = predicate.minimum( *median );
  while( dot( median->v[0], basis[direction] ) > minimum )
  {
    swap( median->v[0], median->v[1] );
    swap( median->v[1], median->v[2] );
  }

  buckysort( begin, median, (direction+1)%DIRECTIONS );
  buckysort( median+1, end, (direction+1)%DIRECTIONS );  
}

// search a sorted array of triangles in O(logN)? time
//
void buckysearch( triangle** results, triangle* begin, triangle* end, hull* query, int direction )
{
  int triangles = end - begin;
  if( triangles < 1 )
    return;
  triangle* median = begin + triangles / 2;

  buckysearch( results, begin, median, query, (direction+1)%DIRECTIONS ); // can't reject near side of plane
  if( dot( median->v[0], basis[direction] ) > query->plane[direction] )
    return; // reject median and far side of plane
  *results++ = median; // we should probably do a real tri-ray check here
  buckysearch( results, median+1, end, query, (direction+1)%DIRECTIONS );
}