diff --git a/vcg/space/index/perfect_spatial_hashing.h b/vcg/space/index/perfect_spatial_hashing.h
new file mode 100644
index 00000000..a7798d5a
--- /dev/null
+++ b/vcg/space/index/perfect_spatial_hashing.h
@@ -0,0 +1,1599 @@
+/****************************************************************************
+* 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 VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
+#define VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
+
+#pragma warning(disable : 4996)
+
+#define _USE_GRID_UTIL_PARTIONING_ 1
+#define _USE_OCTREE_PARTITIONING_  (1-_USE_GRID_UTIL_PARTIONING_)
+
+#include <vector>
+#include <list>
+#include <algorithm>
+
+#include <vcg/space/index/base.h>
+#include <vcg/space/index/grid_util.h>
+
+#include <vcg/space/point2.h>
+#include <vcg/space/point3.h>
+#include <vcg/space/box3.h>
+
+namespace vcg
+{
+
+	// Compute the greatest common divisor between two integers a and b
+	int GreatestCommonDivisor(const int a, const int b)
+	{
+		int m = a;
+		int n = b;
+
+		do 
+		{
+			if (m<n) std::swap(m, n);
+			m = m % n;
+			std::swap(m, n);
+		}
+		while (n!=0);
+		return m;
+	}
+	
+	// Doxygen documentation
+	/** \addtogroup index */
+	/*! @{ */ 
+
+	/*!
+	 * This class implements the perfect spatial hashing by S.Lefebvre and H.Hoppe
+	 * This is an spatial indexing structure such as the uniform grid, but with lower
+	 * memory requirements, since all the empty cells of the uniform grid are removed.
+	 * Access to a non-empty cell is performed looking up in two d-dimensional tables,
+	 * the offset table and the hash table.
+	 * @param OBJECT_TYPE (Template parameter) the type of objects to be indexed
+	 * @param SCALAR_TYPE (Template parameter) the scalar type
+	 */
+	template < class OBJECT_TYPE, class SCALAR_TYPE >
+	class PerfectSpatialHashing : public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE >
+	{
+		// Given an object or a pointer to an object, return the reference to the object
+		template < typename TYPE >
+		struct Dereferencer
+		{
+			static				TYPE& Reference(TYPE				 &t)	{	return  t;	}
+			static				TYPE& Reference(TYPE*			 &t)	{ return *t;  }
+			static const  TYPE& Reference(const TYPE	 &t)	{ return  t;	}
+			static const	TYPE& Reference(const TYPE* &t)	{ return *t;	}
+		};
+
+		// Given a type, holds this type in Type 
+		template < typename TYPE >
+		struct ReferenceType { typedef TYPE Type; };
+
+		// Given as type a "pointer to type", holds the type in Type
+		template < typename TYPE >
+		struct ReferenceType< TYPE * > { typedef typename ReferenceType<TYPE>::Type Type; };
+
+	public:
+		typedef						SCALAR_TYPE													ScalarType;
+		typedef						OBJECT_TYPE													ObjectType;
+		typedef typename	ReferenceType< ObjectType >::Type * ObjectPointer;
+		typedef typename  vcg::Box3< ScalarType >							BoundingBoxType;
+		typedef typename	vcg::Point3< ScalarType >						CoordinateType;
+		
+	protected:
+		/*! \struct NeighboringEntryIterator
+		* This class provides a convenient way to iterate over the six neighboring cells of a given cell.
+		*/
+		struct NeighboringEntryIterator
+		{
+			/*!
+			* Default constructor. 
+			* @param[in] entry				The index of the cell in the UniformGrid around which iterate.
+			* @param[in] table_size  The number of cells in the UniformGrid for each side.
+			*/
+			NeighboringEntryIterator(const vcg::Point3i &entry, const int table_size)
+			{
+				m_Center		= entry;
+				m_TableSize = table_size;
+				m_CurrentNeighbor.X() = (m_Center.X()+m_TableSize-1)%m_TableSize;
+				m_CurrentNeighbor.Y() = m_Center.Y();
+				m_CurrentNeighbor.Z() = m_Center.Z();
+				m_CurrentIteration		= 0;
+			}
+
+			/*!
+			* Increment the iterator to point to the next neighboring entry in the UniformGrid
+			*/
+			void operator++(int)
+			{
+				switch(++m_CurrentIteration)
+				{
+				case 1: m_CurrentNeighbor.X()=(m_Center.X()+1)%m_TableSize; break;
+				case 2: m_CurrentNeighbor.X()=m_Center.X(); m_CurrentNeighbor.Y()=(m_Center.Y()+m_TableSize-1)%m_TableSize; break;
+				case 3: m_CurrentNeighbor.Y()=(m_Center.Y()+1)%m_TableSize; break;
+				case 4: m_CurrentNeighbor.Y()=m_Center.Y(); m_CurrentNeighbor.Z()=(m_Center.Z()+m_TableSize-1)%m_TableSize; break;
+				case 5: m_CurrentNeighbor.Z()=(m_Center.Z()+1)%m_TableSize; break;
+				default: m_CurrentNeighbor = vcg::Point3i(-1, -1, -1); break;
+				}
+			}
+
+			/*!
+			* Dereferencing operator
+			* \return The neighbor of the given cell at the current iteration
+			*/
+			vcg::Point3i operator*() { return m_CurrentNeighbor; }
+
+			/*!
+			* Assignment operator
+			* @param[in] The source neighboring iterator
+			* \return		The reference to this iterator
+			*/
+			NeighboringEntryIterator& operator =(const NeighboringEntryIterator &it)
+			{
+				m_Center						= it.m_Center						;
+				m_CurrentNeighbor		= it.m_CurrentNeighbor	;
+				m_CurrentIteration	= it.m_CurrentIteration	;
+				m_TableSize					= it.m_TableSize				;
+				return *this;
+			}
+
+			/*!
+			* Less than operator. Since each entry in the UniformGrid has only 6 neighbors, 
+			* the iterator over the neighboring entries can be compared with an integer.
+			*/
+			inline bool operator <(const int value) { return m_CurrentIteration<value; }
+
+		protected:
+			vcg::Point3i	m_Center;						/*!< The cell whose neighboring cells are to be looked up.	*/
+			vcg::Point3i	m_CurrentNeighbor;	/*!< The neighboring cell at the current iteration.					*/
+			int						m_CurrentIteration; /*!< The current iteration.																	*/
+			int						m_TableSize;				/*!< The number of cell in the UniformGrid for each side		*/
+		}; // end of class NeighboringEntryIterator
+
+
+
+		/************************************************************************/
+		/*! \class UniformGrid
+		 * This class represent the domain U in the original article. It is used 
+		 * only during the construction of the offset and hash tables. 
+		 */ 
+		/************************************************************************/
+		class UniformGrid
+		{
+		public:
+			typedef vcg::Point3i CellCoordinate;
+
+			/*! \struct EntryIterator
+			 * This class provides a convenient way to iterate over the set of the grid cells.
+			 */
+			struct EntryIterator
+			{
+				friend class UniformGrid;
+
+				/*!
+				 * Default constructor
+				 */
+				EntryIterator(UniformGrid *uniform_grid, const CellCoordinate &position)
+				{
+					m_UniformGrid			= uniform_grid;
+					m_CurrentPosition = position;
+				}
+
+				
+				/*!
+				 * Increment operator. Move the iterator to the next cell in the UniformGrid
+				 */
+				void operator++(int)
+				{
+					if (++m_CurrentPosition.Z()==m_UniformGrid->GetResolution())
+					{
+						m_CurrentPosition.Z() = 0;
+						if (++m_CurrentPosition.Y()==m_UniformGrid->GetResolution())
+						{
+							m_CurrentPosition.Y() = 0;
+							if (++m_CurrentPosition.X()==m_UniformGrid->GetResolution())
+								m_CurrentPosition = CellCoordinate(-1, -1, -1);
+						}
+					}
+				}
+
+
+				/*!
+				 * Copy operator.
+				 * @param[in] it The iterator whose value has to be copied.
+				 */
+				void operator =(const EntryIterator &it)
+				{
+					m_UniformGrid			= it.m_UniformGrid;
+					m_CurrentPosition = it.m_CurrentPosition;
+				}
+
+				/*!
+				 * Equivalence operator
+				 */
+				bool operator==(const EntryIterator &it) const
+				{
+					return m_CurrentPosition==it.m_CurrentPosition;
+				}
+
+				/*!
+				 * Diversity operator
+				 */
+				bool operator!=(const EntryIterator &it) const
+				{
+					return m_CurrentPosition!=it.m_CurrentPosition;
+				}
+
+				/*!
+				 * Dereferencing operator.
+				 * \return The pointer to the vector of the objects contained in the cell pointed to by the iterator.
+				 */
+				std::vector< ObjectPointer >* operator*()
+				{
+					return m_UniformGrid->GetObjects(m_CurrentPosition);
+				}
+
+				/*!
+				 * Return the index of the cell pointed to by the iterator.
+				 */
+				CellCoordinate GetPosition() const
+				{
+					return m_CurrentPosition;
+				}
+
+
+			protected:
+				UniformGrid			* m_UniformGrid;
+				CellCoordinate		m_CurrentPosition;
+			}; // end of struct EntryIterator
+
+
+			/*!
+			 * Default constructor
+			 */
+			UniformGrid() {}
+
+			/*!
+			 * Default destructor
+			 */
+			~UniformGrid() {}
+
+
+			/*!
+			 * These functions return an iterator pointing to the first and the last cell of the grid respectively.
+			 */
+			EntryIterator Begin() { return EntryIterator(this, CellCoordinate( 0,  0,  0)); }
+			EntryIterator End()		{ return EntryIterator(this, CellCoordinate(-1, -1, -1)); }
+
+			
+			/*!
+			 * Return an iterator that iterates over the six adjacent cells of a given cell.
+			 * @param[in] at	The cell around which this iterator takes values.
+			 * \return				The iterator over the neighboring cells of <CODE>at</CODE>.
+			 */
+			NeighboringEntryIterator GetNeighboringEntryIterator(const CellCoordinate &at) { return NeighboringEntryIterator(at, m_CellPerSide); }
+
+
+			/*!
+			 * Allocate the necessary space for the uniform grid.
+			 * @param[in] bounding_box	The bounding box enclosing the whole dataset.
+			 * @param[in] cell_per_side	The resolution of the grid.
+			 */
+			void Allocate(const BoundingBoxType &bounding_box, const int cell_per_side)
+			{
+				m_CellPerSide = cell_per_side;
+				m_BoundingBox = bounding_box;
+				m_CellSize		= (m_BoundingBox.max - m_BoundingBox.min)/ScalarType(cell_per_side);
+				
+				m_Grid.resize(m_CellPerSide);
+				for (int i=0; i<m_CellPerSide; i++)
+				{
+					m_Grid[i].resize(m_CellPerSide);
+					for (int j=0; j<m_CellPerSide; j++)
+						m_Grid[i][j].resize(m_CellPerSide);
+				}
+			}
+
+
+			/*!
+			 * Removes all the reference to the domain data from the UniformGrid cells.
+			 */
+			void Finalize()
+			{
+				m_Grid.clear();
+			}
+
+
+			/*!
+			 * Adds a set of elements to the uniform grid.
+			 * @param[in] begin The iterator addressing the position of the first element in the range to be added.
+			 * @param[in] end		The iterator addressing the position one past the final element in the range to be added.
+			 */
+			template < class OBJECT_ITERATOR >
+			void InsertElements(const OBJECT_ITERATOR &begin, const OBJECT_ITERATOR &end)
+			{
+				typedef OBJECT_ITERATOR ObjectIterator;
+				typedef Dereferencer< typename ReferenceType< typename OBJECT_ITERATOR::value_type >::Type > ObjectDereferencer;
+
+				std::vector< CellCoordinate > cells_occupied;
+				for (ObjectIterator iObject=begin; iObject!=end; iObject++)
+				{
+					ObjectPointer pObject = &ObjectDereferencer::Reference( *iObject );
+					GetCellsIndex( pObject, cells_occupied);
+					for (std::vector< CellCoordinate >::iterator iCell=cells_occupied.begin(), eCell=cells_occupied.end(); iCell!=eCell; iCell++)
+						GetObjects( *iCell )->push_back( pObject );
+					cells_occupied.clear();
+				}
+			}
+
+
+			/*!
+			 * Given a point contained in the UniformGrid, returns the index of the cell where it's contained.
+			 * @param[in] query The 3D point.
+			 * \return					The index of the UniformGrid entry where this point is contained.
+			 */
+			inline CellCoordinate Interize(const CoordinateType &query) const
+			{
+				CellCoordinate result;
+				result.X() = (int) floorf( (query.X()-m_BoundingBox.min.X())/m_CellSize.X() );
+				result.Y() = (int) floorf( (query.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
+				result.Z() = (int) floorf( (query.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
+				return result;
+			}
+
+			/*!
+			 * Given a bounding box contained in the UniformGrid, returns its integer-equivalent bounding box.
+			 * @param[in] bounding_box	The bounding box in the 3D space.
+			 * \return									The integer representation of the bounding box.
+			 */
+			inline vcg::Box3i Interize(const BoundingBoxType &bounding_box) const
+			{
+				vcg::Box3i result;
+				result.min = Interize(bounding_box.min);
+				result.max = Interize(bounding_box.max);
+				return result;
+			}
+
+			
+			/*!
+			 * Given the pointer to an object, returns the set of cells in the uniform grid containing the object.
+			 * @param[in]		pObject					The pointer to the object
+			 * @param[out]	cells_occuppied	The set of cell index containing the object
+			 */
+			void GetCellsIndex(const ObjectPointer pObject, std::vector< CellCoordinate > & cells_occupied)
+			{
+				BoundingBoxType object_bb;
+				(*pObject).GetBBox(object_bb);
+				CoordinateType corner = object_bb.min;
+				
+				while (object_bb.IsIn(corner))
+				{
+					CellCoordinate cell_index;
+					cell_index.X() = (int) floorf( (corner.X()-m_BoundingBox.min.X())/m_CellSize.X() );
+					cell_index.Y() = (int) floorf( (corner.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
+					cell_index.Z() = (int) floorf( (corner.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
+					cells_occupied.push_back( cell_index );
+
+					if ((corner.X()+=m_CellSize.X())>object_bb.max.X())
+					{
+						corner.X() = object_bb.min.X();
+						if ( (corner.Z()+=m_CellSize.Z())>object_bb.max.Z() )
+						{
+							corner.Z() = object_bb.min.Z();
+							corner.Y() += m_CellSize.Y();
+						}
+					}
+				}
+			}
+
+
+			/*!
+			 * Return the number of cells of the uniform grid where there are no item of the input dataset.
+			 * \return The number of cells occupied by at least one item.
+			 */
+			int GetNumberOfNotEmptyCells()
+			{
+				int number_of_not_empty_cell = 0;
+				for (int i=0; i<m_CellPerSide; i++)
+					for (int j=0; j<m_CellPerSide; j++)
+						for (int k=0; k<m_CellPerSide; k++)
+							if (GetObjects(i, j, k)->size()>0)
+								number_of_not_empty_cell++;
+				return number_of_not_empty_cell;
+			}
+
+			/*!
+			 * Returns the number of entries for each side of the grid.
+			 * \return The resolution of the UniformGrid in each dimension.
+			 */
+			inline int GetResolution() const { return m_CellPerSide; }
+
+
+			/*!
+			 * Return the pointer to a vector containing pointers to the objects falling in a given domain cell.
+			 * @param[in] at	The index of the cell of the uniform grid where looking for
+			 * \return				A pointer to a vector of pointers to the objects falling in the cell having index <CODE>at</CODE>.
+			 */
+			std::vector< ObjectPointer >* GetObjects(const int i, const int j, const int k) { return &m_Grid[i][j][k]; }
+			std::vector< ObjectPointer >* GetObjects(const CellCoordinate &at)							{ return &m_Grid[at.X()][at.Y()][at.Z()];}
+			std::vector< ObjectPointer >* operator[](const CellCoordinate &at)							{ return &m_Grid[at.X()][at.Y()][at.Z()];}
+			
+		protected:
+			std::vector< std::vector< std::vector< std::vector< ObjectPointer > > > > 
+											m_Grid;					/*!< The uniform grid												*/
+			BoundingBoxType	m_BoundingBox;	/*!< The bounding box of the uniform grid.	*/
+			int							m_CellPerSide;	/*!< The number of cell per side.						*/
+			CoordinateType 	m_CellSize;			/*!< The dimension of each cell.						*/
+		}; //end of class UniformGrid
+
+
+
+
+		/************************************************************************/
+		/*! \class HashTable
+		 * This class substitutes the uniform grid. 
+		 */
+		/************************************************************************/
+		class HashTable
+		{
+		public:
+			typedef vcg::Point3i EntryCoordinate;
+			
+			// We preferred using the Data structure instead of a pointer 
+			// to the vector of the domain elements just for extensibility
+			struct Data
+			{
+				/*!
+				 * Default constructor
+				 */
+				Data(std::vector< ObjectPointer > *data) 
+				{ 
+					domain_data = data; 
+				}
+				
+				std::vector< ObjectPointer >	*domain_data;
+			};
+
+			/*!
+			 * Default constructor
+			 */
+			HashTable() {}
+
+			/*!
+			 * Default destructor
+			 */
+			~HashTable() { Clear(true); }
+
+
+			/*!
+			 * 
+			 */
+			NeighboringEntryIterator GetNeighborintEntryIterator(const EntryCoordinate &at) { return NeighboringEntryIterator(at, m_EntryPerSide); }
+
+			
+			/*!
+			 * Allocates the space for the hash table; the number of entries created is entry_per_side^3.
+			 * @param[in] entry_per_side The number of entries for each size
+			 */
+			void Allocate(const int entry_per_side)
+			{
+				m_EntryPerSide = entry_per_side;
+				m_Table.resize(m_EntryPerSide);
+				for (int i=0; i<m_EntryPerSide; i++)
+				{
+					m_Table[i].resize(m_EntryPerSide);
+					for (int j=0; j<m_EntryPerSide; j++)
+						m_Table[i][j].resize(m_EntryPerSide, NULL);
+				}
+				
+				BuildFreeEntryList();	
+			}
+
+			
+			/*
+			 * Once the PerfectSpatialHash has been computed, all the unnecessary data can be eliminated.
+			 * This function frees the empyt_list, and substitutes all the pointers to the UniformGrid 
+			 * whit brand new pointers to the input objects.
+			 */
+			void Finalize()
+			{
+				Data *pData;
+				for (int i=0; i<m_EntryPerSide; i++)
+					for (int j=0; j<m_EntryPerSide; j++)
+						for (int k=0; k<m_EntryPerSide; k++)
+							if ((pData=GetData(i, j, k))!=NULL)
+							{
+								std::vector< ObjectPointer >	*domain_data = pData->domain_data;
+								pData->domain_data = new std::vector< ObjectPointer>( *domain_data );
+							}
+
+
+				m_FreeEntries.clear();
+			}
+
+
+			/*!
+			 * Inserts each entry in the hash table in the free entry list.
+			 * When this function is called, each entry in the hash table must be free.
+			 */
+			void BuildFreeEntryList()
+			{
+				m_FreeEntries.clear();
+				for (int i=0; i<m_EntryPerSide; i++)
+					for (int j=0; j<m_EntryPerSide; j++)
+						for (int k=0; k<m_EntryPerSide; k++)
+						{
+							assert(m_Table[i][j][k]==NULL);
+							m_FreeEntries.push_back(EntryCoordinate(i, j, k));
+						}
+			}
+
+			/*!
+			 * Removes all the entries from the table and clears the free entry list
+			 */
+			void Clear(bool delete_vectors=false)
+			{
+				for (int i=0; i<m_EntryPerSide; i++)
+					for (int j=0; j<m_EntryPerSide; j++)
+						for (int k=0; k<m_EntryPerSide; k++)
+							if (m_Table[i][j][k]!=NULL)
+							{
+								if (delete_vectors)
+									delete m_Table[i][j][k]->domain_data;
+
+								delete m_Table[i][j][k];
+								m_Table[i][j][k] = NULL;
+							}
+				
+				m_FreeEntries.clear();
+			}
+
+			/*!
+			 * Returns the reference to the free entry list
+			 * \return The reference to the free entry list
+			 */
+			std::list< EntryCoordinate >* GetFreeEntryList() { return &m_FreeEntries; }
+
+			/*!
+			 * Maps a given domain entry index into a hash table index.
+			 * It corresponds to the \f$f_0\f$ function in the original article.
+			 */
+			EntryCoordinate DomainToHashTable(const typename UniformGrid::CellCoordinate &p)
+			{
+				EntryCoordinate result;
+				result.X() = p.X()%m_EntryPerSide;
+				result.Y() = p.Y()%m_EntryPerSide;
+				result.Z() = p.Z()%m_EntryPerSide;
+				return result;
+			}
+
+			/*!
+			 * Inserts a new element in the hash table at the given position.
+			 * @param[in] at		The position in the hash table where the new element will be created
+			 * @param[in] data	The set of the domain elements contained in this entry
+			 */ 
+			void SetEntry(const EntryCoordinate &at, std::vector< ObjectPointer > *data)
+			{
+				assert(IsFree(at));
+				m_Table[at.X()][at.Y()][at.Z()] = new Data(data);
+				m_FreeEntries.remove(at);
+			}
+
+			/*!
+			 * Given a hash table entry, this function modifies its coordinates in order to guarantee that
+			 * they are in the valid range. Call this function before accessing the hash table.
+			 * @param[in, out] entry The entry whose coordinates have to be checked.
+			 */ 
+			void ValidateEntry(EntryCoordinate &entry)
+			{
+				while (entry.X()<0) entry.X()+=m_EntryPerSide;
+				while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
+				while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
+			}
+
+			/*!
+			 * Check if a given position in the hash table is free.
+			 * @param[in] at The position of the hash table to check.
+			 * \return			 True if and only if the hash table is free at the given position.
+			 */
+			inline bool IsFree(const EntryCoordinate &at) const
+			{
+				return (GetData(at)==NULL);
+			}
+
+			/*!
+			 */
+			inline int GetSize() { return m_EntryPerSide; }
+
+			/*!
+			 * Returns the number of free entries.
+			 */
+			inline int GetNumberOfFreeEntries() 
+			{
+				return int(m_FreeEntries.size());
+			}
+
+			/*!
+			 * Return the number of entries where there is some domain data.
+			 */
+			inline int GetNumberOfNotEmptyEntries()
+			{
+				return (int(powf(float(m_EntryPerSide), 3.0f))-int(m_FreeEntries.size()));
+			}
+
+			/*!
+			 * Return the pointer to the data stored in the hash table at the given position.
+			 * @param[in] at	The position of the hash table where looks for the data.
+			 * \return				A pointer to a valid data only if a valid pointer is stored in the hash 
+			 *								table at the given position; otherwise return NULL.
+			 */
+			inline Data* GetData	 (const int i, const int j, const int k) const { return m_Table[i][j][k]; }
+			inline Data* GetData	 (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
+			inline Data* operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
+
+		protected:
+			int																									m_EntryPerSide; /*!< The number of entries for each side of the hash-table. */
+			std::vector< std::vector< std::vector < Data* > > > m_Table;				/*!< The table. */
+			std::list< EntryCoordinate >												m_FreeEntries;  /*!< The list containing the free entries. */
+		}; //end of class HashTable
+
+		/************************************************************************/
+		/*! \class OffsetTable
+		* This class containts the offsets used for shifting the access to the hash table.
+		*/
+		/************************************************************************/
+		class OffsetTable
+		{
+		public:
+			typedef unsigned char						OffsetType;
+			typedef vcg::Point3<OffsetType>	Offset;
+			typedef Offset								* OffsetPointer;
+			typedef vcg::Point3i						EntryCoordinate;
+			
+			/*! \struct PreImage
+			 * This class represents the pre-image for a given entry in the offset table, that is the set
+			 * \f$h_1^{-1}(q)={p \in S s.t. h_1(p)=q}\f$
+			 */
+			struct PreImage
+			{
+				/*!
+				 * Default constructor.
+				 * @param[in] at				The entry in the offset table where the cells in the preimage are mapped into.
+				 * @param[in] preimage	The set of UniformGrid cells mapping to this entry.
+				 */
+				PreImage(EntryCoordinate &at, std::vector< typename UniformGrid::CellCoordinate > *preimage)
+				{
+					entry_index = at;
+					pre_image		= preimage;
+					cardinality = int(pre_image->size());
+				}
+
+				/*!
+				 * less-than operator: needed for sorting the preimage slots based on their cardinality.
+				 * @param second 
+				 * \return <code>true</code> if and only if the cardinality of this preimage slot is greater than that of <code>second</code>.
+				 */
+				inline bool operator<(const PreImage &second) const { return (cardinality>second.cardinality); }
+				
+
+				std::vector< typename UniformGrid::CellCoordinate >	
+													*	pre_image;		/*!< The set of entries in the uniform grid whose image through \f$h_1\f$ is this entry.*/
+				EntryCoordinate			entry_index;  /*!< The index of the entry inside the offset table.	*/
+				int									cardinality;	/*!< The cardinality of the pre-image.								*/
+			}; // end of struct PreImage
+
+
+			/*!
+			 * Default constructor
+			 */
+			OffsetTable() { m_EntryPerSide=-1; m_NumberOfOccupiedEntries=0;}
+
+			/*!
+			 * Destructor
+			 */
+			~OffsetTable() { Clear(); }
+
+			/*!
+			 * Clear the entries in the offset table and in the preimage table.
+			 */
+			void Clear()
+			{
+				for (int i=0; i<m_EntryPerSide; i++)
+					for (int j=0; j<m_EntryPerSide; j++)
+						for (int k=0; k<m_EntryPerSide; k++)
+							if (m_Table[i][j][k]!=NULL)
+							{
+								delete m_Table[i][j][k];
+								m_Table[i][j][k] = NULL;
+							}
+				m_EntryPerSide = -1;
+				m_H1PreImage.clear();
+				m_NumberOfOccupiedEntries = 0;
+			}
+
+			/*!
+			 * Allocate the space necessary for a offset table containing size entries for 
+			 * each dimension and the necessary space for computing the relative anti-image. 
+			 * @param[in] size The number of entries per side to allocate. 
+			 */
+			void Allocate(int size)
+			{
+				m_NumberOfOccupiedEntries = 0;
+				
+				m_EntryPerSide = size;
+				m_Table.resize(m_EntryPerSide);
+				for (int i=0; i<m_EntryPerSide; i++)
+				{
+					m_Table[i].resize(m_EntryPerSide);
+					for (int j=0; j<m_EntryPerSide; j++)
+						m_Table[i][j].resize(m_EntryPerSide, NULL);
+				}
+
+				m_H1PreImage.resize(m_EntryPerSide);
+				for (int i=0; i<m_EntryPerSide; i++)
+				{
+					m_H1PreImage[i].resize(m_EntryPerSide);
+					for (int j=0; j<m_EntryPerSide; j++)
+						m_H1PreImage[i][j].resize(m_EntryPerSide);
+				}
+			}
+
+
+			/*!
+			 *
+			 */
+			void Finalize()
+			{
+				m_H1PreImage.clear();
+			}
+
+
+			/*!
+			 * Build the pre-image of the \f$h_1\f$ function: the <CODE>m_H1PreImage</CODE> grid contains, for each  
+			 * cell (i, j, k) a list of the domain grid (the UniformGrid) that are mapped through \f$h_1\f$ into that cell.
+			 */
+			void BuildH1PreImage(const typename UniformGrid::EntryIterator &begin, const typename UniformGrid::EntryIterator &end)
+			{
+				for (typename UniformGrid::EntryIterator iter=begin; iter!=end; iter++)
+				{
+					if ((*iter)->size()==0)
+						continue;
+
+					typename UniformGrid::CellCoordinate cell_index = iter.GetPosition();
+					EntryCoordinate at = DomainToOffsetTable(cell_index);
+					m_H1PreImage[at.X()][at.Y()][at.Z()].push_back(cell_index);
+				}
+			}
+
+			/*!
+			 * Sorts the entries of the PreImage table based on their cardinality.
+			 * @param[out] preimage The list containing the entries of the preimage sorted by cardinality
+			 */
+			void GetPreImageSortedPerCardinality(std::list< PreImage > &pre_image)
+			{
+				pre_image.clear();
+				for (int i=0; i<m_EntryPerSide; i++)
+					for (int j=0; j<m_EntryPerSide; j++)
+						for (int k=0; k<m_EntryPerSide; k++)
+						{
+							std::vector< typedef UniformGrid::CellCoordinate > *preimage = &m_H1PreImage[i][j][k];
+							if (preimage->size()>0)
+								pre_image.push_back( PreImage(typename UniformGrid::CellCoordinate(i, j, k), preimage) );
+						}
+				pre_image.sort();
+			}
+			
+				
+			/*!
+			 * Check if the entries in the offset table near the given entry contain a valid offset.
+			 * @param[in]				at The entry of the offset table whose neighboring entries will be checked.
+			 * @param[out] offsets The set of consistent offset found by inspecting the neighboring entries.
+			 * \return a vector containing possible offsets for the given entry
+			 */
+			void SuggestConsistentOffsets(const EntryCoordinate &at, std::vector< Offset > &offsets)
+			{
+				offsets.clear();
+				for (int i=-1; i<2; i++)
+					for (int j=-1; j<2; j++)
+						for (int k=-1; k<2; k++)
+						{
+							if (i==0 && j==0 && k==0)
+								continue;
+
+							int x = (at.X()+i+m_EntryPerSide)%m_EntryPerSide;
+							int y = (at.Y()+j+m_EntryPerSide)%m_EntryPerSide;
+							int z = (at.Z()+k+m_EntryPerSide)%m_EntryPerSide;
+							EntryCoordinate neighboring_entry(x, y, z);
+							if (!IsFree(neighboring_entry))
+								offsets.push_back( *GetOffset(neighboring_entry) );
+						}
+			}
+
+			
+			/*!
+			 * Assures that the given entry can be used to access the offset table without throwing an out-of-bound exception.
+			 * @param[in,out] entry The entry to be checked.
+			 */
+			void ValidateEntryCoordinate(EntryCoordinate &entry)
+			{
+				while (entry.X()<0) entry.X()+=m_EntryPerSide;
+				while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
+				while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
+			}
+
+			/*!
+			 * Converts the coordinate of a given cell in the UniformGrid to a valid entry in the offset table.
+			 * This function corresponds to the \f$h_1\f$ function of the article.
+			 * @param[in] coord The index of a domain cell in the UniformGrid.
+			 * \return					The coordinate of the entry corresponding to <CODE>coord<CODE> through this mapping.
+			 */
+			EntryCoordinate DomainToOffsetTable(const typename UniformGrid::CellCoordinate &coord)
+			{
+				EntryCoordinate result;
+				result.X() = coord.X()%m_EntryPerSide;
+				result.Y() = coord.Y()%m_EntryPerSide;
+				result.Z() = coord.Z()%m_EntryPerSide;
+				return result;
+			}
+
+			/*!
+			 * Adds a new element to the offset table.
+			 * @param[in] coord		The index of the UniformGrid cell whose offset has to be stored.
+			 * @param[in] offset	The offset to associate to the given UniformGrid cell.
+			 */
+			void SetOffset(const typename UniformGrid::CellCoordinate &coord, const Offset &offset)
+			{
+				EntryCoordinate entry = DomainToOffsetTable(coord);
+				assert(IsFree(entry));
+				m_Table[entry.X()][entry.Y()][entry.Z()] = new Offset(offset);
+				m_NumberOfOccupiedEntries++;
+			}
+			
+			/*!
+			 * Return a random offset: this function is used during the first steps of the creation process,
+			 * when the offsets are computed at random.
+			 * @param[out] A random offset
+			 */
+			void GetRandomOffset( Offset &offset )
+			{
+				offset.X() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
+				offset.Y() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
+				offset.Z() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
+			}
+
+			
+			/*!
+			* Return the number of entries of the offset table for each dimension.
+			* \return The number of entries for each side
+			*/
+			inline int	GetSize()			const {return m_EntryPerSide;}
+
+
+			/*!
+			* Checks if the given entry in the offset table is free
+			* @param[in] at The coordinate of the entry to be checked.
+			* \return			 true if and only if the entry with coordinate <CODE>at</CODE> is free.
+			*/
+			inline bool IsFree(const EntryCoordinate &at) const { return GetOffset(at)==NULL;	} 
+			//{ return m_Table[at.X()][at.Y()][at.Z()]==NULL;	}
+
+
+			/*!
+			 * Return the number of entries containing a valid offset.
+			 * \return The number of not empty entries.
+			 */
+			inline int GetNumberOfOccupiedCells() const { return m_NumberOfOccupiedEntries;		}
+
+			/*!
+			 * Return the pointer to the offset stored at the given entry. NULL if that entry doesn't contain a offset
+			 */
+			inline OffsetPointer& GetOffset (const int i, const int j, const int k)				{ return m_Table[i][j][k]; } 
+			inline OffsetPointer  GetOffset (const int i, const int j, const int k) const { return m_Table[i][j][k]; } 
+			
+			inline OffsetPointer& GetOffset (const EntryCoordinate &at)				{ return m_Table[at.X()][at.Y()][at.Z()]; } 
+			inline OffsetPointer  GetOffset (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; } 
+			
+			inline OffsetPointer& operator[](const EntryCoordinate &at)				{ return m_Table[at.X()][at.Y()][at.Z()]; } 
+			inline OffsetPointer  operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; } 
+
+		protected:
+			const static int	m_MAX_VERSOR_LENGTH = 256;	/*!< The maximal length of the single component of each offset. */ 
+			int								m_EntryPerSide;							/*!< The resolution of the offset table.												*/
+			int								m_NumberOfOccupiedEntries;	/*!< The number of entries containing a valid offset.						*/
+			std::vector< std::vector< std::vector< OffsetPointer > > >																				m_Table;			/*!< The offset table.					*/
+			std::vector< std::vector< std::vector< std::vector< typename UniformGrid::CellCoordinate > > > >	m_H1PreImage; /*!< The \f$f1\f$ pre-image.		*/
+		}; //end of class OffsetTable
+
+
+		
+		/*******************************************************************************************************************************/
+		/*! \class BinaryImage
+		* This class is used to encode the sparsity of the dataset. Since the hash table stores data associated with a sparse 
+		* subset of the domain, it may be necessary to determine if an arbitrary query point lies in this defined domain.
+		*/
+		/*******************************************************************************************************************************/
+		 class BinaryImage
+		 {
+		 public:
+			 /*!
+			  * Default constructor
+				*/
+			 BinaryImage() 
+			 {
+				 m_Resolution = -1;
+			 }
+
+			 
+			 /*!
+			  * Destructor
+				*/
+			 ~BinaryImage() {}
+
+
+			 /*!
+			  * Allocate the space necessary to encode the distribution of the dataset over the domain.
+				* @param[in] size The resolution on each dimension of the bitmap.
+				*/
+			 void Allocate(const int size)
+			 {
+				 m_Resolution = size;
+				 m_Mask.resize(m_Resolution);
+				 for (int i=0; i<m_Resolution; i++)
+				 {
+					 m_Mask[i].resize(m_Resolution);
+					 for (int j=0; j<m_Resolution; j++)
+						 m_Mask[i][j].resize(m_Resolution, false);
+				 }
+			 }
+
+
+			 /*!
+			  * Removes all flags from the bitmap. This function is called after a failed attempt to create the offset-table.
+				*/
+			 void Clear()
+			 {
+				 for (int i=0; i<m_Resolution; i++)
+					 for (int j=0; j<m_Resolution; j++)
+						 std::fill(m_Mask[i][j].begin(), m_Mask[i][j].end(), false);
+			 }
+
+
+			 /*!
+			  * Checks if a portion of the dataset fall inside the UniformGrid cell having coordinate <CODE>at</CODE>.
+				* @param[in] at The index of the UniformGrid cell to check.
+				* \return				<code>true</code> if and only if a portion of the dataset in included in this  UniformGrid cell.
+				*/
+			 inline bool ContainsData(const typename UniformGrid::CellCoordinate &at) const { return GetFlag(at)==true;}
+
+			 
+			 /*!
+			  * Returns the number of entries in each dimension.
+				* \return The resolution of the BinaryImage.
+				*/
+			 inline int GetResolution() const { return m_Resolution; }
+
+			 
+			 /*!
+			  * Return the value stored in the d-dimensional bitmap at the given position.
+				* @param[in] i 
+				* @param[in] j 
+				* @param[in] k
+				*	\return
+				*/
+			 inline bool operator()(const int i, const int j, const int k)					{ return m_Mask[i][j][k]; }
+			 
+			 
+			 /*!
+			  * Return the value stored at the given position in the d-dimensional bitmap.
+				* @param[in] at
+				* \return
+				*/
+			 inline bool operator[](const typename UniformGrid::CellCoordinate &at)	{ return m_Mask[at.X()][at.Y()][at.Z()]; }
+			 inline bool& GetFlag(const int i, const int j, const int k)const				{ return m_Mask[i][j][k]; }
+			 inline void  SetFlat(const int i, const int j, const int k)						{ m_Mask[i][j][k] = true; }
+			 
+
+			 inline bool  GetFlag(const typename UniformGrid::CellCoordinate &at)	const			{ return m_Mask[at.X()][at.Y()][at.Z()]; }
+			 inline void  SetFlag(const typename UniformGrid::CellCoordinate &at)						{ m_Mask[at.X()][at.Y()][at.Z()] = true; }
+
+		 protected:
+			 std::vector< std::vector< std::vector< bool > > > 
+						m_Mask;					/*!< The bitmap image.							*/
+			 int	m_Resolution;		/*!< The resolution of the bitmap.	*/
+		 }; // end of class BinaryImage
+
+
+
+		 /*******************************************************************************************************************************/
+		 /*! \struct Neighbor
+		 * This class is used to retrieve the neighboring objects in the spatial queries and to sort them.
+		 */
+		 /*******************************************************************************************************************************/
+		 struct Neighbor
+		 {
+			 /*!
+			  * Default constructor
+				*/
+			 Neighbor() 
+			 { 
+				 object		= NULL;
+				 distance = ScalarType(-1.0);
+				 nearest_point.Zero();
+			 }
+
+			 
+			 /*!
+			  * Constructor
+				* @param[in] pObject	The pointer to the object.
+				* @param[in] dist			The distance between the object and the query point.
+				* @param[in] point		The point on the object having minimal distance from the query point.
+				*/
+			 Neighbor(ObjectPointer pObject, ScalarType dist, CoordinateType point)
+			 {
+				 object = pObject;
+				 distance = dist;
+				 nearest_point(point);
+			 }
+
+			 
+			 /*!
+			  * Less than operator. Needed for sorting a range of neighbor based on their distance from the query object.
+				*/
+			 inline bool operator<(const Neighbor &second)
+			 {
+				 return distance<second.distance;
+			 }
+
+			 ObjectPointer	object;
+			 ScalarType			distance;
+			 CoordinateType nearest_point;
+		 }; // end of struct Neighbor
+
+		///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// 
+		/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+		/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+		/************************************************************************/
+		/* Functions	                                                          */
+		/************************************************************************/
+	public:
+		/*!
+		 * The hash table can be constructed following two different approaches:
+		 * the first is more fast, but might allocate a offset table bigger than the necessary
+		 * the second try to construct the offset table up to m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION times, and then chooses the minimum size for which the construction succeeded. 
+		 */
+		enum		ConstructionApproach { FastConstructionApproach=0, CompactConstructionApproach=1 };
+
+		/*!
+		 * Default constructor
+		 */
+		PerfectSpatialHashing() { srand( (unsigned) time(NULL) ); }
+
+		/*!
+		 * Destructor
+		 */
+		~PerfectSpatialHashing() { /* ... I don't remember if there is something to delete! :D */ } 
+
+		template < class OBJECT_ITERATOR >
+		void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj) 
+		{ Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, NULL); }
+
+		template < class OBJECT_ITERATOR >
+		void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, vcg::CallBackPos *callback) 
+		{ Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, callback); }
+
+		template < class OBJECT_ITERATOR >
+		void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach) 
+		{ Set<OBJECT_ITERATOR>(bObj, eObj, approach, NULL); }
+
+		/*!
+		 * Add the elements to the PerfectSpatialHashing data structure. Since this structure can handle only
+		 * static dataset, the elements mustn't be changed while using this structure.
+		 * @param[in] bObj			The iterator addressing the first element to be included in the hashing.
+		 * @param[in] eObj			The iterator addressing the position after the last element to be included in the hashing.
+		 * @param[in] approach	Either <code>FastConstructionApproach</code> or <code>CompactConstructionApproach</code>.
+		 * @param[in] callback  The callback to call to provide information about the progress of the computation.
+		 */
+		template < class OBJECT_ITERATOR >
+		void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach, vcg::CallBackPos *callback)
+		{
+			BoundingBoxType bounding_box;
+			BoundingBoxType object_bb;
+			bounding_box.SetNull();
+			for (OBJECT_ITERATOR iObj=bObj; iObj!=eObj; iObj++)
+			{
+				(*iObj).GetBBox(object_bb);
+				bounding_box.Add(object_bb);
+			}	
+
+			//...and expand it a bit more
+			BoundingBoxType	resulting_bb(bounding_box);
+			CoordinateType	offset = bounding_box.Dim()*float(m_BOUNDING_BOX_EXPANSION_FACTOR);
+			CoordinateType	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);
+
+			int number_of_objects = int(std::distance(bObj, eObj));
+
+			// Try to find a reasonable space partition
+#ifdef _USE_GRID_UTIL_PARTIONING_
+			vcg::Point3i resolution;
+			vcg::BestDim<ScalarType>(number_of_objects, resulting_bb.Dim(), resolution); 
+			int cells_per_side = resolution.X();
+#else ifdef _USE_OCTREE_PARTITIONING_ // Alternative to find the resolution of the uniform grid:
+			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(number_of_objects)/density);
+				depth++;
+			}
+			while (primitives_per_voxel>16 && depth<15);
+			int cells_per_side = int(powf(2.0f, float(depth)));
+#endif
+
+			m_UniformGrid.Allocate(resulting_bb, cells_per_side);
+			m_UniformGrid.InsertElements(bObj, eObj);
+			m_Bitmap.Allocate(cells_per_side);
+			int number_of_cells_occupied = m_UniformGrid.GetNumberOfNotEmptyCells();
+
+			int hash_table_size = (int) ceilf(powf(float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
+			if (hash_table_size>256)
+				hash_table_size = (int) ceilf(powf(1.01f*float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
+			m_HashTable.Allocate(hash_table_size);
+
+			switch (approach)
+			{
+			case FastConstructionApproach		:	PerformFastConstruction(number_of_cells_occupied, callback)		; break;
+			case CompactConstructionApproach: PerformCompactConstruction(number_of_cells_occupied, callback); break;
+			default: assert(false);
+			}
+			Finalize();
+		} // end of method Set
+
+
+		/*!
+		 * Returns all the objects contained inside a specified sphere
+		 * @param[in]  distance_functor
+		 * @param[in]	 marker
+		 * @param[in]	 sphere_center
+		 * @param[in]	 sphere_radius
+		 * @param[out] objects
+		 * @param[out] distances
+		 * @param[out] points
+		 * \return
+		 */
+		template <class OBJECT_POINT_DISTANCE_FUNCTOR, class OBJECT_MARKER, class OBJECT_POINTER_CONTAINER, class DISTANCE_CONTAINER, class POINT_CONTAINER>
+		unsigned int GetInSphere
+			(
+			OBJECT_POINT_DISTANCE_FUNCTOR		&	distance_functor, 
+			OBJECT_MARKER										&	marker,
+			const CoordType									&	sphere_center, 
+			const ScalarType								&	sphere_radius,
+			OBJECT_POINTER_CONTAINER				&	objects, 
+			DISTANCE_CONTAINER							&	distances, 
+			POINT_CONTAINER									&	points,
+			bool															sort_per_distance   = true,
+			bool															allow_zero_distance = true
+			)
+		{
+			BoundingBoxType query_bb(sphere_center, sphere_radius);
+			vcg::Box3i integer_bb = m_UniformGrid.Interize(query_bb);
+
+			vcg::Point3i index;
+			std::vector< std::vector< ObjectPointer >* > contained_objects;
+			std::vector< ObjectPointer >* tmp;
+			for (index.X()=integer_bb.min.X(); index.X()<=integer_bb.max.X(); index.X()++)
+				for (index.Y()=integer_bb.min.Y(); index.Y()<=integer_bb.max.Y(); index.Y()++)
+					for (index.Z()=integer_bb.min.Z(); index.Z()<=integer_bb.max.Z(); index.Z()++)
+						if ((tmp=(*this)[index])!=NULL)
+							contained_objects.push_back(tmp);
+
+			std::vector< Neighbor > results;
+			for (std::vector< std::vector< ObjectPointer >* >::iterator iVec=contained_objects.begin(), eVec=contained_objects.end(); iVec!=eVec; iVec++)
+				for (std::vector< ObjectPointer >::iterator iObj=(*iVec)->begin(), eObj=(*iVec)->end(); iObj!=eObj; iObj++ )
+				{
+					int r = int(results.size());
+					results.push_back(Neighbor());
+					results[r].object		= *iObj;
+					results[r].distance = sphere_radius;
+					if (!distance_functor(*results[r].object, sphere_center, results[r].distance, results[r].nearest_point) || (results[r].distance==ScalarType(0.0) && !allow_zero_distance) )
+						results.pop_back();
+				}
+
+			if (sort_per_distance)
+				std::sort( results.begin(), results.end() );
+
+			int number_of_objects = int(results.size());
+			distances.resize(number_of_objects);
+			points.resize(number_of_objects);
+			objects.resize(number_of_objects);
+			for (int i=0, size=int(results.size()); i<size; i++)
+			{
+				distances[i]	= results[i].distance;
+				points[i]			= results[i].nearest_point;
+				objects[i]		= results[i].object;
+			}
+			return number_of_objects;
+		} //end of GetInSphere
+
+		
+		/*!
+		* Once the offset table has been built, this function can be used to access the data.
+		* Given a 3D point in the space, this function returns the set of ObjectPointers contained
+		* in the same UniformGrid cell where this point is contained.
+		* @param[in] query The 3D query point.
+		* \return		The pointer to the vector of ObjectPointers contained in the same UG of query. NULL if any.
+		*/		
+		std::vector< ObjectPointer >* operator[](const CoordinateType &query)
+		{
+			typename UniformGrid::CellCoordinate ug_index = m_UniformGrid.Interize(query);
+			if (!m_Bitmap[ug_index])
+				return NULL;
+
+			HashTable::EntryCoordinate ht_index = PerfectHashFunction(ug_index);
+			std::vector< ObjectPointer >* result = m_HashTable[ht_index];
+			return result;
+		}
+
+
+		std::vector< ObjectPointer >* operator[](const typename UniformGrid::CellCoordinate  &query)
+		{
+			if(!m_Bitmap[query])
+				return NULL;
+
+			HashTable::EntryCoordinate ht_index = PerfectHashFunction(query);
+			std::vector< ObjectPointer >* result = m_HashTable[ht_index]->domain_data;
+			return result;
+		}
+
+
+	protected:
+		/*!
+		* The injective mapping from the set of occupied cells to a slot in the hash-table
+		* @param[in]	query		The index of a domain cell whose mapping has to be calculated.
+		* @param[out]	result	The index of a entry in the hash-table where query is mapped into.
+		*/
+		typename HashTable::EntryCoordinate PerfectHashFunction(const typename UniformGrid::CellCoordinate &query)
+		{
+			typename HashTable::EntryCoordinate result;
+			OffsetTable::OffsetPointer offset = m_OffsetTable[ m_OffsetTable.DomainToOffsetTable(query) ];
+			result = m_HashTable.DomainToHashTable( Shift(query, *offset) );
+			return result;
+		} 
+
+		
+		/*!
+		 * Performs the addition between a entry coordinate and an offset.
+		 * @param[in] entry		The index of a given cell.
+		 * @param[in] offset	The offset that must be applied to the entry.
+		 * \return						The entry resulting by the addition of entry and offset.
+		 */
+		typename HashTable::EntryCoordinate Shift(const vcg::Point3i &entry, const typename OffsetTable::Offset &offset)
+		{
+			HashTable::EntryCoordinate result;
+			result.X() = entry.X() + int(offset.X());
+			result.Y() = entry.Y() + int(offset.Y());
+			result.Z() = entry.Z() + int(offset.Z());
+			return result;
+		}
+
+
+		/*!
+		* Finalizes the data structures at the end of the offset-table construction.
+		* This function eliminates all unnecessary data, and encodes sparsity.
+		* TODO At the moment, the sparsity encoding is implemented thought a bitmap, i.e. a boolean grid 
+		*			where each slot tells if the relative UniformGrid has a valid entry in the HashTable.
+		*/
+		void Finalize()
+		{
+#ifdef _DEBUG
+			for (UniformGrid::EntryIterator iUGEntry=m_UniformGrid.Begin(), eUGEntry=m_UniformGrid.End(); iUGEntry!=eUGEntry; iUGEntry++)
+				assert(m_Bitmap.ContainsData(iUGEntry.GetPosition())==((*iUGEntry)->size()>0));
+#endif
+			m_HashTable.Finalize();
+			m_UniformGrid.Finalize();
+			m_OffsetTable.Finalize();
+		}
+
+		
+		/*!
+		 * Check if the given offset is valid for a set of domain cell.
+		 * @param[in] pre_image
+		 * @param[in] offset
+		 * \return
+		 */
+		bool IsAValidOffset(const std::vector< typename UniformGrid::CellCoordinate > *pre_image, const typename OffsetTable::Offset &offset)
+		{
+			int ht_size			= m_HashTable.GetSize();
+			int sqr_ht_size = ht_size*ht_size;
+			std::vector< int > involved_entries;
+			for (int i=0, pre_image_size=int((*pre_image).size()); i<pre_image_size; i++)
+			{
+				typename UniformGrid::CellCoordinate domain_entry = (*pre_image)[i];
+				typename HashTable::EntryCoordinate  hash_entry		= m_HashTable.DomainToHashTable( Shift(domain_entry, offset) );
+				if (!m_HashTable.IsFree(hash_entry))
+					return false;
+				else
+					involved_entries.push_back(hash_entry.X()*sqr_ht_size + hash_entry.Y()*ht_size + hash_entry.Z());
+			}
+
+			// In order to guarantee that the PerfectHashFunction is injective, the image of each domain-entry must be unique in the hash-table
+			std::sort(involved_entries.begin(), involved_entries.end());
+			for (int i=0, j=1, size=int(involved_entries.size()); j<size; i++, j++)
+				if (involved_entries[i]==involved_entries[j])
+					return false;
+
+			return true;
+		}
+
+		
+		/*!
+		 * Given the size of the hash table and an initial seed for the size of the offset table, returns an appropriate size
+		 * for the offset table in order to avoid less effective constructions.
+		 * @param[in] hash_table_size		The number of entries for each side of the hash-table.
+		 * @param[in] offset_table_size The number of entries for each side of the offset-table.
+		 * \return											The next appropriate size for the offset table.
+		 */
+		int GetUnefectiveOffsetTableSize(const int hash_table_size, const int offset_table_size)
+		{
+			int result = offset_table_size;
+			while (GreatestCommonDivisor(hash_table_size, result)!=1 || hash_table_size%result==0)
+				result += (hash_table_size%2==0) ? (result%2)+1 : 1;  //Change the offset table size, otherwise the hash construction should be ineffective
+			return result;
+		}
+
+		
+		/*!
+		* Start the construction of the offset table trying to complete as quickly as possible.
+		* Sometimes the dimension of the offset table will not be minimal.
+		* @param[in] number_of_filled_cells The number of entries in the uniform grid containing some elements of the dataset.
+		*/
+		void PerformFastConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
+		{
+			int offset_table_size = (int) ceilf(powf(m_SIGMA*float(number_of_filled_cells), 1.0f/float(m_DIMENSION)));
+			int hash_table_size   = m_HashTable.GetSize();
+			int failed_construction_count = 0;
+			do 
+			{
+				offset_table_size += failed_construction_count++;
+				offset_table_size  = GetUnefectiveOffsetTableSize(hash_table_size, offset_table_size);
+			} 
+			while(!OffsetTableConstructionSucceded(offset_table_size, callback));
+		}
+
+			
+		/*!
+		* Start the construction of the offset table trying to minimize its dimension.
+		* The offset table size is chosen by performing a binary search over possible values for the offset table size.
+		* For this reason, this method will generally require more time than PerformFastConstruction.
+		*/
+		void PerformCompactConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
+		{
+			int min_successfully_dimension	= std::numeric_limits<int>::max();
+			int hash_table_size							= m_HashTable.GetSize();
+			int half_hash_table_size				= int(float(hash_table_size)/2.0f);
+
+			// According to the original article, a maximum number of trials are to be made in order to select the optimal (i.e. minimal) offset table size
+			for (int t=0; t<m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION; t++)
+			{
+				int lower_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(double(rand())/double(RAND_MAX)*half_hash_table_size + 1) );
+				int upper_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(((double) rand() / (double) RAND_MAX) * hash_table_size + half_hash_table_size));
+
+				// The construction of the offset table using this pair of values surely succeed for max, but not for min
+				// The optimal value for the offset table size is then found via binary search over this pair of values
+				int candidate_offset_table_size;
+				int last_tried_size = -1;
+				while (lower_bound<upper_bound)
+				{
+					candidate_offset_table_size = GetUnefectiveOffsetTableSize(hash_table_size, int(floorf((lower_bound+upper_bound)/2.0f)));
+
+					// If in the previous iteration the offset table has been successfully constructed with the same size,
+					// a minimum for the range [lower_bound, upper_bound] has been found
+					if (last_tried_size==candidate_offset_table_size)
+						break; 
+
+					if ( OffsetTableConstructionSucceded((last_tried_size=candidate_offset_table_size), callback) )
+					{
+						upper_bound = candidate_offset_table_size;
+						min_successfully_dimension = vcg::math::Min(candidate_offset_table_size, min_successfully_dimension);
+					}
+					else
+						lower_bound = candidate_offset_table_size;
+
+					m_HashTable.Clear();
+					m_HashTable.BuildFreeEntryList();
+					m_OffsetTable.Clear();
+					m_Bitmap.Clear();
+				}
+#ifdef _DEBUD
+				printf("\nPerfectSpatialHashing: minimum offset table size found at the %d-th iteration was %d\n", (t+1), min_successfully_dimension);
+#endif
+			}
+
+			// Finally the OffsetTable must be constructed using the minimum valid size
+			while (!OffsetTableConstructionSucceded(min_successfully_dimension, callback)) 
+			{
+				m_HashTable.Clear();
+				m_HashTable.BuildFreeEntryList();
+				m_OffsetTable.Clear();
+				m_Bitmap.Clear();
+			}
+		}
+
+
+
+		/*!
+		* Try to construct the offset table for a given size
+		*	\param[in] offset_table_size  The size of the offset table.
+		*	\return												<CODE>true</CODE> if and only if the construction of the offset table succeeds.
+		*/
+		bool OffsetTableConstructionSucceded(const int offset_table_size, vcg::CallBackPos *callback)
+		{
+			m_OffsetTable.Allocate(offset_table_size); // Create the Offset table
+			m_OffsetTable.BuildH1PreImage(m_UniformGrid.Begin(), m_UniformGrid.End()); // Build the f0 pre-image
+
+			std::list< OffsetTable::PreImage > preimage_slots;
+			m_OffsetTable.GetPreImageSortedPerCardinality(preimage_slots);
+
+			char msg[128];
+			sprintf(msg, "Building offset table of resolution %d", m_OffsetTable.GetSize());
+			int	step = int(preimage_slots.size())/100;
+			int number_of_slots = int(preimage_slots.size());
+			int perc = 0; 
+			int iter = 0;
+			for (std::list< OffsetTable::PreImage >::iterator iPreImage=preimage_slots.begin(), ePreImage=preimage_slots.end(); iPreImage!=ePreImage; iPreImage++, iter++)
+			{
+				if (callback!=NULL && iter%step==0 && (perc=iter*100/number_of_slots)<100) (*callback)(perc, msg);
+
+				bool found_valid_offset = false;
+				typename OffsetTable::Offset candidate_offset;
+
+				// Heuristic #1: try to set the offset value to one stored in a neighboring entry of the offset table
+				std::vector< typename OffsetTable::Offset > consistent_offsets;
+				m_OffsetTable.SuggestConsistentOffsets( (*iPreImage).entry_index, consistent_offsets);
+				for (std::vector< typename OffsetTable::Offset >::iterator iOffset=consistent_offsets.begin(), eOffset=consistent_offsets.end(); iOffset!=eOffset && !found_valid_offset; iOffset++)
+					if (IsAValidOffset(iPreImage->pre_image, *iOffset))
+					{
+						found_valid_offset = true;
+						candidate_offset = *iOffset;
+					}
+
+
+					// Heuristic #2: 
+					if (!found_valid_offset)
+					{
+						std::vector< typename UniformGrid::CellCoordinate > *pre_image = (*iPreImage).pre_image;
+						for (std::vector< typename UniformGrid::CellCoordinate >::iterator iPreImage=pre_image->begin(), ePreImage=pre_image->end(); iPreImage!=ePreImage && !found_valid_offset; iPreImage++)
+							for (NeighboringEntryIterator iUGNeighbourhood=m_UniformGrid.GetNeighboringEntryIterator(*iPreImage); iUGNeighbourhood<6 && !found_valid_offset; iUGNeighbourhood++ )
+								if (!m_OffsetTable.IsFree( m_OffsetTable.DomainToOffsetTable( *iUGNeighbourhood ) ))
+								{
+									HashTable::EntryCoordinate ht_entry = PerfectHashFunction(*iUGNeighbourhood);
+									for (NeighboringEntryIterator iHTNeighbourhood=m_HashTable.GetNeighborintEntryIterator(ht_entry); iHTNeighbourhood<6 && !found_valid_offset; iHTNeighbourhood++)
+										if (m_HashTable.IsFree(*iHTNeighbourhood))
+										{
+											candidate_offset.Import( *iHTNeighbourhood-m_HashTable.DomainToHashTable(*iPreImage) ) ;
+											// m_OffsetTable.ValidateEntry( candidate_offset ); Is'n necessary, becouse the offset type is unsigned char.
+											if (IsAValidOffset(pre_image, candidate_offset))
+												found_valid_offset = true;
+										}
+								}
+					}
+
+					if (!found_valid_offset)
+					{
+						// At the beginning, the offset can be found via random searches
+						for (int i=0; i<m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET && !found_valid_offset; i++)
+						{
+							HashTable::EntryCoordinate base_entry = (*iPreImage).pre_image->at(0);
+							do 
+								m_OffsetTable.GetRandomOffset(candidate_offset);
+							while (!m_HashTable.IsFree( m_HashTable.DomainToHashTable( Shift(base_entry, candidate_offset) ) ));
+
+							if (IsAValidOffset( (*iPreImage).pre_image, candidate_offset))
+								found_valid_offset = true;
+						}
+
+						// The chance to find a valid offset table via random searches decreases toward the end of the offset table construction:
+						// So a exhaustive search over all the free hash table entries is performed.
+						for (std::list< HashTable::EntryCoordinate >::const_iterator iFreeCell=m_HashTable.GetFreeEntryList()->begin(), eFreeCell=m_HashTable.GetFreeEntryList()->end(); iFreeCell!=eFreeCell && !found_valid_offset; iFreeCell++)
+						{
+							UniformGrid::CellCoordinate domain_entry		= (*iPreImage).pre_image->at(0);
+							OffsetTable::EntryCoordinate offset_entry		= m_OffsetTable.DomainToOffsetTable(domain_entry);
+							HashTable::EntryCoordinate hashtable_entry	= m_HashTable.DomainToHashTable(domain_entry);
+							candidate_offset.Import(*iFreeCell - hashtable_entry);
+
+							if ( IsAValidOffset(iPreImage->pre_image, candidate_offset) )
+								found_valid_offset = true;
+						}
+				}
+
+					// If a valid offset has been found, the construction of the offset table continues, 
+					// otherwise the offset table must be enlarged and the construction repeated
+					if (found_valid_offset)
+					{
+						m_OffsetTable.SetOffset( (*iPreImage->pre_image).at(0), candidate_offset);
+						for (int c=0, pre_image_cardinality = iPreImage->cardinality; c<pre_image_cardinality; c++)
+						{
+							HashTable::EntryCoordinate ht_entry = PerfectHashFunction( (*iPreImage->pre_image).at(c));
+							std::vector< ObjectPointer > *domain_data = m_UniformGrid[ (*iPreImage->pre_image).at(c) ];
+							m_HashTable.SetEntry(ht_entry, domain_data /*, (*iPreImage->pre_image).at(c)*/); // might be useful for encoding sparsity
+							m_Bitmap.SetFlag((*iPreImage->pre_image).at(c));
+						}
+					}
+					else
+					{
+						m_OffsetTable.Clear();
+						m_HashTable.Clear();
+						m_HashTable.BuildFreeEntryList();
+						m_Bitmap.Clear();
+						return false; 
+					}
+			}
+
+			if (callback!=NULL) (*callback)(100, msg);
+			return true;
+		} // end of OffsetTableConstructionSucceded
+
+
+		/************************************************************************/
+		/* Data Members                                                         */
+		/************************************************************************/
+	protected: 
+		UniformGrid m_UniformGrid;	/*!< The uniform grid used for partitioning the volume.						*/
+		OffsetTable m_OffsetTable;	/*!< The offset table corresponding to \f$\Phi\f$ in the article. */
+		HashTable		m_HashTable;		/*!< The hash table that will substitute the uniform grid.				*/
+		BinaryImage	m_Bitmap;
+
+		const static float	m_BOUNDING_BOX_EXPANSION_FACTOR;
+		const static float	m_SIGMA;
+		const static int		m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION;
+		const static int		m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET;
+		const static int    m_DIMENSION;
+	}; //end of class PerfectSpatialHashing
+
+	/*! @} */ 
+	//end of Doxygen documentation
+}//end of namespace vcg
+
+template < class OBJECT_TYPE, class SCALAR_TYPE >
+const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET	= 32;
+
+template < class OBJECT_TYPE, class SCALAR_TYPE >
+const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION	= 5;
+
+template < class OBJECT_TYPE, class SCALAR_TYPE >
+const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_DIMENSION	= 3;
+
+template < class OBJECT_TYPE, class SCALAR_TYPE >
+const SCALAR_TYPE vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_BOUNDING_BOX_EXPANSION_FACTOR	= SCALAR_TYPE(0.035);
+
+template < class OBJECT_TYPE, class SCALAR_TYPE >
+const SCALAR_TYPE vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_SIGMA	= SCALAR_TYPE(1.0f/(2.0f*SCALAR_TYPE(m_DIMENSION)));
+
+#endif //VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
\ No newline at end of file