PeriDyno 1.0.0
Loading...
Searching...
No Matches
PointsBehindMesh.cu
Go to the documentation of this file.
1#include "PointsBehindMesh.h"
2
3#include "GLPointVisualModule.h"
4#include "GLSurfaceVisualModule.h"
5#include "GLWireframeVisualModule.h"
6
7#include <thrust/sort.h>
8namespace dyno
9{
10 template<typename TDataType>
11 PointsBehindMesh<TDataType>::PointsBehindMesh()
12 : Sampler<TDataType>()
13 {
14 this->statePointSet()->setDataPtr(std::make_shared<PointSet<TDataType>>());
15 this->statePointSet()->promoteOuput();
16
17 this->statePlane()->setDataPtr(std::make_shared<TriangleSet<TDataType>>());
18
19 this->statePosition()->allocate();
20 this->statePosition()->promoteOuput();
21
22 this->statePointNormal()->allocate();
23 this->statePointNormal()->promoteOuput();
24
25 this->statePointBelongTriangleIndex()->allocate();
26
27
28 //auto tsRender = std::make_shared<GLSurfaceVisualModule>();
29 //tsRender->setColor(Color(0.1f, 0.12f, 0.25f));
30 //tsRender->varAlpha()->setValue(0.5);
31 //tsRender->setVisible(true);
32 //this->statePlane()->connect(tsRender->inTriangleSet());
33 //this->graphicsPipeline()->pushModule(tsRender);
34
35 auto esRender = std::make_shared<GLWireframeVisualModule>();
36 esRender->varBaseColor()->setValue(Color(0, 0, 0));
37 this->statePlane()->connect(esRender->inEdgeSet());
38 this->graphicsPipeline()->pushModule(esRender);
39
40 auto ptRender = std::make_shared<GLPointVisualModule>();
41 ptRender->setColor(Color(0, 0.5, 1));
42 ptRender->varPointSize()->setValue(0.005f);
43 ptRender->setColorMapMode(GLPointVisualModule::PER_VERTEX_SHADER);
44 this->statePointSet()->connect(ptRender->inPointSet());
45 this->graphicsPipeline()->pushModule(ptRender);
46
47 m_NeighborPointQuery = std::make_shared<NeighborPointQuery<TDataType>>();
48 this->statePosition()->connect(m_NeighborPointQuery->inPosition());
49 };
50
51
52 template<typename Coord>
53 __device__ Real WhetherPointInsideTriangle(Coord t_a, Coord t_b, Coord t_c, Coord point_o) {
54
55 /*@Brief : Caculate triangle arear.
56 */
57
58 Coord ab = t_b - t_a;
59 Coord ac = t_c - t_a;
60 Real abc = (ab.cross(ac)).norm() * 0.5;
61
62 Coord ao = point_o - t_a;
63 Coord bo = point_o - t_b;
64 Coord bc = t_c - t_b;
65
66 Real abo = (ab.cross(ao)).norm() * 0.5;
67 Real aco = (ac.cross(ao)).norm() * 0.5;
68 Real bco = (bc.cross(bo)).norm() * 0.5;
69
70 return ((abo + aco + bco) - abc);
71 }
72
73
74
75 template<typename Coord, typename Triangle>
76 __global__ void CalculateSquareFromTriangle(
77 DArray<Coord> square_a,
78 DArray<Coord> square_b,
79 DArray<Coord> square_c,
80 DArray<Coord> square_d,
81 DArray<Triangle> triangles,
82 DArray<Coord> vertices
83 )
84 {
85 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
86 if (tId >= triangles.size()) return;
87
88 Triangle t = triangles[tId];
89 const Coord& v0 = vertices[t[0]];
90 const Coord& v1 = vertices[t[1]];
91 const Coord& v2 = vertices[t[2]];
92
93 Real lv01 = (v0 - v1).norm();
94 Real lv12 = (v1 - v2).norm();
95 Real lv02 = (v0 - v2).norm();
96
97 Coord C(0.0f);
98 if ((lv01 >= lv12) && (lv01 >= lv02))
99 {
100 square_a[tId] = v0;
101 square_b[tId] = v1;
102 C = v2;
103 }
104 else if ((lv12 >= lv01) && (lv12 >= lv02))
105 {
106 square_a[tId] = v1;
107 square_b[tId] = v2;
108 C = v0;
109 }
110 else if ((lv02 >= lv01) && (lv02 >= lv12))
111 {
112 square_a[tId] = v0;
113 square_b[tId] = v2;
114 C = v1;
115 }
116
117 Coord AB = square_b[tId] - square_a[tId];
118 Coord n_AB = AB / AB.norm();
119 Coord AC = C - square_a[tId];
120 Coord proj_O = square_a[tId] + (AC.dot(n_AB)) * (n_AB);
121 Coord trans_vec = C - proj_O;
122
123 square_c[tId] = square_a[tId] + trans_vec;
124 square_d[tId] = square_b[tId] + trans_vec;
125 }
126
127
128 template<typename Coord, typename Triangle>
129 __global__ void CalculateNewTriangle(
130 DArray<Triangle> NewTriangles,
131 DArray<Coord> square_a,
132 DArray<Coord> square_b,
133 DArray<Coord> square_c,
134 DArray<Coord> square_d,
135 DArray<Coord> vertices
136 )
137 {
138 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
139 if (tId >= square_a.size()) return;
140
141 int v_id = tId * 4;
142
143 vertices[v_id] = square_a[tId];
144 vertices[v_id + 1] = square_b[tId];
145 vertices[v_id + 2] = square_c[tId];
146 vertices[v_id + 3] = square_d[tId];
147
148 int t_id = tId * 2;
149
150 Triangle& t0 = NewTriangles[t_id];
151 t0[0] = v_id;
152 t0[1] = v_id + 1;
153 t0[2] = v_id + 2;
154
155 Triangle& t1 = NewTriangles[t_id + 1];
156
157 t1[0] = v_id + 1;
158 t1[1] = v_id + 2;
159 t1[2] = v_id + 3;
160
161 }
162
163
164 template<typename Coord, typename Triangle, typename Real>
165 __global__ void CalculateTriangleBasicVector(
166 DArray<Triangle> triangles,
167 DArray<Coord> vertices,
168 DArray<Coord> basic_x,
169 DArray<Coord> basic_y,
170 DArray<Coord> basic_z,
171 DArray<Coord> square_a,
172 DArray<Coord> square_b,
173 DArray<Coord> square_c,
174 DArray<Coord> square_d,
175 Real dx
176 )
177 {
178 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
179 if (tId >= triangles.size()) return;
180
181 Coord& n_x = basic_x[tId];
182 Coord& n_y = basic_y[tId];
183 Coord& n_z = basic_z[tId];
184
185 n_x = (square_b[tId] - square_a[tId]) / (square_b[tId] - square_a[tId]).norm();
186 n_y = (square_c[tId] - square_a[tId]) / (square_c[tId] - square_a[tId]).norm();
187
188 n_z = n_x.cross(n_y);
189 n_z = n_z / n_z.norm();
190 }
191
192
193 template<typename Coord, typename Real, typename Triangle>
194 __global__ void CalculatePointSizeInSquare(
195 DArray<int> PointSize,
196 DArray<Triangle> triangles,
197 DArray<Coord> vertices,
198 DArray<Coord> basic_x,
199 DArray<Coord> basic_y,
200 DArray<Coord> basic_z,
201 DArray<Coord> square_a,
202 DArray<Coord> square_b,
203 DArray<Coord> square_c,
204 DArray<Coord> square_d,
205 Real dx
206 )
207 {
208 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
209 if (tId >= PointSize.size()) return;
210
211 int& size = PointSize[tId];
212 PointSize[tId] = 0;
213
214 int num_x = (int)((square_a[tId] - square_c[tId]).norm() / dx);
215 int num_y = (int)((square_a[tId] - square_b[tId]).norm() / dx);
216
217 Triangle t = triangles[tId];
218 const Coord& v0 = vertices[t[0]];
219 const Coord& v1 = vertices[t[1]];
220 const Coord& v2 = vertices[t[2]];
221
222 for (int j = 0; j <= num_y + 1; j++)
223 {
224 for (int i = 0; i <= num_x + 1; i++)
225 {
226 int count = i + j * num_x;
227
228 //if (count < size) {
229 Coord candi = square_a[tId] + (Real)(j)*dx * basic_x[tId] + (Real)(i)*dx * basic_y[tId];
230 Real value = WhetherPointInsideTriangle(v0, v1, v2, candi);
231 if (value < 100000000 * dx * dx * EPSILON) {
232 size++;
233 }
234
235 }
236 }
237 }
238
239
240 template<typename Matrix, typename Coord, typename Real, typename Triangle>
241 __global__ void CalculatePointsOnPlane(
242 DArray<Matrix> Rotation,
243 DArray<Coord> Normal,
244 DArray<Triangle> triangles,
245 Real value
246 )
247 {
248 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
249 if (tId >= Normal.size()) return;
250
251 }
252
253
254 template<typename Coord, typename Real, typename Triangle>
255 __global__ void CalculatePointInSquare(
256 DArray<Coord> Points,
257 DArray<int> PointSize,
258 DArray<int> tIndex,
259 DArray<Coord> PointNormal,
260 DArray<int> PointOfTriangle,
261 DArray<Triangle> triangles,
262 DArray<Coord> vertices,
263 DArray<Coord> basic_x,
264 DArray<Coord> basic_y,
265 DArray<Coord> PlaneNormal,
266 DArray<Coord> square_a,
267 DArray<Coord> square_b,
268 DArray<Coord> square_c,
269 DArray<Coord> square_d,
270 Real dx
271 )
272 {
273
274 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
275 if (tId >= PointSize.size()) return;
276
277
278 int num_x = (int)((square_a[tId] - square_c[tId]).norm() / dx);
279 int num_y = (int)((square_a[tId] - square_b[tId]).norm() / dx);
280 int& begin = tIndex[tId];
281 int& size = PointSize[tId];
282
283
284 Triangle t = triangles[tId];
285 const Coord& v0 = vertices[t[0]];
286 const Coord& v1 = vertices[t[1]];
287 const Coord& v2 = vertices[t[2]];
288
289 int counter = 0;
290
291 for (int j = 0; j <= num_y + 1; j++)
292 {
293 for (int i = 0; i <= num_x + 1; i++)
294 {
295 //int count = i + j * num_x;
296 if (counter < size) {
297
298 Coord candi = square_a[tId] + (Real)(j)*dx * basic_x[tId] + (Real)(i)*dx * basic_y[tId];
299 Real value = WhetherPointInsideTriangle(v0, v1, v2, candi);
300
301 if (value < 100000000 * dx * dx * EPSILON)
302 {
303 Points[begin + counter] = candi;
304 PointOfTriangle[begin + counter] = tId;
305 PointNormal[begin + counter] = -1.0 * PlaneNormal[tId];
306 counter++;
307 }
308
309 }
310 }
311 }
312 }
313
314 template<typename Coord, typename Real>
315 __global__ void CalculateGrowingPointSize(
316 DArray<Coord> GrowingPoints,
317 DArray<Coord> SeedPoints,
318 DArray<Coord> Normals,
319 DArray<Coord> outPointNormals,
320 DArray<int> SeedOfTriangleId,
321 DArray<int> PointOfTriangleId,
322 Real thickness,
323 Real dx,
324 bool direction,
325 int pointsLayer
326 )
327 {
328 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
329 if (tId >= SeedPoints.size()) return;
330
331 if (direction == true)
332 {
333 SeedPoints[tId] -= Normals[tId] * dx * 0.5;
334 for (int i = 0; i < pointsLayer; i++)
335 {
336 GrowingPoints[pointsLayer * tId + i] = SeedPoints[tId] - (Real)(i)*dx * Normals[tId];
337 outPointNormals[pointsLayer * tId + i] = Normals[tId];
338 PointOfTriangleId[pointsLayer * tId + i] = SeedOfTriangleId[tId];
339 }
340 }
341 else
342 {
343 SeedPoints[tId] += Normals[tId] * dx * 0.5;
344 for (int i = 0; i < pointsLayer; i++)
345 {
346 GrowingPoints[pointsLayer * tId + i] = SeedPoints[tId] + (Real)(i)*dx * Normals[tId];
347 outPointNormals[pointsLayer * tId + i] = -Normals[tId];
348 PointOfTriangleId[pointsLayer * tId + i] = SeedOfTriangleId[tId];
349 }
350 }
351 }
352
353
354
355 template<typename Coord, typename Real>
356 __global__ void CalculateRemovingPointFlag(
357 DArray<bool> flags,
358 DArray<Coord> positions,
359 DArrayList<int> neighbors,
360 Real threshold_r
361 )
362 {
363 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
364 if (pId >= flags.size()) return;
365
366 flags[pId] = false;
367
368 List<int>& list_i = neighbors[pId];
369 int nbSize = list_i.size();
370 for (int ne = 0; ne < nbSize; ne++)
371 {
372 int j = list_i[ne];
373 if (pId == j) continue;
374 Real r = (positions[pId] - positions[j]).norm();
375 if ((threshold_r > r) && (pId < j)) //
376 {
377 flags[pId] = true;
378 }
379 }
380 }
381
382
383 __global__ void CalculateGhostCounter(
384 DArray<int> counter,
385 DArray<bool> flags
386 )
387 {
388 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
389 if (tId >= flags.size()) return;
390
391 counter[tId] = flags[tId] ? 0 : 1;
392 }
393
394
395 template<typename Coord>
396 __global__ void CalculateLastGhostPoints(
397 DArray<Coord> lastPoints,
398 DArray<Coord> originalPoints,
399 DArray<Coord> lastNormals,
400 DArray<Coord> originalNormals,
401 DArray<int> radix
402 )
403 {
404 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
405 if (pId >= originalPoints.size()) return;
406
407 if (pId == originalPoints.size() - 1 || radix[pId] != radix[pId + 1])
408 {
409 lastPoints[radix[pId]] = originalPoints[pId];
410 lastNormals[radix[pId]] = originalNormals[pId];
411 }
412 }
413
414 template< typename Coord>
415 __global__ void PtsBehindMesh_UpdateTriangleNormal(
416 DArray<TopologyModule::Triangle> d_triangles,
417 DArray<Coord> d_points,
418 DArray<Coord> normal,
419 float length,
420 bool normalization)
421 {
422 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
423 if (pId >= d_triangles.size()) return;
424
425 int a = d_triangles[pId][0];
426 int b = d_triangles[pId][1];
427 int c = d_triangles[pId][2];
428
429 float x = (d_points[a][0] + d_points[b][0] + d_points[c][0]) / 3;
430 float y = (d_points[a][1] + d_points[b][1] + d_points[c][1]) / 3;
431 float z = (d_points[a][2] + d_points[b][2] + d_points[c][2]) / 3;
432
433 Coord ca = d_points[b] - d_points[a];
434 Coord cb = d_points[b] - d_points[c];
435 Coord dirNormal;
436
437 if (normalization)
438 dirNormal = ca.cross(cb).normalize() * -1 * length;
439 else
440 dirNormal = ca.cross(cb) * -1 * length;
441
442 normal[pId] = dirNormal.normalize();
443
444 }
445
446 template<typename TDataType>
447 void PointsBehindMesh<TDataType>::resetStates()
448 {
449 //Normal<TDataType>::resetStates();
450
451 int triangle_num = this->inTriangleSet()->getData().getTriangles().size();
452 //auto& m_triangle_normal = this->stateNormal()->getData();
453
454 if (triangle_num == 0) return;
455
456 this->outPointGrowthDirection()->setValue(this->varGeneratingDirection()->getValue());
457
458 Real dx = this->varSamplingDistance()->getValue();
459
460 if (mTriangleNormal.size() != triangle_num)
461 {
462 mTriangleNormal.resize(triangle_num);
463 }
464
465 if (msquare_1.size() != triangle_num)
466 {
467 msquare_1.resize(triangle_num);
468 msquare_2.resize(triangle_num);
469 msquare_3.resize(triangle_num);
470 msquare_4.resize(triangle_num);
471 }
472
473 if (mBasicVector_x.size() != triangle_num)
474 {
475 mBasicVector_x.resize(triangle_num);
476 mBasicVector_y.resize(triangle_num);
477 mBasicVector_z.resize(triangle_num);
478 mThinPointSize.resize(triangle_num);
479 }
480
481 if (mTriangleTempt.size() != 2 * triangle_num)
482 mTriangleTempt.resize(2 * triangle_num);
483 if (mVerticesTempt.size() != 4 * triangle_num)
484 mVerticesTempt.resize(4 * triangle_num);
485
486 cuExecute(triangle_num,
487 PtsBehindMesh_UpdateTriangleNormal,
488 this->inTriangleSet()->getData().getTriangles(), //DArray<TopologyModule::Triangle> d_triangles,
489 this->inTriangleSet()->getData().getPoints(),
490 //DArray<TopologyModule::Edge> edges,
491 //DArray<Coord> normal_points,
492 mTriangleNormal,
493 //DArray<Coord> triangleCenter,
494 0.2f,
495 true
496 );
497
498
499
500 cuExecute(triangle_num,
501 CalculateSquareFromTriangle,
502 msquare_1,
503 msquare_2,
504 msquare_3,
505 msquare_4,
506 this->inTriangleSet()->getData().getTriangles(),
507 this->inTriangleSet()->getData().getPoints()
508 );
509
510 cuExecute(triangle_num,
511 CalculateNewTriangle,
512 mTriangleTempt,
513 msquare_1,
514 msquare_2,
515 msquare_3,
516 msquare_4,
517 mVerticesTempt
518 );
519
520 auto& temp_triangleSet = this->statePlane()->getData();
521 temp_triangleSet.setTriangles(mTriangleTempt);
522 temp_triangleSet.setPoints(mVerticesTempt);
523
524 cuExecute(triangle_num,
525 CalculateTriangleBasicVector,
526 this->inTriangleSet()->getData().getTriangles(),
527 this->inTriangleSet()->getData().getPoints(),
528 mBasicVector_x,
529 mBasicVector_y,
530 mBasicVector_z,
531 msquare_1,
532 msquare_2,
533 msquare_3,
534 msquare_4,
535 this->varSamplingDistance()->getValue()
536 );
537
538 cuExecute(triangle_num,
539 CalculatePointSizeInSquare,
540 mThinPointSize,
541 this->inTriangleSet()->getData().getTriangles(),
542 this->inTriangleSet()->getData().getPoints(),
543 mBasicVector_x,
544 mBasicVector_y,
545 mBasicVector_z,
546 msquare_1,
547 msquare_2,
548 msquare_3,
549 msquare_4,
550 this->varSamplingDistance()->getValue()
551 );
552
553 Reduction<int> reduce;
554 int total_ghostpoint_number = reduce.accumulate(mThinPointSize.begin(), mThinPointSize.size());
555
556 if (total_ghostpoint_number <= 0)
557 {
558 return;
559 }
560
561 DArray<int> tIndex;
562 tIndex.assign(mThinPointSize);
563
564
565 DArray<Coord> mThinPointNormal;
566 mThinPointNormal.resize(total_ghostpoint_number);
567
568 mThinPoints.resize(total_ghostpoint_number);
569 mSeedOfTriangleId.resize(total_ghostpoint_number);
570
571 thrust::exclusive_scan(thrust::device, tIndex.begin(), tIndex.begin() + tIndex.size(), tIndex.begin());
572
573 cuExecute(triangle_num,
574 CalculatePointInSquare,
575 mThinPoints,
576 mThinPointSize,
577 tIndex,
578 mThinPointNormal,
579 mSeedOfTriangleId,
580 this->inTriangleSet()->getData().getTriangles(),
581 this->inTriangleSet()->getData().getPoints(),
582 mBasicVector_x,
583 mBasicVector_y,
584 mTriangleNormal, //this->stateNormal()->getData(),
585 msquare_1,
586 msquare_2,
587 msquare_3,
588 msquare_4,
589 this->varSamplingDistance()->getValue()
590 );
591
592 int pointLayerCount = (int)(this->varThickness()->getValue() / this->varSamplingDistance()->getValue());
593 pointLayerCount = pointLayerCount < 1 ? 1 : pointLayerCount;
594
595 mThickPoints.resize(mThinPoints.size() * pointLayerCount);
596
597 mPointOfTriangleId.resize(mThinPoints.size() * pointLayerCount);
598
599
600 DArray<Coord> mThickPointNormals;
601 mThickPointNormals.resize(mThinPoints.size() * pointLayerCount);
602
603 cuExecute(mThinPoints.size(),
604 CalculateGrowingPointSize,
605 mThickPoints,
606 mThinPoints,
607 mThinPointNormal,
608 mThickPointNormals,
609 mSeedOfTriangleId,
610 mPointOfTriangleId,
611 this->varThickness()->getValue(),
612 dx,
613 this->varGeneratingDirection()->getValue(),
614 pointLayerCount
615 );
616
617 if (mRemovingFlag.size() != mThickPoints.size())
618 {
619 mRemovingFlag.resize(mThickPoints.size());
620 }
621
622 this->statePosition()->getData().assign(mThickPoints);
623 this->statePosition()->connect(m_NeighborPointQuery->inPosition());
624 m_NeighborPointQuery->inRadius()->setValue(this->varSamplingDistance()->getValue() * 1.0);
625 m_NeighborPointQuery->update();
626
627 cuExecute(mThickPoints.size(),
628 CalculateRemovingPointFlag,
629 mRemovingFlag,
630 mThickPoints,
631 m_NeighborPointQuery->outNeighborIds()->getData(),
632 0.75f * this->varSamplingDistance()->getValue()
633 );
634
635 DArray<int> ghost_counter(mRemovingFlag.size());
636
637 cuExecute(mThickPoints.size(),
638 CalculateGhostCounter,
639 ghost_counter,
640 mRemovingFlag
641 );
642
643 int last_ghostpoint_number = reduce.accumulate(ghost_counter.begin(), ghost_counter.size());
644 scan.exclusive(ghost_counter.begin(), ghost_counter.size());
645
646 DArray<Coord> lastPoints(last_ghostpoint_number);
647 auto& lastPointNormal = this->statePointNormal()->getData();
648 lastPointNormal.resize(last_ghostpoint_number);
649
650 this->statePointNormal()->resize(last_ghostpoint_number);
651
652 if (last_ghostpoint_number > 0)
653 {
654 cuExecute(mThickPoints.size(),
655 CalculateLastGhostPoints,
656 lastPoints,
657 mThickPoints,
658 this->statePointNormal()->getData(),
659 mThickPointNormals,
660 ghost_counter
661 );
662 }
663
664 auto& ghostPointSet = this->statePointSet()->getData();
665 this->statePosition()->getData().clear();
666 this->statePosition()->getData().assign(lastPoints);
667 this->statePointSet()->getData().clear();
668 this->statePointSet()->getData().setPoints(lastPoints);
669
670 this->outSamplingDistance()->setValue(this->varSamplingDistance()->getValue());
671 this->statePointBelongTriangleIndex()->assign(mPointOfTriangleId);
672
673
674
675
676 tIndex.clear();
677 lastPoints.clear();
678 ghost_counter.clear();
679 mThinPointNormal.clear();
680 mThickPointNormals.clear();
681
682 }
683
684
685 DEFINE_CLASS(PointsBehindMesh);
686}
Generated by doxygen 1.13.2