55 vtkPolyData *
const outputCriticalPoints,
56 vtkPolyData *
const outputSeparatrices1,
57 vtkPolyData *
const outputSeparatrices2,
58 const SimplexId *
const inputOffsets,
59 const triangulationType &triangulation) {
61 const int dimensionality = triangulation.getDimensionality();
62 const auto scalars = ttkUtils::GetPointer<scalarType>(inputScalars);
65 criticalPoints_, separatrices1_, separatrices2_, segmentations_, scalars,
66 inputScalars->GetMTime(), inputOffsets, triangulation);
68#ifndef TTK_ENABLE_KAMIKAZE
70 this->
printErr(
"MorseSmaleComplex.execute() error");
77 vtkNew<vtkPoints> points{};
78 vtkNew<vtkSignedCharArray> cellDimensions{};
79 vtkNew<ttkSimplexIdTypeArray> cellIds{};
81 inputScalars->NewInstance()};
82 vtkNew<vtkSignedCharArray> isOnBoundary{};
83 vtkNew<ttkSimplexIdTypeArray> PLVertexIdentifiers{};
84 vtkNew<ttkSimplexIdTypeArray> manifoldSizeScalars{};
85 const auto nPoints = criticalPoints_.
points_.size();
87#ifndef TTK_ENABLE_KAMIKAZE
88 if(!points || !cellDimensions || !cellIds || !cellScalars || !isOnBoundary
89 || !PLVertexIdentifiers || !manifoldSizeScalars) {
90 this->
printErr(
"Critical points vtkDataArray allocation problem.");
95 points->SetNumberOfPoints(nPoints);
97 cellDimensions->SetNumberOfComponents(1);
101 cellIds->SetNumberOfComponents(1);
105 cellScalars->SetNumberOfComponents(1);
106 cellScalars->SetName(inputScalars->GetName());
107 cellScalars->SetNumberOfTuples(nPoints);
108#ifdef TTK_ENABLE_OPENMP
109#pragma omp parallel for num_threads(this->threadNumber_)
111 for(
size_t i = 0; i < nPoints; ++i) {
112 points->SetPoint(i, criticalPoints_.
points_[i].data());
113 cellScalars->SetTuple1(
117 isOnBoundary->SetNumberOfComponents(1);
121 PLVertexIdentifiers->SetNumberOfComponents(1);
125 manifoldSizeScalars->SetNumberOfComponents(1);
136 auto pointData = outputCriticalPoints->GetPointData();
137#ifndef TTK_ENABLE_KAMIKAZE
139 this->
printErr(
"outputCriticalPoints has no point data.");
144 pointData->SetScalars(cellDimensions);
145 pointData->AddArray(cellIds);
146 pointData->AddArray(cellScalars);
147 pointData->AddArray(isOnBoundary);
148 pointData->AddArray(PLVertexIdentifiers);
149 pointData->AddArray(manifoldSizeScalars);
156 vtkNew<vtkFloatArray> pointsCoords{};
157 vtkNew<vtkSignedCharArray> smoothingMask{};
158 vtkNew<vtkSignedCharArray> cellDimensions{};
159 vtkNew<ttkSimplexIdTypeArray> cellIds{};
160 vtkNew<ttkSimplexIdTypeArray> sourceIds{};
161 vtkNew<ttkSimplexIdTypeArray> destinationIds{};
162 vtkNew<ttkSimplexIdTypeArray> separatrixIds{};
163 vtkNew<vtkSignedCharArray> separatrixTypes{};
164 vtkNew<vtkDoubleArray> separatrixFunctionMaxima{};
165 vtkNew<vtkDoubleArray> separatrixFunctionMinima{};
166 vtkNew<vtkDoubleArray> separatrixFunctionDiffs{};
167 vtkNew<vtkSignedCharArray> isOnBoundary{};
169#ifndef TTK_ENABLE_KAMIKAZE
170 if(!pointsCoords || !smoothingMask || !cellDimensions || !cellIds
171 || !sourceIds || !destinationIds || !separatrixIds || !separatrixTypes
172 || !separatrixFunctionMaxima || !separatrixFunctionMinima
173 || !separatrixFunctionDiffs || !isOnBoundary) {
174 this->
printErr(
"1-separatrices vtkDataArray allocation problem.");
179 pointsCoords->SetNumberOfComponents(3);
182 smoothingMask->SetNumberOfComponents(1);
186 cellDimensions->SetNumberOfComponents(1);
190 cellIds->SetNumberOfComponents(1);
194 sourceIds->SetNumberOfComponents(1);
198 destinationIds->SetNumberOfComponents(1);
202 separatrixIds->SetNumberOfComponents(1);
206 separatrixTypes->SetNumberOfComponents(1);
210 separatrixFunctionMaxima->SetNumberOfComponents(1);
212 separatrixFunctionMaxima->SetNumberOfTuples(
215 separatrixFunctionMinima->SetNumberOfComponents(1);
217 separatrixFunctionMinima->SetNumberOfTuples(
220 separatrixFunctionDiffs->SetNumberOfComponents(1);
222 separatrixFunctionDiffs->SetNumberOfTuples(
225#ifdef TTK_ENABLE_OPENMP
226#pragma omp parallel for num_threads(this->threadNumber_)
233 separatrixFunctionMinima->SetTuple1(i, min);
234 separatrixFunctionMaxima->SetTuple1(i, max);
235 separatrixFunctionDiffs->SetTuple1(i, max - min);
238 isOnBoundary->SetNumberOfComponents(1);
242 vtkNew<ttkSimplexIdTypeArray> offsets{}, connectivity{};
243 offsets->SetNumberOfComponents(1);
245 connectivity->SetNumberOfComponents(1);
248#ifdef TTK_ENABLE_OPENMP
249#pragma omp parallel for num_threads(this->threadNumber_)
252 offsets->SetTuple1(i, 2 * i);
255 vtkNew<vtkPoints> points{};
256 points->SetData(pointsCoords);
257 outputSeparatrices1->SetPoints(points);
258 vtkNew<vtkCellArray> cells{};
259#ifndef TTK_ENABLE_64BIT_IDS
260 cells->Use32BitStorage();
262 cells->SetData(offsets, connectivity);
263 outputSeparatrices1->SetLines(cells);
265 auto pointData = outputSeparatrices1->GetPointData();
266 auto cellData = outputSeparatrices1->GetCellData();
268#ifndef TTK_ENABLE_KAMIKAZE
269 if(!pointData || !cellData) {
270 this->
printErr(
"outputSeparatrices1 has no point or no cell data.");
275 pointData->AddArray(smoothingMask);
276 pointData->AddArray(cellDimensions);
277 pointData->AddArray(cellIds);
279 cellData->AddArray(sourceIds);
280 cellData->AddArray(destinationIds);
281 cellData->AddArray(separatrixIds);
282 cellData->SetScalars(separatrixTypes);
283 cellData->AddArray(separatrixFunctionMaxima);
284 cellData->AddArray(separatrixFunctionMinima);
285 cellData->AddArray(separatrixFunctionDiffs);
286 cellData->AddArray(isOnBoundary);
290 if(dimensionality == 3
293 vtkNew<vtkFloatArray> pointsCoords{};
294 vtkNew<ttkSimplexIdTypeArray> sourceIds{};
295 vtkNew<ttkSimplexIdTypeArray> separatrixIds{};
296 vtkNew<vtkSignedCharArray> separatrixTypes{};
297 vtkNew<vtkDoubleArray> separatrixFunctionMaxima{};
298 vtkNew<vtkDoubleArray> separatrixFunctionMinima{};
299 vtkNew<vtkDoubleArray> separatrixFunctionDiffs{};
300 vtkNew<vtkSignedCharArray> isOnBoundary{};
302#ifndef TTK_ENABLE_KAMIKAZE
303 if(!pointsCoords || !sourceIds || !separatrixIds || !separatrixTypes
304 || !separatrixFunctionMaxima || !separatrixFunctionMinima
305 || !separatrixFunctionDiffs || !isOnBoundary) {
306 this->
printErr(
"2-separatrices vtkDataArray allocation problem.");
311 pointsCoords->SetNumberOfComponents(3);
314 sourceIds->SetNumberOfComponents(1);
318 separatrixIds->SetNumberOfComponents(1);
322 separatrixTypes->SetNumberOfComponents(1);
326 separatrixFunctionMaxima->SetNumberOfComponents(1);
328 separatrixFunctionMaxima->SetNumberOfTuples(
331 separatrixFunctionMinima->SetNumberOfComponents(1);
333 separatrixFunctionMinima->SetNumberOfTuples(
336 separatrixFunctionDiffs->SetNumberOfComponents(1);
338 separatrixFunctionDiffs->SetNumberOfTuples(
341#ifdef TTK_ENABLE_OPENMP
342#pragma omp parallel for num_threads(this->threadNumber_)
349 separatrixFunctionMinima->SetTuple1(i, min);
350 separatrixFunctionMaxima->SetTuple1(i, max);
351 separatrixFunctionDiffs->SetTuple1(i, max - min);
354 isOnBoundary->SetNumberOfComponents(1);
358 vtkNew<ttkSimplexIdTypeArray> offsets{}, connectivity{};
359 offsets->SetNumberOfComponents(1);
361 connectivity->SetNumberOfComponents(1);
364 vtkNew<vtkPoints> points{};
365 points->SetData(pointsCoords);
366 outputSeparatrices2->SetPoints(points);
367 vtkNew<vtkCellArray> cells{};
368#ifndef TTK_ENABLE_64BIT_IDS
369 cells->Use32BitStorage();
371 cells->SetData(offsets, connectivity);
372 outputSeparatrices2->SetPolys(cells);
374 auto cellData = outputSeparatrices2->GetCellData();
375#ifndef TTK_ENABLE_KAMIKAZE
377 this->
printErr(
"outputSeparatrices2 has no cell data.");
382 cellData->AddArray(sourceIds);
383 cellData->AddArray(separatrixIds);
384 cellData->AddArray(separatrixTypes);
385 cellData->AddArray(separatrixFunctionMaxima);
386 cellData->AddArray(separatrixFunctionMinima);
387 cellData->AddArray(separatrixFunctionDiffs);
388 cellData->AddArray(isOnBoundary);
395 vtkInformationVector **inputVector,
396 vtkInformationVector *outputVector) {
398 const auto input = vtkDataSet::GetData(inputVector[0]);
399 auto outputCriticalPoints = vtkPolyData::GetData(outputVector, 0);
400 auto outputSeparatrices1 = vtkPolyData::GetData(outputVector, 1);
401 auto outputSeparatrices2 = vtkPolyData::GetData(outputVector, 2);
402 auto outputMorseComplexes = vtkDataSet::GetData(outputVector, 3);
404#ifndef TTK_ENABLE_KAMIKAZE
406 this->
printErr(
"Input pointer is NULL.");
409 if(input->GetNumberOfPoints() == 0) {
410 this->
printErr(
"Input has no point.");
413 if(!outputCriticalPoints or !outputSeparatrices1 or !outputSeparatrices2
414 or !outputMorseComplexes) {
415 this->
printErr(
"Output pointers are NULL.");
421 if(triangulation ==
nullptr) {
422 this->
printErr(
"Triangulation is null");
427 const auto inputScalars = this->GetInputArrayToProcess(0, inputVector);
429#ifndef TTK_ENABLE_KAMIKAZE
430 if(inputScalars ==
nullptr) {
437 input, 0, triangulation,
false, 1, this->ForceInputOffsetScalarField);
439#ifndef TTK_ENABLE_KAMIKAZE
440 if(inputOffsets ==
nullptr) {
444 if(inputOffsets->GetDataType() != VTK_INT
445 and inputOffsets->GetDataType() != VTK_ID_TYPE) {
446 this->
printErr(
"input offset field type not supported.");
451 this->
printMsg(
"Launching computation on field `"
452 + std::string(inputScalars->GetName()) +
"'...");
455 const SimplexId numberOfVertices = triangulation->getNumberOfVertices();
456#ifndef TTK_ENABLE_KAMIKAZE
457 if(!numberOfVertices) {
458 this->
printErr(
"Input has no vertices.");
463 vtkNew<ttkSimplexIdTypeArray> ascendingManifold{};
464 vtkNew<ttkSimplexIdTypeArray> descendingManifold{};
465 vtkNew<ttkSimplexIdTypeArray> morseSmaleManifold{};
466#ifndef TTK_ENABLE_KAMIKAZE
467 if(!ascendingManifold || !descendingManifold || !morseSmaleManifold) {
468 this->
printErr(
"Manifold vtkDataArray allocation problem.");
472 ascendingManifold->SetNumberOfComponents(1);
473 ascendingManifold->SetNumberOfTuples(numberOfVertices);
476 descendingManifold->SetNumberOfComponents(1);
477 descendingManifold->SetNumberOfTuples(numberOfVertices);
480 morseSmaleManifold->SetNumberOfComponents(1);
481 morseSmaleManifold->SetNumberOfTuples(numberOfVertices);
484 this->segmentations_ = {ttkUtils::GetPointer<SimplexId>(ascendingManifold),
485 ttkUtils::GetPointer<SimplexId>(descendingManifold),
486 ttkUtils::GetPointer<SimplexId>(morseSmaleManifold)};
495 inputScalars->GetDataType(), triangulation->getType(),
496 (ret = dispatch<VTK_TT, TTK_TT>(
497 inputScalars, outputCriticalPoints, outputSeparatrices1,
498 outputSeparatrices2, ttkUtils::GetPointer<SimplexId>(inputOffsets),
499 *
static_cast<TTK_TT *
>(triangulation->getData()))));
505 outputMorseComplexes->ShallowCopy(input);
508 vtkPointData *pointData = outputMorseComplexes->GetPointData();
509#ifndef TTK_ENABLE_KAMIKAZE
511 this->
printErr(
"outputMorseComplexes has no point data.");
516 if(ComputeDescendingSegmentation)
517 pointData->AddArray(descendingManifold);
519 pointData->AddArray(ascendingManifold);
522 pointData->AddArray(morseSmaleManifold);