00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
00025 #define VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
00026
00027 #pragma warning(disable : 4996)
00028
00029 #define _USE_GRID_UTIL_PARTIONING_ 1
00030 #define _USE_OCTREE_PARTITIONING_ (1-_USE_GRID_UTIL_PARTIONING_)
00031
00032 #include <vector>
00033 #include <list>
00034 #include <algorithm>
00035
00036 #include <vcg/space/index/base.h>
00037 #include <vcg/space/index/grid_util.h>
00038
00039 #include <vcg/space/point2.h>
00040 #include <vcg/space/point3.h>
00041 #include <vcg/space/box3.h>
00042
00043 namespace vcg
00044 {
00045
00046
00047 int GreatestCommonDivisor(const int a, const int b)
00048 {
00049 int m = a;
00050 int n = b;
00051
00052 do
00053 {
00054 if (m<n) std::swap(m, n);
00055 m = m % n;
00056 std::swap(m, n);
00057 }
00058 while (n!=0);
00059 return m;
00060 }
00061
00062
00075 template < class OBJECT_TYPE, class SCALAR_TYPE >
00076 class PerfectSpatialHashing : public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE >
00077 {
00078
00079 template < typename TYPE >
00080 struct Dereferencer
00081 {
00082 static TYPE& Reference(TYPE &t) { return t; }
00083 static TYPE& Reference(TYPE* &t) { return *t; }
00084 static const TYPE& Reference(const TYPE &t) { return t; }
00085 static const TYPE& Reference(const TYPE* &t) { return *t; }
00086 };
00087
00088
00089 template < typename TYPE >
00090 struct ReferenceType { typedef TYPE Type; };
00091
00092
00093 template < typename TYPE >
00094 struct ReferenceType< TYPE * > { typedef typename ReferenceType<TYPE>::Type Type; };
00095
00096 public:
00097 typedef SCALAR_TYPE ScalarType;
00098 typedef OBJECT_TYPE ObjectType;
00099 typedef typename ReferenceType< ObjectType >::Type * ObjectPointer;
00100 typedef typename vcg::Box3< ScalarType > BoundingBoxType;
00101 typedef typename vcg::Point3< ScalarType > CoordinateType;
00102
00103 protected:
00107 struct NeighboringEntryIterator
00108 {
00114 NeighboringEntryIterator(const vcg::Point3i &entry, const int table_size)
00115 {
00116 m_Center = entry;
00117 m_TableSize = table_size;
00118 m_CurrentNeighbor.X() = (m_Center.X()+m_TableSize-1)%m_TableSize;
00119 m_CurrentNeighbor.Y() = m_Center.Y();
00120 m_CurrentNeighbor.Z() = m_Center.Z();
00121 m_CurrentIteration = 0;
00122 }
00123
00127 void operator++(int)
00128 {
00129 switch(++m_CurrentIteration)
00130 {
00131 case 1: m_CurrentNeighbor.X()=(m_Center.X()+1)%m_TableSize; break;
00132 case 2: m_CurrentNeighbor.X()=m_Center.X(); m_CurrentNeighbor.Y()=(m_Center.Y()+m_TableSize-1)%m_TableSize; break;
00133 case 3: m_CurrentNeighbor.Y()=(m_Center.Y()+1)%m_TableSize; break;
00134 case 4: m_CurrentNeighbor.Y()=m_Center.Y(); m_CurrentNeighbor.Z()=(m_Center.Z()+m_TableSize-1)%m_TableSize; break;
00135 case 5: m_CurrentNeighbor.Z()=(m_Center.Z()+1)%m_TableSize; break;
00136 default: m_CurrentNeighbor = vcg::Point3i(-1, -1, -1); break;
00137 }
00138 }
00139
00144 vcg::Point3i operator*() { return m_CurrentNeighbor; }
00145
00151 NeighboringEntryIterator& operator =(const NeighboringEntryIterator &it)
00152 {
00153 m_Center = it.m_Center ;
00154 m_CurrentNeighbor = it.m_CurrentNeighbor ;
00155 m_CurrentIteration = it.m_CurrentIteration ;
00156 m_TableSize = it.m_TableSize ;
00157 return *this;
00158 }
00159
00164 inline bool operator <(const int value) { return m_CurrentIteration<value; }
00165
00166 protected:
00167 vcg::Point3i m_Center;
00168 vcg::Point3i m_CurrentNeighbor;
00169 int m_CurrentIteration;
00170 int m_TableSize;
00171 };
00172
00173
00174
00175
00180
00181 class UniformGrid
00182 {
00183 public:
00184 typedef vcg::Point3i CellCoordinate;
00185
00189 struct EntryIterator
00190 {
00191 friend class UniformGrid;
00192
00196 EntryIterator(UniformGrid *uniform_grid, const CellCoordinate &position)
00197 {
00198 m_UniformGrid = uniform_grid;
00199 m_CurrentPosition = position;
00200 }
00201
00202
00206 void operator++(int)
00207 {
00208 if (++m_CurrentPosition.Z()==m_UniformGrid->GetResolution())
00209 {
00210 m_CurrentPosition.Z() = 0;
00211 if (++m_CurrentPosition.Y()==m_UniformGrid->GetResolution())
00212 {
00213 m_CurrentPosition.Y() = 0;
00214 if (++m_CurrentPosition.X()==m_UniformGrid->GetResolution())
00215 m_CurrentPosition = CellCoordinate(-1, -1, -1);
00216 }
00217 }
00218 }
00219
00220
00225 void operator =(const EntryIterator &it)
00226 {
00227 m_UniformGrid = it.m_UniformGrid;
00228 m_CurrentPosition = it.m_CurrentPosition;
00229 }
00230
00234 bool operator==(const EntryIterator &it) const
00235 {
00236 return m_CurrentPosition==it.m_CurrentPosition;
00237 }
00238
00242 bool operator!=(const EntryIterator &it) const
00243 {
00244 return m_CurrentPosition!=it.m_CurrentPosition;
00245 }
00246
00251 std::vector< ObjectPointer >* operator*()
00252 {
00253 return m_UniformGrid->GetObjects(m_CurrentPosition);
00254 }
00255
00259 CellCoordinate GetPosition() const
00260 {
00261 return m_CurrentPosition;
00262 }
00263
00264
00265 protected:
00266 UniformGrid * m_UniformGrid;
00267 CellCoordinate m_CurrentPosition;
00268 };
00269
00270
00274 UniformGrid() {}
00275
00279 ~UniformGrid() {}
00280
00281
00285 EntryIterator Begin() { return EntryIterator(this, CellCoordinate( 0, 0, 0)); }
00286 EntryIterator End() { return EntryIterator(this, CellCoordinate(-1, -1, -1)); }
00287
00288
00294 NeighboringEntryIterator GetNeighboringEntryIterator(const CellCoordinate &at) { return NeighboringEntryIterator(at, m_CellPerSide); }
00295
00296
00302 void Allocate(const BoundingBoxType &bounding_box, const int cell_per_side)
00303 {
00304 m_CellPerSide = cell_per_side;
00305 m_BoundingBox = bounding_box;
00306 m_CellSize = (m_BoundingBox.max - m_BoundingBox.min)/ScalarType(cell_per_side);
00307
00308 m_Grid.resize(m_CellPerSide);
00309 for (int i=0; i<m_CellPerSide; i++)
00310 {
00311 m_Grid[i].resize(m_CellPerSide);
00312 for (int j=0; j<m_CellPerSide; j++)
00313 m_Grid[i][j].resize(m_CellPerSide);
00314 }
00315 }
00316
00317
00321 void Finalize()
00322 {
00323 m_Grid.clear();
00324 }
00325
00326
00332 template < class OBJECT_ITERATOR >
00333 void InsertElements(const OBJECT_ITERATOR &begin, const OBJECT_ITERATOR &end)
00334 {
00335 typedef OBJECT_ITERATOR ObjectIterator;
00336 typedef Dereferencer< typename ReferenceType< typename OBJECT_ITERATOR::value_type >::Type > ObjectDereferencer;
00337
00338 std::vector< CellCoordinate > cells_occupied;
00339 for (ObjectIterator iObject=begin; iObject!=end; iObject++)
00340 {
00341 ObjectPointer pObject = &ObjectDereferencer::Reference( *iObject );
00342 GetCellsIndex( pObject, cells_occupied);
00343 for (std::vector< CellCoordinate >::iterator iCell=cells_occupied.begin(), eCell=cells_occupied.end(); iCell!=eCell; iCell++)
00344 GetObjects( *iCell )->push_back( pObject );
00345 cells_occupied.clear();
00346 }
00347 }
00348
00349
00355 inline CellCoordinate Interize(const CoordinateType &query) const
00356 {
00357 CellCoordinate result;
00358 result.X() = (int) floorf( (query.X()-m_BoundingBox.min.X())/m_CellSize.X() );
00359 result.Y() = (int) floorf( (query.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
00360 result.Z() = (int) floorf( (query.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
00361 return result;
00362 }
00363
00369 inline vcg::Box3i Interize(const BoundingBoxType &bounding_box) const
00370 {
00371 vcg::Box3i result;
00372 result.min = Interize(bounding_box.min);
00373 result.max = Interize(bounding_box.max);
00374 return result;
00375 }
00376
00377
00383 void GetCellsIndex(const ObjectPointer pObject, std::vector< CellCoordinate > & cells_occupied)
00384 {
00385 BoundingBoxType object_bb;
00386 (*pObject).GetBBox(object_bb);
00387 CoordinateType corner = object_bb.min;
00388
00389 while (object_bb.IsIn(corner))
00390 {
00391 CellCoordinate cell_index;
00392 cell_index.X() = (int) floorf( (corner.X()-m_BoundingBox.min.X())/m_CellSize.X() );
00393 cell_index.Y() = (int) floorf( (corner.Y()-m_BoundingBox.min.Y())/m_CellSize.Y() );
00394 cell_index.Z() = (int) floorf( (corner.Z()-m_BoundingBox.min.Z())/m_CellSize.Z() );
00395 cells_occupied.push_back( cell_index );
00396
00397 if ((corner.X()+=m_CellSize.X())>object_bb.max.X())
00398 {
00399 corner.X() = object_bb.min.X();
00400 if ( (corner.Z()+=m_CellSize.Z())>object_bb.max.Z() )
00401 {
00402 corner.Z() = object_bb.min.Z();
00403 corner.Y() += m_CellSize.Y();
00404 }
00405 }
00406 }
00407 }
00408
00409
00414 int GetNumberOfNotEmptyCells()
00415 {
00416 int number_of_not_empty_cell = 0;
00417 for (int i=0; i<m_CellPerSide; i++)
00418 for (int j=0; j<m_CellPerSide; j++)
00419 for (int k=0; k<m_CellPerSide; k++)
00420 if (GetObjects(i, j, k)->size()>0)
00421 number_of_not_empty_cell++;
00422 return number_of_not_empty_cell;
00423 }
00424
00429 inline int GetResolution() const { return m_CellPerSide; }
00430
00431
00437 std::vector< ObjectPointer >* GetObjects(const int i, const int j, const int k) { return &m_Grid[i][j][k]; }
00438 std::vector< ObjectPointer >* GetObjects(const CellCoordinate &at) { return &m_Grid[at.X()][at.Y()][at.Z()];}
00439 std::vector< ObjectPointer >* operator[](const CellCoordinate &at) { return &m_Grid[at.X()][at.Y()][at.Z()];}
00440
00441 protected:
00442 std::vector< std::vector< std::vector< std::vector< ObjectPointer > > > >
00443 m_Grid;
00444 BoundingBoxType m_BoundingBox;
00445 int m_CellPerSide;
00446 CoordinateType m_CellSize;
00447 };
00448
00449
00450
00451
00452
00456
00457 class HashTable
00458 {
00459 public:
00460 typedef vcg::Point3i EntryCoordinate;
00461
00462
00463
00464 struct Data
00465 {
00469 Data(std::vector< ObjectPointer > *data)
00470 {
00471 domain_data = data;
00472 }
00473
00474 std::vector< ObjectPointer > *domain_data;
00475 };
00476
00480 HashTable() {}
00481
00485 ~HashTable() { Clear(true); }
00486
00487
00491 NeighboringEntryIterator GetNeighborintEntryIterator(const EntryCoordinate &at) { return NeighboringEntryIterator(at, m_EntryPerSide); }
00492
00493
00498 void Allocate(const int entry_per_side)
00499 {
00500 m_EntryPerSide = entry_per_side;
00501 m_Table.resize(m_EntryPerSide);
00502 for (int i=0; i<m_EntryPerSide; i++)
00503 {
00504 m_Table[i].resize(m_EntryPerSide);
00505 for (int j=0; j<m_EntryPerSide; j++)
00506 m_Table[i][j].resize(m_EntryPerSide, NULL);
00507 }
00508
00509 BuildFreeEntryList();
00510 }
00511
00512
00513
00514
00515
00516
00517
00518 void Finalize()
00519 {
00520 Data *pData;
00521 for (int i=0; i<m_EntryPerSide; i++)
00522 for (int j=0; j<m_EntryPerSide; j++)
00523 for (int k=0; k<m_EntryPerSide; k++)
00524 if ((pData=GetData(i, j, k))!=NULL)
00525 {
00526 std::vector< ObjectPointer > *domain_data = pData->domain_data;
00527 pData->domain_data = new std::vector< ObjectPointer>( *domain_data );
00528 }
00529
00530
00531 m_FreeEntries.clear();
00532 }
00533
00534
00539 void BuildFreeEntryList()
00540 {
00541 m_FreeEntries.clear();
00542 for (int i=0; i<m_EntryPerSide; i++)
00543 for (int j=0; j<m_EntryPerSide; j++)
00544 for (int k=0; k<m_EntryPerSide; k++)
00545 {
00546 assert(m_Table[i][j][k]==NULL);
00547 m_FreeEntries.push_back(EntryCoordinate(i, j, k));
00548 }
00549 }
00550
00554 void Clear(bool delete_vectors=false)
00555 {
00556 for (int i=0; i<m_EntryPerSide; i++)
00557 for (int j=0; j<m_EntryPerSide; j++)
00558 for (int k=0; k<m_EntryPerSide; k++)
00559 if (m_Table[i][j][k]!=NULL)
00560 {
00561 if (delete_vectors)
00562 delete m_Table[i][j][k]->domain_data;
00563
00564 delete m_Table[i][j][k];
00565 m_Table[i][j][k] = NULL;
00566 }
00567
00568 m_FreeEntries.clear();
00569 }
00570
00575 std::list< EntryCoordinate >* GetFreeEntryList() { return &m_FreeEntries; }
00576
00581 EntryCoordinate DomainToHashTable(const typename UniformGrid::CellCoordinate &p)
00582 {
00583 EntryCoordinate result;
00584 result.X() = p.X()%m_EntryPerSide;
00585 result.Y() = p.Y()%m_EntryPerSide;
00586 result.Z() = p.Z()%m_EntryPerSide;
00587 return result;
00588 }
00589
00595 void SetEntry(const EntryCoordinate &at, std::vector< ObjectPointer > *data)
00596 {
00597 assert(IsFree(at));
00598 m_Table[at.X()][at.Y()][at.Z()] = new Data(data);
00599 m_FreeEntries.remove(at);
00600 }
00601
00607 void ValidateEntry(EntryCoordinate &entry)
00608 {
00609 while (entry.X()<0) entry.X()+=m_EntryPerSide;
00610 while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
00611 while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
00612 }
00613
00619 inline bool IsFree(const EntryCoordinate &at) const
00620 {
00621 return (GetData(at)==NULL);
00622 }
00623
00626 inline int GetSize() { return m_EntryPerSide; }
00627
00631 inline int GetNumberOfFreeEntries()
00632 {
00633 return int(m_FreeEntries.size());
00634 }
00635
00639 inline int GetNumberOfNotEmptyEntries()
00640 {
00641 return (int(powf(float(m_EntryPerSide), 3.0f))-int(m_FreeEntries.size()));
00642 }
00643
00650 inline Data* GetData (const int i, const int j, const int k) const { return m_Table[i][j][k]; }
00651 inline Data* GetData (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
00652 inline Data* operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
00653
00654 protected:
00655 int m_EntryPerSide;
00656 std::vector< std::vector< std::vector < Data* > > > m_Table;
00657 std::list< EntryCoordinate > m_FreeEntries;
00658 };
00659
00660
00664
00665 class OffsetTable
00666 {
00667 public:
00668 typedef unsigned char OffsetType;
00669 typedef vcg::Point3<OffsetType> Offset;
00670 typedef Offset * OffsetPointer;
00671 typedef vcg::Point3i EntryCoordinate;
00672
00677 struct PreImage
00678 {
00684 PreImage(EntryCoordinate &at, std::vector< typename UniformGrid::CellCoordinate > *preimage)
00685 {
00686 entry_index = at;
00687 pre_image = preimage;
00688 cardinality = int(pre_image->size());
00689 }
00690
00696 inline bool operator<(const PreImage &second) const { return (cardinality>second.cardinality); }
00697
00698
00699 std::vector< typename UniformGrid::CellCoordinate >
00700 * pre_image;
00701 EntryCoordinate entry_index;
00702 int cardinality;
00703 };
00704
00705
00709 OffsetTable() { m_EntryPerSide=-1; m_NumberOfOccupiedEntries=0;}
00710
00714 ~OffsetTable() { Clear(); }
00715
00719 void Clear()
00720 {
00721 for (int i=0; i<m_EntryPerSide; i++)
00722 for (int j=0; j<m_EntryPerSide; j++)
00723 for (int k=0; k<m_EntryPerSide; k++)
00724 if (m_Table[i][j][k]!=NULL)
00725 {
00726 delete m_Table[i][j][k];
00727 m_Table[i][j][k] = NULL;
00728 }
00729 m_EntryPerSide = -1;
00730 m_H1PreImage.clear();
00731 m_NumberOfOccupiedEntries = 0;
00732 }
00733
00739 void Allocate(int size)
00740 {
00741 m_NumberOfOccupiedEntries = 0;
00742
00743 m_EntryPerSide = size;
00744 m_Table.resize(m_EntryPerSide);
00745 for (int i=0; i<m_EntryPerSide; i++)
00746 {
00747 m_Table[i].resize(m_EntryPerSide);
00748 for (int j=0; j<m_EntryPerSide; j++)
00749 m_Table[i][j].resize(m_EntryPerSide, NULL);
00750 }
00751
00752 m_H1PreImage.resize(m_EntryPerSide);
00753 for (int i=0; i<m_EntryPerSide; i++)
00754 {
00755 m_H1PreImage[i].resize(m_EntryPerSide);
00756 for (int j=0; j<m_EntryPerSide; j++)
00757 m_H1PreImage[i][j].resize(m_EntryPerSide);
00758 }
00759 }
00760
00761
00765 void Finalize()
00766 {
00767 m_H1PreImage.clear();
00768 }
00769
00770
00775 void BuildH1PreImage(const typename UniformGrid::EntryIterator &begin, const typename UniformGrid::EntryIterator &end)
00776 {
00777 for (typename UniformGrid::EntryIterator iter=begin; iter!=end; iter++)
00778 {
00779 if ((*iter)->size()==0)
00780 continue;
00781
00782 typename UniformGrid::CellCoordinate cell_index = iter.GetPosition();
00783 EntryCoordinate at = DomainToOffsetTable(cell_index);
00784 m_H1PreImage[at.X()][at.Y()][at.Z()].push_back(cell_index);
00785 }
00786 }
00787
00792 void GetPreImageSortedPerCardinality(std::list< PreImage > &pre_image)
00793 {
00794 pre_image.clear();
00795 for (int i=0; i<m_EntryPerSide; i++)
00796 for (int j=0; j<m_EntryPerSide; j++)
00797 for (int k=0; k<m_EntryPerSide; k++)
00798 {
00799 std::vector< typedef UniformGrid::CellCoordinate > *preimage = &m_H1PreImage[i][j][k];
00800 if (preimage->size()>0)
00801 pre_image.push_back( PreImage(typename UniformGrid::CellCoordinate(i, j, k), preimage) );
00802 }
00803 pre_image.sort();
00804 }
00805
00806
00813 void SuggestConsistentOffsets(const EntryCoordinate &at, std::vector< Offset > &offsets)
00814 {
00815 offsets.clear();
00816 for (int i=-1; i<2; i++)
00817 for (int j=-1; j<2; j++)
00818 for (int k=-1; k<2; k++)
00819 {
00820 if (i==0 && j==0 && k==0)
00821 continue;
00822
00823 int x = (at.X()+i+m_EntryPerSide)%m_EntryPerSide;
00824 int y = (at.Y()+j+m_EntryPerSide)%m_EntryPerSide;
00825 int z = (at.Z()+k+m_EntryPerSide)%m_EntryPerSide;
00826 EntryCoordinate neighboring_entry(x, y, z);
00827 if (!IsFree(neighboring_entry))
00828 offsets.push_back( *GetOffset(neighboring_entry) );
00829 }
00830 }
00831
00832
00837 void ValidateEntryCoordinate(EntryCoordinate &entry)
00838 {
00839 while (entry.X()<0) entry.X()+=m_EntryPerSide;
00840 while (entry.Y()<0) entry.Y()+=m_EntryPerSide;
00841 while (entry.Z()<0) entry.Z()+=m_EntryPerSide;
00842 }
00843
00850 EntryCoordinate DomainToOffsetTable(const typename UniformGrid::CellCoordinate &coord)
00851 {
00852 EntryCoordinate result;
00853 result.X() = coord.X()%m_EntryPerSide;
00854 result.Y() = coord.Y()%m_EntryPerSide;
00855 result.Z() = coord.Z()%m_EntryPerSide;
00856 return result;
00857 }
00858
00864 void SetOffset(const typename UniformGrid::CellCoordinate &coord, const Offset &offset)
00865 {
00866 EntryCoordinate entry = DomainToOffsetTable(coord);
00867 assert(IsFree(entry));
00868 m_Table[entry.X()][entry.Y()][entry.Z()] = new Offset(offset);
00869 m_NumberOfOccupiedEntries++;
00870 }
00871
00877 void GetRandomOffset( Offset &offset )
00878 {
00879 offset.X() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
00880 offset.Y() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
00881 offset.Z() = OffsetType(rand()%m_MAX_VERSOR_LENGTH);
00882 }
00883
00884
00889 inline int GetSize() const {return m_EntryPerSide;}
00890
00891
00897 inline bool IsFree(const EntryCoordinate &at) const { return GetOffset(at)==NULL; }
00898
00899
00900
00905 inline int GetNumberOfOccupiedCells() const { return m_NumberOfOccupiedEntries; }
00906
00910 inline OffsetPointer& GetOffset (const int i, const int j, const int k) { return m_Table[i][j][k]; }
00911 inline OffsetPointer GetOffset (const int i, const int j, const int k) const { return m_Table[i][j][k]; }
00912
00913 inline OffsetPointer& GetOffset (const EntryCoordinate &at) { return m_Table[at.X()][at.Y()][at.Z()]; }
00914 inline OffsetPointer GetOffset (const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
00915
00916 inline OffsetPointer& operator[](const EntryCoordinate &at) { return m_Table[at.X()][at.Y()][at.Z()]; }
00917 inline OffsetPointer operator[](const EntryCoordinate &at) const { return m_Table[at.X()][at.Y()][at.Z()]; }
00918
00919 protected:
00920 const static int m_MAX_VERSOR_LENGTH = 256;
00921 int m_EntryPerSide;
00922 int m_NumberOfOccupiedEntries;
00923 std::vector< std::vector< std::vector< OffsetPointer > > > m_Table;
00924 std::vector< std::vector< std::vector< std::vector< typename UniformGrid::CellCoordinate > > > > m_H1PreImage;
00925 };
00926
00927
00928
00929
00934
00935 class BinaryImage
00936 {
00937 public:
00941 BinaryImage()
00942 {
00943 m_Resolution = -1;
00944 }
00945
00946
00950 ~BinaryImage() {}
00951
00952
00957 void Allocate(const int size)
00958 {
00959 m_Resolution = size;
00960 m_Mask.resize(m_Resolution);
00961 for (int i=0; i<m_Resolution; i++)
00962 {
00963 m_Mask[i].resize(m_Resolution);
00964 for (int j=0; j<m_Resolution; j++)
00965 m_Mask[i][j].resize(m_Resolution, false);
00966 }
00967 }
00968
00969
00973 void Clear()
00974 {
00975 for (int i=0; i<m_Resolution; i++)
00976 for (int j=0; j<m_Resolution; j++)
00977 std::fill(m_Mask[i][j].begin(), m_Mask[i][j].end(), false);
00978 }
00979
00980
00986 inline bool ContainsData(const typename UniformGrid::CellCoordinate &at) const { return GetFlag(at)==true;}
00987
00988
00993 inline int GetResolution() const { return m_Resolution; }
00994
00995
01003 inline bool operator()(const int i, const int j, const int k) { return m_Mask[i][j][k]; }
01004
01005
01011 inline bool operator[](const typename UniformGrid::CellCoordinate &at) { return m_Mask[at.X()][at.Y()][at.Z()]; }
01012 inline bool& GetFlag(const int i, const int j, const int k)const { return m_Mask[i][j][k]; }
01013 inline void SetFlat(const int i, const int j, const int k) { m_Mask[i][j][k] = true; }
01014
01015
01016 inline bool GetFlag(const typename UniformGrid::CellCoordinate &at) const { return m_Mask[at.X()][at.Y()][at.Z()]; }
01017 inline void SetFlag(const typename UniformGrid::CellCoordinate &at) { m_Mask[at.X()][at.Y()][at.Z()] = true; }
01018
01019 protected:
01020 std::vector< std::vector< std::vector< bool > > >
01021 m_Mask;
01022 int m_Resolution;
01023 };
01024
01025
01026
01027
01031
01032 struct Neighbor
01033 {
01037 Neighbor()
01038 {
01039 object = NULL;
01040 distance = ScalarType(-1.0);
01041 nearest_point.SetZero();
01042 }
01043
01044
01051 Neighbor(ObjectPointer pObject, ScalarType dist, CoordinateType point)
01052 {
01053 object = pObject;
01054 distance = dist;
01055 nearest_point(point);
01056 }
01057
01058
01062 inline bool operator<(const Neighbor &second)
01063 {
01064 return distance<second.distance;
01065 }
01066
01067 ObjectPointer object;
01068 ScalarType distance;
01069 CoordinateType nearest_point;
01070 };
01071
01075
01076
01077
01078
01079
01080 public:
01086 enum ConstructionApproach { FastConstructionApproach=0, CompactConstructionApproach=1 };
01087
01091 PerfectSpatialHashing() { srand( (unsigned) time(NULL) ); }
01092
01096 ~PerfectSpatialHashing() { }
01097
01098 template < class OBJECT_ITERATOR >
01099 void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj)
01100 { Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, NULL); }
01101
01102 template < class OBJECT_ITERATOR >
01103 void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, vcg::CallBackPos *callback)
01104 { Set<OBJECT_ITERATOR>(bObj, eObj, FastConstructionApproach, callback); }
01105
01106 template < class OBJECT_ITERATOR >
01107 void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach)
01108 { Set<OBJECT_ITERATOR>(bObj, eObj, approach, NULL); }
01109
01118 template < class OBJECT_ITERATOR >
01119 void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj, const ConstructionApproach approach, vcg::CallBackPos *callback)
01120 {
01121 BoundingBoxType bounding_box;
01122 BoundingBoxType object_bb;
01123 bounding_box.SetNull();
01124 for (OBJECT_ITERATOR iObj=bObj; iObj!=eObj; iObj++)
01125 {
01126 (*iObj).GetBBox(object_bb);
01127 bounding_box.Add(object_bb);
01128 }
01129
01130
01131 BoundingBoxType resulting_bb(bounding_box);
01132 CoordinateType offset = bounding_box.Dim()*float(m_BOUNDING_BOX_EXPANSION_FACTOR);
01133 CoordinateType center = bounding_box.Center();
01134 resulting_bb.Offset(offset);
01135 float longest_side = vcg::math::Max( resulting_bb.DimX(), vcg::math::Max(resulting_bb.DimY(), resulting_bb.DimZ()) )/2.0f;
01136 resulting_bb.Set(center);
01137 resulting_bb.Offset(longest_side);
01138
01139 int number_of_objects = int(std::distance(bObj, eObj));
01140
01141
01142 #ifdef _USE_GRID_UTIL_PARTIONING_
01143 vcg::Point3i resolution;
01144 vcg::BestDim<ScalarType>(number_of_objects, resulting_bb.Dim(), resolution);
01145 int cells_per_side = resolution.X();
01146 #else ifdef _USE_OCTREE_PARTITIONING_ // Alternative to find the resolution of the uniform grid:
01147 int primitives_per_voxel;
01148 int depth = 4;
01149 do
01150 {
01151 int number_of_voxel = 1<<(3*depth);
01152 float density = float(number_of_voxel)/float(depth);
01153 primitives_per_voxel = int(float(number_of_objects)/density);
01154 depth++;
01155 }
01156 while (primitives_per_voxel>16 && depth<15);
01157 int cells_per_side = int(powf(2.0f, float(depth)));
01158 #endif
01159
01160 m_UniformGrid.Allocate(resulting_bb, cells_per_side);
01161 m_UniformGrid.InsertElements(bObj, eObj);
01162 m_Bitmap.Allocate(cells_per_side);
01163 int number_of_cells_occupied = m_UniformGrid.GetNumberOfNotEmptyCells();
01164
01165 int hash_table_size = (int) ceilf(powf(float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
01166 if (hash_table_size>256)
01167 hash_table_size = (int) ceilf(powf(1.01f*float(number_of_cells_occupied), 1.0f/float(m_DIMENSION)));
01168 m_HashTable.Allocate(hash_table_size);
01169
01170 switch (approach)
01171 {
01172 case FastConstructionApproach : PerformFastConstruction(number_of_cells_occupied, callback) ; break;
01173 case CompactConstructionApproach: PerformCompactConstruction(number_of_cells_occupied, callback); break;
01174 default: assert(false);
01175 }
01176 Finalize();
01177 }
01178
01179
01191 template <class OBJECT_POINT_DISTANCE_FUNCTOR, class OBJECT_MARKER, class OBJECT_POINTER_CONTAINER, class DISTANCE_CONTAINER, class POINT_CONTAINER>
01192 unsigned int GetInSphere
01193 (
01194 OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor,
01195 OBJECT_MARKER & marker,
01196 const CoordType & sphere_center,
01197 const ScalarType & sphere_radius,
01198 OBJECT_POINTER_CONTAINER & objects,
01199 DISTANCE_CONTAINER & distances,
01200 POINT_CONTAINER & points,
01201 bool sort_per_distance = true,
01202 bool allow_zero_distance = true
01203 )
01204 {
01205 BoundingBoxType query_bb(sphere_center, sphere_radius);
01206 vcg::Box3i integer_bb = m_UniformGrid.Interize(query_bb);
01207
01208 vcg::Point3i index;
01209 std::vector< std::vector< ObjectPointer >* > contained_objects;
01210 std::vector< ObjectPointer >* tmp;
01211 for (index.X()=integer_bb.min.X(); index.X()<=integer_bb.max.X(); index.X()++)
01212 for (index.Y()=integer_bb.min.Y(); index.Y()<=integer_bb.max.Y(); index.Y()++)
01213 for (index.Z()=integer_bb.min.Z(); index.Z()<=integer_bb.max.Z(); index.Z()++)
01214 if ((tmp=(*this)[index])!=NULL)
01215 contained_objects.push_back(tmp);
01216
01217 std::vector< Neighbor > results;
01218 for (std::vector< std::vector< ObjectPointer >* >::iterator iVec=contained_objects.begin(), eVec=contained_objects.end(); iVec!=eVec; iVec++)
01219 for (std::vector< ObjectPointer >::iterator iObj=(*iVec)->begin(), eObj=(*iVec)->end(); iObj!=eObj; iObj++ )
01220 {
01221 int r = int(results.size());
01222 results.push_back(Neighbor());
01223 results[r].object = *iObj;
01224 results[r].distance = sphere_radius;
01225 if (!distance_functor(*results[r].object, sphere_center, results[r].distance, results[r].nearest_point) || (results[r].distance==ScalarType(0.0) && !allow_zero_distance) )
01226 results.pop_back();
01227 }
01228
01229 if (sort_per_distance)
01230 std::sort( results.begin(), results.end() );
01231
01232 int number_of_objects = int(results.size());
01233 distances.resize(number_of_objects);
01234 points.resize(number_of_objects);
01235 objects.resize(number_of_objects);
01236 for (int i=0, size=int(results.size()); i<size; i++)
01237 {
01238 distances[i] = results[i].distance;
01239 points[i] = results[i].nearest_point;
01240 objects[i] = results[i].object;
01241 }
01242 return number_of_objects;
01243 }
01244
01245
01253 std::vector< ObjectPointer >* operator[](const CoordinateType &query)
01254 {
01255 typename UniformGrid::CellCoordinate ug_index = m_UniformGrid.Interize(query);
01256 if (!m_Bitmap[ug_index])
01257 return NULL;
01258
01259 HashTable::EntryCoordinate ht_index = PerfectHashFunction(ug_index);
01260 std::vector< ObjectPointer >* result = m_HashTable[ht_index];
01261 return result;
01262 }
01263
01264
01265 std::vector< ObjectPointer >* operator[](const typename UniformGrid::CellCoordinate &query)
01266 {
01267 if(!m_Bitmap[query])
01268 return NULL;
01269
01270 HashTable::EntryCoordinate ht_index = PerfectHashFunction(query);
01271 std::vector< ObjectPointer >* result = m_HashTable[ht_index]->domain_data;
01272 return result;
01273 }
01274
01275
01276 protected:
01282 typename HashTable::EntryCoordinate PerfectHashFunction(const typename UniformGrid::CellCoordinate &query)
01283 {
01284 typename HashTable::EntryCoordinate result;
01285 OffsetTable::OffsetPointer offset = m_OffsetTable[ m_OffsetTable.DomainToOffsetTable(query) ];
01286 result = m_HashTable.DomainToHashTable( Shift(query, *offset) );
01287 return result;
01288 }
01289
01290
01297 typename HashTable::EntryCoordinate Shift(const vcg::Point3i &entry, const typename OffsetTable::Offset &offset)
01298 {
01299 HashTable::EntryCoordinate result;
01300 result.X() = entry.X() + int(offset.X());
01301 result.Y() = entry.Y() + int(offset.Y());
01302 result.Z() = entry.Z() + int(offset.Z());
01303 return result;
01304 }
01305
01306
01313 void Finalize()
01314 {
01315 #ifdef _DEBUG
01316 for (UniformGrid::EntryIterator iUGEntry=m_UniformGrid.Begin(), eUGEntry=m_UniformGrid.End(); iUGEntry!=eUGEntry; iUGEntry++)
01317 assert(m_Bitmap.ContainsData(iUGEntry.GetPosition())==((*iUGEntry)->size()>0));
01318 #endif
01319 m_HashTable.Finalize();
01320 m_UniformGrid.Finalize();
01321 m_OffsetTable.Finalize();
01322 }
01323
01324
01331 bool IsAValidOffset(const std::vector< typename UniformGrid::CellCoordinate > *pre_image, const typename OffsetTable::Offset &offset)
01332 {
01333 int ht_size = m_HashTable.GetSize();
01334 int sqr_ht_size = ht_size*ht_size;
01335 std::vector< int > involved_entries;
01336 for (int i=0, pre_image_size=int((*pre_image).size()); i<pre_image_size; i++)
01337 {
01338 typename UniformGrid::CellCoordinate domain_entry = (*pre_image)[i];
01339 typename HashTable::EntryCoordinate hash_entry = m_HashTable.DomainToHashTable( Shift(domain_entry, offset) );
01340 if (!m_HashTable.IsFree(hash_entry))
01341 return false;
01342 else
01343 involved_entries.push_back(hash_entry.X()*sqr_ht_size + hash_entry.Y()*ht_size + hash_entry.Z());
01344 }
01345
01346
01347 std::sort(involved_entries.begin(), involved_entries.end());
01348 for (int i=0, j=1, size=int(involved_entries.size()); j<size; i++, j++)
01349 if (involved_entries[i]==involved_entries[j])
01350 return false;
01351
01352 return true;
01353 }
01354
01355
01363 int GetUnefectiveOffsetTableSize(const int hash_table_size, const int offset_table_size)
01364 {
01365 int result = offset_table_size;
01366 while (GreatestCommonDivisor(hash_table_size, result)!=1 || hash_table_size%result==0)
01367 result += (hash_table_size%2==0) ? (result%2)+1 : 1;
01368 return result;
01369 }
01370
01371
01377 void PerformFastConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
01378 {
01379 int offset_table_size = (int) ceilf(powf(m_SIGMA*float(number_of_filled_cells), 1.0f/float(m_DIMENSION)));
01380 int hash_table_size = m_HashTable.GetSize();
01381 int failed_construction_count = 0;
01382 do
01383 {
01384 offset_table_size += failed_construction_count++;
01385 offset_table_size = GetUnefectiveOffsetTableSize(hash_table_size, offset_table_size);
01386 }
01387 while(!OffsetTableConstructionSucceded(offset_table_size, callback));
01388 }
01389
01390
01396 void PerformCompactConstruction(const int number_of_filled_cells, vcg::CallBackPos *callback)
01397 {
01398 int min_successfully_dimension = std::numeric_limits<int>::max();
01399 int hash_table_size = m_HashTable.GetSize();
01400 int half_hash_table_size = int(float(hash_table_size)/2.0f);
01401
01402
01403 for (int t=0; t<m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION; t++)
01404 {
01405 int lower_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(double(rand())/double(RAND_MAX)*half_hash_table_size + 1) );
01406 int upper_bound = GetUnefectiveOffsetTableSize(hash_table_size, int(((double) rand() / (double) RAND_MAX) * hash_table_size + half_hash_table_size));
01407
01408
01409
01410 int candidate_offset_table_size;
01411 int last_tried_size = -1;
01412 while (lower_bound<upper_bound)
01413 {
01414 candidate_offset_table_size = GetUnefectiveOffsetTableSize(hash_table_size, int(floorf((lower_bound+upper_bound)/2.0f)));
01415
01416
01417
01418 if (last_tried_size==candidate_offset_table_size)
01419 break;
01420
01421 if ( OffsetTableConstructionSucceded((last_tried_size=candidate_offset_table_size), callback) )
01422 {
01423 upper_bound = candidate_offset_table_size;
01424 min_successfully_dimension = vcg::math::Min(candidate_offset_table_size, min_successfully_dimension);
01425 }
01426 else
01427 lower_bound = candidate_offset_table_size;
01428
01429 m_HashTable.Clear();
01430 m_HashTable.BuildFreeEntryList();
01431 m_OffsetTable.Clear();
01432 m_Bitmap.Clear();
01433 }
01434 #ifdef _DEBUD
01435 printf("\nPerfectSpatialHashing: minimum offset table size found at the %d-th iteration was %d\n", (t+1), min_successfully_dimension);
01436 #endif
01437 }
01438
01439
01440 while (!OffsetTableConstructionSucceded(min_successfully_dimension, callback))
01441 {
01442 m_HashTable.Clear();
01443 m_HashTable.BuildFreeEntryList();
01444 m_OffsetTable.Clear();
01445 m_Bitmap.Clear();
01446 }
01447 }
01448
01449
01450
01456 bool OffsetTableConstructionSucceded(const int offset_table_size, vcg::CallBackPos *callback)
01457 {
01458 m_OffsetTable.Allocate(offset_table_size);
01459 m_OffsetTable.BuildH1PreImage(m_UniformGrid.Begin(), m_UniformGrid.End());
01460
01461 std::list< OffsetTable::PreImage > preimage_slots;
01462 m_OffsetTable.GetPreImageSortedPerCardinality(preimage_slots);
01463
01464 char msg[128];
01465 sprintf(msg, "Building offset table of resolution %d", m_OffsetTable.GetSize());
01466 int step = int(preimage_slots.size())/100;
01467 int number_of_slots = int(preimage_slots.size());
01468 int perc = 0;
01469 int iter = 0;
01470 for (std::list< OffsetTable::PreImage >::iterator iPreImage=preimage_slots.begin(), ePreImage=preimage_slots.end(); iPreImage!=ePreImage; iPreImage++, iter++)
01471 {
01472 if (callback!=NULL && iter%step==0 && (perc=iter*100/number_of_slots)<100) (*callback)(perc, msg);
01473
01474 bool found_valid_offset = false;
01475 typename OffsetTable::Offset candidate_offset;
01476
01477
01478 std::vector< typename OffsetTable::Offset > consistent_offsets;
01479 m_OffsetTable.SuggestConsistentOffsets( (*iPreImage).entry_index, consistent_offsets);
01480 for (std::vector< typename OffsetTable::Offset >::iterator iOffset=consistent_offsets.begin(), eOffset=consistent_offsets.end(); iOffset!=eOffset && !found_valid_offset; iOffset++)
01481 if (IsAValidOffset(iPreImage->pre_image, *iOffset))
01482 {
01483 found_valid_offset = true;
01484 candidate_offset = *iOffset;
01485 }
01486
01487
01488
01489 if (!found_valid_offset)
01490 {
01491 std::vector< typename UniformGrid::CellCoordinate > *pre_image = (*iPreImage).pre_image;
01492 for (std::vector< typename UniformGrid::CellCoordinate >::iterator iPreImage=pre_image->begin(), ePreImage=pre_image->end(); iPreImage!=ePreImage && !found_valid_offset; iPreImage++)
01493 for (NeighboringEntryIterator iUGNeighbourhood=m_UniformGrid.GetNeighboringEntryIterator(*iPreImage); iUGNeighbourhood<6 && !found_valid_offset; iUGNeighbourhood++ )
01494 if (!m_OffsetTable.IsFree( m_OffsetTable.DomainToOffsetTable( *iUGNeighbourhood ) ))
01495 {
01496 HashTable::EntryCoordinate ht_entry = PerfectHashFunction(*iUGNeighbourhood);
01497 for (NeighboringEntryIterator iHTNeighbourhood=m_HashTable.GetNeighborintEntryIterator(ht_entry); iHTNeighbourhood<6 && !found_valid_offset; iHTNeighbourhood++)
01498 if (m_HashTable.IsFree(*iHTNeighbourhood))
01499 {
01500 candidate_offset.Import( *iHTNeighbourhood-m_HashTable.DomainToHashTable(*iPreImage) ) ;
01501
01502 if (IsAValidOffset(pre_image, candidate_offset))
01503 found_valid_offset = true;
01504 }
01505 }
01506 }
01507
01508 if (!found_valid_offset)
01509 {
01510
01511 for (int i=0; i<m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET && !found_valid_offset; i++)
01512 {
01513 HashTable::EntryCoordinate base_entry = (*iPreImage).pre_image->at(0);
01514 do
01515 m_OffsetTable.GetRandomOffset(candidate_offset);
01516 while (!m_HashTable.IsFree( m_HashTable.DomainToHashTable( Shift(base_entry, candidate_offset) ) ));
01517
01518 if (IsAValidOffset( (*iPreImage).pre_image, candidate_offset))
01519 found_valid_offset = true;
01520 }
01521
01522
01523
01524 for (std::list< HashTable::EntryCoordinate >::const_iterator iFreeCell=m_HashTable.GetFreeEntryList()->begin(), eFreeCell=m_HashTable.GetFreeEntryList()->end(); iFreeCell!=eFreeCell && !found_valid_offset; iFreeCell++)
01525 {
01526 UniformGrid::CellCoordinate domain_entry = (*iPreImage).pre_image->at(0);
01527 OffsetTable::EntryCoordinate offset_entry = m_OffsetTable.DomainToOffsetTable(domain_entry);
01528 HashTable::EntryCoordinate hashtable_entry = m_HashTable.DomainToHashTable(domain_entry);
01529 candidate_offset.Import(*iFreeCell - hashtable_entry);
01530
01531 if ( IsAValidOffset(iPreImage->pre_image, candidate_offset) )
01532 found_valid_offset = true;
01533 }
01534 }
01535
01536
01537
01538 if (found_valid_offset)
01539 {
01540 m_OffsetTable.SetOffset( (*iPreImage->pre_image).at(0), candidate_offset);
01541 for (int c=0, pre_image_cardinality = iPreImage->cardinality; c<pre_image_cardinality; c++)
01542 {
01543 HashTable::EntryCoordinate ht_entry = PerfectHashFunction( (*iPreImage->pre_image).at(c));
01544 std::vector< ObjectPointer > *domain_data = m_UniformGrid[ (*iPreImage->pre_image).at(c) ];
01545 m_HashTable.SetEntry(ht_entry, domain_data );
01546 m_Bitmap.SetFlag((*iPreImage->pre_image).at(c));
01547 }
01548 }
01549 else
01550 {
01551 m_OffsetTable.Clear();
01552 m_HashTable.Clear();
01553 m_HashTable.BuildFreeEntryList();
01554 m_Bitmap.Clear();
01555 return false;
01556 }
01557 }
01558
01559 if (callback!=NULL) (*callback)(100, msg);
01560 return true;
01561 }
01562
01563
01564
01565
01566
01567 protected:
01568 UniformGrid m_UniformGrid;
01569 OffsetTable m_OffsetTable;
01570 HashTable m_HashTable;
01571 BinaryImage m_Bitmap;
01572
01573 const static float m_BOUNDING_BOX_EXPANSION_FACTOR;
01574 const static float m_SIGMA;
01575 const static int m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION;
01576 const static int m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET;
01577 const static int m_DIMENSION;
01578 };
01579
01581
01582 }
01583
01584 template < class OBJECT_TYPE, class SCALAR_TYPE >
01585 const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_MAX_NUM_OF_RANDOM_GENERATED_OFFSET = 32;
01586
01587 template < class OBJECT_TYPE, class SCALAR_TYPE >
01588 const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_MAX_TRIALS_IN_COMPACT_CONSTRUCTION = 5;
01589
01590 template < class OBJECT_TYPE, class SCALAR_TYPE >
01591 const int vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_DIMENSION = 3;
01592
01593 template < class OBJECT_TYPE, class SCALAR_TYPE >
01594 const SCALAR_TYPE vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_BOUNDING_BOX_EXPANSION_FACTOR = SCALAR_TYPE(0.035);
01595
01596 template < class OBJECT_TYPE, class SCALAR_TYPE >
01597 const SCALAR_TYPE vcg::PerfectSpatialHashing< OBJECT_TYPE, SCALAR_TYPE >::m_SIGMA = SCALAR_TYPE(1.0f/(2.0f*SCALAR_TYPE(m_DIMENSION)));
01598
01599 #endif //VCG_SPACE_INDEX_PERFECT_SPATIAL_HASHING_H
