b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2aabb_intersection.H
1//------------------------------------------------------------------------
2// b2aabb_intersection.H --
3//
4//
5// written by Mathias Doreille
6//
7// Copyright (c) 2009 SMR Engineering & Development SA
8// 2502 Bienne, Switzerland
9//
10// All Rights Reserved. Proprietary source code. The contents of
11// this file may not be disclosed to third parties, copied or
12// duplicated in any form, in whole or in part, without the prior
13// written permission of SMR.
14//------------------------------------------------------------------------
15
16#ifndef B2AABB_INTERSECTION_H_
17#define B2AABB_INTERSECTION_H_
18
19#include <algorithm>
20#include <iostream>
21#include <limits>
22#include <utility>
23#include <vector>
24
25#include "b2ppconfig.h"
26
27namespace b2000 {
28
29/*
30 struct AABox {
31 static const int ndim;
32 double min[NDIM];
33 double max[NDIM];
34 };
35
36 struct PairIntersectionList {
37 void add(AABOX* a, AABOX* b);
38 };
39 */
40
41// compute the box intersection using streamed Multi-level Segment Tree
42// complexity is n^2(log(n))^NDIM
43// storage is n (we newer store the tree)
44//
45// Afra Zomorodian and Herbert Edelsbrunner. Fast software for box intersection. Int. J. Comput.
46// Geom. Appl., 12:143–172, 2002.
47template <typename AABBOX, typename PAIR_INTERSECTION_LIST>
48class AABBIntersection {
49public:
50 static const int min_size_to_go_to_sort = 64;
51 static const int min_size_to_go_to_scanning = 128; // must be not less than 2
52
53 AABBIntersection(
54 AABBOX* box_begin, AABBOX* box_end, PAIR_INTERSECTION_LIST& pair_intersection_list_)
55 : pair_intersection_list(pair_intersection_list_), level(0), dim(0), rand_gen(0) {
56 std::fill_n(current_box.min, AABBOX::ndim, std::numeric_limits<double>::lowest());
57 std::fill_n(current_box.max, AABBOX::ndim, std::numeric_limits<double>::max());
58 current_list_box.reserve(box_end - box_begin);
59 for (; box_begin != box_end; ++box_begin) { current_list_box.push_back(box_begin); }
60 }
61
62 void compute() { rec_compute(current_list_box.begin(), current_list_box.end()); }
63
64private:
65 PAIR_INTERSECTION_LIST& pair_intersection_list;
66 size_t level;
67 int dim;
68
69 struct AABox {
70 double min[AABBOX::ndim];
71 double max[AABBOX::ndim];
72 };
73
74 AABox current_box;
75 typedef std::vector<AABBOX*> current_list_box_t;
76 current_list_box_t current_list_box;
77 typedef std::vector<std::pair<
78 typename current_list_box_t::const_iterator,
79 typename current_list_box_t::const_iterator> >
80 stack_node_box_t;
81 stack_node_box_t stack_node_box;
82 double buffer_pivot[min_size_to_go_to_sort * 2];
83
84 struct NotInBox {
85 NotInBox(const AABox& box_) : box(box_) {}
86 bool operator()(const AABBOX* b) const {
87 for (int d = 0; d != AABBOX::ndim; ++d) {
88 if (b->min[d] > box.min[d] || b->max[d] < box.max[d]) { return true; }
89 }
90 return false;
91 }
92 const AABox& box;
93 };
94
95 struct NotInInterval {
96 NotInInterval(const double min_, const double max_, int dim_)
97 : min(min_), max(max_), dim(dim_) {}
98 bool operator()(const AABBOX* b) const { return b->min[dim] > min || b->max[dim] < max; }
99 double min;
100 double max;
101 int dim;
102 };
103
104 struct Left {
105 Left(double pivot_, int dim_) : pivot(pivot_), dim(dim_) {}
106 bool operator()(const AABBOX* b) const { return b->min[dim] < pivot; }
107 const double pivot;
108 const int dim;
109 };
110
111 struct Right {
112 Right(double pivot_, int dim_) : pivot(pivot_), dim(dim_) {}
113 bool operator()(const AABBOX* b) const { return b->max[dim] <= pivot; }
114 const double pivot;
115 const int dim;
116 };
117
118 struct ScanningSort {
119 bool operator()(const AABBOX* a, const AABBOX* b) const { return a->min[0] < b->min[0]; }
120 };
121
122 size_t rand_gen;
123 bool get_random(size_t s, size_t& i) {
124#if B2_SIZEOF_SIZE_T == 4
125 rand_gen = 1103515245u * rand_gen + 12345u;
126 i = rand_gen % s;
127 return rand_gen & 2147483649u; // 2**31 + 1
128#elif B2_SIZEOF_SIZE_T == 8
129 rand_gen = 2862933555777941757lu * rand_gen + 3037000493lu;
130 i = rand_gen % s;
131 return rand_gen & 9223372036854775809lu; // 2**63 + 1
132#else
133#error "Unimplemented for this sizeof(size_t)"
134#endif
135 }
136
137 double get_approx_median_pivot(
138 const typename current_list_box_t::iterator begin,
139 const typename current_list_box_t::iterator end, int h = 3) {
140 if (--h == 0) {
141 size_t i;
142 if (get_random(end - begin, i)) {
143 const double r = begin[i]->min[dim];
144 if (r < current_box.min[dim]) {
145 return current_box.max[dim];
146 } else {
147 return r;
148 }
149 } else {
150 const double r = begin[i]->max[dim];
151 if (r > current_box.max[dim]) {
152 return current_box.min[dim];
153 } else {
154 return r;
155 }
156 }
157 } else {
158 const double a1 = get_approx_median_pivot(begin, end, h);
159 const double a2 = get_approx_median_pivot(begin, end, h);
160 const double a3 = get_approx_median_pivot(begin, end, h);
161 return a1 < a2 ? (a2 < a3 ? a2 : a3) : (a1 < a3 ? a1 : a3);
162 }
163 }
164
165 double get_median_pivot(
166 typename current_list_box_t::iterator begin,
167 const typename current_list_box_t::iterator end) {
168 double* ptr = buffer_pivot;
169 while (begin != end) {
170 if ((*begin)->min[dim] > current_box.min[dim]) { *ptr++ = (*begin)->min[dim]; }
171 if ((*begin)->max[dim] < current_box.max[dim]) { *ptr++ = (*begin)->max[dim]; }
172 ++begin;
173 }
174 double* ptr_half = &buffer_pivot[0] + (ptr - buffer_pivot) / 2;
175 std::nth_element(&buffer_pivot[0], ptr_half, ptr);
176 return *ptr_half;
177 }
178
179 static bool aabb_aabb_intersection(const AABBOX& a, const AABBOX& b) {
180 for (int i = 0; i != AABBOX::ndim; ++i) {
181 if (!(a.min[i] < b.max[i] && a.max[i] > b.min[i])) { return false; }
182 }
183 return true;
184 }
185
186 static bool scanning_aabb_aabb_intersection(const AABBOX& a, const AABBOX& b) {
187 for (int i = 1; i != AABBOX::ndim; ++i) {
188 if (!(a.min[i] < b.max[i] && a.max[i] > b.min[i])) { return false; }
189 }
190 return true;
191 }
192
193#define SCANNING 1
194 void rec_compute(
195 typename current_list_box_t::iterator begin,
196 const typename current_list_box_t::iterator end) {
197 /*if (level < 100) {
198 std::cout << "level=" << level << " " << ", size=" << end - begin
199 << ", pair size=" << pair_intersection_list.size() << std::endl;
200 for (int i = 0; i != 3; ++i)
201 std::cout << current_box.min[i] << " " << current_box.max[i] << std::endl;
202 }*/
203 if (end - begin < min_size_to_go_to_scanning) {
204 std::sort(begin, end, ScanningSort());
205 for (; begin != end; ++begin) {
206 for (typename current_list_box_t::const_iterator i = begin + 1; i < end; ++i) {
207#ifdef SCANNING
208 if ((*begin)->max[0] <= (*i)->min[0]) { break; }
209 if (scanning_aabb_aabb_intersection(**begin, **i)) {
210 pair_intersection_list.add(*begin, *i);
211 }
212#else
213 if (aabb_aabb_intersection(**begin, **i)) {
214 pair_intersection_list.add(*begin, *i);
215 }
216#endif
217 }
218 for (typename stack_node_box_t::const_iterator st = stack_node_box.begin();
219 st != stack_node_box.end(); ++st) {
220 for (typename current_list_box_t::const_iterator i = st->first; i != st->second;
221 ++i) {
222 pair_intersection_list.add(*begin, *i);
223 }
224 }
225 }
226 } else {
227 typename current_list_box_t::iterator start_node;
228 if (level == 0) {
229 start_node = end;
230 } else {
231 start_node = std::partition(
232 begin, end, NotInInterval(current_box.min[dim], current_box.max[dim], dim));
233 }
234 if (start_node == begin) {
235 start_node = std::partition(start_node, end, NotInBox(current_box));
236 if (start_node != begin) {
237 if (start_node != end) {
238 stack_node_box.push_back(
239 typename stack_node_box_t::value_type(start_node, end));
240 }
241 dim = ++level % AABBOX::ndim;
242 rec_compute(begin, start_node);
243 dim = --level % AABBOX::ndim;
244 }
245 } else {
246 double pivot;
247 if (start_node - begin < min_size_to_go_to_sort) {
248 pivot = get_median_pivot(begin, start_node);
249 } else {
250 pivot = get_approx_median_pivot(begin, start_node);
251 }
252 // std::cout << "pivot=" << pivot << std::endl;
253
254 start_node = std::partition(start_node, end, NotInBox(current_box));
255 if (start_node != end) {
256 stack_node_box.push_back(
257 typename stack_node_box_t::value_type(start_node, end));
258 }
259
260 typename current_list_box_t::iterator ptr_pivot =
261 std::partition(begin, start_node, Left(pivot, dim));
262
263 if (begin != ptr_pivot) {
264 std::swap(current_box.max[dim], pivot);
265 dim = ++level % AABBOX::ndim;
266 rec_compute(begin, ptr_pivot);
267 dim = --level % AABBOX::ndim;
268 std::swap(current_box.max[dim], pivot);
269 ptr_pivot = std::partition(begin, ptr_pivot, Right(pivot, dim));
270 }
271
272 if (ptr_pivot != start_node) {
273 std::swap(current_box.min[dim], pivot);
274 dim = ++level % AABBOX::ndim;
275 rec_compute(ptr_pivot, start_node);
276 dim = --level % AABBOX::ndim;
277 std::swap(current_box.min[dim], pivot);
278 }
279 }
280 if (start_node != end) {
281 if (start_node != begin) { stack_node_box.pop_back(); }
282 for (; start_node != end; ++start_node) {
283 for (typename current_list_box_t::const_iterator i = start_node + 1; i < end;
284 ++i) {
285 pair_intersection_list.add(*start_node, *i);
286 }
287 for (typename stack_node_box_t::const_iterator st = stack_node_box.begin();
288 st != stack_node_box.end(); ++st) {
289 for (typename current_list_box_t::const_iterator i = st->first;
290 i != st->second; ++i) {
291 pair_intersection_list.add(*start_node, *i);
292 }
293 }
294 }
295 }
296 }
297 }
298
300
301 /*
302 double get_approx_median_pivot(const typename current_list_box_t::iterator begin,
303 const typename current_list_box_t::iterator end, int h = 3) {
304 if (--h == 0) {
305 size_t i;
306 if (get_random(end - begin, i)) {
307 const double r = begin[i]->min[dim];
308 if (r < current_box.min[dim])
309 return current_box.max[dim];
310 else
311 return r;
312 } else {
313 const double r = begin[i]->max[dim];
314 if (r> current_box.max[dim])
315 return current_box.min[dim];
316 else
317 return r;
318 }
319 } else {
320 const double a1 = get_approx_median_pivot(begin, end, h);
321 const double a2 = get_approx_median_pivot(begin, end, h);
322 const double a3 = get_approx_median_pivot(begin, end, h);
323 return a1 < a2 ? (a2 < a3 ? a2 : a3) : (a1 < a3 ? a1 : a3);
324 }
325 }
326
327 double get_pos_pivot(typename current_list_box_t::iterator begin,
328 typename current_list_box_t::iterator& end) {
329
330 }
331
332 double get_median_pivot(typename current_list_box_t::iterator begin,
333 typename current_list_box_t::iterator& end) {
334 double p1, p2, p3;
335 for (;;) {
336 if (end == begin)
337 return 0;
338 if ((*begin)->min[dim] > current_box.min[dim]) {
339 p1 = (*(begin++))->min[dim];
340 break;
341 }
342 if ((*begin)->max[dim] < current_box.max[dim]) {
343 p1 = (*(begin++))->max[dim];
344 break;
345 }
346 std::swap(*begin, *(--end));
347 }
348 for (;;) {
349 if (end == begin)
350 return p1;
351 if (end[-1]->min[dim] > current_box.min[dim]) {
352 p2 = end[-1]->min[dim];
353 break;
354 }
355 if (end[-1]->max[dim] < current_box.max[dim]) {
356 p2 = end[-1]->max[dim];
357 break;
358 }
359 --end;
360 }
361 std::swap(*(begin++), end[-1]);
362 for (;;) {
363 if (end == begin)
364 return p1;
365 typename current_list_box_t::iterator half = begin + (end - begin) / 2;
366 if ((*half)->min[dim] > current_box.min[dim]) {
367 p3 = (*half)->min[dim];
368 break;
369 }
370 if ((*half)->max[dim] < current_box.max[dim]) {
371 p3 = (*(half))->max[dim];
372 break;
373 }
374 std::swap(*half, *(--end));
375 }
376 return p1 < p2 ? (p2 < p3 ? p2 : p3) : (p1 < p3 ? p1 : p3);
377 }
378
379 void rec_compute(const typename current_list_box_t::iterator begin,
380 const typename current_list_box_t::iterator end) {
381 //std::cout << "level=" << level << " " << ", size=" << end - begin
382 //<< ", pair size=" << pair_intersection_list.size() << std::endl;
383 typename current_list_box_t::iterator start_node;
384 if (end - begin < 2)
385 start_node = begin;
386 else {
387 start_node = end;
388 typename current_list_box_t::iterator begin_i = begin;
389 for (;;) {
390 double pivot = get_medina_pivot(begin, start_node);
391 if (begin == start_node) {
392 start_node = std::partition(begin, end, NotInNode(current_box));
393 if (begin != start_node) {
394 dim = ++level % AABBOX::ndim;
395 rec_compute(begin, ptr_pivot);
396 dim = --level % AABBOX::ndim;
397 }
398 goto add_pair;
399 } else {
400 for (;;) {
401 }
402 }
403
404 }
405 }
406
407 if (end - begin < 2)
408 start_node = begin;
409 else {
410 if (end - begin < 3)
411 rand();
412 if (level == 0)
413 start_node = end;
414 else
415 start_node = std::partition(begin, end, NotInNode(current_box));
416 if (start_node != begin) {
417 if (start_node != end)
418 stack_node_box.push_back(typename stack_node_box_t::value_type(start_node, end));
419
420 double pivot;
421 if (start_node - begin < min_size_to_go_to_sort)
422 pivot = get_median_pivot(begin, start_node);
423 else
424 pivot = get_approx_median_pivot(begin, start_node);
425 //std::cout << "pivot=" << pivot << std::endl;
426
427 typename current_list_box_t::iterator ptr_pivot = std::partition(begin, start_node,
428 Left(pivot, dim));
429
430 if (begin != ptr_pivot) {
431 std::swap(current_box.max[dim], pivot);
432 dim = ++level % AABBOX::ndim;
433 rec_compute(begin, ptr_pivot);
434 dim = --level % AABBOX::ndim;
435 std::swap(current_box.max[dim], pivot);
436 ptr_pivot = std::partition(begin, ptr_pivot, Right(pivot, dim));
437 }
438
439 if (ptr_pivot != start_node) {
440 std::swap(current_box.min[dim], pivot);
441 dim = ++level % AABBOX::ndim;
442 rec_compute(ptr_pivot, start_node);
443 dim = --level % AABBOX::ndim;
444 std::swap(current_box.min[dim], pivot);
445 }
446 }
447 }
448
449 add_pair:
450 if (start_node != end) {
451 if (start_node != begin)
452 stack_node_box.pop_back();
453 for (; start_node != end; ++start_node) {
454 for (typename current_list_box_t::const_iterator i = start_node + 1; i
455 < end; ++i)
456 pair_intersection_list.add(*start_node, *i);
457 for (typename stack_node_box_t::const_iterator st = stack_node_box.begin();
458 st != stack_node_box.end(); ++st)
459 for (typename current_list_box_t::const_iterator i = st->first; i != st->second;
460 ++i) pair_intersection_list.add(*start_node, *i);
461 }
462 }
463 }
464 */
465};
466
467template <typename AABBOX, typename PAIR_INTERSECTION_LIST>
468inline void get_aabb_intersection(
469 AABBOX* box_begin, AABBOX* box_end, PAIR_INTERSECTION_LIST& pair_intersection) {
470 AABBIntersection<AABBOX, PAIR_INTERSECTION_LIST>(box_begin, box_end, pair_intersection)
471 .compute();
472}
473
474} // namespace b2000
475
476#endif /* B2AABB_INTERSECTION_H_ */
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32