38 #ifndef COAL_HIERARCHY_TREE_ARRAY_INL_H
39 #define COAL_HIERARCHY_TREE_ARRAY_INL_H
50 namespace implementation_array {
53 template <
typename BV>
55 root_node = NULL_NODE;
58 nodes =
new Node[n_nodes_alloc];
59 for (
size_t i = 0; i < n_nodes_alloc - 1; ++i) nodes[i].next = i + 1;
60 nodes[n_nodes_alloc - 1].next = NULL_NODE;
64 max_lookahead_level = -1;
65 bu_threshold = bu_threshold_;
66 topdown_level = topdown_level_;
70 template <
typename BV>
76 template <
typename BV>
80 init_0(leaves, n_leaves_);
83 init_1(leaves, n_leaves_);
86 init_2(leaves, n_leaves_);
89 init_3(leaves, n_leaves_);
92 init_0(leaves, n_leaves_);
97 template <
typename BV>
101 n_leaves = (size_t)n_leaves_;
102 root_node = NULL_NODE;
103 nodes =
new Node[n_leaves * 2];
104 std::copy(leaves, leaves + n_leaves, nodes);
107 n_nodes_alloc = 2 * n_leaves;
108 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
109 nodes[n_nodes_alloc - 1].next = NULL_NODE;
111 size_t* ids =
new size_t[n_leaves];
112 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
114 root_node = topdown(ids, ids + n_leaves);
118 max_lookahead_level = -1;
122 template <
typename BV>
126 n_leaves = (size_t)n_leaves_;
127 root_node = NULL_NODE;
128 nodes =
new Node[n_leaves * 2];
129 std::copy(leaves, leaves + n_leaves, nodes);
132 n_nodes_alloc = 2 * n_leaves;
133 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
134 nodes[n_nodes_alloc - 1].next = NULL_NODE;
137 if (n_leaves > 0) bound_bv = nodes[0].bv;
138 for (
size_t i = 1; i < n_leaves; ++i) bound_bv += nodes[i].bv;
140 morton_functor<CoalScalar, uint32_t> coder(bound_bv);
141 for (
size_t i = 0; i < n_leaves; ++i)
142 nodes[i].code = coder(nodes[i].bv.center());
144 size_t* ids =
new size_t[n_leaves];
145 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
147 const SortByMorton comp{nodes};
148 std::sort(ids, ids + n_leaves, comp);
149 root_node = mortonRecurse_0(ids, ids + n_leaves, (1 << (coder.bits() - 1)),
156 max_lookahead_level = -1;
160 template <
typename BV>
164 n_leaves = (size_t)n_leaves_;
165 root_node = NULL_NODE;
166 nodes =
new Node[n_leaves * 2];
167 std::copy(leaves, leaves + n_leaves, nodes);
170 n_nodes_alloc = 2 * n_leaves;
171 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
172 nodes[n_nodes_alloc - 1].next = NULL_NODE;
175 if (n_leaves > 0) bound_bv = nodes[0].bv;
176 for (
size_t i = 1; i < n_leaves; ++i) bound_bv += nodes[i].bv;
178 morton_functor<CoalScalar, uint32_t> coder(bound_bv);
179 for (
size_t i = 0; i < n_leaves; ++i)
180 nodes[i].code = coder(nodes[i].bv.center());
182 size_t* ids =
new size_t[n_leaves];
183 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
185 const SortByMorton comp{nodes};
186 std::sort(ids, ids + n_leaves, comp);
187 root_node = mortonRecurse_1(ids, ids + n_leaves, (1 << (coder.bits() - 1)),
194 max_lookahead_level = -1;
198 template <
typename BV>
202 n_leaves = (size_t)n_leaves_;
203 root_node = NULL_NODE;
204 nodes =
new Node[n_leaves * 2];
205 std::copy(leaves, leaves + n_leaves, nodes);
208 n_nodes_alloc = 2 * n_leaves;
209 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
210 nodes[n_nodes_alloc - 1].next = NULL_NODE;
213 if (n_leaves > 0) bound_bv = nodes[0].bv;
214 for (
size_t i = 1; i < n_leaves; ++i) bound_bv += nodes[i].bv;
216 morton_functor<CoalScalar, uint32_t> coder(bound_bv);
217 for (
size_t i = 0; i < n_leaves; ++i)
218 nodes[i].code = coder(nodes[i].bv.center());
220 size_t* ids =
new size_t[n_leaves];
221 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
223 const SortByMorton comp{nodes};
224 std::sort(ids, ids + n_leaves, comp);
225 root_node = mortonRecurse_2(ids, ids + n_leaves);
231 max_lookahead_level = -1;
235 template <
typename BV>
237 size_t node = createNode(NULL_NODE, bv,
data);
238 insertLeaf(root_node, node);
244 template <
typename BV>
252 template <
typename BV>
255 root_node = NULL_NODE;
258 nodes =
new Node[n_nodes_alloc];
259 for (
size_t i = 0; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
260 nodes[n_nodes_alloc - 1].next = NULL_NODE;
264 max_lookahead_level = -1;
268 template <
typename BV>
270 return (n_nodes == 0);
274 template <
typename BV>
276 size_t root = removeLeaf(leaf);
277 if (root != NULL_NODE) {
278 if (lookahead_level > 0) {
280 (i < lookahead_level) && (nodes[root].parent != NULL_NODE); ++i)
281 root = nodes[root].parent;
285 insertLeaf(root, leaf);
289 template <
typename BV>
291 if (nodes[leaf].bv.contain(bv))
return false;
297 template <
typename BV>
304 if (nodes[leaf].bv.contain(bv))
return false;
310 template <
typename BV>
314 if (nodes[leaf].bv.contain(bv))
return false;
320 template <
typename BV>
322 if (root_node == NULL_NODE)
return 0;
324 return getMaxHeight(root_node);
328 template <
typename BV>
330 if (root_node == NULL_NODE)
return 0;
333 getMaxDepth(root_node, 0, max_depth);
338 template <
typename BV>
340 if (root_node != NULL_NODE) {
342 Node* leaves_ = leaves;
343 extractLeaves(root_node, leaves_);
344 root_node = NULL_NODE;
345 std::copy(leaves, leaves + n_leaves, nodes);
348 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
349 nodes[n_nodes_alloc - 1].next = NULL_NODE;
351 size_t* ids =
new size_t[n_leaves];
352 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
354 bottomup(ids, ids + n_leaves);
362 template <
typename BV>
364 if (root_node != NULL_NODE) {
366 Node* leaves_ = leaves;
367 extractLeaves(root_node, leaves_);
368 root_node = NULL_NODE;
369 std::copy(leaves, leaves + n_leaves, nodes);
372 for (
size_t i = n_leaves; i < n_nodes_alloc; ++i) nodes[i].next = i + 1;
373 nodes[n_nodes_alloc - 1].next = NULL_NODE;
375 size_t* ids =
new size_t[n_leaves];
376 for (
size_t i = 0; i < n_leaves; ++i) ids[i] = i;
378 root_node = topdown(ids, ids + n_leaves);
384 template <
typename BV>
386 if (iterations < 0) iterations = (int)n_leaves;
387 if ((root_node != NULL_NODE) && (iterations > 0)) {
388 for (
int i = 0; i < iterations; ++i) {
389 size_t node = root_node;
390 unsigned int bit = 0;
391 while (!nodes[node].isLeaf()) {
392 node = nodes[node].children[(opath >> bit) & 1];
393 bit = (bit + 1) & (
sizeof(
unsigned int) * 8 - 1);
402 template <
typename BV>
404 if (root_node != NULL_NODE) recurseRefit(root_node);
408 template <
typename BV>
410 if (!nodes[root].isLeaf()) {
411 extractLeaves(nodes[root].children[0], leaves);
412 extractLeaves(nodes[root].children[1], leaves);
414 *leaves = nodes[root];
420 template <
typename BV>
426 template <
typename BV>
432 template <
typename BV>
438 template <
typename BV>
440 for (
int i = 0; i < depth; ++i) std::cout <<
" ";
441 Node* n = nodes + root;
442 std::cout <<
" (" << n->
bv.min_[0] <<
", " << n->
bv.min_[1] <<
", "
443 << n->
bv.min_[2] <<
"; " << n->
bv.max_[0] <<
", " << n->
bv.max_[1]
444 <<
", " << n->
bv.max_[2] <<
")" << std::endl;
453 template <
typename BV>
455 if (!nodes[node].isLeaf()) {
456 size_t height1 = getMaxHeight(nodes[node].children[0]);
457 size_t height2 = getMaxHeight(nodes[node].children[1]);
458 return std::max(height1, height2) + 1;
464 template <
typename BV>
466 size_t& max_depth)
const {
467 if (!nodes[node].isLeaf()) {
468 getMaxDepth(nodes[node].children[0], depth + 1, max_depth);
469 getmaxDepth(nodes[node].children[1], depth + 1, max_depth);
471 max_depth = std::max(max_depth, depth);
475 template <
typename BV>
477 size_t* lcur_end = lend;
478 while (lbeg < lcur_end - 1) {
479 size_t *min_it1 =
nullptr, *min_it2 =
nullptr;
480 CoalScalar min_size = (std::numeric_limits<CoalScalar>::max)();
481 for (
size_t* it1 = lbeg; it1 < lcur_end; ++it1) {
482 for (
size_t* it2 = it1 + 1; it2 < lcur_end; ++it2) {
483 CoalScalar cur_size = (nodes[*it1].bv + nodes[*it2].bv).size();
484 if (cur_size < min_size) {
493 createNode(NULL_NODE, nodes[*min_it1].bv, nodes[*min_it2].bv,
nullptr);
494 nodes[p].children[0] = *min_it1;
495 nodes[p].children[1] = *min_it2;
496 nodes[*min_it1].parent = p;
497 nodes[*min_it2].parent = p;
499 size_t tmp = *min_it2;
501 *min_it2 = *lcur_end;
507 template <
typename BV>
509 switch (topdown_level) {
511 return topdown_0(lbeg, lend);
514 return topdown_1(lbeg, lend);
517 return topdown_0(lbeg, lend);
522 template <
typename BV>
524 long num_leaves = lend - lbeg;
525 if (num_leaves > 1) {
526 if (num_leaves > bu_threshold) {
527 BV vol = nodes[*lbeg].bv;
528 for (
size_t* i = lbeg + 1; i < lend; ++i) vol += nodes[*i].bv;
530 size_t best_axis = 0;
531 CoalScalar extent[3] = {vol.width(), vol.height(), vol.depth()};
532 if (extent[1] > extent[0]) best_axis = 1;
533 if (extent[2] > extent[best_axis]) best_axis = 2;
535 nodeBaseLess<BV> comp(nodes, best_axis);
536 size_t* lcenter = lbeg + num_leaves / 2;
537 std::nth_element(lbeg, lcenter, lend, comp);
539 size_t node = createNode(NULL_NODE, vol,
nullptr);
540 nodes[node].children[0] = topdown_0(lbeg, lcenter);
541 nodes[node].children[1] = topdown_0(lcenter, lend);
542 nodes[nodes[node].children[0]].parent = node;
543 nodes[nodes[node].children[1]].parent = node;
546 bottomup(lbeg, lend);
554 template <
typename BV>
556 long num_leaves = lend - lbeg;
557 if (num_leaves > 1) {
558 if (num_leaves > bu_threshold) {
559 Vec3s split_p = nodes[*lbeg].bv.center();
560 BV vol = nodes[*lbeg].bv;
561 for (
size_t* i = lbeg + 1; i < lend; ++i) {
562 split_p += nodes[*i].bv.center();
565 split_p /=
static_cast<CoalScalar>(num_leaves);
567 int bestmidp = (int)num_leaves;
568 int splitcount[3][2] = {{0, 0}, {0, 0}, {0, 0}};
569 for (
size_t* i = lbeg; i < lend; ++i) {
570 Vec3s x = nodes[*i].bv.center() - split_p;
571 for (
int j = 0; j < 3; ++j) ++splitcount[j][x[j] > 0 ? 1 : 0];
574 for (
size_t i = 0; i < 3; ++i) {
575 if ((splitcount[i][0] > 0) && (splitcount[i][1] > 0)) {
576 int midp = std::abs(splitcount[i][0] - splitcount[i][1]);
577 if (midp < bestmidp) {
584 if (best_axis < 0) best_axis = 0;
587 size_t* lcenter = lbeg;
588 for (
size_t* i = lbeg; i < lend; ++i) {
589 if (nodes[*i].bv.center()[best_axis] < split_value) {
597 size_t node = createNode(NULL_NODE, vol,
nullptr);
598 nodes[node].children[0] = topdown_1(lbeg, lcenter);
599 nodes[node].children[1] = topdown_1(lcenter, lend);
600 nodes[nodes[node].children[0]].parent = node;
601 nodes[nodes[node].children[1]].parent = node;
604 bottomup(lbeg, lend);
612 template <
typename BV>
614 const uint32_t& split,
int bits) {
615 long num_leaves = lend - lbeg;
616 if (num_leaves > 1) {
618 const SortByMorton comp{nodes, split};
619 size_t* lcenter = std::lower_bound(lbeg, lend, NULL_NODE, comp);
621 if (lcenter == lbeg) {
622 uint32_t split2 = split | (1 << (bits - 1));
623 return mortonRecurse_0(lbeg, lend, split2, bits - 1);
624 }
else if (lcenter == lend) {
625 uint32_t split1 = (split & (~(1 << bits))) | (1 << (bits - 1));
626 return mortonRecurse_0(lbeg, lend, split1, bits - 1);
628 uint32_t split1 = (split & (~(1 << bits))) | (1 << (bits - 1));
629 uint32_t split2 = split | (1 << (bits - 1));
631 size_t child1 = mortonRecurse_0(lbeg, lcenter, split1, bits - 1);
632 size_t child2 = mortonRecurse_0(lcenter, lend, split2, bits - 1);
633 size_t node = createNode(NULL_NODE,
nullptr);
634 nodes[node].children[0] = child1;
635 nodes[node].children[1] = child2;
636 nodes[child1].parent = node;
637 nodes[child2].parent = node;
641 size_t node = topdown(lbeg, lend);
649 template <
typename BV>
651 const uint32_t& split,
int bits) {
652 long num_leaves = lend - lbeg;
653 if (num_leaves > 1) {
655 const SortByMorton comp{nodes, split};
656 size_t* lcenter = std::lower_bound(lbeg, lend, NULL_NODE, comp);
658 if (lcenter == lbeg) {
659 uint32_t split2 = split | (1 << (bits - 1));
660 return mortonRecurse_1(lbeg, lend, split2, bits - 1);
661 }
else if (lcenter == lend) {
662 uint32_t split1 = (split & (~(1 << bits))) | (1 << (bits - 1));
663 return mortonRecurse_1(lbeg, lend, split1, bits - 1);
665 uint32_t split1 = (split & (~(1 << bits))) | (1 << (bits - 1));
666 uint32_t split2 = split | (1 << (bits - 1));
668 size_t child1 = mortonRecurse_1(lbeg, lcenter, split1, bits - 1);
669 size_t child2 = mortonRecurse_1(lcenter, lend, split2, bits - 1);
670 size_t node = createNode(NULL_NODE,
nullptr);
671 nodes[node].children[0] = child1;
672 nodes[node].children[1] = child2;
673 nodes[child1].parent = node;
674 nodes[child2].parent = node;
678 size_t child1 = mortonRecurse_1(lbeg, lbeg + num_leaves / 2, 0, bits - 1);
679 size_t child2 = mortonRecurse_1(lbeg + num_leaves / 2, lend, 0, bits - 1);
680 size_t node = createNode(NULL_NODE,
nullptr);
681 nodes[node].children[0] = child1;
682 nodes[node].children[1] = child2;
683 nodes[child1].parent = node;
684 nodes[child2].parent = node;
692 template <
typename BV>
694 long num_leaves = lend - lbeg;
695 if (num_leaves > 1) {
696 size_t child1 = mortonRecurse_2(lbeg, lbeg + num_leaves / 2);
697 size_t child2 = mortonRecurse_2(lbeg + num_leaves / 2, lend);
698 size_t node = createNode(NULL_NODE,
nullptr);
699 nodes[node].children[0] = child1;
700 nodes[node].children[1] = child2;
701 nodes[child1].parent = node;
702 nodes[child2].parent = node;
709 template <
typename BV>
711 if (root_node == NULL_NODE) {
713 nodes[leaf].parent = NULL_NODE;
715 if (!nodes[root].isLeaf()) {
717 root = nodes[root].children[
select(leaf, nodes[root].children[0],
718 nodes[root].children[1], nodes)];
719 }
while (!nodes[root].isLeaf());
722 size_t prev = nodes[root].parent;
723 size_t node = createNode(prev, nodes[leaf].bv, nodes[root].bv,
nullptr);
724 if (prev != NULL_NODE) {
725 nodes[prev].children[indexOf(root)] = node;
726 nodes[node].children[0] = root;
727 nodes[root].parent = node;
728 nodes[node].children[1] = leaf;
729 nodes[leaf].parent = node;
731 if (!nodes[prev].bv.contain(nodes[node].bv))
732 nodes[prev].bv = nodes[nodes[prev].children[0]].bv +
733 nodes[nodes[prev].children[1]].bv;
737 }
while (NULL_NODE != (prev = nodes[node].parent));
739 nodes[node].children[0] = root;
740 nodes[root].parent = node;
741 nodes[node].children[1] = leaf;
742 nodes[leaf].parent = node;
749 template <
typename BV>
751 if (leaf == root_node) {
752 root_node = NULL_NODE;
755 size_t parent = nodes[leaf].parent;
756 size_t prev = nodes[parent].parent;
757 size_t sibling = nodes[parent].children[1 - indexOf(leaf)];
759 if (prev != NULL_NODE) {
760 nodes[prev].children[indexOf(parent)] = sibling;
761 nodes[sibling].parent = prev;
763 while (prev != NULL_NODE) {
764 BV new_bv = nodes[nodes[prev].children[0]].bv +
765 nodes[nodes[prev].children[1]].bv;
766 if (!(new_bv == nodes[prev].bv)) {
767 nodes[prev].bv = new_bv;
768 prev = nodes[prev].parent;
773 return (prev != NULL_NODE) ? prev : root_node;
776 nodes[sibling].parent = NULL_NODE;
784 template <
typename BV>
786 size_t root = removeLeaf(leaf);
787 if (root != NULL_NODE) {
788 if (max_lookahead_level >= 0) {
790 (i < max_lookahead_level) && (nodes[root].parent != NULL_NODE); ++i)
791 root = nodes[root].parent;
795 insertLeaf(root, leaf);
800 template <
typename BV>
802 return (nodes[nodes[node].parent].children[1] == node);
806 template <
typename BV>
808 if (freelist == NULL_NODE) {
809 Node* old_nodes = nodes;
811 nodes =
new Node[n_nodes_alloc];
812 std::copy(old_nodes, old_nodes + n_nodes, nodes);
815 for (
size_t i = n_nodes; i < n_nodes_alloc - 1; ++i) nodes[i].next = i + 1;
816 nodes[n_nodes_alloc - 1].next = NULL_NODE;
820 size_t node_id = freelist;
821 freelist = nodes[node_id].next;
822 nodes[node_id].
parent = NULL_NODE;
823 nodes[node_id].children[0] = NULL_NODE;
824 nodes[node_id].children[1] = NULL_NODE;
830 template <
typename BV>
832 const BV& bv2,
void* data) {
833 size_t node = allocateNode();
834 nodes[node].parent = parent;
835 nodes[node].data =
data;
836 nodes[node].bv = bv1 + bv2;
841 template <
typename BV>
843 size_t node = allocateNode();
844 nodes[node].parent = parent;
845 nodes[node].data =
data;
851 template <
typename BV>
853 size_t node = allocateNode();
854 nodes[node].parent = parent;
855 nodes[node].data =
data;
860 template <
typename BV>
862 nodes[node].next = freelist;
868 template <
typename BV>
870 if (!nodes[node].isLeaf()) {
871 recurseRefit(nodes[node].children[0]);
872 recurseRefit(nodes[node].children[1]);
874 nodes[nodes[node].children[0]].bv + nodes[nodes[node].children[1]].bv;
880 template <
typename BV>
882 if ((!nodes[root].isLeaf()) && depth) {
883 fetchLeaves(nodes[root].children[0], leaves, depth - 1);
884 fetchLeaves(nodes[root].children[1], leaves, depth - 1);
887 *leaves = nodes[root];
893 template <
typename BV>
895 : nodes(nodes_),
d(d_) {
900 template <
typename BV>
902 if (nodes[i].bv.center()[(
int)
d] < nodes[j].bv.center()[(
int)
d])
return true;
908 template <
typename S,
typename BV>
910 static bool run(
size_t query,
size_t node1,
size_t node2,
920 static bool run(
const BV& query,
size_t node1,
size_t node2,
932 template <
typename BV>
938 template <
typename BV>
939 size_t select(
const BV& query,
size_t node1,
size_t node2,
945 template <
typename S>
947 static bool run(
size_t query,
size_t node1,
size_t node2,
949 const AABB& bv = nodes[query].
bv;
950 const AABB& bv1 = nodes[node1].
bv;
951 const AABB& bv2 = nodes[node2].
bv;
955 CoalScalar d1 = fabs(v1[0]) + fabs(v1[1]) + fabs(v1[2]);
956 CoalScalar d2 = fabs(v2[0]) + fabs(v2[1]) + fabs(v2[2]);
957 return (
d1 < d2) ? 0 : 1;
960 static bool run(
const AABB& query,
size_t node1,
size_t node2,
962 const AABB& bv = query;
963 const AABB& bv1 = nodes[node1].
bv;
964 const AABB& bv2 = nodes[node2].
bv;
968 CoalScalar d1 = fabs(v1[0]) + fabs(v1[1]) + fabs(v1[2]);
969 CoalScalar d2 = fabs(v2[0]) + fabs(v2[1]) + fabs(v2[2]);
970 return (
d1 < d2) ? 0 : 1;