5dd346007b
clang-format BinPack settings have been disabled to make sure that the VTK-m style guideline is obeyed.
985 lines
38 KiB
C++
985 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 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.
|
|
//============================================================================
|
|
// 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, typename DeviceAdapter>
|
|
class ContourTree
|
|
{
|
|
public:
|
|
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithm;
|
|
|
|
typedef vtkm::cont::ArrayHandle<vtkm::Id> IdArrayType;
|
|
typedef vtkm::cont::ArrayHandle<T> ValueArrayType;
|
|
typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, IdArrayType> PermuteIndexType;
|
|
typedef vtkm::cont::ArrayHandlePermutation<IdArrayType, ValueArrayType> PermuteValueType;
|
|
|
|
// 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, DeviceAdapter> &joinTree, &splitTree;
|
|
|
|
// references to join & split graphs
|
|
ChainGraph<T, StorageType, DeviceAdapter> &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, DeviceAdapter>& JoinTree,
|
|
MergeTree<T, StorageType, DeviceAdapter>& SplitTree,
|
|
ChainGraph<T, StorageType, DeviceAdapter>& JoinGraph,
|
|
ChainGraph<T, StorageType, DeviceAdapter>& 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:
|
|
typedef void ControlSignature(FieldIn<IdType> supernode,
|
|
WholeArrayIn<IdType> superarcs,
|
|
FieldOut<IdType> hasSuperArc);
|
|
typedef _3 ExecutionSignature(_1, _2);
|
|
typedef _1 InputDomain;
|
|
|
|
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, typename DeviceAdapter>
|
|
ContourTree<T, StorageType, DeviceAdapter>::ContourTree(
|
|
const vtkm::cont::ArrayHandle<T, StorageType>& Values,
|
|
MergeTree<T, StorageType, DeviceAdapter>& JoinTree,
|
|
MergeTree<T, StorageType, DeviceAdapter>& SplitTree,
|
|
ChainGraph<T, StorageType, DeviceAdapter>& JoinGraph,
|
|
ChainGraph<T, StorageType, DeviceAdapter>& 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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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);
|
|
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<vtkm::Id> regularToCritical;
|
|
vtkm::cont::ArrayHandleConstant<vtkm::Id> noVertArray(NO_VERTEX_ASSIGNED, nValues);
|
|
DeviceAlgorithm::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;
|
|
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<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
|
|
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<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
|
|
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<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;
|
|
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<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();
|
|
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 <typename T, typename StorageType, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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);
|
|
DeviceAlgorithm::Copy(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)
|
|
{
|
|
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<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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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;
|
|
DeviceAlgorithm::CopyIf(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 <typename T, typename StorageType, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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
|
|
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<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
|
|
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<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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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
|
|
DeviceAlgorithm::Sort(tempArray, SaddlePeakSort());
|
|
DeviceAlgorithm::Copy(tempArray, saddlePeak);
|
|
|
|
#ifdef DEBUG_PRINT
|
|
for (vtkm::UInt32 superarc = 0; superarc < superarcVector.size(); 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, typename DeviceAdapter>
|
|
void ContourTree<T, StorageType, DeviceAdapter>::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;
|
|
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<vtkm::Id> activeUpdegree;
|
|
DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, updegree), activeUpdegree);
|
|
printIndices("Active Up Degree", activeUpdegree);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeDowndegree;
|
|
DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, downdegree), activeDowndegree);
|
|
printIndices("Active Down Degree", activeDowndegree);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeJoinArcs;
|
|
DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, joinArcs), activeJoinArcs);
|
|
printIndices("Active Join Arcs", activeJoinArcs);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeSplitArcs;
|
|
DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, splitArcs), activeSplitArcs);
|
|
printIndices("Active Split Arcs", activeSplitArcs);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> activeSuperarcs;
|
|
DeviceAlgorithm::Copy(PermuteIndexType(activeSupernodes, superarcs), activeSuperarcs);
|
|
printIndices("Active Superarcs", activeSuperarcs);
|
|
std::cout << std::endl;
|
|
} // DebugPrint()
|
|
}
|
|
}
|
|
}
|
|
|
|
#endif
|