diff --git a/Spatial_searching/benchmark/Spatial_searching/CMakeLists.txt b/Spatial_searching/benchmark/Spatial_searching/CMakeLists.txt
index 39144c0325a..79f9dafedc3 100644
--- a/Spatial_searching/benchmark/Spatial_searching/CMakeLists.txt
+++ b/Spatial_searching/benchmark/Spatial_searching/CMakeLists.txt
@@ -33,6 +33,7 @@ if ( CGAL_FOUND )
   create_single_source_cgal_program( "Nearest_neighbor_searching_2D.cpp" )
   create_single_source_cgal_program( "Nearest_neighbor_searching_2D_user_defined.cpp" )
   create_single_source_cgal_program( "Split_data.cpp" )
+create_single_source_cgal_program( "nn3cgal.cpp" )
 
 else()
   
diff --git a/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp b/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp
new file mode 100644
index 00000000000..cc924681225
--- /dev/null
+++ b/Spatial_searching/benchmark/Spatial_searching/nn3cgal.cpp
@@ -0,0 +1,109 @@
+
+#include <CGAL/Simple_cartesian.h>
+#include <CGAL/K_neighbor_search.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Search_traits_3.h>
+#include <CGAL/Search_traits_2.h>
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Timer.h>
+#include <CGAL/Memory_sizer.h>
+
+#include <iostream>
+#include <fstream>
+
+#include <boost/lexical_cast.hpp>
+
+typedef CGAL::Simple_cartesian<double> K;
+typedef CGAL::Search_traits_3<K> TreeTraits_3;
+//typedef CGAL::Search_traits_2<K> TreeTraits_3;
+typedef K::Point_3 Point_3;
+typedef CGAL::Euclidean_distance<TreeTraits_3> Distance;
+typedef CGAL::Sliding_midpoint<TreeTraits_3> Splitter;
+typedef CGAL::Orthogonal_k_neighbor_search<TreeTraits_3,Distance,Splitter> Neighbor_search_3;
+typedef Neighbor_search_3::Tree Tree_3;
+typedef CGAL::Timer Timer;
+
+
+typedef std::vector<Point_3> Points;
+		 
+
+void read(Points& points, char* argv)
+{
+  int d, n;
+  Point_3 p;
+#if 0
+  std::ifstream data(argv);
+  data >> d >> n;
+  assert(d == 3);
+  points.reserve(n);
+  while(data >> p){
+    points.push_back(p);
+  }
+  data.close();
+    
+#else
+ std::ifstream data(argv, std::ios::in | std::ios::binary);
+  CGAL::set_binary_mode(data);
+  CGAL::read(data,d);
+  CGAL::read(data,n);
+
+  points.reserve(n);
+  for(int i=0; i < n; i++){
+    data >> p;
+    points.push_back(p);
+  }
+#endif
+  data.close();
+}
+
+int main(int argc,char *argv[]) 
+{
+  int NN_number = 10;
+  Points query_points_3, data_points_3;
+
+  Point_3 p;
+  Timer t;
+
+  read(data_points_3, argv[1]);
+  read(query_points_3, argv[2]);
+ 
+  int bucketsize = (argc>3) ? boost::lexical_cast<int>(argv[3]) : 10;
+  std::cerr << "bucketsize = "  << bucketsize <<std::endl;
+
+  t.start();
+  // Insert data points in the tree
+  Splitter splitter(bucketsize);
+  Tree_3 tree(data_points_3.begin(), data_points_3.end(),splitter);
+  tree.build();
+  t.stop();
+  std::cerr << "build " << t.time() << " sec\n";
+  int items=0,leafs=0,internals=0;
+  Points result(NN_number);
+
+  t.reset();t.start();
+  bool dump = true;
+  for(Points::iterator it = query_points_3.begin(); it != query_points_3.end(); ++it) {
+    // Initialize the search structure, and search all NN_number neighbors
+    Neighbor_search_3 search(tree, *it, NN_number);
+    int i=0;
+    for (Neighbor_search_3::iterator it = search.begin(); it != search.end();it++, i++) {
+      result[i] = it->first;
+      if(dump){
+        std::cerr << result[i] << std::endl;
+      }
+	}
+	if(dump){
+		tree.statistics(std::cerr);
+		//tree.print();
+	}
+    dump = false;
+//	items += search.items_visited();
+//	leafs += search.leafs_visited();
+//	internals += search.internals_visited();
+  }
+  t.stop();
+  std::cerr << "queries " << t.time() << " sec\n";
+  std::cerr << items <<" items "<<internals<<" internals "<<leafs<< " leafs visited\n";
+  std::cerr << "done\n";
+  return 0;
+}
diff --git a/Spatial_searching/include/CGAL/Euclidean_distance.h b/Spatial_searching/include/CGAL/Euclidean_distance.h
index 4d1031431af..4cbbbd2fdc7 100644
--- a/Spatial_searching/include/CGAL/Euclidean_distance.h
+++ b/Spatial_searching/include/CGAL/Euclidean_distance.h
@@ -46,6 +46,7 @@ namespace CGAL {
     typedef typename SearchTraits::FT    FT;
     typedef typename SearchTraits::Point_d Point_d;
     typedef Point_d Query_item;
+	typedef typename SearchTraits::Dimension D;
 
     	// default constructor
     	Euclidean_distance(const SearchTraits& traits_=SearchTraits()):traits(traits_) {}
@@ -64,7 +65,7 @@ namespace CGAL {
 
 
 	inline FT min_distance_to_rectangle(const Query_item& q,
-					    const Kd_tree_rectangle<FT>& r) const {
+					    const Kd_tree_rectangle<FT,D>& r) const {
 		FT distance = FT(0);
 		typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
                 typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
@@ -82,7 +83,7 @@ namespace CGAL {
 	}
 
 	inline FT max_distance_to_rectangle(const Query_item& q,
-					     const Kd_tree_rectangle<FT>& r) const {
+					     const Kd_tree_rectangle<FT,D>& r) const {
 		FT distance=FT(0);
 		typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
                 typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
diff --git a/Spatial_searching/include/CGAL/K_neighbor_search.h b/Spatial_searching/include/CGAL/K_neighbor_search.h
index 434650a1206..4ec5474a0ad 100644
--- a/Spatial_searching/include/CGAL/K_neighbor_search.h
+++ b/Spatial_searching/include/CGAL/K_neighbor_search.h
@@ -34,6 +34,7 @@ class K_neighbor_search: public internal::K_neighbor_search<SearchTraits,Distanc
   
 public:
   typedef typename Base::FT FT;  
+  typedef typename SearchTraits::DIMENSION D;
 
   K_neighbor_search(const Tree& tree, const typename Base::Query_item& q,  
     unsigned int k=1, FT Eps=FT(0.0), bool Search_nearest=true, const Distance& d=Distance(),bool sorted=true)
@@ -49,7 +50,7 @@ private:
   using Base::branch;
 
   void 
-  compute_neighbors_general(typename Base::Node_const_handle N, const Kd_tree_rectangle<FT>& r) 
+  compute_neighbors_general(typename Base::Node_const_handle N, const Kd_tree_rectangle<FT,D>& r) 
   {
     if (!(N->is_leaf())) {
       this->number_of_internal_nodes_visited++;
diff --git a/Spatial_searching/include/CGAL/Kd_tree.h b/Spatial_searching/include/CGAL/Kd_tree.h
index 2d541378128..9e5df76f019 100644
--- a/Spatial_searching/include/CGAL/Kd_tree.h
+++ b/Spatial_searching/include/CGAL/Kd_tree.h
@@ -60,6 +60,7 @@ public:
   typedef typename std::vector<Point_d>::const_iterator const_iterator;
 
   typedef typename std::vector<Point_d>::size_type size_type;
+  typedef typename SearchTraits::Dimension D;
 
 private:
   SearchTraits traits_;
@@ -68,7 +69,7 @@ private:
 
   Node_handle tree_root;
 
-  Kd_tree_rectangle<FT>* bbox;
+  Kd_tree_rectangle<FT,D>* bbox;
   std::vector<Point_d> pts;
 
   // Instead of storing the points in arrays in the Kd_tree_node
@@ -210,7 +211,7 @@ public:
       data.push_back(&pts[i]);
     }
     Point_container c(dim, data.begin(), data.end(),traits_);
-    bbox = new Kd_tree_rectangle<FT>(c.bounding_box());
+    bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box());
     if (c.size() <= split.bucket_size()){
       tree_root = create_leaf_node(c);
     }else {
@@ -335,7 +336,7 @@ public:
     root()->print();
   }
 
-  const Kd_tree_rectangle<FT>&
+  const Kd_tree_rectangle<FT,D>&
   bounding_box() const
   {
     if(! is_built()){
diff --git a/Spatial_searching/include/CGAL/Kd_tree_node.h b/Spatial_searching/include/CGAL/Kd_tree_node.h
index 9d082536a87..43b6b87884c 100644
--- a/Spatial_searching/include/CGAL/Kd_tree_node.h
+++ b/Spatial_searching/include/CGAL/Kd_tree_node.h
@@ -42,6 +42,7 @@ namespace CGAL {
     typedef typename TreeTraits::FT FT;
     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
     typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator;
+	typedef typename TreeTraits::Dimension D;
 
   private:
 
@@ -264,7 +265,7 @@ namespace CGAL {
     template <class OutputIterator, class FuzzyQueryItem>
     OutputIterator 
     search(OutputIterator it, const FuzzyQueryItem& q,
-	   Kd_tree_rectangle<FT>& b) const
+	   Kd_tree_rectangle<FT,D>& b) const
     {
       if (is_leaf()) { 
 	if (n>0) 
diff --git a/Spatial_searching/include/CGAL/Kd_tree_rectangle.h b/Spatial_searching/include/CGAL/Kd_tree_rectangle.h
index a10ac002749..639647a9d8f 100644
--- a/Spatial_searching/include/CGAL/Kd_tree_rectangle.h
+++ b/Spatial_searching/include/CGAL/Kd_tree_rectangle.h
@@ -25,6 +25,8 @@
 #include <algorithm>
 #include <new>
 #include <CGAL/assertions.h>
+#include <CGAL/array.h>
+#include <CGAL/Dimension.h>
 
 namespace CGAL {
 
@@ -53,8 +55,189 @@ namespace CGAL {
   };
 
 
-  template <class FT_> 
+  template <class FT_, typename D> 
   class Kd_tree_rectangle {
+  public:
+    typedef FT_ FT;
+    typedef FT T;
+    
+  private:
+    
+    //int dim;
+    CGAL::cpp11::array<T,D::value> lower_;
+    CGAL::cpp11::array<T,D::value> upper_;
+    int max_span_coord_;
+    
+  public:
+    
+    inline void 
+    set_upper_bound(int i, const FT& x) 
+    {
+      CGAL_assertion(i >= 0 && i < D);
+      CGAL_assertion(x >= lower_[i]);
+      upper_[i] = x;
+      set_max_span();
+    }
+
+    inline void 
+    set_lower_bound(int i, const FT& x) 
+    {
+      CGAL_assertion(i >= 0 && i < D);
+      CGAL_assertion(x <= upper_[i]);
+      lower_[i] = x;
+      set_max_span();
+    }
+
+    inline void 
+    set_max_span() 
+    {
+      FT span = upper_[0]-lower_[0];
+      max_span_coord_ = 0;
+      for (int i = 1; i < D::value; ++i) {
+	FT tmp = upper_[i] - lower_[i];
+	if (span < tmp) {
+	  span = tmp;
+	  max_span_coord_ = i;
+	}
+      }
+    }
+    
+    Kd_tree_rectangle(int d) 
+      : max_span_coord_(0)
+    {
+      lower_.assign(FT(0));
+      upper_.assign(FT(0));
+    }
+
+    Kd_tree_rectangle() 
+    {}
+    
+    
+    explicit 
+    Kd_tree_rectangle(const Kd_tree_rectangle<FT,D>& r) 
+    : max_span_coord_(r.max_span_coord_) 
+    {
+      lower_ = r.lower_;
+      upper_ = r.upper_;
+    }
+
+    template <class Construct_cartesian_const_iterator_d,class PointPointerIter>
+    void update_from_point_pointers(PointPointerIter begin, 
+                                    PointPointerIter end,
+                                    const Construct_cartesian_const_iterator_d& construct_it
+    ) 
+    {
+      if (begin ==end)
+        return;
+      // initialize with values of first point
+      typename Construct_cartesian_const_iterator_d::result_type bit = construct_it(**begin);
+      
+      for (int i=0; i < D::value; ++i, ++bit) {
+        lower_[i]= *bit; upper_[i]=lower_[i];
+      }
+      begin++;
+      typedef typename std::iterator_traits<PointPointerIter>::value_type P;
+      std::for_each(begin, end,set_bounds_from_pointer<Construct_cartesian_const_iterator_d,P,T>(D::value, &(lower_[0]), &(upper_[0]), construct_it));
+      set_max_span();
+    }
+    
+    template <class Construct_cartesian_const_iterator_d,class PointPointerIter> // was PointIter
+    Kd_tree_rectangle(int d,  PointPointerIter begin,  PointPointerIter end,const Construct_cartesian_const_iterator_d& construct_it)
+    {
+      update_from_point_pointers<Construct_cartesian_const_iterator_d>(begin,end,construct_it);
+    }
+
+    inline int 
+    max_span_coord() const 
+    { 
+      return max_span_coord_; 
+    }
+
+    inline FT 
+    max_span() const 
+    {
+      return upper_[max_span_coord_] - lower_[max_span_coord_];
+    }
+
+    inline FT 
+    min_coord(int i) const 
+    {
+      CGAL_assertion(lower_ != NULL);
+      return lower_[i];
+    }
+    
+    inline FT 
+    max_coord(int i) const 
+    {
+      return upper_[i];
+    }
+    
+    std::ostream& 
+    print(std::ostream& s) const 
+    {
+      s << "Rectangle dimension = " << D;
+      s << "\n lower: ";
+      for (int i=0; i < D; ++i)
+	s << lower_[i] << " ";
+      // std::copy(lower_, lower_ + D,
+      // 	      std::ostream_iterator<FT>(s," "));
+      s << "\n upper: ";
+      for (int j=0; j < D; ++j)
+	s << upper_[j] << " ";
+      // std::copy(upper_, upper_ + D,
+      //	      std::ostream_iterator<FT>(s," "));
+      s << "\n maximum span " << max_span() <<
+	" at coordinate " << max_span_coord() << std::endl;
+      return s;
+    }
+    
+    // Splits rectangle by modifying itself to lower half 
+    // and returns upper half
+    //    Kd_tree_rectangle* 
+    void
+    split(Kd_tree_rectangle& r, int d, FT value)
+    {
+      CGAL_assertion(d >= 0 && d < D);
+      CGAL_assertion(lower_[d] <= value && value <= upper_[d]);
+      
+      //Kd_tree_rectangle* r = new Kd_tree_rectangle(*this);
+      upper_[d]=value;
+      r.lower_[d]=value;
+      //return r;
+    }
+    
+
+    int 
+    dimension() const 
+    {
+      return D::value;
+    }
+
+    const T* lower() const {return lower_;}
+    const T* upper() const {return upper_;}    
+ 
+    Kd_tree_rectangle<FT,D>& 
+    operator=(const Kd_tree_rectangle<FT,D>& r) 
+    {
+      CGAL_assertion(dimension() == r.dimension());
+      if (this != &r) {
+        lower_ = r.lower_;
+        upper_ = r.upper_;
+	set_max_span();
+      }
+      return *this;
+    }
+
+  
+
+  }; // of class Kd_tree_rectangle
+
+
+  // Partial specialization for dimension == 0, which means dimension at runtime
+
+  template <class FT_> 
+  class Kd_tree_rectangle<FT_, Dynamic_dimension_tag> {
+
   public:
     typedef FT_ FT;
     typedef FT T;
@@ -62,8 +245,8 @@ namespace CGAL {
   private:
     
     int dim;
-    T *lower_;
-    T *upper_;
+    T* lower_;
+    T* upper_;
     int max_span_coord_;
     
   public:
@@ -109,11 +292,12 @@ namespace CGAL {
 
     Kd_tree_rectangle() 
       : dim(0), lower_(0), upper_(0) 
-    {}
+    {
+}
     
     
     explicit 
-    Kd_tree_rectangle(const Kd_tree_rectangle<FT>& r) 
+    Kd_tree_rectangle(const Kd_tree_rectangle<FT,Dynamic_dimension_tag>& r) 
       : dim(r.dim), lower_(new FT[dim]), upper_(new FT[dim]), 
 	max_span_coord_(r.max_span_coord_) 
     {
@@ -143,7 +327,7 @@ namespace CGAL {
     
     template <class Construct_cartesian_const_iterator_d,class PointPointerIter> // was PointIter
     Kd_tree_rectangle(int d,  PointPointerIter begin,  PointPointerIter end,const Construct_cartesian_const_iterator_d& construct_it)
-      : dim(d), lower_(new FT[d]), upper_(new FT[d]) 
+      : dim(d), lower_(new FT[d]), upper_(new FT[d])
     {
       update_from_point_pointers<Construct_cartesian_const_iterator_d>(begin,end,construct_it);
     }
@@ -215,7 +399,7 @@ namespace CGAL {
 	if (upper_) delete [] upper_;
       }
     }
-
+    
     int 
     dimension() const 
     {
@@ -225,8 +409,8 @@ namespace CGAL {
     const T* lower() const {return lower_;}
     const T* upper() const {return upper_;}    
  
-    Kd_tree_rectangle<FT>& 
-    operator=(const Kd_tree_rectangle<FT>& r) 
+    Kd_tree_rectangle<FT,Dynamic_dimension_tag>& 
+    operator=(const Kd_tree_rectangle<FT,Dynamic_dimension_tag>& r) 
     {
       CGAL_assertion(dimension() == r.dimension());
       if (this != &r) {
@@ -239,11 +423,11 @@ namespace CGAL {
 
   
 
-  }; // of class Kd_tree_rectangle
+  }; // of partial specialization of class Kd_tree_rectangle<FT,0>
 
-  template <class FT>
+  template <class FT, typename D>
   std::ostream& 
-  operator<<(std::ostream& s, const Kd_tree_rectangle<FT>& r) 
+  operator<<(std::ostream& s, const Kd_tree_rectangle<FT,D>& r) 
   {
     return r.print(s);
   }
diff --git a/Spatial_searching/include/CGAL/Point_container.h b/Spatial_searching/include/CGAL/Point_container.h
index b3f44264756..6ca8100ef6a 100644
--- a/Spatial_searching/include/CGAL/Point_container.h
+++ b/Spatial_searching/include/CGAL/Point_container.h
@@ -45,7 +45,7 @@ public:
   
   typedef typename Point_vector::iterator iterator;
   typedef typename Point_vector::const_iterator const_iterator;
-  
+  typedef typename Traits::Dimension D;
 private:
   Traits traits;
   // the iterator range of the Point_container 
@@ -53,20 +53,20 @@ private:
   boost::optional<iterator> m_e ;
   
   int built_coord;    // a coordinate for which the pointer list is built
-  Kd_tree_rectangle<FT> bbox;       // bounding box, i.e. rectangle of node
-  Kd_tree_rectangle<FT> tbox;       // tight bounding box, 
+  Kd_tree_rectangle<FT,D> bbox;       // bounding box, i.e. rectangle of node
+  Kd_tree_rectangle<FT,D> tbox;       // tight bounding box, 
   // i.e. minimal enclosing bounding
   // box of points
 	                	    
 public:
 
-  inline const Kd_tree_rectangle<FT>& 
+  inline const Kd_tree_rectangle<FT,D>& 
   bounding_box() const 
   { 
     return bbox; 
   }
   
-  inline const Kd_tree_rectangle<FT>&
+  inline const Kd_tree_rectangle<FT,D>&
   tight_bounding_box() const 
   { 
     return tbox; 
diff --git a/Spatial_searching/include/CGAL/Search_traits.h b/Spatial_searching/include/CGAL/Search_traits.h
index fb39d0d28ef..6af99f50e08 100644
--- a/Spatial_searching/include/CGAL/Search_traits.h
+++ b/Spatial_searching/include/CGAL/Search_traits.h
@@ -21,13 +21,17 @@
 
 #ifndef CGAL_KD_TREE_TRAITS_POINT_H
 #define CGAL_KD_TREE_TRAITS_POINT_H
+#include <CGAL\Dimension.h>
 
 
 namespace CGAL {
-  template <class FT_, class Point, class CartesianCoordinateIterator, class ConstructCartesianCoordinateIterator>
+  template <class FT_, class Point, class CartesianCoordinateIterator, class ConstructCartesianCoordinateIterator, int D>
   class Search_traits {
     
   public:
+
+    typedef typename Dimension_tag<D> Dimension;
+
     typedef CartesianCoordinateIterator Cartesian_const_iterator_d;
     typedef ConstructCartesianCoordinateIterator Construct_cartesian_const_iterator_d;
     typedef Point Point_d;
diff --git a/Spatial_searching/include/CGAL/Search_traits_2.h b/Spatial_searching/include/CGAL/Search_traits_2.h
index 666591a8ad3..07e80e3331f 100644
--- a/Spatial_searching/include/CGAL/Search_traits_2.h
+++ b/Spatial_searching/include/CGAL/Search_traits_2.h
@@ -21,6 +21,7 @@
 
 #ifndef CGAL_SEARCH_TRAITS_2_H
 #define CGAL_SEARCH_TRAITS_2_H
+#include <CGAL\Dimension.h>
 
 namespace CGAL {
 
@@ -30,6 +31,7 @@ namespace CGAL {
   class Search_traits_2 {
 
   public:
+    typedef typename Dimension_tag<2> Dimension;
     typedef typename K::Point_2 Point_d;
     typedef typename K::Iso_rectangle_2 Iso_box_d;
     typedef typename K::Circle_2 Sphere_d;
diff --git a/Spatial_searching/include/CGAL/Search_traits_3.h b/Spatial_searching/include/CGAL/Search_traits_3.h
index 9c6288d7825..8cd8bf3d49e 100644
--- a/Spatial_searching/include/CGAL/Search_traits_3.h
+++ b/Spatial_searching/include/CGAL/Search_traits_3.h
@@ -21,6 +21,7 @@
 
 #ifndef CGAL_SEARCH_TRAITS_3_H
 #define CGAL_SEARCH_TRAITS_3_H
+#include <CGAL\Dimension.h>
 
 namespace CGAL {
 
@@ -30,7 +31,8 @@ namespace CGAL {
   class Search_traits_3 {
 
   public:
-    
+
+   typedef typename Dimension_tag<3> Dimension;
     typedef typename K::Cartesian_const_iterator_3 Cartesian_const_iterator_d;
     typedef typename K::Construct_cartesian_const_iterator_3 Construct_cartesian_const_iterator_d;
     typedef typename K::Point_3 Point_d;
diff --git a/Spatial_searching/include/CGAL/Search_traits_d.h b/Spatial_searching/include/CGAL/Search_traits_d.h
index b0aba510310..8c2f3c31de0 100644
--- a/Spatial_searching/include/CGAL/Search_traits_d.h
+++ b/Spatial_searching/include/CGAL/Search_traits_d.h
@@ -21,13 +21,15 @@
 
 #ifndef CGAL_SEARCH_TRAITS_D_H
 #define CGAL_SEARCH_TRAITS_D_H
+#include <CGAL\Dimension.h>
 
 namespace CGAL {
 
-  template <class K>
+  template <class K,typename D = Dynamic_dimension_tag>
   class Search_traits_d {
 
   public:
+	  typedef typename D Dimension;
     typedef typename K::Cartesian_const_iterator_d Cartesian_const_iterator_d;
     typedef typename K::Construct_cartesian_const_iterator_d Construct_cartesian_const_iterator_d;
     typedef typename K::Construct_iso_box_d Construct_iso_box_d;