mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-10-06 10:29:00 +00:00
bddad9b386
Now that the dispatcher does its own TryExecute, filters do not need to do that. This change requires all worklets called by filters to be able to execute without knowing the device a priori.
983 lines
38 KiB
C++
983 lines
38 KiB
C++
//============================================================================
|
|
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
// Copyright 2014 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-NA0003525 with NTESS,
|
|
// 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.
|
|
//============================================================================
|
|
// Copyright (c) 2016, Los Alamos National Security, LLC
|
|
// All rights reserved.
|
|
//
|
|
// Copyright 2016. Los Alamos National Security, LLC.
|
|
// This software was produced under U.S. Government contract DE-AC52-06NA25396
|
|
// for Los Alamos National Laboratory (LANL), which is operated by
|
|
// Los Alamos National Security, LLC for the U.S. Department of Energy.
|
|
// The U.S. Government has rights to use, reproduce, and distribute this
|
|
// software. NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC
|
|
// MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE
|
|
// USE OF THIS SOFTWARE. If software is modified to produce derivative works,
|
|
// such modified software should be clearly marked, so as not to confuse it
|
|
// with the version available from LANL.
|
|
//
|
|
// Additionally, 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 Los Alamos National Security, LLC, Los Alamos
|
|
// National Laboratory, LANL, the U.S. Government, 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 LOS ALAMOS NATIONAL SECURITY, LLC 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 LOS ALAMOS
|
|
// NATIONAL SECURITY, LLC 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 based on 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.
|
|
|
|
//=======================================================================================
|
|
//
|
|
// 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 <vtkm/worklet/contourtree/ChainGraph.h>
|
|
#include <vtkm/worklet/contourtree/CopyJoinSplit.h>
|
|
#include <vtkm/worklet/contourtree/CopyNeighbors.h>
|
|
#include <vtkm/worklet/contourtree/CopySupernodes.h>
|
|
#include <vtkm/worklet/contourtree/DegreeDelta.h>
|
|
#include <vtkm/worklet/contourtree/DegreeSubrangeOffset.h>
|
|
#include <vtkm/worklet/contourtree/FillSupernodes.h>
|
|
#include <vtkm/worklet/contourtree/FindLeaves.h>
|
|
#include <vtkm/worklet/contourtree/MergeTree.h>
|
|
#include <vtkm/worklet/contourtree/PrintVectors.h>
|
|
#include <vtkm/worklet/contourtree/RegularToCandidate.h>
|
|
#include <vtkm/worklet/contourtree/RegularToCriticalDown.h>
|
|
#include <vtkm/worklet/contourtree/RegularToCriticalUp.h>
|
|
#include <vtkm/worklet/contourtree/ResetDegrees.h>
|
|
#include <vtkm/worklet/contourtree/SetJoinAndSplitArcs.h>
|
|
#include <vtkm/worklet/contourtree/SetSupernodeInward.h>
|
|
#include <vtkm/worklet/contourtree/SkipVertex.h>
|
|
#include <vtkm/worklet/contourtree/SubrangeOffset.h>
|
|
#include <vtkm/worklet/contourtree/Types.h>
|
|
#include <vtkm/worklet/contourtree/UpdateOutbound.h>
|
|
|
|
#include <vtkm/Pair.h>
|
|
#include <vtkm/cont/ArrayHandle.h>
|
|
#include <vtkm/cont/ArrayHandleConcatenate.h>
|
|
#include <vtkm/cont/ArrayHandleConcatenate.h>
|
|
#include <vtkm/cont/ArrayHandlePermutation.h>
|
|
#include <vtkm/worklet/WorkletMapField.h>
|
|
|
|
//#define DEBUG_PRINT 1
|
|
//#define DEBUG_TIMING 1
|
|
|
|
namespace vtkm
|
|
{
|
|
namespace worklet
|
|
{
|
|
namespace contourtree
|
|
{
|
|
|
|
template <typename T, typename StorageType>
|
|
class ContourTree
|
|
{
|
|
public:
|
|
using IdArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
|
|
using ValueArrayType = vtkm::cont::ArrayHandle<T>;
|
|
using PermuteIndexType = vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType>;
|
|
using PermuteValueType = vtkm::cont::ArrayHandlePermutation<IdArrayType, ValueArrayType>;
|
|
|
|
// reference to the underlying data
|
|
const vtkm::cont::ArrayHandle<T, StorageType> values;
|
|
|
|
// vector of superarcs in the contour tree (stored as inward-pointing)
|
|
vtkm::cont::ArrayHandle<vtkm::Id> superarcs;
|
|
|
|
// vector of supernodes
|
|
vtkm::cont::ArrayHandle<vtkm::Id> supernodes;
|
|
|
|
// vector of supernodes still unprocessed
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeSupernodes;
|
|
|
|
// references to join & split trees
|
|
MergeTree<T, StorageType> &joinTree, &splitTree;
|
|
|
|
// references to join & split graphs
|
|
ChainGraph<T, StorageType> &joinGraph, &splitGraph;
|
|
|
|
// vectors of up & down degree used during computation
|
|
vtkm::cont::ArrayHandle<vtkm::Id> updegree, downdegree;
|
|
|
|
// vectors for tracking merge arcs
|
|
vtkm::cont::ArrayHandle<vtkm::Id> joinArcs, splitArcs;
|
|
|
|
// counter for how many iterations it took to compute
|
|
vtkm::Id nIterations;
|
|
|
|
// contour tree constructor
|
|
ContourTree(const vtkm::cont::ArrayHandle<T, StorageType>& Values,
|
|
MergeTree<T, StorageType>& JoinTree,
|
|
MergeTree<T, StorageType>& SplitTree,
|
|
ChainGraph<T, StorageType>& JoinGraph,
|
|
ChainGraph<T, StorageType>& 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<vtkm::Pair<vtkm::Id, vtkm::Id>>& saddlePeak);
|
|
|
|
void DebugPrint(const char* message);
|
|
|
|
}; // class ContourTree
|
|
|
|
struct VertexAssigned : vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
using ControlSignature = void(FieldIn<IdType> supernode,
|
|
WholeArrayIn<IdType> superarcs,
|
|
FieldOut<IdType> hasSuperArc);
|
|
using ExecutionSignature = _3(_1, _2);
|
|
using InputDomain = _1;
|
|
|
|
bool vertexIsAssigned;
|
|
|
|
VTKM_EXEC_CONT
|
|
VertexAssigned(bool VertexIsAssigned)
|
|
: vertexIsAssigned(VertexIsAssigned)
|
|
{
|
|
}
|
|
|
|
template <typename InPortalFieldType>
|
|
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 <typename T, typename StorageType>
|
|
ContourTree<T, StorageType>::ContourTree(const vtkm::cont::ArrayHandle<T, StorageType>& Values,
|
|
MergeTree<T, StorageType>& JoinTree,
|
|
MergeTree<T, StorageType>& SplitTree,
|
|
ChainGraph<T, StorageType>& JoinGraph,
|
|
ChainGraph<T, StorageType>& SplitGraph)
|
|
: values(Values)
|
|
, joinTree(JoinTree)
|
|
, splitTree(SplitTree)
|
|
, joinGraph(JoinGraph)
|
|
, splitGraph(SplitGraph)
|
|
{
|
|
|
|
// first we have to get the correct list of supernodes
|
|
// this will also set the degrees of the vertices initially
|
|
FindSupernodes();
|
|
|
|
// 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 <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::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<vtkm::Id> candidates;
|
|
|
|
// take the union of the two sets of vertices
|
|
vtkm::cont::ArrayHandleConcatenate<IdArrayType, IdArrayType> candidateArray(
|
|
joinGraph.valueIndex, splitGraph.valueIndex);
|
|
vtkm::cont::Algorithm::Copy(candidateArray, candidates);
|
|
vtkm::cont::Algorithm::Sort(candidates);
|
|
vtkm::cont::Algorithm::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<vtkm::Id> regularToCritical;
|
|
vtkm::cont::ArrayHandleConstant<vtkm::Id> noVertArray(NO_VERTEX_ASSIGNED, nValues);
|
|
vtkm::cont::Algorithm::Copy(noVertArray, regularToCritical);
|
|
|
|
if (nCandidates > 0)
|
|
{
|
|
RegularToCriticalUp regularToCriticalUp;
|
|
vtkm::worklet::DispatcherMapField<RegularToCriticalUp> 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<vtkm::Id> initCandidateArray(0, nCandidates);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> upCandidate;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> downCandidate;
|
|
vtkm::cont::Algorithm::Copy(initCandidateArray, upCandidate);
|
|
vtkm::cont::Algorithm::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<vtkm::Id> 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<RegularToCandidate> 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
|
|
vtkm::cont::Algorithm::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<vtkm::Id> subsetIndexArray(1, 1, nCandidates - 1);
|
|
if (nCandidates > 0)
|
|
{
|
|
SubrangeOffset subRangeOffset;
|
|
vtkm::worklet::DispatcherMapField<SubrangeOffset> 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<DegreeDelta> 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<RegularToCriticalDown> 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
|
|
vtkm::cont::Algorithm::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<SubrangeOffset> 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<DegreeDelta> 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<vtkm::Id> isSupernode;
|
|
isSupernode.Allocate(nCandidates);
|
|
|
|
// fill the vector in
|
|
if (nCandidates > 0)
|
|
{
|
|
FillSupernodes fillSupernodes;
|
|
vtkm::worklet::DispatcherMapField<FillSupernodes> 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<vtkm::Id> supernodeID;
|
|
vtkm::cont::Algorithm::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<CopySupernodes> 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<SetJoinAndSplitArcs> 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<vtkm::Id> newsuperarcs(NO_VERTEX_ASSIGNED, nSupernodes);
|
|
superarcs.ReleaseResources();
|
|
vtkm::cont::Algorithm::Copy(newsuperarcs, superarcs);
|
|
|
|
// create the active supernode vector
|
|
activeSupernodes.ReleaseResources();
|
|
activeSupernodes.Allocate(nSupernodes);
|
|
vtkm::cont::ArrayHandleIndex supernodeSeq(nSupernodes);
|
|
vtkm::cont::Algorithm::Copy(supernodeSeq, activeSupernodes);
|
|
|
|
#ifdef DEBUG_PRINT
|
|
DebugPrint("Supernodes Found");
|
|
#endif
|
|
} // FindSupernodes()
|
|
|
|
// transfers leaves from join/split trees to contour tree
|
|
template <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::TransferLeaves()
|
|
{
|
|
FindLeaves findLeaves;
|
|
vtkm::worklet::DispatcherMapField<FindLeaves> 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 <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::CollapseRegular(bool isJoin)
|
|
{
|
|
// we'll have a vector for tracking outwards
|
|
vtkm::Id nSupernodes = supernodes.GetNumberOfValues();
|
|
vtkm::cont::ArrayHandleConstant<vtkm::Id> nullArray(0, nSupernodes);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> outbound;
|
|
outbound.Allocate(nSupernodes);
|
|
vtkm::cont::ArrayCopy(nullArray, outbound);
|
|
|
|
// and a reference for the inwards array and to the degrees
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inbound;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> indegree;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> outdegree;
|
|
if (isJoin)
|
|
{
|
|
vtkm::cont::ArrayCopy(joinArcs, inbound);
|
|
vtkm::cont::ArrayCopy(downdegree, indegree);
|
|
vtkm::cont::ArrayCopy(updegree, outdegree);
|
|
}
|
|
else
|
|
{
|
|
vtkm::cont::ArrayCopy(splitArcs, inbound);
|
|
vtkm::cont::ArrayCopy(updegree, indegree);
|
|
vtkm::cont::ArrayCopy(downdegree, outdegree);
|
|
}
|
|
|
|
// loop to copy join/split
|
|
CopyJoinSplit copyJoinSplit;
|
|
vtkm::worklet::DispatcherMapField<CopyJoinSplit> 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<UpdateOutbound> 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<SetSupernodeInward> 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 <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::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<SkipVertex> 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 <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::CompressActiveSupernodes()
|
|
{
|
|
// copy only if the superarc is not set
|
|
vtkm::cont::ArrayHandle<vtkm::Id> noSuperarcArray;
|
|
noSuperarcArray.Allocate(activeSupernodes.GetNumberOfValues());
|
|
|
|
VertexAssigned vertexAssigned(false);
|
|
vtkm::worklet::DispatcherMapField<VertexAssigned> vertexAssignedDispatcher(vertexAssigned);
|
|
|
|
vertexAssignedDispatcher.Invoke(activeSupernodes, superarcs, noSuperarcArray);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> compressSupernodes;
|
|
vtkm::cont::Algorithm::CopyIf(activeSupernodes, noSuperarcArray, compressSupernodes);
|
|
|
|
activeSupernodes.ReleaseResources();
|
|
vtkm::cont::ArrayCopy(compressSupernodes, activeSupernodes);
|
|
|
|
#ifdef DEBUG_PRINT
|
|
DebugPrint("Active Supernodes Compressed");
|
|
#endif
|
|
} // CompressActiveSupernodes()
|
|
|
|
// recomputes the degree of each supernode from the join & split trees
|
|
template <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::FindDegrees()
|
|
{
|
|
if (activeSupernodes.GetNumberOfValues() == 0)
|
|
return;
|
|
|
|
vtkm::Id nActiveSupernodes = activeSupernodes.GetNumberOfValues();
|
|
ResetDegrees resetDegrees;
|
|
vtkm::worklet::DispatcherMapField<ResetDegrees> resetDegreesDispatcher(resetDegrees);
|
|
|
|
resetDegreesDispatcher.Invoke(activeSupernodes, // input
|
|
updegree, // output (whole array)
|
|
downdegree); // output (whole array)
|
|
|
|
// create a temporary sorting array
|
|
vtkm::cont::ArrayHandle<vtkm::Id> 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<CopyNeighbors> copyNeighborsDispatcher(copyNeighbors);
|
|
|
|
copyNeighborsDispatcher.Invoke(activeSupernodeIndexArray, // input
|
|
activeSupernodes, // input (whole array)
|
|
joinArcs, // input (whole array)
|
|
sortVector); // output
|
|
}
|
|
|
|
// 2. Sort the neighbours
|
|
vtkm::cont::Algorithm::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<vtkm::Id> subsetIndexArray(1, 1, nActiveSupernodes - 1);
|
|
if (nActiveSupernodes > 1)
|
|
{
|
|
DegreeSubrangeOffset degreeSubrangeOffset;
|
|
vtkm::worklet::DispatcherMapField<DegreeSubrangeOffset> 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<DegreeDelta> 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<CopyNeighbors> copyNeighborsDispatcher(copyNeighbors);
|
|
|
|
copyNeighborsDispatcher.Invoke(activeSupernodeIndexArray, // input
|
|
activeSupernodes, // input (whole array)
|
|
splitArcs, // input (whole array)
|
|
sortVector); // output
|
|
}
|
|
|
|
// 2. Sort the neighbours
|
|
vtkm::cont::Algorithm::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<DegreeSubrangeOffset> 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<DegreeDelta> 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<vtkm::Id, vtkm::Id>& a,
|
|
const vtkm::Pair<vtkm::Id, vtkm::Id>& 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 <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::CollectSaddlePeak(
|
|
vtkm::cont::ArrayHandle<vtkm::Pair<vtkm::Id, vtkm::Id>>& saddlePeak)
|
|
{
|
|
// Collect the valid saddle peak pairs
|
|
std::vector<vtkm::Pair<vtkm::Id, vtkm::Id>> 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::cont::ArrayHandle<vtkm::Pair<vtkm::Id, vtkm::Id>> tempArray =
|
|
vtkm::cont::make_ArrayHandle(superarcVector);
|
|
|
|
// now sort it
|
|
vtkm::cont::Algorithm::Sort(tempArray, SaddlePeakSort());
|
|
vtkm::cont::Algorithm::Copy(tempArray, saddlePeak);
|
|
|
|
#ifdef DEBUG_PRINT
|
|
const vtkm::Id arcVecSize = static_cast<vtkm::Id>(superarcVector.size());
|
|
for (vtkm::Id superarc = 0; superarc < arcVecSize; superarc++)
|
|
{
|
|
std::cout << std::setw(PRINT_WIDTH) << saddlePeak.GetPortalControl().Get(superarc).first << " ";
|
|
std::cout << std::setw(PRINT_WIDTH) << saddlePeak.GetPortalControl().Get(superarc).second
|
|
<< std::endl;
|
|
}
|
|
#endif
|
|
} // CollectSaddlePeak()
|
|
|
|
// debug routine
|
|
template <typename T, typename StorageType>
|
|
void ContourTree<T, StorageType>::DebugPrint(const char* message)
|
|
{
|
|
std::cout << "---------------------------" << std::endl;
|
|
std::cout << std::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<vtkm::Id> supervalues;
|
|
vtkm::cont::ArrayCopy(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<vtkm::Id> activeUpdegree;
|
|
vtkm::cont::ArrayCopy(PermuteIndexType(activeSupernodes, updegree), activeUpdegree);
|
|
printIndices("Active Up Degree", activeUpdegree);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeDowndegree;
|
|
vtkm::cont::ArrayCopy(PermuteIndexType(activeSupernodes, downdegree), activeDowndegree);
|
|
printIndices("Active Down Degree", activeDowndegree);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeJoinArcs;
|
|
vtkm::cont::ArrayCopy(PermuteIndexType(activeSupernodes, joinArcs), activeJoinArcs);
|
|
printIndices("Active Join Arcs", activeJoinArcs);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeSplitArcs;
|
|
vtkm::cont::ArrayCopy(PermuteIndexType(activeSupernodes, splitArcs), activeSplitArcs);
|
|
printIndices("Active Split Arcs", activeSplitArcs);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeSuperarcs;
|
|
vtkm::cont::ArrayCopy(PermuteIndexType(activeSupernodes, superarcs), activeSuperarcs);
|
|
printIndices("Active Superarcs", activeSuperarcs);
|
|
std::cout << std::endl;
|
|
} // DebugPrint()
|
|
}
|
|
}
|
|
}
|
|
|
|
#endif
|