99 vtkInformationVector **inputVector,
100 vtkInformationVector *outputVector) {
104 auto blockBary = vtkMultiBlockDataSet::GetData(inputVector[0], 0);
105 auto tableCoefficients = vtkTable::GetData(inputVector[1], 0);
106 auto tableVectors = vtkTable::GetData(inputVector[2], 0);
107 auto tableCorrelation = vtkTable::GetData(inputVector[3], 0);
108 auto blockInputTrees = vtkMultiBlockDataSet::GetData(inputVector[4], 0);
113 std::vector<vtkSmartPointer<vtkMultiBlockDataSet>> inputBary, inputTrees;
118 if((baryTreeNodes.size() != 0
119 and inputBary[0]->GetBlock(0) != baryTreeNodes[0])
121 and
int(
allTs_[0].size()) != tableCoefficients->GetNumberOfRows())
122 or tableCoefficientsMTime != tableCoefficients->GetMTime()
123 or tableVectorsMTime != tableVectors->GetMTime()
125 and tableCorrelationMTime != tableCorrelation->GetMTime())
127 and blockInputTreesMTime != blockInputTrees->GetMTime())) {
128 resetDataVisualization();
130 tableCoefficientsMTime = tableCoefficients->GetMTime();
131 tableVectorsMTime = tableVectors->GetMTime();
133 tableCorrelationMTime = tableCorrelation->GetMTime();
135 blockInputTreesMTime = blockInputTrees->GetMTime();
138 printMsg(
"Load parameters from field data.");
139 std::vector<std::string> paramNames;
141 for(
auto paramName : paramNames) {
142 auto array = tableCoefficients->GetFieldData()->GetArray(paramName.c_str());
144 double const value = array->GetTuple1(0);
147 printMsg(
" - " + paramName +
" was not found in the field data.");
152 printMsg(
"Computation with normalized Wasserstein.");
154 printMsg(
"Computation without normalized Wasserstein.");
159 auto tableVNoCols = tableVectors->GetNumberOfColumns() / inputBary.size();
160 unsigned int const numberOfGeodesics = tableVNoCols / 4;
161 auto numberOfInputs = tableCoefficients->GetNumberOfRows();
163 auto baryNoNodes = tableVectors->GetNumberOfRows();
164 auto baryNoNodes2 = baryNoNodes;
165 if(inputBary.size() > 1) {
166 for(
unsigned int i = 0; i < 2; ++i) {
167 auto &baryNoNodesT = (i == 0 ? baryNoNodes : baryNoNodes2);
172 ->GetVariantValue(baryNoNodesT - 1)
182 std::vector<std::string> nonFieldDataNames{
"TreeID"};
183 allTs_.resize(numberOfGeodesics, std::vector<double>(numberOfInputs, 0.0));
184 for(
unsigned int i = 0; i < numberOfGeodesics; ++i)
185 for(
unsigned int j = 0; j < numberOfInputs; ++j) {
187 allTs_[i][j] = tableCoefficients->GetColumnByName(name.c_str())
190 nonFieldDataNames.push_back(name);
191 std::string
const normName
193 nonFieldDataNames.push_back(normName);
198 vtkNew<vtkFieldData> fd{};
199 for(
int i = 0; i < tableCoefficients->GetNumberOfColumns(); ++i) {
200 std::string
const name{tableCoefficients->GetColumnName(i)};
201 if(std::count(nonFieldDataNames.begin(), nonFieldDataNames.end(), name)
203 auto array = tableCoefficients->GetColumn(i);
204 array->SetName(name.c_str());
208 inputFieldData = vtkFieldData::New();
209 inputFieldData->ShallowCopy(fd);
214 pVS_.resize(numberOfGeodesics, std::vector<double *>(2));
216 if(inputBary.size() > 1) {
217 pTrees2Vs_.resize(numberOfGeodesics, std::vector<double *>(2));
220 for(
unsigned int h = 0; h < inputBary.size(); ++h) {
223 bool const secondInput = (h == 1);
224 for(
unsigned int i = 0; i < numberOfGeodesics; ++i) {
225 for(
unsigned int k = 0; k < 2; ++k) {
226 std::string
const name1
228 std::string
const name2
230 pVS[i][k] = ttkUtils::GetPointer<double>(vtkDataArray::SafeDownCast(
231 tableVectors->GetColumnByName(name1.c_str())));
232 pV2s[i][k] = ttkUtils::GetPointer<double>(vtkDataArray::SafeDownCast(
233 tableVectors->GetColumnByName(name2.c_str())));
243 if(tableCorrelation) {
247 auto tableCorrNoRows = tableCorrelation->GetNumberOfRows();
248 auto baryNodeIdArray = tableCorrelation->GetColumnByName(
"BaryNodeId");
250 for(
unsigned int j = 0; j < numberOfGeodesics; ++j) {
252 std::string
const name2
255 = ttkUtils::GetPointer<double>(vtkDataArray::SafeDownCast(
256 tableCorrelation->GetColumnByName(name.c_str())));
258 = ttkUtils::GetPointer<double>(vtkDataArray::SafeDownCast(
259 tableCorrelation->GetColumnByName(name2.c_str())));
265 for(
unsigned int i = 0; i < numberOfInputs; ++i) {
267 auto array = tableCorrelation->GetColumnByName(name.c_str());
270 for(
unsigned int j = 0; j < tableCorrNoRows; ++j) {
271 auto baryNodeId = baryNodeIdArray->GetVariantValue(j).ToInt();
273 baryNodeId, array->GetVariantValue(j).ToDouble(), 0.0);
285 return run<float>(outputVector, inputBary, inputTrees);
301 vtkInformationVector *
ttkNotUsed(outputVector),
309 const int numTrees = inputTrees.size();
311 setDataVisualization(numTrees);
313 std::vector<ttk::ftm::MergeTree<dataType>> baryDTree, inputDTrees;
315 std::vector<bool> useSadMaxPairsVec{
false,
true};
317 useSadMaxPairsVec.erase(useSadMaxPairsVec.begin());
318 ttk::ftm::constructTrees<dataType>(inputBary, baryDTree, baryTreeNodes,
319 baryTreeArcs, baryTreeSegmentation,
323 or (ReconstructInputTrees
325 bool const useSadMaxPairs
327 bool const isInputPD = ttk::ftm::constructTrees<dataType>(
328 inputTrees, inputDTrees, inputTreesNodes, inputTreesArcs,
329 inputTreesSegmentation, useSadMaxPairs);
337 std::string
const preFix{(processFirstInput ?
"" :
"Second Input ")};
338 auto &baryToUse = (processFirstInput ? baryDTree[0] : baryDTree[1]);
343 printMsg(
"Compute Geodesics Distance...");
344 computeGeodesicsDistance<dataType>(baryDTree);
348 if(OutputInputTrees) {
351 processInputTrees<dataType>(inputDTrees);
355 if(ReconstructInputTrees) {
357 printMsg(
"Reconstruct " + preFix +
"Input Trees...");
358 std::vector<ttk::ftm::MergeTree<dataType>> reconstructedTreesT;
359 reconstruction<dataType>(baryToUse, inputDTrees, reconstructedTreesT,
360 reconstructionErrors, recInputMatchings,
361 recBaryMatchings, !processFirstInput);
362 ttk::ftm::mergeTreesTemplateToDouble<dataType>(
363 reconstructedTreesT, reconstructedTrees);
365 reconstructedTrees.clear();
368 if(ConstructGeodesicsTrees) {
370 printMsg(
"Construct Geodesics Trees...");
371 std::vector<std::vector<ttk::ftm::MergeTree<dataType>>> geodesicsTreesT;
372 constructGeodesicsTrees<dataType>(
373 baryToUse, geodesicsTreesT, !processFirstInput);
374 ttk::ftm::mergeTreesTemplateToDouble<dataType>(
375 geodesicsTreesT, geodesicsTrees);
377 geodesicsTrees.clear();
380 if(ConstructEllipses) {
382 printMsg(
"Construct Ellipses Trees...");
383 std::vector<ttk::ftm::MergeTree<dataType>> geodesicsEllipsesT;
384 constructGeodesicsEllipses<dataType>(
385 baryToUse, geodesicsEllipsesT, !processFirstInput);
386 ttk::ftm::mergeTreesTemplateToDouble<dataType>(
387 geodesicsEllipsesT, geodesicsEllipses);
389 geodesicsEllipses.clear();
392 if(ConstructRectangle) {
394 printMsg(
"Construct Rectangle Trees...");
395 std::vector<ttk::ftm::MergeTree<dataType>> geodesicsRectangleT;
396 constructGeodesicsRectangle<dataType>(
397 baryToUse, geodesicsRectangleT, RectangleMultiplier, !processFirstInput);
398 ttk::ftm::mergeTreesTemplateToDouble<dataType>(
399 geodesicsRectangleT, geodesicsRectangle);
401 geodesicsRectangle.clear();
404 if(ConstructSurface) {
406 printMsg(
"Construct " + preFix +
"Surface Trees...");
407 std::vector<ttk::ftm::MergeTree<dataType>> geodesicsSurfaceT;
408 constructGeodesicsSurface<dataType>(
409 baryToUse, geodesicsSurfaceT, !processFirstInput);
410 ttk::ftm::mergeTreesTemplateToDouble<dataType>(
411 geodesicsSurfaceT, geodesicsSurface);
413 geodesicsSurface.clear();
415 ttk::ftm::mergeTreesTemplateToDouble<dataType>(baryDTree, baryMTree);
416 ttk::ftm::mergeTreesTemplateToDouble<dataType>(inputDTrees, inputMTrees);
423 vtkInformationVector *outputVector,
431 auto output = vtkMultiBlockDataSet::GetData(outputVector, 0);
433 unsigned int const noInput = (OutputInputTrees ? inputTrees.size() : 0);
434 unsigned int const noReconst = reconstructedTrees.size();
435 unsigned int const noBary = (OutputBarycenter ? 1 : 0);
436 unsigned int const noGeod
437 = geodesicsTrees.size()
438 * (geodesicsTrees.size() == 0 ? 0 : geodesicsTrees[0].size());
439 unsigned int const noEllipses = geodesicsEllipses.size();
440 unsigned int const noRectangle = geodesicsRectangle.size();
441 unsigned int const noSurface = geodesicsSurface.size();
451 std::vector<unsigned int> allNo{
452 noInput, noReconst, noBary, noGeod, noEllipses, noRectangle, noSurface};
453 std::vector<unsigned int> steps(allNo.size());
454 for(
unsigned int i = 0; i < allNo.size(); ++i)
455 steps[i] = (i == 0 ? 0 : steps[i - 1]) + allNo[i];
457 unsigned int const noTreesTotal = steps[steps.size() - 1];
459 unsigned int const noBlocks
461 output->SetNumberOfBlocks(noBlocks);
464 vtkBlockNodes->SetNumberOfBlocks(noTreesTotal);
465 output->SetBlock(0, vtkBlockNodes);
469 vtkBlockArcs->SetNumberOfBlocks(noTreesTotal);
470 output->SetBlock(1, vtkBlockArcs);
472 if(OutputInputTreesSegmentation) {
475 vtkBlockSegs->SetNumberOfBlocks(noTreesTotal);
479 int const geodesicsTreesOffset = std::max((
int)geodesicsTrees.size() - 1, 0);
484 std::vector<std::vector<ttk::ftm::idNode>> matchingMatrix;
486 getMatchingMatrix<double>(
491 printWrn(
"Please provide input trees and correlation matrix to transfer "
492 "input trees information.");
494#ifdef TTK_ENABLE_OPENMP
495#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
497 for(
unsigned int i = 0; i < noTreesTotal; ++i) {
501 std::vector<double> ts;
505 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(inputMTrees[index], mt);
508 }
else if(i < steps[1]) {
510 index = i - steps[0];
511 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(
512 reconstructedTrees[index], mt);
515 }
else if(i < steps[2]) {
517 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(baryMTree[0], mt);
519 std::array<double, 2> middle;
524 }
else if(i < steps[3]) {
526 index = i - steps[2];
527 int const i0 = index / geodesicsTrees[0].size();
528 int const i1 = index % geodesicsTrees[0].size();
529 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(geodesicsTrees[i0][i1], mt);
532 }
else if(i < steps[4]) {
534 index = i - steps[3];
535 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(
536 geodesicsEllipses[index], mt);
537 treeType = 3 + geodesicsTreesOffset;
539 }
else if(i < steps[5]) {
541 index = i - steps[4];
542 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(
543 geodesicsRectangle[index], mt);
544 treeType = 4 + geodesicsTreesOffset;
546 }
else if(i < steps[6]) {
548 index = i - steps[5];
549 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(
550 geodesicsSurface[index], mt);
551 treeType = 5 + geodesicsTreesOffset;
555 bool const isInputTree = (i < steps[0]);
556 bool const isInputTreeAndGotMesh
557 = (isInputTree and inputTrees[i]->GetNumberOfBlocks() >= 3);
558 bool const isReconstructedTree = (i >= steps[0] and i < steps[1]);
559 bool const isInputOrReconstructedTree
560 = (isInputTree or isReconstructedTree);
575 for(
unsigned int n = 0; n <
vSize_; ++n) {
576 if((
int)matchingMatrix[n][i] != -1) {
591 for(
unsigned int n = 0; n <
vSize_; ++n) {
592 if((
int)matchingMatrix[n][i] != -1) {
593 baryNodeID[matchingMatrix[n][i]] = n;
594 baryNodePers[matchingMatrix[n][i]]
595 = baryMTree[0].tree.template getNodePersistence<double>(n);
598 std::string nameBaryNodePers{
"BaryNodePers"};
600 std::string nameBaryNodeId{
"BaryNodeId"};
611 if(isInputTreeAndGotMesh and OutputInputTreesSegmentation) {
622 if(isReconstructedTree) {
625 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(baryMTree[0], baryMT);
626 std::vector<ttk::ftm::idNode> matchingVector;
628 mt, baryMT, recBaryMatchings[index], matchingVector);
630 for(
unsigned int n = 0; n <
vSize_; ++n) {
631 auto match = matchingVector[n];
633 baryNodeID[match] = n;
635 std::string nameBaryNodeId{
"BaryNodeId"};
641 ttk::ftm::mergeTreeDoubleToTemplate<dataType>(
642 inputMTrees[index], inputMT);
643 std::vector<ttk::ftm::idNode> matchingVector;
645 mt, inputMT, recInputMatchings[index], matchingVector);
647 for(
unsigned int n = 0; n <
vSize_; ++n) {
648 auto match = matchingMatrix[n][index];
649 if((
int)match != -1) {
650 match = matchingVector[match];
652 baryNodeID[match] = n;
655 std::string nameBaryNodeId{
"BaryNodeId"};
675 if(noReconst != 0 or noInput != 0) {
676 vtkNew<vtkFieldData> fd{};
677 if(isInputOrReconstructedTree)
678 fd->DeepCopy(inputFieldData);
680 fd->CopyStructure(inputFieldData);
681 for(
int a = 0; a < inputFieldData->GetNumberOfArrays(); ++a) {
682 auto array = fd->GetAbstractArray(a);
683 array->SetNumberOfTuples(1);
684 if(isInputOrReconstructedTree) {
685 vtkNew<vtkIdList> id;
686 id->InsertNextId(i % (noInput == 0 ? noReconst : noInput));
687 inputFieldData->GetAbstractArray(a)->GetTuples(
id, array);
689 auto downArray = vtkDataArray::SafeDownCast(array);
691 downArray->SetTuple1(0, std::nan(
""));
694 vtkOutputNode->GetFieldData()->ShallowCopy(fd);
698 vtkNew<vtkIntArray> treeTypeArray{};
699 treeTypeArray->SetName(
"TreeType");
700 treeTypeArray->InsertNextTuple1(treeType);
701 vtkOutputNode->GetFieldData()->AddArray(treeTypeArray);
703 vtkNew<vtkIntArray> treeIDArray{};
704 treeIDArray->SetName(
"TreeID");
705 treeIDArray->InsertNextTuple1(i);
706 vtkOutputNode->GetFieldData()->AddArray(treeIDArray);
708 for(
unsigned int j = 0; j < ts.size(); ++j) {
709 vtkNew<vtkDoubleArray> tArray{};
711 tArray->SetName(name.c_str());
712 tArray->InsertNextTuple1(ts[j]);
713 vtkOutputNode->GetFieldData()->AddArray(tArray);
715 vtkNew<vtkDoubleArray> scaledTArray{};
717 scaledTArray->SetName(name2.c_str());
719 vtkOutputNode->GetFieldData()->AddArray(scaledTArray);
722 bool const isSurfaceTree = (treeType == 5 + geodesicsTreesOffset);
723 vtkNew<vtkIntArray> isSurfaceArray{};
724 isSurfaceArray->SetName(
"isSurface");
725 isSurfaceArray->InsertNextTuple1(isSurfaceTree);
726 vtkOutputNode->GetFieldData()->AddArray(isSurfaceArray);
728 vtkNew<vtkIntArray> isBoundaryArray{};
729 isBoundaryArray->SetName(
"isBoundary");
730 isBoundaryArray->InsertNextTuple1(isBoundary);
731 vtkOutputNode->GetFieldData()->AddArray(isBoundaryArray);
733 vtkNew<vtkIntArray> BoundaryIDArray{};
734 BoundaryIDArray->SetName(
"BoundaryID");
735 BoundaryIDArray->InsertNextTuple1(boundaryID);
736 vtkOutputNode->GetFieldData()->AddArray(BoundaryIDArray);
738 vtkNew<vtkDoubleArray> reconstructionErrorArray{};
739 reconstructionErrorArray->SetName(
"ReconstructionError");
740 auto err = (isReconstructedTree and reconstructionErrors.size() != 0
741 ? reconstructionErrors[index]
743 reconstructionErrorArray->InsertNextTuple1(err);
744 vtkOutputNode->GetFieldData()->AddArray(reconstructionErrorArray);
747 vtkMultiBlockDataSet::SafeDownCast(output->GetBlock(0))
748 ->SetBlock(i, vtkOutputNode);
750 vtkMultiBlockDataSet::SafeDownCast(output->GetBlock(1))
751 ->SetBlock(i, vtkOutputArc);
752 if(OutputInputTreesSegmentation)
753 vtkMultiBlockDataSet::SafeDownCast(
755 ->SetBlock(i, vtkOutputSegmentation);
759 std::vector<std::string> paramNames;
761 for(
auto paramName : paramNames) {
762 vtkNew<vtkDoubleArray> array{};
763 array->SetName(paramName.c_str());
765 output->GetFieldData()->AddArray(array);