Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge topic 'working-branch'

Browse Source 
7a179e6a0 Merge branch 'master' into working-branch
c40b98e99 Fix to read the scan inclusive output value directly.
ce30e7c62 Refatored writePortal antipattern.
332532b29 Removed unused Euler Tour header files.
ad497e261 Merged with master.
cd43328ef Fixed signed conversion build warning.
6b82c1ab7 Fixed CUDA build errors.
ce653c78f Fixed contour tree example app.
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <kmorel@sandia.gov>
Merge-request: !2167
This commit is contained in:
Kenneth Moreland 2020-08-10 15:21:32 +00:00 committed by Kitware Robot
commit b60e731b4e
11 changed files with 1352 additions and 1048 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
examples/contour_tree_augmented/ContourTreeApp.cxx
View File

@ -719,12 +719,12 @@ int main(int argc, char* argv[])
ctaug_ns::IdArrayType superarcDependentWeight;
ctaug_ns::IdArrayType supernodeTransferWeight;
ctaug_ns::IdArrayType hyperarcDependentWeight;
ctaug_ns::ProcessContourTree::ComputeVolumeWeights(filter.GetContourTree(),
filter.GetNumIterations(),
superarcIntrinsicWeight, // (output)
superarcDependentWeight, // (output)
supernodeTransferWeight, // (output)
hyperarcDependentWeight); // (output)
ctaug_ns::ProcessContourTree::ComputeVolumeWeightsSerial(filter.GetContourTree(),
filter.GetNumIterations(),
superarcIntrinsicWeight, // (output)
superarcDependentWeight, // (output)
supernodeTransferWeight, // (output)
hyperarcDependentWeight); // (output)
// Record the timings for the branch decomposition
std::stringstream timingsStream; // Use a string stream to log in one message
timingsStream << std::endl;
@ -740,14 +740,14 @@ int main(int argc, char* argv[])
ctaug_ns::IdArrayType branchMaximum;
ctaug_ns::IdArrayType branchSaddle;
ctaug_ns::IdArrayType branchParent;
ctaug_ns::ProcessContourTree::ComputeVolumeBranchDecomposition(filter.GetContourTree(),
superarcDependentWeight,
superarcIntrinsicWeight,
whichBranch, // (output)
branchMinimum, // (output)
branchMaximum, // (output)
branchSaddle, // (output)
branchParent); // (output)
ctaug_ns::ProcessContourTree::ComputeVolumeBranchDecompositionSerial(filter.GetContourTree(),
superarcDependentWeight,
superarcIntrinsicWeight,
whichBranch, // (output)
branchMinimum, // (output)
branchMaximum, // (output)
branchSaddle, // (output)
branchParent); // (output)
// Record and log the branch decompostion timings
timingsStream << " " << std::setw(38) << std::left << "Compute Volume Branch Decomposition"
<< ": " << branchDecompTimer.GetElapsedTime() << " seconds" << std::endl;

831
vtkm/worklet/contourtree_augmented/ProcessContourTree.h
View File

@ -58,6 +58,10 @@
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vtkm/BinaryOperators.h>
#include <vtkm/BinaryPredicates.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
// local includes
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
@ -69,6 +73,8 @@
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleView.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/worklet/contourtree_augmented/ArrayTransforms.h>
#include <vtkm/worklet/contourtree_augmented/ContourTree.h>
#include <vtkm/worklet/contourtree_augmented/PrintVectors.h>
@ -77,13 +83,10 @@
#include <vtkm/worklet/contourtree_augmented/processcontourtree/SuperNodeBranchComparator.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeBestUpDown.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourFirstNext.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourList.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeMinMaxValues.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/EulerTour.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/HypersweepWorklets.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/PointerDoubling.h>
//#define DEBUG_PRINT
namespace process_contourtree_inc_ns =
@ -96,6 +99,7 @@ namespace worklet
namespace contourtree_augmented
{
// TODO Many of the post processing routines still need to be parallelized
// Class with routines for post processing the contour tree
class ProcessContourTree
@ -208,12 +212,12 @@ public:
} // CollectSortedSuperarcs()
// routine to compute the volume for each hyperarc and superarc
void static ComputeVolumeWeights(const ContourTree& contourTree,
const vtkm::Id nIterations,
IdArrayType& superarcIntrinsicWeight,
IdArrayType& superarcDependentWeight,
IdArrayType& supernodeTransferWeight,
IdArrayType& hyperarcDependentWeight)
void static ComputeVolumeWeightsSerial(const ContourTree& contourTree,
const vtkm::Id nIterations,
IdArrayType& superarcIntrinsicWeight,
IdArrayType& superarcDependentWeight,
IdArrayType& supernodeTransferWeight,
IdArrayType& hyperarcDependentWeight)
{ // ContourTreeMaker::ComputeWeights()
// start by storing the first sorted vertex ID for each superarc
IdArrayType firstVertexForSuperparent;
@ -326,12 +330,13 @@ public:
// intrinsic wt 1 2 1 5 2 6 1 1
//
// now, if we do a prefix sum on each of these and add the two sums together, we get:
// sArc s11 s12 s13 s14 s21 s22 s23 s31
// hyperparent sNode ID s11 s11 s11 s11 s21 s21 s21 s31
// transfer weight 0 1 2 1 2 3 1 0
// intrinsic weight 1 2 1 5 2 6 1 1
// sum(xfer + intrinsic) 1 3 3 6 4 9 2 1
// prefix sum (xfer + int) 1 4 7 13 14 26 28 29
// sArc s11 s12 s13 s14 s21 s22 s23 s31
// hyperparent sNode ID s11 s11 s11 s11 s21 s21 s21 s31
// transfer weight 0 1 2 1 2 3 1 0
// intrinsic weight 1 2 1 5 2 6 1 1
// sum(xfer + intrinsic) 1 3 3 6 4 9 2 1
// prefix sum (xfer + int) 1 4 7 13 17 26 28 29
// prefix sum (xfer + int - previous hArc) 1 4 7 13 4 13 15 16
// so, step 1: add xfer + int & store in dependent weight
for (vtkm::Id supernode = firstSupernode; supernode < lastSupernode; supernode++)
@ -395,16 +400,15 @@ public:
} // ContourTreeMaker::ComputeWeights()
// routine to compute the branch decomposition by volume
void static ComputeVolumeBranchDecomposition(const ContourTree& contourTree,
const IdArrayType& superarcDependentWeight,
const IdArrayType& superarcIntrinsicWeight,
IdArrayType& whichBranch,
IdArrayType& branchMinimum,
IdArrayType& branchMaximum,
IdArrayType& branchSaddle,
IdArrayType& branchParent)
void static ComputeVolumeBranchDecompositionSerial(const ContourTree& contourTree,
const IdArrayType& superarcDependentWeight,
const IdArrayType& superarcIntrinsicWeight,
IdArrayType& whichBranch,
IdArrayType& branchMinimum,
IdArrayType& branchMaximum,
IdArrayType& branchSaddle,
IdArrayType& branchParent)
{ // ComputeVolumeBranchDecomposition()
auto superarcsPortal = contourTree.Superarcs.ReadPortal();
auto superarcDependentWeightPortal = superarcDependentWeight.ReadPortal();
auto superarcIntrinsicWeightPortal = superarcIntrinsicWeight.ReadPortal();
@ -441,6 +445,8 @@ public:
#ifdef DEBUG_PRINT
std::cout << "Total Volume: " << totalVolume << std::endl;
#endif
auto superarcsPortal = contourTree.Superarcs.ReadPortal();
// NB: Last element in array is guaranteed to be root superarc to infinity,
// so we can easily skip it by not indexing to the full size
for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++)
@ -565,27 +571,15 @@ public:
vtkm::cont::ArrayHandleConstant<vtkm::Id>((vtkm::Id)NO_SUCH_ELEMENT, nSupernodes);
vtkm::cont::ArrayCopy(noSuchElementArray, whichBranch);
// Set up portals
auto bestUpwardPortal = bestUpward.WritePortal();
auto bestDownwardPortal = bestDownward.WritePortal();
auto whichBranchPortal = whichBranch.WritePortal();
// STAGE III: For each vertex, identify which neighbours are on same branch
// Let v = BestUp(u). Then if u = BestDown(v), copy BestUp(u) to whichBranch(u)
// Otherwise, let whichBranch(u) = BestUp(u) | TERMINAL to mark the end of the side branch
// NB 1: Leaves already have the flag set, but it's redundant so its safe
// NB 2: We don't need to do it downwards because it's symmetric
for (vtkm::Id supernode = 0; supernode != nSupernodes; supernode++)
{ // per supernode
vtkm::Id bestUp = bestUpwardPortal.Get(supernode);
if (NoSuchElement(bestUp))
// flag it as an upper leaf
whichBranchPortal.Set(supernode, TERMINAL_ELEMENT | supernode);
else if (bestDownwardPortal.Get(bestUp) == supernode)
whichBranchPortal.Set(supernode, bestUp);
else
whichBranchPortal.Set(supernode, TERMINAL_ELEMENT | supernode);
} // per supernode
vtkm::cont::Invoker invoke;
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PropagateBestUpDown
propagateBestUpDownWorklet;
invoke(propagateBestUpDownWorklet, bestUpward, bestDownward, whichBranch);
#ifdef DEBUG_PRINT
std::cout << "III. Branch Neighbours Identified" << std::endl;
@ -600,15 +594,16 @@ public:
for (vtkm::Id shifter = nSupernodes; shifter != 0; shifter >>= 1)
numLogSteps++;
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PointerDoubling pointerDoubling(
nSupernodes);
// use pointer-doubling to build the branches
for (vtkm::Id iteration = 0; iteration < numLogSteps; iteration++)
{ // per iteration
// loop through the vertices, updating the chaining array
for (vtkm::Id supernode = 0; supernode < nSupernodes; supernode++)
if (!IsTerminalElement(whichBranchPortal.Get(supernode)))
whichBranchPortal.Set(supernode, whichBranchPortal.Get(whichBranchPortal.Get(supernode)));
invoke(pointerDoubling, whichBranch);
} // per iteration
#ifdef DEBUG_PRINT
std::cout << "IV. Branch Chains Propagated" << std::endl;
PrintHeader(whichBranch.GetNumberOfValues());
@ -616,22 +611,21 @@ public:
std::cout << std::endl;
#endif
// STAGE V: Create an array of branches. To do this, first we need a temporary array storing
// which existing leaf corresponds to which branch. It is possible to estimate the correct number
// by counting leaves, but for parallel, we'll need a compression anyway
// V A. Set up the ID lookup for branches
vtkm::Id nBranches = 0;
// Initialise
IdArrayType chainToBranch;
vtkm::cont::ArrayCopy(noSuchElementArray, chainToBranch);
auto chainToBranchPortal = chainToBranch.WritePortal();
for (vtkm::Id supernode = 0; supernode < nSupernodes; supernode++)
{
// test whether the supernode points to itself to find the top ends
if (MaskedIndex(whichBranchPortal.Get(supernode)) == supernode)
{
chainToBranchPortal.Set(supernode, nBranches++);
}
}
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nSupernodes), chainToBranch);
// Set 1 to every relevant
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PrepareChainToBranch
prepareChainToBranchWorklet;
invoke(prepareChainToBranchWorklet, whichBranch, chainToBranch);
// Prefix scanto get IDs
vtkm::Id nBranches = vtkm::cont::Algorithm::ScanInclusive(chainToBranch, chainToBranch);
vtkm::worklet::contourtree_augmented::process_contourtree_inc::FinaliseChainToBranch
finaliseChainToBranchWorklet;
invoke(finaliseChainToBranchWorklet, whichBranch, chainToBranch);
// V B. Create the arrays for the branches
auto noSuchElementArrayNBranches =
@ -640,10 +634,6 @@ public:
vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchMaximum);
vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchSaddle);
vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchParent);
auto branchMinimumPortal = branchMinimum.WritePortal();
auto branchMaximumPortal = branchMaximum.WritePortal();
auto branchSaddlePortal = branchSaddle.WritePortal();
auto branchParentPortal = branchParent.WritePortal();
#ifdef DEBUG_PRINT
std::cout << "V. Branch Arrays Created" << std::endl;
@ -655,20 +645,14 @@ public:
PrintIndices("Branch Saddle", branchSaddle);
PrintIndices("Branch Parent", branchParent);
#endif
// STAGE VI: Sort all supernodes by [whichBranch, regular index] to get the sequence along the branch
// Assign the upper end of the branch as an ID (for now).
// VI A. Create the sorting array, then sort
IdArrayType supernodeSorter;
supernodeSorter.Allocate(nSupernodes);
auto supernodeSorterPortal = supernodeSorter.WritePortal();
for (vtkm::Id supernode = 0; supernode < nSupernodes; supernode++)
{
supernodeSorterPortal.Set(supernode, supernode);
}
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleIndex(nSupernodes), supernodeSorter);
vtkm::cont::Algorithm::Sort(
supernodeSorter,
process_contourtree_inc_ns::SuperNodeBranchComparator(whichBranch, contourTree.Supernodes));
IdArrayType permutedBranches;
permutedBranches.Allocate(nSupernodes);
PermuteArray<vtkm::Id>(whichBranch, supernodeSorter, permutedBranches);
@ -685,39 +669,13 @@ public:
PrintIndices("Regular ID", permutedRegularID);
#endif
// VI B. And reset the whichBranch to use the new branch IDs
for (vtkm::Id supernode = 0; supernode < nSupernodes; supernode++)
{
whichBranchPortal.Set(supernode,
chainToBranchPortal.Get(MaskedIndex(whichBranchPortal.Get(supernode))));
}
vtkm::worklet::contourtree_augmented::process_contourtree_inc::WhichBranchNewId
whichBranchNewIdWorklet;
invoke(whichBranchNewIdWorklet, chainToBranch, whichBranch);
// VI C. For each segment, the highest element sets up the upper end, the lowest element sets the low end
for (vtkm::Id supernode = 0; supernode < nSupernodes; supernode++)
{ // per supernode
// retrieve supernode & branch IDs
vtkm::Id supernodeID = supernodeSorterPortal.Get(supernode);
vtkm::Id branchID = whichBranchPortal.Get(supernodeID);
// save the branch ID as the owner
// use LHE of segment to set branch minimum
if (supernode == 0)
{ // sn = 0
branchMinimumPortal.Set(branchID, supernodeID);
} // sn = 0
else if (branchID != whichBranchPortal.Get(supernodeSorterPortal.Get(supernode - 1)))
{ // LHE
branchMinimumPortal.Set(branchID, supernodeID);
} // LHE
// use RHE of segment to set branch maximum
if (supernode == nSupernodes - 1)
{ // sn = max
branchMaximumPortal.Set(branchID, supernodeID);
} // sn = max
else if (branchID != whichBranchPortal.Get(supernodeSorterPortal.Get(supernode + 1)))
{ // RHE
branchMaximumPortal.Set(branchID, supernodeID);
} // RHE
} // per supernode
vtkm::worklet::contourtree_augmented::process_contourtree_inc::BranchMinMaxSet
branchMinMaxSetWorklet(nSupernodes);
invoke(branchMinMaxSetWorklet, supernodeSorter, whichBranch, branchMinimum, branchMaximum);
#ifdef DEBUG_PRINT
std::cout << "VI. Branches Set" << std::endl;
@ -728,33 +686,16 @@ public:
PrintIndices("Branch Parent", branchParent);
#endif
// STAGE VII: For each branch, set its parent (initially) to NO_SUCH_ELEMENT
// Then test upper & lower ends of each branch (in their segments) to see whether they are leaves
// At most one is a leaf. In the case of the master branch, both are
// So, non-leaf ends set the parent branch to the branch owned by the BestUp/BestDown corresponding
// while leaf ends do nothing. At the end of this, the master branch still has -1 as the parent,
// while all other branches have their parents correctly set
// BTW: This is inefficient, and we need to compress down the list of branches
for (vtkm::Id branchID = 0; branchID < nBranches; branchID++)
{ // per branch
vtkm::Id branchMax = branchMaximumPortal.Get(branchID);
// check whether the maximum is NOT a leaf
if (!NoSuchElement(bestUpwardPortal.Get(branchMax)))
{ // points to a saddle
branchSaddlePortal.Set(branchID, MaskedIndex(bestUpwardPortal.Get(branchMax)));
// if not, then the bestUp points to a saddle vertex at which we join the parent
branchParentPortal.Set(branchID, whichBranchPortal.Get(bestUpwardPortal.Get(branchMax)));
} // points to a saddle
// now do the same with the branch minimum
vtkm::Id branchMin = branchMinimumPortal.Get(branchID);
// test whether NOT a lower leaf
if (!NoSuchElement(bestDownwardPortal.Get(branchMin)))
{ // points to a saddle
branchSaddlePortal.Set(branchID, MaskedIndex(bestDownwardPortal.Get(branchMin)));
// if not, then the bestUp points to a saddle vertex at which we join the parent
branchParentPortal.Set(branchID, whichBranchPortal.Get(bestDownwardPortal.Get(branchMin)));
} // points to a saddle
} // per branch
vtkm::worklet::contourtree_augmented::process_contourtree_inc::BranchSaddleParentSet
branchSaddleParentSetWorklet;
invoke(branchSaddleParentSetWorklet,
whichBranch,
branchMinimum,
branchMaximum,
bestDownward,
bestUpward,
branchSaddle,
branchParent);
#ifdef DEBUG_PRINT
std::cout << "VII. Branches Constructed" << std::endl;
@ -794,134 +735,9 @@ public:
dataFieldIsSorted);
}
void static findMinMaxParallel(const IdArrayType supernodes,
const cont::ArrayHandle<Vec<Id, 2>> tourEdges,
const bool isMin,
IdArrayType minMaxIndex,
IdArrayType parents)
{
//
// Set up some useful portals
//
auto parentsPortal = parents.WritePortal();
auto tourEdgesPortal = tourEdges.ReadPortal();
//
// Set initial values
//
Id root = tourEdgesPortal.Get(0)[0];
parentsPortal.Set(root, root);
//
// Find what the first and last occurence of a vertex is in the euler tour.
//
cont::ArrayHandle<vtkm::Pair<Id, Id>> firstLastVertex;
firstLastVertex.Allocate(supernodes.GetNumberOfValues());
for (int i = 0; i < firstLastVertex.GetNumberOfValues(); i++)
{
firstLastVertex.WritePortal().Set(i,
{ (vtkm::Id)NO_SUCH_ELEMENT, (vtkm::Id)NO_SUCH_ELEMENT });
}
auto firstLastVertexPortal = firstLastVertex.WritePortal();
for (int i = 0; i < tourEdges.GetNumberOfValues(); i++)
{
// Forward Edge
if (firstLastVertexPortal.Get(tourEdgesPortal.Get(i)[1]).first == NO_SUCH_ELEMENT)
{
//printf("The parent of %d is %d.\n", tourEdges[i][1], tourEdges[i][0]);
parentsPortal.Set(tourEdgesPortal.Get(i)[1], tourEdgesPortal.Get(i)[0]);
//firstLastVertex[tourEdges.Get(i)[1]].first = i;
firstLastVertexPortal.Set(
tourEdgesPortal.Get(i)[1],
{ i, firstLastVertexPortal.Get(tourEdgesPortal.Get(i)[1]).second });
}
//firstLastVertex[tourEdges.Get(i)[0]].second = i;
firstLastVertexPortal.Set(tourEdgesPortal.Get(i)[1],
{ firstLastVertexPortal.Get(tourEdgesPortal.Get(i)[1]).first, i });
}
firstLastVertexPortal.Set(root, { 0, supernodes.GetNumberOfValues() - 1 });
//
// For every vertex look at the subrray between the first and last occurence and find the min/max values in it
//
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeMinMaxValues
computeMinMax(isMin);
vtkm::cont::Invoker Invoke;
Invoke(computeMinMax, supernodes, firstLastVertex, tourEdges, minMaxIndex);
}
void static findMinMax(const IdArrayType::ReadPortalType supernodes,
const cont::ArrayHandle<Vec<Id, 2>>::ReadPortalType tourEdges,
const bool isMin,
IdArrayType::WritePortalType minMaxIndex,
IdArrayType::WritePortalType parents)
{
Id root = tourEdges.Get(0)[0];
struct VertexData
{
long long distance;
unsigned long index;
};
std::vector<VertexData> vertexData(static_cast<unsigned long>(supernodes.GetNumberOfValues()),
{ -1, 0 });
for (unsigned long i = 0; i < vertexData.size(); i++)
{
vertexData[i].index = i;
}
parents.Set(root, root);
vertexData[static_cast<unsigned long>(root)].distance = 0;
for (int i = 0; i < tourEdges.GetNumberOfValues(); i++)
{
const Vec<Id, 2> e = tourEdges.Get(i);
if (-1 == vertexData[static_cast<unsigned long>(e[1])].distance)
{
parents.Set(e[1], e[0]);
vertexData[static_cast<unsigned long>(e[1])].distance =
vertexData[static_cast<unsigned long>(e[0])].distance + 1;
}
}
std::sort(vertexData.begin(), vertexData.end(), [](const VertexData& a, const VertexData& b) {
return a.distance > b.distance;
});
for (int i = 0; i < minMaxIndex.GetNumberOfValues(); i++)
{
minMaxIndex.Set(i, i);
}
for (unsigned int i = 0; i < vertexData.size(); i++)
{
Id vertex = static_cast<Id>(vertexData[i].index);
Id parent = parents.Get(vertex);
Id vertexValue = MaskedIndex(supernodes.Get(minMaxIndex.Get(vertex)));
Id parentValue = MaskedIndex(supernodes.Get(minMaxIndex.Get(parent)));
if ((true == isMin && vertexValue < parentValue) ||
(false == isMin && vertexValue > parentValue))
{
minMaxIndex.Set(parent, minMaxIndex.Get(vertex));
}
}
}
// routine to compute the branch decomposition by volume
void static ComputeHeightBranchDecomposition(const ContourTree& contourTree,
const cont::ArrayHandle<Float64> fieldValues,
const IdArrayType& ctSortOrder,
const bool useParallelMinMaxSearch,
void static ComputeVolumeBranchDecomposition(const ContourTree& contourTree,
const vtkm::Id nIterations,
IdArrayType& whichBranch,
IdArrayType& branchMinimum,
IdArrayType& branchMaximum,
@ -929,94 +745,78 @@ public:
IdArrayType& branchParent)
{ // ComputeHeightBranchDecomposition()
vtkm::cont::Invoker Invoke;
// STEP 1. Compute the number of nodes in every superarc, that's the intrinsic weight
IdArrayType superarcIntrinsicWeight;
superarcIntrinsicWeight.Allocate(contourTree.Superarcs.GetNumberOfValues());
IdArrayType firstVertexForSuperparent;
firstVertexForSuperparent.Allocate(contourTree.Superarcs.GetNumberOfValues());
// Compute the number of regular nodes on every superarcs (the intrinsic weight)
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetFirstVertexForSuperparent
setFirstVertexForSuperparent;
Invoke(setFirstVertexForSuperparent,
contourTree.Nodes,
contourTree.Superparents,
firstVertexForSuperparent);
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeIntrinsicWeight
computeIntrinsicWeight;
Invoke(computeIntrinsicWeight,
contourTree.Arcs,
contourTree.Superarcs,
firstVertexForSuperparent,
superarcIntrinsicWeight);
// Cache the number of non-root supernodes & superarcs
vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
auto noSuchElementArray =
vtkm::cont::ArrayHandleConstant<vtkm::Id>((vtkm::Id)NO_SUCH_ELEMENT, nSupernodes);
// STAGE I: Find the upward and downwards weight for each superarc, and set up arrays
IdArrayType bestUpward;
IdArrayType bestDownward;
// Set up bestUpward and bestDownward array, these are the things we want to compute in this routine.
IdArrayType bestUpward, bestDownward;
vtkm::cont::ArrayCopy(noSuchElementArray, bestUpward);
vtkm::cont::ArrayCopy(noSuchElementArray, bestDownward);
//
// Compute Euler Tours
//
cont::ArrayHandle<Vec<Id, 2>> minTourEdges, maxTourEdges;
// We initiale with the weight of the superarcs, once we sum those up we'll get the hypersweep weight
IdArrayType sumValues;
vtkm::cont::ArrayCopy(superarcIntrinsicWeight, sumValues);
minTourEdges.Allocate(2 * (nSupernodes - 1));
maxTourEdges.Allocate(2 * (nSupernodes - 1));
// This should be 0 here, because we're not changing the root
vtkm::cont::ArrayHandle<vtkm::Id> howManyUsed;
vtkm::cont::ArrayCopy(
vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, contourTree.Hyperarcs.GetNumberOfValues()),
howManyUsed);
// Compute the Euler Tour
process_contourtree_inc_ns::EulerTour tour;
tour.computeEulerTour(contourTree.Superarcs.ReadPortal());
// Perform a sum hypersweep
hyperarcScan<decltype(vtkm::Sum())>(contourTree.Supernodes,
contourTree.Hypernodes,
contourTree.Hyperarcs,
contourTree.Hyperparents,
contourTree.Hyperparents,
contourTree.WhenTransferred,
howManyUsed,
nIterations,
vtkm::Sum(),
sumValues);
// Reroot the Euler Tour at the global min
tour.getTourAtRoot(MaskedIndex(contourTree.Superparents.ReadPortal().Get(0)),
minTourEdges.WritePortal());
// For every directed arc store the volume of it's associate subtree
vtkm::cont::ArrayHandle<vtkm::worklet::contourtree_augmented::EdgeDataVolume> arcs;
arcs.Allocate(contourTree.Superarcs.GetNumberOfValues() * 2 - 2);
// Reroot the Euler Tour at the global max
tour.getTourAtRoot(MaskedIndex(contourTree.Superparents.ReadPortal().Get(
contourTree.Nodes.GetNumberOfValues() - 1)),
maxTourEdges.WritePortal());
vtkm::Id totalVolume = contourTree.Nodes.GetNumberOfValues();
vtkm::worklet::contourtree_augmented::process_contourtree_inc::InitialiseArcsVolume initArcs(
totalVolume);
Invoke(initArcs, sumValues, superarcIntrinsicWeight, contourTree.Superarcs, arcs);
//
// Compute Min/Max per subtree
//
IdArrayType minValues, minParents, maxValues, maxParents;
vtkm::cont::ArrayCopy(noSuchElementArray, minValues);
vtkm::cont::ArrayCopy(noSuchElementArray, minParents);
vtkm::cont::ArrayCopy(noSuchElementArray, maxValues);
vtkm::cont::ArrayCopy(noSuchElementArray, maxParents);
// Sort arcs to obtain the best up and down
vtkm::cont::Algorithm::Sort(arcs, vtkm::SortLess());
// Finding the min/max for every subtree can either be done in parallel or serial
// The parallel implementation is not work efficient and will not work well on a single/few cores
// This is why I have left the option to do a BFS style search in serial instead of doing an a prefix min/max for every subtree in the euler tour
if (false == useParallelMinMaxSearch)
{
ProcessContourTree::findMinMax(contourTree.Supernodes.ReadPortal(),
minTourEdges.ReadPortal(),
true,
minValues.WritePortal(),
minParents.WritePortal());
ProcessContourTree::findMinMax(contourTree.Supernodes.ReadPortal(),
maxTourEdges.ReadPortal(),
false,
maxValues.WritePortal(),
maxParents.WritePortal());
}
else
{
ProcessContourTree::findMinMaxParallel(
contourTree.Supernodes, minTourEdges, true, minValues, minParents);
ProcessContourTree::findMinMaxParallel(
contourTree.Supernodes, maxTourEdges, false, maxValues, maxParents);
}
//
// Compute bestUp and bestDown
//
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeBestUpDown
bestUpDownWorklet;
vtkm::cont::Invoker Invoke;
Invoke(bestUpDownWorklet,
tour.first,
contourTree.Nodes,
contourTree.Supernodes,
minValues,
minParents,
maxValues,
maxParents,
ctSortOrder,
tour.edges,
fieldValues,
bestUpward, // output
bestDownward // output
);
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetBestUpDown setBestUpDown;
Invoke(setBestUpDown, bestUpward, bestDownward, arcs);
ProcessContourTree::ComputeBranchData(contourTree,
whichBranch,
@ -1029,6 +829,367 @@ public:
} // ComputeHeightBranchDecomposition()
// routine to compute the branch decomposition by volume
void static ComputeHeightBranchDecomposition(const ContourTree& contourTree,
const cont::ArrayHandle<Float64> fieldValues,
const IdArrayType& ctSortOrder,
const vtkm::Id nIterations,
IdArrayType& whichBranch,
IdArrayType& branchMinimum,
IdArrayType& branchMaximum,
IdArrayType& branchSaddle,
IdArrayType& branchParent)
{ // ComputeHeightBranchDecomposition()
// Cache the number of non-root supernodes & superarcs
vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
auto noSuchElementArray =
vtkm::cont::ArrayHandleConstant<vtkm::Id>((vtkm::Id)NO_SUCH_ELEMENT, nSupernodes);
// Set up bestUpward and bestDownward array, these are the things we want to compute in this routine.
IdArrayType bestUpward, bestDownward;
vtkm::cont::ArrayCopy(noSuchElementArray, bestUpward);
vtkm::cont::ArrayCopy(noSuchElementArray, bestDownward);
// maxValues and minValues store the values from the max and min hypersweep respectively.
IdArrayType minValues, maxValues;
vtkm::cont::ArrayCopy(contourTree.Supernodes, maxValues);
vtkm::cont::ArrayCopy(contourTree.Supernodes, minValues);
// Store the direction of the superarcs in the min and max hypersweep (find a way to get rid of these, the only differing direction is on the path from the root to the min/max).
IdArrayType minParents, maxParents;
vtkm::cont::ArrayCopy(contourTree.Superarcs, minParents);
vtkm::cont::ArrayCopy(contourTree.Superarcs, maxParents);
auto minParentsPortal = minParents.WritePortal();
auto maxParentsPortal = maxParents.WritePortal();
// Cache the glonal minimum and global maximum (these will be the roots in the min and max hypersweep)
Id minSuperNode = MaskedIndex(contourTree.Superparents.ReadPortal().Get(0));
Id maxSuperNode = MaskedIndex(
contourTree.Superparents.ReadPortal().Get(contourTree.Nodes.GetNumberOfValues() - 1));
// Find the path from the global minimum to the root, not parallelisable (but it's fast, no need to parallelise)
auto minPath = findSuperPathToRoot(contourTree.Superarcs.ReadPortal(), minSuperNode);
// Find the path from the global minimum to the root, not parallelisable (but it's fast, no need to parallelise)
auto maxPath = findSuperPathToRoot(contourTree.Superarcs.ReadPortal(), maxSuperNode);
// Reserve the direction of the superarcs on the min path.
for (std::size_t i = 1; i < minPath.size(); i++)
{
minParentsPortal.Set(minPath[i], minPath[i - 1]);
}
minParentsPortal.Set(minPath[0], 0);
// Reserve the direction of the superarcs on the max path.
for (std::size_t i = 1; i < maxPath.size(); i++)
{
maxParentsPortal.Set(maxPath[i], maxPath[i - 1]);
}
maxParentsPortal.Set(maxPath[0], 0);
vtkm::cont::Invoker Invoke;
vtkm::worklet::contourtree_augmented::process_contourtree_inc::UnmaskArray unmaskArrayWorklet;
Invoke(unmaskArrayWorklet, minValues);
Invoke(unmaskArrayWorklet, maxValues);
// Thse arrays hold the changes hyperarcs in the min and max hypersweep respectively
vtkm::cont::ArrayHandle<vtkm::Id> minHyperarcs, maxHyperarcs;
vtkm::cont::ArrayCopy(contourTree.Hyperarcs, minHyperarcs);
vtkm::cont::ArrayCopy(contourTree.Hyperarcs, maxHyperarcs);
// These arrays hold the changed hyperarcs for the min and max hypersweep
vtkm::cont::ArrayHandle<vtkm::Id> minHyperparents, maxHyperparents;
vtkm::cont::ArrayCopy(contourTree.Hyperparents, minHyperparents);
vtkm::cont::ArrayCopy(contourTree.Hyperparents, maxHyperparents);
auto minHyperparentsPortal = minHyperparents.WritePortal();
auto maxHyperparentsPortal = maxHyperparents.WritePortal();
for (std::size_t i = 0; i < minPath.size(); i++)
{
// Set a unique dummy Id (something that the prefix scan by key will leave alone)
minHyperparentsPortal.Set(minPath[i],
contourTree.Hypernodes.GetNumberOfValues() + minPath[i]);
}
for (std::size_t i = 0; i < maxPath.size(); i++)
{
// Set a unique dummy Id (something that the prefix scan by key will leave alone)
maxHyperparentsPortal.Set(maxPath[i],
contourTree.Hypernodes.GetNumberOfValues() + maxPath[i]);
}
// These arrays hold the number of nodes in each hypearcs that are on the min or max path for the min and max hypersweep respectively.
vtkm::cont::ArrayHandle<vtkm::Id> minHowManyUsed, maxHowManyUsed;
vtkm::cont::ArrayCopy(
vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, maxHyperarcs.GetNumberOfValues()),
minHowManyUsed);
vtkm::cont::ArrayCopy(
vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, maxHyperarcs.GetNumberOfValues()),
maxHowManyUsed);
// Min Hypersweep
const auto minOperator = vtkm::Minimum();
// Cut hyperarcs at the first node on the path from the max to the root
editHyperarcs(contourTree.Hyperparents.ReadPortal(),
minPath,
minHyperarcs.WritePortal(),
minHowManyUsed.WritePortal());
// Perform an ordinary hypersweep on those new hyperarcs
hyperarcScan<decltype(vtkm::Minimum())>(contourTree.Supernodes,
contourTree.Hypernodes,
minHyperarcs,
contourTree.Hyperparents,
minHyperparents,
contourTree.WhenTransferred,
minHowManyUsed,
nIterations,
vtkm::Minimum(),
minValues);
// Prefix sum along the path from the min to the root
fixPath(vtkm::Minimum(), minPath, minValues.WritePortal());
// Max Hypersweep
const auto maxOperator = vtkm::Maximum();
// Cut hyperarcs at the first node on the path from the max to the root
editHyperarcs(contourTree.Hyperparents.ReadPortal(),
maxPath,
maxHyperarcs.WritePortal(),
maxHowManyUsed.WritePortal());
// Perform an ordinary hypersweep on those new hyperarcs
hyperarcScan<decltype(vtkm::Maximum())>(contourTree.Supernodes,
contourTree.Hypernodes,
maxHyperarcs,
contourTree.Hyperparents,
maxHyperparents,
contourTree.WhenTransferred,
maxHowManyUsed,
nIterations,
vtkm::Maximum(),
maxValues);
// Prefix sum along the path from the max to the root
fixPath(vtkm::Maximum(), maxPath, maxValues.WritePortal());
// For every directed edge (a, b) consider that subtree who's root is b and does not contain a.
// We have so far found the min and max in all sub subtrees, now we compare those to a and incorporate a into that.
vtkm::worklet::contourtree_augmented::process_contourtree_inc::IncorporateParent<decltype(
vtkm::Minimum())>
incorporateParentMinimumWorklet(minOperator);
Invoke(incorporateParentMinimumWorklet, minParents, contourTree.Supernodes, minValues);
vtkm::worklet::contourtree_augmented::process_contourtree_inc::IncorporateParent<decltype(
vtkm::Maximum())>
incorporateParentMaximumWorklet(maxOperator);
Invoke(incorporateParentMaximumWorklet, maxParents, contourTree.Supernodes, maxValues);
// Initialise all directed superarcs in the contour tree. Those will correspond to subtrees whos height we need for the branch decomposition.
vtkm::cont::ArrayHandle<vtkm::worklet::contourtree_augmented::EdgeDataHeight> arcs;
arcs.Allocate(contourTree.Superarcs.GetNumberOfValues() * 2 - 2);
vtkm::worklet::contourtree_augmented::process_contourtree_inc::InitialiseArcs initArcs(
0, contourTree.Arcs.GetNumberOfValues() - 1, minPath[minPath.size() - 1]);
Invoke(initArcs, minParents, maxParents, minValues, maxValues, contourTree.Superarcs, arcs);
// Use the min & max to compute the height of all subtrees
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeSubtreeHeight
computeSubtreeHeight;
Invoke(computeSubtreeHeight, fieldValues, ctSortOrder, contourTree.Supernodes, arcs);
// Sort all directed edges based on the height of their subtree
vtkm::cont::Algorithm::Sort(arcs, vtkm::SortLess());
// Select a best up and best down neighbour for every vertex in the contour tree using heights of all subtrees
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetBestUpDown setBestUpDown;
Invoke(setBestUpDown, bestUpward, bestDownward, arcs);
// Having computed the bestUp/Down we can propagte those to obtain the branches of the branch decomposition
ProcessContourTree::ComputeBranchData(contourTree,
whichBranch,
branchMinimum,
branchMaximum,
branchSaddle,
branchParent,
bestUpward,
bestDownward);
} // ComputeHeightBranchDecomposition()
std::vector<Id> static findSuperPathToRoot(
vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType parentsPortal,
vtkm::Id vertex)
{
// Initialise the empty path and starting vertex
std::vector<vtkm::Id> path;
vtkm::Id current = vertex;
// Go up the parent list until we reach the root
while (MaskedIndex(parentsPortal.Get(current)) != 0)
{
path.push_back(current);
current = MaskedIndex(parentsPortal.Get(current));
}
path.push_back(current);
return path;
}
// Given a path from a leaf (the global min/max) to the root of the contour tree and a hypersweep (where all hyperarcs are cut at the path)
// This function performs a prefix scan along that path to obtain the correct hypersweep values (as in the global min/max is the root of the hypersweep)
void static fixPath(const std::function<vtkm::Id(vtkm::Id, vtkm::Id)> operation,
const std::vector<vtkm::Id> path,
vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType minMaxIndex)
{
using vtkm::worklet::contourtree_augmented::MaskedIndex;
// Fix path from the old root to the new root. Parallelisble with a prefix scan, but sufficiently fast for now.
for (auto i = path.size() - 2; i > 0; i--)
{
const auto vertex = path[i + 1];
const auto parent = path[i];
const auto vertexValue = minMaxIndex.Get(vertex);
const auto parentValue = minMaxIndex.Get(parent);
minMaxIndex.Set(parent, operation(vertexValue, parentValue));
}
}
// This function edits all the hyperarcs which contain vertices which are on the supplied path. This path is usually the path between the global min/max to the root of the tree.
// This function effectively cuts hyperarcs at the first node they encounter along that path.
// In addition to this it computed the number of supernodes every hyperarc has on that path. This helps in the function hyperarcScan for choosing the new target of the cut hyperarcs.
// NOTE: It is assumed that the supplied path starts at a leaf and ends at the root of the contour tree.
void static editHyperarcs(
const vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType hyperparentsPortal,
const std::vector<vtkm::Id> path,
vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType hyperarcsPortal,
vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType howManyUsedPortal)
{
using vtkm::worklet::contourtree_augmented::MaskedIndex;
std::size_t i = 0;
while (i < path.size())
{
// Cut the hyperacs at the first point
hyperarcsPortal.Set(MaskedIndex(hyperparentsPortal.Get(path[i])), path[i]);
Id currentHyperparent = MaskedIndex(hyperparentsPortal.Get(path[i]));
// Skip the rest of the supernodes which are on the same hyperarc
while (i < path.size() && MaskedIndex(hyperparentsPortal.Get(path[i])) == currentHyperparent)
{
const auto value = howManyUsedPortal.Get(MaskedIndex(hyperparentsPortal.Get(path[i])));
howManyUsedPortal.Set(MaskedIndex(hyperparentsPortal.Get(path[i])), value + 1);
i++;
}
}
}
template <class BinaryFunctor>
void static hyperarcScan(const vtkm::cont::ArrayHandle<vtkm::Id> supernodes,
const vtkm::cont::ArrayHandle<vtkm::Id> hypernodes,
const vtkm::cont::ArrayHandle<vtkm::Id> hyperarcs,
const vtkm::cont::ArrayHandle<vtkm::Id> hyperparents,
const vtkm::cont::ArrayHandle<vtkm::Id> hyperparentKeys,
const vtkm::cont::ArrayHandle<vtkm::Id> whenTransferred,
const vtkm::cont::ArrayHandle<vtkm::Id> howManyUsed,
const vtkm::Id nIterations,
const BinaryFunctor operation,
vtkm::cont::ArrayHandle<vtkm::Id> minMaxIndex)
{
using vtkm::worklet::contourtree_augmented::MaskedIndex;
vtkm::cont::Invoker invoke;
auto supernodesPortal = supernodes.ReadPortal();
auto hypernodesPortal = hypernodes.ReadPortal();
auto hyperparentsPortal = hyperparents.ReadPortal();
// Set the first supernode per iteration
vtkm::cont::ArrayHandle<vtkm::Id> firstSupernodePerIteration;
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nIterations + 1),
firstSupernodePerIteration);
// The first different from the previous is the first in the iteration
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetFirstSupernodePerIteration
setFirstSupernodePerIteration;
invoke(setFirstSupernodePerIteration, whenTransferred, firstSupernodePerIteration);
auto firstSupernodePerIterationPortal = firstSupernodePerIteration.WritePortal();
for (vtkm::Id iteration = 1; iteration < nIterations; ++iteration)
{
if (firstSupernodePerIterationPortal.Get(iteration) == 0)
{
firstSupernodePerIterationPortal.Set(iteration,
firstSupernodePerIterationPortal.Get(iteration + 1));
}
}
// set the sentinel at the end of the array
firstSupernodePerIterationPortal.Set(nIterations, supernodesPortal.GetNumberOfValues());
// Set the first hypernode per iteration
vtkm::cont::ArrayHandle<vtkm::Id> firstHypernodePerIteration;
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nIterations + 1),
firstHypernodePerIteration);
auto firstHypernodePerIterationPortal = firstHypernodePerIteration.WritePortal();
for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
{
firstHypernodePerIterationPortal.Set(
iteration, hyperparentsPortal.Get(firstSupernodePerIterationPortal.Get(iteration)));
}
// Set the sentinel at the end of the array
firstHypernodePerIterationPortal.Set(nIterations, hypernodesPortal.GetNumberOfValues());
// This workled is used in every iteration of the following loop, so it's initialised outside.
vtkm::worklet::contourtree_augmented::process_contourtree_inc::AddDependentWeightHypersweep<
BinaryFunctor>
addDependentWeightHypersweepWorklet(operation);
// For every iteration do a prefix scan on all hyperarcs in that iteration and then transfer the scanned value to every hyperarc's target supernode
for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
{
// Determine the first and last hypernode in the current iteration (all hypernodes between them are also in the current iteration)
vtkm::Id firstHypernode = firstHypernodePerIterationPortal.Get(iteration);
vtkm::Id lastHypernode = firstHypernodePerIterationPortal.Get(iteration + 1);
lastHypernode = vtkm::Minimum()(lastHypernode, hypernodes.GetNumberOfValues() - 1);
// Determine the first and last supernode in the current iteration (all supernode between them are also in the current iteration)
vtkm::Id firstSupernode = MaskedIndex(hypernodesPortal.Get(firstHypernode));
vtkm::Id lastSupernode = MaskedIndex(hypernodesPortal.Get(lastHypernode));
lastSupernode = vtkm::Minimum()(lastSupernode, hyperparents.GetNumberOfValues() - 1);
// Prefix scan along all hyperarcs in the current iteration
auto subarrayValues = vtkm::cont::make_ArrayHandleView(
minMaxIndex, firstSupernode, lastSupernode - firstSupernode);
auto subarrayKeys = vtkm::cont::make_ArrayHandleView(
hyperparentKeys, firstSupernode, lastSupernode - firstSupernode);
vtkm::cont::Algorithm::ScanInclusiveByKey(
subarrayKeys, subarrayValues, subarrayValues, operation);
// Array containing the Ids of the hyperarcs in the current iteration
vtkm::cont::ArrayHandleCounting<vtkm::Id> iterationHyperarcs(
firstHypernode, 1, lastHypernode - firstHypernode);
// Transfer the value accumulated in the last entry of the prefix scan to the hypernode's targe supernode
invoke(addDependentWeightHypersweepWorklet,
iterationHyperarcs,
hypernodes,
hyperarcs,
howManyUsed,
minMaxIndex);
}
}
}; // class ProcessContourTree
} // namespace contourtree_augmented
} // worklet

94
vtkm/worklet/contourtree_augmented/Types.h
View File

@ -148,6 +148,100 @@ inline std::string FlagString(vtkm::Id flaggedIndex)
return fString;
} // FlagString()
class EdgeDataHeight
{
public:
// RegularNodeID (or sortIndex)
Id I;
// RegularNodeID (or sortIndex)
Id J;
// RegularNodeID (or sortIndex)
Id SubtreeMin;
// RegularNodeID (or sortIndex)
Id SubtreeMax;
bool UpEdge;
Float64 SubtreeHeight;
VTKM_EXEC
bool operator<(const EdgeDataHeight& b) const
{
if (this->I == b.I)
{
if (this->UpEdge == b.UpEdge)
{
if (this->SubtreeHeight == b.SubtreeHeight)
{
if (this->SubtreeMin == b.SubtreeMin)
{
return this->SubtreeMax > b.SubtreeMax;
}
else
{
return this->SubtreeMin < b.SubtreeMin;
}
}
else
{
return this->SubtreeHeight > b.SubtreeHeight;
}
}
else
{
return this->UpEdge < b.UpEdge;
}
}
else
{
return this->I < b.I;
}
}
};
class EdgeDataVolume
{
public:
// RegularNodeID (or sortIndex)
Id I;
// RegularNodeID (or sortIndex)
Id J;
bool UpEdge;
Id SubtreeVolume;
VTKM_EXEC
bool operator<(const EdgeDataVolume& b) const
{
if (this->I == b.I)
{
if (this->UpEdge == b.UpEdge)
{
if (this->SubtreeVolume == b.SubtreeVolume)
{
if (this->UpEdge == true)
{
return this->J > b.J;
}
else
{
return this->J < b.J;
}
}
else
{
return this->SubtreeVolume > b.SubtreeVolume;
}
}
else
{
return this->UpEdge < b.UpEdge;
}
}
else
{
return this->I < b.I;
}
}
};
///
/// Helper struct to collect sizing information from a dataset.

6
vtkm/worklet/contourtree_augmented/processcontourtree/Branch.h
View File

@ -73,6 +73,7 @@ template <typename T>
class Branch
{
public:
vtkm::Id OriginalId; // Index of the extremum in the mesh
vtkm::Id Extremum; // Index of the extremum in the mesh
T ExtremumVal; // Value at the extremum:w
vtkm::Id Saddle; // Index of the saddle in the mesh (or minimum for root branch)
@ -157,7 +158,7 @@ Branch<T>* Branch<T>::ComputeBranchDecomposition(
const IdArrayType& sortOrder,
const vtkm::cont::ArrayHandle<T, StorageType>& dataField,
bool dataFieldIsSorted)
{ // ComputeBranchDecomposition()
{ // C)omputeBranchDecomposition()
auto branchMinimumPortal = branchMinimum.ReadPortal();
auto branchMaximumPortal = branchMaximum.ReadPortal();
auto branchSaddlePortal = branchSaddle.ReadPortal();
@ -176,6 +177,7 @@ Branch<T>* Branch<T>::ComputeBranchDecomposition(
// Reconstruct explicit branch decomposition from array representation
for (std::size_t branchID = 0; branchID < static_cast<std::size_t>(nBranches); ++branchID)
{
branches[branchID]->OriginalId = static_cast<vtkm::Id>(branchID);
if (!NoSuchElement(branchSaddlePortal.Get(static_cast<vtkm::Id>(branchID))))
{
branches[branchID]->Saddle = MaskedIndex(
@ -364,7 +366,7 @@ void Branch<T>::GetRelevantValues(int type, T eps, std::vector<T>& values) const
break;
}
if (Parent)
values.push_back(val);
values.push_back({ val });
for (Branch* c : Children)
c->GetRelevantValues(type, eps, values);
} // GetRelevantValues()

7
vtkm/worklet/contourtree_augmented/processcontourtree/CMakeLists.txt
View File

@ -13,12 +13,9 @@ set(headers
PiecewiseLinearFunction.h
SuperArcVolumetricComparator.h
SuperNodeBranchComparator.h
ComputeBestUpDown.h
ComputeEulerTourFirstNext.h
ComputeEulerTourList.h
ComputeMinMaxValues.h
PointerDoubling.h
HypersweepWorklets.h
SetTriangleSuperarcId.h
EulerTour.h
)
#-----------------------------------------------------------------------------

218
vtkm/worklet/contourtree_augmented/processcontourtree/ComputeBestUpDown.h
View File

@ -1,218 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_compute_best_up_down_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_compute_best_up_down_h
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
/*
* This code is written by Petar Hristov in 09.2019
*
* This worklet is part of the Contour Tree Height Based Simplification Code.
* It selects the bestUp and the bestDown for every supernode in the contour tree.
* The best up and best down are used to construct the branches.
*
* This worklet receives a 1D array of edges, sorter by their first vertex and then by their second vertex.
* Each invocation of the worklet goes through the neighbours of one vertex and looks for the bestUp and bestDown.
*
*/
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class ComputeBestUpDown : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn _1,
WholeArrayIn _2,
WholeArrayIn _3,
WholeArrayIn _4,
WholeArrayIn _5,
WholeArrayIn _6,
WholeArrayIn _7,
WholeArrayIn _8,
WholeArrayIn _9,
WholeArrayIn _10,
WholeArrayOut _11,
WholeArrayOut _12);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12);
using InputDomain = _3;
// Default Constructor
VTKM_EXEC_CONT ComputeBestUpDown() {}
template <typename FirstArrayPortalType,
typename NodesArrayPortalType,
typename SupernodesArrayPortalType,
typename MinValuesArrayPortalType,
typename MinParentsArrayPortalType,
typename MaxValuesArrayPortalType,
typename MaxParentsArrayPortalType,
typename SortOrderArrayPortalType,
typename EdgesLinearArrayPortalType,
typename FieldValuesArrayPortalType,
typename OutputArrayPortalType>
VTKM_EXEC void operator()(const vtkm::Id i,
const FirstArrayPortalType& first,
const NodesArrayPortalType& nodes,
const SupernodesArrayPortalType& supernodes,
const MinValuesArrayPortalType& minValues,
const MinParentsArrayPortalType& minParents,
const MaxValuesArrayPortalType& maxValues,
const MaxParentsArrayPortalType& maxParents,
const SortOrderArrayPortalType& ctSortOrder,
const EdgesLinearArrayPortalType& edgesLinear,
const FieldValuesArrayPortalType& fieldValues,
const OutputArrayPortalType& bestUp, // output
const OutputArrayPortalType& bestDown // output
) const
{
Id k = first.Get(i);
Float64 maxUpSubtreeHeight = 0;
Float64 maxDownSubtreeHeight = 0;
while (k < edgesLinear.GetNumberOfValues() && edgesLinear.Get(k)[0] == i)
{
Id j = edgesLinear.Get(k++)[1];
Id regularVertexValueI = MaskedIndex(supernodes.Get(i));
Id regularVertexValueJ = MaskedIndex(supernodes.Get(j));
//
// Get the minimum of subtree T(j) \cup {i}
//
// If the arc is pointed the right way (according to the rooting of the tree) use the subtree min value
// This is the minimum of T(j)
Id minValueInSubtree = MaskedIndex(supernodes.Get(minValues.Get(j)));
// See if the vertex i has a smaller value,
// This means find the minimum of T(j) \cup {i}
if (minValueInSubtree > regularVertexValueI)
{
minValueInSubtree = MaskedIndex(supernodes.Get(i));
}
// If the dirrection of the arc is not according to the rooting of the tree,
// then the minimum on that subtree must be the global minimum
if (j == minParents.Get(i))
{
minValueInSubtree = 0;
}
//
// Get the maximum of subtree T(j) \cup {i}
//
// See if the vertex i has a bigger value,
// This means find the maximum of T(j) \cup {i}
Id maxValueInSubtree = MaskedIndex(supernodes.Get(maxValues.Get(j)));
// Include the current vertex along with the subtree it points at
if (maxValueInSubtree < regularVertexValueI)
{
maxValueInSubtree = MaskedIndex(supernodes.Get(i));
}
// If the dirrection of the arc is not according to the rooting of the tree,
// then the maximum on that subtree must be the global maximum
if (j == maxParents.Get(i))
{
maxValueInSubtree = nodes.GetNumberOfValues() - 1;
}
// Afte having found the min and the max in T(j) \cup {i} we compute their height difference
Float64 minValue = fieldValues.Get(ctSortOrder.Get(minValueInSubtree));
Float64 maxValue = fieldValues.Get(ctSortOrder.Get(maxValueInSubtree));
Float64 subtreeHeight = maxValue - minValue;
// Downward Edge
if (regularVertexValueI > regularVertexValueJ)
{
if (subtreeHeight > maxDownSubtreeHeight)
{
maxDownSubtreeHeight = subtreeHeight;
bestDown.Set(i, j);
}
}
// UpwardsEdge
else
{
if (subtreeHeight > maxUpSubtreeHeight)
{
maxUpSubtreeHeight = subtreeHeight;
bestUp.Set(i, j);
}
}
}
// Make sure at least one of these was set
assert(false == NoSuchElement(bestUp.Get(i)) || false == NoSuchElement(bestDown.Get(i)));
}
}; // ComputeBestUpDown
} // process_contourtree_inc
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm
#endif

130
vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourList.h
View File

@ -1,130 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_euler_tour_list_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_euler_tour_list_h
#include <vtkm/cont/Algorithm.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
/*
* This code is written by Petar Hristov in 09.2019
*
* This worklet performs the second step in computing an Euler tour.
* That is computing the adjacency list of the euler edge tour.
* See https://en.wikipedia.org/wiki/Euler_tour_technique for details
*
*/
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class ComputeEulerTourList : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn next,
WholeArrayIn first,
WholeArrayIn edges,
WholeArrayOut succ);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _1;
template <typename NextArrayPortalType,
typename FirstArrayPortalType,
typename EdgesArrayPortalType,
typename OutputArrayPortalType>
VTKM_EXEC void operator()(const vtkm::Id i,
const NextArrayPortalType& next,
const FirstArrayPortalType& first,
const EdgesArrayPortalType& edges,
const OutputArrayPortalType& succ) const
{
// For edge (a, b) find the opposite edge (b, a) in the sorted list of edges
// This relies on the fact that in the contour tree nodes have bounded degree.
// This comes from the fact that vertices in a 2D/3D mesh have a bounded degree.
auto currentEdge = edges.Get(i);
auto oppositeIndex = first.Get(currentEdge[1]);
while (oppositeIndex < edges.GetNumberOfValues() &&
currentEdge[1] == edges.Get(oppositeIndex)[0])
{
if (currentEdge[0] == edges.Get(oppositeIndex)[1])
{
break;
}
oppositeIndex++;
}
if (NO_SUCH_ELEMENT == next.Get(oppositeIndex))
{
succ.Set(i, first.Get(edges.Get(i)[1]));
}
else
{
succ.Set(i, next.Get(oppositeIndex));
}
}
}; // ComputeEulerTourList
} // process_contourtree_inc
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm
#endif

134
vtkm/worklet/contourtree_augmented/processcontourtree/ComputeMinMaxValues.h
View File

@ -1,134 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_compute_min_max_values_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_compute_min_max_values_h
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
/*
*
* This code is written by Petar Hristov in 09.2019
*
* This worklet computes a prefix min/max over the subtree of a rooted contour tree
* using the euler tour data structure.
*
* Given an euler tour and an array which has the index of the first and last occurence of every vertex,
* this worklet goes through that subarray of the euler tour to find the min/max value (in terms of regular node Id, or isovalue, does not matter which)
*
* This can be optimised with a vtkm parallel reduce operation
*
*/
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class ComputeMinMaxValues : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn supernodes,
WholeArrayIn firstLast,
WholeArrayIn tourEdges,
WholeArrayOut output);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _1;
bool isMin = true;
// Default Constructor
VTKM_EXEC_CONT ComputeMinMaxValues(bool _isMin)
: isMin(_isMin)
{
}
template <typename SupernodesArrayPortalType,
typename FirstLastArrayPortalType,
typename TourEdgesArrayPortalType,
typename OutputArrayPortalType>
VTKM_EXEC void operator()(const vtkm::Id i,
const SupernodesArrayPortalType& supernodes,
const FirstLastArrayPortalType& firstLast,
const TourEdgesArrayPortalType& tourEdges,
const OutputArrayPortalType& output) const
{
Id optimal = tourEdges.Get(firstLast.Get(i).first)[1];
for (Id j = firstLast.Get(i).first; j < firstLast.Get(i).second; j++)
{
Id vertex = tourEdges.Get(j)[1];
Id vertexValue = MaskedIndex(supernodes.Get(vertex));
Id optimalValue = MaskedIndex(supernodes.Get(optimal));
if ((true == isMin && vertexValue < optimalValue) ||
(false == isMin && vertexValue > optimalValue))
{
optimal = vertex;
}
}
output.Set(i, optimal);
}
}; // ComputeMinMaxValues
} // process_contourtree_inc
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm
#endif

180
vtkm/worklet/contourtree_augmented/processcontourtree/EulerTour.h
View File

@ -1,180 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_euler_tour_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_euler_tour_h
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourFirstNext.h>
#include <vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourList.h>
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class EulerTour
{
public:
vtkm::cont::Invoker Invoke;
// The sequence of the edges in the tours is given by succ array
cont::ArrayHandle<Id> succ;
// The tour consists of the directed edges of the contour tree
cont::ArrayHandle<Vec<Id, 2>> edges;
// The first occurance of an vertex
cont::ArrayHandle<Id> first;
// Compute the Euler Tour with root 0
void computeEulerTour(const IdArrayType::ReadPortalType superarcs)
{
//
// Make a list of all directed edges in the tree
//
edges.Allocate(2 * (superarcs.GetNumberOfValues() - 1));
auto edgesPortal = edges.WritePortal();
int counter = 0;
for (int i = 0; i < superarcs.GetNumberOfValues(); i++)
{
Id j = MaskedIndex(superarcs.Get(i));
if (j != 0)
{
edgesPortal.Set(counter++, { i, j });
edgesPortal.Set(counter++, { j, i });
}
}
vtkm::cont::Algorithm::Sort(edges, vtkm::SortLess());
//
// Initialize first and next arrays. They are used to compute the ciculer linked list that is the euler tour
//
//cont::ArrayHandle<Id> first;
cont::ArrayHandle<Id> next;
first = vtkm::cont::make_ArrayHandleMove(
std::vector<Id>(static_cast<unsigned long>(superarcs.GetNumberOfValues()), NO_SUCH_ELEMENT));
next = vtkm::cont::make_ArrayHandleMove(
std::vector<Id>(static_cast<unsigned long>(edges.GetNumberOfValues()), NO_SUCH_ELEMENT));
//
// Compute First and Next arrays that are needed to compute the euler tour linked list
//
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeEulerTourFirstNext
eulerWorklet;
this->Invoke(eulerWorklet, edges, first, next);
//
// Compute the euler tour as a circular linked list from the first and next arrays
//
succ.Allocate(edges.GetNumberOfValues());
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeEulerTourList
eulerTourListWorklet;
this->Invoke(eulerTourListWorklet, next, first, edges, succ);
}
// Reroot the euler tour at a different root (O(n) for finding the first occurence of the new root and O(1) for rerouting and O(n) for returning it as an array)
void getTourAtRoot(const Id root, const cont::ArrayHandle<Vec<Id, 2>>::WritePortalType tourEdges)
{
auto edgesPortal = edges.WritePortal();
//
// Reroot at the global min/max
//
Id i = 0;
Id start = NO_SUCH_ELEMENT;
do
{
if (edgesPortal.Get(i)[0] == root)
{
start = i;
break;
}
i = succ.WritePortal().Get(i);
} while (i != 0);
//
// Convert linked list to array handle array
//
{
int counter = 0;
i = start;
do
{
tourEdges.Set(counter++, edgesPortal.Get(i));
i = succ.WritePortal().Get(i);
} while (i != start);
}
}
}; // class EulerTour
} // namespace process_contourtree_inc
} // namespace contourtree_augmented
} // worklet
} // vtkm
#endif

718
vtkm/worklet/contourtree_augmented/processcontourtree/HypersweepWorklets.h Normal file
View File

@ -0,0 +1,718 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
// Copyright (c) 2018, The Regents of the University of California, through
// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
// from the U.S. Dept. of Energy). All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
// (1) Redistributions of source code must retain the above copyright notice, this
// list of conditions and the following disclaimer.
//
// (2) Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
//
// (3) Neither the name of the University of California, Lawrence Berkeley National
// Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
// IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
// INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
// OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
// OF THE POSSIBILITY OF SUCH DAMAGE.
//
//=============================================================================
//
// This code is an extension of the algorithm presented in the paper:
// Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
// Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
// Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
// (LDAV), October 2016, Baltimore, Maryland.
//
// The PPP2 algorithm and software were jointly developed by
// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_hypersweep_worklets_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_hypersweep_worklets_h
#include <vtkm/BinaryOperators.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
/**
* Incorporates values of the parent of the current subtree in the subtree for the min and max hypersweeps
*/
namespace vtkm
{
namespace worklet
{
namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class InitialiseArcsVolume : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _3;
vtkm::Id TotalVolume;
VTKM_EXEC_CONT InitialiseArcsVolume(vtkm::Id _totalVolume)
: TotalVolume(_totalVolume)
{
}
template <typename IdWholeArrayInPortalType, typename EdgeWholeArrayInOutPortal>
VTKM_EXEC void operator()(const vtkm::Id currentId,
const IdWholeArrayInPortalType& hypersweepSumValuesPortal,
const IdWholeArrayInPortalType& superarcIntrinsicWeightPortal,
const IdWholeArrayInPortalType& superarcsPortal,
const EdgeWholeArrayInOutPortal& arcsPortal) const
{
Id i = currentId;
Id parent = MaskedIndex(superarcsPortal.Get(i));
if (parent == 0)
{
// We expect the root to the last vertex in the supernodes array
VTKM_ASSERT(i != superarcsPortal.GetNumberOfValues() - 2);
return;
}
EdgeDataVolume edge;
edge.I = i;
edge.J = parent;
edge.UpEdge = IsAscending((superarcsPortal.Get(i)));
edge.SubtreeVolume = (this->TotalVolume - hypersweepSumValuesPortal.Get(i)) +
(superarcIntrinsicWeightPortal.Get(i) - 1);
EdgeDataVolume oppositeEdge;
oppositeEdge.I = parent;
oppositeEdge.J = i;
oppositeEdge.UpEdge = !edge.UpEdge;
oppositeEdge.SubtreeVolume = hypersweepSumValuesPortal.Get(i);
arcsPortal.Set(i * 2, edge);
arcsPortal.Set(i * 2 + 1, oppositeEdge);
}
};
class SetFirstVertexForSuperparent : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3);
using InputDomain = _1;
VTKM_EXEC_CONT SetFirstVertexForSuperparent() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(
const vtkm::Id sortedNode,
const IdWholeArrayInPortalType& nodesPortal,
const IdWholeArrayInPortalType& superparentsPortal,
const IdWholeArrayInOutPortalType& firstVertexForSuperparentPortal) const
{
vtkm::Id sortID = nodesPortal.Get(sortedNode);
vtkm::Id superparent = superparentsPortal.Get(sortID);
if (sortedNode == 0)
{
firstVertexForSuperparentPortal.Set(superparent, sortedNode);
}
else if (superparent != superparentsPortal.Get(nodesPortal.Get(sortedNode - 1)))
{
firstVertexForSuperparentPortal.Set(superparent, sortedNode);
}
}
};
class ComputeIntrinsicWeight : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _2;
VTKM_EXEC_CONT ComputeIntrinsicWeight() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id superarc,
const IdWholeArrayInPortalType& arcsPortal,
const IdWholeArrayInPortalType& superarcsPortal,
const IdWholeArrayInPortalType& firstVertexForSuperparentPortal,
const IdWholeArrayInOutPortalType& superarcIntrinsicWeightPortal) const
{
if (superarc == superarcsPortal.GetNumberOfValues() - 1)
{
superarcIntrinsicWeightPortal.Set(
superarc, arcsPortal.GetNumberOfValues() - firstVertexForSuperparentPortal.Get(superarc));
}
else
{
superarcIntrinsicWeightPortal.Set(superarc,
firstVertexForSuperparentPortal.Get(superarc + 1) -
firstVertexForSuperparentPortal.Get(superarc));
}
}
};
class SetFirstSupernodePerIteration : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2);
using InputDomain = _1;
VTKM_EXEC_CONT SetFirstSupernodePerIteration() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(
const vtkm::Id supernode,
const IdWholeArrayInPortalType& whenTransferredPortal,
const IdWholeArrayInOutPortalType& firstSupernodePerIterationPortal) const
{
vtkm::Id when = MaskedIndex(whenTransferredPortal.Get(supernode));
if (supernode == 0)
{
firstSupernodePerIterationPortal.Set(when, supernode);
}
else if (when != MaskedIndex(whenTransferredPortal.Get(supernode - 1)))
{
firstSupernodePerIterationPortal.Set(when, supernode);
}
}
};
template <typename Operator>
class AddDependentWeightHypersweep : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn iterationHypernodes,
WholeArrayIn hypernodes,
WholeArrayIn hyperarcs,
WholeArrayIn howManyUsed,
AtomicArrayInOut minMaxIndex);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4, _5);
using InputDomain = _1;
Operator Op;
// Default Constructor
VTKM_EXEC_CONT AddDependentWeightHypersweep(Operator _op)
: Op(_op)
{
}
template <typename IdWholeArrayHandleCountingIn,
typename IdWholeArrayIn,
typename IdWholeArrayInOut>
VTKM_EXEC void operator()(const vtkm::Id hyperarcId,
const IdWholeArrayHandleCountingIn& iterationHypernodesPortal,
const IdWholeArrayIn& hypernodesPortal,
const IdWholeArrayIn& hyperarcsPortal,
const IdWholeArrayIn& howManyUsedPortal,
const IdWholeArrayInOut& minMaxIndexPortal) const
{
Id i = iterationHypernodesPortal.Get(hyperarcId);
// If it's the last hyperacs (there's nothing to do it's just the root)
if (i >= hypernodesPortal.GetNumberOfValues() - 1)
{
return;
}
//
// The value of the prefix scan is now accumulated in the last supernode of the hyperarc. Transfer is to the target
//
vtkm::Id lastSupernode = MaskedIndex(hypernodesPortal.Get(i + 1)) - howManyUsedPortal.Get(i);
//
// Transfer the accumulated value to the target supernode
//
Id vertex = lastSupernode - 1;
Id parent = MaskedIndex(hyperarcsPortal.Get(i));
Id vertexValue = minMaxIndexPortal.Get(vertex);
//Id parentValue = minMaxIndexPortal.Get(parent);
//Id writeValue = op(vertexValue, parentValue);
auto cur = minMaxIndexPortal.Get(parent); // Load the current value at idx
vtkm::Id newVal; // will hold the result of the multiplication
vtkm::Id expect; // will hold the expected value before multiplication
do
{
expect = cur; // Used to ensure the value hasn't changed since reading
newVal = this->Op(cur, vertexValue); // the actual multiplication
} while ((cur = minMaxIndexPortal.CompareAndSwap(parent, newVal, expect)) != expect);
//minMaxIndexPortal.Set(parent, writeValue);
}
}; // ComputeMinMaxValues
class InitialiseArcs : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4, _5, _6);
using InputDomain = _1;
vtkm::Id GlobalMinSortedIndex, GlobalMaxSortedIndex, RootSupernodeId;
VTKM_EXEC_CONT InitialiseArcs(vtkm::Id _globalMinSortedIndex,
vtkm::Id _globalMaxSortedIndex,
vtkm::Id _rootSupernodeId)
: GlobalMinSortedIndex(_globalMinSortedIndex)
, GlobalMaxSortedIndex(_globalMaxSortedIndex)
, RootSupernodeId(_rootSupernodeId)
{
}
template <typename IdWholeArrayInPortalType, typename EdgeWholeArrayInOutPortal>
VTKM_EXEC void operator()(const vtkm::Id currentId,
const IdWholeArrayInPortalType& minParentsPortal,
const IdWholeArrayInPortalType& maxParentsPortal,
const IdWholeArrayInPortalType& minValuesPortal,
const IdWholeArrayInPortalType& maxValuesPortal,
const IdWholeArrayInPortalType& superarcsPortal,
const EdgeWholeArrayInOutPortal& arcsPortal) const
{
Id i = currentId;
Id parent = MaskedIndex(superarcsPortal.Get(i));
// The root does not correspond to an arc
if (parent == 0)
return;
EdgeDataHeight edge;
edge.I = i;
edge.J = parent;
edge.UpEdge = IsAscending((superarcsPortal.Get(i)));
EdgeDataHeight oppositeEdge;
oppositeEdge.I = parent;
oppositeEdge.J = i;
oppositeEdge.UpEdge = !edge.UpEdge;
// Is it in the direction of the minRootedTree?
if (MaskedIndex(minParentsPortal.Get(edge.J)) == edge.I)
{
edge.SubtreeMin = minValuesPortal.Get(edge.J);
oppositeEdge.SubtreeMin = this->GlobalMinSortedIndex;
}
else
{
oppositeEdge.SubtreeMin = minValuesPortal.Get(oppositeEdge.J);
edge.SubtreeMin = this->GlobalMinSortedIndex;
}
// Is it in the direction of the maxRootedTree?
if (MaskedIndex(maxParentsPortal.Get(edge.J)) == edge.I)
{
edge.SubtreeMax = maxValuesPortal.Get(edge.J);
oppositeEdge.SubtreeMax = this->GlobalMaxSortedIndex;
}
else
{
oppositeEdge.SubtreeMax = maxValuesPortal.Get(oppositeEdge.J);
edge.SubtreeMax = this->GlobalMaxSortedIndex;
}
// We technically don't need this because the root is supposed to be the last vertex
if (i > this->RootSupernodeId)
{
i--;
}
arcsPortal.Set(i * 2, edge);
arcsPortal.Set(i * 2 + 1, oppositeEdge);
}
};
class ComputeSubtreeHeight : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _4;
VTKM_EXEC_CONT ComputeSubtreeHeight() {}
template <typename Float64WholeArrayInPortalType,
typename IdWholeArrayInPortalType,
typename EdgeWholeArrayInOutPortal>
VTKM_EXEC void operator()(const vtkm::Id currentId,
const Float64WholeArrayInPortalType& fieldValuesPortal,
const IdWholeArrayInPortalType& ctSortOrderPortal,
const IdWholeArrayInPortalType& supernodesPortal,
const EdgeWholeArrayInOutPortal& arcsPortal) const
{
Id i = currentId;
EdgeDataHeight edge = arcsPortal.Get(i);
Float64 minIsoval = fieldValuesPortal.Get(ctSortOrderPortal.Get(edge.SubtreeMin));
Float64 maxIsoval = fieldValuesPortal.Get(ctSortOrderPortal.Get(edge.SubtreeMax));
Float64 vertexIsoval =
fieldValuesPortal.Get(ctSortOrderPortal.Get(supernodesPortal.Get(edge.I)));
// We need to incorporate the value of the vertex into the height of the tree (otherwise leafs edges have 0 persistence)
minIsoval = vtkm::Minimum()(minIsoval, vertexIsoval);
maxIsoval = vtkm::Maximum()(maxIsoval, vertexIsoval);
edge.SubtreeHeight = maxIsoval - minIsoval;
arcsPortal.Set(i, edge);
}
}; // ComputeMinMaxValues
class SetBestUpDown : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayInOut, WholeArrayInOut, WholeArrayIn);
typedef void ExecutionSignature(InputIndex, _1, _2, _3);
using InputDomain = _3;
VTKM_EXEC_CONT SetBestUpDown() {}
template <typename IdWholeArrayInPortalType, typename EdgeWholeArrayInOutPortal>
VTKM_EXEC void operator()(const vtkm::Id currentId,
const IdWholeArrayInPortalType& bestUpwardPortal,
const IdWholeArrayInPortalType& bestDownwardPortal,
const EdgeWholeArrayInOutPortal& arcsPortal) const
{
vtkm::Id i = currentId;
if (i == 0)
{
if (arcsPortal.Get(0).UpEdge == 0)
{
bestDownwardPortal.Set(arcsPortal.Get(0).I, arcsPortal.Get(0).J);
}
else
{
bestUpwardPortal.Set(arcsPortal.Get(0).I, arcsPortal.Get(0).J);
}
}
else
{
if (arcsPortal.Get(i).UpEdge == 0 && arcsPortal.Get(i).I != arcsPortal.Get(i - 1).I)
{
bestDownwardPortal.Set(arcsPortal.Get(i).I, arcsPortal.Get(i).J);
}
if (arcsPortal.Get(i).UpEdge == 1 &&
(arcsPortal.Get(i).I != arcsPortal.Get(i - 1).I || arcsPortal.Get(i - 1).UpEdge == 0))
{
bestUpwardPortal.Set(arcsPortal.Get(i).I, arcsPortal.Get(i).J);
}
}
}
}; // ComputeMinMaxValues
class UnmaskArray : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1);
using InputDomain = _1;
// Default Constructor
VTKM_EXEC_CONT UnmaskArray() {}
template <typename IdWholeArrayInPortalType>
VTKM_EXEC void operator()(const vtkm::Id currentId,
const IdWholeArrayInPortalType& maskedArrayPortal) const
{
const auto currentValue = maskedArrayPortal.Get(currentId);
maskedArrayPortal.Set(currentId, MaskedIndex(currentValue));
}
}; // ComputeMinMaxValues
class PropagateBestUpDown : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3);
using InputDomain = _3;
// Default Constructor
VTKM_EXEC_CONT PropagateBestUpDown() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernodeId,
const IdWholeArrayInPortalType& bestUpwardPortal,
const IdWholeArrayInPortalType& bestDownwardPortal,
const IdWholeArrayOutPortalType& whichBranchPortal) const
{
vtkm::Id bestUp = bestUpwardPortal.Get(supernodeId);
if (NoSuchElement(bestUp))
{
// flag it as an upper leaf
whichBranchPortal.Set(supernodeId, TERMINAL_ELEMENT | supernodeId);
}
else if (bestDownwardPortal.Get(bestUp) == supernodeId)
whichBranchPortal.Set(supernodeId, bestUp);
else
whichBranchPortal.Set(supernodeId, TERMINAL_ELEMENT | supernodeId);
}
}; // ComputeMinMaxValues
class WhichBranchNewId : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2);
using InputDomain = _2;
// Default Constructor
VTKM_EXEC_CONT WhichBranchNewId() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernode,
const IdWholeArrayInPortalType& chainToBranchPortal,
const IdWholeArrayInOutPortalType& whichBranchPortal) const
{
const auto currentValue = MaskedIndex(whichBranchPortal.Get(supernode));
whichBranchPortal.Set(supernode, chainToBranchPortal.Get(currentValue));
}
}; // ComputeMinMaxValues
class BranchMinMaxSet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayInOut, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4);
using InputDomain = _2;
vtkm::Id NumSupernodes;
// Default Constructor
VTKM_EXEC_CONT BranchMinMaxSet(vtkm::Id _NumSupernodes)
: NumSupernodes(_NumSupernodes)
{
}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernode,
const IdWholeArrayInPortalType& supernodeSorterPortal,
const IdWholeArrayInPortalType& whichBranchPortal,
const IdWholeArrayInOutPortalType& branchMinimumPortal,
const IdWholeArrayInOutPortalType& branchMaximumPortal) const
{
// retrieve supernode & branch IDs
vtkm::Id supernodeID = supernodeSorterPortal.Get(supernode);
vtkm::Id branchID = whichBranchPortal.Get(supernodeID);
// save the branch ID as the owner
// use LHE of segment to set branch minimum
if (supernode == 0)
{ // sn = 0
branchMinimumPortal.Set(branchID, supernodeID);
} // sn = 0
else if (branchID != whichBranchPortal.Get(supernodeSorterPortal.Get(supernode - 1)))
{ // LHE
branchMinimumPortal.Set(branchID, supernodeID);
} // LHE
// use RHE of segment to set branch maximum
if (supernode == this->NumSupernodes - 1)
{ // sn = max
branchMaximumPortal.Set(branchID, supernodeID);
} // sn = max
else if (branchID != whichBranchPortal.Get(supernodeSorterPortal.Get(supernode + 1)))
{ // RHE
branchMaximumPortal.Set(branchID, supernodeID);
} // RHE
}
}; // ComputeMinMaxValues
class BranchSaddleParentSet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayIn,
WholeArrayInOut,
WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3, _4, _5, _6, _7);
using InputDomain = _2;
// Default Constructor
VTKM_EXEC_CONT BranchSaddleParentSet() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id branchID,
const IdWholeArrayInPortalType& whichBranchPortal,
const IdWholeArrayInPortalType& branchMinimumPortal,
const IdWholeArrayInPortalType& branchMaximumPortal,
const IdWholeArrayInPortalType& bestDownwardPortal,
const IdWholeArrayInPortalType& bestUpwardPortal,
const IdWholeArrayInOutPortalType& branchSaddlePortal,
const IdWholeArrayInOutPortalType& branchParentPortal) const
{
vtkm::Id branchMax = branchMaximumPortal.Get(branchID);
// check whether the maximum is NOT a leaf
if (!NoSuchElement(bestUpwardPortal.Get(branchMax)))
{ // points to a saddle
branchSaddlePortal.Set(branchID, MaskedIndex(bestUpwardPortal.Get(branchMax)));
// if not, then the bestUp points to a saddle vertex at which we join the parent
branchParentPortal.Set(branchID, whichBranchPortal.Get(bestUpwardPortal.Get(branchMax)));
} // points to a saddle
// now do the same with the branch minimum
vtkm::Id branchMin = branchMinimumPortal.Get(branchID);
// test whether NOT a lower leaf
if (!NoSuchElement(bestDownwardPortal.Get(branchMin)))
{ // points to a saddle
branchSaddlePortal.Set(branchID, MaskedIndex(bestDownwardPortal.Get(branchMin)));
// if not, then the bestUp points to a saddle vertex at which we join the parent
branchParentPortal.Set(branchID, whichBranchPortal.Get(bestDownwardPortal.Get(branchMin)));
} // points to a saddle
}
}; // ComputeMinMaxValues
class PrepareChainToBranch : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2);
using InputDomain = _1;
// Default Constructor
VTKM_EXEC_CONT PrepareChainToBranch() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernode,
const IdWholeArrayInPortalType& whichBranchPortal,
const IdWholeArrayInOutPortalType& chainToBranchPortal) const
{
// test whether the supernode points to itself to find the top ends
if (MaskedIndex(whichBranchPortal.Get(supernode)) == supernode)
{
chainToBranchPortal.Set(supernode, 1);
}
}
}; // ComputeMinMaxValues
class FinaliseChainToBranch : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2);
using InputDomain = _1;
VTKM_EXEC_CONT FinaliseChainToBranch() {}
template <typename IdWholeArrayInPortalType, typename IdWholeArrayInOutPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernode,
const IdWholeArrayInPortalType& whichBranchPortal,
const IdWholeArrayInOutPortalType& chainToBranchPortal) const
{
// test whether the supernode points to itself to find the top ends
if (MaskedIndex(whichBranchPortal.Get(supernode)) == supernode)
{
const auto value = chainToBranchPortal.Get(supernode);
chainToBranchPortal.Set(supernode, value - 1);
}
else
{
// This should just be vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT, but there's the build error - identifier "vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT" is undefined in device code
chainToBranchPortal.Set(supernode, std::numeric_limits<vtkm::Id>::min());
}
}
}; // ComputeMinMaxValues
template <typename Operator>
class IncorporateParent : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn, WholeArrayIn, WholeArrayInOut);
typedef void ExecutionSignature(InputIndex, _1, _2, _3);
using InputDomain = _1;
Operator Op;
// Default Constructor
VTKM_EXEC_CONT IncorporateParent(Operator _op)
: Op(_op)
{
}
template <typename IdWholeArrayIn, typename IdWholeArrayInOut>
VTKM_EXEC void operator()(const vtkm::Id superarcId,
const IdWholeArrayIn& parentsPortal,
const IdWholeArrayIn& supernodesPortal,
const IdWholeArrayInOut& hypersweepValuesPortal) const
{
Id i = superarcId;
Id parent = MaskedIndex(parentsPortal.Get(i));
Id subtreeValue = hypersweepValuesPortal.Get(i);
Id parentValue = MaskedIndex(supernodesPortal.Get(parent));
hypersweepValuesPortal.Set(i, this->Op(subtreeValue, parentValue));
}
}; // ComputeMinMaxValues
} // process_contourtree_inc
} // namespace contourtree_augmented
} // namespace worklet
} // namespace vtkm
#endif

54
vtkm/worklet/contourtree_augmented/processcontourtree/ComputeEulerTourFirstNext.h → vtkm/worklet/contourtree_augmented/processcontourtree/PointerDoubling.h
View File

@ -50,18 +50,18 @@
// Oliver Ruebel (LBNL)
//==============================================================================
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_euler_tour_compute_first_next_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_euler_tour_compute_first_next_h
#ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_inc_pointer_doubling_h
#define vtk_m_worklet_contourtree_augmented_process_contourtree_inc_pointer_doubling_h
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/contourtree_augmented/Types.h>
/*
* This code is written by Petar Hristov in 09.2019
*
* This code is written by Petar Hristov in 03.2020
*
* It does pointer doubling on an array.
*
* This worklet performs the first step in computing an Euler tour.
* That is computing two auxiliary array, from which we can compute the linked edge list of the Euler tour.
* See https://en.wikipedia.org/wiki/Euler_tour_technique for details
*
*/
namespace vtkm
@ -72,40 +72,34 @@ namespace contourtree_augmented
{
namespace process_contourtree_inc
{
class ComputeEulerTourFirstNext : public vtkm::worklet::WorkletMapField
class PointerDoubling : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(WholeArrayIn edges, WholeArrayOut first, WholeArrayOut next);
typedef void ControlSignature(WholeArrayInOut whichBranch);
typedef void ExecutionSignature(InputIndex, _1, _2, _3);
typedef void ExecutionSignature(InputIndex, _1);
using InputDomain = _1;
// Default Constructor
VTKM_EXEC_CONT ComputeEulerTourFirstNext() {}
vtkm::Id NumSupernodes;
template <typename EdgesArrayPortalType, typename OutputArrayPortalType>
VTKM_EXEC void operator()(const vtkm::Id i,
const EdgesArrayPortalType& edges,
const OutputArrayPortalType& first,
const OutputArrayPortalType& next) const
// Default Constructor
VTKM_EXEC_CONT PointerDoubling(vtkm::Id _NumSupernodes)
: NumSupernodes(_NumSupernodes)
{
if (i == 0)
}
template <typename WhichBranchArrayPortalType>
VTKM_EXEC void operator()(const vtkm::Id supernode,
const WhichBranchArrayPortalType& whichBranchPortal) const
{
if (!IsTerminalElement(whichBranchPortal.Get(supernode)))
{
first.Set(0, 0);
}
else
{
if (edges.Get(i)[0] != edges.Get(i - 1)[0])
{
first.Set(edges.Get(i)[0], i);
}
else
{
next.Set(i - 1, i);
}
const auto currentValue =
MaskedIndex(whichBranchPortal.Get(MaskedIndex(whichBranchPortal.Get(supernode))));
whichBranchPortal.Set(supernode, currentValue);
}
}
}; // ComputeEulerTourFirstNext
}; // ComputeMinMaxValues
} // process_contourtree_inc
} // namespace contourtree_augmented
} // namespace worklet