59 vtkPolyData *
const outputCriticalPoints,
60 vtkPolyData *
const outputSeparatrices1,
61 vtkPolyData *
const outputSeparatrices2,
63 const triangulationType &triangulation) {
65 const int dimensionality = triangulation.getDimensionality();
72 = this->
execute(criticalPoints_, separatrices1_, separatrices2_,
73 segmentations_, scalars, inputScalars->GetMTime(),
74 inputOffsets, triangulation, StochasticGradientSeed);
76#ifndef TTK_ENABLE_KAMIKAZE
78 this->
printErr(
"MorseSmaleComplex.execute() error");
85 vtkNew<vtkPoints> points{};
86 vtkNew<vtkSignedCharArray> cellDimensions{};
87 vtkNew<ttkSimplexIdTypeArray> cellIds{};
89 inputScalars->NewInstance()};
90 vtkNew<vtkSignedCharArray> isOnBoundary{};
91 vtkNew<ttkSimplexIdTypeArray> PLVertexIdentifiers{};
92 vtkNew<ttkSimplexIdTypeArray> manifoldSizeScalars{};
93 const auto nPoints = criticalPoints_.points_.size();
95#ifndef TTK_ENABLE_KAMIKAZE
96 if(!points || !cellDimensions || !cellIds || !cellScalars || !isOnBoundary
97 || !PLVertexIdentifiers || !manifoldSizeScalars) {
98 this->
printErr(
"Critical points vtkDataArray allocation problem.");
103 points->SetNumberOfPoints(nPoints);
105 cellDimensions->SetNumberOfComponents(1);
107 setArray(cellDimensions, criticalPoints_.cellDimensions_);
109 cellIds->SetNumberOfComponents(1);
111 setArray(cellIds, criticalPoints_.cellIds_);
113 cellScalars->SetNumberOfComponents(1);
114 cellScalars->SetName(inputScalars->GetName());
115 cellScalars->SetNumberOfTuples(nPoints);
116#ifdef TTK_ENABLE_OPENMP
117#pragma omp parallel for num_threads(this->threadNumber_)
119 for(
size_t i = 0; i < nPoints; ++i) {
120 points->SetPoint(i, criticalPoints_.points_[i].data());
121 cellScalars->SetTuple1(
122 i, scalars[criticalPoints_.PLVertexIdentifiers_[i]]);
125 isOnBoundary->SetNumberOfComponents(1);
127 setArray(isOnBoundary, criticalPoints_.isOnBoundary_);
129 PLVertexIdentifiers->SetNumberOfComponents(1);
131 setArray(PLVertexIdentifiers, criticalPoints_.PLVertexIdentifiers_);
133 manifoldSizeScalars->SetNumberOfComponents(1);
136 criticalPoints_.manifoldSize_.resize(nPoints);
137 std::fill(criticalPoints_.manifoldSize_.begin(),
138 criticalPoints_.manifoldSize_.end(), -1);
140 setArray(manifoldSizeScalars, criticalPoints_.manifoldSize_);
144 auto pointData = outputCriticalPoints->GetPointData();
145#ifndef TTK_ENABLE_KAMIKAZE
147 this->
printErr(
"outputCriticalPoints has no point data.");
152 pointData->SetScalars(cellDimensions);
153 pointData->AddArray(cellIds);
154 pointData->AddArray(cellScalars);
155 pointData->AddArray(isOnBoundary);
156 pointData->AddArray(PLVertexIdentifiers);
157 pointData->AddArray(manifoldSizeScalars);
163 vtkNew<vtkFloatArray> pointsCoords{};
164 vtkNew<vtkSignedCharArray> smoothingMask{};
165 vtkNew<vtkSignedCharArray> cellDimensions{};
166 vtkNew<ttkSimplexIdTypeArray> cellIds{};
167 vtkNew<ttkSimplexIdTypeArray> sourceIds{};
168 vtkNew<ttkSimplexIdTypeArray> destinationIds{};
169 vtkNew<ttkSimplexIdTypeArray> separatrixIds{};
170 vtkNew<vtkSignedCharArray> separatrixTypes{};
171 vtkNew<vtkDoubleArray> separatrixFunctionMaxima{};
172 vtkNew<vtkDoubleArray> separatrixFunctionMinima{};
173 vtkNew<vtkDoubleArray> separatrixFunctionDiffs{};
174 vtkNew<vtkSignedCharArray> isOnBoundary{};
176#ifndef TTK_ENABLE_KAMIKAZE
177 if(!pointsCoords || !smoothingMask || !cellDimensions || !cellIds
178 || !sourceIds || !destinationIds || !separatrixIds || !separatrixTypes
179 || !separatrixFunctionMaxima || !separatrixFunctionMinima
180 || !separatrixFunctionDiffs || !isOnBoundary) {
181 this->
printErr(
"1-separatrices vtkDataArray allocation problem.");
186 pointsCoords->SetNumberOfComponents(3);
187 pointsCoords->SetNumberOfTuples(separatrices1_.pt.numberOfPoints_);
188 for(
int i = 0; i < separatrices1_.pt.numberOfPoints_; i++) {
189 pointsCoords->SetTuple3(i, separatrices1_.pt.points_[3 * i],
190 separatrices1_.pt.points_[3 * i + 1],
191 separatrices1_.pt.points_[3 * i + 2]);
194 smoothingMask->SetNumberOfComponents(1);
196 setArray(smoothingMask, separatrices1_.pt.smoothingMask_);
198 cellDimensions->SetNumberOfComponents(1);
200 setArray(cellDimensions, separatrices1_.pt.cellDimensions_);
202 cellIds->SetNumberOfComponents(1);
204 setArray(cellIds, separatrices1_.pt.cellIds_);
206 sourceIds->SetNumberOfComponents(1);
208 setArray(sourceIds, separatrices1_.cl.sourceIds_);
210 destinationIds->SetNumberOfComponents(1);
212 setArray(destinationIds, separatrices1_.cl.destinationIds_);
214 separatrixIds->SetNumberOfComponents(1);
216 setArray(separatrixIds, separatrices1_.cl.separatrixIds_);
218 separatrixTypes->SetNumberOfComponents(1);
220 setArray(separatrixTypes, separatrices1_.cl.separatrixTypes_);
222 separatrixFunctionMaxima->SetNumberOfComponents(1);
224 separatrixFunctionMaxima->SetNumberOfTuples(
225 separatrices1_.cl.numberOfCells_);
227 separatrixFunctionMinima->SetNumberOfComponents(1);
229 separatrixFunctionMinima->SetNumberOfTuples(
230 separatrices1_.cl.numberOfCells_);
232 separatrixFunctionDiffs->SetNumberOfComponents(1);
234 separatrixFunctionDiffs->SetNumberOfTuples(
235 separatrices1_.cl.numberOfCells_);
237#ifdef TTK_ENABLE_OPENMP
238#pragma omp parallel for num_threads(this->threadNumber_)
240 for(
SimplexId i = 0; i < separatrices1_.cl.numberOfCells_; ++i) {
241 const auto sepId = separatrices1_.cl.separatrixIds_[i];
243 const auto min = scalars[separatrices1_.cl.sepFuncMinId_[sepId]];
244 const auto max = scalars[separatrices1_.cl.sepFuncMaxId_[sepId]];
245 separatrixFunctionMinima->SetTuple1(i, min);
246 separatrixFunctionMaxima->SetTuple1(i, max);
247 separatrixFunctionDiffs->SetTuple1(i, max - min);
250 isOnBoundary->SetNumberOfComponents(1);
252 setArray(isOnBoundary, separatrices1_.cl.isOnBoundary_);
254 vtkNew<ttkSimplexIdTypeArray> offsets{}, connectivity{};
255 offsets->SetNumberOfComponents(1);
256 offsets->SetNumberOfTuples(separatrices1_.cl.numberOfCells_ + 1);
257 connectivity->SetNumberOfComponents(1);
258 setArray(connectivity, separatrices1_.cl.connectivity_);
260#ifdef TTK_ENABLE_OPENMP
261#pragma omp parallel for num_threads(this->threadNumber_)
263 for(
SimplexId i = 0; i < separatrices1_.cl.numberOfCells_ + 1; ++i) {
264 offsets->SetTuple1(i, 2 * i);
267 vtkNew<vtkPoints> points{};
268 points->SetData(pointsCoords);
269 outputSeparatrices1->SetPoints(points);
270 vtkNew<vtkCellArray> cells{};
271#ifndef TTK_ENABLE_64BIT_IDS
272 cells->Use32BitStorage();
274 cells->SetData(offsets, connectivity);
275 outputSeparatrices1->SetLines(cells);
277 auto pointData = outputSeparatrices1->GetPointData();
278 auto cellData = outputSeparatrices1->GetCellData();
280#ifndef TTK_ENABLE_KAMIKAZE
281 if(!pointData || !cellData) {
282 this->
printErr(
"outputSeparatrices1 has no point or no cell data.");
287 pointData->AddArray(smoothingMask);
288 pointData->AddArray(cellDimensions);
289 pointData->AddArray(cellIds);
291 cellData->AddArray(sourceIds);
292 cellData->AddArray(destinationIds);
293 cellData->AddArray(separatrixIds);
294 cellData->SetScalars(separatrixTypes);
295 cellData->AddArray(separatrixFunctionMaxima);
296 cellData->AddArray(separatrixFunctionMinima);
297 cellData->AddArray(separatrixFunctionDiffs);
298 cellData->AddArray(isOnBoundary);
302 if(dimensionality == 3
305 vtkNew<vtkFloatArray> pointsCoords{};
306 vtkNew<ttkSimplexIdTypeArray> sourceIds{};
307 vtkNew<ttkSimplexIdTypeArray> separatrixIds{};
308 vtkNew<vtkSignedCharArray> separatrixTypes{};
309 vtkNew<vtkDoubleArray> separatrixFunctionMaxima{};
310 vtkNew<vtkDoubleArray> separatrixFunctionMinima{};
311 vtkNew<vtkDoubleArray> separatrixFunctionDiffs{};
312 vtkNew<vtkSignedCharArray> isOnBoundary{};
314#ifndef TTK_ENABLE_KAMIKAZE
315 if(!pointsCoords || !sourceIds || !separatrixIds || !separatrixTypes
316 || !separatrixFunctionMaxima || !separatrixFunctionMinima
317 || !separatrixFunctionDiffs || !isOnBoundary) {
318 this->
printErr(
"2-separatrices vtkDataArray allocation problem.");
323 pointsCoords->SetNumberOfComponents(3);
324 pointsCoords->SetNumberOfTuples(separatrices2_.pt.points_.size());
325 for(
int i = 0; i < separatrices2_.pt.numberOfPoints_; i++) {
326 pointsCoords->SetTuple3(i, separatrices2_.pt.points_[3 * i],
327 separatrices2_.pt.points_[3 * i + 1],
328 separatrices2_.pt.points_[3 * i + 2]);
331 sourceIds->SetNumberOfComponents(1);
333 setArray(sourceIds, separatrices2_.cl.sourceIds_);
335 separatrixIds->SetNumberOfComponents(1);
337 setArray(separatrixIds, separatrices2_.cl.separatrixIds_);
339 separatrixTypes->SetNumberOfComponents(1);
341 setArray(separatrixTypes, separatrices2_.cl.separatrixTypes_);
343 separatrixFunctionMaxima->SetNumberOfComponents(1);
345 separatrixFunctionMaxima->SetNumberOfTuples(
346 separatrices2_.cl.numberOfCells_);
348 separatrixFunctionMinima->SetNumberOfComponents(1);
350 separatrixFunctionMinima->SetNumberOfTuples(
351 separatrices2_.cl.numberOfCells_);
353 separatrixFunctionDiffs->SetNumberOfComponents(1);
355 separatrixFunctionDiffs->SetNumberOfTuples(
356 separatrices2_.cl.numberOfCells_);
358#ifdef TTK_ENABLE_OPENMP
359#pragma omp parallel for num_threads(this->threadNumber_)
361 for(
SimplexId i = 0; i < separatrices2_.cl.numberOfCells_; ++i) {
362 const auto sepId = separatrices2_.cl.separatrixIds_[i];
364 const auto min = scalars[separatrices2_.cl.sepFuncMinId_[sepId]];
365 const auto max = scalars[separatrices2_.cl.sepFuncMaxId_[sepId]];
366 separatrixFunctionMinima->SetTuple1(i, min);
367 separatrixFunctionMaxima->SetTuple1(i, max);
368 separatrixFunctionDiffs->SetTuple1(i, max - min);
371 isOnBoundary->SetNumberOfComponents(1);
373 setArray(isOnBoundary, separatrices2_.cl.isOnBoundary_);
375 vtkNew<ttkSimplexIdTypeArray> offsets{}, connectivity{};
376 offsets->SetNumberOfComponents(1);
377 setArray(offsets, separatrices2_.cl.offsets_);
378 connectivity->SetNumberOfComponents(1);
379 setArray(connectivity, separatrices2_.cl.connectivity_);
381 vtkNew<vtkPoints> points{};
382 points->SetData(pointsCoords);
383 outputSeparatrices2->SetPoints(points);
384 vtkNew<vtkCellArray> cells{};
385#ifndef TTK_ENABLE_64BIT_IDS
386 cells->Use32BitStorage();
388 cells->SetData(offsets, connectivity);
389 outputSeparatrices2->SetPolys(cells);
391 auto cellData = outputSeparatrices2->GetCellData();
392#ifndef TTK_ENABLE_KAMIKAZE
394 this->
printErr(
"outputSeparatrices2 has no cell data.");
399 cellData->AddArray(sourceIds);
400 cellData->AddArray(separatrixIds);
401 cellData->AddArray(separatrixTypes);
402 cellData->AddArray(separatrixFunctionMaxima);
403 cellData->AddArray(separatrixFunctionMinima);
404 cellData->AddArray(separatrixFunctionDiffs);
405 cellData->AddArray(isOnBoundary);
411 vtkInformationVector **inputVector,
412 vtkInformationVector *outputVector) {
415 = vtkDataSet::SafeDownCast(vtkDataSet::GetData(inputVector[0]));
416 auto outputCriticalPoints = vtkPolyData::GetData(outputVector, 0);
417 auto outputSeparatrices1 = vtkPolyData::GetData(outputVector, 1);
418 auto outputSeparatrices2 = vtkPolyData::GetData(outputVector, 2);
419 auto outputMorseComplexes = vtkDataSet::GetData(outputVector, 3);
421#ifndef TTK_ENABLE_KAMIKAZE
423 this->
printErr(
"Input pointer is NULL.");
426 if(input->GetNumberOfPoints() == 0) {
427 this->
printErr(
"Input has no point.");
430 if(!outputCriticalPoints or !outputSeparatrices1 or !outputSeparatrices2
431 or !outputMorseComplexes) {
432 this->
printErr(
"Output pointers are NULL.");
438 if(triangulation ==
nullptr) {
439 this->
printErr(
"Triangulation is null");
444 const auto inputScalars = this->GetInputArrayToProcess(0, inputVector);
446#ifndef TTK_ENABLE_KAMIKAZE
447 if(inputScalars ==
nullptr) {
454 input, 0, triangulation,
false, 1, this->ForceInputOffsetScalarField);
456#ifndef TTK_ENABLE_KAMIKAZE
457 if(inputOffsets ==
nullptr) {
461 if(inputOffsets->GetDataType() != VTK_INT
462 and inputOffsets->GetDataType() != VTK_ID_TYPE) {
463 this->
printErr(
"input offset field type not supported.");
468 this->
printMsg(
"Launching computation on field `"
469 + std::string(inputScalars->GetName()) +
"'...");
472 const SimplexId numberOfVertices = triangulation->getNumberOfVertices();
473#ifndef TTK_ENABLE_KAMIKAZE
474 if(!numberOfVertices) {
475 this->
printErr(
"Input has no vertices.");
480 vtkNew<ttkSimplexIdTypeArray> ascendingManifold{};
481 vtkNew<ttkSimplexIdTypeArray> descendingManifold{};
482 vtkNew<ttkSimplexIdTypeArray> morseSmaleManifold{};
483#ifndef TTK_ENABLE_KAMIKAZE
484 if(!ascendingManifold || !descendingManifold || !morseSmaleManifold) {
485 this->
printErr(
"Manifold vtkDataArray allocation problem.");
490 ascendingManifold->SetNumberOfComponents(1);
491 ascendingManifold->SetNumberOfTuples(numberOfVertices);
494 descendingManifold->SetNumberOfComponents(1);
495 descendingManifold->SetNumberOfTuples(numberOfVertices);
498 morseSmaleManifold->SetNumberOfComponents(1);
499 morseSmaleManifold->SetNumberOfTuples(numberOfVertices);
510 const auto imageDataInput = vtkImageData::SafeDownCast(input);
512 if(DiscreteGradientBackend == 0) {
514 DiscreteGradient::BACKEND::CLASSIC_BACKEND);
516 if(DiscreteGradientBackend == 1 && !imageDataInput) {
518 DiscreteGradient::BACKEND::CLASSIC_BACKEND);
519 this->
printWrn(
"The stochastic gradient (IEEE TVCG 2012) can only");
520 this->
printWrn(
"be used on vtkImageData (.vti).");
521 this->
printWrn(
"Defaulting to homotopic expansion (IEEE PAMI 2011)");
523 if(DiscreteGradientBackend == 1 && imageDataInput) {
525 DiscreteGradient::BACKEND::STOCHASTIC_BACKEND);
543 inputScalars->GetDataType(), triangulation->getType(),
545 inputScalars, outputCriticalPoints, outputSeparatrices1,
547 *
static_cast<TTK_TT *
>(triangulation->getData()))));
553 outputMorseComplexes->ShallowCopy(input);
556 vtkPointData *pointData = outputMorseComplexes->GetPointData();
557#ifndef TTK_ENABLE_KAMIKAZE
559 this->
printErr(
"outputMorseComplexes has no point data.");
564 if(ComputeDescendingSegmentation)
565 pointData->AddArray(descendingManifold);
567 pointData->AddArray(ascendingManifold);
570 pointData->AddArray(morseSmaleManifold);