#ifndef __KDTREE_SPLITTERS_HPP
#define __KDTREE_SPLITTERS_HPP

#include <algorithm>

namespace CosmoTool
{

  template<int N, typename ValType, typename CType = ComputePrecision> 
  struct KD_homogeneous_cell_splitter
  {
    typedef typename KDDef<N,CType>::KDCoordinates coords; 
    typedef typename KDDef<N,CType>::CoordType ctype; 

    void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound)
    {
      ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]);
      uint32_t below = 0, above = Ncells-1;
      ctype delta_max = std::abs(cells[0]->coord[axis]-midCoord);
      uint32_t idx_max = 0;

      while (below < above)
	{
          ctype delta = cells[below]->coord[axis]-midCoord;
	  if (delta > 0)
            {
              if (delta < delta_max)
                {
                  delta_max = delta;
                  idx_max = above;
                }
	      std::swap(cells[below], cells[above--]);
            }
	  else
            {
              if (-delta < delta_max)
                {
                  delta_max = -delta;
                  idx_max = below;
                }
	      below++;
            }
	}
      if (idx_max != above)
        {
          bool cond1 = cells[idx_max]->coord[axis] > midCoord;
          bool cond2 = cells[above]->coord[axis] > midCoord;
          if ((cond1 && cond2) || (!cond1 && !cond2))
            {
              split_index = above;
              std::swap(cells[above], cells[idx_max]);
            }
          else if (cond2)
            {
              split_index = above-1;
              std::swap(cells[above-1], cells[idx_max]);
            }
          else
            {
              split_index = above+1;
              std::swap(cells[above+1], cells[idx_max]);
            }
         }
       else split_index = above;
    }
  };


};

#endif