From 64477b3800e587a1b65dfa7fdbaa530c122951f9 Mon Sep 17 00:00:00 2001 From: Mingzhe Li Date: Mon, 8 Jul 2024 17:00:42 -0700 Subject: [PATCH] Moving implementation of SelectTopVolumeContours into SelectTopVolumeContoursBlock.cxx --- .../SelectTopVolumeContoursFilter.cxx | 2070 ++++++++--------- .../SelectTopVolumeContoursFilter.h | 4 + .../internal/SelectTopVolumeContoursBlock.cxx | 1095 ++++++++- .../internal/SelectTopVolumeContoursBlock.h | 12 +- 4 files changed, 2129 insertions(+), 1052 deletions(-) diff --git a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx index 69f956ce1..b9de3166e 100644 --- a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx +++ b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx @@ -16,14 +16,6 @@ #include #include #include -#include -#include -#include -#include -#include -#include - - // vtkm includes #include @@ -130,78 +122,80 @@ VTKM_CONT vtkm::cont::PartitionedDataSet SelectTopVolumeContoursFilter::DoExecut << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); - branch_top_volume_master.foreach ( - [&](SelectTopVolumeContoursBlock* b, const vtkmdiy::Master::ProxyWithLink&) { - using vtkm::worklet::contourtree_augmented::IdArrayType; - const auto& globalSize = firstGlobalPointDimensions; - vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; - const vtkm::cont::DataSet& ds = input.GetPartition(b->LocalBlockNo); + branch_top_volume_master.foreach ([&](SelectTopVolumeContoursBlock* b, + const vtkmdiy::Master::ProxyWithLink&) { + using vtkm::worklet::contourtree_augmented::IdArrayType; + const auto& globalSize = firstGlobalPointDimensions; + vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; + const vtkm::cont::DataSet& ds = input.GetPartition(b->LocalBlockNo); - b->SortBranchByVolume(ds, totalVolume); + b->SortBranchByVolume(ds, totalVolume); + b->SelectLocalTopVolumeBranch(ds, this->GetSavedBranches()); - // copy the top volume branches into a smaller array - // we skip index 0 because it must be the main branch (which has the highest volume) - vtkm::Id nActualSavedBranches = - std::min(this->nSavedBranches, b->SortedBranchByVolume.GetNumberOfValues() - 1); + //// copy the top volume branches into a smaller array + //// we skip index 0 because it must be the main branch (which has the highest volume) + //vtkm::Id nActualSavedBranches = + // std::min(this->nSavedBranches, b->SortedBranchByVolume.GetNumberOfValues() - 1); - vtkm::worklet::contourtree_augmented::IdArrayType topVolumeBranch; - vtkm::cont::Algorithm::CopySubRange( - b->SortedBranchByVolume, 1, nActualSavedBranches, topVolumeBranch); + //vtkm::worklet::contourtree_augmented::IdArrayType topVolumeBranch; + //vtkm::cont::Algorithm::CopySubRange( + // b->SortedBranchByVolume, 1, nActualSavedBranches, topVolumeBranch); - auto branchRootByBranch = ds.GetField("BranchRootByBranch") - .GetData() - .AsArrayHandle>(); + //auto branchRootByBranch = ds.GetField("BranchRootByBranch") + // .GetData() + // .AsArrayHandle>(); - const vtkm::Id nBranches = branchRootByBranch.GetNumberOfValues(); + //const vtkm::Id nBranches = branchRootByBranch.GetNumberOfValues(); - auto branchRootGRId = - ds.GetField("BranchRootGRId").GetData().AsArrayHandle>(); + //auto branchRootGRId = + // ds.GetField("BranchRootGRId").GetData().AsArrayHandle>(); - auto upperEndGRId = ds.GetField("UpperEndGlobalRegularIds") - .GetData() - .AsArrayHandle>(); + //auto upperEndGRId = ds.GetField("UpperEndGlobalRegularIds") + // .GetData() + // .AsArrayHandle>(); - auto lowerEndGRId = ds.GetField("LowerEndGlobalRegularIds") - .GetData() - .AsArrayHandle>(); + //auto lowerEndGRId = ds.GetField("LowerEndGlobalRegularIds") + // .GetData() + // .AsArrayHandle>(); - vtkm::cont::Algorithm::Copy(branchRootByBranch, b->BranchRootByBranch); - vtkm::cont::Algorithm::Copy(branchRootGRId, b->BranchRootGRId); - b->IsParentBranch.AllocateAndFill(nBranches, false); + //vtkm::cont::Algorithm::Copy(branchRootByBranch, b->BranchRootByBranch); + //vtkm::cont::Algorithm::Copy(branchRootGRId, b->BranchRootGRId); + //b->IsParentBranch.AllocateAndFill(nBranches, false); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - branchRootGRId, topVolumeBranch, b->TopVolumeBranchRootGRId); + //vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // branchRootGRId, topVolumeBranch, b->TopVolumeBranchRootGRId); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - upperEndGRId, topVolumeBranch, b->TopVolumeBranchUpperEndGRId); + //vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // upperEndGRId, topVolumeBranch, b->TopVolumeBranchUpperEndGRId); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - lowerEndGRId, topVolumeBranch, b->TopVolumeBranchLowerEndGRId); + //vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // lowerEndGRId, topVolumeBranch, b->TopVolumeBranchLowerEndGRId); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->BranchVolume, topVolumeBranch, b->TopVolumeBranchVolume); + //vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->BranchVolume, topVolumeBranch, b->TopVolumeBranchVolume); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->BranchSaddleEpsilon, topVolumeBranch, b->TopVolumeBranchSaddleEpsilon); + //vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->BranchSaddleEpsilon, topVolumeBranch, b->TopVolumeBranchSaddleEpsilon); - auto resolveArray = [&](const auto& inArray) { - using InArrayHandleType = std::decay_t; - InArrayHandleType topVolBranchSaddleIsoValue; - vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( - inArray, topVolumeBranch, topVolBranchSaddleIsoValue); - b->TopVolumeBranchSaddleIsoValue = topVolBranchSaddleIsoValue; - }; + //auto resolveArray = [&](const auto& inArray) { + // using InArrayHandleType = std::decay_t; + // InArrayHandleType topVolBranchSaddleIsoValue; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + // inArray, topVolumeBranch, topVolBranchSaddleIsoValue); + // b->TopVolumeBranchSaddleIsoValue = topVolBranchSaddleIsoValue; + //}; - b->BranchSaddleIsoValue - .CastAndCallForTypes(resolveArray); - }); + //b->BranchSaddleIsoValue + // .CastAndCallForTypes(resolveArray); + }); timingsStream << " " << std::setw(60) << std::left << "SelectBranchByVolume" << ": " << timer.GetElapsedTime() << " seconds" << std::endl; timer.Start(); // We apply all-to-all broadcast to collect the top nSavedBranches branches by volume + // TODO: change all_to_all to swap with partner? It does not seem necessary to use all_to_all vtkmdiy::all_to_all(branch_top_volume_master, assigner, vtkm::filter::scalar_topology::internal::SelectTopVolumeContoursFunctor( @@ -214,396 +208,395 @@ VTKM_CONT vtkm::cont::PartitionedDataSet SelectTopVolumeContoursFilter::DoExecut // we compute the hierarchy of selected branches adding the root branch for each block branch_top_volume_master.foreach ([&](SelectTopVolumeContoursBlock* b, const vtkmdiy::Master::ProxyWithLink&) { - using vtkm::worklet::contourtree_augmented::IdArrayType; const vtkm::cont::DataSet& ds = input.GetPartition(b->LocalBlockNo); - vtkm::cont::Invoker invoke; - - // we need upper/lower local ends and global ends for hierarchy of branches - auto upperLocalEndIds = - ds.GetField("UpperEndLocalIds").GetData().AsArrayHandle>(); - auto lowerLocalEndIds = - ds.GetField("LowerEndLocalIds").GetData().AsArrayHandle>(); - auto globalRegularIds = ds.GetField("RegularNodeGlobalIds") - .GetData() - .AsArrayHandle>(); - IdArrayType upperEndGRIds = - ds.GetField("UpperEndGlobalRegularIds").GetData().AsArrayHandle(); - IdArrayType lowerEndGRIds = - ds.GetField("LowerEndGlobalRegularIds").GetData().AsArrayHandle(); - - // let's check which top volume branches are known by the block - // We check the branchGRId of top volume branches to see whether there are matches within the block - const vtkm::Id nTopVolBranches = b->TopVolumeBranchLowerEndGRId.GetNumberOfValues(); - - // sortedBranchOrder: the branch order (in the ascending order of branch root) - // the high-level idea is to sort the branch root global regular ids - // and for each top-volume branch, we use binary search to get the original branch index - // if the top-volume branch does not exist in the block, it will be dropped out - IdArrayType sortedBranchGRId; - IdArrayType sortedBranchOrder; - vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(b->BranchRootGRId.GetNumberOfValues()), - sortedBranchOrder); - vtkm::cont::Algorithm::Copy(b->BranchRootGRId, sortedBranchGRId); - vtkm::cont::Algorithm::SortByKey(sortedBranchGRId, sortedBranchOrder); - - b->TopVolBranchKnownByBlockStencil.Allocate(nTopVolBranches); - b->TopVolBranchGROrder.Allocate(nTopVolBranches); - - // We reuse the IdxIfWithinBlockWorklet. - // This worklet searches for given values in a sorted array and returns the stencil & index if the value exists in the array. - // b->TopVolBranchGROrder: the order of the topVolBranch among all known branches if the branch is known by the block. - auto idxIfBranchWithinBlockWorklet = - vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); - invoke(idxIfBranchWithinBlockWorklet, - b->TopVolumeBranchRootGRId, - sortedBranchGRId, - b->TopVolBranchKnownByBlockStencil, - b->TopVolBranchGROrder); - - // Dropping out top-volume branches that are not known by the block. - - // the index of top-volume branches known by the block among all top-volume branches - IdArrayType topVolBranchKnownByBlockIndex; - vtkm::cont::ArrayHandleIndex topVolBranchesIndex(nTopVolBranches); - vtkm::cont::Algorithm::CopyIf( - topVolBranchesIndex, b->TopVolBranchKnownByBlockStencil, topVolBranchKnownByBlockIndex); - - const vtkm::Id nTopVolBranchKnownByBlock = topVolBranchKnownByBlockIndex.GetNumberOfValues(); - - // filtered b->TopVolBranchGROrder, by removing NO_SUCH_ELEMENT - IdArrayType topVolBranchFilteredGROrder; - - // b->TopVolBranchInfoActualIndex: the information index of the top-volume branch - vtkm::cont::Algorithm::CopyIf( - b->TopVolBranchGROrder, b->TopVolBranchKnownByBlockStencil, topVolBranchFilteredGROrder); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - sortedBranchOrder, topVolBranchFilteredGROrder, b->TopVolBranchInfoActualIndex); - - // filtered branch saddle epsilons, global lower/upper end GR ids, - IdArrayType topVolFilteredBranchSaddleEpsilon; - IdArrayType topVolFilteredLowerEndGRId; - IdArrayType topVolFilteredUpperEndGRId; - vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchSaddleEpsilon, - b->TopVolBranchKnownByBlockStencil, - topVolFilteredBranchSaddleEpsilon); - vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchUpperEndGRId, - b->TopVolBranchKnownByBlockStencil, - topVolFilteredUpperEndGRId); - vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchLowerEndGRId, - b->TopVolBranchKnownByBlockStencil, - topVolFilteredLowerEndGRId); - - // for each top-vol branch known by the block - // we get their upper end and lower end local ids - IdArrayType topVolBranchUpperLocalEnd; - IdArrayType topVolBranchLowerLocalEnd; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - upperLocalEndIds, b->TopVolBranchInfoActualIndex, topVolBranchUpperLocalEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - lowerLocalEndIds, b->TopVolBranchInfoActualIndex, topVolBranchLowerLocalEnd); - - IdArrayType topVolLowerLocalEndGRId; - IdArrayType topVolUpperLocalEndGRId; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - globalRegularIds, topVolBranchLowerLocalEnd, topVolLowerLocalEndGRId); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - globalRegularIds, topVolBranchUpperLocalEnd, topVolUpperLocalEndGRId); - - // Below is the code to compute the branch hierarchy of top-volume branches - // We need this information because we not only want to visualize the contour - // on top-volume branches, but also on their parent branches. - // Because we use volume as the metric, the parent branch of a top-volume branch - // is either a top-volume branch or the root branch (where both ends are leaf nodes) - vtkm::worklet::scalar_topology::select_top_volume_contours::BranchSaddleIsKnownWorklet - branchSaddleIsKnownWorklet; - // the branch saddle local ID if the saddle end is known by the block - IdArrayType branchSaddleIsKnown; - branchSaddleIsKnown.Allocate(nTopVolBranchKnownByBlock); - - invoke(branchSaddleIsKnownWorklet, // worklet - topVolFilteredLowerEndGRId, // input - topVolBranchLowerLocalEnd, // input - topVolLowerLocalEndGRId, // input - topVolFilteredUpperEndGRId, // input - topVolBranchUpperLocalEnd, // input - topVolUpperLocalEndGRId, // input - topVolFilteredBranchSaddleEpsilon, // input - branchSaddleIsKnown); // output - - // the order of top volume branches with parents known by the block - IdArrayType topVolChildBranch; - IdArrayType topVolChildBranchSaddle; - - vtkm::cont::Algorithm::CopyIf( - topVolBranchKnownByBlockIndex, - branchSaddleIsKnown, - topVolChildBranch, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); - vtkm::cont::Algorithm::CopyIf( - branchSaddleIsKnown, - branchSaddleIsKnown, - topVolChildBranchSaddle, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); - - const vtkm::Id nChildBranch = topVolChildBranch.GetNumberOfValues(); - - // to compute the parent branch, we need to - // 1. for the branch saddle end, collect all superarcs involving it - // 2. get the branch information for selected superarcs - // 3. eliminate branch information for branches sharing the same saddle end - - auto superarcs = - ds.GetField("Superarcs").GetData().AsArrayHandle>(); - auto branchRoots = - ds.GetField("BranchRoots").GetData().AsArrayHandle>(); - VTKM_ASSERT(superarcs.GetNumberOfValues() == branchRoots.GetNumberOfValues()); - - // we sort all superarcs by target to allow binary search - IdArrayType superarcsByTarget; - vtkm::worklet::scalar_topology::select_top_volume_contours::SuperarcTargetComparator - superarcComparator(superarcs); - vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(superarcs.GetNumberOfValues()), - superarcsByTarget); - vtkm::cont::Algorithm::Sort(superarcsByTarget, superarcComparator); - - IdArrayType permutedSuperarcs; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - superarcs, superarcsByTarget, permutedSuperarcs); - - IdArrayType permutedBranchRoots; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - branchRoots, superarcsByTarget, permutedBranchRoots); - - // the branch root of the superarc of the branch saddle supernode - IdArrayType topVolChildBranchSaddleBranchRoot; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - branchRoots, topVolChildBranchSaddle, topVolChildBranchSaddleBranchRoot); - - // the GR Ids of the superarc of the branch saddle supernode - IdArrayType topVolChildBranchSaddleGRIds; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - globalRegularIds, topVolChildBranchSaddle, topVolChildBranchSaddleGRIds); - - // there is a debate to find all superarcs connect to a supernode - // strategy 1. iterate through saddles and parallelize over superarcs for search - // time complexity: O(nTopVolBranches) - // (nTopVolBranches usually <= 100, based on input parameter setting) + b->ComputeTopVolumeBranchHierarchy(ds); + // using vtkm::worklet::contourtree_augmented::IdArrayType; + // vtkm::cont::Invoker invoke; // - // strategy 2. parallelize over all saddles and use binary search to find superarcs - // time complexity: O(log_2(nSuperarcs)) (nSuperarcs can be considerably large) + // // we need upper/lower local ends and global ends for hierarchy of branches + // auto upperLocalEndIds = + // ds.GetField("UpperEndLocalIds").GetData().AsArrayHandle>(); + // auto lowerLocalEndIds = + // ds.GetField("LowerEndLocalIds").GetData().AsArrayHandle>(); + // auto globalRegularIds = ds.GetField("RegularNodeGlobalIds") + // .GetData() + // .AsArrayHandle>(); + // IdArrayType upperEndGRIds = + // ds.GetField("UpperEndGlobalRegularIds").GetData().AsArrayHandle(); + // IdArrayType lowerEndGRIds = + // ds.GetField("LowerEndGlobalRegularIds").GetData().AsArrayHandle(); // - // here, we choose strategy 2 for better extendability to high nTopVolBranches - // but in tests using nTopVolBranches <= 10, strategy 1 is generally faster - - // note: after getting the branch root superarc, we use binary search to get the branch order - // because BranchRootByBranch is sorted by branch root (superarc) id - -#ifdef DEBUG_PRINT - std::stringstream parentBranchStream; - - parentBranchStream << "Debug for Parent Branch, Block " << b->LocalBlockNo << std::endl; - vtkm::worklet::contourtree_augmented::PrintHeader(nChildBranch, parentBranchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Child Branch Saddle", topVolChildBranchSaddle, -1, parentBranchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Child Saddle Root", topVolChildBranchSaddleBranchRoot, -1, parentBranchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Child Saddle GR Id", topVolChildBranchSaddleGRIds, -1, parentBranchStream); - - vtkm::worklet::contourtree_augmented::PrintHeader(superarcs.GetNumberOfValues(), - parentBranchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Permuted Superarcs", permutedSuperarcs, -1, parentBranchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Permuted Branch roots", permutedBranchRoots, -1, parentBranchStream); - - VTKM_LOG_S(vtkm::cont::LogLevel::Info, parentBranchStream.str()); -#endif // DEBUG_PRINT - - // the corresponding parent branch of child branches - IdArrayType topVolChildBranchParent; - topVolChildBranchParent.AllocateAndFill(nChildBranch, - vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT); - vtkm::worklet::scalar_topology::select_top_volume_contours::GetParentBranchWorklet - getParentBranchWorklet; - invoke(getParentBranchWorklet, - topVolChildBranchSaddle, - topVolChildBranchSaddleBranchRoot, - topVolChildBranchSaddleGRIds, - permutedSuperarcs, - permutedBranchRoots, - b->BranchRootByBranch, - upperEndGRIds, - lowerEndGRIds, - topVolChildBranchParent); - - b->TopVolumeBranchParent.AllocateAndFill( - nTopVolBranches, vtkm::Id(vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)); - - vtkm::worklet::scalar_topology::select_top_volume_contours::AssignValueByIndex - assignParentBranch; - // for each top volume branch, assign the parent branch info id in the block - invoke( - assignParentBranch, topVolChildBranch, topVolChildBranchParent, b->TopVolumeBranchParent); - // for each branch, assign true if it is a parent branch - invoke(assignParentBranch, - topVolChildBranchParent, - vtkm::cont::ArrayHandleConstant(true, nChildBranch), - b->IsParentBranch); - - // Update 06/17/2024: temporarily moving the initialization of IsOuterSaddle outside the function below - // due to a compile error on ubuntu1604+gcc5 - IdArrayType IsOuterSaddle; - IsOuterSaddle.Allocate(nTopVolBranches); - - // sort all top-volume branches based on - // 1. parent branch info id: b->TopVolumeBranchParent - // 2. saddle-end value: b->TopVolumeBranchSaddleIsovalue - // 3. branch root global regular id (anything that can break tie) - auto resolveBranchParent = [&](const auto& inArray) { - using InArrayHandleType = std::decay_t; - using ValueType = typename InArrayHandleType::ValueType; - - vtkm::worklet::scalar_topology::select_top_volume_contours::BranchParentComparator - parentComparator(b->TopVolumeBranchParent, inArray, b->TopVolumeBranchRootGRId); - - // sort index for all top volume branches - IdArrayType topVolSortForOuterSaddleIdx; - vtkm::cont::Algorithm::Copy(topVolBranchesIndex, topVolSortForOuterSaddleIdx); - vtkm::cont::Algorithm::Sort(topVolSortForOuterSaddleIdx, parentComparator); - - IdArrayType parentPermutation; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->TopVolumeBranchParent, topVolSortForOuterSaddleIdx, parentPermutation); - - // warning: when parent is NO_SUCH_ELEMENT, parentSaddleEps obtains 0 - // However, the corresponding element will be discarded in collecting outer saddles - IdArrayType parentSaddleEpsPermutation; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->BranchSaddleEpsilon, parentPermutation, parentSaddleEpsPermutation); - - // Some branches have parent=NO_SUCH_ELEMENT (no parent) - // we collect the isovalue of the first and/or the last branches for each parent branch - // we collect the first if branchSaddleEpsilon(parent) < 0 - // or the last if branchSaddleEpsilon(parent) > 0 - // or both if branchSaddleEpsilon(parent) == 0 - vtkm::worklet::scalar_topology::select_top_volume_contours::CollectOuterSaddle - collectOuterSaddleWorklet; - invoke( - collectOuterSaddleWorklet, parentSaddleEpsPermutation, parentPermutation, IsOuterSaddle); - - // after sorting by index back - // each top volume branch know whether it is the outer saddle of its parent - vtkm::cont::Algorithm::SortByKey(topVolSortForOuterSaddleIdx, IsOuterSaddle); - - // collect branches that need contours on extra minima/maxima - // we store the information of the parent branches (on both directions) - IdArrayType extraMaximaParentBranch; - IdArrayType extraMinimaParentBranch; - IdArrayType extraMaximaParentBranchRootGRId; - IdArrayType extraMinimaParentBranchRootGRId; - - IdArrayType allBranchGRIdByVolume; - IdArrayType branchGRIdByVolumeIdx; - - // we need global branch order including the root branch - // this information should be consistent globally - allBranchGRIdByVolume.Allocate(nTopVolBranches + 1); - vtkm::cont::Algorithm::CopySubRange( - b->TopVolumeBranchRootGRId, 0, nTopVolBranches, allBranchGRIdByVolume, 1); - auto topBranchGRIdWritePortal = allBranchGRIdByVolume.WritePortal(); - auto sortedBranchByVolPortal = b->SortedBranchByVolume.ReadPortal(); - auto branchGRIdReadPortal = b->BranchRootGRId.ReadPortal(); - topBranchGRIdWritePortal.Set(0, branchGRIdReadPortal.Get(sortedBranchByVolPortal.Get(0))); - - vtkm::cont::Algorithm::Copy( - vtkm::cont::ArrayHandleIndex(allBranchGRIdByVolume.GetNumberOfValues()), - branchGRIdByVolumeIdx); - vtkm::cont::Algorithm::SortByKey(allBranchGRIdByVolume, branchGRIdByVolumeIdx); - - vtkm::cont::Algorithm::CopyIf( - b->TopVolumeBranchParent, - IsOuterSaddle, - extraMaximaParentBranch, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); - vtkm::cont::Algorithm::CopyIf( - b->TopVolumeBranchParent, - IsOuterSaddle, - extraMinimaParentBranch, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); - - if (extraMaximaParentBranch.GetNumberOfValues()) - { - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - upperLocalEndIds, extraMaximaParentBranch, b->ExtraMaximaBranchUpperEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - lowerLocalEndIds, extraMaximaParentBranch, b->ExtraMaximaBranchLowerEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->BranchRootGRId, extraMaximaParentBranch, extraMaximaParentBranchRootGRId); - - InArrayHandleType extraMaximaBranchIsoValue; - vtkm::cont::Algorithm::CopyIf( - b->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), - IsOuterSaddle, - extraMaximaBranchIsoValue, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); - b->ExtraMaximaBranchIsoValue = extraMaximaBranchIsoValue; - - // a worklet to binary search a number in a sorted array and return the index - vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet - getParentBranchOrder; - IdArrayType permutedExtraMaximaBranchOrder; - permutedExtraMaximaBranchOrder.Allocate( - extraMaximaParentBranchRootGRId.GetNumberOfValues()); - vtkm::cont::ArrayHandleDiscard discard; - discard.Allocate(extraMaximaParentBranchRootGRId.GetNumberOfValues()); - invoke(getParentBranchOrder, - extraMaximaParentBranchRootGRId, - allBranchGRIdByVolume, - discard, - permutedExtraMaximaBranchOrder); - - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - branchGRIdByVolumeIdx, permutedExtraMaximaBranchOrder, b->ExtraMaximaBranchOrder); - } - - if (extraMinimaParentBranch.GetNumberOfValues()) - { - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - upperLocalEndIds, extraMinimaParentBranch, b->ExtraMinimaBranchUpperEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - lowerLocalEndIds, extraMinimaParentBranch, b->ExtraMinimaBranchLowerEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - b->BranchRootGRId, extraMinimaParentBranch, extraMinimaParentBranchRootGRId); - - InArrayHandleType extraMinimaBranchIsoValue; - vtkm::cont::Algorithm::CopyIf( - b->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), - IsOuterSaddle, - extraMinimaBranchIsoValue, - vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); - b->ExtraMinimaBranchIsoValue = extraMinimaBranchIsoValue; - - vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet - getParentBranchOrder; - IdArrayType permutedExtraMinimaBranchOrder; - permutedExtraMinimaBranchOrder.Allocate( - extraMinimaParentBranchRootGRId.GetNumberOfValues()); - vtkm::cont::ArrayHandleDiscard discard; - discard.Allocate(extraMinimaParentBranchRootGRId.GetNumberOfValues()); - invoke(getParentBranchOrder, - extraMinimaParentBranchRootGRId, - allBranchGRIdByVolume, - discard, - permutedExtraMinimaBranchOrder); - - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - branchGRIdByVolumeIdx, permutedExtraMinimaBranchOrder, b->ExtraMinimaBranchOrder); - } - }; - b->TopVolumeBranchSaddleIsoValue - .CastAndCallForTypes( - resolveBranchParent); + // // let's check which top volume branches are known by the block + // // We check the branchGRId of top volume branches to see whether there are matches within the block + // const vtkm::Id nTopVolBranches = b->TopVolumeBranchLowerEndGRId.GetNumberOfValues(); + // + // // sortedBranchOrder: the branch order (in the ascending order of branch root) + // // the high-level idea is to sort the branch root global regular ids + // // and for each top-volume branch, we use binary search to get the original branch index + // // if the top-volume branch does not exist in the block, it will be dropped out + // IdArrayType sortedBranchGRId; + // IdArrayType sortedBranchOrder; + // vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(b->BranchRootGRId.GetNumberOfValues()), + // sortedBranchOrder); + // vtkm::cont::Algorithm::Copy(b->BranchRootGRId, sortedBranchGRId); + // vtkm::cont::Algorithm::SortByKey(sortedBranchGRId, sortedBranchOrder); + // + // b->TopVolBranchKnownByBlockStencil.Allocate(nTopVolBranches); + // b->TopVolBranchGROrder.Allocate(nTopVolBranches); + // + // // We reuse the IdxIfWithinBlockWorklet. + // // This worklet searches for given values in a sorted array and returns the stencil & index if the value exists in the array. + // // b->TopVolBranchGROrder: the order of the topVolBranch among all known branches if the branch is known by the block. + // auto idxIfBranchWithinBlockWorklet = + // vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); + // invoke(idxIfBranchWithinBlockWorklet, + // b->TopVolumeBranchRootGRId, + // sortedBranchGRId, + // b->TopVolBranchKnownByBlockStencil, + // b->TopVolBranchGROrder); + // + // // Dropping out top-volume branches that are not known by the block. + // + // // the index of top-volume branches known by the block among all top-volume branches + // IdArrayType topVolBranchKnownByBlockIndex; + // vtkm::cont::ArrayHandleIndex topVolBranchesIndex(nTopVolBranches); + // vtkm::cont::Algorithm::CopyIf( + // topVolBranchesIndex, b->TopVolBranchKnownByBlockStencil, topVolBranchKnownByBlockIndex); + // + // const vtkm::Id nTopVolBranchKnownByBlock = topVolBranchKnownByBlockIndex.GetNumberOfValues(); + // + // // filtered b->TopVolBranchGROrder, by removing NO_SUCH_ELEMENT + // IdArrayType topVolBranchFilteredGROrder; + // + // // b->TopVolBranchInfoActualIndex: the information index of the top-volume branch + // vtkm::cont::Algorithm::CopyIf( + // b->TopVolBranchGROrder, b->TopVolBranchKnownByBlockStencil, topVolBranchFilteredGROrder); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // sortedBranchOrder, topVolBranchFilteredGROrder, b->TopVolBranchInfoActualIndex); + // + // // filtered branch saddle epsilons, global lower/upper end GR ids, + // IdArrayType topVolFilteredBranchSaddleEpsilon; + // IdArrayType topVolFilteredLowerEndGRId; + // IdArrayType topVolFilteredUpperEndGRId; + // vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchSaddleEpsilon, + // b->TopVolBranchKnownByBlockStencil, + // topVolFilteredBranchSaddleEpsilon); + // vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchUpperEndGRId, + // b->TopVolBranchKnownByBlockStencil, + // topVolFilteredUpperEndGRId); + // vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchLowerEndGRId, + // b->TopVolBranchKnownByBlockStencil, + // topVolFilteredLowerEndGRId); + // + // // for each top-vol branch known by the block + // // we get their upper end and lower end local ids + // IdArrayType topVolBranchUpperLocalEnd; + // IdArrayType topVolBranchLowerLocalEnd; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // upperLocalEndIds, b->TopVolBranchInfoActualIndex, topVolBranchUpperLocalEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // lowerLocalEndIds, b->TopVolBranchInfoActualIndex, topVolBranchLowerLocalEnd); + // + // IdArrayType topVolLowerLocalEndGRId; + // IdArrayType topVolUpperLocalEndGRId; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // globalRegularIds, topVolBranchLowerLocalEnd, topVolLowerLocalEndGRId); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // globalRegularIds, topVolBranchUpperLocalEnd, topVolUpperLocalEndGRId); + // + // // Below is the code to compute the branch hierarchy of top-volume branches + // // We need this information because we not only want to visualize the contour + // // on top-volume branches, but also on their parent branches. + // // Because we use volume as the metric, the parent branch of a top-volume branch + // // is either a top-volume branch or the root branch (where both ends are leaf nodes) + // vtkm::worklet::scalar_topology::select_top_volume_contours::BranchSaddleIsKnownWorklet + // branchSaddleIsKnownWorklet; + // // the branch saddle local ID if the saddle end is known by the block + // IdArrayType branchSaddleIsKnown; + // branchSaddleIsKnown.Allocate(nTopVolBranchKnownByBlock); + // + // invoke(branchSaddleIsKnownWorklet, // worklet + // topVolFilteredLowerEndGRId, // input + // topVolBranchLowerLocalEnd, // input + // topVolLowerLocalEndGRId, // input + // topVolFilteredUpperEndGRId, // input + // topVolBranchUpperLocalEnd, // input + // topVolUpperLocalEndGRId, // input + // topVolFilteredBranchSaddleEpsilon, // input + // branchSaddleIsKnown); // output + // + // // the order of top volume branches with parents known by the block + // IdArrayType topVolChildBranch; + // IdArrayType topVolChildBranchSaddle; + // + // vtkm::cont::Algorithm::CopyIf( + // topVolBranchKnownByBlockIndex, + // branchSaddleIsKnown, + // topVolChildBranch, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); + // vtkm::cont::Algorithm::CopyIf( + // branchSaddleIsKnown, + // branchSaddleIsKnown, + // topVolChildBranchSaddle, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); + // + // const vtkm::Id nChildBranch = topVolChildBranch.GetNumberOfValues(); + // + // // to compute the parent branch, we need to + // // 1. for the branch saddle end, collect all superarcs involving it + // // 2. get the branch information for selected superarcs + // // 3. eliminate branch information for branches sharing the same saddle end + // + // auto superarcs = + // ds.GetField("Superarcs").GetData().AsArrayHandle>(); + // auto branchRoots = + // ds.GetField("BranchRoots").GetData().AsArrayHandle>(); + // VTKM_ASSERT(superarcs.GetNumberOfValues() == branchRoots.GetNumberOfValues()); + // + // // we sort all superarcs by target to allow binary search + // IdArrayType superarcsByTarget; + // vtkm::worklet::scalar_topology::select_top_volume_contours::SuperarcTargetComparator + // superarcComparator(superarcs); + // vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(superarcs.GetNumberOfValues()), + // superarcsByTarget); + // vtkm::cont::Algorithm::Sort(superarcsByTarget, superarcComparator); + // + // IdArrayType permutedSuperarcs; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // superarcs, superarcsByTarget, permutedSuperarcs); + // + // IdArrayType permutedBranchRoots; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // branchRoots, superarcsByTarget, permutedBranchRoots); + // + // // the branch root of the superarc of the branch saddle supernode + // IdArrayType topVolChildBranchSaddleBranchRoot; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // branchRoots, topVolChildBranchSaddle, topVolChildBranchSaddleBranchRoot); + // + // // the GR Ids of the superarc of the branch saddle supernode + // IdArrayType topVolChildBranchSaddleGRIds; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // globalRegularIds, topVolChildBranchSaddle, topVolChildBranchSaddleGRIds); + // + // // there is a debate to find all superarcs connect to a supernode + // // strategy 1. iterate through saddles and parallelize over superarcs for search + // // time complexity: O(nTopVolBranches) + // // (nTopVolBranches usually <= 100, based on input parameter setting) + // // + // // strategy 2. parallelize over all saddles and use binary search to find superarcs + // // time complexity: O(log_2(nSuperarcs)) (nSuperarcs can be considerably large) + // // + // // here, we choose strategy 2 for better extendability to high nTopVolBranches + // // but in tests using nTopVolBranches <= 10, strategy 1 is generally faster + // + // // note: after getting the branch root superarc, we use binary search to get the branch order + // // because BranchRootByBranch is sorted by branch root (superarc) id + // + //#ifdef DEBUG_PRINT + // std::stringstream parentBranchStream; + // + // parentBranchStream << "Debug for Parent Branch, Block " << b->LocalBlockNo << std::endl; + // vtkm::worklet::contourtree_augmented::PrintHeader(nChildBranch, parentBranchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Child Branch Saddle", topVolChildBranchSaddle, -1, parentBranchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Child Saddle Root", topVolChildBranchSaddleBranchRoot, -1, parentBranchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Child Saddle GR Id", topVolChildBranchSaddleGRIds, -1, parentBranchStream); + // + // vtkm::worklet::contourtree_augmented::PrintHeader(superarcs.GetNumberOfValues(), + // parentBranchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Permuted Superarcs", permutedSuperarcs, -1, parentBranchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Permuted Branch roots", permutedBranchRoots, -1, parentBranchStream); + // + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, parentBranchStream.str()); + //#endif // DEBUG_PRINT + // + // // the corresponding parent branch of child branches + // IdArrayType topVolChildBranchParent; + // topVolChildBranchParent.AllocateAndFill(nChildBranch, + // vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT); + // vtkm::worklet::scalar_topology::select_top_volume_contours::GetParentBranchWorklet + // getParentBranchWorklet; + // invoke(getParentBranchWorklet, + // topVolChildBranchSaddle, + // topVolChildBranchSaddleBranchRoot, + // topVolChildBranchSaddleGRIds, + // permutedSuperarcs, + // permutedBranchRoots, + // b->BranchRootByBranch, + // upperEndGRIds, + // lowerEndGRIds, + // topVolChildBranchParent); + // + // b->TopVolumeBranchParent.AllocateAndFill( + // nTopVolBranches, vtkm::Id(vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)); + // + // vtkm::worklet::scalar_topology::select_top_volume_contours::AssignValueByIndex + // assignParentBranch; + // // for each top volume branch, assign the parent branch info id in the block + // invoke( + // assignParentBranch, topVolChildBranch, topVolChildBranchParent, b->TopVolumeBranchParent); + // // for each branch, assign true if it is a parent branch + // invoke(assignParentBranch, + // topVolChildBranchParent, + // vtkm::cont::ArrayHandleConstant(true, nChildBranch), + // b->IsParentBranch); + // + // // sort all top-volume branches based on + // // 1. parent branch info id: b->TopVolumeBranchParent + // // 2. saddle-end value: b->TopVolumeBranchSaddleIsovalue + // // 3. branch root global regular id (anything that can break tie) + // auto resolveBranchParent = [&](const auto& inArray) + // { + // using InArrayHandleType = std::decay_t; + // using ValueType = typename InArrayHandleType::ValueType; + // + // vtkm::worklet::scalar_topology::select_top_volume_contours::BranchParentComparator + // parentComparator(b->TopVolumeBranchParent, inArray, b->TopVolumeBranchRootGRId); + // + // // sort index for all top volume branches + // IdArrayType topVolSortForOuterSaddleIdx; + // vtkm::cont::Algorithm::Copy(topVolBranchesIndex, topVolSortForOuterSaddleIdx); + // vtkm::cont::Algorithm::Sort(topVolSortForOuterSaddleIdx, parentComparator); + // + // IdArrayType parentPermutation; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->TopVolumeBranchParent, topVolSortForOuterSaddleIdx, parentPermutation); + // + // // warning: when parent is NO_SUCH_ELEMENT, parentSaddleEps obtains 0 + // // However, the corresponding element will be discarded in collecting outer saddles + // IdArrayType parentSaddleEpsPermutation; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->BranchSaddleEpsilon, parentPermutation, parentSaddleEpsPermutation); + // + // // Some branches have parent=NO_SUCH_ELEMENT (no parent) + // // we collect the isovalue of the first and/or the last branches for each parent branch + // // we collect the first if branchSaddleEpsilon(parent) < 0 + // // or the last if branchSaddleEpsilon(parent) > 0 + // // or both if branchSaddleEpsilon(parent) == 0 + // IdArrayType IsOuterSaddle; + // IsOuterSaddle.Allocate(nTopVolBranches); + // vtkm::worklet::scalar_topology::select_top_volume_contours::CollectOuterSaddle + // collectOuterSaddleWorklet; + // invoke( + // collectOuterSaddleWorklet, parentSaddleEpsPermutation, parentPermutation, IsOuterSaddle); + // + // // after sorting by index back + // // each top volume branch know whether it is the outer saddle of its parent + // vtkm::cont::Algorithm::SortByKey(topVolSortForOuterSaddleIdx, IsOuterSaddle); + // + // // collect branches that need contours on extra minima/maxima + // // we store the information of the parent branches (on both directions) + // IdArrayType extraMaximaParentBranch; + // IdArrayType extraMinimaParentBranch; + // IdArrayType extraMaximaParentBranchRootGRId; + // IdArrayType extraMinimaParentBranchRootGRId; + // + // IdArrayType allBranchGRIdByVolume; + // IdArrayType branchGRIdByVolumeIdx; + // + // // we need global branch order including the root branch + // // this information should be consistent globally + // allBranchGRIdByVolume.Allocate(nTopVolBranches + 1); + // vtkm::cont::Algorithm::CopySubRange( + // b->TopVolumeBranchRootGRId, 0, nTopVolBranches, allBranchGRIdByVolume, 1); + // auto topBranchGRIdWritePortal = allBranchGRIdByVolume.WritePortal(); + // auto sortedBranchByVolPortal = b->SortedBranchByVolume.ReadPortal(); + // auto branchGRIdReadPortal = b->BranchRootGRId.ReadPortal(); + // topBranchGRIdWritePortal.Set(0, branchGRIdReadPortal.Get(sortedBranchByVolPortal.Get(0))); + // + // vtkm::cont::Algorithm::Copy( + // vtkm::cont::ArrayHandleIndex(allBranchGRIdByVolume.GetNumberOfValues()), + // branchGRIdByVolumeIdx); + // vtkm::cont::Algorithm::SortByKey(allBranchGRIdByVolume, branchGRIdByVolumeIdx); + // + // vtkm::cont::Algorithm::CopyIf( + // b->TopVolumeBranchParent, + // IsOuterSaddle, + // extraMaximaParentBranch, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); + // vtkm::cont::Algorithm::CopyIf( + // b->TopVolumeBranchParent, + // IsOuterSaddle, + // extraMinimaParentBranch, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); + // + // if (extraMaximaParentBranch.GetNumberOfValues()) + // { + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // upperLocalEndIds, extraMaximaParentBranch, b->ExtraMaximaBranchUpperEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // lowerLocalEndIds, extraMaximaParentBranch, b->ExtraMaximaBranchLowerEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->BranchRootGRId, extraMaximaParentBranch, extraMaximaParentBranchRootGRId); + // + // InArrayHandleType extraMaximaBranchIsoValue; + // vtkm::cont::Algorithm::CopyIf( + // b->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), + // IsOuterSaddle, + // extraMaximaBranchIsoValue, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); + // b->ExtraMaximaBranchIsoValue = extraMaximaBranchIsoValue; + // + // // a worklet to binary search a number in a sorted array and return the index + // vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet + // getParentBranchOrder; + // IdArrayType permutedExtraMaximaBranchOrder; + // permutedExtraMaximaBranchOrder.Allocate( + // extraMaximaParentBranchRootGRId.GetNumberOfValues()); + // vtkm::cont::ArrayHandleDiscard discard; + // discard.Allocate(extraMaximaParentBranchRootGRId.GetNumberOfValues()); + // invoke(getParentBranchOrder, + // extraMaximaParentBranchRootGRId, + // allBranchGRIdByVolume, + // discard, + // permutedExtraMaximaBranchOrder); + // + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // branchGRIdByVolumeIdx, permutedExtraMaximaBranchOrder, b->ExtraMaximaBranchOrder); + // } + // + // if (extraMinimaParentBranch.GetNumberOfValues()) + // { + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // upperLocalEndIds, extraMinimaParentBranch, b->ExtraMinimaBranchUpperEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // lowerLocalEndIds, extraMinimaParentBranch, b->ExtraMinimaBranchLowerEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // b->BranchRootGRId, extraMinimaParentBranch, extraMinimaParentBranchRootGRId); + // + // InArrayHandleType extraMinimaBranchIsoValue; + // vtkm::cont::Algorithm::CopyIf( + // b->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), + // IsOuterSaddle, + // extraMinimaBranchIsoValue, + // vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); + // b->ExtraMinimaBranchIsoValue = extraMinimaBranchIsoValue; + // + // vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet + // getParentBranchOrder; + // IdArrayType permutedExtraMinimaBranchOrder; + // permutedExtraMinimaBranchOrder.Allocate( + // extraMinimaParentBranchRootGRId.GetNumberOfValues()); + // vtkm::cont::ArrayHandleDiscard discard; + // discard.Allocate(extraMinimaParentBranchRootGRId.GetNumberOfValues()); + // invoke(getParentBranchOrder, + // extraMinimaParentBranchRootGRId, + // allBranchGRIdByVolume, + // discard, + // permutedExtraMinimaBranchOrder); + // + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // branchGRIdByVolumeIdx, permutedExtraMinimaBranchOrder, b->ExtraMinimaBranchOrder); + // } + // }; + // b->TopVolumeBranchSaddleIsoValue + // .CastAndCallForTypes( + // resolveBranchParent); }); timingsStream << " " << std::setw(60) << std::left << "ComputeTopVolumeBranchHierarchy" @@ -622,602 +615,605 @@ VTKM_CONT vtkm::cont::PartitionedDataSet SelectTopVolumeContoursFilter::DoExecut timer.Start(); // For each block, we compute the get the extracted isosurface for every selected branch - // Let's create a mesh dataset to extract all the contours first branch_top_volume_master.foreach ([&](SelectTopVolumeContoursBlock* b, const vtkmdiy::Master::ProxyWithLink&) { - using vtkm::worklet::contourtree_augmented::IdArrayType; - - // if no branch to extract from - if (b->TopVolumeBranchRootGRId.GetNumberOfValues() < 1) - return; - const vtkm::cont::DataSet& ds = input.GetPartition(b->LocalBlockNo); - - // branch root global regular ID - // size: nBranches - // usage: identifier of the branch - vtkm::cont::Algorithm::Copy( - ds.GetField("BranchRootGRId").GetData().AsArrayHandle(), b->BranchRootGRId); - - // branch local upper end and lower end - // size: nBranches - // usage: search for the superarc of an arbitrary point (not necessarily on grid) - auto upperEndLocalIds = ds.GetField("UpperEndLocalIds").GetData().AsArrayHandle(); - auto lowerEndLocalIds = ds.GetField("LowerEndLocalIds").GetData().AsArrayHandle(); - - // global regular ids - auto globalRegularIds = - ds.GetField("RegularNodeGlobalIds").GetData().AsArrayHandle(); - - vtkm::Id3 globalPointDimensions; - vtkm::Id3 pointDimensions, globalPointIndexStart; - - ds.GetCellSet().CastAndCallForTypes( - vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), - pointDimensions, - globalPointDimensions, - globalPointIndexStart); - -#ifdef DEBUG_PRINT - VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block size info"); - std::stringstream rs; - { - rs << "globalPointDimensions: " << globalPointDimensions << std::endl; - rs << "pointDimensions: " << pointDimensions << std::endl; - rs << "globalPointIndexStart: " << globalPointIndexStart << std::endl; - rs << "globalRegularIDs: " << globalRegularIds.GetNumberOfValues() << std::endl; - // ds.PrintSummary(rs); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); - } -#endif - - // Tool to relabel local mesh ids to global ids - auto localToGlobalIdRelabeler = vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler( - globalPointIndexStart, pointDimensions, globalPointDimensions); - IdArrayType globalIdsByMesh; - - // Note: the cell set is different from the mesh structure - // Here, we assume that the cell set is structured grid - // A more general way to do this is to use CellSet().GetCellPointIds(i) to extract all the local ids and keep unique ones - // local ids in the mesh - IdArrayType localIdsByMesh; - vtkm::cont::Algorithm::Copy( - vtkm::cont::ArrayHandleIndex(pointDimensions[0] * pointDimensions[1] * pointDimensions[2]), - localIdsByMesh); - // then, we transform the local ids to global ids - auto localTransformToGlobalId = - vtkm::cont::make_ArrayHandleTransform(localIdsByMesh, localToGlobalIdRelabeler); - vtkm::cont::ArrayCopyDevice(localTransformToGlobalId, globalIdsByMesh); - - // detect whether the element in globalRegularIds are in the block - // globalIdsDiscard is just a filler for the worklet format. We do not use it. - // The last slot for the worklet is useful in a later step. - // Here we just reuse the worklet - IdArrayType globalIdsWithinBlockStencil; - vtkm::cont::ArrayHandleDiscard globalIdsDiscard; - globalIdsWithinBlockStencil.Allocate(globalRegularIds.GetNumberOfValues()); - globalIdsDiscard.Allocate(globalRegularIds.GetNumberOfValues()); - - vtkm::cont::Invoker invoke; - // stencil is 1 if the global regular id is within the block, 0 otherwise - auto idxIfWithinBlock = - vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); - invoke(idxIfWithinBlock, - globalRegularIds, - globalIdsByMesh, - globalIdsWithinBlockStencil, - globalIdsDiscard); - - auto resolveArray = [&](const auto& inArray) { - using InArrayHandleType = std::decay_t; - using ValueType = typename InArrayHandleType::ValueType; - - // we need to sort all values based on the global ids - // and remove values of points that do not belong to the local block - auto dataValues = ds.GetField("DataValues").GetData().AsArrayHandle(); - - IdArrayType globalIdsWithinBlock; - IdArrayType localIdsWithinBlock; - InArrayHandleType dataValuesWithinBlock; - - // filter global regular ids, array ids, and data values - vtkm::cont::Algorithm::CopyIf( - globalRegularIds, globalIdsWithinBlockStencil, globalIdsWithinBlock); - vtkm::cont::Algorithm::CopyIf( - vtkm::cont::ArrayHandleIndex(globalRegularIds.GetNumberOfValues()), - globalIdsWithinBlockStencil, - localIdsWithinBlock); - vtkm::cont::Algorithm::CopyIf(dataValues, globalIdsWithinBlockStencil, dataValuesWithinBlock); - - // sorted index based on global regular ids - IdArrayType sortedGlobalIds; - vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(globalIdsByMesh.GetNumberOfValues()), - sortedGlobalIds); - vtkm::cont::Algorithm::SortByKey(globalIdsWithinBlock, sortedGlobalIds); - - // globalIdsWithinBlock (sorted) and globalIdsByMesh should be identical. - // computing globalIdsWithinBlock ensures the input data is correct - bool identical = - globalIdsWithinBlock.GetNumberOfValues() == globalIdsByMesh.GetNumberOfValues(); - if (identical) - { - vtkm::cont::ArrayHandle globalIdsIdentical; - vtkm::cont::Algorithm::Transform( - globalIdsWithinBlock, globalIdsByMesh, globalIdsIdentical, vtkm::Equal()); - identical = vtkm::cont::Algorithm::Reduce(globalIdsIdentical, true, vtkm::LogicalAnd()); - } - if (!identical) - { - vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsByMesh.GetNumberOfValues()); - vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsByMesh", globalIdsByMesh); - vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsWithinBlock.GetNumberOfValues()); - vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsWithinBlock", - globalIdsWithinBlock); - } - VTKM_ASSERT(identical); - - // filtered and sorted local node info ids - // i.e. index of global regular ids, data values, and superparents - // Note: This is not local mesh id! Make sure to distinguish them - IdArrayType sortedLocalNodeInfoIdsWithinBlock; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - localIdsWithinBlock, sortedGlobalIds, sortedLocalNodeInfoIdsWithinBlock); - - // sorted data values - InArrayHandleType sortedDataValuesWithinBlock; - vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( - dataValuesWithinBlock, sortedGlobalIds, sortedDataValuesWithinBlock); - - // create an execution object to find the superarc for an arbitrary point within the mesh - // all information below are required to initialize the execution object - auto superparents = ds.GetField("Superparents").GetData().AsArrayHandle(); - auto supernodes = ds.GetField("Supernodes").GetData().AsArrayHandle(); - auto superarcs = ds.GetField("Superarcs").GetData().AsArrayHandle(); - auto superchildren = ds.GetField("Superchildren").GetData().AsArrayHandle(); - auto whichRound = ds.GetField("WhichRound").GetData().AsArrayHandle(); - auto whichIteration = ds.GetField("WhichIteration").GetData().AsArrayHandle(); - auto hyperparents = ds.GetField("Hyperparents").GetData().AsArrayHandle(); - auto hypernodes = ds.GetField("Hypernodes").GetData().AsArrayHandle(); - auto hyperarcs = ds.GetField("Hyperarcs").GetData().AsArrayHandle(); - - // filtered + sorted superparents of nodes - IdArrayType superparentsWithinBlock; - vtkm::cont::Algorithm::CopyIf( - superparents, globalIdsWithinBlockStencil, superparentsWithinBlock); - IdArrayType sortedSuperparentsWithinBlock; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - superparentsWithinBlock, sortedGlobalIds, sortedSuperparentsWithinBlock); - - // initialize the exec object - // note: terms should include the contour tree as much as possible - // should pass the full arrays to the object instead of the filtered ones - auto findSuperarcForNode = - vtkm::worklet::contourtree_distributed::FindSuperArcForUnknownNode( - superparents, - supernodes, - superarcs, - superchildren, - whichRound, - whichIteration, - hyperparents, - hypernodes, - hyperarcs, - globalRegularIds, - dataValues); - - // let's check which branches are known by the block - // We check the branchGRId of top volume branches to see whether there are matches within the block - vtkm::Id nIsoValues = inArray.GetNumberOfValues(); - vtkm::Id totalNumPoints = - globalPointDimensions[0] * globalPointDimensions[1] * globalPointDimensions[2]; - - // dropping out top-volume branches that are not known by the block - // index of top-volume branches within the block among all top-volume branches - - IdArrayType topVolBranchWithinBlockId; - vtkm::cont::Algorithm::CopyIf(vtkm::cont::ArrayHandleIndex(nIsoValues), - b->TopVolBranchKnownByBlockStencil, - topVolBranchWithinBlockId); - auto topVolBranchWithinBlockIdPortal = topVolBranchWithinBlockId.ReadPortal(); - - vtkm::Id nTopVolBranchWithinBlock = topVolBranchWithinBlockId.GetNumberOfValues(); - - // filtered branch saddle values - InArrayHandleType isoValues; - vtkm::cont::Algorithm::CopyIf(inArray, b->TopVolBranchKnownByBlockStencil, isoValues); - auto isoValuePortal = isoValues.ReadPortal(); - - // filtered branch saddle epsilons - IdArrayType topVolBranchSaddleEpsilons; - vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchSaddleEpsilon, - b->TopVolBranchKnownByBlockStencil, - topVolBranchSaddleEpsilons); - auto topVolBranchSaddleEpsilonPortal = topVolBranchSaddleEpsilons.ReadPortal(); - - // for each top-vol branch in the block - // we get their upper end and lower end local ids - IdArrayType topVolLocalBranchUpperEnd; - IdArrayType topVolLocalBranchLowerEnd; - vtkm::cont::ArrayHandle topVolIsParent; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - upperEndLocalIds, b->TopVolBranchInfoActualIndex, topVolLocalBranchUpperEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - lowerEndLocalIds, b->TopVolBranchInfoActualIndex, topVolLocalBranchLowerEnd); - vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex>( - b->IsParentBranch, b->TopVolBranchInfoActualIndex, topVolIsParent); - auto topVolIsParentPortal = topVolIsParent.ReadPortal(); - - // we compute the superarc of the branch within the block - // around the given isovalue - IdArrayType branchIsoSuperarcs; - branchIsoSuperarcs.Allocate(nTopVolBranchWithinBlock); - - auto branchIsoSuperarcWorklet = - vtkm::worklet::scalar_topology::select_top_volume_contours::GetSuperarcByIsoValueWorklet( - totalNumPoints); - invoke(branchIsoSuperarcWorklet, - topVolLocalBranchUpperEnd, - topVolLocalBranchLowerEnd, - isoValues, - topVolBranchSaddleEpsilons, - branchIsoSuperarcs, - findSuperarcForNode); - auto branchIsoSuperarcsPortal = branchIsoSuperarcs.ReadPortal(); - -#ifdef DEBUG_PRINT - std::stringstream branchStream; - - branchStream << "Debug for branch info, Block " << b->LocalBlockNo << std::endl; - vtkm::worklet::contourtree_augmented::PrintHeader(sortedBranchGRId.GetNumberOfValues(), - branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Raw Branch GR", b->BranchRootGRId, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Raw Upper End", upperLocalEndIds, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Raw Lower End", lowerLocalEndIds, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Sorted Branch GR", sortedBranchGRId, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Sorted Branch Id", sortedBranchOrder, -1, branchStream); - - vtkm::worklet::contourtree_augmented::PrintHeader(nIsoValues, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Top Branch GR", b->TopVolumeBranchRootGRId, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Top Branch Stencil", topVolBranchWithinBlockStencil, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Top Branch Idx", topVolBranchInfoIdx, -1, branchStream); - - vtkm::worklet::contourtree_augmented::PrintHeader(nTopVolBranchKnownByBlock, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Sorted Upper End", topVolBranchUpperLocalEnd, -1, branchStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Sorted Lower End", topVolBranchLowerLocalEnd, -1, branchStream); - - VTKM_LOG_S(vtkm::cont::LogLevel::Info, branchStream.str()); -#endif - - const vtkm::Id nExtraMaximaBranch = b->ExtraMaximaBranchLowerEnd.GetNumberOfValues(); - const vtkm::Id nExtraMinimaBranch = b->ExtraMinimaBranchLowerEnd.GetNumberOfValues(); - InArrayHandleType extraMaximaBranchIsoValue; - InArrayHandleType extraMinimaBranchIsoValue; - - IdArrayType extraMaximaBranchSuperarcs; - IdArrayType extraMinimaBranchSuperarcs; - extraMaximaBranchSuperarcs.Allocate(nExtraMaximaBranch); - extraMinimaBranchSuperarcs.Allocate(nExtraMinimaBranch); - - if (nExtraMaximaBranch) - { - extraMaximaBranchIsoValue = b->ExtraMaximaBranchIsoValue.AsArrayHandle(); - - invoke(branchIsoSuperarcWorklet, - b->ExtraMaximaBranchUpperEnd, - b->ExtraMaximaBranchLowerEnd, - extraMaximaBranchIsoValue, - vtkm::cont::ArrayHandleConstant(1, nExtraMaximaBranch), - extraMaximaBranchSuperarcs, - findSuperarcForNode); - -#ifdef DEBUG_PRINT - std::stringstream extraMaxStream; - extraMaxStream << "Debug for Extra Maxima Branch, Block " << b->LocalBlockNo << std::endl; - vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMaximaBranch, extraMaxStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Upper End", b->ExtraMaximaBranchUpperEnd, -1, extraMaxStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Lower End", b->ExtraMaximaBranchLowerEnd, -1, extraMaxStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "Max Branch IsoValue", extraMaximaBranchIsoValue, -1, extraMaxStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Superarc", extraMaximaBranchSuperarcs, -1, extraMaxStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Order", b->ExtraMaximaBranchOrder, -1, extraMaxStream); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMaxStream.str()); -#endif // DEBUG_PRINT - } - - if (nExtraMinimaBranch) - { - extraMinimaBranchIsoValue = b->ExtraMinimaBranchIsoValue.AsArrayHandle(); - - invoke(branchIsoSuperarcWorklet, - b->ExtraMinimaBranchUpperEnd, - b->ExtraMinimaBranchLowerEnd, - extraMinimaBranchIsoValue, - vtkm::cont::ArrayHandleConstant(-1, nExtraMinimaBranch), - extraMinimaBranchSuperarcs, - findSuperarcForNode); - -#ifdef DEBUG_PRINT - std::stringstream extraMinStream; - extraMaxStream << "Debug for Extra Maxima Branch, Block " << b->LocalBlockNo << std::endl; - vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMinimaBranch, extraMinStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Upper End", b->ExtraMinimaBranchUpperEnd, -1, extraMinStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Lower End", b->ExtraMinimaBranchLowerEnd, -1, extraMinStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "Max Branch IsoValue", extraMinimaBranchIsoValue, -1, extraMinStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Superarc", extraMinimaBranchSuperarcs, -1, extraMinStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Max Branch Order", b->ExtraMinimaBranchOrder, -1, extraMinStream); - VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMinStream.str()); -#endif // DEBUG_PRINT - } - - auto extraMaximaBranchSuperarcPortal = extraMaximaBranchSuperarcs.ReadPortal(); - auto extraMinimaBranchSuperarcPortal = extraMinimaBranchSuperarcs.ReadPortal(); - auto extraMaximaBranchIsoValuePortal = extraMaximaBranchIsoValue.ReadPortal(); - auto extraMinimaBranchIsoValuePortal = extraMinimaBranchIsoValue.ReadPortal(); - auto extraMaximaBranchOrderPortal = b->ExtraMaximaBranchOrder.ReadPortal(); - auto extraMinimaBranchOrderPortal = b->ExtraMinimaBranchOrder.ReadPortal(); - - const vtkm::Id nContours = nTopVolBranchWithinBlock + nExtraMaximaBranch + nExtraMinimaBranch; - b->IsosurfaceEdgesOffset.AllocateAndFill(nContours, 0); - b->IsosurfaceEdgesLabels.AllocateAndFill(nContours, 0); - b->IsosurfaceEdgesOrders.AllocateAndFill(nContours, 0); - InArrayHandleType isosurfaceIsoValue; - isosurfaceIsoValue.AllocateAndFill(nContours, ValueType(0)); - auto edgeOffsetWritePortal = b->IsosurfaceEdgesOffset.WritePortal(); - auto edgeLabelWritePortal = b->IsosurfaceEdgesLabels.WritePortal(); - auto edgeOrderWritePortal = b->IsosurfaceEdgesOrders.WritePortal(); - auto isosurfaceValuePortal = isosurfaceIsoValue.WritePortal(); - - vtkm::Id nContourCandidateMeshes = 0; - // NOTE: nContours denotes the number of isosurfaces for visualization. - // The number is usually small, so linear loop is not too costly. - // NOTE update 06/16/2024: it is hard to pre-compute the superarc of mesh edges, - // because we always need the isovalue of the contour. - // As a result, we have to iterate through nContours (=O(k)). - for (vtkm::Id branchIdx = 0; branchIdx < nContours; branchIdx++) - { - ValueType isoValue; - vtkm::Id currBranchSaddleEpsilon, branchSuperarc, branchOrder, branchLabel = 0; - - using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_SADDLE; - using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_COVER; - using vtkm::worklet::scalar_topology::select_top_volume_contours::MAXIMA_CONTOUR; - - if (branchIdx < nTopVolBranchWithinBlock) - { - isoValue = isoValuePortal.Get(branchIdx); - currBranchSaddleEpsilon = topVolBranchSaddleEpsilonPortal.Get(branchIdx); - branchSuperarc = branchIsoSuperarcsPortal.Get(branchIdx); - branchOrder = topVolBranchWithinBlockIdPortal.Get(branchIdx) + 1; - branchLabel |= BRANCH_SADDLE; - branchLabel |= topVolIsParentPortal.Get(branchIdx) ? BRANCH_COVER : 0; - branchLabel |= currBranchSaddleEpsilon > 0 ? MAXIMA_CONTOUR : 0; - } - else if (branchIdx < nTopVolBranchWithinBlock + nExtraMaximaBranch) - { - const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock; - isoValue = extraMaximaBranchIsoValuePortal.Get(idx); - currBranchSaddleEpsilon = 1; - branchSuperarc = extraMaximaBranchSuperarcPortal.Get(idx); - branchOrder = extraMaximaBranchOrderPortal.Get(idx); - branchLabel |= MAXIMA_CONTOUR; - } - else - { - const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock - nExtraMaximaBranch; - VTKM_ASSERT(idx < nExtraMinimaBranch); - isoValue = extraMinimaBranchIsoValuePortal.Get(idx); - currBranchSaddleEpsilon = -1; - branchSuperarc = extraMinimaBranchSuperarcPortal.Get(idx); - branchOrder = extraMinimaBranchOrderPortal.Get(idx); - } - - edgeOffsetWritePortal.Set(branchIdx, b->IsosurfaceEdgesFrom.GetNumberOfValues()); - edgeLabelWritePortal.Set(branchIdx, branchLabel); - edgeOrderWritePortal.Set(branchIdx, branchOrder); - isosurfaceValuePortal.Set(branchIdx, isoValue); - - if (vtkm::worklet::contourtree_augmented::NoSuchElement(branchSuperarc)) - { - continue; - } - - // Note: by concept, there is no 3D cell if pointDimensions[2] <= 1 - vtkm::Id nCells = globalPointDimensions[2] <= 1 - ? (pointDimensions[0] - 1) * (pointDimensions[1] - 1) - : (pointDimensions[0] - 1) * (pointDimensions[1] - 1) * (pointDimensions[2] - 1); - - // we get the polarity cases of cells - // we use lookup tables to for fast processing - IdArrayType vertexOffset = globalPointDimensions[2] <= 1 - ? vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset2d - : vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset3d; - IdArrayType edgeTable = globalPointDimensions[2] <= 1 - ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTable2d - : (isMarchingCubes - ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableMC3d - : vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableLT3d); - IdArrayType numBoundTable = globalPointDimensions[2] <= 1 - ? vtkm::worklet::scalar_topology::select_top_volume_contours::numLinesTable2d - : (isMarchingCubes - ? vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableMC3d - : vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableLT3d); - IdArrayType boundaryTable = globalPointDimensions[2] <= 1 - ? vtkm::worklet::scalar_topology::select_top_volume_contours::lineTable2d - : (isMarchingCubes - ? vtkm::worklet::scalar_topology::select_top_volume_contours::triTableMC3d - : vtkm::worklet::scalar_topology::select_top_volume_contours::triTableLT3d); - IdArrayType labelEdgeTable = isMarchingCubes - ? vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableMC3d - : vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableLT3d; - - IdArrayType caseCells; - - caseCells.Allocate(nCells); - IdArrayType numEdgesInCell; - numEdgesInCell.Allocate(nCells); - auto caseCellsWorklet = - vtkm::worklet::scalar_topology::select_top_volume_contours::GetCellCasesWorklet< - ValueType>(pointDimensions, currBranchSaddleEpsilon, isoValue); - - invoke(caseCellsWorklet, - vtkm::cont::ArrayHandleIndex(nCells), - sortedDataValuesWithinBlock, - vertexOffset, - caseCells); - - // we compute the number of edges for each cell - // to initialize the array size of edges - IdArrayType numBoundariesInCell; - vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( - numBoundTable, caseCells, numBoundariesInCell); - - vtkm::cont::Algorithm::Transform(numBoundariesInCell, - globalPointDimensions[2] <= 1 - ? vtkm::cont::ArrayHandleConstant(1, nCells) - : vtkm::cont::ArrayHandleConstant(3, nCells), - numEdgesInCell, - vtkm::Multiply()); - - // use prefix sum to get the offset of the starting edge in each cell/cube - vtkm::Id nEdges = - vtkm::cont::Algorithm::Reduce(numEdgesInCell, vtkm::Id(0)); - nContourCandidateMeshes += globalPointDimensions[2] <= 1 ? nEdges : nEdges / 3; - IdArrayType edgesOffset; - vtkm::cont::Algorithm::ScanExclusive(numEdgesInCell, edgesOffset); - - vtkm::cont::ArrayHandle isosurfaceEdgesFrom; - vtkm::cont::ArrayHandle isosurfaceEdgesTo; - // IdArrayType isosurfaceEdgesLabels; - IdArrayType isValidEdges; - isosurfaceEdgesFrom.Allocate(nEdges); - isosurfaceEdgesTo.Allocate(nEdges); - // isosurfaceEdgesLabels.Allocate(nEdges); - isValidEdges.Allocate(nEdges); - - // draw isosurface - auto getEdgesInCellWorklet = - vtkm::worklet::scalar_topology::select_top_volume_contours::GetEdgesInCellWorklet< - ValueType>(pointDimensions, - globalPointIndexStart, - isoValue, - branchSuperarc, - currBranchSaddleEpsilon, - totalNumPoints, - isMarchingCubes); - - invoke(getEdgesInCellWorklet, - edgesOffset, - caseCells, - sortedLocalNodeInfoIdsWithinBlock, - sortedDataValuesWithinBlock, - vertexOffset, - edgeTable, - numBoundTable, - boundaryTable, - labelEdgeTable, - isosurfaceEdgesFrom, - isosurfaceEdgesTo, - isValidEdges, - findSuperarcForNode); - - // isValidEdges: stencil about whether the edge is on the desired superarc - vtkm::cont::ArrayHandle validEdgesFrom; - vtkm::cont::ArrayHandle validEdgesTo; - - // we remove edges that are not on the desired superarc - vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesFrom, isValidEdges, validEdgesFrom); - vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesTo, isValidEdges, validEdgesTo); - - // append edges into the result array - vtkm::Id nValidEdges = validEdgesFrom.GetNumberOfValues(); - vtkm::Id nExistEdges = branchIdx == 0 ? 0 : b->IsosurfaceEdgesFrom.GetNumberOfValues(); - - b->IsosurfaceEdgesFrom.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); - b->IsosurfaceEdgesTo.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); - -#ifdef DEBUG_PRINT - std::stringstream edgeInfoStream; - contourStream << "Debug for Contour Info, Block " << b->LocalBlockNo << std::endl; - vtkm::worklet::contourtree_augmented::PrintHeader(nContours, edgeInfoStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "edgeOffset", b->IsosurfaceEdgesOffset, -1, edgeInfoStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "edgeOrders", b->IsosurfaceEdgesOrders, -1, edgeInfoStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "edgeLabels", b->IsosurfaceEdgesLabels, -1, edgeInfoStream); - - VTKM_LOG_S(vtkm::cont::LogLevel::Info, edgeInfoStream.str()); -#endif // DEBUG_PRINT - - vtkm::cont::Algorithm::CopySubRange( - validEdgesFrom, 0, nValidEdges, b->IsosurfaceEdgesFrom, nExistEdges); - vtkm::cont::Algorithm::CopySubRange( - validEdgesTo, 0, nValidEdges, b->IsosurfaceEdgesTo, nExistEdges); - -#ifdef DEBUG_PRINT - std::stringstream contourStream; - contourStream << "Debug for Contour, Block " << b->LocalBlockNo << std::endl; - contourStream << "Branch Superarc = " << branchSuperarc << ", isoValue = " << isoValue - << std::endl; - - vtkm::worklet::contourtree_augmented::PrintHeader(nCells, contourStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Case of Cells", caseCells, -1, contourStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "# of Edges", numEdgesInCell, -1, contourStream); - - vtkm::worklet::contourtree_augmented::PrintHeader(nEdges, contourStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "Edge Stencil", isValidEdges, -1, contourStream); - - vtkm::worklet::contourtree_augmented::PrintHeader(nValidEdges, contourStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "EdgesFrom", validEdgesFrom, -1, contourStream); - vtkm::worklet::contourtree_augmented::PrintValues( - "EdgesTo", validEdgesTo, -1, contourStream); - vtkm::worklet::contourtree_augmented::PrintIndices( - "EdgesLabels", validEdgesLabels, -1, contourStream); - - VTKM_LOG_S(vtkm::cont::LogLevel::Info, contourStream.str()); -#endif // DEBUG_PRINT - } - b->IsosurfaceIsoValue = isosurfaceIsoValue; - - const vtkm::Id nMeshesOnBranches = globalPointDimensions[2] <= 1 - ? b->IsosurfaceEdgesFrom.GetNumberOfValues() - : b->IsosurfaceEdgesFrom.GetNumberOfValues() / 3; - VTKM_LOG_S(TimingsLogLevel, - std::endl - << "----------- Draw Isosurface (block=" << b->LocalBlockNo << ")------------" - << std::endl - << " " << std::setw(60) << std::left << "Number of Contours: " << nContours - << std::endl - << " " << std::setw(60) << std::left - << "Number of Isosurface Meshes: " << nContourCandidateMeshes << std::endl - << " " << std::setw(60) << std::left - << "Number of Meshes On Branches: " << nMeshesOnBranches << std::endl); - }; - b->TopVolumeBranchSaddleIsoValue - .CastAndCallForTypes(resolveArray); + b->ExtractIsosurfaceOnSelectedBranch(ds, this->GetMarchingCubes(), this->GetTimingsLogLevel()); + // using vtkm::worklet::contourtree_augmented::IdArrayType; + // + // // if no branch to extract from + // if (b->TopVolumeBranchRootGRId.GetNumberOfValues() < 1) + // return; + // + + // + // // branch root global regular ID + // // size: nBranches + // // usage: identifier of the branch + // vtkm::cont::Algorithm::Copy( + // ds.GetField("BranchRootGRId").GetData().AsArrayHandle(), b->BranchRootGRId); + // + // // branch local upper end and lower end + // // size: nBranches + // // usage: search for the superarc of an arbitrary point (not necessarily on grid) + // auto upperEndLocalIds = ds.GetField("UpperEndLocalIds").GetData().AsArrayHandle(); + // auto lowerEndLocalIds = ds.GetField("LowerEndLocalIds").GetData().AsArrayHandle(); + // + // // global regular ids + // auto globalRegularIds = + // ds.GetField("RegularNodeGlobalIds").GetData().AsArrayHandle(); + // + // vtkm::Id3 globalPointDimensions; + // vtkm::Id3 pointDimensions, globalPointIndexStart; + // + // ds.GetCellSet().CastAndCallForTypes( + // vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + // pointDimensions, + // globalPointDimensions, + // globalPointIndexStart); + // + //#ifdef DEBUG_PRINT + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block size info"); + // std::stringstream rs; + // { + // rs << "globalPointDimensions: " << globalPointDimensions << std::endl; + // rs << "pointDimensions: " << pointDimensions << std::endl; + // rs << "globalPointIndexStart: " << globalPointIndexStart << std::endl; + // rs << "globalRegularIDs: " << globalRegularIds.GetNumberOfValues() << std::endl; + // // ds.PrintSummary(rs); + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); + // } + //#endif + // + // // Tool to relabel local mesh ids to global ids + // auto localToGlobalIdRelabeler = vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler( + // globalPointIndexStart, pointDimensions, globalPointDimensions); + // IdArrayType globalIdsByMesh; + // + // // Note: the cell set is different from the mesh structure + // // Here, we assume that the cell set is structured grid + // // A more general way to do this is to use CellSet().GetCellPointIds(i) to extract all the local ids and keep unique ones + // // local ids in the mesh + // IdArrayType localIdsByMesh; + // vtkm::cont::Algorithm::Copy( + // vtkm::cont::ArrayHandleIndex(pointDimensions[0] * pointDimensions[1] * pointDimensions[2]), + // localIdsByMesh); + // // then, we transform the local ids to global ids + // auto localTransformToGlobalId = + // vtkm::cont::make_ArrayHandleTransform(localIdsByMesh, localToGlobalIdRelabeler); + // vtkm::cont::ArrayCopyDevice(localTransformToGlobalId, globalIdsByMesh); + // + // // detect whether the element in globalRegularIds are in the block + // // globalIdsDiscard is just a filler for the worklet format. We do not use it. + // // The last slot for the worklet is useful in a later step. + // // Here we just reuse the worklet + // IdArrayType globalIdsWithinBlockStencil; + // vtkm::cont::ArrayHandleDiscard globalIdsDiscard; + // globalIdsWithinBlockStencil.Allocate(globalRegularIds.GetNumberOfValues()); + // globalIdsDiscard.Allocate(globalRegularIds.GetNumberOfValues()); + // + // vtkm::cont::Invoker invoke; + // // stencil is 1 if the global regular id is within the block, 0 otherwise + // auto idxIfWithinBlock = + // vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); + // invoke(idxIfWithinBlock, + // globalRegularIds, + // globalIdsByMesh, + // globalIdsWithinBlockStencil, + // globalIdsDiscard); + // + // auto resolveArray = [&](const auto& inArray) { + // using InArrayHandleType = std::decay_t; + // using ValueType = typename InArrayHandleType::ValueType; + // + // // we need to sort all values based on the global ids + // // and remove values of points that do not belong to the local block + // auto dataValues = ds.GetField("DataValues").GetData().AsArrayHandle(); + // + // IdArrayType globalIdsWithinBlock; + // IdArrayType localIdsWithinBlock; + // InArrayHandleType dataValuesWithinBlock; + // + // // filter global regular ids, array ids, and data values + // vtkm::cont::Algorithm::CopyIf( + // globalRegularIds, globalIdsWithinBlockStencil, globalIdsWithinBlock); + // vtkm::cont::Algorithm::CopyIf( + // vtkm::cont::ArrayHandleIndex(globalRegularIds.GetNumberOfValues()), + // globalIdsWithinBlockStencil, + // localIdsWithinBlock); + // vtkm::cont::Algorithm::CopyIf(dataValues, globalIdsWithinBlockStencil, dataValuesWithinBlock); + // + // // sorted index based on global regular ids + // IdArrayType sortedGlobalIds; + // vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(globalIdsByMesh.GetNumberOfValues()), + // sortedGlobalIds); + // vtkm::cont::Algorithm::SortByKey(globalIdsWithinBlock, sortedGlobalIds); + // + // // globalIdsWithinBlock (sorted) and globalIdsByMesh should be identical. + // // computing globalIdsWithinBlock ensures the input data is correct + // bool identical = + // globalIdsWithinBlock.GetNumberOfValues() == globalIdsByMesh.GetNumberOfValues(); + // if (identical) + // { + // vtkm::cont::ArrayHandle globalIdsIdentical; + // vtkm::cont::Algorithm::Transform( + // globalIdsWithinBlock, globalIdsByMesh, globalIdsIdentical, vtkm::Equal()); + // identical = vtkm::cont::Algorithm::Reduce(globalIdsIdentical, true, vtkm::LogicalAnd()); + // } + // if (!identical) + // { + // vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsByMesh.GetNumberOfValues()); + // vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsByMesh", globalIdsByMesh); + // vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsWithinBlock.GetNumberOfValues()); + // vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsWithinBlock", + // globalIdsWithinBlock); + // } + // VTKM_ASSERT(identical); + // + // // filtered and sorted local node info ids + // // i.e. index of global regular ids, data values, and superparents + // // Note: This is not local mesh id! Make sure to distinguish them + // IdArrayType sortedLocalNodeInfoIdsWithinBlock; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // localIdsWithinBlock, sortedGlobalIds, sortedLocalNodeInfoIdsWithinBlock); + // + // // sorted data values + // InArrayHandleType sortedDataValuesWithinBlock; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + // dataValuesWithinBlock, sortedGlobalIds, sortedDataValuesWithinBlock); + // + // // create an execution object to find the superarc for an arbitrary point within the mesh + // // all information below are required to initialize the execution object + // auto superparents = ds.GetField("Superparents").GetData().AsArrayHandle(); + // auto supernodes = ds.GetField("Supernodes").GetData().AsArrayHandle(); + // auto superarcs = ds.GetField("Superarcs").GetData().AsArrayHandle(); + // auto superchildren = ds.GetField("Superchildren").GetData().AsArrayHandle(); + // auto whichRound = ds.GetField("WhichRound").GetData().AsArrayHandle(); + // auto whichIteration = ds.GetField("WhichIteration").GetData().AsArrayHandle(); + // auto hyperparents = ds.GetField("Hyperparents").GetData().AsArrayHandle(); + // auto hypernodes = ds.GetField("Hypernodes").GetData().AsArrayHandle(); + // auto hyperarcs = ds.GetField("Hyperarcs").GetData().AsArrayHandle(); + // + // // filtered + sorted superparents of nodes + // IdArrayType superparentsWithinBlock; + // vtkm::cont::Algorithm::CopyIf( + // superparents, globalIdsWithinBlockStencil, superparentsWithinBlock); + // IdArrayType sortedSuperparentsWithinBlock; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // superparentsWithinBlock, sortedGlobalIds, sortedSuperparentsWithinBlock); + // + // // initialize the exec object + // // note: terms should include the contour tree as much as possible + // // should pass the full arrays to the object instead of the filtered ones + // auto findSuperarcForNode = + // vtkm::worklet::contourtree_distributed::FindSuperArcForUnknownNode( + // superparents, + // supernodes, + // superarcs, + // superchildren, + // whichRound, + // whichIteration, + // hyperparents, + // hypernodes, + // hyperarcs, + // globalRegularIds, + // dataValues); + // + // // let's check which branches are known by the block + // // We check the branchGRId of top volume branches to see whether there are matches within the block + // vtkm::Id nIsoValues = inArray.GetNumberOfValues(); + // vtkm::Id totalNumPoints = + // globalPointDimensions[0] * globalPointDimensions[1] * globalPointDimensions[2]; + // + // // dropping out top-volume branches that are not known by the block + // // index of top-volume branches within the block among all top-volume branches + // + // IdArrayType topVolBranchWithinBlockId; + // vtkm::cont::Algorithm::CopyIf(vtkm::cont::ArrayHandleIndex(nIsoValues), + // b->TopVolBranchKnownByBlockStencil, + // topVolBranchWithinBlockId); + // auto topVolBranchWithinBlockIdPortal = topVolBranchWithinBlockId.ReadPortal(); + // + // vtkm::Id nTopVolBranchWithinBlock = topVolBranchWithinBlockId.GetNumberOfValues(); + // + // // filtered branch saddle values + // InArrayHandleType isoValues; + // vtkm::cont::Algorithm::CopyIf(inArray, b->TopVolBranchKnownByBlockStencil, isoValues); + // auto isoValuePortal = isoValues.ReadPortal(); + // + // // filtered branch saddle epsilons + // IdArrayType topVolBranchSaddleEpsilons; + // vtkm::cont::Algorithm::CopyIf(b->TopVolumeBranchSaddleEpsilon, + // b->TopVolBranchKnownByBlockStencil, + // topVolBranchSaddleEpsilons); + // auto topVolBranchSaddleEpsilonPortal = topVolBranchSaddleEpsilons.ReadPortal(); + // + // // for each top-vol branch in the block + // // we get their upper end and lower end local ids + // IdArrayType topVolLocalBranchUpperEnd; + // IdArrayType topVolLocalBranchLowerEnd; + // vtkm::cont::ArrayHandle topVolIsParent; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // upperEndLocalIds, b->TopVolBranchInfoActualIndex, topVolLocalBranchUpperEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // lowerEndLocalIds, b->TopVolBranchInfoActualIndex, topVolLocalBranchLowerEnd); + // vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex>( + // b->IsParentBranch, b->TopVolBranchInfoActualIndex, topVolIsParent); + // auto topVolIsParentPortal = topVolIsParent.ReadPortal(); + // + // // we compute the superarc of the branch within the block + // // around the given isovalue + // IdArrayType branchIsoSuperarcs; + // branchIsoSuperarcs.Allocate(nTopVolBranchWithinBlock); + // + // auto branchIsoSuperarcWorklet = + // vtkm::worklet::scalar_topology::select_top_volume_contours::GetSuperarcByIsoValueWorklet( + // totalNumPoints); + // invoke(branchIsoSuperarcWorklet, + // topVolLocalBranchUpperEnd, + // topVolLocalBranchLowerEnd, + // isoValues, + // topVolBranchSaddleEpsilons, + // branchIsoSuperarcs, + // findSuperarcForNode); + // auto branchIsoSuperarcsPortal = branchIsoSuperarcs.ReadPortal(); + // + //#ifdef DEBUG_PRINT + // std::stringstream branchStream; + // + // branchStream << "Debug for branch info, Block " << b->LocalBlockNo << std::endl; + // vtkm::worklet::contourtree_augmented::PrintHeader(sortedBranchGRId.GetNumberOfValues(), + // branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Raw Branch GR", b->BranchRootGRId, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Raw Upper End", upperLocalEndIds, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Raw Lower End", lowerLocalEndIds, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Sorted Branch GR", sortedBranchGRId, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Sorted Branch Id", sortedBranchOrder, -1, branchStream); + // + // vtkm::worklet::contourtree_augmented::PrintHeader(nIsoValues, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Top Branch GR", b->TopVolumeBranchRootGRId, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Top Branch Stencil", topVolBranchWithinBlockStencil, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Top Branch Idx", topVolBranchInfoIdx, -1, branchStream); + // + // vtkm::worklet::contourtree_augmented::PrintHeader(nTopVolBranchKnownByBlock, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Sorted Upper End", topVolBranchUpperLocalEnd, -1, branchStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Sorted Lower End", topVolBranchLowerLocalEnd, -1, branchStream); + // + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, branchStream.str()); + //#endif + // + // const vtkm::Id nExtraMaximaBranch = b->ExtraMaximaBranchLowerEnd.GetNumberOfValues(); + // const vtkm::Id nExtraMinimaBranch = b->ExtraMinimaBranchLowerEnd.GetNumberOfValues(); + // InArrayHandleType extraMaximaBranchIsoValue; + // InArrayHandleType extraMinimaBranchIsoValue; + // + // IdArrayType extraMaximaBranchSuperarcs; + // IdArrayType extraMinimaBranchSuperarcs; + // extraMaximaBranchSuperarcs.Allocate(nExtraMaximaBranch); + // extraMinimaBranchSuperarcs.Allocate(nExtraMinimaBranch); + // + // if (nExtraMaximaBranch) + // { + // extraMaximaBranchIsoValue = b->ExtraMaximaBranchIsoValue.AsArrayHandle(); + // + // invoke(branchIsoSuperarcWorklet, + // b->ExtraMaximaBranchUpperEnd, + // b->ExtraMaximaBranchLowerEnd, + // extraMaximaBranchIsoValue, + // vtkm::cont::ArrayHandleConstant(1, nExtraMaximaBranch), + // extraMaximaBranchSuperarcs, + // findSuperarcForNode); + // + //#ifdef DEBUG_PRINT + // std::stringstream extraMaxStream; + // extraMaxStream << "Debug for Extra Maxima Branch, Block " << b->LocalBlockNo << std::endl; + // vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMaximaBranch, extraMaxStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Upper End", b->ExtraMaximaBranchUpperEnd, -1, extraMaxStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Lower End", b->ExtraMaximaBranchLowerEnd, -1, extraMaxStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "Max Branch IsoValue", extraMaximaBranchIsoValue, -1, extraMaxStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Superarc", extraMaximaBranchSuperarcs, -1, extraMaxStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Order", b->ExtraMaximaBranchOrder, -1, extraMaxStream); + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMaxStream.str()); + //#endif // DEBUG_PRINT + // } + // + // if (nExtraMinimaBranch) + // { + // extraMinimaBranchIsoValue = b->ExtraMinimaBranchIsoValue.AsArrayHandle(); + // + // invoke(branchIsoSuperarcWorklet, + // b->ExtraMinimaBranchUpperEnd, + // b->ExtraMinimaBranchLowerEnd, + // extraMinimaBranchIsoValue, + // vtkm::cont::ArrayHandleConstant(-1, nExtraMinimaBranch), + // extraMinimaBranchSuperarcs, + // findSuperarcForNode); + // + //#ifdef DEBUG_PRINT + // std::stringstream extraMinStream; + // extraMaxStream << "Debug for Extra Maxima Branch, Block " << b->LocalBlockNo << std::endl; + // vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMinimaBranch, extraMinStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Upper End", b->ExtraMinimaBranchUpperEnd, -1, extraMinStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Lower End", b->ExtraMinimaBranchLowerEnd, -1, extraMinStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "Max Branch IsoValue", extraMinimaBranchIsoValue, -1, extraMinStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Superarc", extraMinimaBranchSuperarcs, -1, extraMinStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Max Branch Order", b->ExtraMinimaBranchOrder, -1, extraMinStream); + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMinStream.str()); + //#endif // DEBUG_PRINT + // } + // + // auto extraMaximaBranchSuperarcPortal = extraMaximaBranchSuperarcs.ReadPortal(); + // auto extraMinimaBranchSuperarcPortal = extraMinimaBranchSuperarcs.ReadPortal(); + // auto extraMaximaBranchIsoValuePortal = extraMaximaBranchIsoValue.ReadPortal(); + // auto extraMinimaBranchIsoValuePortal = extraMinimaBranchIsoValue.ReadPortal(); + // auto extraMaximaBranchOrderPortal = b->ExtraMaximaBranchOrder.ReadPortal(); + // auto extraMinimaBranchOrderPortal = b->ExtraMinimaBranchOrder.ReadPortal(); + // + // const vtkm::Id nContours = nTopVolBranchWithinBlock + nExtraMaximaBranch + nExtraMinimaBranch; + // b->IsosurfaceEdgesOffset.AllocateAndFill(nContours, 0); + // b->IsosurfaceEdgesLabels.AllocateAndFill(nContours, 0); + // b->IsosurfaceEdgesOrders.AllocateAndFill(nContours, 0); + // InArrayHandleType isosurfaceIsoValue; + // isosurfaceIsoValue.AllocateAndFill(nContours, ValueType(0)); + // auto edgeOffsetWritePortal = b->IsosurfaceEdgesOffset.WritePortal(); + // auto edgeLabelWritePortal = b->IsosurfaceEdgesLabels.WritePortal(); + // auto edgeOrderWritePortal = b->IsosurfaceEdgesOrders.WritePortal(); + // auto isosurfaceValuePortal = isosurfaceIsoValue.WritePortal(); + // + // vtkm::Id nContourCandidateMeshes = 0; + // // NOTE: nContours denotes the number of isosurfaces for visualization. + // // The number is usually small, so linear loop is not too costly. + // // NOTE update 06/16/2024: We always need the isovalue of the contour. + // // As a result, we iterate through nContours (=O(k)). + // // This may be parallelizable in future work + // for (vtkm::Id branchIdx = 0; branchIdx < nContours; branchIdx++) + // { + // ValueType isoValue; + // vtkm::Id currBranchSaddleEpsilon, branchSuperarc, branchOrder, branchLabel = 0; + // + // using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_SADDLE; + // using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_COVER; + // using vtkm::worklet::scalar_topology::select_top_volume_contours::MAXIMA_CONTOUR; + // + // if (branchIdx < nTopVolBranchWithinBlock) + // { + // isoValue = isoValuePortal.Get(branchIdx); + // currBranchSaddleEpsilon = topVolBranchSaddleEpsilonPortal.Get(branchIdx); + // branchSuperarc = branchIsoSuperarcsPortal.Get(branchIdx); + // branchOrder = topVolBranchWithinBlockIdPortal.Get(branchIdx) + 1; + // branchLabel |= BRANCH_SADDLE; + // branchLabel |= topVolIsParentPortal.Get(branchIdx) ? BRANCH_COVER : 0; + // branchLabel |= currBranchSaddleEpsilon > 0 ? MAXIMA_CONTOUR : 0; + // } + // else if (branchIdx < nTopVolBranchWithinBlock + nExtraMaximaBranch) + // { + // const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock; + // isoValue = extraMaximaBranchIsoValuePortal.Get(idx); + // currBranchSaddleEpsilon = 1; + // branchSuperarc = extraMaximaBranchSuperarcPortal.Get(idx); + // branchOrder = extraMaximaBranchOrderPortal.Get(idx); + // branchLabel |= MAXIMA_CONTOUR; + // } + // else + // { + // const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock - nExtraMaximaBranch; + // VTKM_ASSERT(idx < nExtraMinimaBranch); + // isoValue = extraMinimaBranchIsoValuePortal.Get(idx); + // currBranchSaddleEpsilon = -1; + // branchSuperarc = extraMinimaBranchSuperarcPortal.Get(idx); + // branchOrder = extraMinimaBranchOrderPortal.Get(idx); + // } + // + // edgeOffsetWritePortal.Set(branchIdx, b->IsosurfaceEdgesFrom.GetNumberOfValues()); + // edgeLabelWritePortal.Set(branchIdx, branchLabel); + // edgeOrderWritePortal.Set(branchIdx, branchOrder); + // isosurfaceValuePortal.Set(branchIdx, isoValue); + // + // if (vtkm::worklet::contourtree_augmented::NoSuchElement(branchSuperarc)) + // { + // continue; + // } + // + // // Note: by concept, there is no 3D cell if pointDimensions[2] <= 1 + // vtkm::Id nCells = globalPointDimensions[2] <= 1 + // ? (pointDimensions[0] - 1) * (pointDimensions[1] - 1) + // : (pointDimensions[0] - 1) * (pointDimensions[1] - 1) * (pointDimensions[2] - 1); + // + // // we get the polarity cases of cells + // // we use lookup tables to for fast processing + // IdArrayType vertexOffset = globalPointDimensions[2] <= 1 + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset2d + // : vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset3d; + // IdArrayType edgeTable = globalPointDimensions[2] <= 1 + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTable2d + // : (isMarchingCubes + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableMC3d + // : vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableLT3d); + // IdArrayType numBoundTable = globalPointDimensions[2] <= 1 + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::numLinesTable2d + // : (isMarchingCubes + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableMC3d + // : vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableLT3d); + // IdArrayType boundaryTable = globalPointDimensions[2] <= 1 + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::lineTable2d + // : (isMarchingCubes + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::triTableMC3d + // : vtkm::worklet::scalar_topology::select_top_volume_contours::triTableLT3d); + // IdArrayType labelEdgeTable = isMarchingCubes + // ? vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableMC3d + // : vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableLT3d; + // + // IdArrayType caseCells; + // + // caseCells.Allocate(nCells); + // IdArrayType numEdgesInCell; + // numEdgesInCell.Allocate(nCells); + // auto caseCellsWorklet = + // vtkm::worklet::scalar_topology::select_top_volume_contours::GetCellCasesWorklet< + // ValueType>(pointDimensions, currBranchSaddleEpsilon, isoValue); + // + // invoke(caseCellsWorklet, + // vtkm::cont::ArrayHandleIndex(nCells), + // sortedDataValuesWithinBlock, + // vertexOffset, + // caseCells); + // + // // we compute the number of edges for each cell + // // to initialize the array size of edges + // IdArrayType numBoundariesInCell; + // vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + // numBoundTable, caseCells, numBoundariesInCell); + // + // vtkm::cont::Algorithm::Transform(numBoundariesInCell, + // globalPointDimensions[2] <= 1 + // ? vtkm::cont::ArrayHandleConstant(1, nCells) + // : vtkm::cont::ArrayHandleConstant(3, nCells), + // numEdgesInCell, + // vtkm::Multiply()); + // + // // use prefix sum to get the offset of the starting edge in each cell/cube + // vtkm::Id nEdges = + // vtkm::cont::Algorithm::Reduce(numEdgesInCell, vtkm::Id(0)); + // nContourCandidateMeshes += globalPointDimensions[2] <= 1 ? nEdges : nEdges / 3; + // IdArrayType edgesOffset; + // vtkm::cont::Algorithm::ScanExclusive(numEdgesInCell, edgesOffset); + // + // vtkm::cont::ArrayHandle isosurfaceEdgesFrom; + // vtkm::cont::ArrayHandle isosurfaceEdgesTo; + // // IdArrayType isosurfaceEdgesLabels; + // IdArrayType isValidEdges; + // isosurfaceEdgesFrom.Allocate(nEdges); + // isosurfaceEdgesTo.Allocate(nEdges); + // // isosurfaceEdgesLabels.Allocate(nEdges); + // isValidEdges.Allocate(nEdges); + // + // // draw isosurface + // auto getEdgesInCellWorklet = + // vtkm::worklet::scalar_topology::select_top_volume_contours::GetEdgesInCellWorklet< + // ValueType>(pointDimensions, + // globalPointIndexStart, + // isoValue, + // branchSuperarc, + // currBranchSaddleEpsilon, + // totalNumPoints, + // isMarchingCubes); + // + // invoke(getEdgesInCellWorklet, + // edgesOffset, + // caseCells, + // sortedLocalNodeInfoIdsWithinBlock, + // sortedDataValuesWithinBlock, + // vertexOffset, + // edgeTable, + // numBoundTable, + // boundaryTable, + // labelEdgeTable, + // isosurfaceEdgesFrom, + // isosurfaceEdgesTo, + // isValidEdges, + // findSuperarcForNode); + // + // // isValidEdges: stencil about whether the edge is on the desired superarc + // vtkm::cont::ArrayHandle validEdgesFrom; + // vtkm::cont::ArrayHandle validEdgesTo; + // + // // we remove edges that are not on the desired superarc + // vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesFrom, isValidEdges, validEdgesFrom); + // vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesTo, isValidEdges, validEdgesTo); + // + // // append edges into the result array + // vtkm::Id nValidEdges = validEdgesFrom.GetNumberOfValues(); + // vtkm::Id nExistEdges = branchIdx == 0 ? 0 : b->IsosurfaceEdgesFrom.GetNumberOfValues(); + // + // b->IsosurfaceEdgesFrom.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); + // b->IsosurfaceEdgesTo.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); + // + //#ifdef DEBUG_PRINT + // std::stringstream edgeInfoStream; + // contourStream << "Debug for Contour Info, Block " << b->LocalBlockNo << std::endl; + // vtkm::worklet::contourtree_augmented::PrintHeader(nContours, edgeInfoStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "edgeOffset", b->IsosurfaceEdgesOffset, -1, edgeInfoStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "edgeOrders", b->IsosurfaceEdgesOrders, -1, edgeInfoStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "edgeLabels", b->IsosurfaceEdgesLabels, -1, edgeInfoStream); + // + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, edgeInfoStream.str()); + //#endif // DEBUG_PRINT + // + // vtkm::cont::Algorithm::CopySubRange( + // validEdgesFrom, 0, nValidEdges, b->IsosurfaceEdgesFrom, nExistEdges); + // vtkm::cont::Algorithm::CopySubRange( + // validEdgesTo, 0, nValidEdges, b->IsosurfaceEdgesTo, nExistEdges); + // + //#ifdef DEBUG_PRINT + // std::stringstream contourStream; + // contourStream << "Debug for Contour, Block " << b->LocalBlockNo << std::endl; + // contourStream << "Branch Superarc = " << branchSuperarc << ", isoValue = " << isoValue + // << std::endl; + // + // vtkm::worklet::contourtree_augmented::PrintHeader(nCells, contourStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Case of Cells", caseCells, -1, contourStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "# of Edges", numEdgesInCell, -1, contourStream); + // + // vtkm::worklet::contourtree_augmented::PrintHeader(nEdges, contourStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "Edge Stencil", isValidEdges, -1, contourStream); + // + // vtkm::worklet::contourtree_augmented::PrintHeader(nValidEdges, contourStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "EdgesFrom", validEdgesFrom, -1, contourStream); + // vtkm::worklet::contourtree_augmented::PrintValues( + // "EdgesTo", validEdgesTo, -1, contourStream); + // vtkm::worklet::contourtree_augmented::PrintIndices( + // "EdgesLabels", validEdgesLabels, -1, contourStream); + // + // VTKM_LOG_S(vtkm::cont::LogLevel::Info, contourStream.str()); + //#endif // DEBUG_PRINT + // } + // b->IsosurfaceIsoValue = isosurfaceIsoValue; + // + // const vtkm::Id nMeshesOnBranches = globalPointDimensions[2] <= 1 + // ? b->IsosurfaceEdgesFrom.GetNumberOfValues() + // : b->IsosurfaceEdgesFrom.GetNumberOfValues() / 3; + // VTKM_LOG_S(TimingsLogLevel, + // std::endl + // << "----------- Draw Isosurface (block=" << b->LocalBlockNo << ")------------" + // << std::endl + // << " " << std::setw(60) << std::left + // << "Number of Contours: " << nContours + // << std::endl + // << " " << std::setw(60) << std::left + // << "Number of Isosurface Meshes: " << nContourCandidateMeshes + // << std::endl + // << " " << std::setw(60) << std::left + // << "Number of Meshes On Branches: " << nMeshesOnBranches << std::endl); + // }; + // b->TopVolumeBranchSaddleIsoValue + // .CastAndCallForTypes(resolveArray); }); timingsStream << " " << std::setw(60) << std::left << "Draw Contours By Branches" diff --git a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h index 41570d996..8b3798411 100644 --- a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h +++ b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h @@ -73,6 +73,10 @@ public: this->isMarchingCubes = marchingCubes; } + VTKM_CONT const vtkm::Id GetSavedBranches() { return this->nSavedBranches; } + VTKM_CONT const bool GetMarchingCubes() { return this->isMarchingCubes; } + VTKM_CONT const vtkm::cont::LogLevel GetTimingsLogLevel() { return this->TimingsLogLevel; } + private: VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet&) override; VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions( diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx index 230316260..9b549e870 100644 --- a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx @@ -52,8 +52,16 @@ #include #include +#include +#include +#include +#include #include +#include #include +#include +#include +#include #include #ifdef DEBUG_PRINT @@ -81,9 +89,8 @@ SelectTopVolumeContoursBlock::SelectTopVolumeContoursBlock(vtkm::Id localBlockNo { } -void SelectTopVolumeContoursBlock::SortBranchByVolume( - const vtkm::cont::DataSet& hierarchicalTreeDataSet, - const vtkm::Id totalVolume) +void SelectTopVolumeContoursBlock::SortBranchByVolume(const vtkm::cont::DataSet& bdDataSet, + const vtkm::Id totalVolume) { /// Pipeline to compute the branch volume /// 1. check both ends of the branch. If both leaves, then main branch, volume = totalVolume @@ -97,26 +104,26 @@ void SelectTopVolumeContoursBlock::SortBranchByVolume( vtkm::cont::ArrayHandle isLowerLeaf; vtkm::cont::ArrayHandle isUpperLeaf; - auto upperEndIntrinsicVolume = hierarchicalTreeDataSet.GetField("UpperEndIntrinsicVolume") + auto upperEndIntrinsicVolume = bdDataSet.GetField("UpperEndIntrinsicVolume") .GetData() .AsArrayHandle>(); - auto upperEndDependentVolume = hierarchicalTreeDataSet.GetField("UpperEndDependentVolume") + auto upperEndDependentVolume = bdDataSet.GetField("UpperEndDependentVolume") .GetData() .AsArrayHandle>(); - auto lowerEndIntrinsicVolume = hierarchicalTreeDataSet.GetField("LowerEndIntrinsicVolume") + auto lowerEndIntrinsicVolume = bdDataSet.GetField("LowerEndIntrinsicVolume") .GetData() .AsArrayHandle>(); - auto lowerEndDependentVolume = hierarchicalTreeDataSet.GetField("LowerEndDependentVolume") + auto lowerEndDependentVolume = bdDataSet.GetField("LowerEndDependentVolume") .GetData() .AsArrayHandle>(); - auto lowerEndSuperarcId = hierarchicalTreeDataSet.GetField("LowerEndSuperarcId") + auto lowerEndSuperarcId = bdDataSet.GetField("LowerEndSuperarcId") .GetData() .AsArrayHandle>(); - auto upperEndSuperarcId = hierarchicalTreeDataSet.GetField("UpperEndSuperarcId") + auto upperEndSuperarcId = bdDataSet.GetField("UpperEndSuperarcId") .GetData() .AsArrayHandle>(); - auto branchRoot = hierarchicalTreeDataSet.GetField("BranchRootByBranch") + auto branchRoot = bdDataSet.GetField("BranchRootByBranch") .GetData() .AsArrayHandle>(); @@ -145,8 +152,7 @@ void SelectTopVolumeContoursBlock::SortBranchByVolume( isLowerLeaf, isUpperLeaf); - vtkm::cont::UnknownArrayHandle upperEndValue = - hierarchicalTreeDataSet.GetField("UpperEndValue").GetData(); + vtkm::cont::UnknownArrayHandle upperEndValue = bdDataSet.GetField("UpperEndValue").GetData(); // Based on the direction info of the branch, store epsilon direction and isovalue of the saddle auto resolveArray = [&](const auto& inArray) { @@ -160,7 +166,7 @@ void SelectTopVolumeContoursBlock::SortBranchByVolume( vtkm::worklet::scalar_topology::select_top_volume_contours::UpdateInfoByBranchDirectionWorklet< ValueType> updateInfoWorklet; - auto lowerEndValue = hierarchicalTreeDataSet.GetField("LowerEndValue") + auto lowerEndValue = bdDataSet.GetField("LowerEndValue") .GetData() .AsArrayHandle>(); @@ -230,6 +236,1069 @@ void SelectTopVolumeContoursBlock::SortBranchByVolume( vtkm::cont::Algorithm::Copy(sortedBranches, this->SortedBranchByVolume); } +// Select the local top K branches by volume +void SelectTopVolumeContoursBlock::SelectLocalTopVolumeBranch(const vtkm::cont::DataSet& bdDataset, + const vtkm::Id nSavedBranches) +{ + using vtkm::worklet::contourtree_augmented::IdArrayType; + // copy the top volume branches into a smaller array + // we skip index 0 because it must be the main branch (which has the highest volume) + vtkm::Id nActualSavedBranches = + std::min(nSavedBranches, this->SortedBranchByVolume.GetNumberOfValues() - 1); + + vtkm::worklet::contourtree_augmented::IdArrayType topVolumeBranch; + vtkm::cont::Algorithm::CopySubRange( + this->SortedBranchByVolume, 1, nActualSavedBranches, topVolumeBranch); + + auto branchRootByBranch = bdDataset.GetField("BranchRootByBranch") + .GetData() + .AsArrayHandle>(); + + const vtkm::Id nBranches = branchRootByBranch.GetNumberOfValues(); + + auto branchRootGRId = bdDataset.GetField("BranchRootGRId") + .GetData() + .AsArrayHandle>(); + + auto upperEndGRId = bdDataset.GetField("UpperEndGlobalRegularIds") + .GetData() + .AsArrayHandle>(); + + auto lowerEndGRId = bdDataset.GetField("LowerEndGlobalRegularIds") + .GetData() + .AsArrayHandle>(); + + vtkm::cont::Algorithm::Copy(branchRootByBranch, this->BranchRootByBranch); + vtkm::cont::Algorithm::Copy(branchRootGRId, this->BranchRootGRId); + + // This seems weird, but we temporarily put the initialization of computing the branch decomposition tree here + this->IsParentBranch.AllocateAndFill(nBranches, false); + + // we permute all branch information to align with the order by volume + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchRootGRId, topVolumeBranch, this->TopVolumeBranchRootGRId); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + upperEndGRId, topVolumeBranch, this->TopVolumeBranchUpperEndGRId); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + lowerEndGRId, topVolumeBranch, this->TopVolumeBranchLowerEndGRId); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->BranchVolume, topVolumeBranch, this->TopVolumeBranchVolume); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->BranchSaddleEpsilon, topVolumeBranch, this->TopVolumeBranchSaddleEpsilon); + + auto resolveArray = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + InArrayHandleType topVolBranchSaddleIsoValue; + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + inArray, topVolumeBranch, topVolBranchSaddleIsoValue); + this->TopVolumeBranchSaddleIsoValue = topVolBranchSaddleIsoValue; + }; + + this->BranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); +} + + +void SelectTopVolumeContoursBlock::ComputeTopVolumeBranchHierarchy( + const vtkm::cont::DataSet& bdDataSet) +{ + using vtkm::worklet::contourtree_augmented::IdArrayType; + vtkm::cont::Invoker invoke; + + // we need upper/lower local ends and global ends for hierarchy of branches + auto upperLocalEndIds = bdDataSet.GetField("UpperEndLocalIds") + .GetData() + .AsArrayHandle>(); + auto lowerLocalEndIds = bdDataSet.GetField("LowerEndLocalIds") + .GetData() + .AsArrayHandle>(); + auto globalRegularIds = bdDataSet.GetField("RegularNodeGlobalIds") + .GetData() + .AsArrayHandle>(); + IdArrayType upperEndGRIds = + bdDataSet.GetField("UpperEndGlobalRegularIds").GetData().AsArrayHandle(); + IdArrayType lowerEndGRIds = + bdDataSet.GetField("LowerEndGlobalRegularIds").GetData().AsArrayHandle(); + + // let's check which top volume branches are known by the block + // We check the branchGRId of top volume branches to see whether there are matches within the block + const vtkm::Id nTopVolBranches = this->TopVolumeBranchLowerEndGRId.GetNumberOfValues(); + + // sortedBranchOrder: the branch order (in the ascending order of branch root) + // the high-level idea is to sort the branch root global regular ids + // and for each top-volume branch, we use binary search to get the original branch index + // if the top-volume branch does not exist in the block, it will be dropped out + IdArrayType sortedBranchGRId; + IdArrayType sortedBranchOrder; + vtkm::cont::Algorithm::Copy( + vtkm::cont::ArrayHandleIndex(this->BranchRootGRId.GetNumberOfValues()), sortedBranchOrder); + vtkm::cont::Algorithm::Copy(this->BranchRootGRId, sortedBranchGRId); + vtkm::cont::Algorithm::SortByKey(sortedBranchGRId, sortedBranchOrder); + + this->TopVolBranchKnownByBlockStencil.Allocate(nTopVolBranches); + this->TopVolBranchGROrder.Allocate(nTopVolBranches); + + // We reuse the IdxIfWithinBlockWorklet. + // This worklet searches for given values in a sorted array and returns the stencil & index if the value exists in the array. + // this->TopVolBranchGROrder: the order of the topVolBranch among all known branches if the branch is known by the block. + auto idxIfBranchWithinBlockWorklet = + vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); + invoke(idxIfBranchWithinBlockWorklet, + this->TopVolumeBranchRootGRId, + sortedBranchGRId, + this->TopVolBranchKnownByBlockStencil, + this->TopVolBranchGROrder); + + // Dropping out top-volume branches that are not known by the block. + + // the index of top-volume branches known by the block among all top-volume branches + IdArrayType topVolBranchKnownByBlockIndex; + vtkm::cont::ArrayHandleIndex topVolBranchesIndex(nTopVolBranches); + vtkm::cont::Algorithm::CopyIf( + topVolBranchesIndex, this->TopVolBranchKnownByBlockStencil, topVolBranchKnownByBlockIndex); + + const vtkm::Id nTopVolBranchKnownByBlock = topVolBranchKnownByBlockIndex.GetNumberOfValues(); + + // filtered this->TopVolBranchGROrder, by removing NO_SUCH_ELEMENT + IdArrayType topVolBranchFilteredGROrder; + + // this->TopVolBranchInfoActualIndex: the information index of the top-volume branch + vtkm::cont::Algorithm::CopyIf( + this->TopVolBranchGROrder, this->TopVolBranchKnownByBlockStencil, topVolBranchFilteredGROrder); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + sortedBranchOrder, topVolBranchFilteredGROrder, this->TopVolBranchInfoActualIndex); + + // filtered branch saddle epsilons, global lower/upper end GR ids, + IdArrayType topVolFilteredBranchSaddleEpsilon; + IdArrayType topVolFilteredLowerEndGRId; + IdArrayType topVolFilteredUpperEndGRId; + vtkm::cont::Algorithm::CopyIf(this->TopVolumeBranchSaddleEpsilon, + this->TopVolBranchKnownByBlockStencil, + topVolFilteredBranchSaddleEpsilon); + vtkm::cont::Algorithm::CopyIf(this->TopVolumeBranchUpperEndGRId, + this->TopVolBranchKnownByBlockStencil, + topVolFilteredUpperEndGRId); + vtkm::cont::Algorithm::CopyIf(this->TopVolumeBranchLowerEndGRId, + this->TopVolBranchKnownByBlockStencil, + topVolFilteredLowerEndGRId); + + // for each top-vol branch known by the block + // we get their upper end and lower end local ids + IdArrayType topVolBranchUpperLocalEnd; + IdArrayType topVolBranchLowerLocalEnd; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + upperLocalEndIds, this->TopVolBranchInfoActualIndex, topVolBranchUpperLocalEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + lowerLocalEndIds, this->TopVolBranchInfoActualIndex, topVolBranchLowerLocalEnd); + + IdArrayType topVolLowerLocalEndGRId; + IdArrayType topVolUpperLocalEndGRId; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + globalRegularIds, topVolBranchLowerLocalEnd, topVolLowerLocalEndGRId); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + globalRegularIds, topVolBranchUpperLocalEnd, topVolUpperLocalEndGRId); + + // Below is the code to compute the branch hierarchy of top-volume branches + // We need this information because we not only want to visualize the contour + // on top-volume branches, but also on their parent branches. + // Because we use volume as the metric, the parent branch of a top-volume branch + // is either a top-volume branch or the root branch (where both ends are leaf nodes) + vtkm::worklet::scalar_topology::select_top_volume_contours::BranchSaddleIsKnownWorklet + branchSaddleIsKnownWorklet; + // the branch saddle local ID if the saddle end is known by the block + IdArrayType branchSaddleIsKnown; + branchSaddleIsKnown.Allocate(nTopVolBranchKnownByBlock); + + invoke(branchSaddleIsKnownWorklet, // worklet + topVolFilteredLowerEndGRId, // input + topVolBranchLowerLocalEnd, // input + topVolLowerLocalEndGRId, // input + topVolFilteredUpperEndGRId, // input + topVolBranchUpperLocalEnd, // input + topVolUpperLocalEndGRId, // input + topVolFilteredBranchSaddleEpsilon, // input + branchSaddleIsKnown); // output + + // the order of top volume branches with parents known by the block + IdArrayType topVolChildBranch; + IdArrayType topVolChildBranchSaddle; + + vtkm::cont::Algorithm::CopyIf( + topVolBranchKnownByBlockIndex, + branchSaddleIsKnown, + topVolChildBranch, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); + vtkm::cont::Algorithm::CopyIf( + branchSaddleIsKnown, + branchSaddleIsKnown, + topVolChildBranchSaddle, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsNonNegative()); + + const vtkm::Id nChildBranch = topVolChildBranch.GetNumberOfValues(); + + // to compute the parent branch, we need to + // 1. for the branch saddle end, collect all superarcs involving it + // 2. get the branch information for selected superarcs + // 3. eliminate branch information for branches sharing the same saddle end + + auto superarcs = + bdDataSet.GetField("Superarcs").GetData().AsArrayHandle>(); + auto branchRoots = + bdDataSet.GetField("BranchRoots").GetData().AsArrayHandle>(); + VTKM_ASSERT(superarcs.GetNumberOfValues() == branchRoots.GetNumberOfValues()); + + // we sort all superarcs by target to allow binary search + IdArrayType superarcsByTarget; + vtkm::worklet::scalar_topology::select_top_volume_contours::SuperarcTargetComparator + superarcComparator(superarcs); + vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(superarcs.GetNumberOfValues()), + superarcsByTarget); + vtkm::cont::Algorithm::Sort(superarcsByTarget, superarcComparator); + + IdArrayType permutedSuperarcs; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + superarcs, superarcsByTarget, permutedSuperarcs); + + IdArrayType permutedBranchRoots; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchRoots, superarcsByTarget, permutedBranchRoots); + + // the branch root of the superarc of the branch saddle supernode + IdArrayType topVolChildBranchSaddleBranchRoot; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchRoots, topVolChildBranchSaddle, topVolChildBranchSaddleBranchRoot); + + // the GR Ids of the superarc of the branch saddle supernode + IdArrayType topVolChildBranchSaddleGRIds; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + globalRegularIds, topVolChildBranchSaddle, topVolChildBranchSaddleGRIds); + + // there is a debate to find all superarcs connect to a supernode + // strategy 1. iterate through saddles and parallelize over superarcs for search + // time complexity: O(nTopVolBranches) + // (nTopVolBranches usually <= 100, based on input parameter setting) + // + // strategy 2. parallelize over all saddles and use binary search to find superarcs + // time complexity: O(log_2(nSuperarcs)) (nSuperarcs can be considerably large) + // + // here, we choose strategy 2 for better extendability to high nTopVolBranches + // but in tests using nTopVolBranches <= 10, strategy 1 is generally faster + + // note: after getting the branch root superarc, we use binary search to get the branch order + // because BranchRootByBranch is sorted by branch root (superarc) id + +#ifdef DEBUG_PRINT + std::stringstream parentBranchStream; + + parentBranchStream << "Debug for Parent Branch, Block " << this->LocalBlockNo << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nChildBranch, parentBranchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Child Branch Saddle", topVolChildBranchSaddle, -1, parentBranchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Child Saddle Root", topVolChildBranchSaddleBranchRoot, -1, parentBranchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Child Saddle GR Id", topVolChildBranchSaddleGRIds, -1, parentBranchStream); + + vtkm::worklet::contourtree_augmented::PrintHeader(superarcs.GetNumberOfValues(), + parentBranchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Permuted Superarcs", permutedSuperarcs, -1, parentBranchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Permuted Branch roots", permutedBranchRoots, -1, parentBranchStream); + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, parentBranchStream.str()); +#endif // DEBUG_PRINT + + // the corresponding parent branch of child branches + IdArrayType topVolChildBranchParent; + topVolChildBranchParent.AllocateAndFill(nChildBranch, + vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT); + vtkm::worklet::scalar_topology::select_top_volume_contours::GetParentBranchWorklet + getParentBranchWorklet; + invoke(getParentBranchWorklet, + topVolChildBranchSaddle, + topVolChildBranchSaddleBranchRoot, + topVolChildBranchSaddleGRIds, + permutedSuperarcs, + permutedBranchRoots, + this->BranchRootByBranch, + upperEndGRIds, + lowerEndGRIds, + topVolChildBranchParent); + + this->TopVolumeBranchParent.AllocateAndFill( + nTopVolBranches, vtkm::Id(vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT)); + + vtkm::worklet::scalar_topology::select_top_volume_contours::AssignValueByIndex assignParentBranch; + // for each top volume branch, assign the parent branch info id in the block + invoke( + assignParentBranch, topVolChildBranch, topVolChildBranchParent, this->TopVolumeBranchParent); + // for each branch, assign true if it is a parent branch + invoke(assignParentBranch, + topVolChildBranchParent, + vtkm::cont::ArrayHandleConstant(true, nChildBranch), + this->IsParentBranch); + + // sort all top-volume branches based on + // 1. parent branch info id: this->TopVolumeBranchParent + // 2. saddle-end value: this->TopVolumeBranchSaddleIsovalue + // 3. branch root global regular id (anything that can break tie) + auto resolveBranchParent = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + using ValueType = typename InArrayHandleType::ValueType; + + vtkm::worklet::scalar_topology::select_top_volume_contours::BranchParentComparator + parentComparator(this->TopVolumeBranchParent, inArray, this->TopVolumeBranchRootGRId); + + // sort index for all top volume branches + IdArrayType topVolSortForOuterSaddleIdx; + vtkm::cont::Algorithm::Copy(topVolBranchesIndex, topVolSortForOuterSaddleIdx); + vtkm::cont::Algorithm::Sort(topVolSortForOuterSaddleIdx, parentComparator); + + IdArrayType parentPermutation; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->TopVolumeBranchParent, topVolSortForOuterSaddleIdx, parentPermutation); + + // warning: when parent is NO_SUCH_ELEMENT, parentSaddleEps obtains 0 + // However, the corresponding element will be discarded in collecting outer saddles + IdArrayType parentSaddleEpsPermutation; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->BranchSaddleEpsilon, parentPermutation, parentSaddleEpsPermutation); + + // Some branches have parent=NO_SUCH_ELEMENT (no parent) + // we collect the isovalue of the first and/or the last branches for each parent branch + // we collect the first if branchSaddleEpsilon(parent) < 0 + // or the last if branchSaddleEpsilon(parent) > 0 + // or both if branchSaddleEpsilon(parent) == 0 + IdArrayType IsOuterSaddle; + IsOuterSaddle.Allocate(nTopVolBranches); + vtkm::worklet::scalar_topology::select_top_volume_contours::CollectOuterSaddle + collectOuterSaddleWorklet; + invoke(collectOuterSaddleWorklet, parentSaddleEpsPermutation, parentPermutation, IsOuterSaddle); + + // after sorting by index back + // each top volume branch know whether it is the outer saddle of its parent + vtkm::cont::Algorithm::SortByKey(topVolSortForOuterSaddleIdx, IsOuterSaddle); + + // collect branches that need contours on extra minima/maxima + // we store the information of the parent branches (on both directions) + IdArrayType extraMaximaParentBranch; + IdArrayType extraMinimaParentBranch; + IdArrayType extraMaximaParentBranchRootGRId; + IdArrayType extraMinimaParentBranchRootGRId; + + IdArrayType allBranchGRIdByVolume; + IdArrayType branchGRIdByVolumeIdx; + + // we need global branch order including the root branch + // this information should be consistent globally + allBranchGRIdByVolume.Allocate(nTopVolBranches + 1); + vtkm::cont::Algorithm::CopySubRange( + this->TopVolumeBranchRootGRId, 0, nTopVolBranches, allBranchGRIdByVolume, 1); + auto topBranchGRIdWritePortal = allBranchGRIdByVolume.WritePortal(); + auto sortedBranchByVolPortal = this->SortedBranchByVolume.ReadPortal(); + auto branchGRIdReadPortal = this->BranchRootGRId.ReadPortal(); + topBranchGRIdWritePortal.Set(0, branchGRIdReadPortal.Get(sortedBranchByVolPortal.Get(0))); + + vtkm::cont::Algorithm::Copy( + vtkm::cont::ArrayHandleIndex(allBranchGRIdByVolume.GetNumberOfValues()), + branchGRIdByVolumeIdx); + vtkm::cont::Algorithm::SortByKey(allBranchGRIdByVolume, branchGRIdByVolumeIdx); + + vtkm::cont::Algorithm::CopyIf( + this->TopVolumeBranchParent, + IsOuterSaddle, + extraMaximaParentBranch, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); + vtkm::cont::Algorithm::CopyIf( + this->TopVolumeBranchParent, + IsOuterSaddle, + extraMinimaParentBranch, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); + + if (extraMaximaParentBranch.GetNumberOfValues()) + { + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + upperLocalEndIds, extraMaximaParentBranch, this->ExtraMaximaBranchUpperEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + lowerLocalEndIds, extraMaximaParentBranch, this->ExtraMaximaBranchLowerEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->BranchRootGRId, extraMaximaParentBranch, extraMaximaParentBranchRootGRId); + + InArrayHandleType extraMaximaBranchIsoValue; + vtkm::cont::Algorithm::CopyIf( + this->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), + IsOuterSaddle, + extraMaximaBranchIsoValue, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMaxima()); + this->ExtraMaximaBranchIsoValue = extraMaximaBranchIsoValue; + + // a worklet to binary search a number in a sorted array and return the index + vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet + getParentBranchOrder; + IdArrayType permutedExtraMaximaBranchOrder; + permutedExtraMaximaBranchOrder.Allocate(extraMaximaParentBranchRootGRId.GetNumberOfValues()); + vtkm::cont::ArrayHandleDiscard discard; + discard.Allocate(extraMaximaParentBranchRootGRId.GetNumberOfValues()); + invoke(getParentBranchOrder, + extraMaximaParentBranchRootGRId, + allBranchGRIdByVolume, + discard, + permutedExtraMaximaBranchOrder); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchGRIdByVolumeIdx, permutedExtraMaximaBranchOrder, this->ExtraMaximaBranchOrder); + } + + if (extraMinimaParentBranch.GetNumberOfValues()) + { + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + upperLocalEndIds, extraMinimaParentBranch, this->ExtraMinimaBranchUpperEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + lowerLocalEndIds, extraMinimaParentBranch, this->ExtraMinimaBranchLowerEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + this->BranchRootGRId, extraMinimaParentBranch, extraMinimaParentBranchRootGRId); + + InArrayHandleType extraMinimaBranchIsoValue; + vtkm::cont::Algorithm::CopyIf( + this->TopVolumeBranchSaddleIsoValue.AsArrayHandle(), + IsOuterSaddle, + extraMinimaBranchIsoValue, + vtkm::worklet::scalar_topology::select_top_volume_contours::IsExtraMinima()); + this->ExtraMinimaBranchIsoValue = extraMinimaBranchIsoValue; + + vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet + getParentBranchOrder; + IdArrayType permutedExtraMinimaBranchOrder; + permutedExtraMinimaBranchOrder.Allocate(extraMinimaParentBranchRootGRId.GetNumberOfValues()); + vtkm::cont::ArrayHandleDiscard discard; + discard.Allocate(extraMinimaParentBranchRootGRId.GetNumberOfValues()); + invoke(getParentBranchOrder, + extraMinimaParentBranchRootGRId, + allBranchGRIdByVolume, + discard, + permutedExtraMinimaBranchOrder); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchGRIdByVolumeIdx, permutedExtraMinimaBranchOrder, this->ExtraMinimaBranchOrder); + } + }; + this->TopVolumeBranchSaddleIsoValue + .CastAndCallForTypes( + resolveBranchParent); +} + + +void SelectTopVolumeContoursBlock::ExtractIsosurfaceOnSelectedBranch( + const vtkm::cont::DataSet& bdDataSet, + const bool isMarchingCubes, + const vtkm::cont::LogLevel timingsLogLevel) +{ + using vtkm::worklet::contourtree_augmented::IdArrayType; + + // if no branch to extract from + if (this->TopVolumeBranchRootGRId.GetNumberOfValues() < 1) + return; + + // Let's create a mesh dataset to extract all the contours first + + // branch root global regular ID + // size: nBranches + // usage: identifier of the branch + vtkm::cont::Algorithm::Copy( + bdDataSet.GetField("BranchRootGRId").GetData().AsArrayHandle(), + this->BranchRootGRId); + + // branch local upper end and lower end + // size: nBranches + // usage: search for the superarc of an arbitrary point (not necessarily on grid) + auto upperEndLocalIds = + bdDataSet.GetField("UpperEndLocalIds").GetData().AsArrayHandle(); + auto lowerEndLocalIds = + bdDataSet.GetField("LowerEndLocalIds").GetData().AsArrayHandle(); + + // global regular ids + auto globalRegularIds = + bdDataSet.GetField("RegularNodeGlobalIds").GetData().AsArrayHandle(); + + vtkm::Id3 globalPointDimensions; + vtkm::Id3 pointDimensions, globalPointIndexStart; + + bdDataSet.GetCellSet().CastAndCallForTypes( + vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + pointDimensions, + globalPointDimensions, + globalPointIndexStart); + +#ifdef DEBUG_PRINT + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Block size info"); + std::stringstream rs; + { + rs << "globalPointDimensions: " << globalPointDimensions << std::endl; + rs << "pointDimensions: " << pointDimensions << std::endl; + rs << "globalPointIndexStart: " << globalPointIndexStart << std::endl; + rs << "globalRegularIDs: " << globalRegularIds.GetNumberOfValues() << std::endl; + // ds.PrintSummary(rs); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); + } +#endif + + // Tool to relabel local mesh ids to global ids + auto localToGlobalIdRelabeler = vtkm::worklet::contourtree_augmented::mesh_dem::IdRelabeler( + globalPointIndexStart, pointDimensions, globalPointDimensions); + IdArrayType globalIdsByMesh; + + // Note: the cell set is different from the mesh structure + // Here, we assume that the cell set is structured grid + // A more general way to do this is to use CellSet().GetCellPointIds(i) to extract all the local ids and keep unique ones + // local ids in the mesh + IdArrayType localIdsByMesh; + vtkm::cont::Algorithm::Copy( + vtkm::cont::ArrayHandleIndex(pointDimensions[0] * pointDimensions[1] * pointDimensions[2]), + localIdsByMesh); + // then, we transform the local ids to global ids + auto localTransformToGlobalId = + vtkm::cont::make_ArrayHandleTransform(localIdsByMesh, localToGlobalIdRelabeler); + vtkm::cont::ArrayCopyDevice(localTransformToGlobalId, globalIdsByMesh); + + // detect whether the element in globalRegularIds are in the block + // globalIdsDiscard is just a filler for the worklet format. We do not use it. + // The last slot for the worklet is useful in a later step. + // Here we just reuse the worklet + IdArrayType globalIdsWithinBlockStencil; + vtkm::cont::ArrayHandleDiscard globalIdsDiscard; + globalIdsWithinBlockStencil.Allocate(globalRegularIds.GetNumberOfValues()); + globalIdsDiscard.Allocate(globalRegularIds.GetNumberOfValues()); + + vtkm::cont::Invoker invoke; + // stencil is 1 if the global regular id is within the block, 0 otherwise + auto idxIfWithinBlock = + vtkm::worklet::scalar_topology::select_top_volume_contours::IdxIfWithinBlockWorklet(); + invoke(idxIfWithinBlock, + globalRegularIds, + globalIdsByMesh, + globalIdsWithinBlockStencil, + globalIdsDiscard); + + auto resolveArray = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + using ValueType = typename InArrayHandleType::ValueType; + + // we need to sort all values based on the global ids + // and remove values of points that do not belong to the local block + auto dataValues = bdDataSet.GetField("DataValues").GetData().AsArrayHandle(); + + IdArrayType globalIdsWithinBlock; + IdArrayType localIdsWithinBlock; + InArrayHandleType dataValuesWithinBlock; + + // filter global regular ids, array ids, and data values + vtkm::cont::Algorithm::CopyIf( + globalRegularIds, globalIdsWithinBlockStencil, globalIdsWithinBlock); + vtkm::cont::Algorithm::CopyIf( + vtkm::cont::ArrayHandleIndex(globalRegularIds.GetNumberOfValues()), + globalIdsWithinBlockStencil, + localIdsWithinBlock); + vtkm::cont::Algorithm::CopyIf(dataValues, globalIdsWithinBlockStencil, dataValuesWithinBlock); + + // sorted index based on global regular ids + IdArrayType sortedGlobalIds; + vtkm::cont::Algorithm::Copy(vtkm::cont::ArrayHandleIndex(globalIdsByMesh.GetNumberOfValues()), + sortedGlobalIds); + vtkm::cont::Algorithm::SortByKey(globalIdsWithinBlock, sortedGlobalIds); + + // globalIdsWithinBlock (sorted) and globalIdsByMesh should be identical. + // computing globalIdsWithinBlock ensures the input data is correct + bool identical = + globalIdsWithinBlock.GetNumberOfValues() == globalIdsByMesh.GetNumberOfValues(); + if (identical) + { + vtkm::cont::ArrayHandle globalIdsIdentical; + vtkm::cont::Algorithm::Transform( + globalIdsWithinBlock, globalIdsByMesh, globalIdsIdentical, vtkm::Equal()); + identical = vtkm::cont::Algorithm::Reduce(globalIdsIdentical, true, vtkm::LogicalAnd()); + } + if (!identical) + { + vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsByMesh.GetNumberOfValues()); + vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsByMesh", globalIdsByMesh); + vtkm::worklet::contourtree_augmented::PrintHeader(globalIdsWithinBlock.GetNumberOfValues()); + vtkm::worklet::contourtree_augmented::PrintIndices("globalIdsWithinBlock", + globalIdsWithinBlock); + } + VTKM_ASSERT(identical); + + // filtered and sorted local node info ids + // i.e. index of global regular ids, data values, and superparents + // Note: This is not local mesh id! Make sure to distinguish them + IdArrayType sortedLocalNodeInfoIdsWithinBlock; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + localIdsWithinBlock, sortedGlobalIds, sortedLocalNodeInfoIdsWithinBlock); + + // sorted data values + InArrayHandleType sortedDataValuesWithinBlock; + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + dataValuesWithinBlock, sortedGlobalIds, sortedDataValuesWithinBlock); + + // create an execution object to find the superarc for an arbitrary point within the mesh + // all information below are required to initialize the execution object + auto superparents = bdDataSet.GetField("Superparents").GetData().AsArrayHandle(); + auto supernodes = bdDataSet.GetField("Supernodes").GetData().AsArrayHandle(); + auto superarcs = bdDataSet.GetField("Superarcs").GetData().AsArrayHandle(); + auto superchildren = bdDataSet.GetField("Superchildren").GetData().AsArrayHandle(); + auto whichRound = bdDataSet.GetField("WhichRound").GetData().AsArrayHandle(); + auto whichIteration = + bdDataSet.GetField("WhichIteration").GetData().AsArrayHandle(); + auto hyperparents = bdDataSet.GetField("Hyperparents").GetData().AsArrayHandle(); + auto hypernodes = bdDataSet.GetField("Hypernodes").GetData().AsArrayHandle(); + auto hyperarcs = bdDataSet.GetField("Hyperarcs").GetData().AsArrayHandle(); + + // filtered + sorted superparents of nodes + IdArrayType superparentsWithinBlock; + vtkm::cont::Algorithm::CopyIf( + superparents, globalIdsWithinBlockStencil, superparentsWithinBlock); + IdArrayType sortedSuperparentsWithinBlock; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + superparentsWithinBlock, sortedGlobalIds, sortedSuperparentsWithinBlock); + + // initialize the exec object + // note: terms should include the contour tree as much as possible + // should pass the full arrays to the object instead of the filtered ones + auto findSuperarcForNode = + vtkm::worklet::contourtree_distributed::FindSuperArcForUnknownNode( + superparents, + supernodes, + superarcs, + superchildren, + whichRound, + whichIteration, + hyperparents, + hypernodes, + hyperarcs, + globalRegularIds, + dataValues); + + // let's check which branches are known by the block + // We check the branchGRId of top volume branches to see whether there are matches within the block + vtkm::Id nIsoValues = inArray.GetNumberOfValues(); + vtkm::Id totalNumPoints = + globalPointDimensions[0] * globalPointDimensions[1] * globalPointDimensions[2]; + + // dropping out top-volume branches that are not known by the block + // index of top-volume branches within the block among all top-volume branches + + IdArrayType topVolBranchWithinBlockId; + vtkm::cont::Algorithm::CopyIf(vtkm::cont::ArrayHandleIndex(nIsoValues), + this->TopVolBranchKnownByBlockStencil, + topVolBranchWithinBlockId); + auto topVolBranchWithinBlockIdPortal = topVolBranchWithinBlockId.ReadPortal(); + + vtkm::Id nTopVolBranchWithinBlock = topVolBranchWithinBlockId.GetNumberOfValues(); + + // filtered branch saddle values + InArrayHandleType isoValues; + vtkm::cont::Algorithm::CopyIf(inArray, this->TopVolBranchKnownByBlockStencil, isoValues); + auto isoValuePortal = isoValues.ReadPortal(); + + // filtered branch saddle epsilons + IdArrayType topVolBranchSaddleEpsilons; + vtkm::cont::Algorithm::CopyIf(this->TopVolumeBranchSaddleEpsilon, + this->TopVolBranchKnownByBlockStencil, + topVolBranchSaddleEpsilons); + auto topVolBranchSaddleEpsilonPortal = topVolBranchSaddleEpsilons.ReadPortal(); + + // for each top-vol branch in the block + // we get their upper end and lower end local ids + IdArrayType topVolLocalBranchUpperEnd; + IdArrayType topVolLocalBranchLowerEnd; + vtkm::cont::ArrayHandle topVolIsParent; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + upperEndLocalIds, this->TopVolBranchInfoActualIndex, topVolLocalBranchUpperEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + lowerEndLocalIds, this->TopVolBranchInfoActualIndex, topVolLocalBranchLowerEnd); + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex>( + this->IsParentBranch, this->TopVolBranchInfoActualIndex, topVolIsParent); + auto topVolIsParentPortal = topVolIsParent.ReadPortal(); + + // we compute the superarc of the branch within the block + // around the given isovalue + IdArrayType branchIsoSuperarcs; + branchIsoSuperarcs.Allocate(nTopVolBranchWithinBlock); + + auto branchIsoSuperarcWorklet = + vtkm::worklet::scalar_topology::select_top_volume_contours::GetSuperarcByIsoValueWorklet( + totalNumPoints); + invoke(branchIsoSuperarcWorklet, + topVolLocalBranchUpperEnd, + topVolLocalBranchLowerEnd, + isoValues, + topVolBranchSaddleEpsilons, + branchIsoSuperarcs, + findSuperarcForNode); + auto branchIsoSuperarcsPortal = branchIsoSuperarcs.ReadPortal(); + +#ifdef DEBUG_PRINT + std::stringstream branchStream; + + branchStream << "Debug for branch info, Block " << this->LocalBlockNo << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(sortedBranchGRId.GetNumberOfValues(), + branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Raw Branch GR", this->BranchRootGRId, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Raw Upper End", upperLocalEndIds, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Raw Lower End", lowerLocalEndIds, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Sorted Branch GR", sortedBranchGRId, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Sorted Branch Id", sortedBranchOrder, -1, branchStream); + + vtkm::worklet::contourtree_augmented::PrintHeader(nIsoValues, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Top Branch GR", this->TopVolumeBranchRootGRId, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Top Branch Stencil", topVolBranchWithinBlockStencil, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Top Branch Idx", topVolBranchInfoIdx, -1, branchStream); + + vtkm::worklet::contourtree_augmented::PrintHeader(nTopVolBranchKnownByBlock, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Sorted Upper End", topVolBranchUpperLocalEnd, -1, branchStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Sorted Lower End", topVolBranchLowerLocalEnd, -1, branchStream); + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, branchStream.str()); +#endif + + const vtkm::Id nExtraMaximaBranch = this->ExtraMaximaBranchLowerEnd.GetNumberOfValues(); + const vtkm::Id nExtraMinimaBranch = this->ExtraMinimaBranchLowerEnd.GetNumberOfValues(); + InArrayHandleType extraMaximaBranchIsoValue; + InArrayHandleType extraMinimaBranchIsoValue; + + IdArrayType extraMaximaBranchSuperarcs; + IdArrayType extraMinimaBranchSuperarcs; + extraMaximaBranchSuperarcs.Allocate(nExtraMaximaBranch); + extraMinimaBranchSuperarcs.Allocate(nExtraMinimaBranch); + + if (nExtraMaximaBranch) + { + extraMaximaBranchIsoValue = + this->ExtraMaximaBranchIsoValue.AsArrayHandle(); + + invoke(branchIsoSuperarcWorklet, + this->ExtraMaximaBranchUpperEnd, + this->ExtraMaximaBranchLowerEnd, + extraMaximaBranchIsoValue, + vtkm::cont::ArrayHandleConstant(1, nExtraMaximaBranch), + extraMaximaBranchSuperarcs, + findSuperarcForNode); + +#ifdef DEBUG_PRINT + std::stringstream extraMaxStream; + extraMaxStream << "Debug for Extra Maxima Branch, Block " << this->LocalBlockNo << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMaximaBranch, extraMaxStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Upper End", this->ExtraMaximaBranchUpperEnd, -1, extraMaxStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Lower End", this->ExtraMaximaBranchLowerEnd, -1, extraMaxStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "Max Branch IsoValue", extraMaximaBranchIsoValue, -1, extraMaxStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Superarc", extraMaximaBranchSuperarcs, -1, extraMaxStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Order", this->ExtraMaximaBranchOrder, -1, extraMaxStream); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMaxStream.str()); +#endif // DEBUG_PRINT + } + + if (nExtraMinimaBranch) + { + extraMinimaBranchIsoValue = + this->ExtraMinimaBranchIsoValue.AsArrayHandle(); + + invoke(branchIsoSuperarcWorklet, + this->ExtraMinimaBranchUpperEnd, + this->ExtraMinimaBranchLowerEnd, + extraMinimaBranchIsoValue, + vtkm::cont::ArrayHandleConstant(-1, nExtraMinimaBranch), + extraMinimaBranchSuperarcs, + findSuperarcForNode); + +#ifdef DEBUG_PRINT + std::stringstream extraMinStream; + extraMaxStream << "Debug for Extra Maxima Branch, Block " << this->LocalBlockNo << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nExtraMinimaBranch, extraMinStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Upper End", this->ExtraMinimaBranchUpperEnd, -1, extraMinStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Lower End", this->ExtraMinimaBranchLowerEnd, -1, extraMinStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "Max Branch IsoValue", extraMinimaBranchIsoValue, -1, extraMinStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Superarc", extraMinimaBranchSuperarcs, -1, extraMinStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Max Branch Order", this->ExtraMinimaBranchOrder, -1, extraMinStream); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, extraMinStream.str()); +#endif // DEBUG_PRINT + } + + auto extraMaximaBranchSuperarcPortal = extraMaximaBranchSuperarcs.ReadPortal(); + auto extraMinimaBranchSuperarcPortal = extraMinimaBranchSuperarcs.ReadPortal(); + auto extraMaximaBranchIsoValuePortal = extraMaximaBranchIsoValue.ReadPortal(); + auto extraMinimaBranchIsoValuePortal = extraMinimaBranchIsoValue.ReadPortal(); + auto extraMaximaBranchOrderPortal = this->ExtraMaximaBranchOrder.ReadPortal(); + auto extraMinimaBranchOrderPortal = this->ExtraMinimaBranchOrder.ReadPortal(); + + const vtkm::Id nContours = nTopVolBranchWithinBlock + nExtraMaximaBranch + nExtraMinimaBranch; + this->IsosurfaceEdgesOffset.AllocateAndFill(nContours, 0); + this->IsosurfaceEdgesLabels.AllocateAndFill(nContours, 0); + this->IsosurfaceEdgesOrders.AllocateAndFill(nContours, 0); + InArrayHandleType isosurfaceIsoValue; + isosurfaceIsoValue.AllocateAndFill(nContours, ValueType(0)); + auto edgeOffsetWritePortal = this->IsosurfaceEdgesOffset.WritePortal(); + auto edgeLabelWritePortal = this->IsosurfaceEdgesLabels.WritePortal(); + auto edgeOrderWritePortal = this->IsosurfaceEdgesOrders.WritePortal(); + auto isosurfaceValuePortal = isosurfaceIsoValue.WritePortal(); + + vtkm::Id nContourCandidateMeshes = 0; + // NOTE: nContours denotes the number of isosurfaces for visualization. + // The number is usually small, so linear loop is not too costly. + // NOTE update 06/16/2024: We always need the isovalue of the contour. + // As a result, we iterate through nContours (=O(k)). + // This may be parallelizable in future work + for (vtkm::Id branchIdx = 0; branchIdx < nContours; branchIdx++) + { + ValueType isoValue; + vtkm::Id currBranchSaddleEpsilon, branchSuperarc, branchOrder, branchLabel = 0; + + using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_SADDLE; + using vtkm::worklet::scalar_topology::select_top_volume_contours::BRANCH_COVER; + using vtkm::worklet::scalar_topology::select_top_volume_contours::MAXIMA_CONTOUR; + + if (branchIdx < nTopVolBranchWithinBlock) + { + isoValue = isoValuePortal.Get(branchIdx); + currBranchSaddleEpsilon = topVolBranchSaddleEpsilonPortal.Get(branchIdx); + branchSuperarc = branchIsoSuperarcsPortal.Get(branchIdx); + branchOrder = topVolBranchWithinBlockIdPortal.Get(branchIdx) + 1; + branchLabel |= BRANCH_SADDLE; + branchLabel |= topVolIsParentPortal.Get(branchIdx) ? BRANCH_COVER : 0; + branchLabel |= currBranchSaddleEpsilon > 0 ? MAXIMA_CONTOUR : 0; + } + else if (branchIdx < nTopVolBranchWithinBlock + nExtraMaximaBranch) + { + const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock; + isoValue = extraMaximaBranchIsoValuePortal.Get(idx); + currBranchSaddleEpsilon = 1; + branchSuperarc = extraMaximaBranchSuperarcPortal.Get(idx); + branchOrder = extraMaximaBranchOrderPortal.Get(idx); + branchLabel |= MAXIMA_CONTOUR; + } + else + { + const vtkm::Id idx = branchIdx - nTopVolBranchWithinBlock - nExtraMaximaBranch; + VTKM_ASSERT(idx < nExtraMinimaBranch); + isoValue = extraMinimaBranchIsoValuePortal.Get(idx); + currBranchSaddleEpsilon = -1; + branchSuperarc = extraMinimaBranchSuperarcPortal.Get(idx); + branchOrder = extraMinimaBranchOrderPortal.Get(idx); + } + + edgeOffsetWritePortal.Set(branchIdx, this->IsosurfaceEdgesFrom.GetNumberOfValues()); + edgeLabelWritePortal.Set(branchIdx, branchLabel); + edgeOrderWritePortal.Set(branchIdx, branchOrder); + isosurfaceValuePortal.Set(branchIdx, isoValue); + + if (vtkm::worklet::contourtree_augmented::NoSuchElement(branchSuperarc)) + { + continue; + } + + // Note: by concept, there is no 3D cell if pointDimensions[2] <= 1 + vtkm::Id nCells = globalPointDimensions[2] <= 1 + ? (pointDimensions[0] - 1) * (pointDimensions[1] - 1) + : (pointDimensions[0] - 1) * (pointDimensions[1] - 1) * (pointDimensions[2] - 1); + + // we get the polarity cases of cells + // we use lookup tables to for fast processing + IdArrayType vertexOffset = globalPointDimensions[2] <= 1 + ? vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset2d + : vtkm::worklet::scalar_topology::select_top_volume_contours::vertexOffset3d; + IdArrayType edgeTable = globalPointDimensions[2] <= 1 + ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTable2d + : (isMarchingCubes + ? vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableMC3d + : vtkm::worklet::scalar_topology::select_top_volume_contours::edgeTableLT3d); + IdArrayType numBoundTable = globalPointDimensions[2] <= 1 + ? vtkm::worklet::scalar_topology::select_top_volume_contours::numLinesTable2d + : (isMarchingCubes + ? vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableMC3d + : vtkm::worklet::scalar_topology::select_top_volume_contours::numTrianglesTableLT3d); + IdArrayType boundaryTable = globalPointDimensions[2] <= 1 + ? vtkm::worklet::scalar_topology::select_top_volume_contours::lineTable2d + : (isMarchingCubes + ? vtkm::worklet::scalar_topology::select_top_volume_contours::triTableMC3d + : vtkm::worklet::scalar_topology::select_top_volume_contours::triTableLT3d); + IdArrayType labelEdgeTable = isMarchingCubes + ? vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableMC3d + : vtkm::worklet::scalar_topology::select_top_volume_contours::labelEdgeTableLT3d; + + IdArrayType caseCells; + + caseCells.Allocate(nCells); + IdArrayType numEdgesInCell; + numEdgesInCell.Allocate(nCells); + auto caseCellsWorklet = + vtkm::worklet::scalar_topology::select_top_volume_contours::GetCellCasesWorklet( + pointDimensions, currBranchSaddleEpsilon, isoValue); + + invoke(caseCellsWorklet, + vtkm::cont::ArrayHandleIndex(nCells), + sortedDataValuesWithinBlock, + vertexOffset, + caseCells); + + // we compute the number of edges for each cell + // to initialize the array size of edges + IdArrayType numBoundariesInCell; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + numBoundTable, caseCells, numBoundariesInCell); + + vtkm::cont::Algorithm::Transform(numBoundariesInCell, + globalPointDimensions[2] <= 1 + ? vtkm::cont::ArrayHandleConstant(1, nCells) + : vtkm::cont::ArrayHandleConstant(3, nCells), + numEdgesInCell, + vtkm::Multiply()); + + // use prefix sum to get the offset of the starting edge in each cell/cube + vtkm::Id nEdges = + vtkm::cont::Algorithm::Reduce(numEdgesInCell, vtkm::Id(0)); + nContourCandidateMeshes += globalPointDimensions[2] <= 1 ? nEdges : nEdges / 3; + IdArrayType edgesOffset; + vtkm::cont::Algorithm::ScanExclusive(numEdgesInCell, edgesOffset); + + vtkm::cont::ArrayHandle isosurfaceEdgesFrom; + vtkm::cont::ArrayHandle isosurfaceEdgesTo; + // IdArrayType isosurfaceEdgesLabels; + IdArrayType isValidEdges; + isosurfaceEdgesFrom.Allocate(nEdges); + isosurfaceEdgesTo.Allocate(nEdges); + // isosurfaceEdgesLabels.Allocate(nEdges); + isValidEdges.Allocate(nEdges); + + // draw isosurface + auto getEdgesInCellWorklet = + vtkm::worklet::scalar_topology::select_top_volume_contours::GetEdgesInCellWorklet< + ValueType>(pointDimensions, + globalPointIndexStart, + isoValue, + branchSuperarc, + currBranchSaddleEpsilon, + totalNumPoints, + isMarchingCubes); + + invoke(getEdgesInCellWorklet, + edgesOffset, + caseCells, + sortedLocalNodeInfoIdsWithinBlock, + sortedDataValuesWithinBlock, + vertexOffset, + edgeTable, + numBoundTable, + boundaryTable, + labelEdgeTable, + isosurfaceEdgesFrom, + isosurfaceEdgesTo, + isValidEdges, + findSuperarcForNode); + + // isValidEdges: stencil about whether the edge is on the desired superarc + vtkm::cont::ArrayHandle validEdgesFrom; + vtkm::cont::ArrayHandle validEdgesTo; + + // we remove edges that are not on the desired superarc + vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesFrom, isValidEdges, validEdgesFrom); + vtkm::cont::Algorithm::CopyIf(isosurfaceEdgesTo, isValidEdges, validEdgesTo); + + // append edges into the result array + vtkm::Id nValidEdges = validEdgesFrom.GetNumberOfValues(); + vtkm::Id nExistEdges = branchIdx == 0 ? 0 : this->IsosurfaceEdgesFrom.GetNumberOfValues(); + + this->IsosurfaceEdgesFrom.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); + this->IsosurfaceEdgesTo.Allocate(nValidEdges + nExistEdges, vtkm::CopyFlag::On); + +#ifdef DEBUG_PRINT + std::stringstream edgeInfoStream; + contourStream << "Debug for Contour Info, Block " << this->LocalBlockNo << std::endl; + vtkm::worklet::contourtree_augmented::PrintHeader(nContours, edgeInfoStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "edgeOffset", this->IsosurfaceEdgesOffset, -1, edgeInfoStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "edgeOrders", this->IsosurfaceEdgesOrders, -1, edgeInfoStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "edgeLabels", this->IsosurfaceEdgesLabels, -1, edgeInfoStream); + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, edgeInfoStream.str()); +#endif // DEBUG_PRINT + + vtkm::cont::Algorithm::CopySubRange( + validEdgesFrom, 0, nValidEdges, this->IsosurfaceEdgesFrom, nExistEdges); + vtkm::cont::Algorithm::CopySubRange( + validEdgesTo, 0, nValidEdges, this->IsosurfaceEdgesTo, nExistEdges); + +#ifdef DEBUG_PRINT + std::stringstream contourStream; + contourStream << "Debug for Contour, Block " << this->LocalBlockNo << std::endl; + contourStream << "Branch Superarc = " << branchSuperarc << ", isoValue = " << isoValue + << std::endl; + + vtkm::worklet::contourtree_augmented::PrintHeader(nCells, contourStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Case of Cells", caseCells, -1, contourStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "# of Edges", numEdgesInCell, -1, contourStream); + + vtkm::worklet::contourtree_augmented::PrintHeader(nEdges, contourStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "Edge Stencil", isValidEdges, -1, contourStream); + + vtkm::worklet::contourtree_augmented::PrintHeader(nValidEdges, contourStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "EdgesFrom", validEdgesFrom, -1, contourStream); + vtkm::worklet::contourtree_augmented::PrintValues( + "EdgesTo", validEdgesTo, -1, contourStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "EdgesLabels", validEdgesLabels, -1, contourStream); + + VTKM_LOG_S(vtkm::cont::LogLevel::Info, contourStream.str()); +#endif // DEBUG_PRINT + } + this->IsosurfaceIsoValue = isosurfaceIsoValue; + + const vtkm::Id nMeshesOnBranches = globalPointDimensions[2] <= 1 + ? this->IsosurfaceEdgesFrom.GetNumberOfValues() + : this->IsosurfaceEdgesFrom.GetNumberOfValues() / 3; + VTKM_LOG_S(timingsLogLevel, + std::endl + << "----------- Draw Isosurface (block=" << this->LocalBlockNo << ")------------" + << std::endl + << " " << std::setw(60) << std::left << "Number of Contours: " << nContours + << std::endl + << " " << std::setw(60) << std::left + << "Number of Isosurface Meshes: " << nContourCandidateMeshes << std::endl + << " " << std::setw(60) << std::left + << "Number of Meshes On Branches: " << nMeshesOnBranches << std::endl); + }; + this->TopVolumeBranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); +} + } // namespace internal } // namespace scalar_topology } // namespace filter diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h index cb855a519..9c5f54015 100644 --- a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h @@ -69,8 +69,16 @@ struct SelectTopVolumeContoursBlock { SelectTopVolumeContoursBlock(vtkm::Id localBlockNo, int globalBlockId); - void SortBranchByVolume(const vtkm::cont::DataSet& hierarchicalTreeDataSet, - const vtkm::Id totalVolume); + void SortBranchByVolume(const vtkm::cont::DataSet& bdDataSet, const vtkm::Id totalVolume); + + void SelectLocalTopVolumeBranch(const vtkm::cont::DataSet& bdDataSet, + const vtkm::Id nSavedBranches); + + void ComputeTopVolumeBranchHierarchy(const vtkm::cont::DataSet& bdDataSet); + + void ExtractIsosurfaceOnSelectedBranch(const vtkm::cont::DataSet& bdDataSet, + const bool isMarchingCubes, + const vtkm::cont::LogLevel timingsLogLevel); // Block metadata vtkm::Id LocalBlockNo;