1/*
2 Copyright (c) 2005-2019 Intel Corporation
3
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
7
8 http://www.apache.org/licenses/LICENSE-2.0
9
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
15*/
16
17/*
18 This file contains a few implementations, so it may look overly complicated.
19 The most efficient implementation is also separated into convex_hull_sample.cpp
20*/
21#include "convex_hull.h"
22
23typedef util::point<double> point_t;
24
25#ifndef USETBB
26 #define USETBB 1
27#endif
28#ifndef USECONCVEC
29 #define USECONCVEC 1
30#endif
31
32#if !USETBB // Serial implementation of Quick Hull algorithm
33
34typedef std::vector< point_t > pointVec_t;
35
36void serial_initialize(pointVec_t &points);
37
38// C++ style serial code
39
40class FindXExtremum : public std::unary_function<const point_t&, void> {
41public:
42 typedef enum {
43 minX, maxX
44 } extremumType;
45
46 FindXExtremum(const point_t& frstPoint, extremumType exType_)
47 : extrXPoint(frstPoint), exType(exType_) {}
48
49 void operator()(const point_t& p) {
50 if(closerToExtremum(p))
51 extrXPoint = p;
52 }
53
54 operator point_t () {
55 return extrXPoint;
56 }
57
58private:
59 const extremumType exType;
60 point_t extrXPoint;
61
62 bool closerToExtremum(const point_t &p) const {
63 switch(exType) {
64 case minX:
65 return p.x<extrXPoint.x; break;
66 case maxX:
67 return p.x>extrXPoint.x; break;
68 }
69 return false; // avoid warning
70 }
71};
72
73template <FindXExtremum::extremumType type>
74point_t extremum(const pointVec_t &points) {
75 assert(!points.empty());
76 return std::for_each(points.begin(), points.end(), FindXExtremum(points[0], type));
77}
78
79class SplitByCP : public std::unary_function<const point_t&, void> {
80 pointVec_t &reducedSet;
81 point_t p1, p2;
82 point_t farPoint;
83 double howFar;
84public:
85
86 SplitByCP( point_t _p1, point_t _p2, pointVec_t &_reducedSet)
87 : p1(_p1), p2(_p2), reducedSet(_reducedSet), howFar(0), farPoint(p1) {}
88
89 void operator()(const point_t& p) {
90 double cp;
91 if( (p != p1) && (p != p2) ) {
92 cp = util::cross_product(p1, p2, p);
93 if(cp>0) {
94 reducedSet.push_back(p);
95 if(cp>howFar) {
96 farPoint = p;
97 howFar = cp;
98 }
99 }
100 }
101 }
102
103 operator point_t (){
104 return farPoint;
105 }
106};
107
108point_t divide(const pointVec_t &P, pointVec_t &P_reduced, const point_t &p1, const point_t &p2) {
109 SplitByCP splitByCP(p1, p2, P_reduced);
110 point_t farPoint = std::for_each(P.begin(), P.end(), splitByCP);
111
112 if(util::verbose) {
113 std::stringstream ss;
114 ss << P.size() << " nodes in bucket"<< ", "
115 << "dividing by: [ " << p1 << ", " << p2 << " ], "
116 << "farthest node: " << farPoint;
117 util::OUTPUT.push_back(ss.str());
118 }
119
120 return farPoint;
121}
122
123void divide_and_conquer(const pointVec_t &P, pointVec_t &H, point_t p1, point_t p2) {
124 assert(P.size() >= 2);
125 pointVec_t P_reduced;
126 pointVec_t H1, H2;
127 point_t p_far = divide(P, P_reduced, p1, p2);
128 if (P_reduced.size()<2) {
129 H.push_back(p1);
130 H.insert(H.end(), P_reduced.begin(), P_reduced.end());
131 }
132 else {
133 divide_and_conquer(P_reduced, H1, p1, p_far);
134 divide_and_conquer(P_reduced, H2, p_far, p2);
135
136 H.insert(H.end(), H1.begin(), H1.end());
137 H.insert(H.end(), H2.begin(), H2.end());
138 }
139}
140
141void quickhull(const pointVec_t &points, pointVec_t &hull) {
142 if (points.size() < 2) {
143 hull.insert(hull.end(), points.begin(), points.end());
144 return;
145 }
146 point_t p_maxx = extremum<FindXExtremum::maxX>(points);
147 point_t p_minx = extremum<FindXExtremum::minX>(points);
148
149 pointVec_t H;
150
151 divide_and_conquer(points, hull, p_maxx, p_minx);
152 divide_and_conquer(points, H, p_minx, p_maxx);
153 hull.insert(hull.end(), H.begin(), H.end());
154}
155
156
157int main(int argc, char* argv[]) {
158 util::ParseInputArgs(argc, argv);
159
160 pointVec_t points;
161 pointVec_t hull;
162 util::my_time_t tm_init, tm_start, tm_end;
163
164 std::cout << "Starting serial version of QUICK HULL algorithm" << std::endl;
165
166 tm_init = util::gettime();
167 serial_initialize(points);
168 tm_start = util::gettime();
169 std::cout << "Init time: " << util::time_diff(tm_init, tm_start) << " Points in input: " << points.size() << "\n";
170 tm_start = util::gettime();
171 quickhull(points, hull);
172 tm_end = util::gettime();
173 std::cout << "Serial time: " << util::time_diff(tm_start, tm_end) << " Points in hull: " << hull.size() << "\n";
174}
175
176#else // USETBB - parallel version of Quick Hull algorithm
177
178#include "tbb/task_scheduler_init.h"
179#include "tbb/parallel_for.h"
180#include "tbb/parallel_reduce.h"
181#include "tbb/blocked_range.h"
182
183typedef tbb::blocked_range<size_t> range_t;
184
185#if USECONCVEC
186#include "tbb/concurrent_vector.h"
187
188typedef tbb::concurrent_vector<point_t> pointVec_t;
189
190void appendVector(const point_t* src, size_t srcSize, pointVec_t& dest) {
191 std::copy(src, src + srcSize, dest.grow_by(srcSize));
192}
193
194void appendVector(const pointVec_t& src, pointVec_t& dest) {
195 std::copy(src.begin(), src.end(), dest.grow_by(src.size()));
196}
197
198void grow_vector_to_at_least(pointVec_t& vect, size_t size) {
199 vect.grow_to_at_least(size);
200}
201#else // USE STD::VECTOR - include spin_mutex.h and lock vector operations
202#include "tbb/spin_mutex.h"
203
204typedef tbb::spin_mutex mutex_t;
205typedef std::vector<point_t> pointVec_t;
206
207void appendVector(mutex_t& insertMutex, const pointVec_t& src, pointVec_t& dest) {
208 mutex_t::scoped_lock lock(insertMutex);
209 dest.insert(dest.end(), src.begin(), src.end());
210}
211
212void appendVector(mutex_t& insertMutex, const point_t* src, size_t srcSize,
213 pointVec_t& dest) {
214 mutex_t::scoped_lock lock(insertMutex);
215 dest.insert(dest.end(), src, src + srcSize);
216}
217
218void grow_vector_to_at_least(mutex_t& mutex, pointVec_t& vect, size_t size) {
219 mutex_t::scoped_lock lock(mutex);
220 if (vect.size()< size){
221 vect.resize(size);
222 }
223}
224#endif // USECONCVEC
225
226class FillRNDPointsVector {
227 pointVec_t &points;
228public:
229 static const size_t grainSize = cfg::generateGrainSize;
230#if !USECONCVEC
231 static mutex_t pushBackMutex;
232#endif // USECONCVEC
233
234 explicit FillRNDPointsVector(pointVec_t& _points)
235 : points(_points){}
236
237 void operator()(const range_t& range) const {
238 util::rng the_rng(range.begin());
239 const size_t i_end = range.end();
240 size_t count = 0;
241#if USECONCVEC
242 points.grow_to_at_least(i_end);
243#else // Locked enlarge to a not thread-safe STD::VECTOR
244 grow_vector_to_at_least(pushBackMutex,points,i_end);
245#endif // USECONCVEC
246
247 for(size_t i = range.begin(); i != i_end; ++i) {
248 points[i]=util::GenerateRNDPoint<double>(count,the_rng,util::rng::max_rand);
249 }
250 }
251};
252
253class FillRNDPointsVector_buf {
254 pointVec_t &points;
255public:
256 static const size_t grainSize = cfg::generateGrainSize;
257#if !USECONCVEC
258 static mutex_t insertMutex;
259#endif // USECONCVEC
260
261 explicit FillRNDPointsVector_buf(pointVec_t& _points)
262 : points(_points){}
263
264 void operator()(const range_t& range) const {
265 util::rng the_rng(range.begin());
266 const size_t i_end = range.end();
267 size_t count = 0, j = 0;
268 point_t tmp_vec[grainSize];
269
270 for(size_t i=range.begin(); i!=i_end; ++i) {
271 tmp_vec[j++] = util::GenerateRNDPoint<double>(count,the_rng,util::rng::max_rand);
272 }
273#if USECONCVEC
274 grow_vector_to_at_least(points,range.end());
275#else // USE STD::VECTOR
276 grow_vector_to_at_least(insertMutex,points,range.end());
277#endif // USECONCVEC
278 std::copy(tmp_vec, tmp_vec+j,points.begin()+range.begin());
279 }
280};
281
282#if !USECONCVEC
283mutex_t FillRNDPointsVector::pushBackMutex = mutex_t();
284mutex_t FillRNDPointsVector_buf::insertMutex = mutex_t();
285#endif
286
287template<typename BodyType>
288void initialize(pointVec_t &points) {
289 //This function generate the same series of point on every call.
290 //Reproducibility is needed for benchmarking to produce reliable results.
291 //It is achieved through the following points:
292 // - FillRNDPointsVector_buf instance has its own local instance
293 // of random number generator, which in turn does not use any global data
294 // - tbb::simple_partitioner produce the same set of ranges on every call to
295 // tbb::parallel_for
296 // - local RNG instances are seeded by the starting indexes of corresponding ranges
297 // - grow_to_at_least() enables putting points into the resulting vector in deterministic order
298 // (unlike concurrent push_back or grow_by).
299
300 // In the buffered version, a temporary storage for as much as grainSize elements
301 // is allocated inside the body. Since auto_partitioner may increase effective
302 // range size which would cause a crash, simple partitioner has to be used.
303
304 tbb::parallel_for(range_t(0, cfg::numberOfPoints, BodyType::grainSize),
305 BodyType(points), tbb::simple_partitioner());
306}
307
308class FindXExtremum {
309public:
310 typedef enum {
311 minX, maxX
312 } extremumType;
313
314 static const size_t grainSize = cfg::findExtremumGrainSize;
315
316 FindXExtremum(const pointVec_t& points_, extremumType exType_)
317 : points(points_), exType(exType_), extrXPoint(points[0]) {}
318
319 FindXExtremum(const FindXExtremum& fxex, tbb::split)
320 : points(fxex.points), exType(fxex.exType), extrXPoint(fxex.extrXPoint) {}
321
322 void operator()(const range_t& range) {
323 const size_t i_end = range.end();
324 if(!range.empty()) {
325 for(size_t i = range.begin(); i != i_end; ++i) {
326 if(closerToExtremum(points[i])) {
327 extrXPoint = points[i];
328 }
329 }
330 }
331 }
332
333 void join(const FindXExtremum &rhs) {
334 if(closerToExtremum(rhs.extrXPoint)) {
335 extrXPoint = rhs.extrXPoint;
336 }
337 }
338
339 point_t extremeXPoint() {
340 return extrXPoint;
341 }
342
343private:
344 const pointVec_t &points;
345 const extremumType exType;
346 point_t extrXPoint;
347 bool closerToExtremum(const point_t &p) const {
348 switch(exType) {
349 case minX:
350 return p.x<extrXPoint.x; break;
351 case maxX:
352 return p.x>extrXPoint.x; break;
353 }
354 return false; // avoid warning
355 }
356};
357
358template <FindXExtremum::extremumType type>
359point_t extremum(const pointVec_t &P) {
360 FindXExtremum fxBody(P, type);
361 tbb::parallel_reduce(range_t(0, P.size(), FindXExtremum::grainSize), fxBody);
362 return fxBody.extremeXPoint();
363}
364
365class SplitByCP {
366 const pointVec_t &initialSet;
367 pointVec_t &reducedSet;
368 point_t p1, p2;
369 point_t farPoint;
370 double howFar;
371public:
372 static const size_t grainSize = cfg::divideGrainSize;
373#if !USECONCVEC
374 static mutex_t pushBackMutex;
375#endif // USECONCVEC
376
377 SplitByCP( point_t _p1, point_t _p2,
378 const pointVec_t &_initialSet, pointVec_t &_reducedSet)
379 : p1(_p1), p2(_p2),
380 initialSet(_initialSet), reducedSet(_reducedSet),
381 howFar(0), farPoint(p1) {
382 }
383
384 SplitByCP( SplitByCP& sbcp, tbb::split )
385 : p1(sbcp.p1), p2(sbcp.p2),
386 initialSet(sbcp.initialSet), reducedSet(sbcp.reducedSet),
387 howFar(0), farPoint(p1) {}
388
389 void operator()( const range_t& range ) {
390 const size_t i_end = range.end();
391 double cp;
392 for(size_t i=range.begin(); i!=i_end; ++i) {
393 if( (initialSet[i] != p1) && (initialSet[i] != p2) ) {
394 cp = util::cross_product(p1, p2, initialSet[i]);
395 if(cp>0) {
396#if USECONCVEC
397 reducedSet.push_back(initialSet[i]);
398#else // Locked push_back to a not thread-safe STD::VECTOR
399 {
400 mutex_t::scoped_lock lock(pushBackMutex);
401 reducedSet.push_back(initialSet[i]);
402 }
403#endif // USECONCVEC
404 if(cp>howFar) {
405 farPoint = initialSet[i];
406 howFar = cp;
407 }
408 }
409 }
410 }
411 }
412
413 void join(const SplitByCP& rhs) {
414 if(rhs.howFar>howFar) {
415 howFar = rhs.howFar;
416 farPoint = rhs.farPoint;
417 }
418 }
419
420 point_t farthestPoint() const {
421 return farPoint;
422 }
423};
424
425class SplitByCP_buf {
426 const pointVec_t &initialSet;
427 pointVec_t &reducedSet;
428 point_t p1, p2;
429 point_t farPoint;
430 double howFar;
431public:
432 static const size_t grainSize = cfg::divideGrainSize;
433#if !USECONCVEC
434 static mutex_t insertMutex;
435#endif // USECONCVEC
436
437 SplitByCP_buf( point_t _p1, point_t _p2,
438 const pointVec_t &_initialSet, pointVec_t &_reducedSet)
439 : p1(_p1), p2(_p2),
440 initialSet(_initialSet), reducedSet(_reducedSet),
441 howFar(0), farPoint(p1) {}
442
443 SplitByCP_buf(SplitByCP_buf& sbcp, tbb::split)
444 : p1(sbcp.p1), p2(sbcp.p2),
445 initialSet(sbcp.initialSet), reducedSet(sbcp.reducedSet),
446 howFar(0), farPoint(p1) {}
447
448 void operator()(const range_t& range) {
449 const size_t i_end = range.end();
450 size_t j = 0;
451 double cp;
452 point_t tmp_vec[grainSize];
453 for(size_t i = range.begin(); i != i_end; ++i) {
454 if( (initialSet[i] != p1) && (initialSet[i] != p2) ) {
455 cp = util::cross_product(p1, p2, initialSet[i]);
456 if(cp>0) {
457 tmp_vec[j++] = initialSet[i];
458 if(cp>howFar) {
459 farPoint = initialSet[i];
460 howFar = cp;
461 }
462 }
463 }
464 }
465
466#if USECONCVEC
467 appendVector(tmp_vec, j, reducedSet);
468#else // USE STD::VECTOR
469 appendVector(insertMutex, tmp_vec, j, reducedSet);
470#endif // USECONCVEC
471 }
472
473 void join(const SplitByCP_buf& rhs) {
474 if(rhs.howFar>howFar) {
475 howFar = rhs.howFar;
476 farPoint = rhs.farPoint;
477 }
478 }
479
480 point_t farthestPoint() const {
481 return farPoint;
482 }
483};
484
485#if !USECONCVEC
486mutex_t SplitByCP::pushBackMutex = mutex_t();
487mutex_t SplitByCP_buf::insertMutex = mutex_t();
488#endif
489
490template <typename BodyType>
491point_t divide(const pointVec_t &P, pointVec_t &P_reduced,
492 const point_t &p1, const point_t &p2) {
493 BodyType body(p1, p2, P, P_reduced);
494 // Must use simple_partitioner (see the comment in initialize() above)
495 tbb::parallel_reduce(range_t(0, P.size(), BodyType::grainSize),
496 body, tbb::simple_partitioner() );
497
498 if(util::verbose) {
499 std::stringstream ss;
500 ss << P.size() << " nodes in bucket"<< ", "
501 << "dividing by: [ " << p1 << ", " << p2 << " ], "
502 << "farthest node: " << body.farthestPoint();
503 util::OUTPUT.push_back(ss.str());
504 }
505
506 return body.farthestPoint();
507}
508
509void divide_and_conquer(const pointVec_t &P, pointVec_t &H,
510 point_t p1, point_t p2, bool buffered) {
511 assert(P.size() >= 2);
512 pointVec_t P_reduced;
513 pointVec_t H1, H2;
514 point_t p_far;
515
516 if(buffered) {
517 p_far = divide<SplitByCP_buf>(P, P_reduced, p1, p2);
518 } else {
519 p_far = divide<SplitByCP>(P, P_reduced, p1, p2);
520 }
521
522 if (P_reduced.size()<2) {
523 H.push_back(p1);
524#if USECONCVEC
525 appendVector(P_reduced, H);
526#else // insert into STD::VECTOR
527 H.insert(H.end(), P_reduced.begin(), P_reduced.end());
528#endif
529 }
530 else {
531 divide_and_conquer(P_reduced, H1, p1, p_far, buffered);
532 divide_and_conquer(P_reduced, H2, p_far, p2, buffered);
533
534#if USECONCVEC
535 appendVector(H1, H);
536 appendVector(H2, H);
537#else // insert into STD::VECTOR
538 H.insert(H.end(), H1.begin(), H1.end());
539 H.insert(H.end(), H2.begin(), H2.end());
540#endif
541 }
542}
543
544void quickhull(const pointVec_t &points, pointVec_t &hull, bool buffered) {
545 if (points.size() < 2) {
546#if USECONCVEC
547 appendVector(points, hull);
548#else // STD::VECTOR
549 hull.insert(hull.end(), points.begin(), points.end());
550#endif // USECONCVEC
551 return;
552 }
553
554 point_t p_maxx = extremum<FindXExtremum::maxX>(points);
555 point_t p_minx = extremum<FindXExtremum::minX>(points);
556
557 pointVec_t H;
558
559 divide_and_conquer(points, hull, p_maxx, p_minx, buffered);
560 divide_and_conquer(points, H, p_minx, p_maxx, buffered);
561#if USECONCVEC
562 appendVector(H, hull);
563#else // STD::VECTOR
564 hull.insert(hull.end(), H.begin(), H.end());
565#endif // USECONCVEC
566}
567
568int main(int argc, char* argv[]) {
569 util::ParseInputArgs(argc, argv);
570
571 int nthreads;
572 util::my_time_t tm_init, tm_start, tm_end;
573
574#if USECONCVEC
575 std::cout << "Starting TBB unbuffered push_back version of QUICK HULL algorithm" << std::endl;
576#else
577 std::cout << "Starting STL locked unbuffered push_back version of QUICK HULL algorithm" << std::endl;
578#endif // USECONCVEC
579
580 for(nthreads=cfg::threads.first; nthreads<=cfg::threads.last; nthreads=cfg::threads.step(nthreads)) {
581 pointVec_t points;
582 pointVec_t hull;
583
584 tbb::task_scheduler_init init(nthreads);
585 tm_init = util::gettime();
586 initialize<FillRNDPointsVector>(points);
587 tm_start = util::gettime();
588 std::cout << "Parallel init time on " << nthreads << " threads: " << util::time_diff(tm_init, tm_start) << " Points in input: " << points.size() << "\n";
589
590 tm_start = util::gettime();
591 quickhull(points, hull, false);
592 tm_end = util::gettime();
593 std::cout << "Time on " << nthreads << " threads: " << util::time_diff(tm_start, tm_end) << " Points in hull: " << hull.size() << "\n";
594 }
595
596#if USECONCVEC
597 std::cout << "Starting TBB buffered version of QUICK HULL algorithm" << std::endl;
598#else
599 std::cout << "Starting STL locked buffered version of QUICK HULL algorithm" << std::endl;
600#endif
601
602 for(nthreads=cfg::threads.first; nthreads<=cfg::threads.last; nthreads=cfg::threads.step(nthreads)) {
603 pointVec_t points;
604 pointVec_t hull;
605
606 tbb::task_scheduler_init init(nthreads);
607
608 tm_init = util::gettime();
609 initialize<FillRNDPointsVector_buf>(points);
610 tm_start = util::gettime();
611 std::cout << "Init time on " << nthreads << " threads: " << util::time_diff(tm_init, tm_start) << " Points in input: " << points.size() << "\n";
612
613 tm_start = util::gettime();
614 quickhull(points, hull, true);
615 tm_end = util::gettime();
616 std::cout << "Time on " << nthreads << " threads: " << util::time_diff(tm_start, tm_end) << " Points in hull: " << hull.size() << "\n";
617 }
618
619 return 0;
620}
621
622#endif // USETBB
623
624void serial_initialize(pointVec_t &points) {
625 points.reserve(cfg::numberOfPoints);
626
627 unsigned int rseed=1;
628 for(size_t i=0, count=0; long(i)<cfg::numberOfPoints; ++i) {
629 points.push_back(util::GenerateRNDPoint<double>(count,&std::rand,RAND_MAX ));
630 }
631}
632