?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
tweak C1-era CARR
[openc2e.git] / rtree.h
blob6c22589c623e3912a2a4a76b4e955bd255a968da
1 /*
2 * rtree.h
3 * openc2e
4 *
5 * Created by Bryan Donlan on Tue 10 Jan 2006.
6 * Copyright (c) 2006 Bryan Donlan. All rights reserved.
7 *
8 * This library is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2 of the License, or (at your option) any later version.
12 *
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 */
19 #ifndef RTREE_H
20 #define RTREE_H 1
21
22 #include <vector>
23 #include <string>
24 #include <climits>
25 #include <sstream>
26 #include <cassert>
27 #include <iostream>
28 #include <cstdio>
29 #include <algorithm>
30 #include <list>
31 #include <deque>
32 #include <boost/lambda/lambda.hpp>
33 #include "slaballoc.h"
34 #include "openc2e.h"
35
36 // REGION_MAX should be odd, or adjust REGION_MIN manually
37 #define REGION_MAX 5
38 #define REGION_MIN (int)(REGION_MAX / 2)
39
40 typedef unsigned long long u64_t;
41 using namespace boost::lambda;
42
43 struct Region {
44 int xmin, xmax;
45 int ymin, ymax;
46
47 Region() { xmin = xmax = ymin = ymax = 0; }
48
49 template<class Archive>
50 void serialize(Archive & ar, const unsigned int version)
51 {
52 ar & xmin & xmax;
53 ar & ymin & ymax;
54 }
55
56 bool contains(int x, int y) const {
57 return (x >= xmin && x <= xmax &&
58 y >= ymin && y <= ymax
59 );
60 }
61
62 bool overlaps(const Region &r2) const {
63 return (xmin <= r2.xmax && xmax >= r2.xmin)
64 && (ymin <= r2.ymax && ymax >= r2.ymin);
65 }
66
67 Region expand(const Region &r2) {
68 Region r1 = *this;
69 if (r1.xmin > r2.xmin)
70 r1.xmin = r2.xmin;
71 if (r1.xmax < r2.xmax)
72 r1.xmax = r2.xmax;
73 if (r1.ymin > r2.ymin)
74 r1.ymin = r2.ymin;
75 if (r1.ymax < r2.ymax)
76 r1.ymax = r2.ymax;
77 return r1;
78 }
79
80 u64_t area() const {
81 return (u64_t)(xmax - xmin) * (u64_t)(ymax - ymin);
82 }
83
84 std::string to_s() const {
85 std::ostringstream oss;
86 oss << "(" << xmin << " - " << xmax << ", " << ymin << " - " << ymax << ")";
87 return oss.str();
88 }
89
90 bool operator==(const Region &r2) const {
91 return (xmin == r2.xmin)
92 && (xmax == r2.xmax)
93 && (ymin == r2.ymin)
94 && (ymax == r2.ymax)
95 ;
96 }
97
98 bool operator!=(const Region &r2) const {
99 return !(*this == r2);
100 }
101
102 Region(const Region &r2)
103 : xmin(r2.xmin), xmax(r2.xmax),
104 ymin(r2.ymin), ymax(r2.ymax)
105 {
106 assert(xmin <= xmax);
107 assert(ymin <= ymax);
108 }
109
110 Region(int xmin_, int xmax_, int ymin_, int ymax_)
111 : xmin(xmin_), xmax(xmax_), ymin(ymin_), ymax(ymax_)
112 {
113 assert(xmin <= xmax);
114 assert(ymin <= ymax);
115 }
116 };
117
118 template <class T> struct RBranch;
119 template <class T> struct RData;
120 template <class T> class RTree;
121
122 template <class T>
123 struct RNode {
124 Region r;
125 RBranch<T> *parent;
126 unsigned int parent_idx;
127 RTree<T> *tree;
128
129 virtual void scan(const Region &target, std::vector<RData<T> *> &list) = 0;
130 virtual RData<T> *scan(const Region &target) = 0;
131
132 virtual RBranch<T> *chooseleaf(const Region &)
133 { assert(false); }
134
135 RNode(const Region &r_, RTree<T> *t) : r(r_), parent(NULL), parent_idx(UINT_MAX), tree(t)
136 {
137 }
138
139 virtual ~RNode() { }
140
141 virtual std::string dump() = 0;
142
143 virtual size_t size() const = 0; // Leaves
144 virtual size_t inner_size() const = 0; // All nodes
145 virtual void expensive_checks() {}
146 };
147
148 template <class T>
149 struct RBranch : public RNode<T> {
150 unsigned int child_count; // : CLD_BITS;
151 bool is_leaf : 1; // if true, all of children are RLeafs
152 RNode<T> *children[REGION_MAX];
153
154 size_t sz, isz;
155
156 virtual size_t size() const { return sz; }
157
158 virtual size_t inner_size() const { return isz; }
159
160 void expensive_checks() {
161 #ifdef RTREE_DEBUG
162 Region r_ = this->r;
163 for (unsigned int i = 0; i < child_count; i++) {
164 assert(children[i]->parent == this);
165 assert(children[i]->parent_idx == i);
166 children[i]->expensive_checks();
167 }
168 minimize();
169 if (this->r != r_) {
170 std::cerr << "expensive_checks() fail, was = " << r_.to_s() << std::endl;
171 std::cerr << dump();
172 abort();
173 }
174 #endif
175 }
176
177 void minimize() {
178 assert(child_count);
179 this->sz = children[0]->size();
180 this->isz = children[0]->inner_size() + 1;
181 this->r = children[0]->r;
182 for (unsigned int i = 1; i < child_count; i++) {
183 this->r = this->r.expand(children[i]->r);
184 this->sz += children[i]->size();
185 this->isz += children[i]->inner_size();
186 }
187 }
188
189 void recalc_up() {
190 RBranch<T> *p = this;
191 while(p) {
192 p->minimize();
193 p = p->parent;
194 }
195 }
196
197 void add_kid(RNode<T> *n) {
198 assert(child_count < REGION_MAX);
199 children[child_count] = n;
200 n->parent = this;
201 n->parent_idx = child_count;
202 child_count++;
203
204 recalc_up();
205 }
206
207 void drop_kid(unsigned int idx, std::deque<RNode<T> *> &reloc, RBranch<T> *&root) { // not yet used
208 assert(child_count && idx < child_count);
209 child_count--;
210
211 for (unsigned int i = idx; i < child_count; i++) {
212 children[i] = children[i + 1];
213 children[i]->parent_idx = i;
214 }
215
216 if (child_count < REGION_MIN) {
217 reloc.push_back(this);
218 if (this->parent)
219 this->parent->drop_kid(this->parent_idx, reloc, root);
220 else {
221 assert(root == this);
222 root = NULL;
223 }
224 return;
225 }
226
227 recalc_up();
228 }
229
230 void add_multi(RNode<T> **n, int count) {
231 assert(child_count + count <= REGION_MAX);
232 for (int i = 0; i < count; i++) {
233 children[child_count + i] = n[i];
234 children[child_count + i]->parent = this;
235 children[child_count + i]->parent_idx = child_count + i;
236 }
237
238 child_count += count;
239
240 recalc_up();
241 }
242
243 RBranch(const Region &r, RTree<T> *t) : RNode<T>(r, t) {
244 // Convenience for the main RTree initializer
245 child_count = 0;
246 is_leaf = true;
247 this->sz = 0;
248 this->isz = 1;
249 }
250
251 virtual void scan(const Region &target, std::vector<RData<T> *> &list) {
252 if (!this->r.overlaps(target))
253 return;
254 for (unsigned int i = 0; i < child_count; i++)
255 if (children[i]->r.overlaps(target))
256 children[i]->scan(target, list);
257 }
258
259 virtual RData<T> *scan(const Region &target) {
260 if (!this->r.overlaps(target))
261 return NULL;
262 for (unsigned int i = 0; i < child_count; i++) {
263 if (children[i]->r.overlaps(target)) {
264 RData<T> *ret = children[i]->scan(target);
265 if (ret) return ret;
266 }
267 }
268 return NULL;
269 }
270
271 virtual RBranch<T> *chooseleaf(const Region &target) {
272 if (is_leaf)
273 return this;
274 assert(child_count);
275 RBranch<T> *best = dynamic_cast<RBranch<T> *>(children[0]);
276 unsigned long expansion = ULONG_MAX;
277
278 for (unsigned int i = 0; i < child_count; i++) {
279 u64_t orig_area = children[i]->r.area();
280 u64_t new_area = children[i]->r.expand(target).area();
281 u64_t exp = new_area - orig_area;
282 if (exp < expansion) {
283 best = dynamic_cast<RBranch<T> *>(children[i]);
284 assert(best);
285 expansion = exp;
286 }
287 }
288 return best->chooseleaf(target);
289 }
290
291 void insert(RNode<T> *node, RBranch<T> *&root) {
292 if (is_leaf)
293 assert(!dynamic_cast<RBranch<T> *>(node));
294 else
295 assert(dynamic_cast<RBranch<T> *>(node));
296
297 if (child_count == REGION_MAX)
298 split(root);
299 add_kid(node);
300 expensive_checks();
301 }
302
303 void divide(std::list<RNode<T> *> &inlist, std::list<RNode<T> *> &node1, std::list<RNode<T> *> &node2) {
304 // quadratic time split
305 // Described in Guttman, A. 1988. R-trees: a dynamic index structure for spatial searching.
306 // section 3.5.2
307
308 // PickSeeds
309 size_t orig_size = inlist.size();
310
311 assert(!node1.size() && !node2.size());
312
313 typename std::list<RNode<T> *>::iterator best1, best2;
314 u64_t best_d = 0;
315 typename std::list<RNode<T> *>::iterator i1, i2;
316
317 for (i1 = inlist.begin(); i1 != inlist.end(); i1++) {
318 for (i2 = inlist.begin(); i2 != inlist.end(); i2++) {
319 if (i1 == i2)
320 continue;
321 Region j = (*i1)->r.expand((*i2)->r);
322 u64_t d = j.area() - (*i1)->r.area() - (*i2)->r.area();
323 if (d >= best_d) {
324 best_d = d;
325 best1 = i1;
326 best2 = i2;
327 }
328 }
329 }
330
331 Region r1((*best1)->r);
332 Region r2((*best2)->r);
333
334 node1.push_back(*best1);
335 node2.push_back(*best2);
336 inlist.erase(best1);
337 inlist.erase(best2);
338
339
340 while (inlist.size()) {
341 if (inlist.size() + node1.size() == REGION_MIN) {
342 std::copy(inlist.begin(), inlist.end(), std::inserter(node1, node1.end()));
343 inlist.clear();
344 return;
345 }
346 if (inlist.size() + node2.size() == REGION_MIN) {
347 std::copy(inlist.begin(), inlist.end(), std::inserter(node2, node2.end()));
348 inlist.clear();
349 return;
350 }
351 // PickNext
352 typename std::list<RNode<T> *>::iterator best;
353 u64_t best_diff = 0;
354
355 for (typename std::list<RNode<T> *>::iterator i = inlist.begin(); i != inlist.end(); i++) {
356 Region &r = (*i)->r;
357 u64_t d1 = r1.expand(r).area() - r1.area();
358 u64_t d2 = r2.expand(r).area() - r2.area();
359 u64_t diff = (d1 < d2) ? d2 - d1 : d1 - d2;
360 if (diff >= best_diff) {
361 best = i;
362 best_diff = diff;
363 }
364 }
365
366 Region &r = (*best)->r;
367 u64_t d1 = r1.expand(r).area() - r1.area();
368 u64_t d2 = r2.expand(r).area() - r2.area();
369 int choice = 0; // node1, node2
370
371 if (d1 > d2)
372 choice = 2;
373 else if (d2 > d2)
374 choice = 1;
375 else if (r1.area() > r2.area())
376 choice = 2;
377 else if (r2.area() > r1.area())
378 choice = 1;
379 else if (node1.size() > node2.size())
380 choice = 2;
381 else if (node2.size() > node1.size())
382 choice = 2;
383 else choice = 1; // arbitrary tiebreaker
384
385 assert(choice);
386 if (choice == 1) {
387 r1 = r1.expand(r);
388 node1.push_back(*best);
389 } else {
390 r2 = r2.expand(r);
391 node2.push_back(*best);
392 }
393
394 inlist.erase(best);
395 }
396
397 assert(orig_size == node1.size() + node2.size());
398 }
399
400 void split(RBranch<T> *&root) {
401 assert(child_count == REGION_MAX);
402
403 unsigned int split = child_count / 2;
404 RBranch<T> *newbranch = new(this->tree->branches) RBranch<T>(children[split]->r, this->tree);
405
406 newbranch->is_leaf = this->is_leaf;
407
408 std::list<RNode<T> *> l, n1, n2;
409 std::copy(children, children + child_count, std::inserter(l, l.end()));
410 assert(l.size() == child_count);
411 divide(l, n1, n2);
412 assert(l.size() == 0);
413 assert(n1.size() + n2.size() == REGION_MAX);
414
415 std::copy(n1.begin(), n1.end(), children);
416 std::copy(n2.begin(), n2.end(), newbranch->children);
417
418 child_count = n1.size();
419 newbranch->child_count = n2.size();
420
421 for (unsigned int i = 0; i < child_count; i++) {
422 children[i]->parent = this;
423 children[i]->parent_idx = i;
424 }
425
426 for (unsigned int i = 0; i < newbranch->child_count; i++) {
427 newbranch->children[i]->parent = newbranch;
428 newbranch->children[i]->parent_idx = i;
429 }
430
431 minimize();
432 newbranch->minimize();
433
434 if (!this->parent) { // Need to make a new root!
435 root = new(this->tree->branches) RBranch<T>(this->r, this->tree);
436 root->is_leaf = false;
437 root->insert(this, root);
438 }
439
440 this->parent->insert(newbranch, root);
441
442 // std::cerr << "Split! After split, myct " << this->child_count <<
443 // " newbranch " << newbranch->child_count << " root " << root->child_count << std::endl;
444
445 expensive_checks();
446 root->expensive_checks();
447 }
448
449 virtual std::string dump() {
450 assert(child_count);
451 Region test(children[0]->r);
452 std::ostringstream oss;
453 oss << "RBranch " << this->r.to_s() << " {" << std::endl;
454 for (unsigned int i = 0; i < child_count; i++) {
455 oss << "#" << i << " " << children[i]->dump() << std::endl;
456 test = test.expand(children[i]->r);
457 }
458 if (test == this->r)
459 oss << "OK ";
460 else {
461 Region test2(test);
462 oss << "NOK " << test.to_s();
463 test2.expand(this->r);
464 if (test2 != test) oss << " NOTMIN ";
465 }
466 oss << "}" << std::endl;
467 return oss.str();
468 }
469
470 SLAB_CLASS(RBranch<T>)
471 };
472
473 template<class T>
474 struct RData : public RNode<T> {
475 T obj;
476
477 virtual size_t size() const { return 1; }
478
479 virtual size_t inner_size() const { return 1; }
480
481 virtual void scan(const Region &target, std::vector<RData<T> *> &list) {
482 if (!this->r.overlaps(target))
483 return;
484 list.push_back(this);
485 }
486 virtual RData<T> *scan(const Region &target) {
487 if (this->r.overlaps(target))
488 return this;
489 return NULL;
490 }
491
492 RData(const Region &r_, RTree<T> *t, const T &obj_) : RNode<T>(r_, t), obj(obj_) { }
493
494 std::string dump() {
495 std::ostringstream oss;
496 oss << "RData " << this->r.to_s() << " -> " << /* obj.to_s() */ (void *)&obj << std::endl;
497 return oss.str();
498 }
499
500 SLAB_CLASS(RData<T>)
501 };
502
503 template <class T>
504 class RTree {
505 protected:
506 RBranch<T> *root;
507 SlabAllocator branches;
508 DestructingSlab leaves;
509 public:
510
511 void free_node(RNode<T> *n) { delete n; }
512
513 friend class ptr;
514 friend class RBranch<T>;
515
516 class ptr {
517 // Wrapper to enforce encapsulation
518 protected:
519 RData<T> *node;
520 RTree<T> *tree;
521 public:
522 const Region &region() const { assert(node); return node->r; }
523 T &data() { assert(node); return node->obj; }
524 protected:
525 friend std::vector<ptr> RTree::find(const Region &r);
526 ptr(RData<T> *n, RTree<T> *t) : node(n), tree(t) { }
527 ptr(RData<T> &n, RTree<T> *t) : node(&n), tree(t) {}
528 public:
529 ptr(const ptr &p) : node(p.node), tree(p.tree) {}
530 void erase() {
531 std::deque<RNode<T> *> reloc;
532 assert(node && node->parent);
533 node->parent->drop_kid(node->parent_idx, reloc, tree->root);
534 tree->free_node(node);
535
536 while (!reloc.empty()) {
537 RNode<T> *node = reloc.front();
538 reloc.pop_front();
539 RBranch<T> *br = dynamic_cast<RBranch<T> *>(node);
540 if (br) {
541 // This is probably slightly inefficient but should work
542 for (unsigned int i = 0; i < br->child_count; i++) {
543 reloc.push_back(br->children[i]);
544 }
545 tree->free_node(br);
546 } else {
547 if (!tree->root)
548 tree->root = new(this->tree->branches) RBranch<T>(node->r, tree);
549 tree->root->chooseleaf(node->r)->insert(node, tree->root);
550 }
551 }
552 }
553 };
554
555
556 RTree() : root(NULL) {}
557
558 void insert(const Region &r, const T &val) {
559 RData<T> *node = new(leaves) RData<T>(r, this, val);
560 if (!root)
561 root = new(branches) RBranch<T>(r, this);
562 root->chooseleaf(r)->insert(node, root);
563 check();
564 }
565 std::vector<ptr> find(const Region &r) {
566 if (!root) return std::vector<ptr>();
567 std::vector<RData<T> *> list;
568 root->scan(r, list);
569
570 typename std::vector<RData<T> *>::iterator i = list.begin();
571 std::vector<ptr> out;
572 while (i != list.end()) {
573 out.push_back(ptr(*i, this));
574 i++;
575 }
576
577 return out;
578 }
579
580 std::vector<ptr> find(int x, int y) {
581 Region r(x, y, x, y);
582 return find(r);
583 }
584
585 T *find_one(const Region &r) {
586 if (!root) return NULL;
587
588 RData<T> *ret = root->scan(r);
589 if (!ret) return NULL;
590 return &ret->obj;
591 }
592
593 std::string dump() {
594 if (!root)
595 return "NULL";
596 return root->dump();
597 }
598
599
600 template<class Archive>
601 void save(Archive & ar, const unsigned int version) const {
602 std::vector<RData<T> *> l;
603 Region za_warudo(INT_MIN, INT_MIN, INT_MAX, INT_MAX);
604 if (root)
605 root->scan(za_warudo, l);
606 size_t len = l.size();
607 assert(len == size());
608 ar & len;
609 std::for_each(l.start(), l.end(),
610 ar & (*_1)->* &RData<T>::r & (*_1)->* &RData<T>::obj);
611 }
612
613 template<class Archive>
614 void load(Archive & ar, const unsigned int version) {
615 leaves.clear();
616 branches.clear();
617 root = NULL;
618
619 size_t len;
620 ar & len;
621
622 for (size_t i = 0; i < len; i++) {
623 Region r;
624 T obj;
625
626 ar & r & obj;
627 insert(r, obj);
628 }
629 }
630
631 void check() { if (root) root->expensive_checks(); }
632 size_t size() const { return root ? root->size() : 0; }
633 size_t inner_size() const { return root ? root->inner_size() : 0; }
634 };
635
636 #endif
637 // vim: set ts=4 tw=4 noet :
Mirror of openc2e SVN