PeriDyno 1.0.0
Loading...
Searching...
No Matches
SparseOctree.cu
Go to the documentation of this file.
1#pragma once
2#include "SparseOctree.h"
3
4#include <thrust/sort.h>
5
6#include "Object.h"
7
8namespace dyno {
9 DYN_FUNC OcKey CalculateMortonCode(Level l, OcIndex x, OcIndex y, OcIndex z)
10 {
11 OcKey key = 0;
12
13 x |= 1U << l;
14 y |= 1U << l;
15 z |= 1U << l;
16
17 for (int i = 0; i < MAX_LEVEL; i++)
18 {
19 key |= (x & 1U << i) << 2 * i | (y & 1U << i) << (2 * i + 1) | (z & 1U << i) << (2 * i + 2);
20 }
21 return key;
22 }
23
24
25 DYN_FUNC void RecoverFromMortonCode(OcKey key, Level& l, OcIndex& x, OcIndex& y, OcIndex& z)
26 {
27
28 }
29
30 void print(DArray<OctreeNode>& d_arr)
31 {
32 CArray<OctreeNode> h_arr;
33 h_arr.resize(d_arr.size());
34
35 h_arr.assign(d_arr);
36
37 for (uint i = 0; i < h_arr.size(); i++)
38 {
39 Level level;
40 OcIndex tnx, tny, tnz;
41 h_arr[i].getCoord(level, tnx, tny, tnz);
42 //std::cout << "Poster order: " << i << " " << h_arr[i].key() << " " << tnx << " " << tny << " " << tnz << std::endl;
43 printf("Node %d: key: %d - %d: %d %d %d : ", i, h_arr[i].m_key, level, tnx, tny, tnz);
44 printf(" Data size: %d ; start loc: %d; data loc: %d; first child: %d ", h_arr[i].getDataSize(), h_arr[i].getStartIndex(), h_arr[i].m_data_loc, h_arr[i].m_first_child_loc);
45
46 printf("\n child: ");
47 for (int j = 0; j < 8; j++)
48 {
49 printf(" %d: %d ", j, h_arr[i].childs[j]);
50 }
51
52 printf("\n");
53 }
54 h_arr.clear();
55
56 printf("\n\n");
57 }
58
59 DYN_FUNC OctreeNode::OctreeNode()
60 : m_key(0)
61 , m_level(0)
62 {
63 for (int i = 0; i < 8; i++)
64 {
65 childs[i] = EMPTY;
66 }
67 }
68
69 DYN_FUNC OctreeNode::OctreeNode(OcKey key)
70 : m_key(key)
71 , m_level(0)
72 {
73 while (key != 0)
74 {
75 key = key >> 3U;
76 m_level++;
77 }
78
79 m_level = m_level == 0 ? 0 : m_level - 1;
80
81 for (int i = 0; i < 8; i++)
82 {
83 childs[i] = EMPTY;
84 }
85 }
86
87 DYN_FUNC OctreeNode::OctreeNode(Level l, OcIndex x, OcIndex y, OcIndex z)
88 : m_key(0)
89 , m_level(l)
90 {
91 m_key = CalculateMortonCode(l, x, y, z);
92
93 for (int i = 0; i < 8; i++)
94 {
95 childs[i] = EMPTY;
96 }
97 }
98
99
100 DYN_FUNC void OctreeNode::getCoord(Level& l, OcIndex& x, OcIndex& y, OcIndex& z)
101 {
102 l = m_level;
103
104 x = 0;
105 y = 0;
106 z = 0;
107 for (int i = 0; i < MAX_LEVEL; i++)
108 {
109 x |= (m_key & (1U << DIM * i)) >> 2 * i;
110 y |= (m_key & (1U << DIM * i + 1)) >> (2 * i + 1);
111 z |= (m_key & (1U << DIM * i + 2)) >> (2 * i + 2);
112 }
113
114 x &= ((1U << l) - 1);
115 y &= ((1U << l) - 1);
116 z &= ((1U << l) - 1);
117 }
118
119
120 DYN_FUNC bool OctreeNode::operator==(const OctreeNode& mc2) const
121 {
122 return m_key == mc2.m_key;
123 }
124
125 DYN_FUNC bool OctreeNode::operator>=(const OctreeNode& mc2) const
126 {
127 if (mc2.isContainedStrictlyIn(*this))
128 return false;
129
130 auto k1 = m_key;
131 auto k2 = mc2.m_key;
132
133 m_level > mc2.m_level ? (k1 = k1 >> DIM * (m_level - mc2.m_level)) : (k2 = k2 >> DIM * (mc2.m_level - m_level));
134
135 return k1 >= k2;
136 }
137
138 DYN_FUNC bool OctreeNode::operator>(const OctreeNode& mc2) const
139 {
140 if (isContainedStrictlyIn(mc2))
141 {
142 return true;
143 }
144
145 auto k1 = m_key;
146 auto k2 = mc2.m_key;
147
148 m_level > mc2.m_level ? (k1 = k1 >> DIM * (m_level - mc2.m_level)) : (k2 = k2 >> DIM * (mc2.m_level - m_level));
149
150 return k1 > k2;
151 }
152
153 DYN_FUNC bool OctreeNode::operator<=(const OctreeNode& mc2) const
154 {
155 return ~(*this > mc2);
156 }
157
158 DYN_FUNC bool OctreeNode::operator<(const OctreeNode& mc2) const
159 {
160 return ~(*this >= mc2);
161 }
162
163 DYN_FUNC bool OctreeNode::isContainedIn(const OctreeNode& mc2) const
164 {
165 if (m_level < mc2.m_level)
166 {
167 return false;
168 }
169
170 auto k1 = m_key >> 3 * (m_level - mc2.m_level);
171 auto k2 = mc2.key();
172
173 return k1 == k2;
174 }
175
176 DYN_FUNC bool OctreeNode::isContainedStrictlyIn(const OctreeNode& mc2) const
177 {
178 if (m_level <= mc2.m_level)
179 {
180 return false;
181 }
182
183 auto k1 = m_key >> DIM * (m_level - mc2.m_level);
184 auto k2 = mc2.key();
185
186 return k1 == k2;
187 }
188
189
190 DYN_FUNC OctreeNode OctreeNode::leastCommonAncestor(const OctreeNode& mc2) const
191 {
192 OcKey k1 = m_key;
193 OcKey k2 = mc2.m_key;
194
195 m_level > mc2.m_level ? (k1 = k1 >> DIM * (m_level - mc2.m_level)) : (k2 = k2 >> DIM * (mc2.m_level - m_level));
196
197 while (k1 != k2)
198 {
199 k1 = k1 >> 3U;
200 k2 = k2 >> 3U;
201 }
202 return OctreeNode(k1);
203 }
204
205 // DYN_FUNC MortonCode3D& MortonCode3D::operator=(const MortonCode3D & mc2)
206 // {
207 // m_key = mc2.m_key;
208 // return *this;
209 // }
210 //
211 // DYN_FUNC MortonCode3D MortonCode3D::operator=(const MortonCode3D & mc2) const
212 // {
213 // return MortonCode3D(mc2.m_key);
214 // }
215
216 DYN_FUNC bool OctreeNode::operator==(const OcKey k) const
217 {
218 return m_key == k;
219 }
220
221 DYN_FUNC bool OctreeNode::operator!=(const OcKey k) const
222 {
223 return m_key != k;
224 }
225
226
227 template<typename TDataType>
228 SparseOctree<TDataType>::SparseOctree()
229 {
230 }
231
232 template<typename TDataType>
233 SparseOctree<TDataType>::~SparseOctree()
234 {
235 }
236
237 template<typename TDataType>
238 void SparseOctree<TDataType>::release()
239 {
240 m_all_nodes.clear();
241 m_post_ordered_nodes.clear();
242
243 node_buffer.clear();
244 node_count.clear();
245 nonRepeatNodes_cpy.clear();
246 aux_nodes.clear();
247 duplicates_count.clear();
248 }
249
250 template<typename TDataType>
251 void SparseOctree<TDataType>::setSpace(Coord lo, Real h, Real L)
252 {
253 m_lo = lo;
254 m_h = h;
255 m_L = L;
256
257 Real segments = m_L / h;
258
259 m_level_max = ceil(log2(segments));
260 //if (m_level_max < 1) m_level_max = 1;
261 }
262
263 template<typename Real, typename Coord>
264 __global__ void SO_ConstructAABB(
265 DArray<AABB> aabb,
266 DArray<Coord> pos,
267 Real radius)
268 {
269 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
270
271 if (tId >= pos.size()) return;
272
273 Coord p = pos[tId];
274
275 aabb[tId].v0 = p - radius;
276 aabb[tId].v1 = p + radius;
277 }
278
279
280 template<typename TDataType>
281 CPU_FUNC OctreeNode SparseOctree<TDataType>::queryNode(Level l, OcIndex x, OcIndex y, OcIndex z)
282 {
283 OcKey key = CalculateMortonCode(l, x, y, z);
284
285 OcKey mask = 7U << 3 * l;
286
287 CArray<OctreeNode> h_ordered_nodes(m_post_ordered_nodes.size());
288 h_ordered_nodes.assign(m_post_ordered_nodes);
289
290 OctreeNode node = h_ordered_nodes[h_ordered_nodes.size() - 1];
291
292 //root: level 0, i should be indexed from 1 to l
293 for (int i = 1; i <= l; i++)
294 {
295 OcKey mask_i = mask >> DIM * i;
296 int child_index = (key & mask_i) >> (DIM * (l - i));
297
298 if (node.childs[child_index] == EMPTY)
299 {
300 break;
301 }
302 else
303 {
304 node = h_ordered_nodes[node.childs[child_index]];
305 }
306 }
307
308 h_ordered_nodes.clear();
309
310 if (key == node.key())
311 {
312 return node;
313 }
314 return OctreeNode();
315 }
316
317 template<typename TDataType>
318 GPU_FUNC Level SparseOctree<TDataType>::requestLevelNumber(const AABB box)
319 {
320 return SO_ComputeLevel(box, m_h, m_level_max);
321 }
322
323
324 template<typename TDataType>
325 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumber(const AABB box)
326 {
327 Level level;
328
329 int nx_lo, ny_lo, nz_lo;
330 int nx_hi, ny_hi, nz_hi;
331
332 int num = SO_ComputeRange(level, nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, m_lo, m_h, m_level_max);
333
334 int ret_num = 0;
335 //box.intersect(box);
336
337 if (num > 0)
338 {
339 for (int i = nx_lo; i <= nx_hi; i++)
340 for (int j = ny_lo; j <= ny_hi; j++)
341 for (int k = nz_lo; k <= nz_hi; k++)
342 {
343 OcKey mc = CalculateMortonCode(level, i, j, k);
344
345 int d_num = this->requestIntersectionNumber(mc, level);
346
347 ret_num += d_num;
348 }
349 }
350
351 return ret_num;
352 }
353
354 template<typename TDataType>
355 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIds(int* ids, const AABB box)
356 {
357 Level level;
358
359 int nx_lo, ny_lo, nz_lo;
360 int nx_hi, ny_hi, nz_hi;
361
362 int num = SO_ComputeRange(level, nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, m_lo, m_h, m_level_max);
363
364 if (num > 0)
365 {
366 int shift = 0;
367 for (int i = nx_lo; i <= nx_hi; i++)
368 for (int j = ny_lo; j <= ny_hi; j++)
369 for (int k = nz_lo; k <= nz_hi; k++)
370 {
371 OcKey mc = CalculateMortonCode(level, i, j, k);
372
373
374
375 this->reqeustIntersectionIds(ids, shift, mc, level);
376 }
377 }
378 }
379
380
381
382 template<typename TDataType>
383 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumberFromLevel(const AABB box, int level)
384 {
385 int nx_lo, ny_lo, nz_lo;
386 int nx_hi, ny_hi, nz_hi;
387
388 int num = SO_ComputeRangeAtLevel(nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, level, m_lo, m_h, m_level_max);
389
390 int ret_num = 0;
391
392 //printf("num = %d\n",num);
393
394 if (num > 0)
395 {
396 for (int i = nx_lo; i <= nx_hi; i++)
397 for (int j = ny_lo; j <= ny_hi; j++)
398 for (int k = nz_lo; k <= nz_hi; k++)
399 {
400 OcKey mc = CalculateMortonCode(level, i, j, k);
401
402 int d_num = this->requestIntersectionNumber(mc, level);
403
404 ret_num += d_num;
405 }
406 }
407 //if(ret_num > 100)
408 //printf("====\n%d\n%d %d %d\n%d %d %d\n============\n",(int)pow(Real(2), int(level)),nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi);
409 return ret_num;
410 }
411
412 template<typename TDataType>
413 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumberFromLevel(const AABB box, AABB* data, int level)
414 {
415 int nx_lo, ny_lo, nz_lo;
416 int nx_hi, ny_hi, nz_hi;
417
418 int num = SO_ComputeRangeAtLevel(nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, level, m_lo, m_h, m_level_max);
419
420 int ret_num = 0;
421
422 //printf("num = %d\n",num);
423
424 if (num > 0)
425 {
426 for (int i = nx_lo; i <= nx_hi; i++)
427 for (int j = ny_lo; j <= ny_hi; j++)
428 for (int k = nz_lo; k <= nz_hi; k++)
429 {
430 OcKey mc = CalculateMortonCode(level, i, j, k);
431
432 int d_num = this->requestIntersectionNumber(mc, level, box, data);
433
434 ret_num += d_num;
435 }
436 }
437 //if(ret_num > 100)
438 //printf("====\n%d\n%d %d %d\n%d %d %d\n============\n",(int)pow(Real(2), int(level)),nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi);
439 return ret_num;
440 }
441
442 template<typename TDataType>
443 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIdsFromLevel(int* ids, const AABB box, int level)
444 {
445 int nx_lo, ny_lo, nz_lo;
446 int nx_hi, ny_hi, nz_hi;
447
448 int num = SO_ComputeRangeAtLevel(nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, level, m_lo, m_h, m_level_max);
449
450
451 if (num > 0)
452 {
453 int shift = 0;
454 for (int i = nx_lo; i <= nx_hi; i++)
455 for (int j = ny_lo; j <= ny_hi; j++)
456 for (int k = nz_lo; k <= nz_hi; k++)
457 {
458 OcKey mc = CalculateMortonCode(level, i, j, k);
459
460 //printf("level of point: %d %d %d %d\n", level, i, j, k);
461 this->reqeustIntersectionIds(ids, shift, mc, level);
462 }
463 }
464 }
465
466 template<typename TDataType>
467 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIdsFromLevel(int* ids, const AABB box, AABB* data, int level)
468 {
469 int nx_lo, ny_lo, nz_lo;
470 int nx_hi, ny_hi, nz_hi;
471
472 int num = SO_ComputeRangeAtLevel(nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, box, level, m_lo, m_h, m_level_max);
473
474
475 if (num > 0)
476 {
477 int shift = 0;
478 for (int i = nx_lo; i <= nx_hi; i++)
479 for (int j = ny_lo; j <= ny_hi; j++)
480 for (int k = nz_lo; k <= nz_hi; k++)
481 {
482 OcKey mc = CalculateMortonCode(level, i, j, k);
483
484 //printf("level of point: %d %d %d %d\n", level, i, j, k);
485 this->reqeustIntersectionIds(ids, shift, mc, level, box, data);
486 }
487 }
488 }
489
490
491 template<typename TDataType>
492 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumberFromBottom(const AABB box)
493 {
494 return requestIntersectionNumberFromLevel(box, m_level_max);
495 }
496
497 template<typename TDataType>
498 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIdsFromBottom(int* ids, const AABB box)
499 {
500 reqeustIntersectionIdsFromLevel(ids, box, m_level_max);
501 }
502
503 template<typename TDataType>
504 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumberFromBottom(const AABB box, AABB* data)
505 {
506 return requestIntersectionNumberFromLevel(box, data, m_level_max);
507 }
508
509 template<typename TDataType>
510 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIdsFromBottom(int* ids, const AABB box, AABB* data)
511 {
512 reqeustIntersectionIdsFromLevel(ids, box, data, m_level_max);
513 }
514
515 template<typename TDataType>
516 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumber(const OcKey key, const Level l)
517 {
518 int ret_num = 0;
519
520 OcKey mask = 7U << DIM * l;
521
522 OctreeNode node = m_post_ordered_nodes[m_post_ordered_nodes.size() - 1];
523
524 ret_num += node.getDataSize();
525
526
527 int lp = node.level();
528 //root: level 0, i should be indexed from 1 to l
529 for (int i = lp + 1; i <= l;)
530 {
531 OcKey mask_i = mask >> DIM * i;
532 int child_index = (key & mask_i) >> (DIM * (l - i));
533
534 if (node.childs[child_index] == EMPTY)
535 {
536 break;
537 }
538 else
539 {
540 node = m_post_ordered_nodes[node.childs[child_index]];
541
542 int lc = node.level();
543 if (lc - lp > 1)
544 {
545 auto k1 = key >> DIM * (l - lc);
546 auto k2 = node.key();
547 if (!(k1 == k2)) break;
548 }
549
550 ret_num += node.getDataSize();
551 //the octree may be incomplete, skip over missing internal nodes.
552
553 i += (lc - lp);
554 lp = lc;
555 }
556 }
557
558 return ret_num;
559 }
560
561 template<typename TDataType>
562 GPU_FUNC int SparseOctree<TDataType>::requestIntersectionNumber(const OcKey key, const Level l, const AABB box, AABB* data)
563 {
564 int ret_num = 0;
565
566 OcKey mask = 7U << DIM * l;
567
568 OctreeNode node = m_post_ordered_nodes[m_post_ordered_nodes.size() - 1];
569
570 //ret_num += node.getDataSize();
571 int nd_size = node.getDataSize();
572 int t_id = node.getStartIndex();
573 AABB tmp_box;
574 for (int t = 0; t < nd_size; t++)
575 {
576 if (box.intersect(data[m_all_nodes[t_id + t].getDataIndex()], tmp_box))
577 {
578 ret_num++;
579 }
580
581 }
582
583
584 int lp = node.level();
585 //root: level 0, i should be indexed from 1 to l
586 for (int i = lp + 1; i <= l;)
587 {
588 OcKey mask_i = mask >> DIM * i;
589 int child_index = (key & mask_i) >> (DIM * (l - i));
590
591 if (node.childs[child_index] == EMPTY)
592 {
593 break;
594 }
595 else
596 {
597 node = m_post_ordered_nodes[node.childs[child_index]];
598
599 int lc = node.level();
600 if (lc - lp > 1)
601 {
602 auto k1 = key >> DIM * (l - lc);
603 auto k2 = node.key();
604 if (!(k1 == k2)) break;
605 }
606 int t_id = node.getStartIndex();
607 int nd_size = node.getDataSize();
608 AABB tmp_box;
609
610 for (int t = 0; t < nd_size; t++)
611 {
612 if (box.intersect(data[m_all_nodes[t_id + t].getDataIndex()], tmp_box))
613 ret_num++;
614 }
615 //the octree may be incomplete, skip over missing internal nodes.
616
617 i += (lc - lp);
618 lp = lc;
619 }
620 }
621
622 return ret_num;
623 }
624
625
626 template<typename TDataType>
627 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIds(int* ids, int& shift, const OcKey key, const Level l)
628 {
629 OcKey mask = 7U << DIM * l;
630
631 OctreeNode node = m_post_ordered_nodes[m_post_ordered_nodes.size() - 1];
632
633 int nd_size = node.getDataSize();
634
635 int t_id = node.getStartIndex();
636 for (int t = 0; t < nd_size; t++)
637 {
638 ids[shift] = m_all_nodes[t_id + t].getDataIndex();
639 shift++;
640 }
641 //printf("level = %d, shift = %d\n", 0, shift);
642 int lp = node.level();
643 //root: level 0, i should be indexed from 1 to l
644 for (int i = lp + 1; i <= l;)
645 {
646
647 OcKey mask_i = mask >> DIM * i;
648 int child_index = (key & mask_i) >> (DIM * (l - i));
649 //printf("child_index = %d\n", child_index);
650 if (node.childs[child_index] == EMPTY)
651 {
652 break;
653 }
654 else
655 {
656 node = m_post_ordered_nodes[node.childs[child_index]];
657 int nd_size = node.getDataSize();
658
659 int lc = node.level();
660
661 if (lc - lp > 1)
662 {
663 auto k1 = key >> DIM * (l - lc);
664 auto k2 = node.key();
665 if (!(k1 == k2)) break;
666 }
667
668 int t_id = node.getStartIndex();
669 for (int t = 0; t < nd_size; t++)
670 {
671 ids[shift] = m_all_nodes[t_id + t].getDataIndex();
672 shift++;
673 }
674
675 //the octree may be incomplete, skip over missing internal nodes.
676
677 i += (lc - lp);
678 lp = lc;
679
680 Level level;
681 OcIndex x, y, z;
682 node.getCoord(level, x, y, z);
683 //printf("level = %d, shift = %d, level = %d, x = %d, y = %d, z = %d\n", i, shift, level, x, y, z);
684 }
685
686
687 }
688 // printf("shift = %d\n", shift);
689 }
690
691 template<typename TDataType>
692 GPU_FUNC void SparseOctree<TDataType>::reqeustIntersectionIds(int* ids, int& shift, const OcKey key, const Level l, const AABB box, AABB* data)
693 {
694 OcKey mask = 7U << DIM * l;
695
696 OctreeNode node = m_post_ordered_nodes[m_post_ordered_nodes.size() - 1];
697
698 int nd_size = node.getDataSize();
699 AABB tmp_box;
700 int t_id = node.getStartIndex();
701 for (int t = 0; t < nd_size; t++)
702 {
703 if (box.intersect(data[m_all_nodes[t_id + t].getDataIndex()], tmp_box))
704 {
705 ids[shift] = m_all_nodes[t_id + t].getDataIndex();
706 shift++;
707 }
708 //ids[shift] = m_all_nodes[t_id + t].getDataIndex();
709 //shift++;
710 }
711 //printf("level = %d, shift = %d\n", 0, shift);
712 int lp = node.level();
713 //root: level 0, i should be indexed from 1 to l
714 for (int i = lp + 1; i <= l;)
715 {
716
717 OcKey mask_i = mask >> DIM * i;
718 int child_index = (key & mask_i) >> (DIM * (l - i));
719 //printf("child_index = %d\n", child_index);
720 if (node.childs[child_index] == EMPTY)
721 {
722 break;
723 }
724 else
725 {
726 node = m_post_ordered_nodes[node.childs[child_index]];
727 int nd_size = node.getDataSize();
728
729 int lc = node.level();
730 if (lc - lp > 1)
731 {
732 auto k1 = key >> DIM * (l - lc);
733 auto k2 = node.key();
734 if (!(k1 == k2)) break;
735 }
736
737 int t_id = node.getStartIndex();
738 AABB tmp_box;
739 for (int t = 0; t < nd_size; t++)
740 {
741 if (box.intersect(data[m_all_nodes[t_id + t].getDataIndex()], tmp_box))
742 {
743 ids[shift] = m_all_nodes[t_id + t].getDataIndex();
744 shift++;
745 }
746 }
747
748 //the octree may be incomplete, skip over missing internal nodes.
749
750 i += (lc - lp);
751 lp = lc;
752
753 Level level;
754 OcIndex x, y, z;
755 node.getCoord(level, x, y, z);
756 //printf("level = %d, shift = %d, level = %d, x = %d, y = %d, z = %d\n", i, shift, level, x, y, z);
757 }
758
759
760 }
761 // printf("shift = %d\n", shift);
762 }
763 template<typename TDataType>
764 void SparseOctree<TDataType>::construct(
765 const DArray<Coord>& points,
766 Real radius)
767 {
768 DArray<AABB> aabb;
769 aabb.resize(points.size());
770
771
772 cuExecute(points.size(),
773 SO_ConstructAABB,
774 aabb,
775 points,
776 radius);
777
778 this->construct(aabb);
779
780
781 aabb.clear();
782 }
783
784 template<typename Real>
785 DYN_FUNC inline Level SO_ComputeLevel(
786 const AABB box,
787 Real h_min,
788 int level_max)
789 {
790 Real len_max = maximum(box.length(0), maximum(box.length(1), box.length(2)));
791 len_max /= 8.0f;
792 return level_max - clamp(int(ceil(log2(len_max / h_min))), 0, level_max);
793 }
794
795 template<typename Coord>
796 DYN_FUNC int SO_ComputeRange(
797 Level& level,
798 int& nx_lo,
799 int& nx_hi,
800 int& ny_lo,
801 int& ny_hi,
802 int& nz_lo,
803 int& nz_hi,
804 AABB box,
805 Coord origin,
806 Real h_min,
807 int level_max)
808 {
809 Coord lo_rel = box.v0 - origin;
810 Coord hi_rel = box.v1 - origin;
811
812 level = SO_ComputeLevel(box, h_min, level_max);
813
814 int grid_size = (int)pow(Real(2), int(level));
815 Real h = h_min * pow(Real(2), level_max - level);
816
817 // nx_lo = clamp(int(floor(lo_rel[0] / h)), 0, grid_size - 1);
818 // ny_lo = clamp(int(floor(lo_rel[1] / h)), 0, grid_size - 1);
819 // nz_lo = clamp(int(floor(lo_rel[2] / h)), 0, grid_size - 1);
820 //
821 // nx_hi = clamp(int(floor(hi_rel[0] / h)), 0, grid_size - 1);
822 // ny_hi = clamp(int(floor(hi_rel[1] / h)), 0, grid_size - 1);
823 // nz_hi = clamp(int(floor(hi_rel[2] / h)), 0, grid_size - 1);
824
825 nx_lo = (int)floor(lo_rel[0] / h);
826 ny_lo = (int)floor(lo_rel[1] / h);
827 nz_lo = (int)floor(lo_rel[2] / h);
828
829 nx_hi = (int)floor(hi_rel[0] / h);
830 ny_hi = (int)floor(hi_rel[1] / h);
831 nz_hi = (int)floor(hi_rel[2] / h);
832
833
834
835
836 if (nx_hi < 0 || nx_lo >= grid_size || ny_hi < 0 || ny_lo >= grid_size || nz_hi < 0 || nz_lo >= grid_size)
837 {
838 return 0;
839 }
840
841
842 nx_lo = clamp(nx_lo, 0, grid_size - 1);
843 ny_lo = clamp(ny_lo, 0, grid_size - 1);
844 nz_lo = clamp(nz_lo, 0, grid_size - 1);
845
846 nx_hi = clamp(nx_hi, 0, grid_size - 1);
847 ny_hi = clamp(ny_hi, 0, grid_size - 1);
848 nz_hi = clamp(nz_hi, 0, grid_size - 1);
849
850 return (nz_hi - nz_lo + 1) * (ny_hi - ny_lo + 1) * (nx_hi - nx_lo + 1);
851 }
852
853 template<typename Coord>
854 DYN_FUNC int SO_ComputeRangeAtLevel(
855 int& nx_lo,
856 int& nx_hi,
857 int& ny_lo,
858 int& ny_hi,
859 int& nz_lo,
860 int& nz_hi,
861 AABB box,
862 int level,
863 Coord origin,
864 Real h_min,
865 int level_max)
866 {
867 Coord lo_rel = box.v0 - origin;
868 Coord hi_rel = box.v1 - origin;
869
870 int grid_size = (int)pow(Real(2), int(level));
871 Real h = h_min * pow(Real(2), level_max - level);
872
873 nx_lo = (int)floor(lo_rel[0] / h);
874 ny_lo = (int)floor(lo_rel[1] / h);
875 nz_lo = (int)floor(lo_rel[2] / h);
876
877 nx_hi = (int)floor(hi_rel[0] / h);
878 ny_hi = (int)floor(hi_rel[1] / h);
879 nz_hi = (int)floor(hi_rel[2] / h);
880
881 if (nx_hi < 0 || nx_lo >= grid_size || ny_hi < 0 || ny_lo >= grid_size || nz_hi < 0 || nz_lo >= grid_size)
882 {
883 return 0;
884 }
885
886 nx_lo = clamp(nx_lo, 0, grid_size - 1);
887 ny_lo = clamp(ny_lo, 0, grid_size - 1);
888 nz_lo = clamp(nz_lo, 0, grid_size - 1);
889
890 nx_hi = clamp(nx_hi, 0, grid_size - 1);
891 ny_hi = clamp(ny_hi, 0, grid_size - 1);
892 nz_hi = clamp(nz_hi, 0, grid_size - 1);
893
894 return (nz_hi - nz_lo + 1) * (ny_hi - ny_lo + 1) * (nx_hi - nx_lo + 1);
895 }
896
897
898 template<typename Coord>
899 __global__ void SO_InitCounter(
900 DArray<int> counter,
901 DArray<AABB> aabb,
902 Coord origin,
903 Real h_min,
904 int level_max)
905 {
906 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
907
908 if (tId >= aabb.size()) return;
909
910 Level level;
911
912 int nx_lo;
913 int ny_lo;
914 int nz_lo;
915
916 int nx_hi;
917 int ny_hi;
918 int nz_hi;
919
920 counter[tId] = SO_ComputeRange(level, nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, aabb[tId], origin, h_min, level_max);
921 }
922
923 template<typename Coord>
924 __global__ void SO_CreateAllNodes(
925 DArray<OctreeNode> nodes,
926 DArray<AABB> aabb,
927 DArray<int> index,
928 Coord origin,
929 Real h_min,
930 int level_max)
931 {
932 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
933
934 if (tId >= aabb.size()) return;
935
936 Level level;
937
938 int nx_lo;
939 int ny_lo;
940 int nz_lo;
941
942 int nx_hi;
943 int ny_hi;
944 int nz_hi;
945
946 int num = SO_ComputeRange(level, nx_lo, nx_hi, ny_lo, ny_hi, nz_lo, nz_hi, aabb[tId], origin, h_min, level_max);
947
948 if (num > 0)
949 {
950 int acc_num = 0;
951 for (int i = nx_lo; i <= nx_hi; i++)
952 for (int j = ny_lo; j <= ny_hi; j++)
953 for (int k = nz_lo; k <= nz_hi; k++)
954 {
955 auto mc = OctreeNode(level, i, j, k);
956 mc.setDataIndex(tId);
957
958 nodes[index[tId] + acc_num] = mc;
959
960 //printf("%d %d %d %d %d \n", index[tId], morton[index[tId] + acc_num].childs[0], i, j, k);
961
962 acc_num++;
963 }
964 }
965 }
966
967 __global__ void SO_CountNonRepeatedNodes(
968 DArray<int> counter,
969 DArray<OctreeNode> morton)
970 {
971 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
972
973 if (tId >= counter.size()) return;
974
975 if ((tId == 0 || morton[tId].key() != morton[tId - 1].key()))
976 {
977 counter[tId] = 1;
978 }
979 else
980 counter[tId] = 0;
981
982 //printf("%d %d \n", tId, counter[tId]);
983 }
984
985 //the size of counter is 1 more than the size of post_ordered_nodes
986 __global__ void SO_RemoveDuplicativeNodes(
987 DArray<OctreeNode> non_repeated_nodes,
988 DArray<OctreeNode> post_ordered_nodes,
989 DArray<int> counter
990 )
991 {
992 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
993
994 if (tId >= post_ordered_nodes.size()) return;
995
996 if (tId == 0 || post_ordered_nodes[tId].key() != post_ordered_nodes[tId - 1].key())
997 {
998 non_repeated_nodes[counter[tId]] = post_ordered_nodes[tId];
999 non_repeated_nodes[counter[tId]].setStartIndex(tId);
1000 }
1001
1002 //printf("%d %d \n", tId, counter[tId]);
1003 }
1004
1005 __global__ void SO_CalculateDataSize(
1006 DArray<OctreeNode> non_repeated_nodes,
1007 int non_repeated_node_num,
1008 int total_node_num)
1009 {
1010 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1011 if (tId >= non_repeated_node_num) return;
1012
1013 if (tId < non_repeated_node_num - 1)
1014 non_repeated_nodes[tId].setDataSize(non_repeated_nodes[tId + 1].getStartIndex() - non_repeated_nodes[tId].getStartIndex());
1015 else
1016 non_repeated_nodes[tId].setDataSize(total_node_num - non_repeated_nodes[tId].getStartIndex());
1017 }
1018
1019 __global__ void SO_GenerateLCA(
1020 DArray<OctreeNode> nodes,
1021 int shift
1022 )
1023 {
1024 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1025 if (tId >= shift - 1) return;
1026
1027 nodes[tId + shift] = nodes[tId].leastCommonAncestor(nodes[tId + 1]);
1028 }
1029
1030 __global__ void SO_RemoveDuplicativeInternalNodes(
1031 DArray<OctreeNode> nodes,
1032 DArray<int> counter,
1033 DArray<OctreeNode> nodes_cpy)
1034 {
1035 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1036
1037 if (tId >= nodes.size()) return;
1038
1039 if (tId > 0 && nodes[tId].key() == nodes[tId - 1].key())
1040 {
1041 counter[tId] = 0;
1042 nodes[tId] = OctreeNode();
1043 }
1044 else
1045 {
1046 counter[tId] = 1;
1047 nodes[tId] = nodes_cpy[tId];
1048 }
1049 }
1050
1051 __global__ void SO_GenerateAuxNodes(
1052 DArray<OctreeNode> aux_nodes,
1053 DArray<OctreeNode> post_ordered_nodes,
1054 int shift)
1055 {
1056 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1057
1058 if (tId >= post_ordered_nodes.size()) return;
1059
1060 aux_nodes[tId] = post_ordered_nodes[tId];
1061 aux_nodes[tId].m_current_loc = tId;
1062 if (tId < shift - 1)
1063 {
1064 auto cp = post_ordered_nodes[tId].leastCommonAncestor(post_ordered_nodes[tId + 1]);
1065 cp.setFirstChildIndex(tId);
1066 cp.m_bCopy = true;
1067
1068 aux_nodes[tId + shift] = cp;
1069 }
1070 }
1071
1072 __global__ void SO_GeneratePostOrderedNodes(
1073 DArray<OctreeNode> post_ordered_nodes,
1074 DArray<OctreeNode> nodes_buffer,
1075 int num)
1076 {
1077 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1078
1079 if (tId >= num) return;
1080
1081 post_ordered_nodes[tId] = nodes_buffer[tId];
1082 }
1083
1084 __global__ void SO_SetupParentChildRelationship(
1085 DArray<OctreeNode> post_ordered_nodes,
1086 DArray<OctreeNode> aux_nodes)
1087 {
1088 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1089
1090 if (tId >= aux_nodes.size() - 1)
1091 return;
1092
1093 int aux_num = aux_nodes.size();
1094
1095 int t = tId + 1;
1096 if ((aux_nodes[tId].m_bCopy != aux_nodes[t].m_bCopy))
1097 {
1098 int aux_level = aux_nodes[tId].m_level;
1099 int aux_key = aux_nodes[tId].key();
1100 int aux_cur_loc = aux_nodes[tId].m_current_loc;
1101 while (t < aux_num && (aux_nodes[tId].key() == aux_nodes[t].key()))
1102 {
1103 int tLoc = aux_nodes[t].getFirstChildIndex();
1104
1105 OctreeNode child = post_ordered_nodes[tLoc];
1106
1107 auto cKey = child.key();
1108
1109 int cIndex = (cKey >> 3 * (child.m_level - aux_level - 1)) & 7U;
1110
1111 post_ordered_nodes[aux_cur_loc].childs[cIndex] = tLoc;
1112
1113 t++;
1114 }
1115 }
1116 }
1117
1118 template<typename TDataType>
1119 void SparseOctree<TDataType>::construct(const DArray<AABB>& aabb)
1120 {
1121
1122 data_count.resize(aabb.size());
1123
1124
1125 /*************** step 1: identify nodes that containing data ****************/
1126 //step 1.1: calculate the cell number each point covers
1127 cuExecute(data_count.size(),
1128 SO_InitCounter,
1129 data_count,
1130 aabb,
1131 m_lo,
1132 m_h,
1133 m_level_max);
1134
1135 //step 1.2: allocate enough space to store all cells
1136 int total_node_num = thrust::reduce(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), (int)0, thrust::plus<int>());
1137 thrust::exclusive_scan(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), data_count.begin());
1138
1139 m_all_nodes.resize(total_node_num);
1140
1141
1142 //std::cout << MAX_LEVEL << ' ' << m_level_max << std::endl;
1143
1144 //step 1.3: create all octree nodes, record the data location using node.setDataIndex().
1145 cuExecute(aabb.size(),
1146 SO_CreateAllNodes,
1147 m_all_nodes,
1148 aabb,
1149 data_count,
1150 m_lo,
1151 m_h,
1152 m_level_max);
1153
1154 //print(m_all_nodes);
1155 //std::cout << total_node_num << std::endl;
1156 thrust::sort(thrust::device, m_all_nodes.begin(), m_all_nodes.begin() + m_all_nodes.size(), NodeCmp());
1157
1158 //print(m_all_nodes);
1159
1160 construct(m_all_nodes);
1161
1162 //data_count.clear();
1163 //thrust::sort_by_key(thrust::device, m_key.getDataPtr(), m_key.getDataPtr() + m_key.size(), m_morton.getDataPtr());
1164 }
1165
1166 template<typename TDataType>
1167 void SparseOctree<TDataType>::construct(const DArray<OctreeNode>& nodes)
1168 {
1169 /*************** step 2: remove duplicative nodes ****************/
1170 int total_node_num = nodes.size();
1171
1172 duplicates_count.resize(total_node_num);
1173
1174 cuExecute(duplicates_count.size(),
1175 SO_CountNonRepeatedNodes,
1176 duplicates_count,
1177 nodes);
1178
1179 int non_duplicative_num = thrust::reduce(thrust::device, duplicates_count.begin(), duplicates_count.begin() + duplicates_count.size(), (int)0, thrust::plus<int>());
1180 thrust::exclusive_scan(thrust::device, duplicates_count.begin(), duplicates_count.begin() + duplicates_count.size(), duplicates_count.begin());
1181
1182 node_buffer.resize(2 * non_duplicative_num - 1);
1183
1184 //Remove duplicative nodes, record the first location using node.setStartIndex();
1185 cuExecute(nodes.size(),
1186 SO_RemoveDuplicativeNodes,
1187 node_buffer,
1188 nodes,
1189 duplicates_count);
1190
1191 //print(node_buffer);
1192
1193 cuExecute(non_duplicative_num,
1194 SO_CalculateDataSize,
1195 node_buffer,
1196 non_duplicative_num,
1197 total_node_num);
1198
1199 //print(node_buffer);
1200
1201
1202 //thrust::sort(thrust::device, nonRepeatMorton.getDataPtr(), nonRepeatMorton.getDataPtr() + non_repeated_num, NodeCmp());
1203
1204 /*************** step 3: step up parent-child relationship * ****************/
1205
1206 //step 3.1: Generate internal nodes & sort
1207 cuExecute(non_duplicative_num - 1,
1208 SO_GenerateLCA,
1209 node_buffer,
1210 non_duplicative_num);
1211
1212 thrust::sort(thrust::device, node_buffer.begin(), node_buffer.begin() + node_buffer.size(), NodeCmp());
1213
1214 //print(node_buffer);
1215 nonRepeatNodes_cpy.resize(node_buffer.size());
1216 nonRepeatNodes_cpy.assign(node_buffer);
1217
1218 //leaf + LCA
1219 node_count.resize(node_buffer.size());
1220
1221 //step 3.2: remove duplicates
1222 cuExecute(node_buffer.size(),
1223 SO_RemoveDuplicativeInternalNodes,
1224 node_buffer,
1225 node_count,
1226 nonRepeatNodes_cpy);
1227
1228 thrust::sort(thrust::device, node_buffer.begin(), node_buffer.begin() + node_buffer.size(), NodeCmp());
1229 non_duplicative_num = thrust::reduce(thrust::device, node_count.begin(), node_count.begin() + node_count.size(), (int)0, thrust::plus<int>());
1230
1231
1232 m_post_ordered_nodes.resize(non_duplicative_num);
1233
1234 cuExecute(m_post_ordered_nodes.size(),
1235 SO_GeneratePostOrderedNodes,
1236 m_post_ordered_nodes,
1237 node_buffer,
1238 non_duplicative_num);
1239
1240 //print(m_post_ordered_nodes);
1241 aux_nodes.resize(2 * m_post_ordered_nodes.size() - 1);
1242
1243 cuExecute(m_post_ordered_nodes.size(),
1244 SO_GenerateAuxNodes,
1245 aux_nodes,
1246 m_post_ordered_nodes,
1247 m_post_ordered_nodes.size());
1248
1249 //print(morton_b);
1250
1251 thrust::sort(thrust::device, aux_nodes.begin(), aux_nodes.begin() + aux_nodes.size(), NodeCmp());
1252
1253 //print(aux_nodes);
1254
1255 cuExecute(aux_nodes.size(),
1256 SO_SetupParentChildRelationship,
1257 m_post_ordered_nodes,
1258 aux_nodes);
1259
1260 //print(m_post_ordered_nodes);
1261 }
1262
1263 template<typename TDataType>
1264 void SparseOctree<TDataType>::printPostOrderedTree()
1265 {
1266 print(m_post_ordered_nodes);
1267 }
1268
1269
1270 template<typename TDataType>
1271 void SparseOctree<TDataType>::printAllNodes()
1272 {
1273 print(m_all_nodes);
1274 }
1275
1276 DEFINE_CLASS(SparseOctree);
1277}
Generated by doxygen 1.13.2