VCG Library
trimesh_indexing.cpp
1 #include <iostream>
2 #ifdef _OPENMP
3 #include <omp.h>
4 #endif
5 
6 #include "nanoflann.hpp"
7 
8 #include <vcg/complex/complex.h>
9 
10 #include <wrap/io_trimesh/import.h>
11 #include <wrap/io_trimesh/export.h>
12 
13 
14 #include <vcg/space/index/kdtree/kdtree.h>
15 #include <vcg/space/index/grid_static_ptr.h>
16 #include <vcg/space/index/perfect_spatial_hashing.h>
17 #include <vcg/space/index/spatial_hashing.h>
18 #include <vcg/space/index/octree.h>
19 
20 int num_test = 1000;
21 int kNearest = 256;
22 float queryDist = 0.0037;
23 float bboxratio = 1000.0f;
24 
25 
26 class CVertex;
27 class CFace;
28 class CEdge;
29 
30 class CUsedTypes : public vcg::UsedTypes < vcg::Use< CVertex >::AsVertexType, vcg::Use< CFace >::AsFaceType>{};
31 class CVertex : public vcg::Vertex < CUsedTypes, vcg::vertex::Coord3f, vcg::vertex::Normal3f, vcg::vertex::Radiusf, vcg::vertex::BitFlags, vcg::vertex::Qualityf, vcg::vertex::Color4b>{};
32 class CFace : public vcg::Face < CUsedTypes, vcg::face::VertexRef>{};
33 
34 class CMesh : public vcg::tri::TriMesh < std::vector< CVertex >, std::vector< CFace > > {};
35 
36 int elapsed(int t)
37 {
38  return ((clock()-t)*1000.0)/CLOCKS_PER_SEC;
39 }
40 
41 template <typename T>
42 struct PointCloud
43 {
44  struct Point
45  {
46  T x,y,z;
47  };
48 
49  std::vector<Point> pts;
50 
51  inline size_t kdtree_get_point_count() const { return pts.size(); }
52 
53  inline T kdtree_distance(const T *p1, const size_t idx_p2,size_t size) const
54  {
55  const T d0=p1[0]-pts[idx_p2].x;
56  const T d1=p1[1]-pts[idx_p2].y;
57  const T d2=p1[2]-pts[idx_p2].z;
58  return d0*d0+d1*d1+d2*d2;
59  }
60 
61  inline T kdtree_get_pt(const size_t idx, int dim) const
62  {
63  if (dim==0) return pts[idx].x;
64  else if (dim==1) return pts[idx].y;
65  else return pts[idx].z;
66  }
67 
68  template <class BBOX>
69  bool kdtree_get_bbox(BBOX &bb) const { return false; }
70 
71 };
72 
73 
74 void testKDTree(CMesh& mesh, std::vector<unsigned int>& test_indeces, std::vector<vcg::Point3f>& randomSamples)
75 {
76  std::cout << "==================================================="<< std::endl;
77  std::cout << "KDTree" << std::endl;
78  int t0=clock();
79  // Construction of the kdTree
80  vcg::ConstDataWrapper<CMesh::VertexType::CoordType> wrapperVcg(&mesh.vert[0].P(), mesh.vert.size(), size_t(mesh.vert[1].P().V()) - size_t(mesh.vert[0].P().V()));
81  vcg::KdTree<CMesh::ScalarType> kdTreeVcg(wrapperVcg);
82  std::cout << "Build: " << elapsed(t0) << " ms" << std::endl;
83  int nn=1;
84  // Computation of the point radius
85  float mAveragePointSpacing = 0;
86  t0=clock();
87  #pragma omp parallel for reduction(+: mAveragePointSpacing) schedule(dynamic, 10)
88  for (int i = 0; i < mesh.vert.size(); i++)
89  {
90 #ifdef _OPENMP
91  nn =omp_get_num_threads();
92 #endif
93  vcg::KdTree<CMesh::ScalarType>::PriorityQueue queue;
94  kdTreeVcg.doQueryK(mesh.vert[i].cP(), 16, queue);
95  float newRadius = 2.0f * sqrt(queue.getWeight(0)/ queue.getNofElements());
96  mesh.vert[i].R() -= newRadius;
97  mAveragePointSpacing += newRadius;
98  }
99  std::cout << "Num trhread " << nn << std::endl;
100  mAveragePointSpacing /= mesh.vert.size();
101  std::cout << "Average point radius (OpenMP with" << nn << " threads) " << mAveragePointSpacing << std::endl;
102  std::cout << "Time (OpenMP): " << elapsed(t0) << " ms" << std::endl;
103 
104  queryDist = mAveragePointSpacing * 150;
105 
106  // Test with the radius search
107  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
108  float avgTime = 0.0f;
109  for (int ii = 0; ii < num_test; ii++)
110  {
111  int t0=clock();
112  std::vector<unsigned int> indeces;
113  std::vector<float> dists;
114  kdTreeVcg.doQueryDist(mesh.vert[test_indeces[ii]].cP(), queryDist, indeces, dists);
115  avgTime += elapsed(t0);
116  }
117  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl;
118 
119  // Test with the k-nearest search
120  std::cout << "k-Nearest search (" << num_test*10 << " tests)"<< std::endl;
121  avgTime = 0.0f;
122  for (int ii = 0; ii < num_test * 10; ii++)
123  {
124  int t0=clock();
125  vcg::KdTree<CMesh::ScalarType>::PriorityQueue queue;
126  kdTreeVcg.doQueryK(mesh.vert[test_indeces[ii]].cP(), kNearest, queue);
127  avgTime += elapsed(t0);
128  }
129  std::cout << "Time (k = " << kNearest << "): " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl;
130 
131  // Test with the closest search
132  std::cout << "Closest search (" << num_test*10 << " tests)"<< std::endl;
133  avgTime = 0.0f;
134  for (int ii = 0; ii < num_test * 10; ii++)
135  {
136  int t0=clock();
137  unsigned int index;
138  float minDist;
139  kdTreeVcg.doQueryClosest(randomSamples[ii], index, minDist);
140  avgTime += elapsed(t0);
141  }
142  std::cout << "Time : " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl << std::endl;
143 }
144 
145 
146 void testNanoFLANN(CMesh& mesh, std::vector<unsigned int>& test_indeces, std::vector<vcg::Point3f> randomSamples)
147 {
148  std::cout << "==================================================="<< std::endl;
149  std::cout << "nanoFLANN" << std::endl;
150 
151  PointCloud<float> cloud;
152  cloud.pts.resize(mesh.vert.size());
153  for (size_t i=0; i < mesh.vert.size(); i++)
154  {
155  cloud.pts[i].x = mesh.vert[i].P().X();
156  cloud.pts[i].y = mesh.vert[i].P().Y();
157  cloud.pts[i].z = mesh.vert[i].P().Z();
158  }
159 
160  typedef nanoflann::KDTreeSingleIndexAdaptor<
161  nanoflann::L2_Simple_Adaptor<float, PointCloud<float> > ,
162  PointCloud<float>,
163  3 /* dim */
164  > my_kd_tree_t;
165 
166  // Construction of the nanoFLANN KDtree
167  int t0=clock();
168  my_kd_tree_t index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(16) );
169  index.buildIndex();
170  std::cout << "Build nanoFlann: " << elapsed(t0) << " ms" << std::endl;
171 
172  // Test with the radius search
173  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
174  float avgTime = 0.0f;
175  std::vector<std::pair<size_t,float> > ret_matches;
176  nanoflann::SearchParams params;
177  for (int ii = 0; ii < num_test; ii++)
178  {
179  t0=clock();
180  const size_t nMatches = index.radiusSearch(mesh.vert[test_indeces[ii]].P().V(), queryDist, ret_matches, params);
181  avgTime += elapsed(t0);
182  }
183  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl;
184 
185  // Test with the k-nearest search
186  std::cout << "k-Nearest search (" << num_test*10 << " tests)"<< std::endl;
187  avgTime = 0.0f;
188  std::vector<size_t> ret_index(kNearest);
189  std::vector<float> out_dist_sqr(kNearest);
190  for (int ii = 0; ii < num_test * 10; ii++)
191  {
192  t0=clock();
193  index.knnSearch(mesh.vert[test_indeces[ii]].P().V(), kNearest, &ret_index[0], &out_dist_sqr[0]);
194  avgTime += elapsed(t0);
195  }
196  std::cout << "Time (k = " << kNearest << "): " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl;
197 
198  // Test with the closest search
199  std::cout << "Closest search (" << num_test*10 << " tests)"<< std::endl;
200  avgTime = 0.0f;
201  std::vector<size_t> ret_index_clos(1);
202  std::vector<float> out_dist_sqr_clos(1);
203  for (int ii = 0; ii < num_test * 10; ii++)
204  {
205  t0=clock();
206  index.knnSearch(randomSamples[ii].V(), 1, &ret_index_clos[0], &out_dist_sqr_clos[0]);
207  avgTime += elapsed(t0);
208  }
209  std::cout << "Time : " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl << std::endl;
210 }
211 
212 
213 void testUniformGrid(CMesh& mesh, std::vector<unsigned int>& test_indeces, std::vector<vcg::Point3f>& randomSamples)
214 {
215  std::cout << "==================================================="<< std::endl;
216  std::cout << "Uniform Grid" << std::endl;
217  int t0=clock();
218 
219  // Construction of the uniform grid
220  typedef vcg::GridStaticPtr<CMesh::VertexType, CMesh::VertexType::ScalarType> MeshGrid;
221  MeshGrid uniformGrid;
222  uniformGrid.Set(mesh.vert.begin(), mesh.vert.end());
223  std::cout << "Build: " << elapsed(t0) << " ms" << std::endl;
224 
225  // Test with the radius search
226  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
227  float avgTime = 0.0f;
228  for (int ii = 0; ii < num_test; ii++)
229  {
230  t0=clock();
231  std::vector<CMesh::VertexPointer> vertexPtr;
232  std::vector<CMesh::VertexType::CoordType> points;
233  std::vector<float> dists;
234  vcg::tri::GetInSphereVertex(mesh, uniformGrid, mesh.vert[test_indeces[ii]].cP(), queryDist, vertexPtr, dists, points);
235  avgTime += elapsed(t0);
236  }
237  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl;
238 
239  // Test with the k-nearest search
240  std::cout << "k-Nearest search (" << num_test*10 << " tests)"<< std::endl;
241  avgTime = 0.0f;
242  for (int ii = 0; ii < num_test * 10; ii++)
243  {
244  t0=clock();
245  std::vector<CMesh::VertexPointer> vertexPtr;
246  std::vector<CMesh::VertexType::CoordType> points;
247  std::vector<float> dists;
248  vcg::tri::GetKClosestVertex(mesh, uniformGrid, kNearest, mesh.vert[test_indeces[ii]].cP(), mesh.bbox.Diag(), vertexPtr, dists, points);
249  avgTime += elapsed(t0);
250  }
251  std::cout << "Time (k = " << kNearest << "): " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl;
252 
253  // Test with the Closest search
254  std::cout << "Closest search (" << num_test*10 << " tests)"<< std::endl;
255  avgTime = 0.0f;
256  for (int ii = 0; ii < num_test * 10; ii++)
257  {
258  t0=clock();
259  float minDist;
260  vcg::tri::GetClosestVertex(mesh, uniformGrid, randomSamples[ii], mesh.bbox.Diag(), minDist);
261  avgTime += elapsed(t0);
262  }
263  std::cout << "Time : " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl << std::endl;
264 }
265 
266 
267 
268 void testSpatialHashing(CMesh& mesh, std::vector<unsigned int>& test_indeces, std::vector<vcg::Point3f>& randomSamples)
269 {
270  std::cout << "==================================================="<< std::endl;
271  std::cout << "Spatial Hashing" << std::endl;
272  int t0=clock();
273 
274  // Construction of the uniform grid
275  typedef vcg::SpatialHashTable<CMesh::VertexType, CMesh::VertexType::ScalarType> MeshGrid;
276  MeshGrid uniformGrid;
277  uniformGrid.Set(mesh.vert.begin(), mesh.vert.end());
278  std::cout << "Build: " << elapsed(t0) << " ms" << std::endl;
279 
280  // Test with the radius search
281  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
282  float avgTime = 0.0f;
283  for (int ii = 0; ii < num_test; ii++)
284  {
285  t0=clock();
286  std::vector<CMesh::VertexPointer> vertexPtr;
287  std::vector<CMesh::VertexType::CoordType> points;
288  std::vector<float> dists;
289  vcg::tri::GetInSphereVertex(mesh, uniformGrid, mesh.vert[test_indeces[ii]].cP(), queryDist, vertexPtr, dists, points);
290  avgTime += elapsed(t0);
291  }
292  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl;
293 
294  // Test with the k-nearest search
295  std::cout << "k-Nearest search (" << num_test*10 << " tests)"<< std::endl;
296  avgTime = 0.0f;
297  for (int ii = 0; ii < num_test * 10; ii++)
298  {
299  t0=clock();
300  std::vector<CMesh::VertexPointer> vertexPtr;
301  std::vector<CMesh::VertexType::CoordType> points;
302  std::vector<float> dists;
303  vcg::tri::GetKClosestVertex(mesh, uniformGrid, kNearest, mesh.vert[test_indeces[ii]].cP(), mesh.bbox.Diag(), vertexPtr, dists, points);
304  avgTime += elapsed(t0);
305  }
306  std::cout << "Time (k = " << kNearest << "): " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl;
307 
308  // Test with the Closest search
309  std::cout << "Closest search (" << num_test*10 << " tests)"<< std::endl;
310  avgTime = 0.0f;
311  for (int ii = 0; ii < num_test * 10; ii++)
312  {
313  t0=clock();
314  float minDist;
315  vcg::tri::GetClosestVertex(mesh, uniformGrid, randomSamples[ii], mesh.bbox.Diag(), minDist);
316  avgTime += elapsed(t0);
317  }
318  std::cout << "Time : " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl << std::endl;
319 }
320 
321 
322 
323 void testPerfectSpatialHashing(CMesh& mesh, std::vector<unsigned int>& test_indeces)
324 {
325  std::cout << "==================================================="<< std::endl;
326  std::cout << "Perfect Spatial Hashing" << std::endl;
327  int t0=clock();
328 
329  // Construction of the uniform grid
330  typedef vcg::SpatialHashTable<CMesh::VertexType, CMesh::VertexType::ScalarType> MeshGrid;
331  MeshGrid uniformGrid;
332  uniformGrid.Set(mesh.vert.begin(), mesh.vert.end());
333  std::cout << "Build: " << elapsed(t0) << " ms" << std::endl;
334 
335  // Test with the radius search
336  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
337  float avgTime = 0.0f;
338  for (int ii = 0; ii < num_test; ii++)
339  {
340  t0=clock();
341  std::vector<CMesh::VertexPointer> vertexPtr;
342  std::vector<CMesh::VertexType::CoordType> points;
343  std::vector<float> dists;
344  vcg::tri::GetInSphereVertex(mesh, uniformGrid, mesh.vert[test_indeces[ii]].cP(), queryDist, vertexPtr, dists, points);
345  avgTime += elapsed(t0);
346  }
347  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl << std::endl;
348 }
349 
350 
351 void testOctree(CMesh& mesh, std::vector<unsigned int>& test_indeces, std::vector<vcg::Point3f>& randomSamples)
352 {
353  std::cout << "==================================================="<< std::endl;
354  std::cout << "Octree" << std::endl;
355  int t0=clock();
356 
357  // Construction of the uniform grid
358  typedef vcg::Octree<CMesh::VertexType, CMesh::VertexType::ScalarType> MeshGrid;
359  MeshGrid uniformGrid;
360  uniformGrid.Set(mesh.vert.begin(), mesh.vert.end());
361  std::cout << "Build: " << elapsed(t0) << " ms" << std::endl;
362 
363  // Test with the radius search
364  std::cout << "Radius search (" << num_test << " tests)"<< std::endl;
365  float avgTime = 0.0f;
366  for (int ii = 0; ii < num_test; ii++)
367  {
368  t0=clock();
369  std::vector<CMesh::VertexPointer> vertexPtr;
370  std::vector<CMesh::VertexType::CoordType> points;
371  std::vector<float> dists;
372  vcg::tri::GetInSphereVertex(mesh, uniformGrid, mesh.vert[test_indeces[ii]].cP(), queryDist, vertexPtr, dists, points);
373  avgTime += elapsed(t0);
374  }
375  std::cout << "Time (radius = " << queryDist << "): " << avgTime << " ms (mean " << avgTime / num_test << "ms)" << std::endl;
376 
377  // Test with the k-nearest search
378  std::cout << "k-Nearest search (" << num_test*10 << " tests)"<< std::endl;
379  avgTime = 0.0f;
380  for (int ii = 0; ii < num_test * 10; ii++)
381  {
382  t0=clock();
383  std::vector<CMesh::VertexPointer> vertexPtr;
384  std::vector<CMesh::VertexType::CoordType> points;
385  std::vector<float> dists;
386  vcg::tri::GetKClosestVertex(mesh, uniformGrid, kNearest, mesh.vert[test_indeces[ii]].cP(), mesh.bbox.Diag(), vertexPtr, dists, points);
387  avgTime += elapsed(t0);
388  }
389  std::cout << "Time (k = " << kNearest << "): " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl;
390 
391  // Test with the Closest search
392  std::cout << "Closest search (" << num_test*10 << " tests)"<< std::endl;
393  avgTime = 0.0f;
394  for (int ii = 0; ii < num_test * 10; ii++)
395  {
396  t0=clock();
397  float minDist;
398  vcg::tri::GetClosestVertex(mesh, uniformGrid, randomSamples[ii], mesh.bbox.Diag(), minDist);
399  avgTime += elapsed(t0);
400  }
401  std::cout << "Time : " << avgTime << " ms (mean " << avgTime / (num_test * 10) << "ms)" << std::endl << std::endl;
402 }
403 
404 
405 
406 int main( int argc, char * argv[] )
407 {
408  if (argc < 2) {
409  std::cout << "Invalid arguments" << std::endl;
410  exit(-1);
411  }
412  CMesh mesh;
413  if (vcg::tri::io::Importer<CMesh>::Open(mesh, argv[1]) != 0)
414  std::cout << "Invalid filename" << std::endl;
415 
416  std::cout << "Mesh BBox diagonal: " << mesh.bbox.Diag() << std::endl;
417  std::cout << "Max point random offset: " << mesh.bbox.Diag() / 1000.0f << std::endl << std::endl;
418 
419  vcg::math::MarsenneTwisterRNG randGen;
420  randGen.initialize(0);
421  std::vector<vcg::Point3f> randomSamples;
422  for (int i = 0; i < num_test * 10; i++)
423  randomSamples.push_back(vcg::math::GeneratePointOnUnitSphereUniform<float>(randGen) * randGen.generate01() * mesh.bbox.Diag() / bboxratio);
424 
425  std::vector<unsigned int> test_indeces;
426  for (int i = 0; i < num_test * 10; i++)
427  {
428  int index = randGen.generate01() * (mesh.vert.size() - 1);
429  test_indeces.push_back(index);
430  randomSamples[i] += mesh.vert[i].P();
431  }
432 
433  testKDTree(mesh, test_indeces, randomSamples);
434  testNanoFLANN(mesh, test_indeces, randomSamples);
435  testUniformGrid(mesh, test_indeces, randomSamples);
436  testSpatialHashing(mesh, test_indeces, randomSamples);
437  testPerfectSpatialHashing(mesh, test_indeces);
438  testOctree(mesh, test_indeces, randomSamples);
439 }