GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <geogram/mesh/boxes_intersections.h>
41			#include <geogram/basic/process.h>
42			#include <geogram/basic/algorithm.h>
43			#include <random>
44
45			/*
46			* The algorithm implemented here for fast vector<Box> x vector<Box> intersection
47			* in the function hybrid() is described in:
48			* Fast software for box intersections
49			* Afra Zoromodian and Herbert Edelsbrunner
50			* International Journal of Computational Geometry & Applications
51			* 2002
52			*
53			* Internally is uses a fast median approximation implemented in the
54			* approximate_median() function, and described in:
55			* Approximating center points with iterative radon points
56			* K. L. Clarkson, D. Eppstein, G. L. Miller, C. Sturtivant, and S- H. Teng.
57			* International Journal of Computational Geometry & Applications
58			* 6 (1996) 357–377
59			*/
60
61
62			namespace {
63			using namespace GEO;
64
65			/**
66			* \brief random number generator.
67			*/
68			#ifdef GARGANTUA
69			typedef std::mt19937_64 RNG;
70			#else
71			typedef std::mt19937 RNG;
72			#endif
73
74			/**
75			* \brief comparison between two coordinates with simulation-of-simplicity
76			* symbolic perturbation.
77			* \details used to uniquely assign an interval or a point to a given
78			* subrange when splitting ranges of intervals or points. Also used for
79			* sorting boxes.
80			* \param[in] a , b the two coordinates to be compared
81			* \param[in] ida , idb two ids or indices corresponding to \p a and \p b
82			* \retval true if \p a is smaller than \p b. In case of equality, the ids
83			* \p ida and \p idb are used to disambiguate
84			* \retval false otherwise
85			*/
86		✗	inline bool lt_sos(double a, double b, index_t ida, index_t idb) {
87		✗	return (a < b || (a == b && ida < idb));
88			}
89
90			/**
91			* \brief Comparator class for sorting boxes along various dimensions
92			*/
93			struct BoxesCompare {
94			/**
95			* \brief Compares two boxes
96			* \param[in] ida , idb the indices of the two boxes to be compared
97			* \retval true if box \p ida is before \p idb
98			* \retval false otherwise
99			* \details uses lt_sos()
100			*/
101		✗	bool operator() (index_t ida, index_t idb) const {
102		✗	double a = boxes[ida].xyz_min[d];
103		✗	double b = boxes[idb].xyz_min[d];
104		✗	return lt_sos(a,b,ida,idb);
105			}
106
107			const Box3d* boxes;
108			index_t d;
109			};
110
111			/************************************************************************/
112
113			/**
114			* \brief a range of boxes used to implement box-box intersection algorithm
115			* \details Stores a pointer to an array of boxes and a begin-end
116			* range of pointers to box indices
117			*/
118			struct BoxesRange {
119			/**
120			* \brief Tests whether a BoxRange is empty
121			* \retval true if this BoxRange contains no box
122			* \retval false otherwise
123			*/
124		✗	bool empty() const {
125		✗	return e == b;
126			}
127
128			/**
129			* \brief Gets the size of a BoxRange
130			* \return the number of boxes in this BoxRange
131			*/
132		✗	index_t size() const {
133		✗	return index_t(e-b);
134			}
135
136			/**
137			* \brief Splits this range into two subranges
138			* \details first subrange contains boxes e for which
139			* (precond && predicate(e)) evaluates as true, and second
140			* subrange contains boxes e for which (precond && predicate(e))
141			* evaluates as false
142			* \param[in] predicate a function taking the index of a box and
143			* returning true if box should be moved to first subrange and false
144			* otherwise
145			* \param[in] precond a precondition anded with the predicate
146			* \return a pair of ranges
147			*/
148		✗	std::pair<BoxesRange, BoxesRange> split(
149			std::function<bool(index_t)> predicate, bool precond = true
150			) {
151		✗	if(!precond) {
152		✗	return std::make_pair(
153		✗	BoxesRange{boxes,b,b},
154		✗	BoxesRange{boxes,b,e}
155		✗	);
156			}
157		✗	index_t* m = std::partition(b,e,predicate);
158			// a BoxesRange is just three pointers, so there is no harm returning
159			// two of them by value.
160		✗	return std::make_pair(
161		✗	BoxesRange{boxes,b,m},
162		✗	BoxesRange{boxes,m,e}
163		✗	);
164			}
165
166
167			/**
168			* \brief gets the lower bound of a box for a given coordinate
169			* \param[in] i the index of the box
170			* \param[in] d the coordinate
171			* \return the coordinate \p d of the lower bound of box \p i
172			*/
173		✗	double xmin(index_t i, index_t d) const {
174		✗	geo_debug_assert(d < 3);
175		✗	return boxes[i].xyz_min[d];
176			}
177
178			/**
179			* \brief gets the lower bound of a box for a given coordinate
180			* \param[in] i the index of the box
181			* \param[in] d the coordinate
182			* \return the coordinate \p d of the lower bound of box \p i
183			*/
184		✗	double xmax(index_t i, index_t d) const {
185		✗	geo_debug_assert(d < 3);
186		✗	return boxes[i].xyz_max[d];
187			}
188
189			/**
190			* \brief gets the lower bound of a box for a given coordinate
191			* \param[in] i a pointer to the index of the box
192			* \param[in] d the coordinate
193			* \return the coordinate \p d of the lower bound of box \p i
194			*/
195		✗	double xmin(index_t* i, index_t d) const {
196		✗	geo_debug_assert(i >= b && i < e);
197		✗	return xmin(*i, d);
198			}
199
200			/**
201			* \brief gets the upper bound of a box for a given coordinate
202			* \param[in] i a pointer to the index of the box
203			* \param[in] d the coordinate
204			* \return the coordinate \p d of the upper bound of box \p i
205			*/
206		✗	double xmax(index_t* i, index_t d) const {
207		✗	geo_debug_assert(d < 3);
208		✗	return xmax(*i,d);
209			}
210
211			/**
212			* \brief Tests whether boxes as interval intersect along all coordinates
213			* from 1 to a given upper bound
214			* \param[in] i pointer to first box index
215			* \param[in] J the BoxesRange in which \p j resides
216			* \param[in] j pointer to second box index
217			* \param[in] d upper bound of coordinates to be tested
218			* \retval true if all coordinate invervals intersect from 1 to d
219			* \retval false otherwise
220			*/
221		✗	bool II_isect_1tod(
222			index_t* i, const BoxesRange& J, index_t* j, index_t d
223			) {
224		✗	if(boxes == J.boxes && *i == *j) { // do not report self-intersection
225		✗	return false;
226			}
227		✗	for(index_t dim=1; dim<=d; ++dim) {
228		✗	if(xmax(i,dim) < J.xmin(j,dim) || xmin(i,dim) > J.xmax(j,dim)) {
229		✗	return false;
230			}
231			}
232		✗	return true;
233			}
234
235			/**
236			* \brief Tests whether a box seen as an interval contains a
237			* (box seen as a) point for a given coordinate
238			* \details does symbolic perturbation on the left bound of the interval,
239			* like in BoxesCompare
240			* \param[in] i pointer to index of first box
241			* \param[in] P the BoxesRange in which \p p resides
242			* \param[in] p pointer to index of second box
243			* \param[in] d coordinate to be tested
244			* \retval true if \p i seen as an interval contains \p p seen as a point
245			* for coord \p d
246			* \retval false otherwise
247			*/
248		✗	bool I_contains_P(
249			index_t* i, const BoxesRange& P, index_t* p, index_t d
250			) {
251		✗	double ixmin = xmin(i,d);
252		✗	double ixmax = xmax(i,d);
253		✗	double px = P.xmin(p,d);
254		✗	return lt_sos(ixmin, px, *i, *p) && (ixmax >= px);
255			}
256
257			const Box3d* boxes; /**< a pointer to the array of 3d boxes */
258			index_t* b; /**< a pointer to the first index */
259			index_t* e; /**< a pointer to one position past the last index */
260			};
261
262			/************************************************************************/
263
264			/**
265			* \brief Reports all pairs (i,p) such that box i contains point p
266			* in a set of boxes I seen as intervals and a set of boxes P seen as
267			* points (stabbing)
268			* \details see Fast software for box intersections,
269			* Afra Zoromodian and Herbert Edelsbrunner,
270			* International Journal of Computational Geometry & Applications, 2002
271			* \param[in] I a set of boxes seen as intervals
272			* \param[in] P a set of boxes seen as points
273			* \param[in] d maximum dimension to be tested
274			* \param[in] report_isect the callback used to report intersections, takes
275			* two integers, i and p
276			* \param[in] swap_ip if set, the parameters i,p of the callback are swapped
277			*/
278		✗	void one_way_scan(
279			BoxesRange I, BoxesRange P, index_t d,
280			std::function<void(index_t,index_t)> report_isect, bool swap_ip = false
281			) {
282		✗	GEO::sort(I.b, I.e, BoxesCompare{I.boxes,0});
283		✗	GEO::sort(P.b, P.e, BoxesCompare{P.boxes,0});
284		✗	for(index_t* i = I.b; i != I.e; ++i) {
285		✗	while(P.b != P.e && lt_sos(P.xmin(P.b,0), I.xmin(i,0), *P.b, *i)) {
286		✗	++P.b;
287			}
288		✗	for(index_t* p=P.b; p!=P.e && P.xmin(p,0) <= I.xmax(i,0); ++p) {
289		✗	if(P.II_isect_1tod(p, I, i, d)) {
290		✗	report_isect(swap_ip ? *p : *i, swap_ip ? *i : *p);
291			}
292			}
293			}
294		✗	}
295
296			/**
297			* \brief Reports all pairs (i,p) such that boxes i and p have intersections
298			* in two sets of boxes I and P
299			* \details reports the same intersections as obtained by calling both
300			* one_way_scan(I,P) and one_way_scan(P,I).
301			* see Fast software for box intersections,
302			* Afra Zoromodian and Herbert Edelsbrunner,
303			* International Journal of Computational Geometry & Applications, 2002
304			* \param[in] I a set of boxes
305			* \param[in] P a set of boxes
306			* \param[in] d maximum dimension to be tested
307			* \param[in] report_isect the callback used to report intersections, takes
308			* two integers, i and p
309			* \param[in] swap_ip if set, the parameters i,p of the callback are swapped
310			*/
311		✗	void modified_two_way_scan(
312			BoxesRange I, BoxesRange P, index_t d,
313			std::function<void(index_t,index_t)> report_isect, bool swap_ip = false
314			) {
315		✗	GEO::sort(I.b, I.e, BoxesCompare{I.boxes,0});
316		✗	GEO::sort(P.b, P.e, BoxesCompare{P.boxes,0});
317		✗	while(I.b != I.e && P.b != P.e) {
318		✗	if(lt_sos(I.xmin(I.b,0), P.xmin(P.b,0), *I.b, *P.b)) {
319		✗	for(index_t* p=P.b; p!=P.e && P.xmin(p,0)<=I.xmax(I.b,0); ++p) {
320		✗	if(P.II_isect_1tod(p,I,I.b,d) && I.I_contains_P(I.b,P,p,d)) {
321		✗	report_isect(swap_ip ? *p : *I.b, swap_ip ? *I.b : *p);
322			}
323			}
324		✗	++I.b;
325			} else {
326		✗	for(index_t* i=I.b; i!=I.e && I.xmin(i,0)<=P.xmax(P.b,0); ++i) {
327		✗	if(I.II_isect_1tod(i,P,P.b,d) && I.I_contains_P(i,P,P.b,d)) {
328		✗	report_isect(swap_ip ? *P.b : *i, swap_ip ? *i : *P.b);
329			}
330			}
331		✗	++P.b;
332			}
333			}
334		✗	}
335
336			/************************************************************************/
337
338			/**
339			* \brief given three boxes seen as points, gets the median one for a
340			* given coordinate
341			* \param[in] boxes pointer to an array of boxes
342			* \param[in] a , b , c pointers to three indices referring to three boxes
343			* in \p boxes
344			* \param[in] dim the dimension along which boxes are sorted
345			* \return a pointer to the index of the median box relative to the order
346			* of their lower points along dimension \p dim
347			*/
348		✗	index_t* median_of_three(
349			const Box3d* boxes, index_t* a, index_t* b, index_t* c, index_t dim
350			) {
351		✗	BoxesCompare C{boxes, dim};
352		✗	if(C(*a,*b)) {
353		✗	if(C(*b,*c)) {
354		✗	return b;
355		✗	} else if(C(*a,*c)) {
356		✗	return c;
357			} else {
358		✗	return a;
359			}
360		✗	} else if(C(*a,*c)) {
361		✗	return a;
362		✗	} else if(C(*b,*c)) {
363		✗	return c;
364			} else {
365		✗	return b;
366			}
367			}
368
369			/************************************************************************/
370
371			/**
372			* \brief computes the approximate median of a set of boxes along a given
373			* dimension
374			* \param[in] boxes pointer to an array of boxes
375			* \param[in] b , e pointers to a sequence of box indices
376			* \param[in] d dimension along with boxes are stored
377			* \param[in] levels number of levels in the approximation of the median
378			* \param[in] rng a random number generator
379			* \details see
380			* Approximating center points with iterative radon points
381			* K. L. Clarkson, D. Eppstein, G. L. Miller, C. Sturtivant
382			* and S- H. Teng.
383			* International Journal of Computational Geometry & Applications
384			* 6 (1996) 357–377
385			*/
386		✗	index_t* approximate_median(
387			const Box3d* boxes, index_t* b, index_t* e, index_t d, int levels,
388			RNG& rng
389			) {
390		✗	if(levels < 0) {
391		✗	auto N = std::distance(b,e);
392		✗	geo_assert(N >= 1);
393		✗	return b + size_t(
394		✗	std::uniform_int_distribution<index_t>(0,index_t(N-1))(rng)
395		✗	);
396			}
397
398		✗	return median_of_three(
399			boxes,
400			approximate_median(boxes, b, e, d, levels-1, rng),
401			approximate_median(boxes, b, e, d, levels-1, rng),
402			approximate_median(boxes, b, e, d, levels-1, rng),
403			d
404		✗	);
405			}
406
407			/***************************************************************************/
408
409			/**
410			* \brief Splits a BoxesRange seen as points along a given dimension
411			* \details this reorders the elements in the input BoxesRange
412			* \param[in] P the BoxesRange seen as points to be split
413			* \param[in] d the dimension along which to split \p P
414			* \param[in] rng a random number generator
415			* \return a triple (P1, x, P2) where P1 and P2 are the constructed
416			* subranges and x the approximate median coordinate.
417			*/
418		✗	std::tuple<BoxesRange, double, BoxesRange> split_points(
419			BoxesRange& P, index_t d, RNG& rng
420			) {
421		✗	index_t N = P.size();
422		✗	int levels = int(0.91 * std::log(double(N)/137.035999206)+1.0);
423		✗	levels = (levels <= 0) ? 1 : levels;
424		✗	index_t* m = approximate_median(P.boxes, P.b, P.e, d, levels, rng);
425		✗	double px_m = P.xmin(m,d);
426			// a BoxesRange is just three pointers, so there is no harm returning
427			// two of them by value.
428		✗	std::pair<BoxesRange, BoxesRange> P1P2 = P.split(
429		✗	[&P, d, px_m](index_t i)->bool {
430		✗	return P.xmin(i,d) < px_m;
431			}
432			);
433		✗	return std::make_tuple(P1P2.first, px_m, P1P2.second);
434			}
435
436			/***************************************************************************/
437
438			struct JobsGroup;
439
440			void hybrid(
441			BoxesRange I, BoxesRange P, index_t d,
442			std::function<void(index_t,index_t)> report_isect,
443			bool swap_ip, double lo, double hi, RNG& rng,
444			JobsGroup* jobs = nullptr
445			);
446
447
448			/**
449			* \brief Stores information related with an invocation of the
450			* hybrid() function that implements the box-box intersection
451			* algorithm.
452			* \details This makes it possible to schedule parallel invocations
453			* of hybrid() during its recursive evaluation.
454			*/
455			struct Job {
456		✗	Job() = default;
457			/**
458			* \brief Default constructor
459			* \details Parameters are the same as hybrid(). Index sequences
460			* pointed at by \p I_in and \p P_in are copied into local vectors.
461			* \see hybrid()
462			*/
463		✗	Job(
464			BoxesRange& I_in, BoxesRange& P_in, index_t d_in,
465			bool swap_ip_in, double lo_in, double hi_in, const RNG& rng_in
466		✗	) :
467		✗	I(I_in), P(P_in), d(d_in),
468		✗	swap_ip(swap_ip_in), lo(lo_in), hi(hi_in),
469		✗	rng(rng_in) {
470		✗	idx.reserve(I.size());
471		✗	pdx.reserve(P.size());
472		✗	idx.insert(idx.end(), I.b, I.e);
473		✗	pdx.insert(pdx.end(), P.b, P.e);
474		✗	I.b = idx.data(); I.e = idx.data() + idx.size();
475		✗	P.b = pdx.data(); P.e = pdx.data() + pdx.size();
476		✗	}
477
478			// TODO: test invoking the job directly and having a lock instead.
479			/**
480			* \brief Runs the jobs
481			* \details Calls the true hybrid() function with the stored parameters
482			*/
483		✗	void run() {
484		✗	::hybrid(
485			I, P, d,
486		✗	[this](index_t a, index_t b) {
487		✗	intersections.emplace_back(a,b);
488		✗	},
489		✗	swap_ip, lo, hi, rng
490			);
491		✗	}
492
493			#ifdef GEO_DEBUG
494			/**
495			* \brief Sanity check
496			* \details Checks that idx and pdx vectors were not moved
497			* in memory between constructor invocation and now.
498			*/
499		✗	void check() {
500		✗	geo_assert(I.b == idx.data());
501		✗	geo_assert(I.e == idx.data() + idx.size());
502		✗	geo_assert(P.b == pdx.data());
503		✗	geo_assert(P.e == pdx.data() + pdx.size());
504		✗	}
505			#endif
506
507			BoxesRange I;
508			BoxesRange P;
509			index_t d;
510			bool swap_ip;
511			double lo;
512			double hi;
513			RNG rng;
514			vector<index_t> idx; /**< local copy of I index range */
515			vector<index_t> pdx; /**< local copy of P index range */
516			vector<std::pair<index_t, index_t>> intersections; /**< local output */
517			};
518
519			/**
520			* \brief A group of Job objects
521			* \details JobGroup has the same interface as hybrid(), but gathers
522			* invocations into a vector of Job objects, and invokes them in
523			* parallel chunk-by-chunk.
524			*/
525			struct JobsGroup {
526			public:
527			/** \brief Maximum number of jobs to be created in parallel mode */
528			static constexpr index_t max_jobs = 128;
529			/** \brief Maximum size of I and P in a job */
530			static constexpr index_t job_cutoff = 131072;
531
532			/**
533			* \brief JobsGroup constructor
534			* \param[in] report_isect_in user callback to report intersections
535			* \param[in] rd_in random device used to initiazize random number
536			* generator in each Job.
537			*/
538		✗	JobsGroup(
539			std::function<void(index_t,index_t)> report_isect_in,
540			std::random_device& rd_in
541		✗	) : report_isect(report_isect_in), rd(rd_in) {
542		✗	jobs.reserve(max_jobs);
543			// To my great surprise, I discovered that reallocating in a
544			// vector of things that have vector members may change the
545			// addresses in the vector members, I did not expect that (or
546			// maybe there is something I did not understand somewhere else,
547			// I should create a minimal example to make sure). It seems that
548			// vector<Job>::resize() calls the copy constructor rather than
549			// the move constructor. Reserving memory in advance avoids this
550			// (unwanted) behavior.
551		✗	}
552
553			/**
554			* \brief JobsGroup destructor
555			* \details Invokes all queued Job objects
556			*/
557		✗	~JobsGroup() {
558		✗	if(jobs.size() != 0) {
559		✗	run_and_flush();
560			}
561		✗	}
562
563			/**
564			* \brief Queues an invocation to hybrid()
565			* \details Parameters are the same as hybrid(). Copies all parameters
566			* and index sequences into a local Job object. If stored vector of
567			* Job object has already of size max_jobs, invokes all stored jobs
568			* in parallel and flushes job queue.
569			* \retval true if I and P were small enough and a job was queued
570			* \retval false otherwise
571			* \see hybrid()
572			*/
573		✗	bool hybrid(
574			BoxesRange& I, BoxesRange& P, index_t d,
575			bool swap_ip, double lo, double hi
576			) {
577		✗	if(I.size() > job_cutoff || P.size() > job_cutoff) {
578		✗	return false;
579			}
580		✗	if(jobs.size() == max_jobs) {
581		✗	run_and_flush();
582			}
583		✗	queue_hybrid(I,P,d,swap_ip,lo,hi);
584		✗	return true;
585			}
586
587			private:
588			/**
589			* \brief Queues a call to hybrid()
590			* \pre queue contains less jobs than max_jobs
591			* \details Parameters are the same as hybrid()
592			* \see hybrid()
593			*/
594		✗	void queue_hybrid(
595			BoxesRange& I, BoxesRange& P, index_t d,
596			bool swap_ip, double lo, double hi
597			) {
598		✗	geo_debug_assert(jobs.size() < max_jobs);
599		✗	jobs.emplace_back(I, P, d, swap_ip, lo, hi, RNG(rd()));
600		✗	}
601
602			/**
603			* \brief Runs all queued jobs in parallel and flushes the queue
604			* \details Uses the callback and random device passed to the constructor
605			*/
606		✗	void run_and_flush() {
607			// Run all job in parallel by calling the real hybrid() function.
608		✗	parallel_for(
609			0,jobs.size(),
610		✗	[this](index_t j) {
611		✗	Job& J = jobs[j];
612			// Make sure nothing weird happened with reallocations (see
613			// remark near jobs.reserve() above).
614			#ifdef GEO_DEBUG
615		✗	J.check();
616			#endif
617		✗	J.run();
618		✗	}
619			);
620
621			// Call user callback for each detected box intersection
622		✗	for(auto& J: jobs) {
623		✗	for(auto& ip: J.intersections) {
624		✗	report_isect(ip.first, ip.second);
625			}
626			}
627
628		✗	jobs.resize(0);
629		✗	}
630
631			vector<Job> jobs;
632			std::function<void(index_t,index_t)> report_isect;
633			std::random_device& rd;
634			};
635
636
637			/***************************************************************************/
638
639			/**
640			* \brief Optimized algorithm that reports all pairs (i,p) such that
641			* boxes i contains (boxes seen as) points p in two sets of
642			* boxes I and P. Much faster than one_way_scan().
643			* \details the recursion is initiated as follows
644			* \code
645			* hybrid(
646			* I,P,2,callback,false,
647			* -std::numeric_limits<double>::max(),
648			* std::numeric_limits<double>::max()
649			* );
650			* \endcode
651			* This is the main algorithm described in:
652			* Fast software for box intersections,
653			* Afra Zoromodian and Herbert Edelsbrunner,
654			* International Journal of Computational Geometry & Applications, 2002
655			* \pre each p in P belongs to [lo,hi) and each i in I intersects [lo,hi)
656			* \param[in] I a set of boxes
657			* \param[in] P a set of boxes
658			* \param[in] d maximum dimension to be tested
659			* \param[in] report_isect the callback used to report intersections, takes
660			* two integers, i and p
661			* \param[in] swap_ip if set, the parameters i,p of the callback are swapped
662			* \param[in] lo , hi range of coordinates
663			* \param[in] rng a random number generator
664			*/
665
666		✗	void hybrid(
667			BoxesRange I, BoxesRange P, index_t d,
668			std::function<void(index_t,index_t)> report_isect,
669			bool swap_ip, double lo, double hi, RNG& rng,
670			JobsGroup* jobs
671			) {
672			static constexpr index_t scanning_cutoff = 1024;
673			static constexpr double inf = -std::numeric_limits<double>::max();
674			static constexpr double sup = std::numeric_limits<double>::max();
675
676			#ifdef GEO_DEBUG // check preconditions
677		✗	for(index_t* p = P.b; p != P.e; ++p) { // Each p belongs to [lo,hi)
678		✗	geo_debug_assert(P.xmin(p,d) >= lo && P.xmin(p,d) < hi);
679			}
680		✗	for(index_t* i = I.b; i != I.e; ++i) { // Each i intersects [lo,hi)
681		✗	geo_debug_assert(I.xmin(i,d) < hi && I.xmax(i,d) >= lo);
682			}
683			#endif
684
685		✗	if( I.empty() || P.empty() || lo >= hi ) {
686		✗	return;
687			}
688
689			// Tentatively send the job to the optional JobsGroup if present
690		✗	if(jobs != nullptr && jobs->hybrid(I,P,d,swap_ip,lo,hi)) {
691		✗	return;
692			}
693
694			// First hybridization: scan instead of third level of segment tree
695		✗	if(d == 0) {
696		✗	one_way_scan(I, P, d, report_isect, swap_ip);
697		✗	return;
698			}
699
700			// Second hybridization: cutoffs to switch to scanning
701		✗	if(I.size() <= scanning_cutoff || P.size() <= scanning_cutoff) {
702		✗	modified_two_way_scan(I, P, d, report_isect, swap_ip);
703		✗	return;
704			}
705
706			// Split I into the parts Ispan that span [lo,hi] and the rest Inonspan
707			// Ispan correspond to the list of segments that would be stored in
708			// current node if using a standard representation of a segment tree.
709			// note: not using structured binding auto[Ispan, Inonspan] because
710			// they are later captured (requires c++20)
711			BoxesRange Ispan; BoxesRange Inonspan;
712		✗	std::tie(Ispan, Inonspan) = I.split(
713		✗	[&I,d,lo,hi](index_t i)->bool{
714		✗	return (I.xmin(i,d) < lo && I.xmax(i,d) > hi);
715			},
716		✗	lo != inf && hi != sup
717		✗	);
718
719			// Two calls for roots of the segment tree at the next level
720		✗	if(!Ispan.empty()) {
721		✗	hybrid(Ispan, P, d-1, report_isect, swap_ip, inf, sup, rng, jobs);
722		✗	hybrid(P, Ispan, d-1, report_isect, !swap_ip, inf, sup, rng, jobs);
723			}
724
725			// divide [lo,hi) into [lo, mi) and [mi, hi)
726			// Pl is the sef of points contained in the left subsegment [lo,mi)
727			// Pr is the sef of points contained in the left subsegment [mi,hi)
728			// note: not using structured binding auto[Pl, mi, Pr] because
729			// they are later captured (requires c++20)
730			BoxesRange Pl; double mi; BoxesRange Pr;
731		✗	std::tie(Pl, mi, Pr) = split_points(P, d, rng);
732
733			// Special case: unable to split points (fallback: modified_two_way_scan)
734		✗	if(Pl.empty() || Pr.empty()) {
735		✗	modified_two_way_scan(Inonspan, P, d, report_isect, swap_ip);
736		✗	return;
737			}
738
739			// Il is the sef of intervals that intersect the left subsegment [lo,mi)
740			// but that do not span the entire segment [lo,hi)
741		✗	BoxesRange Il = Inonspan.split(
742		✗	[&Inonspan, d, mi](index_t i)->bool {
743		✗	return (Inonspan.xmin(i,d) < mi);
744			}
745		✗	).first;
746		✗	hybrid(Il, Pl, d, report_isect, swap_ip, lo, mi, rng, jobs);
747
748			// Ir is the sef of intervals that intersect the right subsegment [mi,hi)
749			// but that do not span the entire segment [lo,hi)
750			// Note that Il and Ir are usually not disjoint.
751		✗	BoxesRange Ir = Inonspan.split(
752		✗	[&Inonspan, d, mi](index_t i)->bool {
753		✗	return (Inonspan.xmax(i,d) >= mi);
754			}
755		✗	).first;
756		✗	hybrid(Ir, Pr, d, report_isect, swap_ip, mi, hi, rng, jobs);
757			}
758
759			/*******************************************************************************/
760
761			#ifdef GEO_ALTERNATIVE_IMPLEMENTATION_KEPT_FOR_REFERENCE
762			/**
763			* \brief Parallel version of boxes_intersections()
764			* \details Used by boxes_intersections() for large datasets
765			* I and P will be subdivided into nI and nP subsets respectively,
766			* then intersections will be computed in parallel on each nI*nP couple
767			* Kept for reference, boxes_intersections_parallel() is faster in general
768			* \see boxes_intersections()
769			*/
770			void boxes_intersections_parallel_split_I_P(
771			const vector<Box3d>& boxes,
772			std::function<void(index_t, index_t)> callback,
773			index_t nI = 6, index_t nP = 6
774			) {
775			vector<index_t> idx(boxes.size());
776			for(index_t i=0; i<boxes.size(); ++i) {
777			idx[i] = i;
778			}
779			vector<index_t> pdx(idx);
780
781			// Create nP copies of I and nI copies of P, so that
782			// each nI*nP thread can manipulate its own copy of I and P
783			vector<index_t> idx_copies;
784			idx_copies.reserve(idx.size()*nP);
785			vector<index_t> pdx_copies;
786			pdx_copies.reserve(pdx.size()*nI);
787			for(index_t p=0; p<nP; ++p) {
788			idx_copies.insert(idx_copies.end(), idx.begin(), idx.end());
789			}
790			for(index_t i=0; i<nI; ++i) {
791			pdx_copies.insert(pdx_copies.end(), pdx.begin(), pdx.end());
792			}
793
794			index_t I_size = index_t(idx.size());
795			index_t I_batch_size = index_t(I_size/nI);
796			index_t P_size = index_t(pdx.size());
797			index_t P_batch_size = index_t(P_size/nP);
798
799			// Job (i,p) I indices are [ I_ptr(i,p) ... I_ptr(i+1,p) )
800			auto I_ptr = [&](index_t i, index_t p)->index_t* {
801			return idx_copies.data() + (
802			p*I_size + ((i==nI) ? I_size : (i * I_batch_size))
803			);
804			};
805
806			// Job (i,p) P indices are [ P_ptr(i,p) ... P_ptr(i,p+1) )
807			auto P_ptr = [&](index_t i, index_t p)->index_t* {
808			return pdx_copies.data() + (
809			i*P_size + ((p==nP) ? P_size : (p * P_batch_size))
810			); // :-) -------^
811			};
812
813			// Initialize jobs (note: not using local storage idx and pdx)
814			std::random_device rd;
815			vector<Job> jobs(nI*nP);
816			for(index_t p=0; p<nP; ++p) {
817			for(index_t i=0; i<nI; ++i) {
818			Job& J = jobs[p*nI+i];
819			J.I=BoxesRange{boxes.data(),I_ptr(i,p),I_ptr(i+1,p)};
820			J.P=BoxesRange{boxes.data(),P_ptr(i,p),P_ptr(i,p+1)};
821			J.d=2;
822			J.swap_ip=false;
823			J.lo=-std::numeric_limits<double>::max();
824			J.hi= std::numeric_limits<double>::max();
825			J.rng = RNG(rd());
826			}
827			}
828
829			// Let's rock and roll !
830			parallel_for(
831			0, nI*nP, [&jobs](index_t j) {
832			jobs[j].run();
833			}
834			);
835
836			// We could also have called callback() asynchronously and handled
837			// concurrency with a lock but it seems to be faster to store
838			// intersections in one vector per thread.
839			for(const auto& job: jobs) {
840			for(const auto& ij: job.intersections) {
841			callback(ij.first, ij.second);
842			}
843			}
844			}
845			#endif
846
847			/**
848			* \brief Reports all pairs (i,p) such that boxes i and p have intersections
849			* in two sets of boxes I and P.
850			* \details Wrapper around the lower-level hybrid() function. Goes parallel
851			* if one of the sets has more than 1024 elements.
852			* \pre each p in P belongs to [lo,hi) and each i in I intersects [lo,hi)
853			* \param[in] I a set of boxes
854			* \param[in] P a set of boxes
855			* \param[in] callback the callback used to report intersections, takes
856			* two integers, i and p
857			*/
858		✗	void hybrid(
859			BoxesRange& I, BoxesRange& P,
860			std::function<void(index_t, index_t)> callback
861			) {
862		✗	std::random_device rd;
863		✗	RNG rng(rd());
864		✗	if(
865		✗	(I.size() >= 1024 || P.size() >= 1024) &&
866		✗	GEO::uses_parallel_algorithm()
867			) {
868			// Parallel mode:
869			// jobs with both I and P smaller than JobsGroup:job_cutoff are
870			// stored in the JobsGroup and later executed in parallel.
871		✗	JobsGroup jobs(callback,rd);
872		✗	hybrid(
873			I,P,2,callback,false,
874		✗	-std::numeric_limits<double>::max(),
875			std::numeric_limits<double>::max(),
876			rng, &jobs
877			);
878		✗	} else {
879		✗	hybrid(
880			I,P,2,callback,false,
881		✗	-std::numeric_limits<double>::max(),
882			std::numeric_limits<double>::max(),
883			rng
884			);
885			}
886		✗	}
887			}
888
889			/******************************************************************************/
890
891			namespace GEO {
892
893		✗	void boxes_intersections(
894			const vector<Box3d>& boxes,
895			std::function<void(index_t, index_t)> callback
896			) {
897		✗	vector<index_t> idx(boxes.size());
898		✗	vector<index_t> pdx(boxes.size());
899		✗	for(index_t i=0; i<boxes.size(); ++i) {
900		✗	idx[i] = i;
901		✗	pdx[i] = i;
902			}
903		✗	BoxesRange I{boxes.data(), idx.data(), idx.data()+idx.size()};
904		✗	BoxesRange P{boxes.data(), pdx.data(), pdx.data()+pdx.size()};
905		✗	hybrid(I, P, callback);
906		✗	}
907
908		✗	void boxes_intersections_hybrid_impl(
909			const Box3d* Iboxes,
910			index_t* Ib,
911			index_t* Ie,
912			const Box3d* Pboxes,
913			index_t* Pb,
914			index_t* Pe,
915			std::function<void(index_t, index_t)> report_intersection
916			) {
917		✗	if((Pb >= Ib && Ib <= Ie) || (Ib >= Pb && Ib <= Pe)) {
918		✗	vector<index_t> pdx;
919		✗	pdx.insert(pdx.end(), Pb, Pe);
920		✗	BoxesRange I{Iboxes, Ib, Ie};
921		✗	BoxesRange P{Pboxes, pdx.data(), pdx.data()+pdx.size()};
922		✗	hybrid(I,P,report_intersection);
923		✗	} else {
924		✗	BoxesRange I{Iboxes, Ib, Ie};
925		✗	BoxesRange P{Pboxes, Pb, Pe};
926		✗	hybrid(I,P,report_intersection);
927			}
928		✗	}
929
930			}
931