TTK
Loading...
Searching...
No Matches
ttkTriangulationRequest.cpp
Go to the documentation of this file.
1#include <ttkTriangulationRequest.h>
2#include <ttkUtils.h>
3
4#include <vtkCellData.h>
5#include <vtkDataArray.h>
6#include <vtkDataSet.h>
7#include <vtkInformation.h>
8#include <vtkNew.h>
9#include <vtkPointData.h>
10#include <vtkSignedCharArray.h>
11#include <vtkUnstructuredGrid.h>
12
13vtkStandardNewMacro(ttkTriangulationRequest);
14
15ttkTriangulationRequest::ttkTriangulationRequest() {
16 this->setDebugMsgPrefix("TriangulationRequest");
17 this->SetNumberOfInputPorts(1);
18 this->SetNumberOfOutputPorts(1);
19}
20
21int ttkTriangulationRequest::FillInputPortInformation(int port,
22 vtkInformation *info) {
23 if(port == 0) {
24 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
25 return 1;
26 }
27 return 0;
28}
29
30int ttkTriangulationRequest::FillOutputPortInformation(int port,
31 vtkInformation *info) {
32 if(port == 0) {
33 info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkUnstructuredGrid");
34 return 1;
35 }
36 return 0;
37}
38
39int ttkTriangulationRequest::RequestData(vtkInformation *ttkNotUsed(request),
40 vtkInformationVector **inputVector,
41 vtkInformationVector *outputVector) {
42
43 ttk::Timer timer;
44
45 vtkDataSet *input = vtkDataSet::GetData(inputVector[0]);
46 vtkUnstructuredGrid *output = vtkUnstructuredGrid::GetData(outputVector);
47
48 const auto triangulation = ttkAlgorithm::GetTriangulation(input);
49 if(!triangulation)
50 return 0;
51
52 const int dimensionality = triangulation->getDimensionality();
53
54 vtkNew<vtkPoints> points{};
55 vtkNew<vtkUnstructuredGrid> cells{};
56 vtkNew<ttkSimplexIdTypeArray> cellIds{};
57 cellIds->SetNumberOfComponents(1);
58 cellIds->SetName("CellId");
59 vtkNew<vtkSignedCharArray> cellDims{};
60 cellDims->SetNumberOfComponents(1);
61 cellDims->SetName("CellDimension");
62
63 using ttk::SimplexId;
64
65 std::vector<SimplexId> vertices;
66 const SimplexId numberOfVertices = triangulation->getNumberOfVertices();
67 std::vector<SimplexId> isVisited(numberOfVertices, -1);
68
69 this->printMsg(ttk::debug::Separator::L1);
70 this->printMsg({
71 {"#Threads", std::to_string(this->threadNumber_)},
72 {"#Vertices", std::to_string(numberOfVertices)},
73 });
74
75 auto addVertex = [&](const SimplexId vertexId) {
76 if(vertexId == -1)
77 return vtkIdType(-1);
78
79 std::array<float, 3> p{};
80 triangulation->getVertexPoint(vertexId, p[0], p[1], p[2]);
81 vertices.push_back(vertexId);
82 return points->InsertNextPoint(p.data());
83 };
84
85 auto addEdge = [&](const SimplexId edgeId) {
86 std::array<vtkIdType, 2> pointIds{};
87 std::array<SimplexId, 2> vertexIds{};
88
89 for(int i = 0; i < 2; ++i) {
90 triangulation->getEdgeVertex(edgeId, i, vertexIds[i]);
91
92 const SimplexId vertexId = vertexIds[i];
93 if(vertexId == -1)
94 return -1;
95
96 if(isVisited[vertexId] == -1) {
97 pointIds[i] = addVertex(vertexId);
98 isVisited[vertexId] = pointIds[i];
99 } else
100 pointIds[i] = isVisited[vertexId];
101 }
102
103 cells->InsertNextCell(VTK_LINE, 2, pointIds.data());
104 cellIds->InsertNextTuple1(edgeId);
105 cellDims->InsertNextTuple1(1);
106
107 return 0;
108 };
109
110 auto addTriangle = [&](const SimplexId triangleId) {
111 std::array<vtkIdType, 3> pointIds{};
112 std::array<SimplexId, 3> vertexIds{};
113
114 for(int i = 0; i < 3; ++i) {
115 if(dimensionality == 3)
116 triangulation->getTriangleVertex(triangleId, i, vertexIds[i]);
117 else
118 triangulation->getCellVertex(triangleId, i, vertexIds[i]);
119
120 const SimplexId vertexId = vertexIds[i];
121 if(vertexId == -1)
122 return -1;
123
124 if(isVisited[vertexId] == -1) {
125 pointIds[i] = addVertex(vertexId);
126 isVisited[vertexId] = pointIds[i];
127 } else
128 pointIds[i] = isVisited[vertexId];
129 }
130
131 cells->InsertNextCell(VTK_TRIANGLE, 3, pointIds.data());
132 cellIds->InsertNextTuple1(triangleId);
133 cellDims->InsertNextTuple1(2);
134
135 return 0;
136 };
137
138 auto addTetra = [&](const SimplexId tetraId) {
139 std::array<vtkIdType, 4> pointIds{};
140 std::array<SimplexId, 4> vertexIds{};
141
142 for(int i = 0; i < 4; ++i) {
143 triangulation->getCellVertex(tetraId, i, vertexIds[i]);
144
145 const SimplexId vertexId = vertexIds[i];
146 if(vertexId == -1)
147 return -1;
148
149 if(isVisited[vertexId] == -1) {
150 pointIds[i] = addVertex(vertexId);
151 isVisited[vertexId] = pointIds[i];
152 } else
153 pointIds[i] = isVisited[vertexId];
154 }
155
156 cells->InsertNextCell(VTK_TETRA, 4, pointIds.data());
157 cellIds->InsertNextTuple1(tetraId);
158 cellDims->InsertNextTuple1(3);
159
160 return 0;
161 };
162
163 auto addStar = [&](const SimplexId starId) {
164 if(dimensionality == 2)
165 addTriangle(starId);
166 else if(dimensionality == 3)
167 addTetra(starId);
168 };
169
170 // parse SimplexIdentifier into a vector of identifiers
171 std::vector<SimplexId> ids{};
172 std::istringstream iss(this->SimplexIdentifier);
173 for(SimplexId i; iss >> i;) {
174 ids.emplace_back(i);
175 if(iss.peek() == ',') {
176 iss.ignore();
177 }
178 }
179
180 // remove duplicates and negative values
181 {
182 std::sort(ids.rbegin(), ids.rend()); // decreasing order
183 const auto last = std::unique(ids.begin(), ids.end());
184 ids.erase(last, ids.end());
185
186 // first negative value
187 const auto it = std::find_if(
188 ids.begin(), ids.end(), [](const SimplexId a) { return a < 0; });
189 // remove negative values
190 ids.erase(it, ids.end());
191 }
192
193 const auto detectOutOfBounds = [this, numberOfVertices, dimensionality,
194 &triangulation](const SimplexId si) {
195 switch(this->SimplexType) {
196 case SIMPLEX::VERTEX:
197 if(si >= numberOfVertices) {
198 this->printWrn("Vertex ID " + std::to_string(si)
199 + " out of bounds (max. "
200 + std::to_string(numberOfVertices) + ").");
201 return true;
202 }
203 break;
204
205 case SIMPLEX::EDGE:
206 triangulation->preconditionEdges();
207 if(si >= triangulation->getNumberOfEdges()) {
208 this->printWrn(
209 "Edge ID " + std::to_string(si) + " out of bounds (max. "
210 + std::to_string(triangulation->getNumberOfEdges()) + ").");
211 return true;
212 }
213 break;
214
215 case SIMPLEX::TRIANGLE: {
216 if(dimensionality == 3) {
217 triangulation->preconditionTriangles();
218 }
219 const auto numberOfTriangles
220 = dimensionality == 2 ? triangulation->getNumberOfCells()
221 : triangulation->getNumberOfTriangles();
222 if(si >= numberOfTriangles) {
223 this->printWrn("Triangle ID " + std::to_string(si)
224 + " out of bounds (max. "
225 + std::to_string(numberOfTriangles) + ").");
226 return true;
227 }
228 } break;
229
230 case SIMPLEX::TETRA:
231 if(dimensionality == 3)
232 if(si >= triangulation->getNumberOfCells()) {
233 this->printWrn(
234 "Tetrahedron ID " + std::to_string(si) + " out of bounds (max. "
235 + std::to_string(triangulation->getNumberOfCells()) + ").");
236 return true;
237 }
238 break;
239 }
240 return false;
241 };
242
243 // remove out-of-bounds values
244 {
245 const auto last = std::remove_if(ids.begin(), ids.end(), detectOutOfBounds);
246 ids.erase(last, ids.end());
247 }
248
249 // sanity check
250 if(ids.empty()) {
251 this->printErr("Invalid simplex indices");
252 return 0;
253 }
254
255 // do minimum preprocess and put watchdog on SimplexIdentifier
256 for(const auto si : ids) {
257 switch(this->RequestType) {
258 case REQUEST::COMPUTE_SIMPLEX:
259 switch(this->SimplexType) {
260 case SIMPLEX::VERTEX: {
261 const auto vid = addVertex(si);
262 cells->InsertNextCell(VTK_VERTEX, 1, &vid);
263 cellIds->InsertNextTuple1(vid);
264 cellDims->InsertNextTuple1(0);
265 } break;
266
267 case SIMPLEX::EDGE:
268 addEdge(si);
269 break;
270
271 case SIMPLEX::TRIANGLE:
272 addTriangle(si);
273 break;
274
275 case SIMPLEX::TETRA:
276 if(dimensionality == 3)
277 addTetra(si);
278 break;
279 }
280 break;
281
282 case REQUEST::COMPUTE_FACET:
283 switch(this->SimplexType) {
284 case SIMPLEX::VERTEX:
285 break;
286
287 case SIMPLEX::EDGE:
288 for(int i = 0; i < 2; ++i) {
289 SimplexId vertexId;
290 triangulation->getEdgeVertex(si, i, vertexId);
291 addVertex(vertexId);
292 }
293 break;
294
295 case SIMPLEX::TRIANGLE:
296 if(dimensionality == 2) {
297 triangulation->preconditionCellEdges();
298 for(int i = 0; i < 3; ++i) {
299 SimplexId edgeId;
300 triangulation->getCellEdge(si, i, edgeId);
301 addEdge(edgeId);
302 }
303 } else if(dimensionality == 3) {
304 triangulation->preconditionTriangleEdges();
305 for(int i = 0; i < 3; ++i) {
306 SimplexId edgeId;
307 triangulation->getTriangleEdge(si, i, edgeId);
308 addEdge(edgeId);
309 }
310 }
311 break;
312
313 case SIMPLEX::TETRA:
314 if(dimensionality == 3) {
315 triangulation->preconditionCellTriangles();
316 for(int i = 0; i < 4; ++i) {
317 SimplexId triangleId;
318 triangulation->getCellTriangle(si, i, triangleId);
319 addTriangle(triangleId);
320 }
321 }
322 break;
323 }
324 break;
325
326 case REQUEST::COMPUTE_COFACET:
327 switch(this->SimplexType) {
328 case SIMPLEX::VERTEX:
329 triangulation->preconditionVertexNeighbors();
330 triangulation->preconditionVertexEdges();
331 {
332 const SimplexId edgeNumber
333 = triangulation->getVertexEdgeNumber(si);
334 for(SimplexId i = 0; i < edgeNumber; ++i) {
335 SimplexId edgeId;
336 triangulation->getVertexEdge(si, i, edgeId);
337 addEdge(edgeId);
338 }
339 }
340 break;
341
342 case SIMPLEX::EDGE:
343 if(dimensionality == 2) {
344 triangulation->preconditionEdgeStars();
345 const SimplexId starNumber = triangulation->getEdgeStarNumber(si);
346 for(SimplexId i = 0; i < starNumber; ++i) {
347 SimplexId starId;
348 triangulation->getEdgeStar(si, i, starId);
349 addStar(starId);
350 }
351 } else if(dimensionality == 3) {
352 triangulation->preconditionEdgeTriangles();
353 const SimplexId triangleNumber
354 = triangulation->getEdgeTriangleNumber(si);
355 for(SimplexId i = 0; i < triangleNumber; ++i) {
356 SimplexId triangleId;
357 triangulation->getEdgeTriangle(si, i, triangleId);
358 addTriangle(triangleId);
359 }
360 }
361 break;
362
363 case SIMPLEX::TRIANGLE:
364 if(dimensionality == 3) {
365 triangulation->preconditionTriangleStars();
366 const SimplexId starNumber
367 = triangulation->getTriangleStarNumber(si);
368 for(SimplexId i = 0; i < starNumber; ++i) {
369 SimplexId starId;
370 triangulation->getTriangleStar(si, i, starId);
371 addStar(starId);
372 }
373 }
374 break;
375
376 case SIMPLEX::TETRA:
377 break;
378 }
379 break;
380
381 case REQUEST::COMPUTE_STAR:
382 switch(this->SimplexType) {
383 case SIMPLEX::VERTEX:
384 triangulation->preconditionVertexStars();
385 {
386 const SimplexId starNumber
387 = triangulation->getVertexStarNumber(si);
388 for(SimplexId i = 0; i < starNumber; ++i) {
389 SimplexId starId;
390 triangulation->getVertexStar(si, i, starId);
391 addStar(starId);
392 }
393 }
394 break;
395
396 case SIMPLEX::EDGE:
397 triangulation->preconditionEdgeStars();
398 {
399 const SimplexId starNumber = triangulation->getEdgeStarNumber(si);
400 for(SimplexId i = 0; i < starNumber; ++i) {
401 SimplexId starId;
402 triangulation->getEdgeStar(si, i, starId);
403 addStar(starId);
404 }
405 }
406 break;
407
408 case SIMPLEX::TRIANGLE:
409 if(dimensionality == 3) {
410 triangulation->preconditionTriangleStars();
411 const SimplexId starNumber
412 = triangulation->getTriangleStarNumber(si);
413 for(SimplexId i = 0; i < starNumber; ++i) {
414 SimplexId starId;
415 triangulation->getTriangleStar(si, i, starId);
416 addStar(starId);
417 }
418 }
419 break;
420
421 case SIMPLEX::TETRA:
422 break;
423 }
424 break;
425
426 case REQUEST::COMPUTE_LINK:
427 switch(this->SimplexType) {
428 case SIMPLEX::VERTEX:
429 triangulation->preconditionVertexLinks();
430 if(dimensionality == 2) {
431 triangulation->preconditionEdges();
432 const SimplexId linkNumber
433 = triangulation->getVertexLinkNumber(si);
434 for(SimplexId i = 0; i < linkNumber; ++i) {
435 SimplexId linkId;
436 triangulation->getVertexLink(si, i, linkId);
437 addEdge(linkId);
438 }
439 } else if(dimensionality == 3) {
440 triangulation->preconditionVertexTriangles();
441 const SimplexId linkNumber
442 = triangulation->getVertexLinkNumber(si);
443 for(SimplexId i = 0; i < linkNumber; ++i) {
444 SimplexId linkId;
445 triangulation->getVertexLink(si, i, linkId);
446 addTriangle(linkId);
447 }
448 }
449 break;
450
451 case SIMPLEX::EDGE:
452 triangulation->preconditionEdgeLinks();
453 if(dimensionality == 2) {
454 const SimplexId linkNumber = triangulation->getEdgeLinkNumber(si);
455 for(SimplexId i = 0; i < linkNumber; ++i) {
456 SimplexId linkId;
457 triangulation->getEdgeLink(si, i, linkId);
458 addVertex(linkId);
459 }
460 } else if(dimensionality == 3) {
461 const SimplexId linkNumber = triangulation->getEdgeLinkNumber(si);
462 for(SimplexId i = 0; i < linkNumber; ++i) {
463 SimplexId linkId;
464 triangulation->getEdgeLink(si, i, linkId);
465 addEdge(linkId);
466 }
467 }
468 break;
469
470 case SIMPLEX::TRIANGLE:
471 if(dimensionality == 3) {
472 triangulation->preconditionTriangleLinks();
473 const SimplexId linkNumber
474 = triangulation->getTriangleLinkNumber(si);
475 for(SimplexId i = 0; i < linkNumber; ++i) {
476 SimplexId linkId;
477 triangulation->getTriangleLink(si, i, linkId);
478 addVertex(linkId);
479 }
480 }
481 break;
482
483 case SIMPLEX::TETRA:
484 break;
485 }
486 break;
487 case REQUEST::COMPUTE_BOUNDARY:
488 switch(this->SimplexType) {
489 case SIMPLEX::VERTEX:
490 triangulation->preconditionBoundaryVertices();
491 for(SimplexId v = 0; v < triangulation->getNumberOfVertices();
492 ++v) {
493 if(triangulation->isVertexOnBoundary(v)) {
494 const auto vid = addVertex(v);
495 cells->InsertNextCell(VTK_VERTEX, 1, &vid);
496 cellIds->InsertNextTuple1(vid);
497 cellDims->InsertNextTuple1(0);
498 }
499 }
500 break;
501 case SIMPLEX::EDGE:
502 if(dimensionality == 1) {
503 triangulation->preconditionBoundaryVertices();
504 triangulation->preconditionVertexStars();
505 for(SimplexId v = 0; v < triangulation->getNumberOfVertices();
506 ++v) {
507 if(triangulation->isVertexOnBoundary(v)) {
508 SimplexId e{-1};
509 triangulation->getVertexStar(v, 0, e);
510 addEdge(e);
511 }
512 }
513 } else {
514 triangulation->preconditionBoundaryEdges();
515 for(SimplexId e = 0; e < triangulation->getNumberOfEdges(); ++e) {
516 if(triangulation->isEdgeOnBoundary(e)) {
517 addEdge(e);
518 }
519 }
520 }
521 break;
522 case SIMPLEX::TRIANGLE:
523 if(dimensionality == 2) {
524 triangulation->preconditionBoundaryEdges();
525 triangulation->preconditionEdgeStars();
526 for(SimplexId e = 0; e < triangulation->getNumberOfEdges(); ++e) {
527 if(triangulation->isEdgeOnBoundary(e)) {
528 SimplexId t{-1};
529 triangulation->getEdgeStar(e, 0, t);
530 addTriangle(t);
531 }
532 }
533 } else if(dimensionality == 3) {
534 triangulation->preconditionBoundaryTriangles();
535 for(SimplexId t = 0; t < triangulation->getNumberOfTriangles();
536 ++t) {
537 if(triangulation->isTriangleOnBoundary(t)) {
538 addTriangle(t);
539 }
540 }
541 } else {
542 this->printErr("No triangles on the boundary of 1D data-sets");
543 }
544 break;
545 case SIMPLEX::TETRA:
546 if(dimensionality == 3) {
547 triangulation->preconditionBoundaryTriangles();
548 triangulation->preconditionTriangleStars();
549 for(SimplexId t = 0; t < triangulation->getNumberOfTriangles();
550 ++t) {
551 if(triangulation->isTriangleOnBoundary(t)) {
552 SimplexId T{-1};
553 triangulation->getTriangleStar(t, 0, T);
554 addTetra(T);
555 }
556 }
557 } else {
558 this->printErr(
559 "No tetrahedron on the boundary of 1D or 2D data-sets");
560 }
561 break;
562 }
563 break;
564 }
565 }
566
567 cells->SetPoints(points);
568 cells->GetCellData()->AddArray(cellIds);
569 cells->GetCellData()->AddArray(cellDims);
570
571 output->ShallowCopy(cells);
572
573 if(KeepAllDataArrays) {
574 vtkPointData *inputPointData = input->GetPointData();
575#ifndef TTK_ENABLE_KAMIKAZE
576 if(!inputPointData) {
577 this->printErr("Input has no point data.");
578 return 0;
579 }
580#endif
581 const int numberOfInputArrays = inputPointData->GetNumberOfArrays();
582
583 vtkPointData *outputPointData = output->GetPointData();
584#ifndef TTK_ENABLE_KAMIKAZE
585 if(!outputPointData) {
586 this->printErr("Output has no point data.");
587 return 0;
588 }
589#endif
590
591 for(int i = 0; i < numberOfInputArrays; ++i) {
592 vtkDataArray *arr = inputPointData->GetArray(i);
593
594 if(arr and arr->GetNumberOfComponents() == 1) {
595 vtkDataArray *newArr = arr->NewInstance();
596 if(newArr == nullptr) {
597 continue;
598 }
599 newArr->SetName(arr->GetName());
600 newArr->SetNumberOfComponents(1);
601
602 for(SimplexId const v : vertices)
603 newArr->InsertNextTuple1(arr->GetTuple1(v));
604
605 outputPointData->AddArray(newArr);
606 }
607 }
608 }
609
610 {
611 this->printMsg("Complete", 1, timer.getElapsedTime());
612 this->printMsg(ttk::debug::Separator::L1);
613 }
614 return 1;
615}
ttkNotUsed
#define ttkNotUsed(x)
Mark function/method parameters that are not used in the function body at all.
Definition BaseClass.h:47
ttkAlgorithm::GetTriangulation
ttk::Triangulation * GetTriangulation(vtkDataSet *dataSet)
Definition ttkAlgorithm.cpp:41
ttkTriangulationRequest
TTK VTK-filter that wraps the triangulation processing package.
Definition ttkTriangulationRequest.h:33
ttkTriangulationRequest::FillOutputPortInformation
int FillOutputPortInformation(int port, vtkInformation *info) override
Definition ttkTriangulationRequest.cpp:30
ttkTriangulationRequest::RequestData
int RequestData(vtkInformation *request, vtkInformationVector **inputVector, vtkInformationVector *outputVector) override
Definition ttkTriangulationRequest.cpp:39
ttkTriangulationRequest::FillInputPortInformation
int FillInputPortInformation(int port, vtkInformation *info) override
Definition ttkTriangulationRequest.cpp:21
ttk::Debug::printWrn
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:159
ttk::Debug::setDebugMsgPrefix
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
ttk::Debug::printErr
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition Debug.h:149
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
std::to_string
std::string to_string(__int128)
Definition ripserpy.cpp:99
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
vtkStandardNewMacro
vtkStandardNewMacro(ttkTriangulationRequest)
printMsg
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Mon May 6 2024 07:44:31 for TTK by doxygen 1.9.8