PeriDyno 1.0.0
Loading...
Searching...
No Matches
FastSweepingMethodGPU.cu
Go to the documentation of this file.
1#include "FastSweepingMethodGPU.h"
2
3#include "Algorithm/Reduction.h"
4
5#include "Collision/Distance3D.h"
6
7namespace dyno
8{
9 IMPLEMENT_TCLASS(FastSweepingMethodGPU, TDataType)
10
11 template<typename TDataType>
12 FastSweepingMethodGPU<TDataType>::FastSweepingMethodGPU()
13 : ComputeModule()
14 {
15 }
16
17 template<typename TDataType>
18 FastSweepingMethodGPU<TDataType>::~FastSweepingMethodGPU()
19 {
20 }
21
22 enum GridType
23 {
24 Accepted = 0,
25 Tentative = 1,
26 Infinite = 2
27 };
28
29 template<typename TDataType>
30 void FastSweepingMethodGPU<TDataType>::compute()
31 {
32 this->makeLevelSet();
33 }
34
35#define TRI_MIN(x, y, z) glm::min(x, glm::min(y, z))
36#define TRI_MAX(x, y, z) glm::max(x, glm::max(y, z))
37
38 template<typename Real>
39 __global__ void FSM_AssignInifity(
40 DArray3D<Real> phi,
41 DArray3D<int> closest_tri_id,
42 DArray3D<GridType> types)
43 {
44 int i = threadIdx.x + (blockIdx.x * blockDim.x);
45 int j = threadIdx.y + (blockIdx.y * blockDim.y);
46 int k = threadIdx.z + (blockIdx.z * blockDim.z);
47
48 uint nx = phi.nx();
49 uint ny = phi.ny();
50 uint nz = phi.nz();
51
52 if (i >= nx || j >= ny || k >= nz) return;
53
54 phi(i, j, k) = REAL_MAX;
55 closest_tri_id(i, j, k) = -1;
56 types(i, j, k) = GridType::Infinite;
57 }
58
59 template<typename Real, typename Coord>
60 __global__ void FSM_FindNeighboringGrids(
61 DArray3D<uint> counter,
62 DArray<Coord> vertices,
63 DArray<TopologyModule::Triangle> indices,
64 Coord origin,
65 Real dx)
66 {
67 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
68
69 if (tId >= indices.size()) return;
70
71 int ni = counter.nx();
72 int nj = counter.ny();
73 int nk = counter.nz();
74
75 auto t = indices[tId];
76
77 Coord vp = vertices[t[0]];
78 Coord vq = vertices[t[1]];
79 Coord vr = vertices[t[2]];
80
81 // coordinates in grid to high precision
82 Real fip = (vp[0] - origin[0]) / dx;
83 Real fjp = (vp[1] - origin[1]) / dx;
84 Real fkp = (vp[2] - origin[2]) / dx;
85 Real fiq = (vq[0] - origin[0]) / dx;
86 Real fjq = (vq[1] - origin[1]) / dx;
87 Real fkq = (vq[2] - origin[2]) / dx;
88 Real fir = (vr[0] - origin[0]) / dx;
89 Real fjr = (vr[1] - origin[1]) / dx;
90 Real fkr = (vr[2] - origin[2]) / dx;
91
92 // do distances nearby
93 const int exact_band = 1;
94
95 int i0 = glm::clamp(int(TRI_MIN(fip, fiq, fir)) - exact_band, 0, ni - 1);
96 int i1 = glm::clamp(int(TRI_MAX(fip, fiq, fir)) + exact_band + 1, 0, ni - 1);
97 int j0 = glm::clamp(int(TRI_MIN(fjp, fjq, fjr)) - exact_band, 0, nj - 1);
98 int j1 = glm::clamp(int(TRI_MAX(fjp, fjq, fjr)) + exact_band + 1, 0, nj - 1);
99 int k0 = glm::clamp(int(TRI_MIN(fkp, fkq, fkr)) - exact_band, 0, nk - 1);
100 int k1 = glm::clamp(int(TRI_MAX(fkp, fkq, fkr)) + exact_band + 1, 0, nk - 1);
101
102 TTriangle3D<Real> tri(vp, vq, vr);
103 for (int k = k0; k <= k1; ++k) for (int j = j0; j <= j1; ++j) for (int i = i0; i <= i1; ++i) {
104 TPoint3D<Real> point(origin + Coord(i * dx, j * dx, k * dx));
105 if (glm::abs(point.distance(tri)) < dx * 1.74) //1.74 is chosen to be slightly larger than sqrt(3)
106 {
107 atomicAdd(&counter(i, j, k), 1);
108 }
109 }
110 }
111
112 template<typename Real, typename Coord>
113 __global__ void FSM_StoreTriangleIds(
114 DArrayList<uint> triIds,
115 DArray3D<uint> counter,
116 DArray<Coord> vertices,
117 DArray<TopologyModule::Triangle> indices,
118 Coord origin,
119 Real dx)
120 {
121 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
122
123 if (tId >= indices.size()) return;
124
125 int ni = counter.nx();
126 int nj = counter.ny();
127 int nk = counter.nz();
128
129 auto t = indices[tId];
130
131 Coord vp = vertices[t[0]];
132 Coord vq = vertices[t[1]];
133 Coord vr = vertices[t[2]];
134
135 // coordinates in grid to high precision
136 Real fip = (vp[0] - origin[0]) / dx;
137 Real fjp = (vp[1] - origin[1]) / dx;
138 Real fkp = (vp[2] - origin[2]) / dx;
139 Real fiq = (vq[0] - origin[0]) / dx;
140 Real fjq = (vq[1] - origin[1]) / dx;
141 Real fkq = (vq[2] - origin[2]) / dx;
142 Real fir = (vr[0] - origin[0]) / dx;
143 Real fjr = (vr[1] - origin[1]) / dx;
144 Real fkr = (vr[2] - origin[2]) / dx;
145
146 // do distances nearby
147 const int exact_band = 1;
148
149 int i0 = glm::clamp(int(TRI_MIN(fip, fiq, fir)) - exact_band, 0, ni - 1);
150 int i1 = glm::clamp(int(TRI_MAX(fip, fiq, fir)) + exact_band + 1, 0, ni - 1);
151 int j0 = glm::clamp(int(TRI_MIN(fjp, fjq, fjr)) - exact_band, 0, nj - 1);
152 int j1 = glm::clamp(int(TRI_MAX(fjp, fjq, fjr)) + exact_band + 1, 0, nj - 1);
153 int k0 = glm::clamp(int(TRI_MIN(fkp, fkq, fkr)) - exact_band, 0, nk - 1);
154 int k1 = glm::clamp(int(TRI_MAX(fkp, fkq, fkr)) + exact_band + 1, 0, nk - 1);
155
156 TTriangle3D<Real> tri(vp, vq, vr);
157 for (int k = k0; k <= k1; ++k) for (int j = j0; j <= j1; ++j) for (int i = i0; i <= i1; ++i) {
158 TPoint3D<Real> point(origin + Coord(i * dx, j * dx, k * dx));
159 if (glm::abs(point.distance(tri)) < dx * 1.74) //1.74 is chosen to be slightly larger than sqrt(3)
160 {
161 triIds[counter.index(i, j, k)].atomicInsert(tId);
162 }
163 }
164 }
165
166 __global__ void FSM_Array3D_To_Array1d(
167 DArray<uint> arr1d,
168 DArray3D<uint> arr3d)
169 {
170 int i = threadIdx.x + (blockIdx.x * blockDim.x);
171 int j = threadIdx.y + (blockIdx.y * blockDim.y);
172 int k = threadIdx.z + (blockIdx.z * blockDim.z);
173
174 uint nx = arr3d.nx();
175 uint ny = arr3d.ny();
176 uint nz = arr3d.nz();
177
178 if (i >= nx || j >= ny || k >= nz) return;
179
180 uint index1d = arr3d.index(i, j, k);
181
182 arr1d[index1d] = arr3d(i, j, k);
183 }
184
185 template<typename Real, typename Coord, typename Tri2Edg>
186 __global__ void FSM_InitializeSDFNearMesh(
187 DArray3D<Real> phi,
188 DArray3D<int> closestId,
189 DArray3D<GridType> gridType,
190 DArrayList<uint> triIds,
191 DArray<Coord> vertices,
192 DArray<TopologyModule::Triangle> indices,
193 DArray<TopologyModule::Edge> edges,
194 DArray<Tri2Edg> t2e,
195 DArray<Coord> edgeN,
196 DArray<Coord> vertexN,
197 Coord origin,
198 Real dx)
199 {
200 int i = threadIdx.x + (blockIdx.x * blockDim.x);
201 int j = threadIdx.y + (blockIdx.y * blockDim.y);
202 int k = threadIdx.z + (blockIdx.z * blockDim.z);
203
204 uint nx = phi.nx();
205 uint ny = phi.ny();
206 uint nz = phi.nz();
207
208 if (i >= nx || j >= ny || k >= nz) return;
209
210 Coord gx(i * dx + origin[0], j * dx + origin[1], k * dx + origin[2]);
211
212 auto& list = triIds[phi.index(i, j, k)];
213
214 if (list.size() > 0)
215 {
216 ProjectedPoint3D<Real> p3d;
217
218 bool valid = calculateSignedDistance2TriangleSetFromNormal(p3d, gx, vertices, edges, indices, t2e, edgeN, vertexN, list);
219 if (valid) {
220 phi(i, j, k) = p3d.signed_distance;
221 closestId(i, j, k) = p3d.id;
222 gridType(i, j, k) = GridType::Accepted;
223 }
224 }
225 }
226
227 template<typename Real, typename Coord>
228 GPU_FUNC void FSM_CheckNeighbour(
229 DArray3D<Real>& phi,
230 DArray3D<int>& closest_tri,
231 const DArray<Coord>& vert,
232 const DArray<TopologyModule::Triangle>& tri,
233 const Coord gx,
234 int i0, int j0, int k0,
235 int i1, int j1, int k1)
236 {
237 if (closest_tri(i1, j1, k1) >= 0) {
238 unsigned int p, q, r;
239 auto trijk = tri[closest_tri(i1, j1, k1)];
240 p = trijk[0];
241 q = trijk[1];
242 r = trijk[2];
243
244 auto t = TTriangle3D<Real>(vert[p], vert[q], vert[r]);
245 Real d = TPoint3D<Real>(gx).distance(t);
246
247 Real phi0 = phi(i0, j0, k0);
248 Real phi1 = phi(i1, j1, k1);
249
250 if (glm::abs(d) < glm::abs(phi0)) {
251 phi(i0, j0, k0) = phi1 > 0 ? glm::abs(d) : -glm::abs(d);
252 closest_tri(i0, j0, k0) = closest_tri(i1, j1, k1);
253 }
254 }
255 }
256
257 template<typename Real, typename Coord>
258 __global__ void FSM_SweepX(
259 DArray3D<Real> phi,
260 DArray3D<int> closestId,
261 DArray3D<GridType> gridType,
262 DArray<Coord> vertices,
263 DArray<TopologyModule::Triangle> indices,
264 Coord origin,
265 Real dx,
266 int di)
267 {
268 int j = threadIdx.x + (blockIdx.x * blockDim.x);
269 int k = threadIdx.y + (blockIdx.y * blockDim.y);
270
271 uint ny = phi.ny();
272 uint nz = phi.nz();
273
274 if (j >= ny || k >= nz) return;
275
276 int i0 = di > 0 ? 1 : phi.nx() - 2;
277 int i1 = di > 0 ? phi.nx() : -1;
278
279 for (int i = i0; i != i1; i += di) {
280 if (gridType(i, j, k) != GridType::Accepted)
281 {
282 Coord gx(i * dx + origin[0], j * dx + origin[1], k * dx + origin[2]);
283 FSM_CheckNeighbour(phi, closestId, vertices, indices, gx, i, j, k, i - di, j, k);
284 }
285 }
286 }
287
288 template<typename Real, typename Coord>
289 __global__ void FSM_SweepY(
290 DArray3D<Real> phi,
291 DArray3D<int> closestId,
292 DArray3D<GridType> gridType,
293 DArray<Coord> vertices,
294 DArray<TopologyModule::Triangle> indices,
295 Coord origin,
296 Real dx,
297 int dj)
298 {
299 int i = threadIdx.x + (blockIdx.x * blockDim.x);
300 int k = threadIdx.y + (blockIdx.y * blockDim.y);
301
302 uint nx = phi.nx();
303 uint nz = phi.nz();
304
305 if (i >= nx || k >= nz) return;
306
307 int j0 = dj > 0 ? 1 : phi.ny() - 2;
308 int j1 = dj > 0 ? phi.ny() : -1;
309
310 for (int j = j0; j != j1; j += dj) {
311 if (gridType(i, j, k) != GridType::Accepted)
312 {
313 Coord gx(i * dx + origin[0], j * dx + origin[1], k * dx + origin[2]);
314 FSM_CheckNeighbour(phi, closestId, vertices, indices, gx, i, j, k, i, j - dj, k);
315 }
316 }
317 }
318
319 template<typename Real, typename Coord>
320 __global__ void FSM_SweepZ(
321 DArray3D<Real> phi,
322 DArray3D<int> closestId,
323 DArray3D<GridType> gridType,
324 DArray<Coord> vertices,
325 DArray<TopologyModule::Triangle> indices,
326 Coord origin,
327 Real dx,
328 int dk)
329 {
330 int i = threadIdx.x + (blockIdx.x * blockDim.x);
331 int j = threadIdx.y + (blockIdx.y * blockDim.y);
332
333 uint nx = phi.nx();
334 uint ny = phi.ny();
335
336 if (i >= nx || j >= ny) return;
337
338 int k0 = dk > 0 ? 1 : phi.nz() - 2;
339 int k1 = dk > 0 ? phi.nz() : -1;
340
341 for (int k = k0; k != k1; k += dk) {
342 if (gridType(i, j, k) != GridType::Accepted)
343 {
344 Coord gx(i * dx + origin[0], j * dx + origin[1], k * dx + origin[2]);
345 FSM_CheckNeighbour(phi, closestId, vertices, indices, gx, i, j, k, i, j, k - dk);
346 }
347 }
348 }
349
350 template<typename TDataType>
351 void FastSweepingMethodGPU<TDataType>::makeLevelSet()
352 {
353 if (this->outLevelSet()->isEmpty()) {
354 this->outLevelSet()->allocate();
355 }
356
357 DistanceField3D<TDataType>& sdf = this->outLevelSet()->getDataPtr()->getSDF();
358
359 auto ts = this->inTriangleSet()->constDataPtr();
360
361 DArray<Coord> edgeNormal, vertexNormal;
362 ts->updateEdgeNormal(edgeNormal);
363 ts->updateAngleWeightedVertexNormal(vertexNormal);
364
365 ts->updateTriangle2Edge();
366 auto& tri2edg = ts->getTriangle2Edge();
367
368 auto& vertices = ts->getPoints();
369 auto& edges = ts->getEdges();
370 auto& triangles = ts->getTriangles();
371
372 Reduction<Coord> reduce;
373 Coord min_box = reduce.minimum(vertices.begin(), vertices.size());
374 Coord max_box = reduce.maximum(vertices.begin(), vertices.size());
375
376 Real dx = this->varSpacing()->getValue();
377 uint padding = this->varPadding()->getValue();
378
379 Coord unit(1, 1, 1);
380 min_box -= padding * dx * unit;
381 max_box += padding * dx * unit;
382
383 origin = min_box;
384 maxPoint = max_box;
385
386 sdf.setSpace(min_box, max_box, dx);
387 auto& phi = sdf.distances();
388
389 ni = phi.nx();
390 nj = phi.ny();
391 nk = phi.nz();
392
393 //initialize distances near the mesh
394 DArray3D<uint> counter3d(ni, nj, nk);
395 DArray<uint> counter1d(ni * nj * nk);
396 counter3d.reset();
397
398 DArray3D<GridType> gridType(ni, nj, nk);
399 DArray3D<int> closestTriId(ni, nj, nk);
400
401 cuExecute3D(make_uint3(ni, nj, nk),
402 FSM_AssignInifity,
403 phi,
404 closestTriId,
405 gridType);
406
407 cuExecute(triangles.size(),
408 FSM_FindNeighboringGrids,
409 counter3d,
410 vertices,
411 triangles,
412 origin,
413 dx);
414
415 cuExecute3D(make_uint3(ni, nj, nk),
416 FSM_Array3D_To_Array1d,
417 counter1d,
418 counter3d);
419
420 DArrayList<uint> neighbors;
421 neighbors.resize(counter1d);
422
423 cuExecute(triangles.size(),
424 FSM_StoreTriangleIds,
425 neighbors,
426 counter3d,
427 vertices,
428 triangles,
429 origin,
430 dx);
431
432 cuExecute3D(make_uint3(ni, nj, nk),
433 FSM_InitializeSDFNearMesh,
434 phi,
435 closestTriId,
436 gridType,
437 neighbors,
438 vertices,
439 triangles,
440 edges,
441 tri2edg,
442 edgeNormal,
443 vertexNormal,
444 origin,
445 dx);
446
447 // fill in the rest of the distances with fast sweeping
448 uint passMax = this->varPassNumber()->getValue();
449 for (unsigned int pass = 0; pass < passMax; ++pass) {
450
451 // +x direction
452 cuExecute2D(make_uint2(nj, nk),
453 FSM_SweepX,
454 phi,
455 closestTriId,
456 gridType,
457 vertices,
458 triangles,
459 origin,
460 dx,
461 1);
462
463 // -x direction
464 cuExecute2D(make_uint2(nj, nk),
465 FSM_SweepX,
466 phi,
467 closestTriId,
468 gridType,
469 vertices,
470 triangles,
471 origin,
472 dx,
473 -1);
474
475 // +y direction
476 cuExecute2D(make_uint2(ni, nk),
477 FSM_SweepY,
478 phi,
479 closestTriId,
480 gridType,
481 vertices,
482 triangles,
483 origin,
484 dx,
485 1);
486
487 // -y direction
488 cuExecute2D(make_uint2(ni, nk),
489 FSM_SweepY,
490 phi,
491 closestTriId,
492 gridType,
493 vertices,
494 triangles,
495 origin,
496 dx,
497 -1);
498
499 // +z direction
500 cuExecute2D(make_uint2(ni, nj),
501 FSM_SweepZ,
502 phi,
503 closestTriId,
504 gridType,
505 vertices,
506 triangles,
507 origin,
508 dx,
509 1);
510
511 // -z direction
512 cuExecute2D(make_uint2(ni, nj),
513 FSM_SweepZ,
514 phi,
515 closestTriId,
516 gridType,
517 vertices,
518 triangles,
519 origin,
520 dx,
521 -1);
522 }
523
524 counter3d.clear();
525 counter1d.clear();
526 gridType.clear();
527 neighbors.clear();
528 edgeNormal.clear();
529 vertexNormal.clear();
530 }
531
532 DEFINE_CLASS(FastSweepingMethodGPU);
533}
Generated by doxygen 1.13.2