PeriDyno 1.0.0
Loading...
Searching...
No Matches
TriangleSet.cu
Go to the documentation of this file.
1#include "TriangleSet.h"
2#include <fstream>
3#include <iostream>
4#include <sstream>
5
6#include <thrust/sort.h>
7#define TINYOBJLOADER_IMPLEMENTATION
8#include "tinyobjloader/tiny_obj_loader.h"
9
10namespace dyno
11{
12 template<typename TDataType>
13 TriangleSet<TDataType>::TriangleSet()
14 : EdgeSet<TDataType>()
15 {
16 }
17
18 template<typename TDataType>
19 TriangleSet<TDataType>::~TriangleSet()
20 {
21 mTriangleIndex.clear();
22 mVertexNormal.clear();
23 mVer2Tri.clear();
24 mEdg2Tri.clear();
25 mTri2Edg.clear();
26 }
27
28 template<typename Triangle>
29 __global__ void TS_CountTriangles(
30 DArray<uint> counter,
31 DArray<Triangle> triangles)
32 {
33 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
34 if (tId >= triangles.size()) return;
35
36 Triangle t = triangles[tId];
37
38 atomicAdd(&counter[t[0]], 1);
39 atomicAdd(&counter[t[1]], 1);
40 atomicAdd(&counter[t[2]], 1);
41 }
42
43 template<typename Triangle>
44 __global__ void TS_SetupTriIds(
45 DArrayList<int> triIds,
46 DArray<Triangle> triangles)
47 {
48 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
49 if (tId >= triangles.size()) return;
50
51 Triangle t = triangles[tId];
52
53 triIds[t[0]].atomicInsert(tId);
54 triIds[t[1]].atomicInsert(tId);
55 triIds[t[2]].atomicInsert(tId);
56 }
57
58 template<typename TDataType>
59 DArrayList<int>& TriangleSet<TDataType>::getVertex2Triangles()
60 {
61 DArray<uint> counter(this->mCoords.size());
62 counter.reset();
63
64 cuExecute(mTriangleIndex.size(),
65 TS_CountTriangles,
66 counter,
67 mTriangleIndex);
68
69 mVer2Tri.resize(counter);
70
71 counter.reset();
72 cuExecute(mTriangleIndex.size(),
73 TS_SetupTriIds,
74 mVer2Tri,
75 mTriangleIndex);
76
77 counter.clear();
78
79 return mVer2Tri;
80 }
81
82 template<typename Edg2Tri>
83 __global__ void TS_setupIds(
84 DArray<int> edgIds,
85 DArray<int> triIds,
86 DArray<Edg2Tri> edg2Tri)
87 {
88 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
89 if (tId >= edg2Tri.size()) return;
90
91 Edg2Tri e = edg2Tri[tId];
92
93 triIds[2 * tId] = e[0];
94 triIds[2 * tId + 1] = e[1];
95
96 edgIds[2 * tId] = tId;
97 edgIds[2 * tId + 1] = tId;
98 }
99
100 template<typename Tri2Edg, typename Edge, typename Triangle>
101 __global__ void TS_SetupTri2Edg(
102 DArray<Tri2Edg> tri2Edg,
103 DArray<int> triIds,
104 DArray<int> edgIds,
105 DArray<Edge> edges,
106 DArray<Triangle> triangles)
107 {
108 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
109 if (tId >= triIds.size()) return;
110
111 if (tId == 0 || triIds[tId] != triIds[tId - 1])
112 {
113 Tri2Edg t2E(EMPTY, EMPTY, EMPTY);
114
115 EKey te0(triangles[triIds[tId]][0], triangles[triIds[tId]][1]);
116 EKey te1(triangles[triIds[tId]][1], triangles[triIds[tId]][2]);
117 EKey te2(triangles[triIds[tId]][2], triangles[triIds[tId]][0]);
118
119 EKey e0(edges[edgIds[tId]][0], edges[edgIds[tId]][1]);
120 EKey e1(edges[edgIds[tId + 1]][0], edges[edgIds[tId + 1]][1]);
121 EKey e2(edges[edgIds[tId + 2]][0], edges[edgIds[tId + 2]][1]);
122
123 if (te0 == e0)
124 t2E[0] = edgIds[tId];
125 else if (te0 == e1)
126 t2E[0] = edgIds[tId + 1];
127 else if (te0 == e2)
128 t2E[0] = edgIds[tId + 2];
129
130 if (te1 == e0)
131 t2E[1] = edgIds[tId];
132 else if (te1 == e1)
133 t2E[1] = edgIds[tId + 1];
134 else if (te1 == e2)
135 t2E[1] = edgIds[tId + 2];
136
137 if (te2 == e0)
138 t2E[2] = edgIds[tId];
139 else if (te2 == e1)
140 t2E[2] = edgIds[tId + 1];
141 else if (te2 == e2)
142 t2E[2] = edgIds[tId + 2];
143
144 int shift = tId / 3;
145 tri2Edg[shift] = t2E;
146
147 //printf("tri2Edg: %d, %d %d %d, %d %d %d \n", shift, triangles[triIds[tId]][0], triangles[triIds[tId]][1], triangles[triIds[tId]][2],tri2Edg[shift][0], tri2Edg[shift][1], tri2Edg[shift][2]);
148 }
149 }
150
151 template<typename TDataType>
152 void TriangleSet<TDataType>::updateTriangle2Edge()
153 {
154 if (mEdg2Tri.size() == 0)
155 this->updateEdges();
156
157 uint edgSize = mEdg2Tri.size();
158
159 DArray<int> triIds, edgIds;
160 triIds.resize(2 * edgSize);
161 edgIds.resize(2 * edgSize);
162
163 cuExecute(edgSize,
164 TS_setupIds,
165 edgIds,
166 triIds,
167 mEdg2Tri);
168
169 thrust::sort_by_key(thrust::device, triIds.begin(), triIds.begin() + triIds.size(), edgIds.begin());
170
171 auto& pEdges = this->getEdges();
172
173 mTri2Edg.resize(mTriangleIndex.size());
174 cuExecute(triIds.size(),
175 TS_SetupTri2Edg,
176 mTri2Edg,
177 triIds,
178 edgIds,
179 pEdges,
180 mTriangleIndex);
181
182 triIds.clear();
183 edgIds.clear();
184 }
185
186 template<typename EKey, typename Triangle>
187 __global__ void TS_SetupKeys(
188 DArray<EKey> keys,
189 DArray<int> ids,
190 DArray<Triangle> triangles)
191 {
192 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
193 if (tId >= triangles.size()) return;
194
195 Triangle tri = triangles[tId];
196 keys[3 * tId] = EKey(tri[0], tri[1]);
197 keys[3 * tId + 1] = EKey(tri[1], tri[2]);
198 keys[3 * tId + 2] = EKey(tri[2], tri[0]);
199
200 ids[3 * tId] = tId;
201 ids[3 * tId + 1] = tId;
202 ids[3 * tId + 2] = tId;
203 }
204
205 template<typename EKey>
206 __global__ void TS_CountEdgeNumber(
207 DArray<int> counter,
208 DArray<EKey> keys)
209 {
210 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
211 if (tId >= keys.size()) return;
212
213 if (tId == 0 || keys[tId] != keys[tId - 1])
214 counter[tId] = 1;
215 else
216 counter[tId] = 0;
217 }
218
219 template<typename Edge, typename Edg2Tri, typename EKey>
220 __global__ void TS_SetupEdges(
221 DArray<Edge> edges,
222 DArray<Edg2Tri> edg2Tri,
223 DArray<EKey> keys,
224 DArray<int> counter,
225 DArray<int> triIds)
226 {
227 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
228 if (tId >= keys.size()) return;
229
230 int shift = counter[tId];
231 if (tId == 0 || keys[tId] != keys[tId - 1])
232 {
233 EKey key = keys[tId];
234 edges[shift] = Edge(key[0], key[1]);
235
236 Edg2Tri e2T(EMPTY, EMPTY);
237 e2T[0] = triIds[tId];
238
239 if (tId + 1 < keys.size() && keys[tId + 1] == key)
240 e2T[1] = triIds[tId + 1];
241
242 edg2Tri[shift] = e2T;
243 }
244 }
245
246 template<typename TDataType>
247 void TriangleSet<TDataType>::updateEdges()
248 {
249 uint triSize = mTriangleIndex.size();
250
251 DArray<EKey> keys;
252 DArray<int> triIds;
253
254 keys.resize(3 * triSize);
255 triIds.resize(3 * triSize);
256
257 cuExecute(triSize,
258 TS_SetupKeys,
259 keys,
260 triIds,
261 mTriangleIndex);
262
263 thrust::sort_by_key(thrust::device, keys.begin(), keys.begin() + keys.size(), triIds.begin());
264
265 DArray<int> counter;
266 counter.resize(3 * triSize);
267
268 cuExecute(keys.size(),
269 TS_CountEdgeNumber,
270 counter,
271 keys);
272
273 int edgeNum = thrust::reduce(thrust::device, counter.begin(), counter.begin() + counter.size());
274 thrust::exclusive_scan(thrust::device, counter.begin(), counter.begin() + counter.size(), counter.begin());
275
276 mEdg2Tri.resize(edgeNum);
277
278 auto& pEdges = this->getEdges();
279 pEdges.resize(edgeNum);
280 cuExecute(keys.size(),
281 TS_SetupEdges,
282 pEdges,
283 mEdg2Tri,
284 keys,
285 counter,
286 triIds);
287
288 counter.clear();
289 triIds.clear();
290 keys.clear();
291 }
292
293 template<typename TDataType>
294 void TriangleSet<TDataType>::setTriangles(std::vector<Triangle>& triangles)
295 {
296 mTriangleIndex.resize(triangles.size());
297 mTriangleIndex.assign(triangles);
298 }
299
300 template<typename TDataType>
301 void TriangleSet<TDataType>::setTriangles(DArray<Triangle>& triangles)
302 {
303 mTriangleIndex.resize(triangles.size());
304 mTriangleIndex.assign(triangles);
305 }
306
307 template<typename TDataType>
308 bool TriangleSet<TDataType>::loadObjFile(std::string filename)
309 {
310 std::vector<Coord> vertList;
311 std::vector<Triangle> faceList;
312
313 tinyobj::attrib_t myattrib;
314 std::vector <tinyobj::shape_t> myshape;
315 std::vector <tinyobj::material_t> mymat;
316 std::string mywarn;
317 std::string myerr;
318
319 char* fname = (char*)filename.c_str();
320
321 bool succeed = tinyobj::LoadObj(&myattrib, &myshape, &mymat, &mywarn, &myerr, fname, nullptr ,true, true);
322 if (!succeed)
323 return false;
324
325 for (int i = 0; i < myattrib.GetVertices().size() / 3; i++)
326 {
327 vertList.push_back(Coord(myattrib.GetVertices()[3 * i], myattrib.GetVertices()[3 * i + 1], myattrib.GetVertices()[3 * i + 2]));
328 }
329
330 for (int i = 0;i < myshape.size();i++)
331 {
332 for (int s = 0;s < myshape[i].mesh.indices.size()/3; s++)
333 {
334 faceList.push_back(Triangle(myshape[i].mesh.indices[3 * s].vertex_index, myshape[i].mesh.indices[3 * s + 1].vertex_index, myshape[i].mesh.indices[3 * s + 2].vertex_index));
335 }
336 }
337 this->setPoints(vertList);
338 this->setTriangles(faceList);
339 this->update();
340
341 vertList.clear();
342 faceList.clear();
343 myshape.clear();
344 mymat.clear();
345
346 return true;
347 }
348
349 template<typename TDataType>
350 void TriangleSet<TDataType>::copyFrom(TriangleSet<TDataType>& triangleSet)
351 {
352 mVer2Tri.assign(triangleSet.mVer2Tri);
353
354 mTriangleIndex.resize(triangleSet.mTriangleIndex.size());
355 mTriangleIndex.assign(triangleSet.mTriangleIndex);
356
357 mEdg2Tri.resize(triangleSet.mEdg2Tri.size());
358 mEdg2Tri.assign(triangleSet.mEdg2Tri);
359
360 EdgeSet<TDataType>::copyFrom(triangleSet);
361 }
362
363 template<typename Triangle>
364 __global__ void TS_UpdateIndex(
365 DArray<Triangle> indices,
366 uint vertexOffset,
367 uint indexOffset)
368 {
369 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
370 if (tId >= indices.size() - indexOffset) return;
371
372 Triangle t = indices[indexOffset + tId];
373 t[0] += vertexOffset;
374 t[1] += vertexOffset;
375 t[2] += vertexOffset;
376
377 indices[indexOffset + tId] = t;
378 }
379
380 template<typename TDataType>
381 std::shared_ptr<TriangleSet<TDataType>> TriangleSet<TDataType>::merge(TriangleSet<TDataType>& ts)
382 {
383 auto ret = std::make_shared<TriangleSet<TDataType>>();
384
385 auto& vertices = ret->getPoints();
386 auto& indices = ret->getTriangles();
387
388 uint vNum0 = PointSet<TDataType>::mCoords.size();
389 uint vNum1 = ts.getPoints().size();
390
391 uint tNum0 = mTriangleIndex.size();
392 uint tNum1 = ts.getTriangles().size();
393
394 vertices.resize(vNum0 + vNum1);
395 indices.resize(tNum0 + tNum1);
396
397 vertices.assign(PointSet<TDataType>::mCoords, vNum0, 0, 0);
398 vertices.assign(ts.getPoints(), vNum1, vNum0, 0);
399
400 indices.assign(mTriangleIndex, tNum0, 0, 0);
401 indices.assign(ts.getTriangles(), tNum1, tNum0, 0);
402
403 cuExecute(tNum1,
404 TS_UpdateIndex,
405 indices,
406 vNum0,
407 tNum0);
408
409 return ret;
410 }
411
412 template<typename TDataType>
413 bool TriangleSet<TDataType>::isEmpty()
414 {
415 return mTriangleIndex.size() == 0 && EdgeSet<TDataType>::isEmpty();
416 }
417
418 template<typename TDataType>
419 void TriangleSet<TDataType>::clear()
420 {
421 mTriangleIndex.clear();
422 mVer2Tri.clear();
423 mEdg2Tri.clear();
424 mTri2Edg.clear();
425
426 mVertexNormal.clear();
427
428 EdgeSet<TDataType>::clear();
429 }
430
431 template<typename Coord, typename Triangle>
432 __global__ void TS_SetupVertexNormals(
433 DArray<Coord> normals,
434 DArray<Coord> vertices,
435 DArray<Triangle> triangles,
436 DArrayList<int> triIds)
437 {
438 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
439 if (tId >= normals.size()) return;
440
441 List<int>& list_i = triIds[tId];
442 int triSize = list_i.size();
443
444 Coord N = Coord(0);
445 for (int ne = 0; ne < triSize; ne++)
446 {
447 int j = list_i[ne];
448 Triangle t = triangles[j];
449
450 Coord v0 = vertices[t[0]];
451 Coord v1 = vertices[t[1]];
452 Coord v2 = vertices[t[2]];
453
454 N += (v1 - v0).cross(v2 - v0);
455 }
456
457 N.normalize();
458
459 normals[tId] = N;
460 }
461
462 template<typename TDataType>
463 void TriangleSet<TDataType>::setNormals(DArray<Coord>& normals)
464 {
465 mVertexNormal.assign(normals);
466 }
467
468 template<typename TDataType>
469 void TriangleSet<TDataType>::updateVertexNormal()
470 {
471 uint vertSize = this->mCoords.size();
472
473 if (vertSize <= 0)
474 return;
475
476 if (mVertexNormal.size() != vertSize) {
477 mVertexNormal.resize(vertSize);
478 }
479
480 auto& vert2Tri = getVertex2Triangles();
481 cuExecute(vertSize,
482 TS_SetupVertexNormals,
483 mVertexNormal,
484 this->mCoords,
485 mTriangleIndex,
486 vert2Tri);
487 }
488
489 template<typename Coord, typename Triangle>
490 __global__ void TS_SetupAngleWeightedVertexNormals(
491 DArray<Coord> normals,
492 DArray<Coord> vertices,
493 DArray<Triangle> triangles,
494 DArrayList<int> triIds)
495 {
496 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
497 if (tId >= normals.size()) return;
498
499 List<int>& list_i = triIds[tId];
500 int triSize = list_i.size();
501
502 Coord N = Coord(0);
503 for (int ne = 0; ne < triSize; ne++)
504 {
505 int j = list_i[ne];
506 Triangle t = triangles[j];
507
508 Coord v0 = vertices[t[0]];
509 Coord v1 = vertices[t[1]];
510 Coord v2 = vertices[t[2]];
511
512 Real e0 = (v1 - v2).norm();
513 Real e1 = (v2 - v0).norm();
514 Real e2 = (v1 - v0).norm();
515
516 Real cosangle = 0;
517 if (t[0] == tId)
518 cosangle = (e1 * e1 + e2 * e2 - e0 * e0) / (2.0 * e1 * e2);
519 else if (t[1] == tId)
520 cosangle = (e0 * e0 + e2 * e2 - e1 * e1) / (2.0 * e0 * e2);
521 else if (t[2] == tId)
522 cosangle = (e1 * e1 + e0 * e0 - e2 * e2) / (2.0 * e1 * e0);
523
524 Real angle = acos(cosangle);
525 Coord norm = (v1 - v0).cross(v2 - v0);
526 norm.normalize();
527 N += angle * norm;
528
529 //printf("vertex normal: %d, %f %f %f, %d %f, %f %f %f, %f %f %f \n", tId, vertices[tId][0], vertices[tId][1], vertices[tId][2],
530 // j, angle, norm[0], norm[1], norm[2], N[0], N[1], N[2]);
531 }
532
533 N.normalize();
534
535 normals[tId] = N;
536 }
537
538 template<typename TDataType>
539 void TriangleSet<TDataType>::updateAngleWeightedVertexNormal(DArray<Coord>& vertexNormal)
540 {
541 uint vertSize = this->mCoords.size();
542
543 vertexNormal.resize(vertSize);
544
545 auto& vert2Tri = getVertex2Triangles();
546
547 cuExecute(vertSize,
548 TS_SetupAngleWeightedVertexNormals,
549 vertexNormal,
550 this->mCoords,
551 mTriangleIndex,
552 vert2Tri);
553 }
554
555 template<typename Coord, typename Triangle, typename Edg2Tri>
556 __global__ void TS_SetupEdgeNormals(
557 DArray<Coord> normals,
558 DArray<Coord> vertices,
559 DArray<Triangle> triangles,
560 DArray<Edg2Tri> triIds)
561 {
562 int tId = threadIdx.x + (blockIdx.x * blockDim.x);
563 if (tId >= normals.size()) return;
564
565 Edg2Tri& edge = triIds[tId];
566
567 Coord N = Coord(0);
568 for (int ne = 0; ne < 2; ne++)
569 {
570 int j = edge[ne];
571 Triangle t = triangles[j];
572
573 Coord v0 = vertices[t[0]];
574 Coord v1 = vertices[t[1]];
575 Coord v2 = vertices[t[2]];
576
577 Coord norm = (v1 - v0).cross(v2 - v0);
578 norm.normalize();
579 N += norm;
580 }
581
582 N.normalize();
583
584 normals[tId] = N;
585 }
586
587 template<typename TDataType>
588 void TriangleSet<TDataType>::updateEdgeNormal(DArray<Coord>& edgeNormal)
589 {
590 if (mEdg2Tri.size() == 0)
591 updateEdges();
592
593 edgeNormal.resize(mEdg2Tri.size());
594
595 cuExecute(mEdg2Tri.size(),
596 TS_SetupEdgeNormals,
597 edgeNormal,
598 this->mCoords,
599 mTriangleIndex,
600 mEdg2Tri);
601 }
602
603 template<typename TDataType>
604 void TriangleSet<TDataType>::updateTopology()
605 {
606 this->updateTriangles();
607
608 if(bAutoUpdateNormal)
609 this->updateVertexNormal();
610
611 this->EdgeSet<TDataType>::updateTopology();
612 }
613
614 template <typename Real, typename Coord>
615 __global__ void PS_RotateNormal(
616 DArray<Coord> normals,
617 Quat<Real> q)
618 {
619 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
620 if (pId >= normals.size()) return;
621 SquareMatrix<Real, 3> rot = q.toMatrix3x3();
622
623 normals[pId] = rot * normals[pId];
624 }
625
626 template<typename TDataType>
627 void TriangleSet<TDataType>::rotate(const Coord angle)
628 {
629 EdgeSet<TDataType>::rotate(angle);
630 }
631
632 template<typename TDataType>
633 void TriangleSet<TDataType>::rotate(const Quat<Real> q)
634 {
635 EdgeSet<TDataType>::rotate(q);
636
637 cuExecute(mVertexNormal.size(), PS_RotateNormal, mVertexNormal, q);
638 }
639
640 DEFINE_CLASS(TriangleSet);
641}
Generated by doxygen 1.13.2