diff --git a/sample/testkd2.cpp b/sample/testkd2.cpp
index 0226ba0..2b79408 100644
--- a/sample/testkd2.cpp
+++ b/sample/testkd2.cpp
@@ -5,6 +5,7 @@
 #include <fstream>
 #define __KD_TREE_NUMNODES
 #include "mykdtree.hpp"
+#include "kdtree_splitters.hpp"
 
 #define NTRY 100 
 #define ND 3
@@ -12,7 +13,7 @@
 using namespace std;
 using namespace CosmoTool;
 
-typedef KDTree<ND,char> MyTree;
+typedef KDTree<ND,char,ComputePrecision,KD_homogeneous_cell_splitter<ND, char> > MyTree;
 typedef KDCell<ND,char> MyCell;
 
 MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells)
diff --git a/src/dinterpolate.tcc b/src/dinterpolate.tcc
index ecd04ae..58e7dfb 100644
--- a/src/dinterpolate.tcc
+++ b/src/dinterpolate.tcc
@@ -217,7 +217,7 @@ namespace CosmoTool {
        kdc[i] = c[i];
 
     // It may happen that we are unlucky and have to iterate to farther
-    // neighbors. It should happen, especially on the boundaries.
+    // neighbors. It is bound to happen, especially on the boundaries.
     do
       {
 	uint32_t i;
diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp
index d661376..c2087b5 100644
--- a/src/mykdtree.hpp
+++ b/src/mykdtree.hpp
@@ -59,6 +59,7 @@ namespace CosmoTool {
     uint32_t currentRank;
     uint32_t numCells;
   };
+  
 
   template<int N, typename ValType, typename CType = ComputePrecision>
   class RecursionMultipleInfo
@@ -76,7 +77,13 @@ namespace CosmoTool {
     }
   };
 
-  template<int N, typename ValType, typename CType = ComputePrecision>
+  template<int N, typename ValType, typename CType = ComputePrecision> 
+  struct KD_default_cell_splitter
+  {
+    void operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound);
+  };
+
+  template<int N, typename ValType, typename CType = ComputePrecision, typename CellSplitter = KD_default_cell_splitter<N,ValType,CType> >
   class KDTree
   {
   public:
@@ -84,6 +91,8 @@ namespace CosmoTool {
     typedef typename KDDef<N>::KDCoordinates coords;
     typedef KDCell<N,ValType,CType> Cell;
     typedef KDTreeNode<N,ValType,CType> Node;
+    
+    CellSplitter splitter;
 
     KDTree(Cell *cells, uint32_t Ncells);
     ~KDTree();
diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc
index 9899d60..6f5eebf 100644
--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@@ -23,13 +23,13 @@ namespace CosmoTool {
     int rank;
   };
 
-  template<int N, typename ValType, typename CType>
-  KDTree<N,ValType,CType>::~KDTree()
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::~KDTree()
   {
   }
 
-  template<int N, typename ValType, typename CType>
-  KDTree<N,ValType,CType>::KDTree(Cell *cells, uint32_t Ncells)
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::KDTree(Cell *cells, uint32_t Ncells)
   {
     
     base_cell = cells;
@@ -43,8 +43,8 @@ namespace CosmoTool {
     optimize();
   }
 
-  template<int N, typename ValType, typename CType>
-  void KDTree<N,ValType,CType>::optimize()
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void KDTree<N,ValType,CType,CellSplitter>::optimize()
   {
     coords absoluteMin, absoluteMax;
 
@@ -76,9 +76,9 @@ namespace CosmoTool {
     std::cout << "   done." << std::endl;
   }
 
-  template<int N, typename ValType, typename CType>
-  uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r, 
-						    KDTree<N,ValType,CType>::Cell **cells,
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
+						    KDTree<N,ValType,CType,CellSplitter>::Cell **cells,
 						    uint32_t numCells)
     throw (NotEnoughCells)
   {
@@ -96,8 +96,8 @@ namespace CosmoTool {
     return info.currentRank;
   }
 
-  template<int N, typename ValType, typename CType>
-  uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r, 
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::getIntersection(const coords& x, CoordType r, 
 						    Cell **cells,
 						    CoordType *distances,
 						    uint32_t numCells)
@@ -117,8 +117,8 @@ namespace CosmoTool {
     return info.currentRank;
   }
 
-  template<int N, typename ValType, typename CType>
-  uint32_t KDTree<N,ValType,CType>::countCells(const coords& x, CoordType r)
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::countCells(const coords& x, CoordType r)
   {
     RecursionInfoCells<N,ValType> info;
 
@@ -134,9 +134,9 @@ namespace CosmoTool {
     return info.currentRank;
   }
 
-  template<int N, typename ValType, typename CType>
+  template<int N, typename ValType, typename CType, typename CellSplitter>
   template<bool justCount>
-  void KDTree<N,ValType,CType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
+  void KDTree<N,ValType,CType,CellSplitter>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
 							   Node *node,
 							   int level)
     throw (NotEnoughCells)
@@ -208,25 +208,31 @@ namespace CosmoTool {
   }
 
   template<int N, typename ValType, typename CType>
-  KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType>::buildTree(Cell **cell0,
-								    uint32_t Ncells,
-								    uint32_t depth,
-								    coords minBound,
-								    coords maxBound)
+  void  KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound)
+  {
+    CellCompare<N,ValType,CType> compare(axis);
+    std::sort(cells, cells+Ncells, compare);
+    split_index = Ncells/2;
+  }
+
+
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType,CellSplitter>::buildTree(Cell **cell0,
+								  uint32_t Ncells,
+								  uint32_t depth,
+								  coords minBound,
+								  coords maxBound)
   {
     if (Ncells == 0)
       return 0;
 
     int axis = depth % N;
     Node *node = &nodes[lastNode++];
-    uint32_t mid = Ncells/2;
+    uint32_t mid;
     coords tmpBound;
     
     // Isolate the environment
-    {
-      CellCompare<N,ValType,CType> compare(axis);
-      std::sort(cell0, cell0+Ncells, compare);
-    }
+    splitter(cell0, Ncells, mid, axis, minBound, maxBound);
 
     node->value = *(cell0+mid);
     memcpy(&node->minBound[0], &minBound[0], sizeof(coords));
@@ -252,8 +258,8 @@ namespace CosmoTool {
     return node;
   }
 
-  template<int N, typename ValType, typename CType>
-  uint32_t KDTree<N,ValType,CType>::countActives() const
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  uint32_t KDTree<N,ValType,CType,CellSplitter>::countActives() const
   {
     uint32_t numActive = 0;
     for (uint32_t i = 0; i < lastNode; i++)
@@ -264,9 +270,9 @@ namespace CosmoTool {
     return numActive;
   }
 
-  template<int N, typename ValType, typename CType>
+  template<int N, typename ValType, typename CType, typename CellSplitter>
   typename KDDef<N,CType>::CoordType
-  KDTree<N,ValType,CType>::computeDistance(const Cell *cell, const coords& x) const
+  KDTree<N,ValType,CType,CellSplitter>::computeDistance(const Cell *cell, const coords& x) const
   {
     CoordType d2 = 0;
 
@@ -278,9 +284,9 @@ namespace CosmoTool {
     return d2;
   }
 
-  template<int N, typename ValType, typename CType>
+  template<int N, typename ValType, typename CType, typename CellSplitter>
   void
-  KDTree<N,ValType,CType>::recursiveNearest(
+  KDTree<N,ValType,CType,CellSplitter>::recursiveNearest(
 				       Node *node,
 				       int level,
 				       const coords& x,
@@ -352,9 +358,9 @@ namespace CosmoTool {
       }
   }
 
- template<int N, typename ValType, typename CType>
+ template<int N, typename ValType, typename CType, typename CellSplitter>
  KDCell<N,ValType,CType> *
- KDTree<N,ValType,CType>::getNearestNeighbour(const coords& x)
+ KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbour(const coords& x)
  {
    CoordType R2 = INFINITY;
    Cell *best = 0;
@@ -364,9 +370,9 @@ namespace CosmoTool {
    return best;
  }
 
-  template<int N, typename ValType, typename CType>
+  template<int N, typename ValType, typename CType, typename CellSplitter>
   void
-  KDTree<N,ValType,CType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
+  KDTree<N,ValType,CType,CellSplitter>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
 					       int level)
   {
     CoordType d2 = 0;
@@ -421,8 +427,8 @@ namespace CosmoTool {
       }    
   }
 
- template<int N, typename ValType, typename CType>
- void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
+ template<int N, typename ValType, typename CType, typename CellSplitter>
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
 						    Cell **cells)
  {
    RecursionMultipleInfo<N,ValType> info(x, cells, N2);
@@ -435,8 +441,8 @@ namespace CosmoTool {
    //   std::cout << "Traversed = " << info.traversed << std::endl;
  }
 
- template<int N, typename ValType, typename CType>
- void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
+ template<int N, typename ValType, typename CType, typename CellSplitter>
+ void KDTree<N,ValType,CType,CellSplitter>::getNearestNeighbours(const coords& x, uint32_t N2,
 						    Cell **cells,
 						    CoordType *distances)
  {
@@ -470,8 +476,8 @@ namespace CosmoTool {
     long rootId;
   };
 
-  template<int N, typename ValType, typename CType>
-  void KDTree<N,ValType,CType>::saveTree(std::ostream& o) const
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  void KDTree<N,ValType,CType,CellSplitter>::saveTree(std::ostream& o) const
   {
     KDTreeHeader h;
 
@@ -505,8 +511,8 @@ namespace CosmoTool {
       }
   }
 
-  template<int N, typename ValType, typename CType>
-  KDTree<N,ValType,CType>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells)
+  template<int N, typename ValType, typename CType, typename CellSplitter>
+  KDTree<N,ValType,CType,CellSplitter>::KDTree(std::istream& in, Cell *cells, uint32_t Ncells)
     throw (InvalidOnDiskKDTree)
   {
     KDTreeHeader h;