VTK-m  2.0
worklet/ContourTreeUniformAugmented.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 // Copyright (c) 2018, The Regents of the University of California, through
11 // Lawrence Berkeley National Laboratory (subject to receipt of any required approvals
12 // from the U.S. Dept. of Energy). All rights reserved.
13 //
14 // Redistribution and use in source and binary forms, with or without modification,
15 // are permitted provided that the following conditions are met:
16 //
17 // (1) Redistributions of source code must retain the above copyright notice, this
18 // list of conditions and the following disclaimer.
19 //
20 // (2) Redistributions in binary form must reproduce the above copyright notice,
21 // this list of conditions and the following disclaimer in the documentation
22 // and/or other materials provided with the distribution.
23 //
24 // (3) Neither the name of the University of California, Lawrence Berkeley National
25 // Laboratory, U.S. Dept. of Energy nor the names of its contributors may be
26 // used to endorse or promote products derived from this software without
27 // specific prior written permission.
28 //
29 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
30 // ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
31 // WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
32 // IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
33 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
34 // BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
35 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
36 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
37 // OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
38 // OF THE POSSIBILITY OF SUCH DAMAGE.
39 //
40 //=============================================================================
41 //
42 // This code is an extension of the algorithm presented in the paper:
43 // Parallel Peak Pruning for Scalable SMP Contour Tree Computation.
44 // Hamish Carr, Gunther Weber, Christopher Sewell, and James Ahrens.
45 // Proceedings of the IEEE Symposium on Large Data Analysis and Visualization
46 // (LDAV), October 2016, Baltimore, Maryland.
47 //
48 // The PPP2 algorithm and software were jointly developed by
49 // Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and
50 // Oliver Ruebel (LBNL)
51 //==============================================================================
52 
53 #ifndef vtk_m_worklet_ContourTreeUniformAugmented_h
54 #define vtk_m_worklet_ContourTreeUniformAugmented_h
55 
56 
57 #include <sstream>
58 #include <utility>
59 
60 // VTKM includes
61 #include <vtkm/Math.h>
62 #include <vtkm/Types.h>
63 #include <vtkm/cont/ArrayHandle.h>
65 #include <vtkm/cont/Field.h>
66 #include <vtkm/cont/Timer.h>
69 
70 // Contour tree worklet includes
82 
83 namespace vtkm
84 {
85 namespace worklet
86 {
87 
90 {
91 public:
98 
100  std::string TimingsLogString;
101 
102 
123  template <typename FieldType,
124  typename StorageType,
125  typename MeshType,
126  typename MeshBoundaryMeshExecType>
128  MeshType& mesh,
131  vtkm::Id& nIterations,
132  unsigned int computeRegularStructure,
133  const MeshBoundaryMeshExecType& meshBoundary)
134  {
136  fieldArray, // Just a place-holder to fill the required field. Used when calling SortData on the contour tree which is a no-op
137  contourTree,
138  sortOrder,
139  nIterations,
140  mesh,
141  computeRegularStructure,
142  meshBoundary);
143  return;
144  }
145 
168  template <typename FieldType, typename StorageType>
172  vtkm::Id& nIterations,
173  const vtkm::Id3 meshSize,
174  bool useMarchingCubes = false,
175  unsigned int computeRegularStructure = 1)
176  {
177  using namespace vtkm::worklet::contourtree_augmented;
178  // 2D Contour Tree
179  if (meshSize[2] == 1)
180  {
181  // Build the mesh and fill in the values
182  DataSetMeshTriangulation2DFreudenthal mesh(vtkm::Id2{ meshSize[0], meshSize[1] });
183  // Run the contour tree on the mesh
184  RunContourTree(fieldArray,
185  contourTree,
186  sortOrder,
187  nIterations,
188  mesh,
189  computeRegularStructure,
190  mesh.GetMeshBoundaryExecutionObject());
191  return;
192  }
193  // 3D Contour Tree using marching cubes
194  else if (useMarchingCubes)
195  {
196  // Build the mesh and fill in the values
198  // Run the contour tree on the mesh
199  RunContourTree(fieldArray,
200  contourTree,
201  sortOrder,
202  nIterations,
203  mesh,
204  computeRegularStructure,
206  return;
207  }
208  // 3D Contour Tree with Freudenthal
209  else
210  {
211  // Build the mesh and fill in the values
213  // Run the contour tree on the mesh
214  RunContourTree(fieldArray,
215  contourTree,
216  sortOrder,
217  nIterations,
218  mesh,
219  computeRegularStructure,
221  return;
222  }
223  }
224 
225 
226 private:
249  template <typename FieldType,
250  typename StorageType,
251  typename MeshClass,
252  typename MeshBoundaryClass>
256  vtkm::Id& nIterations,
257  MeshClass& mesh,
258  unsigned int computeRegularStructure,
259  const MeshBoundaryClass& meshBoundary)
260  {
261  using namespace vtkm::worklet::contourtree_augmented;
262  // Stage 1: Load the data into the mesh. This is done in the Run() method above and accessible
263  // here via the mesh parameter. The actual data load is performed outside of the
264  // worklet in the example contour tree app (or whoever uses the worklet)
265 
266  // Stage 2 : Sort the data on the mesh to initialize sortIndex & indexReverse on the mesh
267  // Start the timer for the mesh sort
268  vtkm::cont::Timer timer;
269  timer.Start();
270  std::stringstream timingsStream; // Use a string stream to log in one message
271 
272  // Sort the mesh data
273  mesh.SortData(fieldArray);
274  timingsStream << " " << std::setw(38) << std::left << "Sort Data"
275  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
276  timer.Start();
277 
278  // Stage 3: Assign every mesh vertex to a peak
279  MeshExtrema extrema(mesh.NumVertices);
280  extrema.SetStarts(mesh, true);
281  extrema.BuildRegularChains(true);
282  timingsStream << " " << std::setw(38) << std::left << "Join Tree Regular Chains"
283  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
284  timer.Start();
285 
286  // Stage 4: Identify join saddles & construct Active Join Graph
287  MergeTree joinTree(mesh.NumVertices, true);
288  ActiveGraph joinGraph(true);
289  joinGraph.Initialise(mesh, extrema);
290  timingsStream << " " << std::setw(38) << std::left << "Join Tree Initialize Active Graph"
291  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
292 
293 #ifdef DEBUG_PRINT
294  joinGraph.DebugPrint("Active Graph Instantiated", __FILE__, __LINE__);
295 #endif
296  timer.Start();
297 
298  // Stage 5: Compute Join Tree Hyperarcs from Active Join Graph
299  joinGraph.MakeMergeTree(joinTree, extrema);
300  timingsStream << " " << std::setw(38) << std::left << "Join Tree Compute"
301  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
302 #ifdef DEBUG_PRINT
303  joinTree.DebugPrint("Join tree Computed", __FILE__, __LINE__);
304  joinTree.DebugPrintTree("Join tree", __FILE__, __LINE__, mesh);
305 #endif
306  timer.Start();
307 
308  // Stage 6: Assign every mesh vertex to a pit
309  extrema.SetStarts(mesh, false);
310  extrema.BuildRegularChains(false);
311  timingsStream << " " << std::setw(38) << std::left << "Split Tree Regular Chains"
312  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
313  timer.Start();
314 
315  // Stage 7: Identify split saddles & construct Active Split Graph
316  MergeTree splitTree(mesh.NumVertices, false);
317  ActiveGraph splitGraph(false);
318  splitGraph.Initialise(mesh, extrema);
319  timingsStream << " " << std::setw(38) << std::left << "Split Tree Initialize Active Graph"
320  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
321 #ifdef DEBUG_PRINT
322  splitGraph.DebugPrint("Active Graph Instantiated", __FILE__, __LINE__);
323 #endif
324  timer.Start();
325 
326  // Stage 8: Compute Split Tree Hyperarcs from Active Split Graph
327  splitGraph.MakeMergeTree(splitTree, extrema);
328  timingsStream << " " << std::setw(38) << std::left << "Split Tree Compute"
329  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
330 #ifdef DEBUG_PRINT
331  splitTree.DebugPrint("Split tree Computed", __FILE__, __LINE__);
332  // Debug split and join tree
333  joinTree.DebugPrintTree("Join tree", __FILE__, __LINE__, mesh);
334  splitTree.DebugPrintTree("Split tree", __FILE__, __LINE__, mesh);
335 #endif
336  timer.Start();
337 
338  // Stage 9: Join & Split Tree are Augmented, then combined to construct Contour Tree
339  contourTree.Init(mesh.NumVertices);
340  ContourTreeMaker treeMaker(contourTree, joinTree, splitTree);
341  // 9.1 First we compute the hyper- and super- structure
342  treeMaker.ComputeHyperAndSuperStructure();
343  timingsStream << " " << std::setw(38) << std::left
344  << "Contour Tree Hyper and Super Structure"
345  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
346  timer.Start();
347 
348  // 9.2 Then we compute the regular structure
349  if (computeRegularStructure == 1) // augment with all vertices
350  {
351  treeMaker.ComputeRegularStructure(extrema);
352  timingsStream << " " << std::setw(38) << std::left << "Contour Tree Regular Structure"
353  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
354  }
355  else if (computeRegularStructure == 2) // augment by the mesh boundary
356  {
357  treeMaker.ComputeBoundaryRegularStructure(extrema, mesh, meshBoundary);
358  timingsStream << " " << std::setw(38) << std::left
359  << "Contour Tree Boundary Regular Structure"
360  << ": " << timer.GetElapsedTime() << " seconds" << std::endl;
361  }
362  timer.Start();
363 
364  // Collect the output data
365  nIterations = treeMaker.ContourTreeResult.NumIterations;
366  // Need to make a copy of sortOrder since ContourTreeMesh uses a smart array handle
367  // TODO: Check if we can just make sortOrder a return array with variable type or if we can make the SortOrder return optional
368  // TODO/FIXME: According to Ken Moreland the short answer is no. We may need to go back and refactor this when we
369  // improve the contour tree API. https://gitlab.kitware.com/vtk/vtk-m/-/merge_requests/2263#note_831128 for more details.
370  vtkm::cont::Algorithm::Copy(mesh.SortOrder, sortOrder);
371  // ProcessContourTree::CollectSortedSuperarcs<DeviceAdapter>(contourTree, mesh.SortOrder, saddlePeak);
372  // contourTree.SortedArcPrint(mesh.SortOrder);
373  // contourTree.PrintDotSuperStructure();
374 
375  // Log the collected timing results in one coherent log entry
376  this->TimingsLogString = timingsStream.str();
377  if (this->TimingsLogLevel != vtkm::cont::LogLevel::Off)
378  {
379  VTKM_LOG_S(this->TimingsLogLevel,
380  std::endl
381  << " ------------------- Contour Tree Worklet Timings ----------------------"
382  << std::endl
383  << this->TimingsLogString);
384  }
385  }
386 };
387 
388 } // namespace vtkm
389 } // namespace vtkm::worklet
390 
391 #endif // vtk_m_worklet_ContourTreeUniformAugmented_h
vtkm::cont::ArrayHandle
Manages an array-worth of data.
Definition: ArrayHandle.h:283
ArrayHandle.h
MeshBoundary2D.h
vtkm::worklet::ContourTreeAugmented::TimingsLogString
std::string TimingsLogString
Remember the results from our time-keeping so we can customize our logging.
Definition: worklet/ContourTreeUniformAugmented.h:100
vtkm::worklet::contourtree_augmented::MeshExtrema::BuildRegularChains
VTKM_CONT void BuildRegularChains(bool isMaximal)
Definition: MeshExtrema.h:129
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
Types.h
WorkletMapField.h
vtkm::worklet::contourtree_augmented::ContourTree
Definition: augmented/ContourTree.h:106
vtkm::cont::Timer::Start
VTKM_CONT void Start()
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation2DFreudenthal
Class representing a 2D dataset mesh with freudenthal triangulation connectivity for contour tree com...
Definition: DataSetMeshTriangulation2DFreudenthal.h:73
vtkm::cont::LogLevel
LogLevel
Log levels for use with the logging macros.
Definition: Logging.h:298
vtkm::worklet::contourtree_augmented::MeshExtrema::SetStarts
void SetStarts(MeshType &mesh, bool isMaximal)
Definition: MeshExtrema.h:151
vtkm::worklet::contourtree_augmented::MergeTree::DebugPrintTree
void DebugPrintTree(const char *message, const char *fileName, long lineNum, const ContourTreeMesh< FieldType > &mesh)
Definition: augmented/MergeTree.h:203
vtkm::worklet::contourtree_augmented
Definition: BuildChainsWorklet.h:63
MergeTree.h
vtkm::worklet::contourtree_augmented::MeshExtrema
Definition: MeshExtrema.h:79
vtkm::cont::Algorithm::Copy
static VTKM_CONT bool Copy(vtkm::cont::DeviceAdapterId devId, const vtkm::cont::ArrayHandle< T, CIn > &input, vtkm::cont::ArrayHandle< U, COut > &output)
Definition: Algorithm.h:410
vtkm::worklet::contourtree_augmented::ContourTree::Init
void Init(vtkm::Id dataSize)
Definition: augmented/ContourTree.h:196
vtkm::worklet::contourtree_augmented::ContourTreeMaker::ComputeHyperAndSuperStructure
void ComputeHyperAndSuperStructure()
Definition: ContourTreeMaker.h:193
vtkm::worklet::contourtree_augmented::ActiveGraph::DebugPrint
void DebugPrint(const char *message, const char *fileName, long lineNum)
Definition: ActiveGraph.h:917
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
DispatcherMapField.h
vtkm::worklet::contourtree_augmented::ContourTreeMaker::ComputeBoundaryRegularStructure
void ComputeBoundaryRegularStructure(MeshExtrema &meshExtrema, const Mesh &mesh, const MeshBoundaryExecObj &meshBoundary)
Definition: ContourTreeMaker.h:525
vtkm::worklet::contourtree_augmented::ContourTree::NumIterations
vtkm::Id NumIterations
Definition: augmented/ContourTree.h:153
ContourTreeMaker.h
ActiveGraph.h
vtkm::worklet::ContourTreeAugmented::TimingsLogLevel
vtkm::cont::LogLevel TimingsLogLevel
Definition: worklet/ContourTreeUniformAugmented.h:97
Math.h
ContourTreeMesh.h
MeshBoundary3D.h
Timer.h
vtkm::worklet::contourtree_augmented::MergeTree::DebugPrint
void DebugPrint(const char *message, const char *fileName, long lineNum)
Definition: augmented/MergeTree.h:168
vtkm::worklet::contourtree_augmented::ContourTreeMaker::ComputeRegularStructure
void ComputeRegularStructure(MeshExtrema &meshExtrema)
Definition: ContourTreeMaker.h:455
vtkm::worklet::contourtree_augmented::ContourTreeMaker
Definition: ContourTreeMaker.h:117
vtkm::worklet::contourtree_augmented::ContourTreeMaker::ContourTreeResult
ContourTree & ContourTreeResult
Definition: ContourTreeMaker.h:123
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation3DFreudenthal
Definition: DataSetMeshTriangulation3DFreudenthal.h:73
DataSetMesh.h
Types.h
MeshBoundaryContourTreeMesh.h
vtkm::cont::Timer
A class that can be used to time operations in VTK-m that might be occuring in parallel.
Definition: Timer.h:43
vtkm::worklet::ContourTreeAugmented
Compute the contour tree for 2d and 3d uniform grids and arbitrary topology graphs.
Definition: worklet/ContourTreeUniformAugmented.h:89
MeshExtrema.h
VTKM_LOG_S
#define VTKM_LOG_S(level,...)
Writes a message using stream syntax to the indicated log level.
Definition: Logging.h:261
vtkm::Vec< vtkm::Id, 3 >
vtkm::worklet::contourtree_augmented::MergeTree
Definition: augmented/MergeTree.h:77
vtkm::cont::LogLevel::Off
@ Off
Used with SetStderrLogLevel to silence the log.
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation3DMarchingCubes
Definition: DataSetMeshTriangulation3DMarchingCubes.h:72
vtkm::worklet::ContourTreeAugmented::Run
void Run(const vtkm::cont::ArrayHandle< FieldType, StorageType > fieldArray, contourtree_augmented::ContourTree &contourTree, contourtree_augmented::IdArrayType &sortOrder, vtkm::Id &nIterations, const vtkm::Id3 meshSize, bool useMarchingCubes=false, unsigned int computeRegularStructure=1)
Definition: worklet/ContourTreeUniformAugmented.h:169
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation3DMarchingCubes::GetMeshBoundaryExecutionObject
MeshBoundary3DExec GetMeshBoundaryExecutionObject() const
Definition: DataSetMeshTriangulation3DMarchingCubes.h:183
ArrayHandleCounting.h
Field.h
vtkm::worklet::contourtree_augmented::DataSetMeshTriangulation3DFreudenthal::GetMeshBoundaryExecutionObject
MeshBoundary3DExec GetMeshBoundaryExecutionObject() const
Definition: DataSetMeshTriangulation3DFreudenthal.h:151
vtkm::worklet::ContourTreeAugmented::Run
void Run(const vtkm::cont::ArrayHandle< FieldType, StorageType > fieldArray, MeshType &mesh, contourtree_augmented::ContourTree &contourTree, contourtree_augmented::IdArrayType &sortOrder, vtkm::Id &nIterations, unsigned int computeRegularStructure, const MeshBoundaryMeshExecType &meshBoundary)
Definition: worklet/ContourTreeUniformAugmented.h:127
vtkm::worklet::contourtree_augmented::ActiveGraph::Initialise
void Initialise(Mesh &mesh, const MeshExtrema &meshExtrema)
Definition: ActiveGraph.h:244
vtkm::worklet::contourtree_augmented::ActiveGraph
Definition: ActiveGraph.h:116
vtkm::worklet::contourtree_augmented::ActiveGraph::MakeMergeTree
void MakeMergeTree(MergeTree &tree, MeshExtrema &meshExtrema)
Definition: ActiveGraph.h:376
ContourTree.h
vtkm::cont::Timer::GetElapsedTime
VTKM_CONT vtkm::Float64 GetElapsedTime() const
Get the elapsed time measured by the given device adapter.
vtkm::worklet::ContourTreeAugmented::RunContourTree
void RunContourTree(const vtkm::cont::ArrayHandle< FieldType, StorageType > fieldArray, contourtree_augmented::ContourTree &contourTree, contourtree_augmented::IdArrayType &sortOrder, vtkm::Id &nIterations, MeshClass &mesh, unsigned int computeRegularStructure, const MeshBoundaryClass &meshBoundary)
Definition: worklet/ContourTreeUniformAugmented.h:253
vtkm::cont::LogLevel::Perf
@ Perf
General timing data and algorithm flow information, such as filter execution, worklet dispatches,...