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
elapsed(int t)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
kdtree_get_point_countPointCloud51 inline size_t kdtree_get_point_count() const { return pts.size(); }
52
kdtree_distancePointCloud53 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
kdtree_get_ptPointCloud61 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>
kdtree_get_bboxPointCloud69 bool kdtree_get_bbox(BBOX &bb) const { return false; }
70
71 };
72
73
testKDTree(CMesh & mesh,std::vector<unsigned int> & test_indeces,std::vector<vcg::Point3f> & randomSamples)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
testNanoFLANN(CMesh & mesh,std::vector<unsigned int> & test_indeces,std::vector<vcg::Point3f> randomSamples)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
testUniformGrid(CMesh & mesh,std::vector<unsigned int> & test_indeces,std::vector<vcg::Point3f> & randomSamples)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
testSpatialHashing(CMesh & mesh,std::vector<unsigned int> & test_indeces,std::vector<vcg::Point3f> & randomSamples)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
testPerfectSpatialHashing(CMesh & mesh,std::vector<unsigned int> & test_indeces)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
testOctree(CMesh & mesh,std::vector<unsigned int> & test_indeces,std::vector<vcg::Point3f> & randomSamples)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
main(int argc,char * argv[])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 }
440