PeriDyno 1.0.0
Loading...
Searching...
No Matches
VolumeHelper.cu
Go to the documentation of this file.
1#include "VolumeHelper.h"
2
3#include <thrust/sort.h>
4
5namespace dyno
6{
7 template<typename TCoord>
8 struct NodeCmp
9 {
10 DYN_FUNC bool operator()(const VoxelOctreeNode<TCoord>& A, const VoxelOctreeNode<TCoord>& B)
11 {
12 return A > B;
13 }
14 };
15
16 template <typename Real, typename Coord>
17 GPU_FUNC int SO_ComputeGrid(
18 int& nx_lo,
19 int& ny_lo,
20 int& nz_lo,
21 int& nx_hi,
22 int& ny_hi,
23 int& nz_hi,
24 Coord surf_p,
25 Coord surf_q,
26 Coord surf_r,
27 Coord origin_,
28 Real dx_,
29 int nx_,
30 int ny_,
31 int nz_,
32 int extend_band)
33 {
34 Coord fp = (surf_p - origin_) / dx_;
35 Coord fq = (surf_q - origin_) / dx_;
36 Coord fr = (surf_r - origin_) / dx_;
37
38 nx_hi = clamp(int(maximum(fp[0], maximum(fq[0], fr[0]))) + extend_band, 0, nx_ - 1);
39 ny_hi = clamp(int(maximum(fp[1], maximum(fq[1], fr[1]))) + extend_band, 0, ny_ - 1);
40 nz_hi = clamp(int(maximum(fp[2], maximum(fq[2], fr[2]))) + extend_band, 0, nz_ - 1);
41
42 nx_lo = clamp(int(minimum(fp[0], minimum(fq[0], fr[0]))) - extend_band, 0, nx_ - 1);
43 ny_lo = clamp(int(minimum(fp[1], minimum(fq[1], fr[1]))) - extend_band, 0, ny_ - 1);
44 nz_lo = clamp(int(minimum(fp[2], minimum(fq[2], fr[2]))) - extend_band, 0, nz_ - 1);
45
46 if ((nx_hi % 2) != 1) nx_hi++;
47 if ((ny_hi % 2) != 1) ny_hi++;
48 if ((nz_hi % 2) != 1) nz_hi++;
49 if ((nx_lo % 2) != 0) nx_lo--;
50 if ((ny_lo % 2) != 0) ny_lo--;
51 if ((nz_lo % 2) != 0) nz_lo--;
52
53 return (nz_hi - nz_lo + 1) * (ny_hi - ny_lo + 1) * (nx_hi - nx_lo + 1);
54 }
55
56 template <typename Real, typename Coord, typename Triangle>
57 __global__ void SO_SurfaceCount(
58 DArray<int> counter,
59 DArray<Triangle> surf_triangles,
60 DArray<Coord> surf_points,
61 Coord origin_,
62 Real dx_,
63 int nx_,
64 int ny_,
65 int nz_,
66 int extend_band)
67 {
68 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
69
70 if (tId >= counter.size()) return;
71
72 int p = surf_triangles[tId][0];
73 int q = surf_triangles[tId][1];
74 int r = surf_triangles[tId][2];
75
76 int nx_lo;
77 int ny_lo;
78 int nz_lo;
79
80 int nx_hi;
81 int ny_hi;
82 int nz_hi;
83
84 counter[tId] = SO_ComputeGrid(nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi, surf_points[p], surf_points[q], surf_points[r], origin_, dx_, nx_, ny_, nz_, extend_band);
85 }
86
87 template <typename Real, typename Coord, typename Triangle>
88 __global__ void SO_SurfaceInit(
89 DArray<PositionNode> nodes,
90 DArray<int> counter,
91 DArray<Triangle> surf_triangles,
92 DArray<Coord> surf_points,
93 Coord origin_,
94 Real dx_,
95 int nx_,
96 int ny_,
97 int nz_,
98 int extend_band)
99 {
100 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
101
102 if (tId >= counter.size()) return;
103
104 int p = surf_triangles[tId][0];
105 int q = surf_triangles[tId][1];
106 int r = surf_triangles[tId][2];
107
108 int nx_lo;
109 int ny_lo;
110 int nz_lo;
111
112 int nx_hi;
113 int ny_hi;
114 int nz_hi;
115
116 int num = SO_ComputeGrid(nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi, surf_points[p], surf_points[q], surf_points[r], origin_, dx_, nx_, ny_, nz_, extend_band);
117
118 if (num > 0)
119 {
120 int acc_num = 0;
121 for (int k = nz_lo; k <= nz_hi; k++) {
122 for (int j = ny_lo; j <= ny_hi; j++) {
123 for (int i = nx_lo; i <= nx_hi; i++)
124 {
125 OcKey index = CalculateMortonCode(i, j, k);
126 nodes[counter[tId] + acc_num] = PositionNode(tId, index);
127
128 acc_num++;
129 }
130 }
131 }
132 }
133 }
134
135 __global__ void SO_CountNonRepeatedPosition(
136 DArray<int> counter,
137 DArray<PositionNode> nodes)
138 {
139 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
140
141 if (tId >= counter.size()) return;
142
143 if ((tId == 0 || nodes[tId].position_index != nodes[tId - 1].position_index))
144 {
145 counter[tId] = 1;
146 }
147 }
148
149 template <typename Coord, typename Triangle, typename Tri2Edg, typename Edge>
150 GPU_FUNC void SO_ComputeObjectAndNormal(
151 Coord& pobject,
152 Coord& pnormal,
153 DArray<Tri2Edg>& t2e,
154 DArray<Edge>& edge,
155 DArray<Coord>& edgeN,
156 DArray<Coord>& vertexN,
157 DArray<Triangle>& surf_triangles,
158 DArray<Coord>& surf_points,
159 Coord ppos,
160 int surf_id)
161 {
162 int p = surf_triangles[surf_id][0];
163 int q = surf_triangles[surf_id][1];
164 int r = surf_triangles[surf_id][2];
165 Coord p0 = surf_points[p];
166 Coord p1 = surf_points[q];
167 Coord p2 = surf_points[r];
168
169 int eid0 = t2e[surf_id][0];
170 int eid1 = t2e[surf_id][1];
171 int eid2 = t2e[surf_id][2];
172
173 Coord dir = p0 - ppos;
174 Coord e0 = p1 - p0;
175 Coord e1 = p2 - p0;
176 Coord e2 = p2 - p1;
177 Real a = e0.dot(e0);
178 Real b = e0.dot(e1);
179 Real c = e1.dot(e1);
180 Real d = e0.dot(dir);
181 Real e = e1.dot(dir);
182 Real f = dir.dot(dir);
183
184 Real det = a * c - b * b;
185 Real s = b * e - c * d;
186 Real t = b * d - a * e;
187
188 Real maxL = maximum(maximum(e0.norm(), e1.norm()), e2.norm());
189 //handle degenerate triangles
190 if (det < REAL_EPSILON * maxL * maxL)
191 {
192 Real g = e2.normSquared();
193 Real l_max = a;
194
195 Coord op0 = p0;
196 Coord op1 = p1;
197 EKey oe(p, q);
198 if (c > l_max)
199 {
200 op0 = p0;
201 op1 = p2;
202 oe = EKey(p, r);
203
204 l_max = c;
205 }
206 if (g > l_max)
207 {
208 op0 = p1;
209 op1 = p2;
210 oe = EKey(q, r);
211 }
212
213 Coord el = ppos - op0;
214 Coord edir = op1 - op0;
215 if (edir.normSquared() < REAL_EPSILON_SQUARED)
216 {
217 pobject = surf_points[oe[0]];
218 pnormal = vertexN[oe[0]];
219 return;
220 }
221
222 Real et = el.dot(edir) / edir.normSquared();
223
224 if (et <= 0)
225 {
226 pobject = surf_points[oe[0]];
227 pnormal = vertexN[oe[0]];
228 return;
229 }
230 else if (et >= 1)
231 {
232 pobject = surf_points[oe[1]];
233 pnormal = vertexN[oe[1]];
234 return;
235 }
236 else
237 {
238 Coord eq = op0 + et * edir;
239 pobject = eq;
240 if (oe == EKey(edge[eid0][0], edge[eid0][1]))
241 {
242 pnormal = edgeN[eid0];
243 return;
244 }
245 else if (oe == EKey(edge[eid1][0], edge[eid1][1]))
246 {
247 pnormal = edgeN[eid1];
248 return;
249 }
250 else if (oe == EKey(edge[eid2][0], edge[eid2][1]))
251 {
252 pnormal = edgeN[eid2];
253 return;
254 }
255 }
256 }
257 if (s + t <= det)
258 {
259 if (s < 0)
260 {
261 if (t < 0)
262 {
263 //region 4
264 s = 0;
265 t = 0;
266 }
267 else
268 {
269 // region 3
270 s = 0;
271 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e / c));
272 }
273 }
274 else
275 {
276 if (t < 0)
277 {
278 //region 5
279 s = (d >= 0 ? 0 : (-d >= a ? 1 : -d / a));
280 t = 0;
281 }
282 else
283 {
284 //region 0
285 Real invDet = 1 / det;
286 s *= invDet;
287 t *= invDet;
288 }
289 }
290 }
291 else
292 {
293 if (s < 0)
294 {
295 //region 2
296 s = 0;
297 t = 1;
298 }
299 else if (t < 0)
300 {
301 //region 6
302 s = 1;
303 t = 0;
304 }
305 else
306 {
307 //region 1
308 Real numer = c + e - b - d;
309 if (numer <= 0) {
310 s = 0;
311 }
312 else {
313 Real denom = a - 2 * b + c; // positive quantity
314 s = (numer >= denom ? 1 : numer / denom);
315 }
316 t = 1 - s;
317 }
318 }
319 pobject = (p0 + s * e0 + t * e1);
320 if (s == 0 && t == 0)
321 {
322 pnormal = vertexN[p];
323 return;
324 }
325 else if (s == 0 && t == 1)
326 {
327 pnormal = vertexN[r];
328 return;
329 }
330 else if (s == 1 && t == 0)
331 {
332 pnormal = vertexN[q];
333 return;
334 }
335 else if (s == 0 && t < 1)
336 {
337 pnormal = edgeN[eid2];
338 return;
339 }
340 else if (s < 1 && t == 0)
341 {
342 pnormal = edgeN[eid0];
343 return;
344 }
345 else if (s + t == 1)
346 {
347 pnormal = edgeN[eid1];
348 return;
349 }
350 else
351 {
352 pnormal = (p1 - p0).cross(p2 - p0);
353 pnormal.normalize();
354 return;
355 }
356 }
357
358 //注意counter中第i+1个元素存的是0-i元素的和
359 template <typename Real, typename Coord, typename Triangle, typename Tri2Edg, typename Edge>
360 __global__ void SO_FetchNonRepeatedPosition(
361 DArray<VoxelOctreeNode<Coord>> nodes,
362 DArray<Real> nodes_value,
363 DArray<Coord> nodes_object,
364 DArray<Coord> nodes_normal,
365 DArray<IndexNode> x_index,
366 DArray<IndexNode> y_index,
367 DArray<IndexNode> z_index,
368 DArray<PositionNode> all_nodes,
369 DArray<int> counter,
370 DArray<Tri2Edg> t2e,
371 DArray<Edge> edge,
372 DArray<Coord> edgeN,
373 DArray<Coord> vertexN,
374 DArray<Triangle> surf_triangles,
375 DArray<Coord> surf_points,
376 Coord origin_,
377 Real dx_,
378 int nx_,
379 int ny_,
380 int nz_)
381 {
382 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
383
384 if (tId >= counter.size()) return;
385
386 if ((tId == 0 || all_nodes[tId].position_index != all_nodes[tId - 1].position_index))
387 {
388 Coord pnormal(0), pobject(0);
389 int surf_g = all_nodes[tId].surface_index;
390
391 OcIndex gnx, gny, gnz;
392 RecoverFromMortonCode((OcKey)(all_nodes[tId].position_index), gnx, gny, gnz);
393 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
394
395 SO_ComputeObjectAndNormal(pobject, pnormal, t2e, edge, edgeN, vertexN, surf_triangles, surf_points, pos, surf_g);
396
397 Real sign = (pos - pobject).dot(pnormal) < Real(0) ? Real(-1) : Real(1);
398 Real dist = sign * (pos - pobject).norm();
399
400 int acc_num = 1;
401 while (((tId + acc_num) < counter.size()) && (all_nodes[tId + acc_num].position_index == all_nodes[tId].position_index))
402 {
403 int surf_g_i = all_nodes[tId + acc_num].surface_index;
404
405 Coord pnormal_i(0), pobject_i(0);
406 SO_ComputeObjectAndNormal(pobject_i, pnormal_i, t2e, edge, edgeN, vertexN, surf_triangles, surf_points, pos, surf_g_i);
407
408 Real sign_i = (pos - pobject_i).dot(pnormal_i) < Real(0) ? Real(-1) : Real(1);
409 Real dist_i = sign_i * (pos - pobject_i).norm();
410
411 if (std::abs(dist_i) < std::abs(dist))
412 {
413 pnormal = pnormal_i;
414 pobject = pobject_i;
415 dist = dist_i;
416 }
417 acc_num++;
418 }
419
420 VoxelOctreeNode<Coord> mc((Level)0, gnx, gny, gnz, pos);
421 mc.setValueLocation(counter[tId]);
422
423 nodes[counter[tId]] = mc;
424 nodes_value[counter[tId]] = dist;
425 nodes_object[counter[tId]] = pobject;
426 nodes_normal[counter[tId]] = pnormal;
427
428 x_index[counter[tId]] = IndexNode((gnz * nx_ * ny_ + gny * nx_ + gnx), counter[tId]);
429 y_index[counter[tId]] = IndexNode((gnx * ny_ * nz_ + gnz * ny_ + gny), counter[tId]);
430 z_index[counter[tId]] = IndexNode((gny * nz_ * nx_ + gnx * nz_ + gnz), counter[tId]);
431 }
432 }
433
434 template <typename Coord>
435 __global__ void SO_UpdateBottomNeighbors(
436 DArray<VoxelOctreeNode<Coord>> nodes,
437 DArray<IndexNode> x_index,
438 DArray<IndexNode> y_index,
439 DArray<IndexNode> z_index,
440 int nx_,
441 int ny_,
442 int nz_)
443 {
444 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
445 if (tId >= nodes.size()) return;
446
447 int xnx = (x_index[tId].xyz_index) % nx_;
448 if ((xnx != 0 && tId > 0) && (x_index[tId - 1].xyz_index == (x_index[tId].xyz_index - 1)))
449 nodes[x_index[tId].node_index].m_neighbor[0] = x_index[tId - 1].node_index;
450 if ((xnx != (nx_ - 1) && tId < (nodes.size() - 1)) && (x_index[tId + 1].xyz_index == (x_index[tId].xyz_index + 1)))
451 nodes[x_index[tId].node_index].m_neighbor[1] = x_index[tId + 1].node_index;
452
453 int yny = (y_index[tId].xyz_index) % ny_;
454 if ((yny != 0 && tId > 0) && (y_index[tId - 1].xyz_index == (y_index[tId].xyz_index - 1)))
455 nodes[y_index[tId].node_index].m_neighbor[2] = y_index[tId - 1].node_index;
456 if ((yny != (ny_ - 1) && tId < (nodes.size() - 1)) && (y_index[tId + 1].xyz_index == (y_index[tId].xyz_index + 1)))
457 nodes[y_index[tId].node_index].m_neighbor[3] = y_index[tId + 1].node_index;
458
459 int znz = (z_index[tId].xyz_index) % nz_;
460 if ((znz != 0 && tId > 0) && (z_index[tId - 1].xyz_index == (z_index[tId].xyz_index - 1)))
461 nodes[z_index[tId].node_index].m_neighbor[4] = z_index[tId - 1].node_index;
462 if ((znz != (nz_ - 1) && tId < (nodes.size() - 1)) && (z_index[tId + 1].xyz_index == (z_index[tId].xyz_index + 1)))
463 nodes[z_index[tId].node_index].m_neighbor[5] = z_index[tId + 1].node_index;
464 }
465
466 template <typename Real, typename Coord>
467 __global__ void SO_ComputeUpLevelGrids(
468 DArray<VoxelOctreeNode<Coord>> up_level_nodes,
469 DArray<Real> up_level_value,
470 DArray<Coord> up_level_object,
471 DArray<Coord> up_level_normal,
472 DArray<VoxelOctreeNode<Coord>> nodes,
473 DArray<Coord> nodes_object,
474 DArray<Coord> nodes_normal,
475 Coord origin_,
476 Real dx_)
477 {
478 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
479
480 if (tId >= nodes.size()) return;
481
482 Level g_level = nodes[tId].level();
483 g_level++;
484
485 OcIndex fnx, fny, fnz;
486 OcKey cg_key = nodes[tId].key();
487 OcKey fg_key = cg_key >> 6;
488 RecoverFromMortonCode(fg_key, fnx, fny, fnz);
489
490 OcIndex gnx, gny, gnz;
491 int gtId = 7 - (tId % 8);
492 RecoverFromMortonCode((OcKey)gtId, gnx, gny, gnz);
493 Real dx_l = Real(pow(Real(2), int(g_level)) * dx_);
494 Coord pos(origin_[0] + (fnx * 2 + gnx + 0.5) * dx_l, origin_[1] + (fny * 2 + gny + 0.5) * dx_l, origin_[2] + (fnz * 2 + gnz + 0.5) * dx_l);
495
496 int ctId = tId - (tId % 8);
497
498 Real sign = (pos - nodes_object[ctId]).dot(nodes_normal[ctId]) < Real(0) ? Real(-1) : Real(1);
499 Real node_value = sign * (pos - nodes_object[ctId]).norm();
500 int node_index = 0;
501 for (int i = 1; i < 8; i++)
502 {
503 Real sign_i = (pos - nodes_object[ctId + i]).dot(nodes_normal[ctId + i]) < Real(0) ? Real(-1) : Real(1);
504 Real node_value_i = sign_i * (pos - nodes_object[ctId + i]).norm();
505
506 if (abs(node_value_i) < abs(node_value))
507 {
508 node_value = node_value_i;
509 node_index = i;
510 }
511 }
512
513 up_level_nodes[tId] = VoxelOctreeNode<Coord>(g_level, (2 * fnx + gnx), (2 * fny + gny), (2 * fnz + gnz), pos);
514 up_level_nodes[tId].setValueLocation(tId);
515 up_level_value[tId] = node_value;
516 up_level_object[tId] = nodes_object[ctId + node_index];
517 up_level_normal[tId] = nodes_normal[ctId + node_index];
518
519 int gIndex = (cg_key >> 3) & 7U;
520 if (gIndex == gtId)
521 {
522 up_level_nodes[tId].setMidsideNode();
523 up_level_nodes[tId].setChildIndex(ctId);
524 }
525 }
526
527 template <typename Real, typename Coord>
528 __global__ void SO_CountNonRepeatedGrids(
529 DArray<int> counter,
530 DArray<VoxelOctreeNode<Coord>> nodes,
531 DArray<Real> nodes_value)
532 {
533 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
534
535 if (tId >= counter.size()) return;
536
537 if ((tId == 0 || nodes[tId].key() != nodes[tId - 1].key()))
538 {
539 int child_id;
540 bool ismidside = nodes[tId].midside();
541 if (ismidside) child_id = nodes[tId].child();
542
543 int counter_index = tId;
544 Real counter_dist = std::abs(nodes_value[nodes[tId].value()]);
545 int rep_num = 1;
546 while (((tId + rep_num) < counter.size()) && (nodes[tId].key() == nodes[tId + rep_num].key()))
547 {
548
549 ismidside = ismidside || (nodes[tId + rep_num].midside());
550 if (nodes[tId + rep_num].midside())
551 child_id = nodes[tId + rep_num].child();
552
553 if (abs(nodes_value[nodes[tId + rep_num].value()]) < counter_dist)
554 {
555 counter_index = tId + rep_num;
556 counter_dist = abs(nodes_value[nodes[tId + rep_num].value()]);
557 }
558 rep_num++;
559 }
560 if (ismidside)
561 {
562 nodes[counter_index].setMidsideNode();
563 nodes[counter_index].setChildIndex(child_id);
564 }
565 counter[counter_index] = 1;
566 }
567 }
568
569 //注意counter中第i+1个元素存的是0-i元素的和
570 template <typename Real, typename Coord>
571 __global__ void SO_FetchNonRepeatedGrids(
572 DArray<VoxelOctreeNode<Coord>> nodes,
573 DArray<Real> nodes_value,
574 DArray<Coord> nodes_object,
575 DArray<Coord> nodes_normal,
576 DArray<IndexNode> x_index,
577 DArray<IndexNode> y_index,
578 DArray<IndexNode> z_index,
579 DArray<VoxelOctreeNode<Coord>> all_nodes,
580 DArray<Real> all_nodes_value,
581 DArray<Coord> all_nodes_object,
582 DArray<Coord> all_nodes_normal,
583 DArray<int> counter,
584 int nx_,
585 int ny_,
586 int nz_)
587 {
588 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
589
590 if (tId >= counter.size()) return;
591
592 if (tId < (counter.size() - 1) && counter[tId] != counter[tId + 1])
593 {
594 nodes[counter[tId]] = all_nodes[tId];
595 nodes[counter[tId]].setValueLocation(counter[tId]);
596 nodes_value[counter[tId]] = all_nodes_value[all_nodes[tId].value()];
597 nodes_object[counter[tId]] = all_nodes_object[all_nodes[tId].value()];
598 nodes_normal[counter[tId]] = all_nodes_normal[all_nodes[tId].value()];
599
600 OcIndex gnx, gny, gnz;
601 OcKey g_key = all_nodes[tId].key();
602 RecoverFromMortonCode(g_key, gnx, gny, gnz);
603 x_index[counter[tId]] = IndexNode((gnz * nx_ * ny_ + gny * nx_ + gnx), counter[tId]);
604 y_index[counter[tId]] = IndexNode((gnx * ny_ * nz_ + gnz * ny_ + gny), counter[tId]);
605 z_index[counter[tId]] = IndexNode((gny * nz_ * nx_ + gnx * nz_ + gnz), counter[tId]);
606 }
607 else if (tId == (counter.size() - 1) && (counter[tId] < nodes.size()))
608 {
609 nodes[counter[tId]] = all_nodes[tId];
610 nodes[counter[tId]].setValueLocation(counter[tId]);
611 nodes_value[counter[tId]] = all_nodes_value[all_nodes[tId].value()];
612 nodes_object[counter[tId]] = all_nodes_object[all_nodes[tId].value()];
613 nodes_normal[counter[tId]] = all_nodes_normal[all_nodes[tId].value()];
614
615 OcIndex gnx, gny, gnz;
616 OcKey g_key = all_nodes[tId].key();
617 RecoverFromMortonCode(g_key, gnx, gny, gnz);
618 x_index[counter[tId]] = IndexNode((gnz * nx_ * ny_ + gny * nx_ + gnx), counter[tId]);
619 y_index[counter[tId]] = IndexNode((gnx * ny_ * nz_ + gnz * ny_ + gny), counter[tId]);
620 z_index[counter[tId]] = IndexNode((gny * nz_ * nx_ + gnx * nz_ + gnz), counter[tId]);
621 }
622 }
623
624 template <typename Real, typename Coord>
625 __global__ void SO_FIMUpLevelGrids(
626 DArray<VoxelOctreeNode<Coord>> nodes,
627 DArray<Real> nodes_value,
628 DArray<IndexNode> x_index,
629 DArray<IndexNode> y_index,
630 DArray<IndexNode> z_index,
631 DArray<Coord> nodes_object,
632 DArray<Coord> nodes_normal,
633 int nx_,
634 int ny_,
635 int nz_)
636 {
637 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
638 if (tId >= nodes.size()) return;
639
640 int xnx = (x_index[tId].xyz_index) % nx_;
641 if ((xnx != 0 && tId > 0) && (x_index[tId - 1].xyz_index == (x_index[tId].xyz_index - 1)))
642 nodes[x_index[tId].node_index].m_neighbor[0] = x_index[tId - 1].node_index;
643 if ((xnx != (nx_ - 1) && tId < (nodes.size() - 1)) && (x_index[tId + 1].xyz_index == (x_index[tId].xyz_index + 1)))
644 nodes[x_index[tId].node_index].m_neighbor[1] = x_index[tId + 1].node_index;
645
646 int yny = (y_index[tId].xyz_index) % ny_;
647 if ((yny != 0 && tId > 0) && (y_index[tId - 1].xyz_index == (y_index[tId].xyz_index - 1)))
648 nodes[y_index[tId].node_index].m_neighbor[2] = y_index[tId - 1].node_index;
649 if ((yny != (ny_ - 1) && tId < (nodes.size() - 1)) && (y_index[tId + 1].xyz_index == (y_index[tId].xyz_index + 1)))
650 nodes[y_index[tId].node_index].m_neighbor[3] = y_index[tId + 1].node_index;
651
652 int znz = (z_index[tId].xyz_index) % nz_;
653 if ((znz != 0 && tId > 0) && (z_index[tId - 1].xyz_index == (z_index[tId].xyz_index - 1)))
654 nodes[z_index[tId].node_index].m_neighbor[4] = z_index[tId - 1].node_index;
655 if ((znz != (nz_ - 1) && tId < (nodes.size() - 1)) && (z_index[tId + 1].xyz_index == (z_index[tId].xyz_index + 1)))
656 nodes[z_index[tId].node_index].m_neighbor[5] = z_index[tId + 1].node_index;
657
658 __syncthreads();
659
660 for (int i = 0; i < 6; i++)
661 {
662 if (nodes[tId].m_neighbor[i] != EMPTY)
663 {
664 int nb_id = nodes[tId].m_neighbor[i];
665 Real sign = (nodes[tId].position() - nodes_object[nb_id]).dot(nodes_normal[nb_id]) < Real(0) ? Real(-1) : Real(1);
666 Real dist = sign * (nodes[tId].position() - nodes_object[nb_id]).norm();
667
668 if (abs(dist) < abs(nodes_value[tId]))
669 {
670 nodes_value[tId] = dist;
671 nodes_object[tId] = nodes_object[nb_id];
672 nodes_normal[tId] = nodes_normal[nb_id];
673 }
674 }
675 }
676 }
677
678 template <typename Real, typename Coord>
679 __global__ void SO_ComputeTopLevelGrids(
680 DArray<VoxelOctreeNode<Coord>> top_level_nodes,
681 DArray<Coord> top_level_object,
682 DArray<Coord> top_level_normal,
683 DArray<int> top_count,
684 DArray<VoxelOctreeNode<Coord>> nodes,
685 DArray<Coord> nodes_object,
686 DArray<Coord> nodes_normal,
687 Coord origin_,
688 Real dx_,
689 int tnx,
690 int tny,
691 int tnz,
692 Level grid_level)
693 {
694 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
695
696 if (tId >= nodes.size()) return;
697
698 Level cg_level = nodes[tId].level();
699 if (cg_level != (grid_level - 1)) return;
700
701 OcIndex gnx, gny, gnz;
702 OcKey cg_key = nodes[tId].key();
703 OcKey g_key = cg_key >> 3;
704 RecoverFromMortonCode(g_key, gnx, gny, gnz);
705
706 int index = gnx + gny * tnx + gnz * tnx * tny;
707
708 Coord gpos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
709
710 auto mc = VoxelOctreeNode<Coord>(grid_level, g_key);
711
712 if ((tId % 8) == 0)
713 {
714 top_level_nodes[index] = mc;
715 top_level_nodes[index].setMidsideNode();
716 top_level_nodes[index].setChildIndex(tId);
717 top_level_nodes[index].setPosition(gpos);
718 }
719
720 int gIndex = cg_key & 7U;
721 top_level_object[8 * index + gIndex] = nodes_object[tId];
722 top_level_normal[8 * index + gIndex] = nodes_normal[tId];
723 top_count[index] = 1;
724 }
725
726 template <typename Real, typename Coord>
727 __global__ void SO_UpdateTopLevelGrids(
728 DArray<VoxelOctreeNode<Coord>> top_level_nodes,
729 DArray<Real> top_level_val,
730 DArray<Coord> top_level_object,
731 DArray<Coord> top_level_normal,
732 DArray<int> node_ind,
733 DArray<Coord> top_object_buf,
734 DArray<Coord> top_normal_buf,
735 Coord origin_,
736 Real dx_,
737 int tnx,
738 int tny,
739 int tnz,
740 Level grid_level)
741 {
742 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
743
744 if (tId >= top_level_nodes.size()) return;
745
746 int gnz = (int)(tId / (tnx * tny));
747 int gny = (int)((tId % (tnx * tny)) / tnx);
748 int gnx = (tId % (tnx * tny)) % tnx;
749
750 if (node_ind[tId] == 1)
751 {
752 int index = 0;
753 Real sign = (top_level_nodes[tId].position() - top_object_buf[8 * tId]).dot(top_normal_buf[8 * tId]) < Real(0) ? Real(-1) : Real(1);
754 Real min_value = sign * (top_level_nodes[tId].position() - top_object_buf[8 * tId]).norm();
755
756 for (int i = 1; i <= 7; i++)
757 {
758 Real sign_i = (top_level_nodes[tId].position() - top_object_buf[8 * tId + i]).dot(top_normal_buf[8 * tId + i]) < Real(0) ? Real(-1) : Real(1);
759 Real min_value_i = sign_i * (top_level_nodes[tId].position() - top_object_buf[8 * tId + i]).norm();
760
761 if (abs(min_value) > abs(min_value_i))
762 {
763 index = i;
764 min_value = min_value_i;
765 }
766 }
767 top_level_val[tId] = min_value;
768 top_level_object[tId] = top_object_buf[8 * tId + index];
769 top_level_normal[tId] = top_normal_buf[8 * tId + index];
770 top_level_nodes[tId].setValueLocation(tId);
771 }
772 else
773 {
774 Coord gpos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
775
776 auto mc = VoxelOctreeNode<Coord>(grid_level, gnx, gny, gnz);
777 mc.setPosition(gpos);
778 mc.setValueLocation(tId);
779 top_level_nodes[tId] = mc;
780 }
781
782 //update m_neighbor
783 if (gnx > 0)
784 top_level_nodes[tId].m_neighbor[0] = tId - 1;
785 if (gnx < (tnx - 1))
786 top_level_nodes[tId].m_neighbor[1] = tId + 1;
787 if (gny > 0)
788 top_level_nodes[tId].m_neighbor[2] = tId - tnx;
789 if (gny < (tny - 1))
790 top_level_nodes[tId].m_neighbor[3] = tId + tnx;
791 if (gnz > 0)
792 top_level_nodes[tId].m_neighbor[4] = tId - tnx * tny;
793 if (gnz < (tnz - 1))
794 top_level_nodes[tId].m_neighbor[5] = tId + tnx * tny;
795 }
796
797 template <typename Real, typename Coord>
798 DYN_FUNC void SO_ComputeGridWithNeighbor(
799 bool& update,
800 int& update_id,
801 Real& value_id,
802 int index_id,
803 DArray<Coord>& nodes_object,
804 DArray<Coord>& nodes_normal,
805 Coord grid_pos)
806 {
807 Real sign = (grid_pos - nodes_object[index_id]).dot(nodes_normal[index_id]) < Real(0) ? Real(-1) : Real(1);
808 Real dist_value = sign * (grid_pos - nodes_object[index_id]).norm();
809
810 if (abs(dist_value) < abs(value_id))
811 {
812 update_id = index_id;
813 value_id = dist_value;
814 update = true;
815 }
816 }
817
818
819 template <typename Real, typename Coord>
820 __global__ void SO_FIMComputeTopLevelGrids(
821 DArray<VoxelOctreeNode<Coord>> top_level_nodes,
822 DArray<Real> top_level_value,
823 DArray<Coord> top_level_object,
824 DArray<Coord> top_level_normal,
825 DArray<int> node_ind,
826 DArray<Coord> object_temp,
827 DArray<Coord> normal_temp,
828 DArray<int> node_ind_temp,
829 int tnx,
830 int tny,
831 int tnz)
832 {
833 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
834
835 if (tId >= top_level_nodes.size()) return;
836
837 int gnz = (int)(tId / (tnx * tny));
838 int gny = (int)((tId % (tnx * tny)) / tnx);
839 int gnx = (tId % (tnx * tny)) % tnx;
840
841 Coord gpos = top_level_nodes[tId].position();
842
843 bool update = false;
844 int update_id;
845 Real value = std::numeric_limits<Real>::max();
846 if (node_ind_temp[tId] == 1)
847 {
848 value = top_level_value[tId];
849 update_id = tId;
850 }
851 if (gnx > 0)
852 if (node_ind_temp[tId - 1] == 1)
853 {
854 SO_ComputeGridWithNeighbor(update, update_id, value, (tId - 1), object_temp, normal_temp, gpos);
855 }
856 if (gnx < (tnx - 1))
857 if (node_ind_temp[tId + 1] == 1)
858 {
859 SO_ComputeGridWithNeighbor(update, update_id, value, (tId + 1), object_temp, normal_temp, gpos);
860 }
861 if (gny > 0)
862 if (node_ind_temp[tId - tnx] == 1)
863 {
864 SO_ComputeGridWithNeighbor(update, update_id, value, (tId - tnx), object_temp, normal_temp, gpos);
865 }
866 if (gny < (tny - 1))
867 if (node_ind_temp[tId + tnx] == 1)
868 {
869 SO_ComputeGridWithNeighbor(update, update_id, value, (tId + tnx), object_temp, normal_temp, gpos);
870 }
871 if (gnz > 0)
872 if (node_ind_temp[tId - tnx * tny] == 1)
873 {
874 SO_ComputeGridWithNeighbor(update, update_id, value, (tId - tnx * tny), object_temp, normal_temp, gpos);
875 }
876 if (gnz < (tnz - 1))
877 if (node_ind_temp[tId + tnx * tny] == 1)
878 {
879 SO_ComputeGridWithNeighbor(update, update_id, value, (tId + tnx * tny), object_temp, normal_temp, gpos);
880 }
881
882 if (update)
883 {
884 top_level_value[tId] = value;
885 node_ind[tId] = 1;
886 top_level_object[tId] = object_temp[update_id];
887 top_level_normal[tId] = normal_temp[update_id];
888 }
889 }
890
891 template <typename Real, typename Coord>
892 __global__ void SO_CollectionGrids(
893 DArray<VoxelOctreeNode<Coord>> total_nodes,
894 DArray<Real> total_value,
895 DArray<Coord> total_object,
896 DArray<Coord> total_normal,
897 DArray<VoxelOctreeNode<Coord>> level0_nodes,
898 DArray<Real> level0_value,
899 DArray<Coord> level0_object,
900 DArray<Coord> level0_normal,
901 DArray<VoxelOctreeNode<Coord>> level1_nodes,
902 DArray<Real> level1_value,
903 DArray<Coord> level1_object,
904 DArray<Coord> level1_normal,
905 int level0_num,
906 int level0_child_num)
907 {
908 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
909
910 if (tId >= total_nodes.size()) return;
911
912 if (tId >= level0_num)
913 {
914 total_nodes[tId] = level1_nodes[(tId - level0_num)];
915 total_nodes[tId].setValueLocation(tId);
916 if (total_nodes[tId].midside() == true)
917 total_nodes[tId].plusChildIndex(level0_num - level0_child_num);
918 for (int i = 0; i < 6; i++)
919 {
920 if (total_nodes[tId].m_neighbor[i] != EMPTY)
921 total_nodes[tId].m_neighbor[i] += level0_num;
922 }
923
924 total_value[tId] = level1_value[(tId - level0_num)];
925 total_object[tId] = level1_object[(tId - level0_num)];
926 total_normal[tId] = level1_normal[(tId - level0_num)];
927 }
928 else
929 {
930 total_nodes[tId] = level0_nodes[tId];
931 total_value[tId] = level0_value[tId];
932 total_object[tId] = level0_object[tId];
933 total_normal[tId] = level0_normal[tId];
934 }
935 }
936
937 template <typename Real, typename Coord>
938 __global__ void SO_CollectionGridsThree(
939 DArray<VoxelOctreeNode<Coord>> total_nodes,
940 DArray<Real> total_value,
941 DArray<Coord> total_object,
942 DArray<Coord> total_normal,
943 DArray<VoxelOctreeNode<Coord>> level0_nodes,
944 DArray<Real> level0_value,
945 DArray<Coord> level0_object,
946 DArray<Coord> level0_normal,
947 DArray<VoxelOctreeNode<Coord>> level1_nodes,
948 DArray<Real> level1_value,
949 DArray<Coord> level1_object,
950 DArray<Coord> level1_normal,
951 DArray<VoxelOctreeNode<Coord>> levelT_nodes,
952 DArray<Real> levelT_value,
953 DArray<Coord> levelT_object,
954 DArray<Coord> levelT_normal,
955 int level0_num,
956 int level1_num)
957 {
958 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
959
960 if (tId >= total_nodes.size()) return;
961
962 if (tId >= level1_num)
963 {
964 total_nodes[tId] = levelT_nodes[(tId - level1_num)];
965 total_nodes[tId].setValueLocation(tId);
966 if (total_nodes[tId].midside() == true)
967 total_nodes[tId].plusChildIndex(level0_num);
968 for (int i = 0; i < 6; i++)
969 {
970 if (total_nodes[tId].m_neighbor[i] != EMPTY)
971 total_nodes[tId].m_neighbor[i] += level1_num;
972 }
973
974 total_value[tId] = levelT_value[(tId - level1_num)];
975 total_object[tId] = levelT_object[(tId - level1_num)];
976 total_normal[tId] = levelT_normal[(tId - level1_num)];
977 }
978 else if (tId >= level0_num)
979 {
980 total_nodes[tId] = level1_nodes[(tId - level0_num)];
981 total_nodes[tId].setValueLocation(tId);
982 for (int i = 0; i < 6; i++)
983 {
984 if (total_nodes[tId].m_neighbor[i] != EMPTY)
985 total_nodes[tId].m_neighbor[i] += level0_num;
986 }
987
988 total_value[tId] = level1_value[(tId - level0_num)];
989 total_object[tId] = level1_object[(tId - level0_num)];
990 total_normal[tId] = level1_normal[(tId - level0_num)];
991 }
992 else
993 {
994 total_nodes[tId] = level0_nodes[tId];
995 total_nodes[tId].setValueLocation(tId);
996
997 total_value[tId] = level0_value[tId];
998 total_object[tId] = level0_object[tId];
999 total_normal[tId] = level0_normal[tId];
1000 }
1001 }
1002
1003 template <typename Real, typename Coord>
1004 __global__ void SO_CollectionGridsFour(
1005 DArray<VoxelOctreeNode<Coord>> total_nodes,
1006 DArray<Real> total_value,
1007 DArray<Coord> total_object,
1008 DArray<Coord> total_normal,
1009 DArray<VoxelOctreeNode<Coord>> level0_nodes,
1010 DArray<Real> level0_value,
1011 DArray<Coord> level0_object,
1012 DArray<Coord> level0_normal,
1013 DArray<VoxelOctreeNode<Coord>> level1_nodes,
1014 DArray<Real> level1_value,
1015 DArray<Coord> level1_object,
1016 DArray<Coord> level1_normal,
1017 DArray<VoxelOctreeNode<Coord>> level2_nodes,
1018 DArray<Real> level2_value,
1019 DArray<Coord> level2_object,
1020 DArray<Coord> level2_normal,
1021 DArray<VoxelOctreeNode<Coord>> levelT_nodes,
1022 DArray<Real> levelT_value,
1023 DArray<Coord> levelT_object,
1024 DArray<Coord> levelT_normal,
1025 int level0_num,
1026 int level1_num,
1027 int level2_num)
1028 {
1029 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1030
1031 if (tId >= total_nodes.size()) return;
1032
1033 if (tId >= level2_num)
1034 {
1035 total_nodes[tId] = levelT_nodes[(tId - level2_num)];
1036 total_nodes[tId].setValueLocation(tId);
1037 if (total_nodes[tId].midside() == true)
1038 total_nodes[tId].plusChildIndex(level1_num);
1039 for (int i = 0; i < 6; i++)
1040 {
1041 if (total_nodes[tId].m_neighbor[i] != EMPTY)
1042 total_nodes[tId].m_neighbor[i] += level2_num;
1043 }
1044
1045 total_value[tId] = levelT_value[(tId - level2_num)];
1046 total_object[tId] = levelT_object[(tId - level2_num)];
1047 total_normal[tId] = levelT_normal[(tId - level2_num)];
1048 }
1049 else if (tId >= level1_num)
1050 {
1051 total_nodes[tId] = level2_nodes[(tId - level1_num)];
1052 total_nodes[tId].setValueLocation(tId);
1053 if (total_nodes[tId].midside() == true)
1054 total_nodes[tId].plusChildIndex(level0_num);
1055 for (int i = 0; i < 6; i++)
1056 {
1057 if (total_nodes[tId].m_neighbor[i] != EMPTY)
1058 total_nodes[tId].m_neighbor[i] += level1_num;
1059 }
1060
1061 total_value[tId] = level2_value[(tId - level1_num)];
1062 total_object[tId] = level2_object[(tId - level1_num)];
1063 total_normal[tId] = level2_normal[(tId - level1_num)];
1064 }
1065 else if (tId >= level0_num)
1066 {
1067 total_nodes[tId] = level1_nodes[(tId - level0_num)];
1068 total_nodes[tId].setValueLocation(tId);
1069 for (int i = 0; i < 6; i++)
1070 {
1071 if (total_nodes[tId].m_neighbor[i] != EMPTY)
1072 total_nodes[tId].m_neighbor[i] += level0_num;
1073 }
1074
1075 total_value[tId] = level1_value[(tId - level0_num)];
1076 total_object[tId] = level1_object[(tId - level0_num)];
1077 total_normal[tId] = level1_normal[(tId - level0_num)];
1078 }
1079 else
1080 {
1081 total_nodes[tId] = level0_nodes[tId];
1082 total_value[tId] = level0_value[tId];
1083 total_object[tId] = level0_object[tId];
1084 total_normal[tId] = level0_normal[tId];
1085 }
1086 }
1087
1088 template <typename Real, typename Coord>
1089 __global__ void SO_CollectionGridsFive(
1090 DArray<VoxelOctreeNode<Coord>> total_nodes,
1091 DArray<Real> total_value,
1092 DArray<Coord> total_object,
1093 DArray<Coord> total_normal,
1094 DArray<VoxelOctreeNode<Coord>> level0_nodes,
1095 DArray<Real> level0_value,
1096 DArray<Coord> level0_object,
1097 DArray<Coord> level0_normal,
1098 DArray<VoxelOctreeNode<Coord>> level1_nodes,
1099 DArray<Real> level1_value,
1100 DArray<Coord> level1_object,
1101 DArray<Coord> level1_normal,
1102 DArray<VoxelOctreeNode<Coord>> level2_nodes,
1103 DArray<Real> level2_value,
1104 DArray<Coord> level2_object,
1105 DArray<Coord> level2_normal,
1106 DArray<VoxelOctreeNode<Coord>> level3_nodes,
1107 DArray<Real> level3_value,
1108 DArray<Coord> level3_object,
1109 DArray<Coord> level3_normal,
1110 DArray<VoxelOctreeNode<Coord>> levelT_nodes,
1111 DArray<Real> levelT_value,
1112 DArray<Coord> levelT_object,
1113 DArray<Coord> levelT_normal,
1114 int level0_num,
1115 int level1_num,
1116 int level2_num,
1117 int level3_num)
1118 {
1119 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1120
1121 if (tId >= total_nodes.size()) return;
1122
1123 if (tId >= level3_num)
1124 {
1125 total_nodes[tId] = levelT_nodes[(tId - level3_num)];
1126 total_value[tId] = levelT_value[(tId - level3_num)];
1127 total_object[tId] = levelT_object[(tId - level3_num)];
1128 total_normal[tId] = levelT_normal[(tId - level3_num)];
1129 }
1130 else if (tId >= level2_num)
1131 {
1132 total_nodes[tId] = level3_nodes[(tId - level2_num)];
1133 total_value[tId] = level3_value[(tId - level2_num)];
1134 total_object[tId] = level3_object[(tId - level2_num)];
1135 total_normal[tId] = level3_normal[(tId - level2_num)];
1136 }
1137 else if (tId >= level1_num)
1138 {
1139 total_nodes[tId] = level2_nodes[(tId - level1_num)];
1140 total_value[tId] = level2_value[(tId - level1_num)];
1141 total_object[tId] = level2_object[(tId - level1_num)];
1142 total_normal[tId] = level2_normal[(tId - level1_num)];
1143 }
1144 else if (tId >= level0_num)
1145 {
1146 total_nodes[tId] = level1_nodes[(tId - level0_num)];
1147 total_value[tId] = level1_value[(tId - level0_num)];
1148 total_object[tId] = level1_object[(tId - level0_num)];
1149 total_normal[tId] = level1_normal[(tId - level0_num)];
1150 }
1151 else
1152 {
1153 total_nodes[tId] = level0_nodes[tId];
1154 total_value[tId] = level0_value[tId];
1155 total_object[tId] = level0_object[tId];
1156 total_normal[tId] = level0_normal[tId];
1157 }
1158 total_nodes[tId].setValueLocation(tId);
1159 }
1160
1161 template<typename TDataType>
1162 void VolumeHelper<TDataType>::levelBottom(
1163 DArray<VoxelOctreeNode<Coord>>& grid0,
1164 DArray<Real>& grid0_value,
1165 DArray<Coord>& grid0_object,
1166 DArray<Coord>& grid0_normal,
1167 std::shared_ptr<TriangleSet<TDataType>> triSet,
1168 Coord m_origin,
1169 int m_nx,
1170 int m_ny,
1171 int m_nz,
1172 int padding,
1173 int& m_level0,
1174 Real m_dx)
1175 {
1176 // initialize data
1177 auto& triangles = triSet->getTriangles();
1178 auto& points = triSet->getPoints();
1179 auto& edges = triSet->getEdges();
1180 std::printf("the num of points edges surfaces is: %d %d %d \n", points.size(), edges.size(), triangles.size());
1181
1182 triSet->updateTriangle2Edge();
1183 auto& tri2edg = triSet->getTriangle2Edge();
1184
1185 DArray<Coord> edgeNormal, vertexNormal;
1186 triSet->updateEdgeNormal(edgeNormal);
1187 triSet->updateAngleWeightedVertexNormal(vertexNormal);
1188
1189 DArray<int> data_count;
1190 data_count.resize(triangles.size());
1191
1192 //数一下level_0中active grid的数目
1193 cuExecute(data_count.size(),
1194 SO_SurfaceCount,
1195 data_count,
1196 triangles,
1197 points,
1198 m_origin,
1199 m_dx,
1200 m_nx,
1201 m_ny,
1202 m_nz,
1203 padding);
1204
1205 int grid_num = thrust::reduce(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), (int)0, thrust::plus<int>());
1206 thrust::exclusive_scan(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), data_count.begin());
1207
1208 DArray<PositionNode> grid_buf;
1209 grid_buf.resize(grid_num);
1210
1211 //将level_0中的active grid取出
1212 cuExecute(data_count.size(),
1213 SO_SurfaceInit,
1214 grid_buf,
1215 data_count,
1216 triangles,
1217 points,
1218 m_origin,
1219 m_dx,
1220 m_nx,
1221 m_ny,
1222 m_nz,
1223 padding);
1224
1225 thrust::sort(thrust::device, grid_buf.begin(), grid_buf.begin() + grid_buf.size(), PositionCmp());
1226
1227 data_count.resize(grid_num);
1228 data_count.reset();
1229 //数不重复的grid
1230 cuExecute(data_count.size(),
1231 SO_CountNonRepeatedPosition,
1232 data_count,
1233 grid_buf);
1234
1235 int grid0_num = thrust::reduce(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), (int)0, thrust::plus<int>());
1236 thrust::exclusive_scan(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), data_count.begin());
1237
1238 DArray<IndexNode> xIndex, yIndex, zIndex;
1239 xIndex.resize(grid0_num);
1240 yIndex.resize(grid0_num);
1241 zIndex.resize(grid0_num);
1242 grid0.resize(grid0_num);
1243 grid0_value.resize(grid0_num);
1244 grid0_object.resize(grid0_num);
1245 grid0_normal.resize(grid0_num);
1246 //去重
1247 cuExecute(data_count.size(),
1248 SO_FetchNonRepeatedPosition,
1249 grid0,
1250 grid0_value,
1251 grid0_object,
1252 grid0_normal,
1253 xIndex,
1254 yIndex,
1255 zIndex,
1256 grid_buf,
1257 data_count,
1258 tri2edg,
1259 edges,
1260 edgeNormal,
1261 vertexNormal,
1262 triangles,
1263 points,
1264 m_origin,
1265 m_dx,
1266 m_nx,
1267 m_ny,
1268 m_nz);
1269
1270 m_level0 = grid0_num;
1271 //std::printf("the grid num of level 0 is: %d \n", grid0_num);
1272
1273 thrust::sort(thrust::device, xIndex.begin(), xIndex.begin() + xIndex.size(), IndexCmp());
1274 thrust::sort(thrust::device, yIndex.begin(), yIndex.begin() + yIndex.size(), IndexCmp());
1275 thrust::sort(thrust::device, zIndex.begin(), zIndex.begin() + zIndex.size(), IndexCmp());
1276
1277 cuExecute(grid0_num,
1278 SO_UpdateBottomNeighbors,
1279 grid0,
1280 xIndex,
1281 yIndex,
1282 zIndex,
1283 m_nx,
1284 m_ny,
1285 m_nz);
1286
1287 xIndex.clear();
1288 yIndex.clear();
1289 zIndex.clear();
1290 data_count.clear();
1291 grid_buf.clear();
1292 edgeNormal.clear();
1293 vertexNormal.clear();
1294 }
1295
1296 template<typename TDataType>
1297 void VolumeHelper<TDataType>::levelMiddle(
1298 DArray<VoxelOctreeNode<Coord>>& grid1,
1299 DArray<Real>& grid1_value,
1300 DArray<Coord>& grid1_object,
1301 DArray<Coord>& grid1_normal,
1302 DArray<VoxelOctreeNode<Coord>>& grid0,
1303 DArray<Coord>& grid0_object,
1304 DArray<Coord>& grid0_normal,
1305 Coord m_origin,
1306 Level multi_level,
1307 int m_nx,
1308 int m_ny,
1309 int m_nz,
1310 Real m_dx)
1311 {
1312 Real coef = Real(pow(Real(2), int(multi_level)));
1313 int up_nx = m_nx / coef;
1314 int up_ny = m_ny / coef;
1315 int up_nz = m_nz / coef;
1316
1317 int grid0_num = grid0.size();
1318 DArray<VoxelOctreeNode<Coord>> grid_buf1;
1319 DArray<Real> grid_value_buf1;
1320 DArray<Coord> grid_object_buf1, grid_normal_buf1;
1321 grid_buf1.resize(grid0_num);
1322 grid_value_buf1.resize(grid0_num);
1323 grid_object_buf1.resize(grid0_num);
1324 grid_normal_buf1.resize(grid0_num);
1325 //生成父节点及其兄弟节点
1326 cuExecute(grid0_num,
1327 SO_ComputeUpLevelGrids,
1328 grid_buf1,
1329 grid_value_buf1,
1330 grid_object_buf1,
1331 grid_normal_buf1,
1332 grid0,
1333 grid0_object,
1334 grid0_normal,
1335 m_origin,
1336 m_dx);
1337
1338 thrust::sort(thrust::device, grid_buf1.begin(), grid_buf1.begin() + grid_buf1.size(), NodeCmp<Coord>());
1339
1340 DArray<int> data_count;
1341 data_count.resize(grid0_num);
1342 data_count.reset();
1343 //数不重复的grid
1344 cuExecute(data_count.size(),
1345 SO_CountNonRepeatedGrids,
1346 data_count,
1347 grid_buf1,
1348 grid_value_buf1);
1349
1350 int grid1_num = thrust::reduce(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), (int)0, thrust::plus<int>());
1351 thrust::exclusive_scan(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), data_count.begin());
1352
1353 DArray<IndexNode> xIndex, yIndex, zIndex;
1354 xIndex.resize(grid1_num);
1355 yIndex.resize(grid1_num);
1356 zIndex.resize(grid1_num);
1357 grid1.resize(grid1_num);
1358 grid1_value.resize(grid1_num);
1359 grid1_object.resize(grid1_num);
1360 grid1_normal.resize(grid1_num);
1361 //去重
1362 cuExecute(data_count.size(),
1363 SO_FetchNonRepeatedGrids,
1364 grid1,
1365 grid1_value,
1366 grid1_object,
1367 grid1_normal,
1368 xIndex,
1369 yIndex,
1370 zIndex,
1371 grid_buf1,
1372 grid_value_buf1,
1373 grid_object_buf1,
1374 grid_normal_buf1,
1375 data_count,
1376 up_nx,
1377 up_ny,
1378 up_nz);
1379
1380 thrust::sort(thrust::device, xIndex.begin(), xIndex.begin() + xIndex.size(), IndexCmp());
1381 thrust::sort(thrust::device, yIndex.begin(), yIndex.begin() + yIndex.size(), IndexCmp());
1382 thrust::sort(thrust::device, zIndex.begin(), zIndex.begin() + zIndex.size(), IndexCmp());
1383
1384 //FIM更新
1385 for (int i = 0; i < 3; i++)
1386 {
1387 cuExecute(grid1_num,
1388 SO_FIMUpLevelGrids,
1389 grid1,
1390 grid1_value,
1391 xIndex,
1392 yIndex,
1393 zIndex,
1394 grid1_object,
1395 grid1_normal,
1396 up_nx,
1397 up_ny,
1398 up_nz);
1399 }
1400
1401 xIndex.clear();
1402 yIndex.clear();
1403 zIndex.clear();
1404 grid_buf1.clear();
1405 grid_value_buf1.clear();
1406 grid_object_buf1.clear();
1407 grid_normal_buf1.clear();
1408 data_count.clear();
1409 }
1410
1411 template<typename TDataType>
1412 void VolumeHelper<TDataType>::levelTop(
1413 DArray<VoxelOctreeNode<Coord>>& grid2,
1414 DArray<Real>& grid2_value,
1415 DArray<Coord>& grid2_object,
1416 DArray<Coord>& grid2_normal,
1417 DArray<VoxelOctreeNode<Coord>>& grid1,
1418 DArray<Coord>& grid1_object,
1419 DArray<Coord>& grid1_normal,
1420 Coord m_origin,
1421 Level multi_level,
1422 int m_nx,
1423 int m_ny,
1424 int m_nz,
1425 Real m_dx)
1426 {
1427 Real coef = Real(pow(Real(2), int(multi_level)));
1428 int up_nx = m_nx / coef;
1429 int up_ny = m_ny / coef;
1430 int up_nz = m_nz / coef;
1431 Real dx_top = m_dx * coef;
1432 int grids_num = up_nx * up_ny * up_nz;
1433
1434 grid2.resize(grids_num);
1435 grid2_value.resize(grids_num);
1436 grid2_object.resize(grids_num);
1437 grid2_normal.resize(grids_num);
1438
1439 DArray<Coord> grid_object_buf, grid_normal_buf;
1440 grid_object_buf.resize(8 * grids_num);
1441 grid_normal_buf.resize(8 * grids_num);
1442 DArray<int> data_count;
1443 data_count.resize(grids_num);
1444 data_count.reset();
1445 int grid1_num = grid1.size();
1446 //生成topside节点:主要构建父子关系
1447 cuExecute(grid1_num,
1448 SO_ComputeTopLevelGrids,
1449 grid2,
1450 grid_object_buf,
1451 grid_normal_buf,
1452 data_count,
1453 grid1,
1454 grid1_object,
1455 grid1_normal,
1456 m_origin,
1457 dx_top,
1458 up_nx,
1459 up_ny,
1460 up_nz,
1461 multi_level);
1462
1463 //更新topside节点中的值
1464 cuExecute(grids_num,
1465 SO_UpdateTopLevelGrids,
1466 grid2,
1467 grid2_value,
1468 grid2_object,
1469 grid2_normal,
1470 data_count,
1471 grid_object_buf,
1472 grid_normal_buf,
1473 m_origin,
1474 dx_top,
1475 up_nx,
1476 up_ny,
1477 up_nz,
1478 multi_level);
1479
1480 Reduction<int> reduce;
1481 DArray<int> data_count_temp;
1482 int total_num = reduce.accumulate(data_count.begin(), data_count.size());
1483 while (total_num < (grids_num))
1484 {
1485 grid_object_buf.assign(grid2_object);
1486 grid_normal_buf.assign(grid2_normal);
1487 data_count_temp.assign(data_count);
1488 //topside节点中的值FIM更新
1489 cuExecute(grids_num,
1490 SO_FIMComputeTopLevelGrids,
1491 grid2,
1492 grid2_value,
1493 grid2_object,
1494 grid2_normal,
1495 data_count,
1496 grid_object_buf,
1497 grid_normal_buf,
1498 data_count_temp,
1499 up_nx,
1500 up_ny,
1501 up_nz);
1502
1503 total_num = reduce.accumulate(data_count.begin(), data_count.size());
1504 }
1505 data_count_temp.clear();
1506 data_count.clear();
1507 grid_object_buf.clear();
1508 grid_normal_buf.clear();
1509 }
1510
1511
1512 template<typename TDataType>
1513 void VolumeHelper<TDataType>::levelCollection(
1514 DArray<VoxelOctreeNode<Coord>>& grids,
1515 DArray<Real>& grids_value,
1516 DArray<Coord>& grids_object,
1517 DArray<Coord>& grids_normal,
1518 DArray<VoxelOctreeNode<Coord>>& grid1,
1519 DArray<Real>& grid1_value,
1520 DArray<Coord>& grid1_object,
1521 DArray<Coord>& grid1_normal,
1522 int uplevel_num)
1523 {
1524 int grids_num_temp = grids.size() + grid1.size();
1525 DArray<VoxelOctreeNode<Coord>> grids_temp;
1526 DArray<Real> grids_value_temp;
1527 DArray<Coord> grids_object_temp, grids_normal_temp;
1528 grids_temp.resize(grids_num_temp);
1529 grids_value_temp.resize(grids_num_temp);
1530 grids_object_temp.resize(grids_num_temp);
1531 grids_normal_temp.resize(grids_num_temp);
1532 cuExecute(grids_num_temp,
1533 SO_CollectionGrids,
1534 grids_temp,
1535 grids_value_temp,
1536 grids_object_temp,
1537 grids_normal_temp,
1538 grids,
1539 grids_value,
1540 grids_object,
1541 grids_normal,
1542 grid1,
1543 grid1_value,
1544 grid1_object,
1545 grid1_normal,
1546 grids.size(),
1547 uplevel_num);
1548 grids.assign(grids_temp);
1549 grids_value.assign(grids_value_temp);
1550 grids_object.assign(grids_object_temp);
1551 grids_normal.assign(grids_normal_temp);
1552 grids_temp.clear();
1553 grids_value_temp.clear();
1554 grids_object_temp.clear();
1555 grids_normal_temp.clear();
1556 }
1557
1558 template <typename Real, typename Coord>
1559 __global__ void SO_CollectionGridsTwo(
1560 DArray<VoxelOctreeNode<Coord>> total_nodes,
1561 DArray<Real> total_value,
1562 DArray<Coord> total_object,
1563 DArray<Coord> total_normal,
1564 DArray<VoxelOctreeNode<Coord>> level0_nodes,
1565 DArray<Real> level0_value,
1566 DArray<Coord> level0_object,
1567 DArray<Coord> level0_normal,
1568 DArray<VoxelOctreeNode<Coord>> levelT_nodes,
1569 DArray<Real> levelT_value,
1570 DArray<Coord> levelT_object,
1571 DArray<Coord> levelT_normal,
1572 int level0_num)
1573 {
1574 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1575
1576 if (tId >= total_nodes.size()) return;
1577
1578 if (tId >= level0_num)
1579 {
1580 total_nodes[tId] = levelT_nodes[(tId - level0_num)];
1581 total_nodes[tId].setValueLocation(tId);
1582 for (int i = 0; i < 6; i++)
1583 {
1584 if (total_nodes[tId].m_neighbor[i] != EMPTY)
1585 total_nodes[tId].m_neighbor[i] += level0_num;
1586 }
1587
1588 total_value[tId] = levelT_value[(tId - level0_num)];
1589 total_object[tId] = levelT_object[(tId - level0_num)];
1590 total_normal[tId] = levelT_normal[(tId - level0_num)];
1591 }
1592 else
1593 {
1594 total_nodes[tId] = level0_nodes[tId];
1595 total_nodes[tId].setValueLocation(tId);
1596
1597 total_value[tId] = level0_value[tId];
1598 total_object[tId] = level0_object[tId];
1599 total_normal[tId] = level0_normal[tId];
1600 }
1601 }
1602
1603 template<typename TDataType>
1604 void VolumeHelper<TDataType>::collectionGridsTwo(
1605 DArray<VoxelOctreeNode<Coord>>& total_nodes,
1606 DArray<Real>& total_value,
1607 DArray<Coord>& total_object,
1608 DArray<Coord>& total_normal,
1609 DArray<VoxelOctreeNode<Coord>>& level0_nodes,
1610 DArray<Real>& level0_value,
1611 DArray<Coord>& level0_object,
1612 DArray<Coord>& level0_normal,
1613 DArray<VoxelOctreeNode<Coord>>& levelT_nodes,
1614 DArray<Real>& levelT_value,
1615 DArray<Coord>& levelT_object,
1616 DArray<Coord>& levelT_normal,
1617 int level0_num,
1618 int grid_total_num)
1619 {
1620 cuExecute(grid_total_num,
1621 SO_CollectionGridsTwo,
1622 total_nodes,
1623 total_value,
1624 total_object,
1625 total_normal,
1626 level0_nodes,
1627 level0_value,
1628 level0_object,
1629 level0_normal,
1630 levelT_nodes,
1631 levelT_value,
1632 levelT_object,
1633 levelT_normal,
1634 level0_nodes.size());
1635 }
1636
1637 template<typename TDataType>
1638 void VolumeHelper<TDataType>::collectionGridsThree(
1639 DArray<VoxelOctreeNode<Coord>>& total_nodes,
1640 DArray<Real>& total_value,
1641 DArray<Coord>& total_object,
1642 DArray<Coord>& total_normal,
1643 DArray<VoxelOctreeNode<Coord>>& level0_nodes,
1644 DArray<Real>& level0_value,
1645 DArray<Coord>& level0_object,
1646 DArray<Coord>& level0_normal,
1647 DArray<VoxelOctreeNode<Coord>>& level1_nodes,
1648 DArray<Real>& level1_value,
1649 DArray<Coord>& level1_object,
1650 DArray<Coord>& level1_normal,
1651 DArray<VoxelOctreeNode<Coord>>& levelT_nodes,
1652 DArray<Real>& levelT_value,
1653 DArray<Coord>& levelT_object,
1654 DArray<Coord>& levelT_normal,
1655 int level0_num,
1656 int level1_num,
1657 int grid_total_num)
1658 {
1659 cuExecute(grid_total_num,
1660 SO_CollectionGridsThree,
1661 total_nodes,
1662 total_value,
1663 total_object,
1664 total_normal,
1665 level0_nodes,
1666 level0_value,
1667 level0_object,
1668 level0_normal,
1669 level1_nodes,
1670 level1_value,
1671 level1_object,
1672 level1_normal,
1673 levelT_nodes,
1674 levelT_value,
1675 levelT_object,
1676 levelT_normal,
1677 level0_nodes.size(),
1678 (level0_nodes.size() + level1_nodes.size()));
1679 }
1680
1681
1682 template<typename TDataType>
1683 void VolumeHelper<TDataType>::collectionGridsFour(
1684 DArray<VoxelOctreeNode<Coord>>& total_nodes,
1685 DArray<Real>& total_value,
1686 DArray<Coord>& total_object,
1687 DArray<Coord>& total_normal,
1688 DArray<VoxelOctreeNode<Coord>>& level0_nodes,
1689 DArray<Real>& level0_value,
1690 DArray<Coord>& level0_object,
1691 DArray<Coord>& level0_normal,
1692 DArray<VoxelOctreeNode<Coord>>& level1_nodes,
1693 DArray<Real>& level1_value,
1694 DArray<Coord>& level1_object,
1695 DArray<Coord>& level1_normal,
1696 DArray<VoxelOctreeNode<Coord>>& level2_nodes,
1697 DArray<Real>& level2_value,
1698 DArray<Coord>& level2_object,
1699 DArray<Coord>& level2_normal,
1700 DArray<VoxelOctreeNode<Coord>>& levelT_nodes,
1701 DArray<Real>& levelT_value,
1702 DArray<Coord>& levelT_object,
1703 DArray<Coord>& levelT_normal,
1704 int level0_num,
1705 int level1_num,
1706 int level2_num,
1707 int grid_total_num)
1708 {
1709 cuExecute(grid_total_num,
1710 SO_CollectionGridsFour,
1711 total_nodes,
1712 total_value,
1713 total_object,
1714 total_normal,
1715 level0_nodes,
1716 level0_value,
1717 level0_object,
1718 level0_normal,
1719 level1_nodes,
1720 level1_value,
1721 level1_object,
1722 level1_normal,
1723 level2_nodes,
1724 level2_value,
1725 level2_object,
1726 level2_normal,
1727 levelT_nodes,
1728 levelT_value,
1729 levelT_object,
1730 levelT_normal,
1731 level0_nodes.size(),
1732 (level0_nodes.size() + level1_nodes.size()),
1733 (level0_nodes.size() + level1_nodes.size() + level2_nodes.size()));
1734 }
1735
1736 template<typename TDataType>
1737 void VolumeHelper<TDataType>::collectionGridsFive(
1738 DArray<VoxelOctreeNode<Coord>>& total_nodes,
1739 DArray<Real>& total_value,
1740 DArray<Coord>& total_object,
1741 DArray<Coord>& total_normal,
1742 DArray<VoxelOctreeNode<Coord>>& level0_nodes,
1743 DArray<Real>& level0_value,
1744 DArray<Coord>& level0_object,
1745 DArray<Coord>& level0_normal,
1746 DArray<VoxelOctreeNode<Coord>>& level1_nodes,
1747 DArray<Real>& level1_value,
1748 DArray<Coord>& level1_object,
1749 DArray<Coord>& level1_normal,
1750 DArray<VoxelOctreeNode<Coord>>& level2_nodes,
1751 DArray<Real>& level2_value,
1752 DArray<Coord>& level2_object,
1753 DArray<Coord>& level2_normal,
1754 DArray<VoxelOctreeNode<Coord>>& level3_nodes,
1755 DArray<Real>& level3_value,
1756 DArray<Coord>& level3_object,
1757 DArray<Coord>& level3_normal,
1758 DArray<VoxelOctreeNode<Coord>>& levelT_nodes,
1759 DArray<Real>& levelT_value,
1760 DArray<Coord>& levelT_object,
1761 DArray<Coord>& levelT_normal,
1762 int level0_num,
1763 int level1_num,
1764 int level2_num,
1765 int level3_num,
1766 int grid_total_num)
1767 {
1768 cuExecute(grid_total_num,
1769 SO_CollectionGridsFive,
1770 total_nodes,
1771 total_value,
1772 total_object,
1773 total_normal,
1774 level0_nodes,
1775 level0_value,
1776 level0_object,
1777 level0_normal,
1778 level1_nodes,
1779 level1_value,
1780 level1_object,
1781 level1_normal,
1782 level2_nodes,
1783 level2_value,
1784 level2_object,
1785 level2_normal,
1786 level3_nodes,
1787 level3_value,
1788 level3_object,
1789 level3_normal,
1790 levelT_nodes,
1791 levelT_value,
1792 levelT_object,
1793 levelT_normal,
1794 level0_nodes.size(),
1795 (level0_nodes.size() + level1_nodes.size()),
1796 (level0_nodes.size() + level1_nodes.size() + level2_nodes.size()),
1797 (level0_nodes.size() + level1_nodes.size() + level2_nodes.size() + level3_nodes.size()));
1798 }
1799
1800
1801 DYN_FUNC static void kernel3(Real& val, Real val_x)
1802 {
1803 if (std::abs(val_x) < 1)
1804 val = ((0.5*abs(val_x)*abs(val_x)*abs(val_x)) - (abs(val_x)*abs(val_x)) + (2.0f / 3.0f));
1805 else if (std::abs(val_x) < 2)
1806 val = (1.0f / 6.0f)*(2.0f - (std::abs(val_x)))*(2.0f - (std::abs(val_x)))*(2.0f - (std::abs(val_x)));
1807 else
1808 val = 0.0f;
1809 }
1810
1811 template <typename Coord>
1812 __global__ void SO_NodesInitialTopology(
1813 DArray<PositionNode> nodes,
1814 DArray<Coord> pos,
1815 DArray<VoxelOctreeNode<Coord>> nodes_a,
1816 DArray<VoxelOctreeNode<Coord>> nodes_b,
1817 int off_ax,
1818 int off_ay,
1819 int off_az,
1820 int off_bx,
1821 int off_by,
1822 int off_bz,
1823 int level0_a)
1824 {
1825 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1826
1827 if (tId >= nodes.size()) return;
1828
1829 if (tId < level0_a)
1830 {
1831 OcIndex gnx0, gny0, gnz0, gnx1, gny1, gnz1;
1832 OcKey old_key = nodes_a[tId].key();
1833 RecoverFromMortonCode(old_key, gnx0, gny0, gnz0);
1834
1835 gnx1 = gnx0 + off_ax;
1836 gny1 = gny0 + off_ay;
1837 gnz1 = gnz0 + off_az;
1838 OcKey new_key = CalculateMortonCode(gnx1, gny1, gnz1);
1839
1840 nodes[tId] = PositionNode(tId, new_key);
1841 pos[tId] = nodes_a[tId].position();
1842 }
1843 else
1844 {
1845 OcIndex gnx0, gny0, gnz0, gnx1, gny1, gnz1;
1846 OcKey old_key = nodes_b[tId - level0_a].key();
1847 RecoverFromMortonCode(old_key, gnx0, gny0, gnz0);
1848
1849 gnx1 = gnx0 + off_bx;
1850 gny1 = gny0 + off_by;
1851 gnz1 = gnz0 + off_bz;
1852 OcKey new_key = CalculateMortonCode(gnx1, gny1, gnz1);
1853
1854 nodes[tId] = PositionNode(tId, new_key);
1855 pos[tId] = nodes_b[tId - level0_a].position();
1856 }
1857 }
1858
1859 template <typename Real>
1860 __global__ void SO_CountNonRepeatedUnion(
1861 DArray<int> counter,
1862 DArray<PositionNode> nodes,
1863 DArray<Real> sdf_a,
1864 DArray<Real> sdf_b,
1865 DArray<Real> value_a,
1866 DArray<Real> value_b,
1867 int level0_a)
1868 {
1869 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1870
1871 if (tId >= counter.size()) return;
1872
1873 if (tId == 0 || (nodes[tId].position_index != nodes[tId - 1].position_index))
1874 {
1875 if (tId < counter.size() - 1 && (nodes[tId].position_index == nodes[tId + 1].position_index))
1876 {
1877 Real sdf_v1, sdf_v2;
1878 if (nodes[tId].surface_index < level0_a)
1879 {
1880 sdf_v1 = sdf_a[nodes[tId].surface_index];
1881 sdf_v2 = sdf_b[nodes[tId + 1].surface_index - level0_a];
1882 }
1883 else
1884 {
1885 sdf_v1 = sdf_b[nodes[tId].surface_index - level0_a];
1886 sdf_v2 = sdf_a[nodes[tId + 1].surface_index];
1887 }
1888
1889 if (sdf_v1 < sdf_v2)
1890 counter[tId] = 1;
1891 else
1892 counter[tId + 1] = 1;
1893 }
1894 else if (tId == counter.size() - 1 || (nodes[tId].position_index != nodes[tId + 1].position_index))
1895 {
1896 if (nodes[tId].surface_index < level0_a)
1897 {
1898 if (value_b[nodes[tId].surface_index] > 0)
1899 counter[tId] = 1;
1900 }
1901 else
1902 {
1903 if (value_a[nodes[tId].surface_index] > 0)
1904 counter[tId] = 1;
1905 }
1906 }
1907 }
1908 }
1909
1910 template <typename Real>
1911 __global__ void SO_CountNonRepeatedIntersection(
1912 DArray<int> counter,
1913 DArray<PositionNode> nodes,
1914 DArray<Real> sdf_a,
1915 DArray<Real> sdf_b,
1916 DArray<Real> value_a,
1917 DArray<Real> value_b,
1918 int level0_a)
1919 {
1920 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1921
1922 if (tId >= counter.size()) return;
1923
1924 if (tId == 0 || (nodes[tId].position_index != nodes[tId - 1].position_index))
1925 {
1926 if (tId < counter.size() - 1 && (nodes[tId].position_index == nodes[tId + 1].position_index))
1927 {
1928 Real sdf_v1, sdf_v2;
1929 if (nodes[tId].surface_index < level0_a)
1930 {
1931 sdf_v1 = sdf_a[nodes[tId].surface_index];
1932 sdf_v2 = sdf_b[nodes[tId + 1].surface_index - level0_a];
1933 }
1934 else
1935 {
1936 sdf_v1 = sdf_b[nodes[tId].surface_index - level0_a];
1937 sdf_v2 = sdf_a[nodes[tId + 1].surface_index];
1938 }
1939
1940 if (sdf_v1 > sdf_v2)
1941 counter[tId] = 1;
1942 else
1943 counter[tId + 1] = 1;
1944 }
1945 else if (tId == counter.size() - 1 || (nodes[tId].position_index != nodes[tId + 1].position_index))
1946 {
1947 if (nodes[tId].surface_index < level0_a)
1948 {
1949 if (value_b[nodes[tId].surface_index] < 0)
1950 counter[tId] = 1;
1951 }
1952 else
1953 {
1954 if (value_a[nodes[tId].surface_index] < 0)
1955 counter[tId] = 1;
1956 }
1957 }
1958 }
1959 }
1960
1961 template <typename Real>
1962 __global__ void SO_CountNonRepeatedSubtractionA(
1963 DArray<int> counter,
1964 DArray<PositionNode> nodes,
1965 DArray<Real> sdf_a,
1966 DArray<Real> sdf_b,
1967 DArray<Real> value_a,
1968 DArray<Real> value_b,
1969 int level0_a)
1970 {
1971 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
1972
1973 if (tId >= counter.size()) return;
1974
1975 if (tId == 0 || (nodes[tId].position_index != nodes[tId - 1].position_index))
1976 {
1977 if (tId < counter.size() - 1 && (nodes[tId].position_index == nodes[tId + 1].position_index))
1978 {
1979 Real sdf_v1, sdf_v2;
1980 if (nodes[tId].surface_index < level0_a)
1981 {
1982 sdf_v1 = sdf_a[nodes[tId].surface_index];
1983 sdf_v2 = sdf_b[nodes[tId + 1].surface_index - level0_a];
1984 if (sdf_v1 > -sdf_v2)
1985 counter[tId] = 1;
1986 else
1987 counter[tId + 1] = 1;
1988 }
1989 else
1990 {
1991 sdf_v1 = sdf_b[nodes[tId].surface_index - level0_a];
1992 sdf_v2 = sdf_a[nodes[tId + 1].surface_index];
1993 if (-sdf_v1 > sdf_v2)
1994 counter[tId] = 1;
1995 else
1996 counter[tId + 1] = 1;
1997 }
1998 }
1999 else if (tId == counter.size() - 1 || (nodes[tId].position_index != nodes[tId + 1].position_index))
2000 {
2001 if (nodes[tId].surface_index < level0_a)
2002 {
2003 if (value_b[nodes[tId].surface_index] > 0)
2004 counter[tId] = 1;
2005 }
2006 else
2007 {
2008 if (value_a[nodes[tId].surface_index] < 0)
2009 counter[tId] = 1;
2010 }
2011 }
2012 }
2013 }
2014
2015 template <typename Real>
2016 __global__ void SO_CountNonRepeatedSubtractionB(
2017 DArray<int> counter,
2018 DArray<PositionNode> nodes,
2019 DArray<Real> sdf_a,
2020 DArray<Real> sdf_b,
2021 DArray<Real> value_a,
2022 DArray<Real> value_b,
2023 int level0_a)
2024 {
2025 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2026
2027 if (tId >= counter.size()) return;
2028
2029 if (tId == 0 || (nodes[tId].position_index != nodes[tId - 1].position_index))
2030 {
2031 if (tId < counter.size() - 1 && (nodes[tId].position_index == nodes[tId + 1].position_index))
2032 {
2033 Real sdf_v1, sdf_v2;
2034 if (nodes[tId].surface_index < level0_a)
2035 {
2036 sdf_v1 = sdf_a[nodes[tId].surface_index];
2037 sdf_v2 = sdf_b[nodes[tId + 1].surface_index - level0_a];
2038 if (-sdf_v1 > sdf_v2)
2039 counter[tId] = 1;
2040 else
2041 counter[tId + 1] = 1;
2042 }
2043 else
2044 {
2045 sdf_v1 = sdf_b[nodes[tId].surface_index - level0_a];
2046 sdf_v2 = sdf_a[nodes[tId + 1].surface_index];
2047 if (sdf_v1 > -sdf_v2)
2048 counter[tId] = 1;
2049 else
2050 counter[tId + 1] = 1;
2051 }
2052 }
2053 else if (tId == counter.size() - 1 || (nodes[tId].position_index != nodes[tId + 1].position_index))
2054 {
2055 if (nodes[tId].surface_index < level0_a)
2056 {
2057 if (value_b[nodes[tId].surface_index] < 0)
2058 counter[tId] = 1;
2059 }
2060 else
2061 {
2062 if (value_a[nodes[tId].surface_index] > 0)
2063 counter[tId] = 1;
2064 }
2065 }
2066 }
2067 }
2068
2069 //注意counter中第i+1个元素存的是0-i元素的和
2070 template <typename Real, typename Coord>
2071 __global__ void SO_FetchNonRepeatedGridsTopology(
2072 DArray<VoxelOctreeNode<Coord>> nodes,
2073 DArray<Real> nodes_value,
2074 DArray<Coord> nodes_object,
2075 DArray<Coord> nodes_normal,
2076 DArray<IndexNode> x_index,
2077 DArray<IndexNode> y_index,
2078 DArray<IndexNode> z_index,
2079 DArray<PositionNode> all_nodes,
2080 DArray<int> counter,
2081 DArray<VoxelOctreeNode<Coord>> nodes_a,
2082 DArray<Real> nodes_value_a,
2083 DArray<Coord> nodes_object_a,
2084 DArray<Coord> nodes_normal_a,
2085 DArray<VoxelOctreeNode<Coord>> nodes_b,
2086 DArray<Real> nodes_value_b,
2087 DArray<Coord> nodes_object_b,
2088 DArray<Coord> nodes_normal_b,
2089 int level0_a,
2090 int op,
2091 Coord origin_,
2092 Real dx_,
2093 int nx_,
2094 int ny_,
2095 int nz_)
2096 {
2097 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2098
2099 if (tId >= counter.size()) return;
2100
2101 if (tId < (counter.size() - 1) && counter[tId] != counter[tId + 1])
2102 {
2103 OcIndex gnx, gny, gnz;
2104 RecoverFromMortonCode((OcKey)(all_nodes[tId].position_index), gnx, gny, gnz);
2105 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2106
2107 VoxelOctreeNode<Coord> mc((Level)0, gnx, gny, gnz, pos);
2108 mc.setValueLocation(counter[tId]);
2109
2110 if (all_nodes[tId].surface_index < level0_a)
2111 {
2112 nodes[counter[tId]] = mc;
2113 nodes_value[counter[tId]] = nodes_value_a[all_nodes[tId].surface_index];
2114 nodes_object[counter[tId]] = nodes_object_a[all_nodes[tId].surface_index];
2115 //if (op == VolumeOctreeUnion<DataType3f>::BooleanOperation::SUBTRACTION_SETB)
2116 // nodes_normal[counter[tId]] = -nodes_normal_a[all_nodes[tId].surface_index];
2117 //else
2118 nodes_normal[counter[tId]] = nodes_normal_a[all_nodes[tId].surface_index];
2119 }
2120 else
2121 {
2122 nodes[counter[tId]] = mc;
2123 nodes_value[counter[tId]] = nodes_value_b[all_nodes[tId].surface_index - level0_a];
2124 nodes_object[counter[tId]] = nodes_object_b[all_nodes[tId].surface_index - level0_a];
2125 if (op == 2)
2126 nodes_normal[counter[tId]] = -nodes_normal_b[all_nodes[tId].surface_index - level0_a];
2127 else
2128 nodes_normal[counter[tId]] = nodes_normal_b[all_nodes[tId].surface_index - level0_a];
2129 }
2130
2131 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
2132 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
2133 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
2134 }
2135 else if (tId == (counter.size() - 1) && (counter[tId] < nodes.size()))
2136 {
2137 OcIndex gnx, gny, gnz;
2138 RecoverFromMortonCode((OcKey)(all_nodes[tId].position_index), gnx, gny, gnz);
2139 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2140
2141 VoxelOctreeNode<Coord> mc((Level)0, gnx, gny, gnz, pos);
2142 mc.setValueLocation(counter[tId]);
2143
2144 if (all_nodes[tId].surface_index < level0_a)
2145 {
2146 nodes[counter[tId]] = mc;
2147 nodes_value[counter[tId]] = nodes_value_a[all_nodes[tId].surface_index];
2148 nodes_object[counter[tId]] = nodes_object_a[all_nodes[tId].surface_index];
2149 //if (op == VolumeOctreeUnion<DataType3f>::BooleanOperation::SUBTRACTION_SETB)
2150 // nodes_normal[counter[tId]] = -nodes_normal_a[all_nodes[tId].surface_index];
2151 //else
2152 nodes_normal[counter[tId]] = nodes_normal_a[all_nodes[tId].surface_index];
2153 }
2154 else
2155 {
2156 nodes[counter[tId]] = mc;
2157 nodes_value[counter[tId]] = nodes_value_b[all_nodes[tId].surface_index - level0_a];
2158 nodes_object[counter[tId]] = nodes_object_b[all_nodes[tId].surface_index - level0_a];
2159 if (op == 2)
2160 nodes_normal[counter[tId]] = -nodes_normal_b[all_nodes[tId].surface_index - level0_a];
2161 else
2162 nodes_normal[counter[tId]] = nodes_normal_b[all_nodes[tId].surface_index - level0_a];
2163 }
2164
2165 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
2166 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
2167 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
2168 }
2169 }
2170
2171 template<typename TDataType>
2172 void VolumeHelper<TDataType>::finestLevelBoolean(
2173 DArray<VoxelOctreeNode<Coord>>& grid0,
2174 DArray<Real>& grid0_value,
2175 DArray<Coord>& grid0_object,
2176 DArray<Coord>& grid0_normal,
2177 DArray<VoxelOctreeNode<Coord>>& sdfOctreeNode_a,
2178 DArray<Real>& sdfValue_a,
2179 DArray<Coord>& object_a,
2180 DArray<Coord>& normal_a,
2181 DArray<VoxelOctreeNode<Coord>>& sdfOctreeNode_b,
2182 DArray<Real>& sdfValue_b,
2183 DArray<Coord>& object_b,
2184 DArray<Coord>& normal_b,
2185 std::shared_ptr<VoxelOctree<TDataType>> sdfOctree_a,
2186 std::shared_ptr<VoxelOctree<TDataType>> sdfOctree_b,
2187 int offset_ax,
2188 int offset_ay,
2189 int offset_az,
2190 int offset_bx,
2191 int offset_by,
2192 int offset_bz,
2193 int level0_a,
2194 int level0_b,
2195 int& m_level0,
2196 Coord m_origin,
2197 Real m_dx,
2198 int m_nx,
2199 int m_ny,
2200 int m_nz,
2201 int boolean)//boolean=0:union, boolean=1:intersection, boolean=2:subtraction
2202 {
2203
2204 int level0_num = level0_a + level0_b;
2205
2206 DArray<PositionNode> grid_buf;
2207 grid_buf.resize(level0_num);
2208 DArray<Coord> grid_pos;
2209 grid_pos.resize(level0_num);
2210 //将level0放到一起
2211 cuExecute(level0_num,
2212 SO_NodesInitialTopology,
2213 grid_buf,
2214 grid_pos,
2215 sdfOctreeNode_a,
2216 sdfOctreeNode_b,
2217 offset_ax,
2218 offset_ay,
2219 offset_az,
2220 offset_bx,
2221 offset_by,
2222 offset_bz,
2223 level0_a);
2224
2225 thrust::sort(thrust::device, grid_buf.begin(), grid_buf.begin() + grid_buf.size(), PositionCmp());
2226
2227 DArray<Real> pos_value_a, pos_value_b;
2228 DArray<Coord> pos_normal_a, pos_normal_b;
2229 sdfOctree_a->getSignDistanceMLS(grid_pos, pos_value_a, pos_normal_a);
2230 sdfOctree_b->getSignDistanceMLS(grid_pos, pos_value_b, pos_normal_b);
2231 pos_normal_a.clear();
2232 pos_normal_b.clear();
2233
2234 DArray<int> data_count;
2235 data_count.resize(level0_num);
2236 data_count.reset();
2237 //数不重复的grid
2238 if (boolean == 0)
2239 {
2240 cuExecute(data_count.size(),
2241 SO_CountNonRepeatedUnion,
2242 data_count,
2243 grid_buf,
2244 sdfValue_a,
2245 sdfValue_b,
2246 pos_value_a,
2247 pos_value_b,
2248 level0_a);
2249 }
2250 else if (boolean == 1)
2251 {
2252 cuExecute(data_count.size(),
2253 SO_CountNonRepeatedIntersection,
2254 data_count,
2255 grid_buf,
2256 sdfValue_a,
2257 sdfValue_b,
2258 pos_value_a,
2259 pos_value_b,
2260 level0_a);
2261 }
2262 else if (boolean == 2)
2263 {
2264 cuExecute(data_count.size(),
2265 SO_CountNonRepeatedSubtractionA,
2266 data_count,
2267 grid_buf,
2268 sdfValue_a,
2269 sdfValue_b,
2270 pos_value_a,
2271 pos_value_b,
2272 level0_a);
2273 }
2274 int grid0_num = thrust::reduce(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), (int)0, thrust::plus<int>());
2275 thrust::exclusive_scan(thrust::device, data_count.begin(), data_count.begin() + data_count.size(), data_count.begin());
2276
2277 if (grid0_num == 0) { return; }
2278
2279 m_level0 = grid0_num;
2280
2281 DArray<IndexNode> xIndex, yIndex, zIndex;
2282 xIndex.resize(grid0_num);
2283 yIndex.resize(grid0_num);
2284 zIndex.resize(grid0_num);
2285 grid0.resize(grid0_num);
2286 grid0_value.resize(grid0_num);
2287 grid0_object.resize(grid0_num);
2288 grid0_normal.resize(grid0_num);
2289 //去重
2290 cuExecute(data_count.size(),
2291 SO_FetchNonRepeatedGridsTopology,
2292 grid0,
2293 grid0_value,
2294 grid0_object,
2295 grid0_normal,
2296 xIndex,
2297 yIndex,
2298 zIndex,
2299 grid_buf,
2300 data_count,
2301 sdfOctreeNode_a,
2302 sdfValue_a,
2303 object_a,
2304 normal_a,
2305 sdfOctreeNode_b,
2306 sdfValue_b,
2307 object_b,
2308 normal_b,
2309 level0_a,
2310 boolean,
2311 m_origin,
2312 m_dx,
2313 m_nx,
2314 m_ny,
2315 m_nz);
2316
2317 thrust::sort(thrust::device, xIndex.begin(), xIndex.begin() + xIndex.size(), IndexCmp());
2318 thrust::sort(thrust::device, yIndex.begin(), yIndex.begin() + yIndex.size(), IndexCmp());
2319 thrust::sort(thrust::device, zIndex.begin(), zIndex.begin() + zIndex.size(), IndexCmp());
2320
2321 cuExecute(grid0_num,
2322 SO_UpdateBottomNeighbors,
2323 grid0,
2324 xIndex,
2325 yIndex,
2326 zIndex,
2327 m_nx,
2328 m_ny,
2329 m_nz);
2330
2331 xIndex.clear();
2332 yIndex.clear();
2333 zIndex.clear();
2334 data_count.clear();
2335 grid_buf.clear();
2336 grid_pos.clear();
2337 pos_value_a.clear();
2338 pos_value_b.clear();
2339 }
2340
2341 template <typename Real, typename Coord>
2342 GPU_FUNC int SO_ComputeGrid(
2343 int& nx_lo,
2344 int& ny_lo,
2345 int& nz_lo,
2346 int& nx_hi,
2347 int& ny_hi,
2348 int& nz_hi,
2349 VoxelOctreeNode<Coord>& node,
2350 Coord p0,
2351 Coord p1,
2352 Coord origin_,
2353 Real dx_,
2354 int nx_,
2355 int ny_,
2356 int nz_)
2357 {
2358 nx_lo = std::floor((p0[0] - origin_[0]) / dx_);
2359 ny_lo = std::floor((p0[1] - origin_[1]) / dx_);
2360 nz_lo = std::floor((p0[2] - origin_[2]) / dx_);
2361
2362 nx_hi = std::ceil((p1[0] - origin_[0]) / dx_);
2363 ny_hi = std::ceil((p1[1] - origin_[1]) / dx_);
2364 nz_hi = std::ceil((p1[2] - origin_[2]) / dx_);
2365
2366 if ((node.m_neighbor[0] == EMPTY) && ((nx_lo % 2) != 0)) nx_lo--;
2367 if ((node.m_neighbor[2] == EMPTY) && ((ny_lo % 2) != 0)) ny_lo--;
2368 if ((node.m_neighbor[4] == EMPTY) && ((nz_lo % 2) != 0)) nz_lo--;
2369 if ((node.m_neighbor[1] == EMPTY) && ((nx_hi % 2) != 1)) nx_hi++;
2370 if ((node.m_neighbor[3] == EMPTY) && ((ny_hi % 2) != 1)) ny_hi++;
2371 if ((node.m_neighbor[5] == EMPTY) && ((nz_hi % 2) != 1)) nz_hi++;
2372
2373 return (nz_hi - nz_lo + 1) * (ny_hi - ny_lo + 1) * (nx_hi - nx_lo + 1);
2374 }
2375
2376 template <typename Real, typename Coord>
2377 __global__ void SO_CountReconstructedGrids(
2378 DArray<int> counter,
2379 DArray<VoxelOctreeNode<Coord>> nodes,
2380 Coord origin_,
2381 Real dx_,
2382 Real olddx_,
2383 int nx_,
2384 int ny_,
2385 int nz_)
2386 {
2387 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2388
2389 if (tId >= counter.size()) return;
2390
2391 Coord p0 = nodes[tId].position() - 0.5*olddx_;
2392 Coord p1 = nodes[tId].position() + 0.5*olddx_;
2393
2394 int nx_lo;
2395 int ny_lo;
2396 int nz_lo;
2397
2398 int nx_hi;
2399 int ny_hi;
2400 int nz_hi;
2401
2402 counter[tId] = SO_ComputeGrid(nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi, nodes[tId], p0, p1, origin_, dx_, nx_, ny_, nz_);
2403 }
2404
2405 template <typename Real, typename Coord>
2406 __global__ void SO_FetchReconstructedGrids(
2407 DArray<PositionNode> nodes_buf,
2408 DArray<int> counter,
2409 DArray<VoxelOctreeNode<Coord>> nodes,
2410 Coord origin_,
2411 Real dx_,
2412 Real olddx_,
2413 int nx_,
2414 int ny_,
2415 int nz_)
2416 {
2417 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2418
2419 if (tId >= counter.size()) return;
2420
2421 Coord p0 = nodes[tId].position() - 0.5*olddx_;
2422 Coord p1 = nodes[tId].position() + 0.5*olddx_;
2423
2424 int nx_lo;
2425 int ny_lo;
2426 int nz_lo;
2427
2428 int nx_hi;
2429 int ny_hi;
2430 int nz_hi;
2431
2432 int num = SO_ComputeGrid(nx_lo, ny_lo, nz_lo, nx_hi, ny_hi, nz_hi, nodes[tId], p0, p1, origin_, dx_, nx_, ny_, nz_);
2433
2434 if (num > 0)
2435 {
2436 int acc_num = 0;
2437 for (int k = nz_lo; k <= nz_hi; k++) {
2438 for (int j = ny_lo; j <= ny_hi; j++) {
2439 for (int i = nx_lo; i <= nx_hi; i++)
2440 {
2441 OcKey index = CalculateMortonCode(i, j, k);
2442 nodes_buf[counter[tId] + acc_num] = PositionNode(tId, index);
2443
2444 acc_num++;
2445 }
2446 }
2447 }
2448 }
2449 }
2450
2451 __global__ void SO_CountNoRepeatedReconstructedGrids(
2452 DArray<int> counter,
2453 DArray<PositionNode> nodes)
2454 {
2455 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2456
2457 if (tId >= counter.size()) return;
2458
2459 if (tId == 0 || nodes[tId].position_index != nodes[tId - 1].position_index)
2460 {
2461 counter[tId] = 1;
2462 }
2463 }
2464
2465 template <typename Real, typename Coord>
2466 __global__ void SO_FetchNoRepeatedReconstructedGrids(
2467 DArray<PositionNode> nodes_new,
2468 DArray<Real> nodes_sdf,
2469 DArray<Coord> nodes_object,
2470 DArray<Coord> nodes_normal,
2471 DArray<int> counter,
2472 DArray<PositionNode> nodes,
2473 DArray<VoxelOctreeNode<Coord>> nodes_old,
2474 DArray<Coord> object,
2475 DArray<Coord> normal,
2476 Coord origin_,
2477 Real dx_)
2478 {
2479 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2480
2481 if (tId >= counter.size()) return;
2482
2483 if (tId == 0 || nodes[tId].position_index != nodes[tId - 1].position_index)
2484 {
2485 OcIndex gnx, gny, gnz;
2486 RecoverFromMortonCode((OcKey)(nodes[tId].position_index), gnx, gny, gnz);
2487 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2488 Coord pos_old(nodes_old[nodes[tId].surface_index].position());
2489
2490 //直接取最近的投影点
2491 Real dist = (pos - object[nodes[tId].surface_index]).norm();
2492 int index = tId;
2493 int num = 1;
2494 while (nodes[tId + num].position_index == nodes[tId].position_index)
2495 {
2496 pos_old = nodes_old[nodes[tId + num].surface_index].position();
2497
2498 Real dist_num = (pos - object[nodes[tId + num].surface_index]).norm();
2499 if (dist_num < dist)
2500 {
2501 index = tId + num;
2502 dist = dist_num;
2503 }
2504 num++;
2505 }
2506 nodes_new[counter[tId]] = nodes[index];
2507 nodes_object[counter[tId]] = object[nodes[index].surface_index];
2508 nodes_normal[counter[tId]] = normal[nodes[index].surface_index];
2509 Real sign = (pos - nodes_object[counter[tId]]).dot(nodes_normal[counter[tId]]) < Real(0) ? Real(-1) : Real(1);
2510 nodes_sdf[counter[tId]] = sign * dist;
2511
2512
2513 ////直接求平均
2514 //Coord object_sum(object[nodes[tId].surface_index]);
2515 //Coord normal_sum(normal[nodes[tId].surface_index]);
2516 //int num = 1;
2517 //while (nodes[tId + num].position_index == nodes[tId].position_index)
2518 //{
2519 // object_sum += object[nodes[tId + num].surface_index];
2520 // normal_sum += normal[nodes[tId + num].surface_index];
2521 // num++;
2522 //}
2523 //nodes_new[counter[tId]] = nodes[tId];
2524 //nodes_object[counter[tId]] = object_sum / num;
2525 //nodes_normal[counter[tId]] = normal_sum / num;
2526 //Real sign = (pos - nodes_object[counter[tId]]).dot(nodes_normal[counter[tId]]) < Real(0) ? Real(-1) : Real(1);
2527 //Real dist = (pos - nodes_object[counter[tId]]).norm();
2528 //nodes_sdf[counter[tId]] = sign * dist;
2529
2530 ////线性插值
2531 //int nnx = 0, nny = 0, nnz = 0;
2532 //Coord object_linter, normal_linter;
2533 //Coord obj_old = object[nodes[tId].surface_index];
2534 //Coord nor_old = normal[nodes[tId].surface_index];
2535 //int num = 1;
2536 //while (nodes[tId + num].position_index == nodes[tId].position_index)
2537 //{
2538 // Coord pos_num(nodes_old[nodes[tId + num].surface_index].position());
2539 // Coord obj_num(object[nodes[tId + num].surface_index]);
2540 // Coord nor_num(normal[nodes[tId + num].surface_index]);
2541 // if (nnx == 0 && abs(pos_num[0] - pos_old[0]) > (10 * REAL_EPSILON))
2542 // {
2543 // object_linter[0] = obj_num[0] * (pos[0] - pos_old[0]) / (pos_num[0] - pos_old[0]) + obj_old[0] * (pos_num[0] - pos[0]) / (pos_num[0] - pos_old[0]);
2544 // normal_linter[0] = nor_num[0] * (pos[0] - pos_old[0]) / (pos_num[0] - pos_old[0]) + nor_old[0] * (pos_num[0] - pos[0]) / (pos_num[0] - pos_old[0]);
2545 // nnx = 1;
2546 // }
2547 // if (nny == 0 && abs(pos_num[1] - pos_old[1]) > (10 * REAL_EPSILON))
2548 // {
2549 // object_linter[1] = obj_num[1] * (pos[1] - pos_old[1]) / (pos_num[1] - pos_old[1]) + obj_old[1] * (pos_num[1] - pos[1]) / (pos_num[1] - pos_old[1]);
2550 // normal_linter[1] = nor_num[1] * (pos[1] - pos_old[1]) / (pos_num[1] - pos_old[1]) + nor_old[1] * (pos_num[1] - pos[1]) / (pos_num[1] - pos_old[1]);
2551 // nny = 1;
2552 // }
2553 // if (nnz == 0 && abs(pos_num[2] - pos_old[2]) > (10 * REAL_EPSILON))
2554 // {
2555 // object_linter[2] = obj_num[2] * (pos[2] - pos_old[2]) / (pos_num[2] - pos_old[2]) + obj_old[2] * (pos_num[2] - pos[2]) / (pos_num[2] - pos_old[2]);
2556 // normal_linter[2] = nor_num[2] * (pos[2] - pos_old[2]) / (pos_num[2] - pos_old[2]) + nor_old[2] * (pos_num[2] - pos[2]) / (pos_num[2] - pos_old[2]);
2557 // nnz = 1;
2558 // }
2559 // if (nnx == 1 && (nny == 1 && nnz == 1))
2560 // {
2561 // break;
2562 // }
2563 // num++;
2564 //}
2565 //nodes_new[counter[tId]] = nodes[tId];
2566 //nodes_object[counter[tId]] = object_linter;
2567 //nodes_normal[counter[tId]] = normal_linter;
2568 //Real sign = (pos - nodes_object[counter[tId]]).dot(nodes_normal[counter[tId]]) < Real(0) ? Real(-1) : Real(1);
2569 //Real dist = (pos - nodes_object[counter[tId]]).norm();
2570 //nodes_sdf[counter[tId]] = sign * dist;
2571
2572 //printf("Boolean neighbor: %d %d \n", tId, num);
2573 }
2574 }
2575
2576 //根据现有的最底层网格进行插值（采用权重与距离相关的非线性插值）
2577 template <typename Real, typename Coord>
2578 __global__ void SO_FetchAndInterpolateReconstructedGrids(
2579 DArray<PositionNode> nodes_new,
2580 DArray<Real> nodes_sdf,
2581 DArray<Coord> nodes_object,
2582 DArray<Coord> nodes_normal,
2583 DArray<IndexNode> x_index,
2584 DArray<IndexNode> y_index,
2585 DArray<IndexNode> z_index,
2586 DArray<int> counter,
2587 DArray<PositionNode> nodes,
2588 DArray<VoxelOctreeNode<Coord>> nodes_old,
2589 DArray<Coord> object,
2590 DArray<Coord> normal,
2591 Coord origin_,
2592 Real dx_,
2593 int nx_,
2594 int ny_,
2595 int nz_)
2596 {
2597 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2598 if (tId >= counter.size()) return;
2599
2600 if (tId == 0 || nodes[tId].position_index != nodes[tId - 1].position_index)
2601 {
2602 OcIndex gnx, gny, gnz;
2603 RecoverFromMortonCode((OcKey)(nodes[tId].position_index), gnx, gny, gnz);
2604 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2605 Coord pos_old(nodes_old[nodes[tId].surface_index].position());
2606
2607 int multiple = 1.5;
2608
2609 Real dist = (pos - pos_old).norm();
2610 Real weight_sum(1.0f), weight(1.0f);
2611 kernel3(weight_sum, dist / (multiple*dx_));
2612
2613 int num = 1;
2614 while (nodes[tId + num].position_index == nodes[tId].position_index)
2615 {
2616 pos_old = nodes_old[nodes[tId + num].surface_index].position();
2617
2618 dist = (pos - pos_old).norm();
2619 kernel3(weight, dist / (multiple*dx_));
2620
2621 weight_sum += weight;
2622
2623 num++;
2624 }
2625
2626 Coord pobject(object[nodes[tId].surface_index]);
2627 Coord pnormal(normal[nodes[tId].surface_index]);
2628
2629 if (weight_sum > 10 * REAL_EPSILON)
2630 {
2631 pos_old = nodes_old[nodes[tId].surface_index].position();
2632 dist = (pos - pos_old).norm();
2633 kernel3(weight, dist / (multiple*dx_));
2634
2635 pobject = pobject * (weight / weight_sum);
2636 pnormal = pnormal * (weight / weight_sum);
2637
2638 num = 1;
2639 while (nodes[tId + num].position_index == nodes[tId].position_index)
2640 {
2641 pos_old = nodes_old[nodes[tId + num].surface_index].position();
2642 dist = (pos - pos_old).norm();
2643 kernel3(weight, dist / (multiple*dx_));
2644 pobject += ((weight / weight_sum)*object[(nodes[tId + num].surface_index)]);
2645 pnormal += ((weight / weight_sum)*normal[(nodes[tId + num].surface_index)]);
2646
2647 num++;
2648 }
2649 }
2650 nodes_new[counter[tId]] = nodes[tId];
2651 nodes_object[counter[tId]] = pobject;
2652 nodes_normal[counter[tId]] = pnormal;
2653
2654 dist = (pos - nodes_object[counter[tId]]).norm();
2655 Real sign = (pos - nodes_object[counter[tId]]).dot(nodes_normal[counter[tId]]) < Real(0) ? Real(-1) : Real(1);
2656 nodes_sdf[counter[tId]] = sign * dist;
2657
2658 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
2659 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
2660 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
2661
2662 }
2663 }
2664
2665 //根据现有的最底层网格进行矩阵求解来计算新网格的投影点和法向（根据法线垂直于表面来求解）
2666 template <typename Real, typename Coord>
2667 __global__ void SO_FetchAndInterpolateReconstructedGrids2(
2668 DArray<PositionNode> nodes_new,
2669 DArray<Real> nodes_sdf,
2670 DArray<Coord> nodes_object,
2671 DArray<Coord> nodes_normal,
2672 DArray<IndexNode> x_index,
2673 DArray<IndexNode> y_index,
2674 DArray<IndexNode> z_index,
2675 DArray<int> counter,
2676 DArray<PositionNode> nodes,
2677 DArray<VoxelOctreeNode<Coord>> nodes_old,
2678 DArray<Coord> object,
2679 DArray<Coord> normal,
2680 Coord origin_,
2681 Real dx_,
2682 int nx_,
2683 int ny_,
2684 int nz_)
2685 {
2686 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2687 if (tId >= counter.size()) return;
2688
2689 if (tId == 0 || nodes[tId].position_index != nodes[tId - 1].position_index)
2690 {
2691 OcIndex gnx, gny, gnz;
2692 RecoverFromMortonCode((OcKey)(nodes[tId].position_index), gnx, gny, gnz);
2693 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2694 Coord pos_old(nodes_old[nodes[tId].surface_index].position());
2695
2696 int multiple = 4;
2697
2698 Real dist = (pos - pos_old).norm();
2699 Real weight_sum(1.0f), weight(1.0f);
2700 kernel3(weight_sum, dist / (multiple*dx_));
2701
2702 int num = 1;
2703 while (nodes[tId + num].position_index == nodes[tId].position_index)
2704 {
2705 pos_old = nodes_old[nodes[tId + num].surface_index].position();
2706 dist = (pos - pos_old).norm();
2707 kernel3(weight, dist / (multiple*dx_));
2708
2709 weight_sum += weight;
2710 num++;
2711 }
2712
2713 Coord pobject(object[nodes[tId].surface_index]);
2714 Coord pnormal(normal[nodes[tId].surface_index]);
2715
2716 if (weight_sum > 10 * REAL_EPSILON)
2717 {
2718 pos_old = nodes_old[nodes[tId].surface_index].position();
2719 dist = (pos - pos_old).norm();
2720 kernel3(weight, dist / (multiple*dx_));
2721 weight = weight / weight_sum;
2722 pnormal = pnormal * weight;
2723
2724 Vec3d b(pnormal[0], pnormal[1], pnormal[2]);
2725 b *= (weight*(pobject.dot(pnormal)));
2726 //b *= (pobject.dot(pnormal));
2727
2728 Mat3d M(pnormal[0] * pnormal[0], pnormal[0] * pnormal[1], pnormal[0] * pnormal[2],
2729 pnormal[1] * pnormal[0], pnormal[1] * pnormal[1], pnormal[1] * pnormal[2],
2730 pnormal[2] * pnormal[0], pnormal[2] * pnormal[1], pnormal[2] * pnormal[2]);
2731 M *= weight;
2732
2733 num = 1;
2734 while (nodes[tId + num].position_index == nodes[tId].position_index)
2735 {
2736 pos_old = nodes_old[nodes[tId + num].surface_index].position();
2737 dist = (pos - pos_old).norm();
2738 kernel3(weight, dist / (multiple*dx_));
2739 weight = weight / weight_sum;
2740
2741 Coord pobj(object[(nodes[tId + num].surface_index)]);
2742 Coord pnor(normal[(nodes[tId + num].surface_index)]);
2743 pnormal += (weight*pnor);
2744
2745 Vec3d b_i(pnor[0], pnor[1], pnor[2]);
2746 b += (weight*b_i*(pobj.dot(pnor)));
2747
2748 M(0, 0) += weight * pnor[0] * pnor[0]; M(0, 1) += weight * pnor[0] * pnor[1]; M(0, 2) += weight * pnor[0] * pnor[2];
2749 M(1, 0) += weight * pnor[1] * pnor[0]; M(1, 1) += weight * pnor[1] * pnor[1]; M(1, 2) += weight * pnor[1] * pnor[2];
2750 M(2, 0) += weight * pnor[2] * pnor[0]; M(2, 1) += weight * pnor[2] * pnor[1]; M(2, 2) += weight * pnor[2] * pnor[2];
2751
2752 num++;
2753 }
2754 Mat3d M_inv = M.inverse();
2755 Vec3d x = M_inv * b;
2756 pobject = Coord(x[0], x[1], x[2]);
2757 }
2758
2759 nodes_new[counter[tId]] = nodes[tId];
2760 nodes_object[counter[tId]] = pobject;
2761 nodes_normal[counter[tId]] = pnormal;
2762
2763 dist = (pos - nodes_object[counter[tId]]).norm();
2764 Real sign = (pos - nodes_object[counter[tId]]).dot(nodes_normal[counter[tId]]) < Real(0) ? Real(-1) : Real(1);
2765 nodes_sdf[counter[tId]] = sign * dist;
2766
2767 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
2768 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
2769 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
2770 }
2771 }
2772
2773 template <typename Real, typename Coord>
2774 GPU_FUNC void SO_UpdateReconstructionGrid(
2775 DArray<Real>& nodes_sdf,
2776 DArray<Coord>& nodes_object,
2777 DArray<Coord>& nodes_normal,
2778 Coord ppos,
2779 int id,
2780 int nid)
2781 {
2782 Real dist = (ppos - nodes_object[nid]).norm();
2783 if (abs(dist) < abs(nodes_sdf[id]))
2784 {
2785 nodes_object[id] = nodes_object[nid];
2786 nodes_normal[id] = nodes_normal[nid];
2787
2788 Real sign = (ppos - nodes_object[id]).dot(nodes_normal[id]) < Real(0) ? Real(-1) : Real(1);
2789
2790 nodes_sdf[id] = sign * dist;
2791 }
2792 }
2793
2794 template <typename Real, typename Coord>
2795 __global__ void SO_FIMUpdateReconstructionGrids(
2796 DArray<PositionNode> nodes,
2797 DArray<Real> nodes_sdf,
2798 DArray<Coord> nodes_object,
2799 DArray<Coord> nodes_normal,
2800 DArray<IndexNode> x_index,
2801 DArray<IndexNode> y_index,
2802 DArray<IndexNode> z_index,
2803 Coord origin_,
2804 Real dx_,
2805 int nx_,
2806 int ny_,
2807 int nz_)
2808 {
2809 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2810 if (tId >= nodes.size()) return;
2811
2812 OcIndex gnx, gny, gnz;
2813 RecoverFromMortonCode((OcKey)(nodes[tId].position_index), gnx, gny, gnz);
2814 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2815
2816 int xnx = (x_index[tId].xyz_index) % nx_;
2817 if ((xnx != 0 && tId > 0) && (x_index[tId - 1].xyz_index == (x_index[tId].xyz_index - 1)))
2818 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (x_index[tId].node_index), (x_index[tId - 1].node_index));
2819 if ((xnx != (nx_ - 1) && tId < (nodes.size() - 1)) && (x_index[tId + 1].xyz_index == (x_index[tId].xyz_index + 1)))
2820 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (x_index[tId].node_index), (x_index[tId + 1].node_index));
2821
2822 int yny = (y_index[tId].xyz_index) % ny_;
2823 if ((yny != 0 && tId > 0) && (y_index[tId - 1].xyz_index == (y_index[tId].xyz_index - 1)))
2824 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (y_index[tId].node_index), (y_index[tId - 1].node_index));
2825 if ((yny != (ny_ - 1) && tId < (nodes.size() - 1)) && (y_index[tId + 1].xyz_index == (y_index[tId].xyz_index + 1)))
2826 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (y_index[tId].node_index), (y_index[tId + 1].node_index));
2827
2828 int znz = (z_index[tId].xyz_index) % nz_;
2829 if ((znz != 0 && tId > 0) && (z_index[tId - 1].xyz_index == (z_index[tId].xyz_index - 1)))
2830 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (z_index[tId].node_index), (z_index[tId - 1].node_index));
2831 if ((znz != (nz_ - 1) && tId < (nodes.size() - 1)) && (z_index[tId + 1].xyz_index == (z_index[tId].xyz_index + 1)))
2832 SO_UpdateReconstructionGrid(nodes_sdf, nodes_object, nodes_normal, pos, (z_index[tId].node_index), (z_index[tId + 1].node_index));
2833 }
2834
2835 template<typename TDataType>
2836 void VolumeHelper<TDataType>::finestLevelReconstruction(
2837 DArray<PositionNode>& recon_node,
2838 DArray<Real>& recon_sdf,
2839 DArray<Coord>& recon_object,
2840 DArray<Coord>& recon_normal,
2841 DArray<VoxelOctreeNode<Coord>>& grid,
2842 DArray<Coord>& grid_object,
2843 DArray<Coord>& grid_normal,
2844 Coord m_origin,
2845 Real m_dx,
2846 int m_nx,
2847 int m_ny,
2848 int m_nz,
2849 Real m_dx_old,
2850 int level_0,
2851 int& level_0_recon)
2852 {
2853 DArray<int> count;
2854 count.resize(level_0);
2855 //数重建网格的数目(with repeat)
2856 cuExecute(count.size(),
2857 SO_CountReconstructedGrids,
2858 count,
2859 grid,
2860 m_origin,
2861 m_dx,
2862 m_dx_old,
2863 m_nx,
2864 m_ny,
2865 m_nz);
2866
2867 int grid_num = thrust::reduce(thrust::device, count.begin(), count.begin() + count.size(), (int)0, thrust::plus<int>());
2868 thrust::exclusive_scan(thrust::device, count.begin(), count.begin() + count.size(), count.begin());
2869 //std::printf("the reconstructed grids(with repeat) is: %d %d \n", sdfOctree_a->getLevel0(), grid_num);
2870
2871 DArray<PositionNode> grid_buf;
2872 grid_buf.resize(grid_num);
2873 //取出重建的网格(with repeat)
2874 cuExecute(count.size(),
2875 SO_FetchReconstructedGrids,
2876 grid_buf,
2877 count,
2878 grid,
2879 m_origin,
2880 m_dx,
2881 m_dx_old,
2882 m_nx,
2883 m_ny,
2884 m_nz);
2885
2886 thrust::sort(thrust::device, grid_buf.begin(), grid_buf.begin() + grid_buf.size(), PositionCmp());
2887
2888 count.resize(grid_num);
2889 count.reset();
2890 //数重建网格不重复的网格
2891 cuExecute(grid_num,
2892 SO_CountNoRepeatedReconstructedGrids,
2893 count,
2894 grid_buf);
2895
2896 level_0_recon = thrust::reduce(thrust::device, count.begin(), count.begin() + count.size(), (int)0, thrust::plus<int>());
2897 thrust::exclusive_scan(thrust::device, count.begin(), count.begin() + count.size(), count.begin());
2898
2899 DArray<IndexNode> xIndex, yIndex, zIndex;
2900 xIndex.resize(level_0_recon);
2901 yIndex.resize(level_0_recon);
2902 zIndex.resize(level_0_recon);
2903 recon_node.resize(level_0_recon);
2904 recon_sdf.resize(level_0_recon);
2905 recon_object.resize(level_0_recon);
2906 recon_normal.resize(level_0_recon);
2907 //取出重建网格不重复的网格
2908 //cuExecute(grid_num,
2909 // SO_FetchNoRepeatedReconstructedGrids,
2910 // recon_node,
2911 // recon_sdf,
2912 // recon_object,
2913 // recon_normal,
2914 // count,
2915 // grid_buf,
2916 // grid,
2917 // grid_object,
2918 // grid_normal,
2919 // m_origin,
2920 // m_dx);
2921 cuExecute(grid_num,
2922 SO_FetchAndInterpolateReconstructedGrids,
2923 recon_node,
2924 recon_sdf,
2925 recon_object,
2926 recon_normal,
2927 xIndex,
2928 yIndex,
2929 zIndex,
2930 count,
2931 grid_buf,
2932 grid,
2933 grid_object,
2934 grid_normal,
2935 m_origin,
2936 m_dx,
2937 m_nx,
2938 m_ny,
2939 m_nz);
2940 cuExecute(grid_num,
2941 SO_FIMUpdateReconstructionGrids,
2942 recon_node,
2943 recon_sdf,
2944 recon_object,
2945 recon_normal,
2946 xIndex,
2947 yIndex,
2948 zIndex,
2949 m_origin,
2950 m_dx,
2951 m_nx,
2952 m_ny,
2953 m_nz);
2954
2955 count.clear();
2956 grid_buf.clear();
2957 xIndex.clear();
2958 yIndex.clear();
2959 zIndex.clear();
2960 }
2961
2962 template <typename Coord>
2963 __global__ void SO_NodesInitialReconstruction(
2964 DArray<PositionNode> nodes,
2965 DArray<Coord> pos,
2966 DArray<PositionNode> nodes_a,
2967 DArray<VoxelOctreeNode<Coord>> nodes_b,
2968 int off_bx,
2969 int off_by,
2970 int off_bz,
2971 int level0_a,
2972 Coord origin_,
2973 Real dx_)
2974 {
2975 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
2976
2977 if (tId >= nodes.size()) return;
2978
2979 if (tId < level0_a)
2980 {
2981 OcIndex gnx, gny, gnz;
2982 RecoverFromMortonCode((OcKey)(nodes_a[tId].position_index), gnx, gny, gnz);
2983 Coord pos0(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
2984
2985 nodes[tId] = PositionNode(tId, nodes_a[tId].position_index);
2986 pos[tId] = pos0;
2987 }
2988 else
2989 {
2990 OcIndex gnx0, gny0, gnz0, gnx1, gny1, gnz1;
2991 OcKey old_key = nodes_b[tId - level0_a].key();
2992 RecoverFromMortonCode(old_key, gnx0, gny0, gnz0);
2993
2994 gnx1 = gnx0 + off_bx;
2995 gny1 = gny0 + off_by;
2996 gnz1 = gnz0 + off_bz;
2997 OcKey new_key = CalculateMortonCode(gnx1, gny1, gnz1);
2998
2999 nodes[tId] = PositionNode(tId, new_key);
3000 pos[tId] = nodes_b[tId - level0_a].position();
3001 }
3002 }
3003
3004 //注意counter中第i+1个元素存的是0-i元素的和
3005 template <typename Real, typename Coord>
3006 __global__ void SO_FetchNonRepeatedGridsReconstruction(
3007 DArray<VoxelOctreeNode<Coord>> nodes,
3008 DArray<Real> nodes_value,
3009 DArray<Coord> nodes_object,
3010 DArray<Coord> nodes_normal,
3011 DArray<IndexNode> x_index,
3012 DArray<IndexNode> y_index,
3013 DArray<IndexNode> z_index,
3014 DArray<PositionNode> all_nodes,
3015 DArray<int> counter,
3016 DArray<Real> nodes_value_a,
3017 DArray<Coord> nodes_object_a,
3018 DArray<Coord> nodes_normal_a,
3019 DArray<VoxelOctreeNode<Coord>> nodes_b,
3020 DArray<Real> nodes_value_b,
3021 DArray<Coord> nodes_object_b,
3022 DArray<Coord> nodes_normal_b,
3023 int level0_a,
3024 int op,
3025 Coord origin_,
3026 Real dx_,
3027 int nx_,
3028 int ny_,
3029 int nz_)
3030 {
3031 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
3032
3033 if (tId >= counter.size()) return;
3034
3035 if (tId < (counter.size() - 1) && counter[tId] != counter[tId + 1])
3036 {
3037 OcIndex gnx, gny, gnz;
3038 RecoverFromMortonCode((OcKey)(all_nodes[tId].position_index), gnx, gny, gnz);
3039 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
3040
3041 VoxelOctreeNode<Coord> mc((Level)0, gnx, gny, gnz, pos);
3042 mc.setValueLocation(counter[tId]);
3043
3044 if (all_nodes[tId].surface_index < level0_a)
3045 {
3046 nodes[counter[tId]] = mc;
3047 nodes_value[counter[tId]] = nodes_value_a[all_nodes[tId].surface_index];
3048 nodes_object[counter[tId]] = nodes_object_a[all_nodes[tId].surface_index];
3049 if (op == 3)
3050 nodes_normal[counter[tId]] = -nodes_normal_a[all_nodes[tId].surface_index];
3051 else
3052 nodes_normal[counter[tId]] = nodes_normal_a[all_nodes[tId].surface_index];
3053 }
3054 else
3055 {
3056 nodes[counter[tId]] = mc;
3057 nodes_value[counter[tId]] = nodes_value_b[all_nodes[tId].surface_index - level0_a];
3058 nodes_object[counter[tId]] = nodes_object_b[all_nodes[tId].surface_index - level0_a];
3059 if (op == 2)
3060 nodes_normal[counter[tId]] = -nodes_normal_b[all_nodes[tId].surface_index - level0_a];
3061 else
3062 nodes_normal[counter[tId]] = nodes_normal_b[all_nodes[tId].surface_index - level0_a];
3063 }
3064
3065 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
3066 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
3067 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
3068 }
3069 else if (tId == (counter.size() - 1) && (counter[tId] < nodes.size()))
3070 {
3071 OcIndex gnx, gny, gnz;
3072 RecoverFromMortonCode((OcKey)(all_nodes[tId].position_index), gnx, gny, gnz);
3073 Coord pos(origin_[0] + (gnx + 0.5) * dx_, origin_[1] + (gny + 0.5) * dx_, origin_[2] + (gnz + 0.5) * dx_);
3074
3075 VoxelOctreeNode<Coord> mc((Level)0, gnx, gny, gnz, pos);
3076 mc.setValueLocation(counter[tId]);
3077
3078 if (all_nodes[tId].surface_index < level0_a)
3079 {
3080 nodes[counter[tId]] = mc;
3081 nodes_value[counter[tId]] = nodes_value_a[all_nodes[tId].surface_index];
3082 nodes_object[counter[tId]] = nodes_object_a[all_nodes[tId].surface_index];
3083 if (op == 3)
3084 nodes_normal[counter[tId]] = -nodes_normal_a[all_nodes[tId].surface_index];
3085 else
3086 nodes_normal[counter[tId]] = nodes_normal_a[all_nodes[tId].surface_index];
3087 }
3088 else
3089 {
3090 nodes[counter[tId]] = mc;
3091 nodes_value[counter[tId]] = nodes_value_b[all_nodes[tId].surface_index - level0_a];
3092 nodes_object[counter[tId]] = nodes_object_b[all_nodes[tId].surface_index - level0_a];
3093 if (op == 2)
3094 nodes_normal[counter[tId]] = -nodes_normal_b[all_nodes[tId].surface_index - level0_a];
3095 else
3096 nodes_normal[counter[tId]] = nodes_normal_b[all_nodes[tId].surface_index - level0_a];
3097 }
3098
3099 x_index[counter[tId]] = IndexNode((gnz*nx_*ny_ + gny * nx_ + gnx), counter[tId]);
3100 y_index[counter[tId]] = IndexNode((gnx*ny_*nz_ + gnz * ny_ + gny), counter[tId]);
3101 z_index[counter[tId]] = IndexNode((gny*nz_*nx_ + gnx * nz_ + gnz), counter[tId]);
3102 }
3103 }
3104
3105
3106 template <typename Real, typename Coord>
3107 __global__ void SO_FIMUpdateGrids(
3108 DArray<VoxelOctreeNode<Coord>> nodes,
3109 DArray<Real> nodes_value,
3110 DArray<Coord> nodes_object,
3111 DArray<Coord> nodes_normal)
3112 {
3113 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
3114 if (tId >= nodes.size()) return;
3115
3116 for (int i = 0; i < 6; i++)
3117 {
3118 if (nodes[tId].m_neighbor[i] != EMPTY)
3119 {
3120 int nb_id = nodes[tId].m_neighbor[i];
3121 Real sign = (nodes[tId].position() - nodes_object[nb_id]).dot(nodes_normal[nb_id]) < Real(0) ? Real(-1) : Real(1);
3122 Real dist = sign * (nodes[tId].position() - nodes_object[nb_id]).norm();
3123
3124 if (abs(dist) < abs(nodes_value[tId]))
3125 {
3126 nodes_value[tId] = dist;
3127 nodes_object[tId] = nodes_object[nb_id];
3128 nodes_normal[tId] = nodes_normal[nb_id];
3129 }
3130 }
3131 }
3132 }
3133
3134 template<typename TDataType>
3135 void VolumeHelper<TDataType>::finestLevelReconstBoolean(
3136 DArray<VoxelOctreeNode<Coord>>& grid0,
3137 DArray<Real>& grid0_value,
3138 DArray<Coord>& grid0_object,
3139 DArray<Coord>& grid0_normal,
3140 DArray<PositionNode>& recon_node,
3141 DArray<Real>& recon_sdf,
3142 DArray<Coord>& recon_object,
3143 DArray<Coord>& recon_normal,
3144 DArray<VoxelOctreeNode<Coord>>& sdfOctreeNode_2,
3145 DArray<Real>& sdfValue_2,
3146 DArray<Coord>& object_2,
3147 DArray<Coord>& normal_2,
3148 std::shared_ptr<VoxelOctree<TDataType>> sdfOctree_1,
3149 std::shared_ptr<VoxelOctree<TDataType>> sdfOctree_2,
3150 int level0_1,
3151 int level0_2,
3152 int& m_level0,
3153 int offset_nx,
3154 int offset_ny,
3155 int offset_nz,
3156 Coord m_origin,
3157 Real m_dx,
3158 int m_nx,
3159 int m_ny,
3160 int m_nz,
3161 int boolean)//boolean=0:union, boolean=1:intersection, boolean=2:subtraction, boolean=3:subtraction with model b reconstruction
3162 {
3163 int level0_num = level0_1 + level0_2;
3164 DArray<PositionNode> grid_buf;
3165 grid_buf.resize(level0_num);
3166 DArray<Coord> grid_pos;
3167 grid_pos.resize(level0_num);
3168 //重建网格中不重复的网格与level0的网格放在一起
3169 cuExecute(level0_num,
3170 SO_NodesInitialReconstruction,
3171 grid_buf,
3172 grid_pos,
3173 recon_node,
3174 sdfOctreeNode_2,
3175 offset_nx,
3176 offset_ny,
3177 offset_nz,
3178 level0_1,
3179 m_origin,
3180 m_dx);
3181
3182 thrust::sort(thrust::device, grid_buf.begin(), grid_buf.begin() + grid_buf.size(), PositionCmp());
3183
3184 DArray<Real> pos_value_1, pos_value_2;
3185 DArray<Coord> pos_normal_1, pos_normal_2;
3186 sdfOctree_1->getSignDistanceMLS(grid_pos, pos_value_1, pos_normal_1);
3187 sdfOctree_2->getSignDistanceMLS(grid_pos, pos_value_2, pos_normal_2);
3188 pos_normal_1.clear();
3189 pos_normal_2.clear();
3190
3191 DArray<int> count;
3192 count.resize(level0_num);
3193 count.reset();
3194 //数不重复的grid
3195 if (boolean == 0)
3196 {
3197 cuExecute(count.size(),
3198 SO_CountNonRepeatedUnion,
3199 count,
3200 grid_buf,
3201 recon_sdf,
3202 sdfValue_2,
3203 pos_value_1,
3204 pos_value_2,
3205 level0_1);
3206 }
3207 else if (boolean == 1)
3208 {
3209 cuExecute(count.size(),
3210 SO_CountNonRepeatedIntersection,
3211 count,
3212 grid_buf,
3213 recon_sdf,
3214 sdfValue_2,
3215 pos_value_1,
3216 pos_value_2,
3217 level0_1);
3218 }
3219 else if (boolean == 2)
3220 {
3221 cuExecute(count.size(),
3222 SO_CountNonRepeatedSubtractionA,
3223 count,
3224 grid_buf,
3225 recon_sdf,
3226 sdfValue_2,
3227 pos_value_1,
3228 pos_value_2,
3229 level0_1);
3230 }
3231 else if (boolean == 3)
3232 {
3233 cuExecute(count.size(),
3234 SO_CountNonRepeatedSubtractionB,
3235 count,
3236 grid_buf,
3237 recon_sdf,
3238 sdfValue_2,
3239 pos_value_1,
3240 pos_value_2,
3241 level0_1);
3242 }
3243 int grid0_num = thrust::reduce(thrust::device, count.begin(), count.begin() + count.size(), (int)0, thrust::plus<int>());
3244 thrust::exclusive_scan(thrust::device, count.begin(), count.begin() + count.size(), count.begin());
3245
3246 if (grid0_num == 0) { return; }
3247
3248 m_level0 = grid0_num;
3249
3250 DArray<IndexNode> xIndex, yIndex, zIndex;
3251 xIndex.resize(grid0_num);
3252 yIndex.resize(grid0_num);
3253 zIndex.resize(grid0_num);
3254 grid0.resize(grid0_num);
3255 grid0_value.resize(grid0_num);
3256 grid0_object.resize(grid0_num);
3257 grid0_normal.resize(grid0_num);
3258 //去重
3259 cuExecute(count.size(),
3260 SO_FetchNonRepeatedGridsReconstruction,
3261 grid0,
3262 grid0_value,
3263 grid0_object,
3264 grid0_normal,
3265 xIndex,
3266 yIndex,
3267 zIndex,
3268 grid_buf,
3269 count,
3270 recon_sdf,
3271 recon_object,
3272 recon_normal,
3273 sdfOctreeNode_2,
3274 sdfValue_2,
3275 object_2,
3276 normal_2,
3277 level0_1,
3278 boolean,
3279 m_origin,
3280 m_dx,
3281 m_nx,
3282 m_ny,
3283 m_nz);
3284
3285 thrust::sort(thrust::device, xIndex.begin(), xIndex.begin() + xIndex.size(), IndexCmp());
3286 thrust::sort(thrust::device, yIndex.begin(), yIndex.begin() + yIndex.size(), IndexCmp());
3287 thrust::sort(thrust::device, zIndex.begin(), zIndex.begin() + zIndex.size(), IndexCmp());
3288
3289 cuExecute(grid0_num,
3290 SO_UpdateBottomNeighbors,
3291 grid0,
3292 xIndex,
3293 yIndex,
3294 zIndex,
3295 m_nx,
3296 m_ny,
3297 m_nz);
3298
3299 for (int i = 0; i < 3; i++)
3300 {
3301 cuExecute(grid0_num,
3302 SO_FIMUpdateGrids,
3303 grid0,
3304 grid0_value,
3305 grid0_object,
3306 grid0_normal);
3307 }
3308
3309 xIndex.clear();
3310 yIndex.clear();
3311 zIndex.clear();
3312 count.clear();
3313 grid_buf.clear();
3314 grid_pos.clear();
3315 pos_value_1.clear();
3316 pos_value_2.clear();
3317 }
3318
3319 template <typename Real>
3320 __global__ void SO_IntersectionSet(
3321 DArray<Real> inter_value,
3322 DArray<int> inter_index,
3323 DArray<Real> inter_value_a,
3324 DArray<Real> inter_value_b)
3325 {
3326 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
3327
3328 if (tId >= inter_index.size()) return;
3329
3330 if (inter_value_a[tId] < 0 && inter_value_b[tId] < 0)
3331 inter_value[inter_index[tId]] = -std::abs(inter_value[inter_index[tId]]);
3332 else
3333 inter_value[inter_index[tId]] = std::abs(inter_value[inter_index[tId]]);
3334 }
3335
3336 template <typename Real>
3337 __global__ void SO_UnionSet(
3338 DArray<Real> inter_value,
3339 DArray<int> inter_index,
3340 DArray<Real> inter_value_a,
3341 DArray<Real> inter_value_b)
3342 {
3343 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
3344
3345 if (tId >= inter_index.size()) return;
3346
3347 if (inter_value_a[tId] < 0 || inter_value_b[tId] < 0)
3348 inter_value[inter_index[tId]] = -std::abs(inter_value[inter_index[tId]]);
3349 else
3350 inter_value[inter_index[tId]] = std::abs(inter_value[inter_index[tId]]);
3351 }
3352
3353 template <typename Real>
3354 __global__ void SO_SubtractionSetA(
3355 DArray<Real> inter_value,
3356 DArray<int> inter_index,
3357 DArray<Real> inter_value_a,
3358 DArray<Real> inter_value_b)
3359 {
3360 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
3361
3362 if (tId >= inter_index.size()) return;
3363
3364 if (inter_value_a[tId] <= 0 && inter_value_b[tId] >= 0)
3365 inter_value[inter_index[tId]] = -std::abs(inter_value[inter_index[tId]]);
3366 else
3367 inter_value[inter_index[tId]] = std::abs(inter_value[inter_index[tId]]);
3368 }
3369
3370 template<typename TDataType>
3371 void VolumeHelper<TDataType>::updateBooleanSigned(
3372 DArray<Real>& leaf_value,
3373 DArray<int>& leaf_index,
3374 DArray<Real>& leaf_value_a,
3375 DArray<Real>& leaf_value_b,
3376 int boolean) //boolean = 0:union, boolean = 1 : intersection, boolean = 2 : subtraction
3377 {
3378 if (boolean == 0)
3379 {
3380 cuExecute(leaf_index.size(),
3381 SO_UnionSet,
3382 leaf_value,
3383 leaf_index,
3384 leaf_value_a,
3385 leaf_value_b);
3386
3387 }
3388 else if (boolean == 1)
3389 {
3390 cuExecute(leaf_index.size(),
3391 SO_IntersectionSet,
3392 leaf_value,
3393 leaf_index,
3394 leaf_value_a,
3395 leaf_value_b);
3396 }
3397 else if (boolean == 2)
3398 {
3399 cuExecute(leaf_index.size(),
3400 SO_SubtractionSetA,
3401 leaf_value,
3402 leaf_index,
3403 leaf_value_a,
3404 leaf_value_b);
3405 }
3406 }
3407 DEFINE_CLASS(VolumeHelper);
3408}
Generated by doxygen 1.13.2