//============================================================================ // 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 2014 Sandia Corporation. // Copyright 2014 UT-Battelle, LLC. // Copyright 2014 Los Alamos National Security. // // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, // the U.S. Government retains certain rights in this software. // // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National // Laboratory (LANL), the U.S. Government retains certain rights in // this software. //============================================================================ //======================================================================================= // // Second Attempt to Compute Contour Tree in Data-Parallel Mode // // Started August 19, 2015 // // Copyright Hamish Carr, University of Leeds // // ContourTree.h - class representing the contour tree // //======================================================================================= // // COMMENTS: // // // i.e. based on PeakPitPruningCriticalSerial // // Under the old merge approach, we had an essentially breadth-first queue for transferring // leaves from the merge trees to the contour tree. // // Most of these leaves are completely independent of each other, and can (on principle) // be processed simultaneously. However, the interior of the tree is dependent on them // having been dealt with already. This version, therefore, will make multiple passes, // in each pass pruning all maxima then all minima, interspersed with updating the merge // and split trees. To understand this, consider what happens in the merge algorithm when // a maximum is added: // // 1. The vertex v is removed from the queue: it has one join neighbour, w // 2. Edge (v,w) is removed from the join tree, along with vertex v // 3. Edge (v,w) is added to the contour tree, with v, w if necessary // 4. Vertex v is removed from the split tree, bridging edges past it if necessary // 5. Vertex w is added to the queue iff it is now a leaf // // To parallelise this: // For all vertices v // Set contourArc[v] = NO_VERTEX_ASSIGNED // Set nContourArcs = 0; // While (nContourArcs) > 0 // might be one, or something else - base case isn't clear // a. Use reduction to compute updegree from join tree, downdegree from split tree // b. For each vertex v // // omit previously processed vertices // if (contourArc[v] == NO_VERTEX_ASSIGNED) // continue; // // Test for extremality // i. If ((updegree[v] == 0) && (downdegree[v] == 1)) // { // Maximum // contourArc[v] = joinArc[v]; // } // Maximum // ii. Else if ((updegree[v] = 1) && (downdegree[v] == 0)) // { // Minimum // contourArc[v] = splitArc[v]; // } // Minimum // c. For (log n iterations) // i. For each vertex v // retrieve it's join neighbour j // retrieve it's split neighbour s // if v has no join neighbour (i.e. j == -1) // skip (i.e. v is the root) // else if j has a contour arc assigned // set v's neighbour to j's neighbour // if v has no split neighbour (i.e. s == -1) // skip (i.e. v is the root) // else if s has a contour arc assigned // set v's neighbour to s's neighbour // // Initially, we will do this with all vertices, regular or otherwise, then restrict to // the critical points. Number of iterations - regular vertices will slow this down, so // the worst case is O(n) passes. Even if we restrict to critical points, W's in the tree // will serialise, so O(n) still applies. I believe that the W edges can be suppressed, // but let's leave that to optimisation for now. // //======================================================================================= #ifndef vtkm_worklet_contourtree_contourtree_h #define vtkm_worklet_contourtree_contourtree_h // local includes #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define DEBUG_PRINT 1 //#define DEBUG_TIMING 1 namespace vtkm { namespace worklet { namespace contourtree { template class ContourTree { public: typedef typename vtkm::cont::DeviceAdapterAlgorithm DeviceAlgorithm; typedef vtkm::cont::ArrayHandle IdArrayType; typedef vtkm::cont::ArrayHandle ValueArrayType; typedef vtkm::cont::ArrayHandlePermutation PermuteIndexType; typedef vtkm::cont::ArrayHandlePermutation PermuteValueType; // device DeviceAdapter device; // reference to the underlying data const vtkm::cont::ArrayHandle values; // vector of superarcs in the contour tree (stored as inward-pointing) vtkm::cont::ArrayHandle superarcs; // vector of supernodes vtkm::cont::ArrayHandle supernodes; // vector of supernodes still unprocessed vtkm::cont::ArrayHandle activeSupernodes; // references to join & split trees MergeTree &joinTree, &splitTree; // references to join & split graphs ChainGraph &joinGraph, &splitGraph; // vectors of up & down degree used during computation vtkm::cont::ArrayHandle updegree, downdegree; // vectors for tracking merge arcs vtkm::cont::ArrayHandle joinArcs, splitArcs; // counter for how many iterations it took to compute vtkm::Id nIterations; // contour tree constructor ContourTree(const vtkm::cont::ArrayHandle &Values, DeviceAdapter Device, MergeTree &JoinTree, ChainGraph &JoinGraph, MergeTree &SplitTree, ChainGraph &SplitGraph); // routines for computing the contour tree // combines the list of active vertices for join & split trees // then reduces them to eliminate regular vertices & non-connectivity critical points void FindSupernodes(); // transfers leaves from join/split trees to contour tree void TransferLeaves(); // collapses regular edges along leaf superarcs void CollapseRegular(bool isJoin); // compresses trees to remove transferred vertices void CompressTrees(); // compresses active set of supernodes void CompressActiveSupernodes(); // finds the degree of each supernode from the join & split trees void FindDegrees(); // collect the resulting saddle peaks in sort pairs void CollectSaddlePeak(vtkm::cont::ArrayHandle > &saddlePeak); void DebugPrint(const char *message); }; // class ContourTree struct VertexAssigned : vtkm::worklet::WorkletMapField { public: typedef void ControlSignature(FieldIn supernode, WholeArrayIn superarcs, FieldOut hasSuperArc); typedef _3 ExecutionSignature(_1, _2); typedef _1 InputDomain; bool vertexIsAssigned; VTKM_EXEC_CONT VertexAssigned(bool VertexIsAssigned) : vertexIsAssigned(VertexIsAssigned) {} template VTKM_EXEC vtkm::Id operator()(const vtkm::Id supernode, const InPortalFieldType& superarcs) const { if (vertexIsAssigned == false) { if (superarcs.Get(supernode) == NO_VERTEX_ASSIGNED) return vtkm::Id(1); else return vtkm::Id(0); } else { if (superarcs.Get(supernode) != NO_VERTEX_ASSIGNED) return vtkm::Id(1); else return vtkm::Id(0); } } }; // creates contour tree template ContourTree::ContourTree( const vtkm::cont::ArrayHandle &Values, DeviceAdapter Device, MergeTree &JoinTree, ChainGraph &JoinGraph, MergeTree &SplitTree, ChainGraph &SplitGraph) : values(Values), device(Device), joinTree(JoinTree), joinGraph(JoinGraph), splitTree(SplitTree), splitGraph(SplitGraph) { // first we have to get the correct list of supernodes // this will also set the degrees of the vertices initially FindSupernodes(); // set the active vertices to something greater than 0 for the loop vtkm::Id nActiveVertices = supernodes.GetNumberOfValues(); // and track how many iterations it takes nIterations = 0; // loop until no arcs remaining to be found // tree can end with either 0 or 1 vertices unprocessed // 0 means the last edge was pruned from both ends // 1 means that there were two final edges meeting at a vertex while (activeSupernodes.GetNumberOfValues() > 1) { // loop until no active vertices remaining #ifdef DEBUG_PRINT std::cout << "========================================" << std::endl; std::cout << " " << std::endl; std::cout << "Iteration " << nIterations << " Size " << activeSupernodes.GetNumberOfValues() << std::endl; std::cout << " " << std::endl; std::cout << "========================================" << std::endl; #endif // transfer all leaves to the contour tree TransferLeaves(); // collapse regular vertices from leaves, upper then lower CollapseRegular(true); CollapseRegular(false); // compress the join and split trees CompressTrees(); // compress the active list of supernodes CompressActiveSupernodes(); // recompute the vertex degrees FindDegrees(); nIterations++; } } // constructor // combines the list of active vertices for join & split trees // then reduces them to eliminate regular vertices & non-connectivity critical points template void ContourTree::FindSupernodes() { // both trees may have non-connectivity critical points, so we first make a joint list // here, we will explicitly assume that the active lists are in numerical order // which is how we are currently constructing them vtkm::Id nCandidates = joinGraph.valueIndex.GetNumberOfValues() + splitGraph.valueIndex.GetNumberOfValues(); vtkm::cont::ArrayHandle candidates; // take the union of the two sets of vertices vtkm::cont::ArrayHandleConcatenate candidateArray(joinGraph.valueIndex, splitGraph.valueIndex); DeviceAlgorithm::Copy(candidateArray, candidates); DeviceAlgorithm::Sort(candidates); DeviceAlgorithm::Unique(candidates); nCandidates = candidates.GetNumberOfValues(); vtkm::cont::ArrayHandleIndex candidateIndexArray(nCandidates); // we need an array lookup to convert vertex ID's vtkm::Id nValues = values.GetNumberOfValues(); vtkm::cont::ArrayHandle regularToCritical; vtkm::cont::ArrayHandleConstant noVertArray(NO_VERTEX_ASSIGNED, nValues); DeviceAlgorithm::Copy(noVertArray, regularToCritical); if (nCandidates > 0) { RegularToCriticalUp regularToCriticalUp; vtkm::worklet::DispatcherMapField regularToCriticalUpDispatcher(regularToCriticalUp); regularToCriticalUpDispatcher.Invoke(candidateIndexArray, // input candidates, // input regularToCritical); // output (whole array) } // now that we have a complete list of active nodes from each, we can call the trees // to connect them properly joinTree.ComputeAugmentedSuperarcs(); joinTree.ComputeAugmentedArcs(candidates); splitTree.ComputeAugmentedSuperarcs(); splitTree.ComputeAugmentedArcs(candidates); // we create up & down degree arrays vtkm::cont::ArrayHandleConstant initCandidateArray(0, nCandidates); vtkm::cont::ArrayHandle upCandidate; vtkm::cont::ArrayHandle downCandidate; DeviceAlgorithm::Copy(initCandidateArray, upCandidate); DeviceAlgorithm::Copy(initCandidateArray, downCandidate); // This next chunk changes in parallel - it has to count the up & down degree for each // vertex. It's a simple loop in serial, but in parallel, what we have to do is: // 1. Copy the lower ends of the edges, converting from regular ID to candidate ID // 2. Sort the lower ends of the edges // 3. For each value, store the beginning of the range // 4. Compute the delta to get the degree. // create a sorting vector vtkm::cont::ArrayHandle sortVector; sortVector.Allocate(nCandidates); // 1. Copy the lower ends of the edges, converting from regular ID to candidate ID if (nCandidates > 0) { RegularToCandidate regularToCandidate; vtkm::worklet::DispatcherMapField regularToCandidateDispatcher(regularToCandidate); regularToCandidateDispatcher.Invoke(candidates, // input joinTree.mergeArcs, // input (whole array) regularToCritical, // input (whole array) sortVector); // output } // 2. Sort the lower ends of the edges DeviceAlgorithm::Sort(sortVector); // 3. For each value, store the beginning & end of the range (in parallel) // The 0th element is guaranteed to be NO_VERTEX_ASSIGNED, & can be skipped // Otherwise, if the i-1th element is different, we are the offset for the subrange // and store into the ith place vtkm::cont::ArrayHandleCounting subsetIndexArray(1, 1, nCandidates-1); if (nCandidates > 0) { SubrangeOffset subRangeOffset; vtkm::worklet::DispatcherMapField subrangeOffsetDispatcher(subRangeOffset); subrangeOffsetDispatcher.Invoke(subsetIndexArray, // index sortVector, // input upCandidate); // output } // 4. Compute the delta to get the degree. if (nCandidates > 0) { DegreeDelta degreeDelta(nCandidates); vtkm::worklet::DispatcherMapField degreeDeltaDispatcher(degreeDelta); degreeDeltaDispatcher.Invoke(subsetIndexArray, // input sortVector, // input (whole array) upCandidate); // output (whole array) } // Now repeat the same steps for the downdegree // 1. Copy the upper ends of the edges, converting from regular ID to candidate ID if (nCandidates > 0) { RegularToCriticalDown regularToCriticalDown; vtkm::worklet::DispatcherMapField regularToCriticalDownDispatcher(regularToCriticalDown); regularToCriticalDownDispatcher.Invoke(candidates, // input splitTree.mergeArcs, // input (whole array) regularToCritical, // input (whole array) sortVector); // output } // 2. Sort the lower ends of the edges DeviceAlgorithm::Sort(sortVector); // 3. For each value, store the beginning & end of the range (in parallel) // The 0th element is guaranteed to be NO_VERTEX_ASSIGNED, & can be skipped // Otherwise, if the i-1th element is different, we are the offset for the subrange // and store into the ith place if (nCandidates > 0) { SubrangeOffset subRangeOffset; vtkm::worklet::DispatcherMapField subrangeOffsetDispatcher(subRangeOffset); subrangeOffsetDispatcher.Invoke(subsetIndexArray, // index sortVector, // input downCandidate); // output } // 4. Compute the delta to get the degree. if (nCandidates > 0) { DegreeDelta degreeDelta(nCandidates); vtkm::worklet::DispatcherMapField degreeDeltaDispatcher(degreeDelta); degreeDeltaDispatcher.Invoke(subsetIndexArray, // index sortVector, // input downCandidate); // in out } // create an index vector for whether the vertex is to be kept vtkm::cont::ArrayHandle isSupernode; isSupernode.Allocate(nCandidates); // fill the vector in if (nCandidates > 0) { FillSupernodes fillSupernodes; vtkm::worklet::DispatcherMapField fillSupernodesDispatcher(fillSupernodes); fillSupernodesDispatcher.Invoke(upCandidate, // input downCandidate, // input isSupernode); // output } // do a compaction to find the new index for each // We end with 0 in position 0, and need one extra position to find the new size vtkm::cont::ArrayHandle supernodeID; DeviceAlgorithm::ScanExclusive(isSupernode, supernodeID); // size is the position of the last element + the size of the last element (0/1) vtkm::Id nSupernodes = supernodeID.GetPortalConstControl().Get(nCandidates-1) + isSupernode.GetPortalConstControl().Get(nCandidates-1); // allocate memory for our arrays supernodes.ReleaseResources(); updegree.ReleaseResources(); downdegree.ReleaseResources(); supernodes.Allocate(nSupernodes); updegree.Allocate(nSupernodes); downdegree.Allocate(nSupernodes); // now copy over the positions to compact if (nCandidates > 0) { CopySupernodes copySupernodes; vtkm::worklet::DispatcherMapField copySupernodesDispatcher(copySupernodes); copySupernodesDispatcher.Invoke(isSupernode, // input candidates, // input supernodeID, // input upCandidate, // input downCandidate, // input regularToCritical, // output (whole array) supernodes, // output (whole array) updegree, // output (whole array) downdegree); // output (whole array) } // now we call the merge tree again to reset the merge arcs joinTree.ComputeAugmentedArcs(supernodes); splitTree.ComputeAugmentedArcs(supernodes); // next we create the working arrays of merge arcs nSupernodes = supernodes.GetNumberOfValues(); vtkm::cont::ArrayHandleIndex supernodeIndexArray(nSupernodes); joinArcs.ReleaseResources(); splitArcs.ReleaseResources(); joinArcs.Allocate(nSupernodes); splitArcs.Allocate(nSupernodes); // and copy them across, setting IDs for both ends SetJoinAndSplitArcs setJoinAndSplitArcs; vtkm::worklet::DispatcherMapField setJoinAndSplitArcsDispatcher(setJoinAndSplitArcs); setJoinAndSplitArcsDispatcher.Invoke(supernodes, // input joinTree.mergeArcs, // input (whole array) splitTree.mergeArcs, // input (whole array) regularToCritical, // input (whole array) joinArcs, // output splitArcs); // output vtkm::cont::ArrayHandleConstant newsuperarcs(NO_VERTEX_ASSIGNED, nSupernodes); superarcs.ReleaseResources(); DeviceAlgorithm::Copy(newsuperarcs, superarcs); // create the active supernode vector activeSupernodes.ReleaseResources(); activeSupernodes.Allocate(nSupernodes); vtkm::cont::ArrayHandleIndex supernodeSeq(nSupernodes); DeviceAlgorithm::Copy(supernodeSeq, activeSupernodes); #ifdef DEBUG_PRINT DebugPrint("Supernodes Found"); #endif } // FindSupernodes() // transfers leaves from join/split trees to contour tree template void ContourTree::TransferLeaves() { FindLeaves findLeaves; vtkm::worklet::DispatcherMapField findLeavesDispatcher(findLeaves); findLeavesDispatcher.Invoke(activeSupernodes, // input updegree, // input (whole array) downdegree, // input (whole array) joinArcs, // input (whole array) splitArcs, // input (whole array) superarcs); // i/o (whole array) #ifdef DEBUG_PRINT DebugPrint("Leaves Transferred"); #endif } // TransferLeaves() // collapses regular edges along leaf superarcs template void ContourTree::CollapseRegular(bool isJoin) { // we'll have a vector for tracking outwards vtkm::Id nSupernodes = supernodes.GetNumberOfValues(); vtkm::cont::ArrayHandleConstant nullArray(0, nSupernodes); vtkm::cont::ArrayHandle outbound; outbound.Allocate(nSupernodes); DeviceAlgorithm::Copy(nullArray, outbound); // and a reference for the inwards array and to the degrees vtkm::cont::ArrayHandle inbound; vtkm::cont::ArrayHandle indegree; vtkm::cont::ArrayHandle outdegree; if (isJoin) { DeviceAlgorithm::Copy(joinArcs, inbound); DeviceAlgorithm::Copy(downdegree, indegree); DeviceAlgorithm::Copy(updegree, outdegree); } else { DeviceAlgorithm::Copy(splitArcs, inbound); DeviceAlgorithm::Copy(updegree, indegree); DeviceAlgorithm::Copy(downdegree, outdegree); } // loop to copy join/split CopyJoinSplit copyJoinSplit; vtkm::worklet::DispatcherMapField copyJoinSplitDispatcher(copyJoinSplit); copyJoinSplitDispatcher.Invoke(activeSupernodes, // input inbound, // input (whole array) indegree, // input (whole array) outdegree, // input (whole array) outbound); // output (whole array) // Compute the number of log steps required in this pass vtkm::Id nLogSteps = 1; vtkm::Id nActiveSupernodes = activeSupernodes.GetNumberOfValues(); for (vtkm::Id shifter = nActiveSupernodes; shifter != 0; shifter >>= 1) nLogSteps++; // loop to find the now-regular vertices and collapse past them without altering // the existing join & split arcs for (vtkm::Id iteration = 0; iteration < nLogSteps; iteration++) { UpdateOutbound updateOutbound; vtkm::worklet::DispatcherMapField updateOutboundDispatcher(updateOutbound); updateOutboundDispatcher.Invoke(activeSupernodes, // input outbound); // i/o (whole array) } // at this point, the outbound vector chains everything outwards to the leaf // any vertices on the last outbound leaf superarc point to the leaf // Now, any regular leaf vertex points out to a leaf, so the condition we test is // a. outbound is not -1 (i.e. vertex is regular) // b. superarc[outbound] is not -1 (i.e. outbound is a leaf) SetSupernodeInward setSupernodeInward; vtkm::worklet::DispatcherMapField setSupernodeInwardDispatcher(setSupernodeInward); setSupernodeInwardDispatcher.Invoke(activeSupernodes, // input inbound, // input (whole array) outbound, // input (whole array) indegree, // input (whole array) outdegree, // input (whole array) superarcs); // i/o (whole array) outbound.ReleaseResources(); #ifdef DEBUG_PRINT DebugPrint(isJoin ? "Upper Regular Nodes Collapsed" : "Lower Regular Nodes Collapsed"); #endif } // CollapseRegular() // compresses trees to remove transferred vertices template void ContourTree::CompressTrees() { // Compute the number of log steps required in this pass vtkm::Id nActiveSupernodes = activeSupernodes.GetNumberOfValues(); vtkm::Id nLogSteps = 1; for (vtkm::Id shifter = nActiveSupernodes; shifter != 0; shifter >>= 1) nLogSteps++; // loop to update the merge trees for (vtkm::Id logStep = 0; logStep < nLogSteps; logStep++) { SkipVertex skipVertex; vtkm::worklet::DispatcherMapField skipVertexDispatcher(skipVertex); skipVertexDispatcher.Invoke(activeSupernodes, // input superarcs, // input (whole array) joinArcs, // i/o (whole array) splitArcs); // i/o (whole array) } #ifdef DEBUG_PRINT DebugPrint("Trees Compressed"); #endif } // CompressTrees() // compresses active set of supernodes template void ContourTree::CompressActiveSupernodes() { // copy only if the superarc is not set vtkm::cont::ArrayHandle noSuperarcArray; noSuperarcArray.Allocate(activeSupernodes.GetNumberOfValues()); VertexAssigned vertexAssigned(false); vtkm::worklet::DispatcherMapField vertexAssignedDispatcher(vertexAssigned); vertexAssignedDispatcher.Invoke(activeSupernodes, superarcs, noSuperarcArray); vtkm::cont::ArrayHandle compressSupernodes; DeviceAlgorithm::StreamCompact(activeSupernodes, noSuperarcArray, compressSupernodes); activeSupernodes.ReleaseResources(); DeviceAlgorithm::Copy(compressSupernodes, activeSupernodes); #ifdef DEBUG_PRINT DebugPrint("Active Supernodes Compressed"); #endif } // CompressActiveSupernodes() // recomputes the degree of each supernode from the join & split trees template void ContourTree::FindDegrees() { if (activeSupernodes.GetNumberOfValues() == 0) return; vtkm::Id nActiveSupernodes = activeSupernodes.GetNumberOfValues(); ResetDegrees resetDegrees; vtkm::worklet::DispatcherMapField resetDegreesDispatcher(resetDegrees); resetDegreesDispatcher.Invoke(activeSupernodes, // input updegree, // output (whole array) downdegree); // output (whole array) // create a temporary sorting array vtkm::cont::ArrayHandle sortVector; sortVector.Allocate(nActiveSupernodes); vtkm::cont::ArrayHandleIndex activeSupernodeIndexArray(nActiveSupernodes); // 1. Copy the neighbours for each active edge if (nActiveSupernodes > 0) { CopyNeighbors copyNeighbors; vtkm::worklet::DispatcherMapField copyNeighborsDispatcher(copyNeighbors); copyNeighborsDispatcher.Invoke(activeSupernodeIndexArray, // input activeSupernodes, // input (whole array) joinArcs, // input (whole array) sortVector); // output } // 2. Sort the neighbours DeviceAlgorithm::Sort(sortVector); // 3. For each value, store the beginning & end of the range (in parallel) // The 0th element is guaranteed to be NO_VERTEX_ASSIGNED, & can be skipped // Otherwise, if the i-1th element is different, we are the offset for the subrange // and store into the ith place vtkm::cont::ArrayHandleCounting subsetIndexArray(1, 1, nActiveSupernodes-1); if (nActiveSupernodes > 1) { DegreeSubrangeOffset degreeSubrangeOffset; vtkm::worklet::DispatcherMapField degreeSubrangeOffsetDispatcher(degreeSubrangeOffset); degreeSubrangeOffsetDispatcher.Invoke(subsetIndexArray, // input sortVector, // input (whole array) updegree); // output (whole array) } // 4. Compute the delta to get the degree. if (nActiveSupernodes > 1) { DegreeDelta degreeDelta(nActiveSupernodes); vtkm::worklet::DispatcherMapField degreeDeltaDispatcher(degreeDelta); degreeDeltaDispatcher.Invoke(subsetIndexArray, // input sortVector, // input updegree); // in out } // Now repeat the same steps for the downdegree // 1. Copy the neighbours for each active edge if (nActiveSupernodes > 0) { CopyNeighbors copyNeighbors; vtkm::worklet::DispatcherMapField copyNeighborsDispatcher(copyNeighbors); copyNeighborsDispatcher.Invoke(activeSupernodeIndexArray, // input activeSupernodes, // input (whole array) splitArcs, // input (whole array) sortVector); // output } // 2. Sort the neighbours DeviceAlgorithm::Sort(sortVector); // 3. For each value, store the beginning & end of the range (in parallel) // The 0th element is guaranteed to be NO_VERTEX_ASSIGNED, & can be skipped // Otherwise, if the i-1th element is different, we are the offset for the subrange // and store into the ith place if (nActiveSupernodes > 1) { DegreeSubrangeOffset degreeSubrangeOffset; vtkm::worklet::DispatcherMapField degreeSubrangeOffsetDispatcher(degreeSubrangeOffset); degreeSubrangeOffsetDispatcher.Invoke(subsetIndexArray, // input sortVector, // input (whole array) downdegree); // output (whole array) } // 4. Compute the delta to get the degree. if (nActiveSupernodes > 1) { DegreeDelta degreeDelta(nActiveSupernodes); vtkm::worklet::DispatcherMapField degreeDeltaDispatcher(degreeDelta); degreeDeltaDispatcher.Invoke(subsetIndexArray, // input sortVector, // input (whole array) downdegree); // in out (whole array) } #ifdef DEBUG_PRINT DebugPrint("Degrees Recomputed"); #endif } // FindDegrees() // small class for storing the contour arcs class EdgePair { public: vtkm::Id low, high; // constructor - defaults to -1 EdgePair(vtkm::Id Low = -1, vtkm::Id High = -1) : low(Low), high(High) {} }; // comparison operator < bool operator < (const EdgePair LHS, const EdgePair RHS) { if (LHS.low < RHS.low) return true; if (LHS.low > RHS.low) return false; if (LHS.high < RHS.high) return true; if (LHS.high > RHS.high) return false; return false; } struct SaddlePeakSort { VTKM_EXEC_CONT bool operator()(const vtkm::Pair& a, const vtkm::Pair& b) const { if (a.first < b.first) return true; if (a.first > b.first) return false; if (a.second < b.second) return true; if (a.second > b.second) return false; return false; } }; // sorted print routine template void ContourTree::CollectSaddlePeak( vtkm::cont::ArrayHandle > &saddlePeak) { // Collect the valid saddle peak pairs std::vector > superarcVector; for (vtkm::Id supernode = 0; supernode < supernodes.GetNumberOfValues(); supernode++) { // ID of regular node vtkm::Id regularID = supernodes.GetPortalConstControl().Get(supernode); // retrieve ID of target supernode vtkm::Id superTo = superarcs.GetPortalConstControl().Get(supernode); // if this is true, it is the last pruned vertex if (superTo == NO_VERTEX_ASSIGNED) continue; // retrieve the regular ID for it vtkm::Id regularTo = supernodes.GetPortalConstControl().Get(superTo); // how we print depends on which end has lower ID if (regularID < regularTo) { // from is lower // extra test to catch duplicate edge if (superarcs.GetPortalConstControl().Get(superTo) != supernode) superarcVector.push_back(vtkm::make_Pair(regularID, regularTo)); } // from is lower else superarcVector.push_back(vtkm::make_Pair(regularTo, regularID)); } // per vertex // Setting saddlePeak reference to the make_ArrayHandle directly does not work vtkm::Id nSuperarcs = superarcVector.size(); vtkm::cont::ArrayHandle > tempArray = vtkm::cont::make_ArrayHandle(superarcVector); // now sort it DeviceAlgorithm::Sort(tempArray, SaddlePeakSort()); DeviceAlgorithm::Copy(tempArray, saddlePeak); #ifdef DEBUG_PRINT for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++) { std::cout << setw(PRINT_WIDTH) << saddlePeak.GetPortalControl().Get(superarc).first << " "; std::cout << setw(PRINT_WIDTH) << saddlePeak.GetPortalControl().Get(superarc).second << std::endl; } #endif } // CollectSaddlePeak() // debug routine template void ContourTree::DebugPrint(const char *message) { std::cout << "---------------------------" << std::endl; std::cout << string(message) << std::endl; std::cout << "---------------------------" << std::endl; std::cout << std::endl; // print out the supernode arrays vtkm::Id nSupernodes = supernodes.GetNumberOfValues(); printHeader(nSupernodes); printIndices("Supernodes", supernodes); vtkm::cont::ArrayHandle supervalues; DeviceAlgorithm::Copy(PermuteValueType(supernodes, values), supervalues); printValues("Value", supervalues); printIndices("Up degree", updegree); printIndices("Down degree", downdegree); printIndices("Join arc", joinArcs); printIndices("Split arc", splitArcs); printIndices("Superarcs", superarcs); std::cout << std::endl; // print out the active supernodes vtkm::Id nActiveSupernodes = activeSupernodes.GetNumberOfValues(); printHeader(nActiveSupernodes); printIndices("Active Supernodes", activeSupernodes); vtkm::cont::ArrayHandle activeUpdegree; DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, updegree), activeUpdegree); printIndices("Active Up Degree", activeUpdegree); vtkm::cont::ArrayHandle activeDowndegree; DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, downdegree), activeDowndegree); printIndices("Active Down Degree", activeDowndegree); vtkm::cont::ArrayHandle activeJoinArcs; DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, joinArcs), activeJoinArcs); printIndices("Active Join Arcs", activeJoinArcs); vtkm::cont::ArrayHandle activeSplitArcs; DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, splitArcs), activeSplitArcs); printIndices("Active Split Arcs", activeSplitArcs); vtkm::cont::ArrayHandle activeSuperarcs; DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, superarcs), activeSuperarcs); printIndices("Active Superarcs", activeSuperarcs); std::cout << std::endl; } // DebugPrint() } } } #endif