diff --git a/vcg/space/index/octree.h b/vcg/space/index/octree.h
new file mode 100644
index 00000000..9f944bfd
--- /dev/null
+++ b/vcg/space/index/octree.h
@@ -0,0 +1,552 @@
+/****************************************************************************
+* VCGLib                                                            o o     *
+* Visual and Computer Graphics Library                            o     o   *
+*                                                                _   O  _   *
+* Copyright(C) 2004                                                \/)\/    *
+* Visual Computing Lab                                            /\/|      *
+* ISTI - Italian National Research Council                           |      *
+*                                                                    \      *
+* All rights reserved.                                                      *
+*                                                                           *
+* This program is free software; you can redistribute it and/or modify      *   
+* it under the terms of the GNU General Public License as published by      *
+* the Free Software Foundation; either version 2 of the License, or         *
+* (at your option) any later version.                                       *
+*                                                                           *
+* This program is distributed in the hope that it will be useful,           *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+* GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
+* for more details.                                                         *
+*                                                                           *
+****************************************************************************/
+
+#ifndef OCTREE_H
+#define OCTREE_H
+
+#include <stdlib.h>
+#include <algorithm>
+#include <vector>
+#include <iterator>
+
+#include <vcg/space/color4.h>
+#include <vcg/space/index/base.h>
+#include <vcg/space/index/octree_template.h>
+#include <vcg/space/box3.h>
+#include <wrap/callback.h>
+
+namespace vcg
+{
+	/*!
+	* Given an object or an object pointer, return the reference to the object
+	*/
+	template <typename TYPE>
+	struct Dereferencer
+	{
+		static				TYPE& Ref(TYPE &t)				{ return ( t);	}
+		static				TYPE& Ref(TYPE* &t)				{ return (*t);	}
+		static const	TYPE& Ref(const TYPE &t)	{ return ( t);	}
+		static const	TYPE& Ref(const TYPE* &t) { return (*t);	}
+	};
+
+	
+	/*!
+	* Given a type, return the type
+	*/
+	template <typename T>
+	class ReferenceType
+	{
+	public:
+		typedef T Type;
+	};
+
+	
+	/*!
+	* Given as type a pointer to type, return the type
+	*/
+	template <typename T>
+	class ReferenceType<T *>
+	{
+	public:
+		typedef typename ReferenceType<T>::Type Type;
+	};
+
+
+	
+	/*!
+	* The type of the octree voxels
+	*/
+	struct Voxel
+	{
+		Voxel() { count = begin = end = -1; }
+
+		void SetRange(const int begin, const int end)
+		{
+			this->begin = begin;
+			this->end		= end;
+			count				= end-begin;
+		};
+
+		void AddRange(const Voxel *voxel)
+		{
+			assert(voxel->end>end);
+
+			count += voxel->count;
+			end		 = voxel->end;
+		};
+
+		int begin;
+		int end;
+		int count;
+	};
+	
+	
+	template < class OBJECT_TYPE, class SCALAR_TYPE>
+	class Octree : public OctreeTemplate< Voxel >, public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE >
+	{
+	public: 
+		typedef						OBJECT_TYPE													ObjectType;
+		typedef						vcg::Box3<ScalarType>								BoundingBoxType;	
+		typedef typename	Octree::Leaf											* LeafPointer;
+		typedef typename	Octree::InnerNode									* InnerNodePointer;
+		typedef typename	ReferenceType<OBJECT_TYPE>::Type	* ObjectPointer;
+
+		/*!
+		* Structure which holds the rendering settings
+		*/
+		struct OcreeRenderingSetting
+		{
+			OcreeRenderingSetting()
+			{
+				color						= vcg::Color4b(155, 155, 155, 255);
+				isVisible				= false;
+				minVisibleDepth	= 1;
+				maxVisibleDepth	= 4;
+			};
+
+			int						isVisible;
+			int						minVisibleDepth;
+			int						maxVisibleDepth;
+			vcg::Color4b	color;
+		};
+
+
+
+
+		
+	protected:
+		/***********************************************
+		*     INNER DATA STRUCTURES AND PREDICATES     *
+		***********************************************/
+		/*!
+		* Structure used during the sorting of the dataset
+		*/
+		template < typename LEAF_TYPE >
+		struct ObjectPlaceholder
+		{
+			typedef LEAF_TYPE* LeafPointer;
+
+			ObjectPlaceholder() { z_order = object_index = -1, leaf_pointer = NULL;}
+
+			ObjectPlaceholder(unsigned long long z_order, void* leaf_pointer, unsigned int object_index)
+			{
+				this->z_order				= z_order;
+				this->leaf_pointer	= leaf_pointer;
+				this->object_index	= object_index;
+			}
+
+			unsigned long long	z_order;
+			LeafPointer					leaf_pointer;
+			unsigned int				object_index;
+		};
+
+		
+		/*!
+		* Predicate used during the sorting of the dataset
+		*/
+		template <typename LEAF_TYPE >
+		struct ObjectSorter
+		{
+			inline bool operator()(const ObjectPlaceholder< LEAF_TYPE > &first, const ObjectPlaceholder< LEAF_TYPE > &second)
+			{
+				return (first.z_order<second.z_order);
+			}
+		};
+
+		/*!
+		* Structure which holds the reference to the object and the position of the mark for that object
+		*/
+		struct ObjectReference
+		{
+			ObjectReference() {mark_position=-1; pObject=NULL;}
+
+			int						mark_position;
+			ObjectPointer pObject;
+		};
+
+		/*
+		* The generic item in the neighbors vector computed by GetNearestNeighbors;
+		*/
+		struct Neighbour
+		{
+			Neighbour() 
+			{
+				this->object	 = NULL;
+				this->distance = -1.0f;
+			};
+
+			Neighbour(ObjectPointer &object, CoordType &point, float distance)
+			{
+				this->object	 = object;
+				this->point		 = point;
+				this->distance = distance;
+			}
+
+
+			ObjectPointer		object;
+			CoordType				point;
+			float						distance;
+		};
+
+		/*
+		* The operator used for sorting the items in neighbors based on the distances
+		*/
+		struct DistanceCompare
+		{
+			inline bool operator()( const Neighbour &p1, const Neighbour &p2) const { return p1.distance<p2.distance; }
+		}; //end of DistanceCompare
+
+
+
+public:
+		~Octree()
+		{
+			delete []marks;
+			int node_count = NodeCount();
+			for (int i=0; i<node_count; i++)
+				delete nodes[i];
+			nodes.clear();
+		}
+
+
+		/*!
+		* Populate the octree
+		*/
+		template <class OBJECT_ITERATOR, class BOUNDING_BOX_FUNCTOR>
+		void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const BoundingBoxType &bounding_box, const BOUNDING_BOX_FUNCTOR &bb_functor/*, vcg::CallBackPos *callback=NULL*/) 
+		{
+			typedef Dereferencer<typename ReferenceType<typename OBJECT_ITERATOR::value_type>::Type >	DereferencerType;
+
+			// Compute the bounding-box enclosing the whole dataset
+			BoundingBoxType resulting_bb(bounding_box);
+			CoordType offset = bounding_box.Dim()*Octree::EXPANSION_FACTOR;
+			CoordType center = bounding_box.Center();
+			resulting_bb.Offset(offset);
+			float longest_side = vcg::math::Max( resulting_bb.DimX(), vcg::math::Max(resulting_bb.DimY(), resulting_bb.DimZ()) )/2.0f;
+			resulting_bb.Set(center);
+			resulting_bb.Offset(longest_side);
+			boundingBox = resulting_bb;
+
+			// Try to find a reasonable octree depth 
+			int dataset_dimension = std::distance(bObj, eObj);
+
+			// Allocate the mark array
+			global_mark				= 1;
+			marks							= new unsigned char[dataset_dimension];
+			memset(&marks[0], 0, sizeof(unsigned char)*dataset_dimension);
+
+			int primitives_per_voxel;
+			int depth = 4;
+			do 
+			{
+				int		number_of_voxel = 1<<(3*depth); // i.e. 8^depth
+				float density					= float(number_of_voxel)/float(depth);
+				primitives_per_voxel	= int(float(dataset_dimension)/density);
+				depth++;
+			}
+			while (primitives_per_voxel>25 && depth<15);
+			Initialize(depth);
+
+			// Sort the dataset (using the lebesgue space filling curve...)
+			std::string message("Indexing dataset...");
+			Octree::NodePointer *route = new Octree::NodePointer[depth+1];
+			OBJECT_ITERATOR iObj = bObj;
+
+			//if (callback!=NULL) callback(int((i+1)*100/dataset_dimension), message.c_str());
+
+			std::vector< ObjectPlaceholder< Octree::Node > > placeholders/*(dataset_dimension)*/;
+			vcg::Box3<ScalarType>		object_bb;
+			vcg::Point3<ScalarType> hit_leaf;
+			for (int i=0; i<dataset_dimension; i++, iObj++)
+			{
+				object_bb = bb_functor(*iObj);
+				hit_leaf  = object_bb.min;
+
+				while (object_bb.IsIn(hit_leaf))
+				{
+					int placeholder_index = placeholders.size();
+					placeholders.push_back( ObjectPlaceholder< Octree::Node >() ); 
+					placeholders[placeholder_index].z_order			 = BuildRoute(hit_leaf, route);
+					placeholders[placeholder_index].leaf_pointer = route[depth];
+					placeholders[placeholder_index].object_index = i;					
+					
+					hit_leaf.X() += leafDimension.X();
+					if (hit_leaf.X()>object_bb.max.X())
+					{
+						hit_leaf.X() = object_bb.min.X();
+						hit_leaf.Z()+= leafDimension.Z();
+						if (hit_leaf.Z()>object_bb.max.Z())
+						{
+							hit_leaf.Z() = object_bb.min.Z();
+							hit_leaf.Y()+= leafDimension.Y();
+						}
+					}
+				}
+			}
+			delete []route;
+
+			int placeholder_count = int(placeholders.size());
+			std::sort(placeholders.begin(), placeholders.end(), ObjectSorter< Octree::Node >());
+			std::vector< Octree::Node* > filled_leaves(placeholder_count);
+			sorted_dataset.resize( dataset_dimension );
+			for (int i=0; i<placeholder_count; i++)
+			{
+				std::advance((iObj=bObj), placeholders[i].object_index);
+				sorted_dataset[i].pObject				= &DereferencerType::Ref(*iObj);
+				sorted_dataset[i].mark_position = i;
+				filled_leaves[i]								= placeholders[i].leaf_pointer;
+			}
+
+			// The dataset is sorted and the octree is built, but the indexing information aren't stored yet in the octree:
+			// we assign to each leaf the range inside the sorted dataset of the primitives contained inside the leaf
+			int begin									= -1;
+			NodePointer initial_leaf	= NULL;
+			for (int end=0; end<placeholder_count; )
+			{
+				begin = end;
+				initial_leaf = filled_leaves[begin];
+				do end++;
+				while (end<placeholder_count && initial_leaf==filled_leaves[end]);
+
+				VoxelType *voxel = Voxel(initial_leaf);
+				voxel->SetRange(begin, end);	
+			}
+
+			// The octree is built, the dataset is sorted but only the leaves are indexed: 
+			// we propaget the indexing information bottom-up to the root
+			IndexInnerNodes( Root() );
+		} //end of Set
+
+		/*!
+		* Retrive the k closest element to the query point
+		*/
+		template <class OBJECT_POINT_DISTANCE_FUNCTOR, class OBJECT_MARKER, class OBJECT_POINTER_CONTAINER, class DISTANCE_CONTAINER, class POINT_POINTER_CONTAINER>
+		unsigned int GetKClosest
+			(
+			OBJECT_POINT_DISTANCE_FUNCTOR &distance_functor, 
+			OBJECT_MARKER									&/*marker*/, 
+			const unsigned int						 k, 
+			const CoordType								&query_point, 
+			const ScalarType							&max_distance,
+			OBJECT_POINTER_CONTAINER			&objects, 
+			DISTANCE_CONTAINER						&distances, 
+			POINT_POINTER_CONTAINER				&points,
+			bool													 sort_per_distance   = false,
+			bool													 allow_zero_distance = false
+			)
+		{
+			// costruisco una bounging box centrata in query di dimensione pari a quella di una foglia.
+			// e controllo se in tale bounging box sono contenute un numero di elementi >= a k.
+			// Altrimenti espando il bounding box.
+			BoundingBoxType query_bb;
+			query_bb.Set(query_point);
+
+			// Il raggio della sfera centrata nel punto query
+			float sphere_radius = 0.0f;
+
+			// se il bounding-box non interseca il bounding-box dell'octree, lo espando subito
+			if (!query_bb.IsIn(query_point))
+			{
+				do
+				{ 
+					query_bb.Offset(leafDiagonal); 
+					sphere_radius += leafDiagonal;
+				}
+				while ( !boundingBox.Collide(query_bb) || sphere_radius>max_distance); // TODO check the termination condintion
+			}
+
+			std::vector< NodePointer > leaves;
+			std::vector< Neighbour	 > neighbors;
+
+			unsigned int object_count;
+			int					 leaves_count;
+			float				 k_distance = leafDiagonal;
+			do
+			{
+				// update the marks
+				global_mark = (global_mark+1)%255;
+				if (global_mark == 0)
+				{
+					memset(&marks[0], 0, sizeof(unsigned char)*int(sorted_dataset.size()));
+					global_mark++;
+				}
+
+				do
+				{
+					leaves.clear();
+					object_count = 0;
+					query_bb.Offset(k_distance);
+					sphere_radius += vcg::math::Max<float>(leafDiagonal, k_distance);
+
+					ContainedLeaves(query_bb, leaves, Root(), boundingBox);
+
+					leaves_count = int(leaves.size());
+					object_count = 0;
+					for (int i=0; i<leaves_count; i++)
+						object_count += Voxel( leaves[i] )->count;
+				}
+				while (object_count<k && sphere_radius>max_distance); // TODO check the termination condintion
+
+				// i punti contenuti nelle foglie sono sufficienti. 
+				// seleziono solamente i k punti più vicini.
+				CoordType				 closest_point;
+				VoxelType				*voxel;
+				ObjectReference *ref;
+				int							 begin;
+				int							 end;
+				float						 distance;
+
+				neighbors.clear();
+				for (int i=0; i<leaves_count; i++)
+				{
+					voxel = Voxel(leaves[i]);
+					begin = voxel->begin;
+					end		= voxel->end;
+					for ( ; begin<end; begin++)
+					{
+						ref				= &sorted_dataset[begin];
+						if (marks[ref->mark_position]==global_mark)
+							continue;
+						
+						distance = max_distance;
+						if (!distance_functor(*ref->pObject, query_point, distance, closest_point))
+							continue;
+
+						marks[ref->mark_position]=global_mark;
+						if ((distance!=0.0f || allow_zero_distance)/* && distance<max_distance*/)
+							neighbors.push_back( Neighbour(ref->pObject, closest_point, distance) );
+					} //end of for ( ; begin<end; begin++)
+				} // end of for (int i=0; i<leavesCount; i++)
+
+
+				if (sphere_radius<max_distance)
+				{
+					if (int(neighbors.size())<k)
+					{
+						k_distance = 2.0f*sphere_radius;
+						continue;
+					}
+					else
+						object_count = k;
+				}
+				else 
+					object_count=int(neighbors.size());
+
+				std::vector< Neighbour >::iterator first = neighbors.begin();
+				std::vector< Neighbour >::iterator last	 = neighbors.end();
+
+				if (sort_per_distance)  std::partial_sort(first, first+object_count, last, DistanceCompare());
+				else										std::nth_element (first, first+object_count, last, DistanceCompare()); 
+				k_distance = neighbors[object_count-1].distance;
+			} 
+			while (k_distance>sphere_radius && sphere_radius<max_distance);
+
+			// copy the nearest object into  
+			if (points.size()!=object_count) 
+			{
+				points.resize(object_count);
+				distances.resize(object_count);
+				objects.resize(object_count);
+			}
+
+			POINT_POINTER_CONTAINER::iterator  iPoint			= points.begin();
+			DISTANCE_CONTAINER::iterator			 iDistance	= distances.begin();
+			OBJECT_POINTER_CONTAINER::iterator iObject		= objects.begin();
+			for (unsigned int n=0; n<object_count; n++, iPoint++, iDistance++, iObject++)
+			{
+				(*iPoint)		 = neighbors[n].point;
+				(*iDistance) = neighbors[n].distance;
+				(*iObject)	 = neighbors[n].object;
+			}
+			return object_count;
+		}; //end of GetKClosest
+
+
+		/*!
+		* Draw the octree in a valid OpenGL context according to the rendering settings
+		*/
+		void Octree::DrawOctree(vcg::Box3f &boundingBox, NodePointer n)
+		{
+			char level = Level(n);
+			NodePointer son;
+			if (rendering_settings.minVisibleDepth>level)
+			{
+				for (int s=0; s<8; s++)
+					if ((son=Son(n, s))!=0)
+						DrawOctree(SubBox(boundingBox, s), son);
+			}
+			else
+			{
+				vcg::glBoxWire(boundingBox);
+				if (level<rendering_settings.maxVisibleDepth)
+					for (int s=0; s<8; s++)
+						if ((son=Son(n, s))!=0)
+							DrawOctree(SubBox(boundingBox, s), son);
+			}
+		};
+
+	protected:
+		/*!
+		* Contains pointer to the objects in the dataset. 
+		* The pointers are sorted so that the object pointed to result ordered in the space
+		*/
+		std::vector< ObjectReference >	sorted_dataset;
+
+		/*!
+		* Markers used to avoid duplication of the same result during a query
+		*/
+		unsigned char	*marks;
+		unsigned char  global_mark;
+
+	public:
+		OcreeRenderingSetting					rendering_settings;
+
+	protected:
+		/*!
+		* When all the leaves are indexed, this procedure is called in order to propagate the indexing information to the inner nodes
+		*/
+		void IndexInnerNodes(NodePointer n)
+		{
+			assert(n!=NULL);
+
+			VoxelPointer current_voxel = Voxel(n);
+			VoxelPointer son_voxel;
+			for (int s=0; s<8; s++)
+			{
+				NodePointer son_index = Son(n, s);
+				if (son_index!=NULL)
+				{
+					if (Level(son_index)!=maximumDepth)
+						IndexInnerNodes(son_index);
+
+					son_voxel = Voxel(son_index);
+					current_voxel->AddRange( son_voxel );
+				}
+			}
+		}; // end of IndexInnerNodes
+	};
+
+
+} //end of namespace vcg
+
+#endif //OCTREE_H
diff --git a/vcg/space/index/octree_template.h b/vcg/space/index/octree_template.h
new file mode 100644
index 00000000..8bc97a93
--- /dev/null
+++ b/vcg/space/index/octree_template.h
@@ -0,0 +1,701 @@
+/****************************************************************************
+* VCGLib                                                            o o     *
+* Visual and Computer Graphics Library                            o     o   *
+*                                                                _   O  _   *
+* Copyright(C) 2004                                                \/)\/    *
+* Visual Computing Lab                                            /\/|      *
+* ISTI - Italian National Research Council                           |      *
+*                                                                    \      *
+* All rights reserved.                                                      *
+*                                                                           *
+* This program is free software; you can redistribute it and/or modify      *   
+* it under the terms of the GNU General Public License as published by      *
+* the Free Software Foundation; either version 2 of the License, or         *
+* (at your option) any later version.                                       *
+*                                                                           *
+* This program is distributed in the hope that it will be useful,           *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+* GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
+* for more details.                                                         *
+*                                                                           *
+****************************************************************************/
+
+
+#ifndef OCTREETEMPLATE_H
+#define OCTREETEMPLATE_H
+
+#include <vcg/space/point3.h>
+#include <vcg/space/box3.h>
+#include <vector>
+
+
+namespace vcg
+{
+
+/* Octree Template
+
+Tiene un dataset volumetrico come un octree
+Assunzione che la grandezza sia una potenza di due
+La prof max e' fissa.
+E' un octree in cui il dato e' nella cella dell'octree.
+Anche i nodi non foglia hanno il dato Voxel
+
+Assunzioni sul tipo voxel:
+che abbia definiti gli operatori per poterci fare sopra pushpull.
+
+Si tiene int invece di puntatori per garantirsi reallocazione dinamica.
+I dati veri e propri stanno in un vettore di nodi
+*/
+
+template <typename VOXEL_TYPE>
+class OctreeTemplate
+{
+protected:
+	struct Node;
+
+public:
+	// Octree Type Definitions
+	typedef typename VOXEL_TYPE		VoxelType;
+	typedef typename VOXEL_TYPE  *VoxelPointer;
+	typedef vcg::Point3i					CenterType;	
+	static const float						EXPANSION_FACTOR;
+	typedef	int										NodeIndex;
+	typedef Node								 *NodePointer;
+
+protected:
+	/*
+	* Inner structures:
+	* Contains the information related to the octree node
+	*/
+	struct Node
+	{
+		// Default constructor: fill the data members with non-meaningful values
+		Node()
+		{
+			parent = NULL;
+			level  = -1;
+		}
+
+		// Constructor: create a new Node
+		Node(NodePointer parent, int level)
+		{
+			this->parent	= parent;
+			this->level		= (char) level;
+		}
+
+		inline virtual NodePointer &Son(int sonIndex) = 0;
+
+		inline virtual bool	IsLeaf() = 0;
+
+		// The position of the center of the node in integer coords in the 0..2^(2*sz) -1 range
+		// The root has position (lsz/2,lsz/2,lsz/2)
+		CenterType	center;   
+		char				level;	
+		NodePointer parent;
+		VoxelType		voxel;
+	};
+
+	/*
+	* Inner struct: Node 
+	*/
+	struct InnerNode : public Node
+	{
+		InnerNode() : Node() {};
+		InnerNode(NodePointer parent, int level) : Node(parent, level) 
+		{
+			memset(&sons[0], NULL, 8*sizeof(Node*));
+		}
+
+		inline NodePointer &Son(int sonIndex)
+		{
+			assert(0<=sonIndex && sonIndex<=8);
+			return sons[sonIndex];
+		}
+
+		inline bool IsLeaf()
+		{
+			return false;
+		}
+
+		NodePointer sons[8];
+	};
+
+	/*
+	* Inner struct: Leaf
+	*/
+	struct Leaf : public Node
+	{
+		Leaf() : Node() {};
+		Leaf(NodePointer parent, int level) : Node(parent, level) {}
+
+		inline NodePointer &Son(int /*sonIndex*/)
+		{
+			assert(false);
+			return parent;
+		}
+
+		inline bool IsLeaf()
+		{
+			return true;
+		}
+	};
+
+public:
+	// Inizializza l'octree
+	void Initialize(int maximumDepth)
+	{
+		this->maximumDepth	= maximumDepth;
+		size	= 1<< maximumDepth;			// e.g. 1*2^maxDepth
+		lSize = 1<<(maximumDepth+1);	// e.g. 1*2^(maxDepth+1)
+		
+		InnerNode *root = new InnerNode(NULL,0);
+		nodes.clear();
+		nodes.push_back( root );
+		root->center	= CenterType(size, size, size);
+		
+		float szf = (float) size;
+		leafDimension  = boundingBox.Dim();
+		leafDimension /= szf;
+		leafDiagonal	 = leafDimension.Norm();
+	};
+
+	// Return the octree bounding-box
+	inline vcg::Box3f		BoundingBox()											{ return boundingBox;				}
+	
+	// Return the Voxel of the n-th node
+	inline VoxelPointer	Voxel(const NodePointer n)				{ return &(n->voxel);				}
+
+	// Return the octree node count
+	inline int					NodeCount()									const { return int(nodes.size()); }
+
+	// Return the root index
+	inline const				NodePointer Root()					const	{ return nodes[0];					}
+
+	// Return the level of the n-th node
+	inline char					Level(const NodePointer n)	const	{ return n->level;					}
+
+	// Return the referente to the i-th son of the n-th node
+	inline NodePointer&	Son(NodePointer n, int i)		const	{ return n->Son(i);					}
+
+	// Return the parent index of the n-th node
+	inline NodePointer	Parent(const NodePointer n) const	{ return n->parent;					}
+
+	// Return the index of the current node in its father
+	int WhatSon(NodePointer n) const
+	{
+		if(n==Root()) 
+			assert(false);
+		
+		NodePointer parent = Parent(n);
+		for(int i=0;i<8;++i)
+			if(parent->Son(i)==n) 
+				return i;
+		
+		return -1;
+	}
+
+	// Return the center of the n-th node 
+	inline CenterType CenterInOctreeCoordinates(const NodePointer n) const { return n->center;}
+	
+	// Return the center of the n-th node expressed in world-coordinate
+	inline vcg::Point3f CenterInWorldCoordinates(const NodePointer n) const
+	{
+		assert(0<=n && n<NodeCount());
+
+		int level = n->level;
+		int shift = maxDepth-level+1;
+		vcg::Point3f wcCenter; 
+		vcg::Point3f ocCenter = CenterInOctreeCoordinates(n);
+		vcg::Point3f nodeSize = boundingBox.Dim()/float(1<<level);
+		wcCenter.X() = boundingBox.min.X() + (nodeSize.X()*(0.5f+(ocCenter.X()>>shift)));
+		wcCenter.Y() = boundingBox.min.Y() + (nodeSize.Y()*(0.5f+(ocCenter.Y()>>shift)));
+		wcCenter.Z() = boundingBox.min.Z() + (nodeSize.Z()*(0.5f+(ocCenter.Z()>>shift)));
+		return wcCenter
+	}
+
+	// Given a node (even not leaf) it returns the center of the box it represent.
+	// the center is expressed not in world-coordinates.
+	// e.g. the root is (sz/2,sz/2,sz/2);
+	// and the finest element in the grid in lower left corner has center (.5, .5, .5)
+	/*
+
+	4----------------   4----------------   4----------------
+	|   |   |   |   |   |       |       |   |               |
+	3---+---+---+---|   3       |       |   3               |
+	|   |   |   |   |   |       |       |   |               |
+	2---+---+---+---|   2---+---+---+---|   2       c       |
+	|   |   |   |   |   |       |       |   |               |
+	1---+---+---+---|   1   b   +       |   1               |
+	| a |   |   |   |   |       |       |   |               |
+	0---1---2---3---4   0---1---2---3---4   0---1---2---3---4
+
+	This is a tree with maxdepth==2, so sz is 2^2=4
+
+	a) a leaf at the deepest level 2 has position (.5,.5)
+	b) a mid node (lev 1) has position (1,1)
+	c) root has level 0 and position (sz/2,sz/2) = (2,2) 
+
+	The center of a node has integer coords in the 2^(MaxDepth+1) range.
+
+	The other approach is to use position as a bit string 
+	codifying the tree path, but in this case you have to 
+	supply also the level (e.g. the string lenght) 
+	you desire. The lower left corner node is always 0 ( (,) for the root (0,0) level 1, and (00,00) for level 2)
+
+	|              ~~~              |
+	|      0~~      |      1~~      |
+	|  00~  |  01~  |  10~  |  11~  |
+	|000|001|010|011|100|101|110|111|      
+
+	The interesting properties is that
+	if your octree represent a space [minv,maxv] and you want 
+	to find the octree cell containing a point p in [minv,maxv] 
+	you just have to convert p in the range [0,sz) truncate it to an integer and use it as a path.
+	For example, consider an octree of depth 3, representing a range [0..100)
+	sz=8 (each cell contains form 0 to 12.5
+	the point 
+	5		-> 0.4	-> path is 000
+	45	-> 3.6	-> path is 011 
+	50  -> 4.0	-> path is 100  
+	100 -> 8		-> ERROR the interval is right open!!!
+
+	Note how each cell is meant to contains a right open interval (e.g. the first cell contains [0,12.5) and the second [12.5,25) and so on)
+
+	The center of each cell can simply be obtained by adding .5 to the path of the leaves.
+	*/
+	vcg::Point3f Center(NodePointer n) const
+	{
+		vcg::Point3f center;
+		center.Import(GetPath(n));
+		center+=Point3f(.5f,.5f,.5f);
+
+		//TODO verify the assert
+		assert(center==nodes[n]->center);
+
+		return center;
+	}
+
+	// Return the bounding-box of the n-th node expressed in world-coordinate
+	vcg::Box3f BoundingBoxInWorldCoordinates(const NodePointer n)
+	{
+		assert(0<=n && n<Size());
+
+		char level = Level(n);
+		int shift  = maximumDepth-level+1;
+		vcg::Point3f nodeDim = boundingBox.Dim()/float(1<<level);
+		CenterType   center	 = CenterInOctreeCoordinates(n);
+		vcg::Box3f	 nodeBB;
+		nodeBB.min.X() = boundingBox.min.X() + (nodeDim.X()*(center.X()>>shift));
+		nodeBB.min.Y() = boundingBox.min.Y() + (nodeDim.Y()*(center.Y()>>shift));
+		nodeBB.min.Z() = boundingBox.min.Z() + (nodeDim.Z()*(center.Z()>>shift));
+		nodeBB.max = nodeBB.min+nodeDim;
+		return nodeBB;
+	};
+
+	// Return one of the 8 subb box of a given box.
+	vcg::Box3f SubBox(vcg::Box3f &lbb, int i)
+	{
+		vcg::Box3f bs;
+		if (i&1)	bs.min.X()=(lbb.min.X()+(bs.max.X()=lbb.max.X()))/2.0f;
+		else			bs.max.X()=((bs.min.X()=lbb.min.X())+lbb.max.X())/2.0f;
+		if (i&2)	bs.min.Y()=(lbb.min.Y()+(bs.max.Y()=lbb.max.Y()))/2.0f;
+		else			bs.max.Y()=((bs.min.Y()=lbb.min.Y())+lbb.max.Y())/2.0f;
+		if (i&4)	bs.min.Z()=(lbb.min.Z()+(bs.max.Z()=lbb.max.Z()))/2.0f;
+		else			bs.max.Z()=((bs.min.Z()=lbb.min.Z())+lbb.max.Z())/2.0f;
+		
+		return bs;
+	}
+
+	// Given the bounding-box and the center (both in world-coordinates)
+	// of a node, return the bounding-box (in world-coordinats) of the i-th son 
+	vcg::Box3f SubBoxAndCenterInWorldCoordinates(vcg::Box3f &lbb, vcg::Point3f &center, int i)
+	{
+		vcg::Box3f bs;
+		if (i&1)	
+		{
+			bs.min[0]=center[0];
+			bs.max[0]=lbb.max[0];
+		}
+		else			
+		{
+			bs.min[0]=lbb.min[0];
+			bs.max[0]=center[0];
+		}
+		if (i&2)
+		{
+			bs.min[1]=center[1];
+			bs.max[1]=lbb.max[1];
+		}
+		else			
+		{
+			bs.max[1]=center[1];
+			bs.min[1]=lbb.min[1];
+		}
+		if (i&4)	
+		{
+			bs.min[2]=center[2];
+			bs.max[2]=lbb.max[2];
+		}
+		else			
+		{
+			bs.max[2]=center[2];
+			bs.min[2]=lbb.min[2];
+		}
+		return bs;
+	};
+
+	/*
+	* Add a new Node to the octree. 
+	*	The created node is the i-th son of the node pointed to by parent.
+	*	Return the pointer to the new node
+	*/
+	NodePointer NewNode(NodePointer parent, int i)
+	{
+		assert(0<=i && i<8);
+		assert(Son(parent, i)==NULL);
+
+		//int index  = NodeCount();
+		char level = Level(parent)+1;
+				
+		Node *node = (level<maximumDepth)? (Node*) new InnerNode(parent, level) : (Node*) new Leaf(parent, level);
+		nodes.push_back( node );
+		Son(parent, i) = node;
+		
+		CenterType *parentCenter = &(parent->center);
+		int displacement = 1<<(maximumDepth-level);
+		node->center.X() = parentCenter->X() + ((i&1)? displacement : -displacement);
+		node->center.Y() = parentCenter->Y() + ((i&2)? displacement : -displacement);
+		node->center.Z() = parentCenter->Z() + ((i&4)? displacement : -displacement);
+
+		return node;
+	}
+
+	// Aggiunge un nodo all'octree nella posizione specificata e al livello specificato.
+	// Vengono inoltre inseriti tutti gli antenati mancanti per andare dalla radice
+	// al nodo ed al livello specificato seguendo path.
+	NodePointer AddNode(CenterType path)
+	{
+		//the input coordinates must be in the range 0..2^maxdepth
+		assert(path[0]>=0 && path[0]<size);    
+		assert(path[1]>=0 && path[1]<size);    
+		assert(path[2]>=0 && path[2]<size);
+
+		NodePointer curNode	= Root();
+		int rootLevel				= 0;
+		int shiftLevel			= maximumDepth-1;
+		
+		while(shiftLevel >= rootLevel)
+		{
+			int nextSon=0;
+			if((path[0]>>shiftLevel)%2) nextSon +=1; 
+			if((path[1]>>shiftLevel)%2) nextSon +=2; 
+			if((path[2]>>shiftLevel)%2) nextSon +=4; 
+			NodePointer nextNode = Son(curNode, nextSon);	
+			if(nextNode!=NULL) // nessun nodo può aver Root() per figlio 
+				curNode = nextNode;
+			else 
+			{
+				NodePointer newNode = NewNode(curNode, nextSon);
+				assert(Son(curNode, nextSon)==newNode); // TODO delete an assignment
+				curNode=newNode;
+			}
+			--shiftLevel;
+		}
+		return curNode;
+	}
+
+	// Convert the point p coordinates to the integer based representation 
+	// in the range 0..size, where size is 2^maxdepth
+	vcg::Point3i Interize(const vcg::Point3f &pf) const  
+	{
+		vcg::Point3i pi;
+
+		assert(pf.X()>=boundingBox.min.X() &&  pf.X()<=boundingBox.max.X());
+		assert(pf.Y()>=boundingBox.min.Y() &&  pf.Y()<=boundingBox.max.Y());
+		assert(pf.Z()>=boundingBox.min.Z() &&  pf.Z()<=boundingBox.max.Z());
+
+		pi.X() = int((pf.X() - boundingBox.min.X()) * size / (boundingBox.max.X() - boundingBox.min.X()));
+		pi.Y() = int((pf.Y() - boundingBox.min.Y()) * size / (boundingBox.max.Y() - boundingBox.min.Y()));
+		pi.Z() = int((pf.Z() - boundingBox.min.Z()) * size / (boundingBox.max.Z() - boundingBox.min.Z()));
+		
+		return pi;
+	}
+
+	// Inverse function of Interize; 
+	// Return to the original coords space (not to the original values!!)
+	vcg::Point3f DeInterize(const vcg::Point3i &pi ) const 
+	{
+		vcg::Point3f pf;
+
+		assert(pi.X()>=0 && pi.X()<size); 
+		assert(pi.Y()>=0 && pi.Y()<size); 
+		assert(pi.Z()>=0 && pi.Z()<size); 
+
+		pf.X() = vi.X() * (boundingBox.max.X() - boundingBox.min.X()) / size + boundingBox.min.X();
+		pf.Y() = vi.Y() * (boundingBox.max.Y() - boundingBox.min.Y()) / size + boundingBox.min.Y();
+		pf.Z() = vi.Z() * (boundingBox.max.Z() - boundingBox.min.Z()) / size + boundingBox.min.Z();
+		
+		return pf;
+	}
+
+	// Compute the z-ordering integer value for a given node;
+	// this value can be used to compute a complete ordering of the nodes of a given level of the octree. 
+	// It assumes that the octree has a max depth of 10.
+	unsigned long long ZOrder(NodePointer n)
+	{
+		unsigned long long finalPosition = 0;
+		unsigned long long currentPosition;
+		char level				= Level(n);
+
+		CenterType path = GetPath(n); 
+		
+		for(int i=0; i<level; ++i)
+		{
+			currentPosition = 0;
+			int mask=1<<i;
+			if(path[0]&mask) currentPosition|=1;
+			if(path[1]&mask) currentPosition|=2;
+			if(path[2]&mask) currentPosition|=4;
+			currentPosition = currentPosition<<(i*3);
+			finalPosition	 |= currentPosition;
+		}
+		return finalPosition;
+	}
+
+	// Funzione principale di accesso secondo un path; 
+	// restituisce l'indice del voxel di profondita' massima 
+	// che contiene il punto espresso in range 0..2^maxk	
+	NodePointer DeepestNode(CenterType path, int MaxLev)
+	{
+		assert(path[0]>=0 && path[0]<sz);
+		assert(path[1]>=0 && path[1]<sz);
+		assert(path[2]>=0 && path[2]<sz);
+
+		NodePointer curNode	= Root();
+		int shift						= maximumDepth-1;
+		
+		while(shift && Level(curNode) < MaxLev)
+		{
+			int son = 0;
+			if((path[0]>>shift)%2) son +=1; 
+			if((path[1]>>shift)%2) son +=2; 
+			if((path[2]>>shift)%2) son +=4; 
+			NodePointer nextNode = Son(curNode, son);
+			if(nextNode!=NULL) 
+				CurNode=nextNode;
+			else 
+				break;
+
+			--shift;
+		}
+		return CurNode;
+	}
+
+
+	// Return the 'path' from root to the specified node;
+	// the path is coded as a point3s; each bit of each component code the direction in one level
+	// only the last 'level' bits of the returned value are meaningful
+	// for example for the root no bit are meaningfull (path is 0)
+	// for the first level only one bit of each one of the three components are maninguful;
+	CenterType GetPath(NodePointer n) const
+	{
+		if(n==Root()) 
+			return CenterType(0,0,0);
+		
+		CenterType	path(0,0,0);
+
+		int shift, mask, son;
+		int startingLevel = int(Level(n));
+		while (n!=Root())
+		{
+			shift = startingLevel-Level(n); //nodes[n].level
+			mask = 1 << shift;							// e.g. 1*2^shift
+			son = WhatSon(n);
+			if(son&1) path[0] |= mask;
+			if(son&2) path[1] |= mask;
+			if(son&4) path[2] |= mask;
+			n = Parent(n);									// nodes[n].parent
+		}       
+		return path;
+	}
+
+	// Dato un punto 3D nello spazio restituisce un array contenente 
+	// i puntatori ai nodi che lo contengono, dalla radice fino alle foglie. 
+	// I nodi mancanti dalla radice fino a profondità maxDepth vengono aggiunti. 
+	// In posizione i ci sarà il nodo di livello i.
+	// Restituisce lo z-order del punto p
+	unsigned long long BuildRoute(const vcg::Point3f &p, NodePointer *&route)
+	{
+		assert( boundingBox.min.X()<=p.X() && p.X()<=boundingBox.max.X() );
+		assert( boundingBox.min.Y()<=p.Y() && p.Y()<=boundingBox.max.Y() );
+		assert( boundingBox.min.Z()<=p.Z() && p.Z()<=boundingBox.max.Z() );
+
+		route[0]						= Root();
+		NodePointer curNode = Root();
+		int finalLevel			= 0;
+		int shift						= maximumDepth-1;
+
+		CenterType path = CenterType::Construct(Interize(p));
+		while(shift >= finalLevel)
+		{
+			int son = 0;
+			if((path[0]>>shift)%2) son +=1; 
+			if((path[1]>>shift)%2) son +=2; 
+			if((path[2]>>shift)%2) son +=4; 
+			NodePointer nextNode = Son(curNode, son);
+
+			if(nextNode!=NULL)
+			{
+				route[maximumDepth-shift] = nextNode;
+				curNode										= nextNode;
+			}
+			else 
+			{
+				NodePointer newNode = NewNode(curNode, son);
+				route[maximumDepth-shift] = newNode;
+				curNode = newNode;
+			}
+			--shift;
+		}
+		return ZOrder(route[maximumDepth]);
+	}; //end of BuildRoute
+
+
+	// Restituisce il percorso dalla radice fino al nodo di profondità
+	// massima presente nell'octree contenente il nodo p. Nessun nuovo nodo è aggiunto
+	// all'octree. In route sono inseriti gli indici dei nodi contenti p, dalla radice
+	// fino al nodo di profontidà massima presente; nelle eventuali posizioni rimaste
+	// libere è inserito il valore -1. Restituisce true se il punto p cade in una foglia
+	// dell'otree, false altrimenti
+	bool GetRoute(const vcg::Point3f &p, NodePointer *&route)
+	{
+		assert( boundingBox.min.X()<=p.X() && p.X()<=boundingBox.max.X() );
+		assert( boundingBox.min.Y()<=p.Y() && p.Y()<=boundingBox.max.Y() );
+		assert( boundingBox.min.Z()<=p.Z() && p.Z()<=boundingBox.max.Z() );
+
+		memset(route, NULL, maxDepth*sizeof(NodePointer));
+
+		CenterType path				 = CenterType::Construct(Interize(p));
+		int shift							 = maxDepth-1;
+		NodePointer finalLevel = Root();
+		NodePointer curNode		 = Root();
+		
+		while(shift >= finalLevel)
+		{
+			int son=0;
+			if((path[0]>>shift)%2) son +=1; 
+			if((path[1]>>shift)%2) son +=2; 
+			if((path[2]>>shift)%2) son +=4; 
+			NodePointer nextNode = Son(curNode, son); 
+			if(nextNode!=NULL) 
+			{
+				route[maxDepth-shift] = nextNode;
+				curNode = nextNode;
+			}
+			else 
+				return false;
+			
+			--shift;
+		}
+		return true;
+	}; //end of GetReoute
+
+	// Data una bounding-box bb_query, calcola l'insieme dei nodi di 
+	// profondità depth il cui bounding-box ha intersezione non nulla con 
+	// bb (la bounding-box dell'octree); i puntatori a tali nodi sono 
+	// inseriti progressivamente in contained_nodes.
+	// The vector nodes must be cleared before calling this method.
+	void ContainedNodes(
+		vcg::Box3f								 &query,
+		std::vector< NodePointer > &nodes,
+		int													depth, 
+		NodePointer									n,
+		vcg::Box3f								 &nodeBB)
+	{
+		if (!query.Collide(nodeBB))
+			return;
+
+		if (Level(n)==depth)
+			nodes.push_back(n);
+		else
+		{
+			NodePointer	son;
+			vcg::Box3f	sonBB;
+			vcg::Point3f nodeCenter = nodeBB.Center();
+			for (int s=0; s<8; s++)
+			{
+				son = Son(n, s);
+				if (son!=NULL)
+				{
+					sonBB = SubBoxAndCenterInWorldCoordinates(nodeBB, nodeCenter, s);
+					ContainedNodes(query, nodes, depth, son, sonBB);
+				}
+			}
+		}
+	}; //end of ContainedNodes
+	
+
+	// Data una bounding-box bb, calcola l'insieme delle foglie il cui 
+	// bounding-box ha intersezione non nulla con bb; i loro indici
+	// sono inseriti all'interno di leaves.
+	void ContainedLeaves(
+		vcg::Box3f								 &query, 
+		std::vector< NodePointer > &leaves,
+		NodePointer									node,
+		vcg::Box3f								 &nodeBB
+		)
+	{
+		NodePointer	son;
+		vcg::Box3f	sonBB;
+		vcg::Point3f nodeCenter = nodeBB.Center();
+		for (int s=0; s<8; s++)
+		{
+			son = Son(node, s); //nodes[nodeIndex].sonIndex[s]
+			if (son!=NULL)
+			{
+				sonBB = SubBoxAndCenterInWorldCoordinates(nodeBB, nodeCenter, s);
+				if ( query.Collide(sonBB) )
+				{
+					if ( son->IsLeaf() )
+						leaves.push_back(son);
+					else
+						ContainedLeaves(query, leaves, son, sonBB);
+				}
+			}
+		}
+	}; //end of ContainedLeaves
+
+	
+	/*
+	* Octree Data Members
+	*/
+public:
+	// the size of the finest grid available (2^maxDepth)
+	int						size;
+
+	// double the size(2^maxDepth)
+	int						lSize;
+
+	// The allowed maximum depth
+	int						maximumDepth;
+
+	// The dimension of a leaf
+	vcg::Point3f	leafDimension; 
+
+	// The diagonal of a leaf
+	float					leafDiagonal;
+	
+	// The Octree nodes
+	std::vector< Node* > nodes;
+
+	// The bounding box containing the octree (in world coordinate)
+	vcg::Box3f				boundingBox;
+}; //end of class OctreeTemplate
+
+template <typename VOXEL_TYPE>
+const float OctreeTemplate<VOXEL_TYPE>::EXPANSION_FACTOR  = 0.035f;
+
+}
+
+#endif //OCTREETEMPLATE_H