1 /****************************************************************************
2 * VCGLib                                                            o o     *
3 * Visual and Computer Graphics Library                            o     o   *
4 *                                                                _   O  _   *
5 * Copyright(C) 2004-2016                                           \/)\/    *
6 * Visual Computing Lab                                            /\/|      *
7 * ISTI - Italian National Research Council                           |      *
8 *                                                                    \      *
9 * All rights reserved.                                                      *
10 *                                                                           *
11 * This program is free software; you can redistribute it and/or modify      *
12 * it under the terms of the GNU General Public License as published by      *
13 * the Free Software Foundation; either version 2 of the License, or         *
14 * (at your option) any later version.                                       *
15 *                                                                           *
16 * This program is distributed in the hope that it will be useful,           *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
20 * for more details.                                                         *
21 *                                                                           *
22 ****************************************************************************/
23 
24 #ifndef VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
25 #define VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
26 
27 #pragma warning(disable : 4996)
28 
29 #define _USE_GRID_UTIL_PARTIONING_ 1
30 #define _USE_OCTREE_PARTITIONING_  (1-_USE_GRID_UTIL_PARTIONING_)
31 
32 #include <vector>
33 #include <list>
34 #include <algorithm>
35 
36 #include <vcg/space/index/base.h>
37 #include <vcg/space/index/grid_util.h>
38 
39 #include <vcg/space/point2.h>
40 #include <vcg/space/point3.h>
41 #include <vcg/space/box3.h>
42 
43 namespace vcg
44 {
45 
46     // Compute the greatest common divisor between two integers a and b
GreatestCommonDivisor(const int a,const int b)47     int GreatestCommonDivisor(const int a, const int b)
48     {
49         int m = a;
50         int n = b;
51 
52         do
53         {
54             if (m<n) std::swap(m, n);
55             m = m % n;
56             std::swap(m, n);
57         }
58         while (n!=0);
59         return m;
60     }
61 
62     // Doxygen documentation
63     /** \addtogroup index */
64     /*! @{ */
65 
66     /*!
67      * This class implements the perfect spatial hashing by S.Lefebvre and H.Hoppe
68      * This is an spatial indexing structure such as the uniform grid, but with lower
69      * memory requirements, since all the empty cells of the uniform grid are removed.
70      * Access to a non-empty cell is performed looking up in two d-dimensional tables,
71      * the offset table and the hash table.
72      * @param OBJECT_TYPE (Template parameter) the type of objects to be indexed
73      * @param SCALAR_TYPE (Template parameter) the scalar type
74      */
75     template < class OBJECT_TYPE, class SCALAR_TYPE >
76     class PerfectSpatialHashing : public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE >
77     {
78         // Given an object or a pointer to an object, return the reference to the object
79         template < typename TYPE >
80         struct Dereferencer
81         {
ReferenceDereferencer82             static				TYPE& Reference(TYPE				 &t)	{	return  t;	}
ReferenceDereferencer83             static				TYPE& Reference(TYPE*			 &t)	{ return *t;  }
ReferenceDereferencer84             static const  TYPE& Reference(const TYPE	 &t)	{ return  t;	}
ReferenceDereferencer85             static const	TYPE& Reference(const TYPE* &t)	{ return *t;	}
86         };
87 
88         // Given a type, holds this type in Type
89         template < typename TYPE >
90         struct ReferenceType { typedef TYPE Type; };
91 
92         // Given as type a "pointer to type", holds the type in Type
93         template < typename TYPE >
94         struct ReferenceType< TYPE * > { typedef typename ReferenceType<TYPE>::Type Type; };
95 
96     public:
97         typedef						SCALAR_TYPE													ScalarType;
98         typedef						OBJECT_TYPE													ObjectType;
99         typedef typename	ReferenceType< ObjectType >::Type * ObjectPointer;
100         typedef typename  vcg::Box3< ScalarType >							BoundingBoxType;
101         typedef typename	vcg::Point3< ScalarType >						CoordinateType;
102 
103     protected:
104         /*! \struct NeighboringEntryIterator
105         * This class provides a convenient way to iterate over the six neighboring cells of a given cell.
106         */
107         struct NeighboringEntryIterator
108         {
109             /*!
110             * Default constructor.
111             * @param[in] entry				The index of the cell in the UniformGrid around which iterate.
112             * @param[in] table_size  The number of cells in the UniformGrid for each side.
113             */
114             NeighboringEntryIterator(const vcg::Point3i &entry, const int table_size)
115             {
116                 m_Center		= entry;
117                 m_TableSize = table_size;
118                 m_CurrentNeighbor.X() = (m_Center.X()+m_TableSize-1)%m_TableSize;
119                 m_CurrentNeighbor.Y() = m_Center.Y();
120                 m_CurrentNeighbor.Z() = m_Center.Z();
121                 m_CurrentIteration		= 0;
122             }
123 
124             /*!
125             * Increment the iterator to point to the next neighboring entry in the UniformGrid
126             */
127             void operator++(int)
128             {
129                 switch(++m_CurrentIteration)
130                 {
131                 case 1: m_CurrentNeighbor.X()=(m_Center.X()+1)%m_TableSize; break;
132                 case 2: m_CurrentNeighbor.X()=m_Center.X(); m_CurrentNeighbor.Y()=(m_Center.Y()+m_TableSize-1)%m_TableSize; break;
133                 case 3: m_CurrentNeighbor.Y()=(m_Center.Y()+1)%m_TableSize; break;
134                 case 4: m_CurrentNeighbor.Y()=m_Center.Y(); m_CurrentNeighbor.Z()=(m_Center.Z()+m_TableSize-1)%m_TableSize; break;
135                 case 5: m_CurrentNeighbor.Z()=(m_Center.Z()+1)%m_TableSize; break;
136                 default: m_CurrentNeighbor = vcg::Point3i(-1, -1, -1); break;
137                 }
138             }
139 
140             /*!
141             * Dereferencing operator
142             * \return The neighbor of the given cell at the current iteration
143             */
144             vcg::Point3i operator*() { return m_CurrentNeighbor; }
145 
146             /*!
147             * Assignment operator
148             * @param[in] The source neighboring iterator
149             * \return		The reference to this iterator
150             */
151             NeighboringEntryIterator& operator =(const NeighboringEntryIterator &it)
152             {
153                 m_Center						= it.m_Center						;
154                 m_CurrentNeighbor		= it.m_CurrentNeighbor	;
155                 m_CurrentIteration	= it.m_CurrentIteration	;
156                 m_TableSize					= it.m_TableSize				;
157                 return *this;
158             }
159 
160             /*!
161             * Less than operator. Since each entry in the UniformGrid has only 6 neighbors,
162             * the iterator over the neighboring entries can be compared with an integer.
163             */
164             inline bool operator <(const int value) { return m_CurrentIteration<value; }
165 
166         protected:
167             vcg::Point3i	m_Center;						/*!< The cell whose neighboring cells are to be looked up.	*/
168             vcg::Point3i	m_CurrentNeighbor;	/*!< The neighboring cell at the current iteration.					*/
169             int						m_CurrentIteration; /*!< The current iteration.																	*/
170             int						m_TableSize;				/*!< The number of cell in the UniformGrid for each side		*/
171         }; // end of class NeighboringEntryIterator
172 
173 
174 
175         /************************************************************************/
176         /*! \class UniformGrid
177          * This class represent the domain U in the original article. It is used
178          * only during the construction of the offset and hash tables.
179          */
180         /************************************************************************/
181         class UniformGrid
182         {
183         public:
184             typedef vcg::Point3i CellCoordinate;
185 
186             /*! \struct EntryIterator
187              * This class provides a convenient way to iterate over the set of the grid cells.
188              */
189             struct EntryIterator
190             {
191                 friend class UniformGrid;
192 
193                 /*!
194                  * Default constructor
195                  */
196                 EntryIterator(UniformGrid *uniform_grid, const CellCoordinate &position)
197                 {
198                     m_UniformGrid			= uniform_grid;
199                     m_CurrentPosition = position;
200                 }
201 
202 
203                 /*!
204                  * Increment operator. Move the iterator to the next cell in the UniformGrid
205                  */
206                 void operator++(int)
207                 {
208                     if (++m_CurrentPosition.Z()==m_UniformGrid->GetResolution())
209                     {
210                         m_CurrentPosition.Z() = 0;
211                         if (++m_CurrentPosition.Y()==m_UniformGrid->GetResolution())
212                         {
213                             m_CurrentPosition.Y() = 0;
214                             if (++m_CurrentPosition.X()==m_UniformGrid->GetResolution())
215                                 m_CurrentPosition = CellCoordinate(-1, -1, -1);
216                         }
217                     }
218                 }
219 
220 
221                 /*!
222                  * Copy operator.
223                  * @param[in] it The iterator whose value has to be copied.
224                  */
225                 void operator =(const EntryIterator &it)
226                 {
227                     m_UniformGrid			= it.m_UniformGrid;
228                     m_CurrentPosition = it.m_CurrentPosition;
229                 }
230 
231                 /*!
232                  * Equivalence operator
233                  */
234                 bool operator==(const EntryIterator &it) const
235                 {
236                     return m_CurrentPosition==it.m_CurrentPosition;
237                 }
238 
239                 /*!
240                  * Diversity operator
241                  */
242                 bool operator!=(const EntryIterator &it) const
243                 {
244                     return m_CurrentPosition!=it.m_CurrentPosition;
245                 }
246 
247                 /*!
248                  * Dereferencing operator.
249                  * \return The pointer to the vector of the objects contained in the cell pointed to by the iterator.
250                  */
251                 std::vector< ObjectPointer >* operator*()
252                 {
253                     return m_UniformGrid->GetObjects(m_CurrentPosition);
254                 }
255 
256                 /*!
257                  * Return the index of the cell pointed to by the iterator.
258                  */
259                 CellCoordinate GetPosition() const
260                 {
261                     return m_CurrentPosition;
262                 }
263 
264 
265             protected:
266                 UniformGrid			* m_UniformGrid;
267                 CellCoordinate		m_CurrentPosition;
268             }; // end of struct EntryIterator
269 
270 
271             /*!
272              * Default constructor
273              */
274             UniformGrid() {}
275 
276             /*!
277              * Default destructor
278              */
279             ~UniformGrid() {}
280 
281 
282             /*!
283              * These functions return an iterator pointing to the first and the last cell of the grid respectively.
284              */
285             EntryIterator Begin() { return EntryIterator(this, CellCoordinate( 0,  0,  0)); }
286             EntryIterator End()		{ return EntryIterator(this, CellCoordinate(-1, -1, -1)); }
287 
288 
289             /*!
290              * Return an iterator that iterates over the six adjacent cells of a given cell.
291              * @param[in] at	The cell around which this iterator takes values.
292              * \return				The iterator over the neighboring cells of <CODE>at</CODE>.
293              */
294             NeighboringEntryIterator GetNeighboringEntryIterator(const CellCoordinate &at) { return NeighboringEntryIterator(at, m_CellPerSide); }
295 
296 
297             /*!
298              * Allocate the necessary space for the uniform grid.
299              * @param[in] bounding_box	The bounding box enclosing the whole dataset.
300              * @param[in] cell_per_side	The resolution of the grid.
301              */
302             void Allocate(const BoundingBoxType &bounding_box, const int cell_per_side)
303             {
304                 m_CellPerSide = cell_per_side;
305                 m_BoundingBox = bounding_box;
306                 m_CellSize		= (m_BoundingBox.max - m_BoundingBox.min)/ScalarType(cell_per_side);
307 
308                 m_Grid.resize(m_CellPerSide);
309                 for (int i=0; i<m_CellPerSide; i++)
310                 {
311                     m_Grid[i].resize(m_CellPerSide);
312                     for (int j=0; j<m_CellPerSide; j++)
313                         m_Grid[i][j].resize(m_CellPerSide);
314                 }
315             }
316 
317 
318             /*!
319              * Removes all the reference to the domain data from the UniformGrid cells.
320              */
321             void Finalize()
322             {
323                 m_Grid.clear();
324             }
325 
326 
327             /*!
328              * Adds a set of elements to the uniform grid.
329              * @param[in] begin The iterator addressing the position of the first element in the range to be added.
330              * @param[in] end		The iterator addressing the position one past the final element in the range to be added.
331              */
332             template < class OBJECT_ITERATOR >
333             void InsertElements(const OBJECT_ITERATOR &begin, const OBJECT_ITERATOR &end)
334             {
335                 typedef OBJECT_ITERATOR ObjectIterator;
336                 typedef Dereferencer< typename ReferenceType< typename OBJECT_ITERATOR::value_type >::Type > ObjectDereferencer;
337 
338                 std::vector< CellCoordinate > cells_occupied;
339                 for (ObjectIterator iObject=begin; iObject!=end; iObject++)
340                 {
341                     ObjectPointer pObject = &ObjectDereferencer::Reference( *iObject );
342                     GetCellsIndex( pObject, cells_occupied);
343                     for (std::vector< CellCoordinate >::iterator iCell=cells_occupied.begin(), eCell=cells_occupied.end(); iCell!=eCell; iCell++)
344                         GetObjects( *iCell )->push_back( pObject );
345                     cells_occupied.clear();
346                 }
347             }
348 
349 
350             /*!
351              * Given a point contained in the UniformGrid, returns the index of the cell where it's contained.
352              * @param[in] query The 3D point.
353              * \return					The index of the UniformGrid entry where this point is contained.
354              */
355             inline CellCoordinate Interize(const CoordinateType &query) const
356             {
357                 CellCoordinate result;
358                 result.X() = (int) floorf( (query.X()-m_BoundingBox.min.X())/m_CellSize.X() );
359                 result.Y() = (int) floorf( (query.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
360                 result.Z() = (int) floorf( (query.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
361                 return result;
362             }
363 
364             /*!
365              * Given a bounding box contained in the UniformGrid, returns its integer-equivalent bounding box.
366              * @param[in] bounding_box	The bounding box in the 3D space.
367              * \return									The integer representation of the bounding box.
368              */
369             inline vcg::Box3i Interize(const BoundingBoxType &bounding_box) const
370             {
371                 vcg::Box3i result;
372                 result.min = Interize(bounding_box.min);
373                 result.max = Interize(bounding_box.max);
374                 return result;
375             }
376 
377 
378             /*!
379              * Given the pointer to an object, returns the set of cells in the uniform grid containing the object.
380              * @param[in]		pObject					The pointer to the object
381              * @param[out]	cells_occuppied	The set of cell index containing the object
382              */
383             void GetCellsIndex(const ObjectPointer pObject, std::vector< CellCoordinate > & cells_occupied)
384             {
385                 BoundingBoxType object_bb;
386                 (*pObject).GetBBox(object_bb);
387                 CoordinateType corner = object_bb.min;
388 
389                 while (object_bb.IsIn(corner))
390                 {
391                     CellCoordinate cell_index;
392                     cell_index.X() = (int) floorf( (corner.X()-m_BoundingBox.min.X())/m_CellSize.X() );
393                     cell_index.Y() = (int) floorf( (corner.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
394                     cell_index.Z() = (int) floorf( (corner.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
395                     cells_occupied.push_back( cell_index );
396 
397                     if ((corner.X()+=m_CellSize.X())>object_bb.max.X())
398                     {
399                         corner.X() = object_bb.min.X();
400                         if ( (corner.Z()+=m_CellSize.Z())>object_bb.max.Z() )
401                         {
402                             corner.Z() = object_bb.min.Z();
403                             corner.Y() += m_CellSize.Y();
404                         }
405                     }
406                 }
407             }
408 
409 
410             /*!
411              * Return the number of cells of the uniform grid where there are no item of the input dataset.
412              * \return The number of cells occupied by at least one item.
413              */
414             int GetNumberOfNotEmptyCells()
415             {
416                 int number_of_not_empty_cell = 0;
417                 for (int i=0; i<m_CellPerSide; i++)
418                     for (int j=0; j<m_CellPerSide; j++)
419                         for (int k=0; k<m_CellPerSide; k++)
420                             if (GetObjects(i, j, k)->size()>0)
421                                 number_of_not_empty_cell++;
422                 return number_of_not_empty_cell;
423             }
424 
425             /*!
426              * Returns the number of entries for each side of the grid.
427              * \return The resolution of the UniformGrid in each dimension.
428              */
429             inline int GetResolution() const { return m_CellPerSide; }
430 
431 
432             /*!
433              * Return the pointer to a vector containing pointers to the objects falling in a given domain cell.
434              * @param[in] at	The index of the cell of the uniform grid where looking for
435              * \return				A pointer to a vector of pointers to the objects falling in the cell having index <CODE>at</CODE>.
436              */
437             std::vector< ObjectPointer >* GetObjects(const int i, const int j, const int k) { return &m_Grid[i][j][k]; }
438             std::vector< ObjectPointer >* GetObjects(const CellCoordinate &at)							{ return &m_Grid[at.X()][at.Y()][at.Z()];}
439             std::vector< ObjectPointer >* operator[](const CellCoordinate &at)							{ return &m_Grid[at.X()][at.Y()][at.Z()];}
440 
441         protected:
442             std::vector< std::vector< std::vector< std::vector< ObjectPointer > > > >
443                                             m_Grid;					/*!< The uniform grid												*/
444             BoundingBoxType	m_BoundingBox;	/*!< The bounding box of the uniform grid.	*/
445             int							m_CellPerSide;	/*!< The number of cell per side.						*/
446             CoordinateType 	m_CellSize;			/*!< The dimension of each cell.						*/
447         }; //end of class UniformGrid
448 
449 
450 
451 
452         /************************************************************************/
453         /*! \class HashTable
454          * This class substitutes the uniform grid.
455          */
456         /************************************************************************/
457         class HashTable
458         {
459         public:
460             typedef vcg::Point3i EntryCoordinate;
461 
462             // We preferred using the Data structure instead of a pointer
463             // to the vector of the domain elements just for extensibility
464             struct Data
465             {
466                 /*!
467                  * Default constructor
468                  */
469                 Data(std::vector< ObjectPointer > *data)
470                 {
471                     domain_data = data;
472                 }
473 
474                 std::vector< ObjectPointer >	*domain_data;
475             };
476 
477             /*!
478              * Default constructor
479              */
480             HashTable() {}
481 
482             /*!
483              * Default destructor
484              */
485             ~HashTable() { Clear(true); }
486 
487 
488             /*!
489              *
490              */
491             NeighboringEntryIterator GetNeighborintEntryIterator(const EntryCoordinate &at) { return NeighboringEntryIterator(at, m_EntryPerSide); }
492 
493 
494             /*!
495              * Allocates the space for the hash table; the number of entries created is entry_per_side^3.
496              * @param[in] entry_per_side The number of entries for each size
497              */
498             void Allocate(const int entry_per_side)
499             {
500                 m_EntryPerSide = entry_per_side;
501                 m_Table.resize(m_EntryPerSide);
502                 for (int i=0; i<m_EntryPerSide; i++)
503                 {
504                     m_Table[i].resize(m_EntryPerSide);
505                     for (int j=0; j<m_EntryPerSide; j++)
506                         m_Table[i][j].resize(m_EntryPerSide, NULL);
507                 }
508 
509                 BuildFreeEntryList();
510             }
511 
512 
513             /*
514              * Once the PerfectSpatialHash has been computed, all the unnecessary data can be eliminated.
515              * This function frees the empyt_list, and substitutes all the pointers to the UniformGrid
516              * whit brand new pointers to the input objects.
517              */
518             void Finalize()
519             {
520                 Data *pData;
521                 for (int i=0; i<m_EntryPerSide; i++)
522                     for (int j=0; j<m_EntryPerSide; j++)
523                         for (int k=0; k<m_EntryPerSide; k++)
524                             if ((pData=GetData(i, j, k))!=NULL)
525                             {
526                                 std::vector< ObjectPointer >	*domain_data = pData->domain_data;
527                                 pData->domain_data = new std::vector< ObjectPointer>( *domain_data );
528                             }
529 
530 
531                 m_FreeEntries.clear();
532             }
533 
534 
535             /*!
536              * Inserts each entry in the hash table in the free entry list.
537              * When this function is called, each entry in the hash table must be free.
538              */
539             void BuildFreeEntryList()
540             {
541                 m_FreeEntries.clear();
542                 for (int i=0; i<m_EntryPerSide; i++)
543                     for (int j=0; j<m_EntryPerSide; j++)
544                         for (int k=0; k<m_EntryPerSide; k++)
545                         {
546                             assert(m_Table[i][j][k]==NULL);
547                             m_FreeEntries.push_back(EntryCoordinate(i, j, k));
548                         }
549             }
550 
551             /*!
552              * Removes all the entries from the table and clears the free entry list
553              */
554             void Clear(bool delete_vectors=false)
555             {
556                 for (int i=0; i<m_EntryPerSide; i++)
557                     for (int j=0; j<m_EntryPerSide; j++)
558                         for (int k=0; k<m_EntryPerSide; k++)
559                             if (m_Table[i][j][k]!=NULL)
560                             {
561                                 if (delete_vectors)
562                                     delete m_Table[i][j][k]->domain_data;
563 
564                                 delete m_Table[i][j][k];
565                                 m_Table[i][j][k] = NULL;
566                             }
567 
568                 m_FreeEntries.clear();
569             }
570 
571             /*!
572              * Returns the reference to the free entry list
573              * \return The reference to the free entry list
574              */
575             std::list< EntryCoordinate >* GetFreeEntryList() { return &m_FreeEntries; }
576 
577             /*!
578              * Maps a given domain entry index into a hash table index.
579              * It corresponds to the \f$f_0\f$ function in the original article.
580              */
581             EntryCoordinate DomainToHashTable(const typename UniformGrid::CellCoordinate &p)
582             {
583                 EntryCoordinate result;
584                 result.X() = p.X()%m_EntryPerSide;
585                 result.Y() = p.Y()%m_EntryPerSide;
586                 result.Z() = p.Z()%m_EntryPerSide;
587                 return result;
588             }
589 
590             /*!
591              * Inserts a new element in the hash table at the given position.
592              * @param[in] at		The position in the hash table where the new element will be created
593              * @param[in] data	The set of the domain elements contained in this entry
594              */
595             void SetEntry(const EntryCoordinate &at, std::vector< ObjectPointer > *data)
596             {
597                 assert(IsFree(at));
598                 m_Table[at.X()][at.Y()][at.Z()] = new Data(data);
599                 m_FreeEntries.remove(at);
600             }
601 
602             /*!
603              * Given a hash table entry, this function modifies its coordinates in order to guarantee that
604              * they are in the valid range. Call this function before accessing the hash table.
605              * @param[in, out] entry The entry whose coordinates have to be checked.
606              */
607             void ValidateEntry(EntryCoordinate &entry)
608             {
609                 while (entry.X()<0) entry.X()+=m_EntryPerSide;
610                 while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
611                 while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
612             }
613 
614             /*!
615              * Check if a given position in the hash table is free.
616              * @param[in] at The position of the hash table to check.
617              * \return			 True if and only if the hash table is free at the given position.
618              */
619             inline bool IsFree(const EntryCoordinate &at) const
620             {
621                 return (GetData(at)==NULL);
622             }
623 
624             /*!
625              */
626             inline int GetSize() { return m_EntryPerSide; }
627 
628             /*!
629              * Returns the number of free entries.
630              */
631             inline int GetNumberOfFreeEntries()
632             {
633                 return int(m_FreeEntries.size());
634             }
635 
636             /*!
637              * Return the number of entries where there is some domain data.
638              */
639             inline int GetNumberOfNotEmptyEntries()
640             {
641                 return (int(powf(float(m_EntryPerSide), 3.0f))-int(m_FreeEntries.size()));
642             }
643 
644             /*!
645              * Return the pointer to the data stored in the hash table at the given position.
646              * @param[in] at	The position of the hash table where looks for the data.
647              * \return				A pointer to a valid data only if a valid pointer is stored in the hash
648              *								table at the given position; otherwise return NULL.
649              */
650             inline Data* GetData	 (const int i, const int j, const int k) const { return m_Table[i][j][k]; }
651             inline Data* GetData	 (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
652             inline Data* operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
653 
654         protected:
655             int																									m_EntryPerSide; /*!< The number of entries for each side of the hash-table. */
656             std::vector< std::vector< std::vector < Data* > > > m_Table;				/*!< The table. */
657             std::list< EntryCoordinate >												m_FreeEntries;  /*!< The list containing the free entries. */
658         }; //end of class HashTable
659 
660         /************************************************************************/
661         /*! \class OffsetTable
662         * This class containts the offsets used for shifting the access to the hash table.
663         */
664         /************************************************************************/
665         class OffsetTable
666         {
667         public:
668             typedef unsigned char						OffsetType;
669             typedef vcg::Point3<OffsetType>	Offset;
670             typedef Offset								* OffsetPointer;
671             typedef vcg::Point3i						EntryCoordinate;
672 
673             /*! \struct PreImage
674              * This class represents the pre-image for a given entry in the offset table, that is the set
675              * \f$h_1^{-1}(q)={p \in S s.t. h_1(p)=q}\f$
676              */
677             struct PreImage
678             {
679                 /*!
680                  * Default constructor.
681                  * @param[in] at				The entry in the offset table where the cells in the preimage are mapped into.
682                  * @param[in] preimage	The set of UniformGrid cells mapping to this entry.
683                  */
684                 PreImage(EntryCoordinate &at, std::vector< typename UniformGrid::CellCoordinate > *preimage)
685                 {
686                     entry_index = at;
687                     pre_image		= preimage;
688                     cardinality = int(pre_image->size());
689                 }
690 
691                 /*!
692                  * less-than operator: needed for sorting the preimage slots based on their cardinality.
693                  * @param second
694                  * \return <code>true</code> if and only if the cardinality of this preimage slot is greater than that of <code>second</code>.
695                  */
696                 inline bool operator<(const PreImage &second) const { return (cardinality>second.cardinality); }
697 
698 
699                 std::vector< typename UniformGrid::CellCoordinate >
700                                                     *	pre_image;		/*!< The set of entries in the uniform grid whose image through \f$h_1\f$ is this entry.*/
701                 EntryCoordinate			entry_index;  /*!< The index of the entry inside the offset table.	*/
702                 int									cardinality;	/*!< The cardinality of the pre-image.								*/
703             }; // end of struct PreImage
704 
705 
706             /*!
707              * Default constructor
708              */
709             OffsetTable() { m_EntryPerSide=-1; m_NumberOfOccupiedEntries=0;}
710 
711             /*!
712              * Destructor
713              */
714             ~OffsetTable() { Clear(); }
715 
716             /*!
717              * Clear the entries in the offset table and in the preimage table.
718              */
719             void Clear()
720             {
721                 for (int i=0; i<m_EntryPerSide; i++)
722                     for (int j=0; j<m_EntryPerSide; j++)
723                         for (int k=0; k<m_EntryPerSide; k++)
724                             if (m_Table[i][j][k]!=NULL)
725                             {
726                                 delete m_Table[i][j][k];
727                                 m_Table[i][j][k] = NULL;
728                             }
729                 m_EntryPerSide = -1;
730                 m_H1PreImage.clear();
731                 m_NumberOfOccupiedEntries = 0;
732             }
733 
734             /*!
735              * Allocate the space necessary for a offset table containing size entries for
736              * each dimension and the necessary space for computing the relative anti-image.
737              * @param[in] size The number of entries per side to allocate.
738              */
739             void Allocate(int size)
740             {
741                 m_NumberOfOccupiedEntries = 0;
742 
743                 m_EntryPerSide = size;
744                 m_Table.resize(m_EntryPerSide);
745                 for (int i=0; i<m_EntryPerSide; i++)
746                 {
747                     m_Table[i].resize(m_EntryPerSide);
748                     for (int j=0; j<m_EntryPerSide; j++)
749                         m_Table[i][j].resize(m_EntryPerSide, NULL);
750                 }
751 
752                 m_H1PreImage.resize(m_EntryPerSide);
753                 for (int i=0; i<m_EntryPerSide; i++)
754                 {
755                     m_H1PreImage[i].resize(m_EntryPerSide);
756                     for (int j=0; j<m_EntryPerSide; j++)
757                         m_H1PreImage[i][j].resize(m_EntryPerSide);
758                 }
759             }
760 
761 
762             /*!
763              *
764              */
765             void Finalize()
766             {
767                 m_H1PreImage.clear();
768             }
769 
770 
771             /*!
772              * Build the pre-image of the \f$h_1\f$ function: the <CODE>m_H1PreImage</CODE> grid contains, for each
773              * cell (i, j, k) a list of the domain grid (the UniformGrid) that are mapped through \f$h_1\f$ into that cell.
774              */
775             void BuildH1PreImage(const typename UniformGrid::EntryIterator &begin, const typename UniformGrid::EntryIterator &end)
776             {
777                 for (typename UniformGrid::EntryIterator iter=begin; iter!=end; iter++)
778                 {
779                     if ((*iter)->size()==0)
780                         continue;
781 
782                     typename UniformGrid::CellCoordinate cell_index = iter.GetPosition();
783                     EntryCoordinate at = DomainToOffsetTable(cell_index);
784                     m_H1PreImage[at.X()][at.Y()][at.Z()].push_back(cell_index);
785                 }
786             }
787 
788             /*!
789              * Sorts the entries of the PreImage table based on their cardinality.
790              * @param[out] preimage The list containing the entries of the preimage sorted by cardinality
791              */
792             void GetPreImageSortedPerCardinality(std::list< PreImage > &pre_image)
793             {
794                 pre_image.clear();
795                 for (int i=0; i<m_EntryPerSide; i++)
796                     for (int j=0; j<m_EntryPerSide; j++)
797                         for (int k=0; k<m_EntryPerSide; k++)
798                         {
799                             std::vector< typename UniformGrid::CellCoordinate > *preimage = &m_H1PreImage[i][j][k];
800                             if (preimage->size()>0)
801                                 pre_image.push_back( PreImage(typename UniformGrid::CellCoordinate(i, j, k), preimage) );
802                         }
803                 pre_image.sort();
804             }
805 
806 
807             /*!
808              * Check if the entries in the offset table near the given entry contain a valid offset.
809              * @param[in]				at The entry of the offset table whose neighboring entries will be checked.
810              * @param[out] offsets The set of consistent offset found by inspecting the neighboring entries.
811              * \return a vector containing possible offsets for the given entry
812              */
813             void SuggestConsistentOffsets(const EntryCoordinate &at, std::vector< Offset > &offsets)
814             {
815                 offsets.clear();
816                 for (int i=-1; i<2; i++)
817                     for (int j=-1; j<2; j++)
818                         for (int k=-1; k<2; k++)
819                         {
820                             if (i==0 && j==0 && k==0)
821                                 continue;
822 
823                             int x = (at.X()+i+m_EntryPerSide)%m_EntryPerSide;
824                             int y = (at.Y()+j+m_EntryPerSide)%m_EntryPerSide;
825                             int z = (at.Z()+k+m_EntryPerSide)%m_EntryPerSide;
826                             EntryCoordinate neighboring_entry(x, y, z);
827                             if (!IsFree(neighboring_entry))
828                                 offsets.push_back( *GetOffset(neighboring_entry) );
829                         }
830             }
831 
832 
833             /*!
834              * Assures that the given entry can be used to access the offset table without throwing an out-of-bound exception.
835              * @param[in,out] entry The entry to be checked.
836              */
837             void ValidateEntryCoordinate(EntryCoordinate &entry)
838             {
839                 while (entry.X()<0) entry.X()+=m_EntryPerSide;
840                 while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
841                 while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
842             }
843 
844             /*!
845              * Converts the coordinate of a given cell in the UniformGrid to a valid entry in the offset table.
846              * This function corresponds to the \f$h_1\f$ function of the article.
847              * @param[in] coord The index of a domain cell in the UniformGrid.
848              * \return					The coordinate of the entry corresponding to <CODE>coord<CODE> through this mapping.
849              */
850             EntryCoordinate DomainToOffsetTable(const typename UniformGrid::CellCoordinate &coord)
851             {
852                 EntryCoordinate result;
853                 result.X() = coord.X()%m_EntryPerSide;
854                 result.Y() = coord.Y()%m_EntryPerSide;
855                 result.Z() = coord.Z()%m_EntryPerSide;
856                 return result;
857             }
858 
859             /*!
860              * Adds a new element to the offset table.
861              * @param[in] coord		The index of the UniformGrid cell whose offset has to be stored.
862              * @param[in] offset	The offset to associate to the given UniformGrid cell.
863              */
864             void SetOffset(const typename UniformGrid::CellCoordinate &coord, const Offset &offset)
865             {
866                 EntryCoordinate entry = DomainToOffsetTable(coord);
867                 assert(IsFree(entry));
868                 m_Table[entry.X()][entry.Y()][entry.Z()] = new Offset(offset);
869                 m_NumberOfOccupiedEntries++;
870             }
871 
872             /*!
873              * Return a random offset: this function is used during the first steps of the creation process,
874              * when the offsets are computed at random.
875              * @param[out] A random offset
876              */
877             void GetRandomOffset( Offset &offset )
878             {
879                 offset.X() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
880                 offset.Y() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
881                 offset.Z() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
882             }
883 
884 
885             /*!
886             * Return the number of entries of the offset table for each dimension.
887             * \return The number of entries for each side
888             */
889             inline int	GetSize()			const {return m_EntryPerSide;}
890 
891 
892             /*!
893             * Checks if the given entry in the offset table is free
894             * @param[in] at The coordinate of the entry to be checked.
895             * \return			 true if and only if the entry with coordinate <CODE>at</CODE> is free.
896             */
897             inline bool IsFree(const EntryCoordinate &at) const { return GetOffset(at)==NULL;	}
898             //{ return m_Table[at.X()][at.Y()][at.Z()]==NULL;	}
899 
900 
901             /*!
902              * Return the number of entries containing a valid offset.
903              * \return The number of not empty entries.
904              */
905             inline int GetNumberOfOccupiedCells() const { return m_NumberOfOccupiedEntries;		}
906 
907             /*!
908              * Return the pointer to the offset stored at the given entry. NULL if that entry doesn't contain a offset
909              */
910             inline OffsetPointer& GetOffset (const int i, const int j, const int k)				{ return m_Table[i][j][k]; }
911             inline OffsetPointer  GetOffset (const int i, const int j, const int k) const { return m_Table[i][j][k]; }
912 
913             inline OffsetPointer& GetOffset (const EntryCoordinate &at)				{ return m_Table[at.X()][at.Y()][at.Z()]; }
914             inline OffsetPointer  GetOffset (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
915 
916             inline OffsetPointer& operator[](const EntryCoordinate &at)				{ return m_Table[at.X()][at.Y()][at.Z()]; }
917             inline OffsetPointer  operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
918 
919         protected:
920             const static int	m_MAX_VERSOR_LENGTH = 256;	/*!< The maximal length of the single component of each offset. */
921             int								m_EntryPerSide;							/*!< The resolution of the offset table.												*/
922             int								m_NumberOfOccupiedEntries;	/*!< The number of entries containing a valid offset.						*/
923             std::vector< std::vector< std::vector< OffsetPointer > > >																				m_Table;			/*!< The offset table.					*/
924             std::vector< std::vector< std::vector< std::vector< typename UniformGrid::CellCoordinate > > > >	m_H1PreImage; /*!< The \f$f1\f$ pre-image.		*/
925         }; //end of class OffsetTable
926 
927 
928 
929         /*******************************************************************************************************************************/
930         /*! \class BinaryImage
931         * This class is used to encode the sparsity of the dataset. Since the hash table stores data associated with a sparse
932         * subset of the domain, it may be necessary to determine if an arbitrary query point lies in this defined domain.
933         */
934         /*******************************************************************************************************************************/
935          class BinaryImage
936          {
937          public:
938              /*!
939               * Default constructor
940                 */
941              BinaryImage()
942              {
943                  m_Resolution = -1;
944              }
945 
946 
947              /*!
948               * Destructor
949                 */
950              ~BinaryImage() {}
951 
952 
953              /*!
954               * Allocate the space necessary to encode the distribution of the dataset over the domain.
955                 * @param[in] size The resolution on each dimension of the bitmap.
956                 */
957              void Allocate(const int size)
958              {
959                  m_Resolution = size;
960                  m_Mask.resize(m_Resolution);
961                  for (int i=0; i<m_Resolution; i++)
962                  {
963                      m_Mask[i].resize(m_Resolution);
964                      for (int j=0; j<m_Resolution; j++)
965                          m_Mask[i][j].resize(m_Resolution, false);
966                  }
967              }
968 
969 
970              /*!
971               * Removes all flags from the bitmap. This function is called after a failed attempt to create the offset-table.
972                 */
973              void Clear()
974              {
975                  for (int i=0; i<m_Resolution; i++)
976                      for (int j=0; j<m_Resolution; j++)
977                          std::fill(m_Mask[i][j].begin(), m_Mask[i][j].end(), false);
978              }
979 
980 
981              /*!
982               * Checks if a portion of the dataset fall inside the UniformGrid cell having coordinate <CODE>at</CODE>.
983                 * @param[in] at The index of the UniformGrid cell to check.
984                 * \return				<code>true</code> if and only if a portion of the dataset in included in this  UniformGrid cell.
985                 */
986              inline bool ContainsData(const typename UniformGrid::CellCoordinate &at) const { return GetFlag(at)==true;}
987 
988 
989              /*!
990               * Returns the number of entries in each dimension.
991                 * \return The resolution of the BinaryImage.
992                 */
993              inline int GetResolution() const { return m_Resolution; }
994 
995 
996              /*!
997               * Return the value stored in the d-dimensional bitmap at the given position.
998                 * @param[in] i
999                 * @param[in] j
1000                 * @param[in] k
1001                 *	\return
1002                 */
1003              inline bool operator()(const int i, const int j, const int k)					{ return m_Mask[i][j][k]; }
1004 
1005 
1006              /*!
1007               * Return the value stored at the given position in the d-dimensional bitmap.
1008                 * @param[in] at
1009                 * \return
1010                 */
1011              inline bool operator[](const typename UniformGrid::CellCoordinate &at)	{ return m_Mask[at.X()][at.Y()][at.Z()]; }
1012              inline bool& GetFlag(const int i, const int j, const int k)const				{ return m_Mask[i][j][k]; }
1013              inline void  SetFlat(const int i, const int j, const int k)						{ m_Mask[i][j][k] = true; }
1014 
1015 
1016              inline bool  GetFlag(const typename UniformGrid::CellCoordinate &at)	const			{ return m_Mask[at.X()][at.Y()][at.Z()]; }
1017              inline void  SetFlag(const typename UniformGrid::CellCoordinate &at)						{ m_Mask[at.X()][at.Y()][at.Z()] = true; }
1018 
1019          protected:
1020              std::vector< std::vector< std::vector< bool > > >
1021                         m_Mask;					/*!< The bitmap image.							*/
1022              int	m_Resolution;		/*!< The resolution of the bitmap.	*/
1023          }; // end of class BinaryImage
1024 
1025 
1026 
1027          /*******************************************************************************************************************************/
1028          /*! \struct Neighbor
1029          * This class is used to retrieve the neighboring objects in the spatial queries and to sort them.
1030          */
1031          /*******************************************************************************************************************************/
1032          struct Neighbor
1033          {
1034              /*!
1035               * Default constructor
1036                 */
1037              Neighbor()
1038              {
1039                  object		= NULL;
1040                  distance = ScalarType(-1.0);
1041                  nearest_point.SetZero();
1042              }
1043 
1044 
1045              /*!
1046               * Constructor
1047                 * @param[in] pObject	The pointer to the object.
1048                 * @param[in] dist			The distance between the object and the query point.
1049                 * @param[in] point		The point on the object having minimal distance from the query point.
1050                 */
1051              Neighbor(ObjectPointer pObject, ScalarType dist, CoordinateType point)
1052              {
1053                  object = pObject;
1054                  distance = dist;
1055                  nearest_point(point);
1056              }
1057 
1058 
1059              /*!
1060               * Less than operator. Needed for sorting a range of neighbor based on their distance from the query object.
1061                 */
1062              inline bool operator<(const Neighbor &second)
1063              {
1064                  return distance<second.distance;
1065              }
1066 
1067              ObjectPointer	object;
1068              ScalarType			distance;
1069              CoordinateType nearest_point;
1070          }; // end of struct Neighbor
1071 
1072         /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1073         /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1074         /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1075 
1076 
1077         /************************************************************************/
1078         /* Functions	                                                          */
1079         /************************************************************************/
1080     public:
1081         /*!
1082          * The hash table can be constructed following two different approaches:
1083          * the first is more fast, but might allocate a offset table bigger than the necessary
1084          * 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.
1085          */
1086         enum		ConstructionApproach { FastConstructionApproach=0, CompactConstructionApproach=1 };
1087 
1088         /*!
1089          * Default constructor
1090          */
1091         PerfectSpatialHashing() { srand( (unsigned) time(NULL) ); }
1092 
1093         /*!
1094          * Destructor
1095          */
1096         ~PerfectSpatialHashing() { /* ... I don't remember if there is something to delete! :D */ }
1097 
1098         template < class OBJECT_ITERATOR >
1099         void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj)
1100         { Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, NULL); }
1101 
1102         template < class OBJECT_ITERATOR >
1103         void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, vcg::CallBackPos *callback)
1104         { Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, callback); }
1105 
1106         template < class OBJECT_ITERATOR >
1107         void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach)
1108         { Set<OBJECT_ITERATOR>(bObj, eObj, approach, NULL); }
1109 
1110         /*!
1111          * Add the elements to the PerfectSpatialHashing data structure. Since this structure can handle only
1112          * static dataset, the elements mustn't be changed while using this structure.
1113          * @param[in] bObj			The iterator addressing the first element to be included in the hashing.
1114          * @param[in] eObj			The iterator addressing the position after the last element to be included in the hashing.
1115          * @param[in] approach	Either <code>FastConstructionApproach</code> or <code>CompactConstructionApproach</code>.
1116          * @param[in] callback  The callback to call to provide information about the progress of the computation.
1117          */
1118         template < class OBJECT_ITERATOR >
1119         void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach, vcg::CallBackPos *callback)
1120         {
1121             BoundingBoxType bounding_box;
1122             BoundingBoxType object_bb;
1123             bounding_box.SetNull();
1124             for (OBJECT_ITERATOR iObj=bObj; iObj!=eObj; iObj++)
1125             {
1126                 (*iObj).GetBBox(object_bb);
1127                 bounding_box.Add(object_bb);
1128             }
1129 
1130             //...and expand it a bit more
1131             BoundingBoxType	resulting_bb(bounding_box);
1132             CoordinateType	offset = bounding_box.Dim()*float(m_BOUNDING_BOX_EXPANSION_FACTOR);
1133             CoordinateType	center = bounding_box.Center();
1134             resulting_bb.Offset(offset);
1135             float longest_side = vcg::math::Max( resulting_bb.DimX(), vcg::math::Max(resulting_bb.DimY(), resulting_bb.DimZ()) )/2.0f;
1136             resulting_bb.Set(center);
1137             resulting_bb.Offset(longest_side);
1138 
1139             int number_of_objects = int(std::distance(bObj, eObj));
1140 
1141             // Try to find a reasonable space partition
1142 #ifdef _USE_GRID_UTIL_PARTIONING_
1143             vcg::Point3i resolution;
1144             vcg::BestDim<ScalarType>(number_of_objects, resulting_bb.Dim(), resolution);
1145             int cells_per_side = resolution.X();
1146 #else ifdef _USE_OCTREE_PARTITIONING_ // Alternative to find the resolution of the uniform grid:
1147             int primitives_per_voxel;
1148             int depth = 4;
1149             do
1150             {
1151                 int		number_of_voxel = 1<<(3*depth); // i.e. 8^depth
1152                 float density					= float(number_of_voxel)/float(depth);
1153                 primitives_per_voxel	= int(float(number_of_objects)/density);
1154                 depth++;
1155             }
1156             while (primitives_per_voxel>16 && depth<15);
1157             int cells_per_side = int(powf(2.0f, float(depth)));
1158 #endif
1159 
1160             m_UniformGrid.Allocate(resulting_bb, cells_per_side);
1161             m_UniformGrid.InsertElements(bObj, eObj);
1162             m_Bitmap.Allocate(cells_per_side);
1163             int number_of_cells_occupied = m_UniformGrid.GetNumberOfNotEmptyCells();
1164 
1165             int hash_table_size = (int) ceilf(powf(float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
1166             if (hash_table_size>256)
1167                 hash_table_size = (int) ceilf(powf(1.01f*float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
1168             m_HashTable.Allocate(hash_table_size);
1169 
1170             switch (approach)
1171             {
1172             case FastConstructionApproach		:	PerformFastConstruction(number_of_cells_occupied, callback)		; break;
1173             case CompactConstructionApproach: PerformCompactConstruction(number_of_cells_occupied, callback); break;
1174             default: assert(false);
1175             }
1176             Finalize();
1177         } // end of method Set
1178 
1179 
1180         /*!
1181          * Returns all the objects contained inside a specified sphere
1182          * @param[in]  distance_functor
1183          * @param[in]	 marker
1184          * @param[in]	 sphere_center
1185          * @param[in]	 sphere_radius
1186          * @param[out] objects
1187          * @param[out] distances
1188          * @param[out] points
1189          * \return
1190          */
1191         template <class OBJECT_POINT_DISTANCE_FUNCTOR, class OBJECT_MARKER, class OBJECT_POINTER_CONTAINER, class DISTANCE_CONTAINER, class POINT_CONTAINER>
1192         unsigned int GetInSphere
1193             (
1194             OBJECT_POINT_DISTANCE_FUNCTOR		&	distance_functor,
1195             OBJECT_MARKER										&	marker,
1196             const CoordinateType									&	sphere_center,
1197             const ScalarType								&	sphere_radius,
1198             OBJECT_POINTER_CONTAINER				&	objects,
1199             DISTANCE_CONTAINER							&	distances,
1200             POINT_CONTAINER									&	points,
1201             bool															sort_per_distance   = true,
1202             bool															allow_zero_distance = true
1203             )
1204         {
1205             BoundingBoxType query_bb(sphere_center, sphere_radius);
1206             vcg::Box3i integer_bb = m_UniformGrid.Interize(query_bb);
1207 
1208             vcg::Point3i index;
1209             std::vector< std::vector< ObjectPointer >* > contained_objects;
1210             std::vector< ObjectPointer >* tmp;
1211             for (index.X()=integer_bb.min.X(); index.X()<=integer_bb.max.X(); index.X()++)
1212                 for (index.Y()=integer_bb.min.Y(); index.Y()<=integer_bb.max.Y(); index.Y()++)
1213                     for (index.Z()=integer_bb.min.Z(); index.Z()<=integer_bb.max.Z(); index.Z()++)
1214                         if ((tmp=(*this)[index])!=NULL)
1215                             contained_objects.push_back(tmp);
1216 
1217             std::vector< Neighbor > results;
1218             for (typename std::vector< typename std::vector< ObjectPointer >* >::iterator iVec=contained_objects.begin(), eVec=contained_objects.end(); iVec!=eVec; iVec++)
1219                 for (typename std::vector< ObjectPointer >::iterator iObj=(*iVec)->begin(), eObj=(*iVec)->end(); iObj!=eObj; iObj++ )
1220                 {
1221                     int r = int(results.size());
1222                     results.push_back(Neighbor());
1223                     results[r].object		= *iObj;
1224                     results[r].distance = sphere_radius;
1225                     if (!distance_functor(*results[r].object, sphere_center, results[r].distance, results[r].nearest_point) || (results[r].distance==ScalarType(0.0) && !allow_zero_distance) )
1226                         results.pop_back();
1227                 }
1228 
1229             if (sort_per_distance)
1230                 std::sort( results.begin(), results.end() );
1231 
1232             int number_of_objects = int(results.size());
1233             distances.resize(number_of_objects);
1234             points.resize(number_of_objects);
1235             objects.resize(number_of_objects);
1236             for (int i=0, size=int(results.size()); i<size; i++)
1237             {
1238                 distances[i]	= results[i].distance;
1239                 points[i]			= results[i].nearest_point;
1240                 objects[i]		= results[i].object;
1241             }
1242             return number_of_objects;
1243         } //end of GetInSphere
1244 
1245 
1246         /*!
1247         * Once the offset table has been built, this function can be used to access the data.
1248         * Given a 3D point in the space, this function returns the set of ObjectPointers contained
1249         * in the same UniformGrid cell where this point is contained.
1250         * @param[in] query The 3D query point.
1251         * \return		The pointer to the vector of ObjectPointers contained in the same UG of query. NULL if any.
1252         */
1253         std::vector< ObjectPointer >* operator[](const CoordinateType &query)
1254         {
1255             typename UniformGrid::CellCoordinate ug_index = m_UniformGrid.Interize(query);
1256             if (!m_Bitmap[ug_index])
1257                 return NULL;
1258 
1259             typename HashTable::EntryCoordinate ht_index = PerfectHashFunction(ug_index);
1260             std::vector< ObjectPointer >* result = m_HashTable[ht_index];
1261             return result;
1262         }
1263 
1264 
1265         std::vector< ObjectPointer >* operator[](const typename UniformGrid::CellCoordinate  &query)
1266         {
1267             if(!m_Bitmap[query])
1268                 return NULL;
1269 
1270             typename HashTable::EntryCoordinate ht_index = PerfectHashFunction(query);
1271             std::vector< ObjectPointer >* result = m_HashTable[ht_index]->domain_data;
1272             return result;
1273         }
1274 
1275 
1276     protected:
1277         /*!
1278         * The injective mapping from the set of occupied cells to a slot in the hash-table
1279         * @param[in]	query		The index of a domain cell whose mapping has to be calculated.
1280         * @param[out]	result	The index of a entry in the hash-table where query is mapped into.
1281         */
1282         typename HashTable::EntryCoordinate PerfectHashFunction(const typename UniformGrid::CellCoordinate &query)
1283         {
1284             typename HashTable::EntryCoordinate result;
1285             typename OffsetTable::OffsetPointer offset = m_OffsetTable[ m_OffsetTable.DomainToOffsetTable(query) ];
1286             result = m_HashTable.DomainToHashTable( Shift(query, *offset) );
1287             return result;
1288         }
1289 
1290 
1291         /*!
1292          * Performs the addition between a entry coordinate and an offset.
1293          * @param[in] entry		The index of a given cell.
1294          * @param[in] offset	The offset that must be applied to the entry.
1295          * \return						The entry resulting by the addition of entry and offset.
1296          */
1297         typename HashTable::EntryCoordinate Shift(const vcg::Point3i &entry, const typename OffsetTable::Offset &offset)
1298         {
1299             typename HashTable::EntryCoordinate result;
1300             result.X() = entry.X() + int(offset.X());
1301             result.Y() = entry.Y() + int(offset.Y());
1302             result.Z() = entry.Z() + int(offset.Z());
1303             return result;
1304         }
1305 
1306 
1307         /*!
1308         * Finalizes the data structures at the end of the offset-table construction.
1309         * This function eliminates all unnecessary data, and encodes sparsity.
1310         * TODO At the moment, the sparsity encoding is implemented thought a bitmap, i.e. a boolean grid
1311         *			where each slot tells if the relative UniformGrid has a valid entry in the HashTable.
1312         */
1313         void Finalize()
1314         {
1315 #ifdef _DEBUG
1316             for (UniformGrid::EntryIterator iUGEntry=m_UniformGrid.Begin(), eUGEntry=m_UniformGrid.End(); iUGEntry!=eUGEntry; iUGEntry++)
1317                 assert(m_Bitmap.ContainsData(iUGEntry.GetPosition())==((*iUGEntry)->size()>0));
1318 #endif
1319             m_HashTable.Finalize();
1320             m_UniformGrid.Finalize();
1321             m_OffsetTable.Finalize();
1322         }
1323 
1324 
1325         /*!
1326          * Check if the given offset is valid for a set of domain cell.
1327          * @param[in] pre_image
1328          * @param[in] offset
1329          * \return
1330          */
1331         bool IsAValidOffset(const std::vector< typename UniformGrid::CellCoordinate > *pre_image, const typename OffsetTable::Offset &offset)
1332         {
1333             int ht_size			= m_HashTable.GetSize();
1334             int sqr_ht_size = ht_size*ht_size;
1335             std::vector< int > involved_entries;
1336             for (int i=0, pre_image_size=int((*pre_image).size()); i<pre_image_size; i++)
1337             {
1338                 typename UniformGrid::CellCoordinate domain_entry = (*pre_image)[i];
1339                 typename HashTable::EntryCoordinate  hash_entry		= m_HashTable.DomainToHashTable( Shift(domain_entry, offset) );
1340                 if (!m_HashTable.IsFree(hash_entry))
1341                     return false;
1342                 else
1343                     involved_entries.push_back(hash_entry.X()*sqr_ht_size + hash_entry.Y()*ht_size + hash_entry.Z());
1344             }
1345 
1346             // In order to guarantee that the PerfectHashFunction is injective, the image of each domain-entry must be unique in the hash-table
1347             std::sort(involved_entries.begin(), involved_entries.end());
1348             for (int i=0, j=1, size=int(involved_entries.size()); j<size; i++, j++)
1349                 if (involved_entries[i]==involved_entries[j])
1350                     return false;
1351 
1352             return true;
1353         }
1354 
1355 
1356         /*!
1357          * Given the size of the hash table and an initial seed for the size of the offset table, returns an appropriate size
1358          * for the offset table in order to avoid less effective constructions.
1359          * @param[in] hash_table_size		The number of entries for each side of the hash-table.
1360          * @param[in] offset_table_size The number of entries for each side of the offset-table.
1361          * \return											The next appropriate size for the offset table.
1362          */
1363         int GetUnefectiveOffsetTableSize(const int hash_table_size, const int offset_table_size)
1364         {
1365             int result = offset_table_size;
1366             while (GreatestCommonDivisor(hash_table_size, result)!=1 || hash_table_size%result==0)
1367                 result += (hash_table_size%2==0) ? (result%2)+1 : 1;  //Change the offset table size, otherwise the hash construction should be ineffective
1368             return result;
1369         }
1370 
1371 
1372         /*!
1373         * Start the construction of the offset table trying to complete as quickly as possible.
1374         * Sometimes the dimension of the offset table will not be minimal.
1375         * @param[in] number_of_filled_cells The number of entries in the uniform grid containing some elements of the dataset.
1376         */
1377         void PerformFastConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
1378         {
1379             int offset_table_size = (int) ceilf(powf(m_SIGMA()*float(number_of_filled_cells), 1.0f/float(m_DIMENSION())));
1380             int hash_table_size   = m_HashTable.GetSize();
1381             int failed_construction_count = 0;
1382             do
1383             {
1384                 offset_table_size += failed_construction_count++;
1385                 offset_table_size  = GetUnefectiveOffsetTableSize(hash_table_size, offset_table_size);
1386             }
1387             while(!OffsetTableConstructionSucceded(offset_table_size, callback));
1388         }
1389 
1390 
1391         /*!
1392         * Start the construction of the offset table trying to minimize its dimension.
1393         * The offset table size is chosen by performing a binary search over possible values for the offset table size.
1394         * For this reason, this method will generally require more time than PerformFastConstruction.
1395         */
1396         void PerformCompactConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
1397         {
1398             int min_successfully_dimension	= std::numeric_limits<int>::max();
1399             int hash_table_size							= m_HashTable.GetSize();
1400             int half_hash_table_size				= int(float(hash_table_size)/2.0f);
1401 
1402             // 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
1403             for (int t=0; t<m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION(); t++)
1404             {
1405                 int lower_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(double(rand())/double(RAND_MAX)*half_hash_table_size + 1) );
1406                 int upper_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(((double) rand() / (double) RAND_MAX) * hash_table_size + half_hash_table_size));
1407 
1408                 // The construction of the offset table using this pair of values surely succeed for max, but not for min
1409                 // The optimal value for the offset table size is then found via binary search over this pair of values
1410                 int candidate_offset_table_size;
1411                 int last_tried_size = -1;
1412                 while (lower_bound<upper_bound)
1413                 {
1414                     candidate_offset_table_size = GetUnefectiveOffsetTableSize(hash_table_size, int(floorf((lower_bound+upper_bound)/2.0f)));
1415 
1416                     // If in the previous iteration the offset table has been successfully constructed with the same size,
1417                     // a minimum for the range [lower_bound, upper_bound] has been found
1418                     if (last_tried_size==candidate_offset_table_size)
1419                         break;
1420 
1421                     if ( OffsetTableConstructionSucceded((last_tried_size=candidate_offset_table_size), callback) )
1422                     {
1423                         upper_bound = candidate_offset_table_size;
1424                         min_successfully_dimension = std::min(candidate_offset_table_size, min_successfully_dimension);
1425                     }
1426                     else
1427                         lower_bound = candidate_offset_table_size;
1428 
1429                     m_HashTable.Clear();
1430                     m_HashTable.BuildFreeEntryList();
1431                     m_OffsetTable.Clear();
1432                     m_Bitmap.Clear();
1433                 }
1434 #ifdef _DEBUD
1435                 printf("\nPerfectSpatialHashing: minimum offset table size found at the %d-th iteration was %d\n", (t+1), min_successfully_dimension);
1436 #endif
1437             }
1438 
1439             // Finally the OffsetTable must be constructed using the minimum valid size
1440             while (!OffsetTableConstructionSucceded(min_successfully_dimension, callback))
1441             {
1442                 m_HashTable.Clear();
1443                 m_HashTable.BuildFreeEntryList();
1444                 m_OffsetTable.Clear();
1445                 m_Bitmap.Clear();
1446             }
1447         }
1448 
1449 
1450 
1451         /*!
1452         * Try to construct the offset table for a given size
1453         *	\param[in] offset_table_size  The size of the offset table.
1454         *	\return												<CODE>true</CODE> if and only if the construction of the offset table succeeds.
1455         */
1456         bool OffsetTableConstructionSucceded(const int offset_table_size, vcg::CallBackPos *callback)
1457         {
1458             m_OffsetTable.Allocate(offset_table_size); // Create the Offset table
1459             m_OffsetTable.BuildH1PreImage(m_UniformGrid.Begin(), m_UniformGrid.End()); // Build the f0 pre-image
1460 
1461             std::list< typename OffsetTable::PreImage > preimage_slots;
1462             m_OffsetTable.GetPreImageSortedPerCardinality(preimage_slots);
1463 
1464             char msg[128];
1465             sprintf(msg, "Building offset table of resolution %d", m_OffsetTable.GetSize());
1466             int	step = int(preimage_slots.size())/100;
1467             int number_of_slots = int(preimage_slots.size());
1468             int perc = 0;
1469             int iter = 0;
1470             for (typename std::list< typename OffsetTable::PreImage >::iterator iPreImage=preimage_slots.begin(), ePreImage=preimage_slots.end(); iPreImage!=ePreImage; iPreImage++, iter++)
1471             {
1472                 if (callback!=NULL && iter%step==0 && (perc=iter*100/number_of_slots)<100) (*callback)(perc, msg);
1473 
1474                 bool found_valid_offset = false;
1475                 typename OffsetTable::Offset candidate_offset;
1476 
1477                 // Heuristic #1: try to set the offset value to one stored in a neighboring entry of the offset table
1478                 std::vector< typename OffsetTable::Offset > consistent_offsets;
1479                 m_OffsetTable.SuggestConsistentOffsets( (*iPreImage).entry_index, consistent_offsets);
1480                 for (typename std::vector< typename OffsetTable::Offset >::iterator iOffset=consistent_offsets.begin(), eOffset=consistent_offsets.end(); iOffset!=eOffset && !found_valid_offset; iOffset++)
1481                     if (IsAValidOffset(iPreImage->pre_image, *iOffset))
1482                     {
1483                         found_valid_offset = true;
1484                         candidate_offset = *iOffset;
1485                     }
1486 
1487 
1488                     // Heuristic #2:
1489                     if (!found_valid_offset)
1490                     {
1491                         std::vector< typename UniformGrid::CellCoordinate > *pre_image = (*iPreImage).pre_image;
1492                         for (typename std::vector< typename UniformGrid::CellCoordinate >::iterator iPreImage=pre_image->begin(), ePreImage=pre_image->end(); iPreImage!=ePreImage && !found_valid_offset; iPreImage++)
1493                             for (NeighboringEntryIterator iUGNeighbourhood=m_UniformGrid.GetNeighboringEntryIterator(*iPreImage); iUGNeighbourhood<6 && !found_valid_offset; iUGNeighbourhood++ )
1494                                 if (!m_OffsetTable.IsFree( m_OffsetTable.DomainToOffsetTable( *iUGNeighbourhood ) ))
1495                                 {
1496                                     typename HashTable::EntryCoordinate ht_entry = PerfectHashFunction(*iUGNeighbourhood);
1497                                     for (NeighboringEntryIterator iHTNeighbourhood=m_HashTable.GetNeighborintEntryIterator(ht_entry); iHTNeighbourhood<6 && !found_valid_offset; iHTNeighbourhood++)
1498                                         if (m_HashTable.IsFree(*iHTNeighbourhood))
1499                                         {
1500                                             candidate_offset.Import( *iHTNeighbourhood-m_HashTable.DomainToHashTable(*iPreImage) ) ;
1501                                             // m_OffsetTable.ValidateEntry( candidate_offset ); Is'n necessary, becouse the offset type is unsigned char.
1502                                             if (IsAValidOffset(pre_image, candidate_offset))
1503                                                 found_valid_offset = true;
1504                                         }
1505                                 }
1506                     }
1507 
1508                     if (!found_valid_offset)
1509                     {
1510                         // At the beginning, the offset can be found via random searches
1511                         for (int i=0; i<m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET() && !found_valid_offset; i++)
1512                         {
1513                             typename HashTable::EntryCoordinate base_entry = (*iPreImage).pre_image->at(0);
1514                             do
1515                                 m_OffsetTable.GetRandomOffset(candidate_offset);
1516                             while (!m_HashTable.IsFree( m_HashTable.DomainToHashTable( Shift(base_entry, candidate_offset) ) ));
1517 
1518                             if (IsAValidOffset( (*iPreImage).pre_image, candidate_offset))
1519                                 found_valid_offset = true;
1520                         }
1521 
1522                         // The chance to find a valid offset table via random searches decreases toward the end of the offset table construction:
1523                         // So a exhaustive search over all the free hash table entries is performed.
1524                         for (typename std::list< typename HashTable::EntryCoordinate >::const_iterator iFreeCell=m_HashTable.GetFreeEntryList()->begin(), eFreeCell=m_HashTable.GetFreeEntryList()->end(); iFreeCell!=eFreeCell && !found_valid_offset; iFreeCell++)
1525                         {
1526                             typename UniformGrid::CellCoordinate domain_entry		= (*iPreImage).pre_image->at(0);
1527                             typename OffsetTable::EntryCoordinate offset_entry		= m_OffsetTable.DomainToOffsetTable(domain_entry);
1528                             typename HashTable::EntryCoordinate hashtable_entry	= m_HashTable.DomainToHashTable(domain_entry);
1529                             candidate_offset.Import(*iFreeCell - hashtable_entry);
1530 
1531                             if ( IsAValidOffset(iPreImage->pre_image, candidate_offset) )
1532                                 found_valid_offset = true;
1533                         }
1534                 }
1535 
1536                     // If a valid offset has been found, the construction of the offset table continues,
1537                     // otherwise the offset table must be enlarged and the construction repeated
1538                     if (found_valid_offset)
1539                     {
1540                         m_OffsetTable.SetOffset( (*iPreImage->pre_image).at(0), candidate_offset);
1541                         for (int c=0, pre_image_cardinality = iPreImage->cardinality; c<pre_image_cardinality; c++)
1542                         {
1543                             typename HashTable::EntryCoordinate ht_entry = PerfectHashFunction( (*iPreImage->pre_image).at(c));
1544                             std::vector< ObjectPointer > *domain_data = m_UniformGrid[ (*iPreImage->pre_image).at(c) ];
1545                             m_HashTable.SetEntry(ht_entry, domain_data /*, (*iPreImage->pre_image).at(c)*/); // might be useful for encoding sparsity
1546                             m_Bitmap.SetFlag((*iPreImage->pre_image).at(c));
1547                         }
1548                     }
1549                     else
1550                     {
1551                         m_OffsetTable.Clear();
1552                         m_HashTable.Clear();
1553                         m_HashTable.BuildFreeEntryList();
1554                         m_Bitmap.Clear();
1555                         return false;
1556                     }
1557             }
1558 
1559             if (callback!=NULL) (*callback)(100, msg);
1560             return true;
1561         } // end of OffsetTableConstructionSucceded
1562 
1563 
1564         /************************************************************************/
1565         /* Data Members                                                         */
1566         /************************************************************************/
1567     protected:
1568         UniformGrid m_UniformGrid;	/*!< The uniform grid used for partitioning the volume.						*/
1569         OffsetTable m_OffsetTable;	/*!< The offset table corresponding to \f$\Phi\f$ in the article. */
1570         HashTable		m_HashTable;		/*!< The hash table that will substitute the uniform grid.				*/
1571         BinaryImage	m_Bitmap;
1572 
1573         const static float  m_BOUNDING_BOX_EXPANSION_FACTOR() { return SCALAR_TYPE(0.035); }
1574         const static float  m_SIGMA() {return SCALAR_TYPE(1.0f/(2.0f*SCALAR_TYPE(m_DIMENSION)));}
1575         const static int    m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION () { return 5; }
1576         const static int    m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET() {return 32;}
1577         const static int    m_DIMENSION() {return 3;}
1578     }; //end of class PerfectSpatialHashing
1579 
1580     /*! @} */
1581     //end of Doxygen documentation
1582 }//end of namespace vcg
1583 
1584 
1585 #endif //VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
1586