VTK-m  2.0
ProcessContourTree.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 
54 #ifndef vtk_m_worklet_contourtree_augmented_process_contourtree_h
55 #define vtk_m_worklet_contourtree_augmented_process_contourtree_h
56 
57 // global includes
58 #include <algorithm>
59 #include <iomanip>
60 #include <iostream>
61 
62 // local includes
65 
66 //VTKM includes
67 #include <vtkm/Pair.h>
68 #include <vtkm/Types.h>
69 #include <vtkm/cont/Algorithm.h>
70 #include <vtkm/cont/ArrayCopy.h>
73 #include <vtkm/cont/Timer.h>
80 
81 #include <vtkm/cont/Invoker.h>
84 
85 //#define DEBUG_PRINT
86 
87 
90 
91 namespace vtkm
92 {
93 namespace worklet
94 {
95 namespace contourtree_augmented
96 {
97 
98 
99 // TODO Many of the post processing routines still need to be parallelized
100 // Class with routines for post processing the contour tree
102 { // class ProcessContourTree
103 public:
104  // initialises contour tree arrays - rest is done by another class
106  { // ProcessContourTree()
107  } // ProcessContourTree()
108 
109  // collect the sorted arcs
110  void static CollectSortedArcs(const ContourTree& contourTree,
111  const IdArrayType& sortOrder,
112  EdgePairArray& sortedArcs)
113  { // CollectSortedArcs
114  // create an array for sorting the arcs
115  std::vector<EdgePair> arcSorter;
116 
117  // fill it up
118  auto arcsPortal = contourTree.Arcs.ReadPortal();
119  auto sortOrderPortal = sortOrder.ReadPortal();
120 
121  for (vtkm::Id node = 0; node < contourTree.Arcs.GetNumberOfValues(); node++)
122  { // per node
123  // retrieve ID of target supernode
124  vtkm::Id arcTo = arcsPortal.Get(node);
125 
126  // if this is true, it is the last pruned vertex & is omitted
127  if (NoSuchElement(arcTo))
128  continue;
129 
130  // otherwise, strip out the flags
131  arcTo = MaskedIndex(arcTo);
132 
133  // now convert to mesh IDs from sort IDs
134  // otherwise, we need to convert the IDs to regular mesh IDs
135  vtkm::Id regularID = sortOrderPortal.Get(node);
136 
137  // retrieve the regular ID for it
138  vtkm::Id regularTo = sortOrderPortal.Get(arcTo);
139 
140  // how we print depends on which end has lower ID
141  if (regularID < regularTo)
142  arcSorter.push_back(EdgePair(regularID, regularTo));
143  else
144  arcSorter.push_back(EdgePair(regularTo, regularID));
145  } // per vertex
146 
147  // now sort it
148  // Setting saddlePeak reference to the make_ArrayHandle directly does not work
149  sortedArcs = vtkm::cont::make_ArrayHandle(arcSorter, vtkm::CopyFlag::On);
151  } // CollectSortedArcs
152 
153  // collect the sorted superarcs
154  void static CollectSortedSuperarcs(const ContourTree& contourTree,
155  const IdArrayType& sortOrder,
156  EdgePairArray& saddlePeak)
157  { // CollectSortedSuperarcs()
158  // create an array for sorting the arcs
159  std::vector<EdgePair> superarcSorter;
160 
161  // fill it up
162  auto supernodesPortal = contourTree.Supernodes.ReadPortal();
163  auto superarcsPortal = contourTree.Superarcs.ReadPortal();
164  auto sortOrderPortal = sortOrder.ReadPortal();
165 
166  for (vtkm::Id supernode = 0; supernode < contourTree.Supernodes.GetNumberOfValues();
167  supernode++)
168  { // per supernode
169  // sort ID of the supernode
170  vtkm::Id sortID = supernodesPortal.Get(supernode);
171 
172  // retrieve ID of target supernode
173  vtkm::Id superTo = superarcsPortal.Get(supernode);
174 
175  // if this is true, it is the last pruned vertex & is omitted
176  if (NoSuchElement(superTo))
177  continue;
178 
179  // otherwise, strip out the flags
180  superTo = MaskedIndex(superTo);
181 
182  // otherwise, we need to convert the IDs to regular mesh IDs
183  vtkm::Id regularID = sortOrderPortal.Get(MaskedIndex(sortID));
184 
185  // retrieve the regular ID for it
186  vtkm::Id regularTo = sortOrderPortal.Get(MaskedIndex(supernodesPortal.Get(superTo)));
187 
188  // how we print depends on which end has lower ID
189  if (regularID < regularTo)
190  { // from is lower
191  // extra test to catch duplicate edge
192  if (superarcsPortal.Get(superTo) != supernode)
193  {
194  superarcSorter.push_back(EdgePair(regularID, regularTo));
195  }
196  } // from is lower
197  else
198  {
199  superarcSorter.push_back(EdgePair(regularTo, regularID));
200  }
201  } // per vertex
202 
203  // Setting saddlePeak reference to the make_ArrayHandle directly does not work
204  saddlePeak = vtkm::cont::make_ArrayHandle(superarcSorter, vtkm::CopyFlag::On);
205 
206  // now sort it
208  } // CollectSortedSuperarcs()
209 
210  // routine to compute the volume for each hyperarc and superarc
211  void static ComputeVolumeWeightsSerial(const ContourTree& contourTree,
212  const vtkm::Id nIterations,
213  IdArrayType& superarcIntrinsicWeight,
214  IdArrayType& superarcDependentWeight,
215  IdArrayType& supernodeTransferWeight,
216  IdArrayType& hyperarcDependentWeight)
217  { // ContourTreeMaker::ComputeWeights()
218  // start by storing the first sorted vertex ID for each superarc
219  IdArrayType firstVertexForSuperparent;
220  firstVertexForSuperparent.Allocate(contourTree.Superarcs.GetNumberOfValues());
221  superarcIntrinsicWeight.Allocate(contourTree.Superarcs.GetNumberOfValues());
222  auto superarcIntrinsicWeightPortal = superarcIntrinsicWeight.WritePortal();
223  auto firstVertexForSuperparentPortal = firstVertexForSuperparent.WritePortal();
224  auto superparentsPortal = contourTree.Superparents.ReadPortal();
225  auto hyperparentsPortal = contourTree.Hyperparents.ReadPortal();
226  auto hypernodesPortal = contourTree.Hypernodes.ReadPortal();
227  auto hyperarcsPortal = contourTree.Hyperarcs.ReadPortal();
228  // auto superarcsPortal = contourTree.Superarcs.ReadPortal();
229  auto nodesPortal = contourTree.Nodes.ReadPortal();
230  // auto whenTransferredPortal = contourTree.WhenTransferred.ReadPortal();
231  for (vtkm::Id sortedNode = 0; sortedNode < contourTree.Arcs.GetNumberOfValues(); sortedNode++)
232  { // per node in sorted order
233  vtkm::Id sortID = nodesPortal.Get(sortedNode);
234  vtkm::Id superparent = superparentsPortal.Get(sortID);
235  if (sortedNode == 0)
236  firstVertexForSuperparentPortal.Set(superparent, sortedNode);
237  else if (superparent != superparentsPortal.Get(nodesPortal.Get(sortedNode - 1)))
238  firstVertexForSuperparentPortal.Set(superparent, sortedNode);
239  } // per node in sorted order
240  // now we use that to compute the intrinsic weights
241  for (vtkm::Id superarc = 0; superarc < contourTree.Superarcs.GetNumberOfValues(); superarc++)
242  if (superarc == contourTree.Superarcs.GetNumberOfValues() - 1)
243  superarcIntrinsicWeightPortal.Set(superarc,
244  contourTree.Arcs.GetNumberOfValues() -
245  firstVertexForSuperparentPortal.Get(superarc));
246  else
247  superarcIntrinsicWeightPortal.Set(superarc,
248  firstVertexForSuperparentPortal.Get(superarc + 1) -
249  firstVertexForSuperparentPortal.Get(superarc));
250 
251  // now initialise the arrays for transfer & dependent weights
254  superarcDependentWeight);
257  supernodeTransferWeight);
260  hyperarcDependentWeight);
261 
262  // set up the array which tracks which supernodes to deal with on which iteration
263  auto firstSupernodePerIterationPortal = contourTree.FirstSupernodePerIteration.ReadPortal();
264  auto firstHypernodePerIterationPortal = contourTree.FirstHypernodePerIteration.ReadPortal();
265  auto supernodeTransferWeightPortal = supernodeTransferWeight.WritePortal();
266  auto superarcDependentWeightPortal = superarcDependentWeight.WritePortal();
267  auto hyperarcDependentWeightPortal = hyperarcDependentWeight.WritePortal();
268 
269  /*
270  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, nIterations + 1),
271  firstSupernodePerIteration);
272  auto firstSupernodePerIterationPortal = firstSupernodePerIteration.WritePortal();
273  for (vtkm::Id supernode = 0; supernode < contourTree.Supernodes.GetNumberOfValues();
274  supernode++)
275  { // per supernode
276  vtkm::Id when = MaskedIndex(whenTransferredPortal.Get(supernode));
277  if (supernode == 0)
278  { // zeroth supernode
279  firstSupernodePerIterationPortal.Set(when, supernode);
280  } // zeroth supernode
281  else if (when != MaskedIndex(whenTransferredPortal.Get(supernode - 1)))
282  { // non-matching supernode
283  firstSupernodePerIterationPortal.Set(when, supernode);
284  } // non-matching supernode
285  } // per supernode
286  for (vtkm::Id iteration = 1; iteration < nIterations; ++iteration)
287  if (firstSupernodePerIterationPortal.Get(iteration) == 0)
288  firstSupernodePerIterationPortal.Set(iteration,
289  firstSupernodePerIterationPortal.Get(iteration + 1));
290 
291  // set the sentinel at the end of the array
292  firstSupernodePerIterationPortal.Set(nIterations, contourTree.Supernodes.GetNumberOfValues());
293 
294  // now use that array to construct a similar array for hypernodes
295  IdArrayType firstHypernodePerIteration;
296  firstHypernodePerIteration.Allocate(nIterations + 1);
297  auto firstHypernodePerIterationPortal = firstHypernodePerIteration.WritePortal();
298  auto supernodeTransferWeightPortal = supernodeTransferWeight.WritePortal();
299  auto superarcDependentWeightPortal = superarcDependentWeight.WritePortal();
300  auto hyperarcDependentWeightPortal = hyperarcDependentWeight.WritePortal();
301  for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
302  firstHypernodePerIterationPortal.Set(
303  iteration, hyperparentsPortal.Get(firstSupernodePerIterationPortal.Get(iteration)));
304  firstHypernodePerIterationPortal.Set(nIterations, contourTree.Hypernodes.GetNumberOfValues());
305  */
306 
307  // now iterate, propagating weights inwards
308  for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
309  { // per iteration
310  // pull the array bounds into register
311  vtkm::Id firstSupernode = firstSupernodePerIterationPortal.Get(iteration);
312  vtkm::Id lastSupernode = firstSupernodePerIterationPortal.Get(iteration + 1);
313  vtkm::Id firstHypernode = firstHypernodePerIterationPortal.Get(iteration);
314  vtkm::Id lastHypernode = firstHypernodePerIterationPortal.Get(iteration + 1);
315 
316  // Recall that the superarcs are sorted by (iteration, hyperarc), & that all superarcs for a given hyperarc are processed
317  // in the same iteration. Assume therefore that:
318  // i. we now have the intrinsic weight assigned for each superarc, and
319  // ii. we also have the transfer weight assigned for each supernode.
320  //
321  // Suppose we have a sequence of superarcs
322  // s11 s12 s13 s14 s21 s22 s23 s31
323  // with transfer weights at their origins and intrinsic weights along them
324  // sArc s11 s12 s13 s14 s21 s22 s23 s31
325  // transfer wt 0 1 2 1 2 3 1 0
326  // intrinsic wt 1 2 1 5 2 6 1 1
327  //
328  // now, if we do a prefix sum on each of these and add the two sums together, we get:
329  // sArc s11 s12 s13 s14 s21 s22 s23 s31
330  // hyperparent sNode ID s11 s11 s11 s11 s21 s21 s21 s31
331  // transfer weight 0 1 2 1 2 3 1 0
332  // intrinsic weight 1 2 1 5 2 6 1 1
333  // sum(xfer + intrinsic) 1 3 3 6 4 9 2 1
334  // prefix sum (xfer + int) 1 4 7 13 17 26 28 29
335  // prefix sum (xfer + int - previous hArc) 1 4 7 13 4 13 15 16
336 
337  // so, step 1: add xfer + int & store in dependent weight
338  for (vtkm::Id supernode = firstSupernode; supernode < lastSupernode; supernode++)
339  {
340  superarcDependentWeightPortal.Set(supernode,
341  supernodeTransferWeightPortal.Get(supernode) +
342  superarcIntrinsicWeightPortal.Get(supernode));
343  }
344 
345  // step 2: perform prefix sum on the dependent weight range
346  for (vtkm::Id supernode = firstSupernode + 1; supernode < lastSupernode; supernode++)
347  superarcDependentWeightPortal.Set(supernode,
348  superarcDependentWeightPortal.Get(supernode) +
349  superarcDependentWeightPortal.Get(supernode - 1));
350 
351  // step 3: subtract out the dependent weight of the prefix to the entire hyperarc. This will be a transfer, but for now, it's easier
352  // to show it in serial. NB: Loops backwards so that computation uses the correct value
353  // As a bonus, note that we test > firstsupernode, not >=. This is because we've got unsigned integers, & otherwise it will not terminate
354  // But the first is always correct anyway (same reason as the short-cut termination on hyperparent), so we're fine
355  for (vtkm::Id supernode = lastSupernode - 1; supernode > firstSupernode; supernode--)
356  { // per supernode
357  // retrieve the hyperparent & convert to a supernode ID
358  vtkm::Id hyperparent = hyperparentsPortal.Get(supernode);
359  vtkm::Id hyperparentSuperID = hypernodesPortal.Get(hyperparent);
360 
361  // if the hyperparent is the first in the sequence, dependent weight is already correct
362  if (hyperparent == firstHypernode)
363  continue;
364 
365  // otherwise, subtract out the dependent weight *immediately* before the hyperparent's supernode
366  superarcDependentWeightPortal.Set(
367  supernode,
368  superarcDependentWeightPortal.Get(supernode) -
369  superarcDependentWeightPortal.Get(hyperparentSuperID - 1));
370  } // per supernode
371 
372  // step 4: transfer the dependent weight to the hyperarc's target supernode
373  for (vtkm::Id hypernode = firstHypernode; hypernode < lastHypernode; hypernode++)
374  { // per hypernode
375  // last superarc for the hyperarc
376  vtkm::Id lastSuperarc;
377  // special case for the last hyperarc
378  if (hypernode == contourTree.Hypernodes.GetNumberOfValues() - 1)
379  // take the last superarc in the array
380  lastSuperarc = contourTree.Supernodes.GetNumberOfValues() - 1;
381  else
382  // otherwise, take the next hypernode's ID and subtract 1
383  lastSuperarc = hypernodesPortal.Get(hypernode + 1) - 1;
384 
385  // now, given the last superarc for the hyperarc, transfer the dependent weight
386  hyperarcDependentWeightPortal.Set(hypernode,
387  superarcDependentWeightPortal.Get(lastSuperarc));
388 
389  // note that in parallel, this will have to be split out as a sort & partial sum in another array
390  vtkm::Id hyperarcTarget = MaskedIndex(hyperarcsPortal.Get(hypernode));
391  supernodeTransferWeightPortal.Set(hyperarcTarget,
392  supernodeTransferWeightPortal.Get(hyperarcTarget) +
393  hyperarcDependentWeightPortal.Get(hypernode));
394  } // per hypernode
395  } // per iteration
396  } // ContourTreeMaker::ComputeWeights()
397 
398  // routine to compute the branch decomposition by volume
399  void static ComputeVolumeBranchDecompositionSerial(const ContourTree& contourTree,
400  const IdArrayType& superarcDependentWeight,
401  const IdArrayType& superarcIntrinsicWeight,
402  IdArrayType& whichBranch,
403  IdArrayType& branchMinimum,
404  IdArrayType& branchMaximum,
405  IdArrayType& branchSaddle,
406  IdArrayType& branchParent)
407  { // ComputeVolumeBranchDecomposition()
408  auto superarcDependentWeightPortal = superarcDependentWeight.ReadPortal();
409  auto superarcIntrinsicWeightPortal = superarcIntrinsicWeight.ReadPortal();
410 
411  // cache the number of non-root supernodes & superarcs
412  vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
413  vtkm::Id nSuperarcs = nSupernodes - 1;
414 
415  // STAGE I: Find the upward and downwards weight for each superarc, and set up arrays
416  IdArrayType upWeight;
417  upWeight.Allocate(nSuperarcs);
418  auto upWeightPortal = upWeight.WritePortal();
419  IdArrayType downWeight;
420  downWeight.Allocate(nSuperarcs);
421  auto downWeightPortal = downWeight.WritePortal();
422  IdArrayType bestUpward;
423  auto noSuchElementArray =
425  vtkm::cont::ArrayCopy(noSuchElementArray, bestUpward);
426  IdArrayType bestDownward;
427  vtkm::cont::ArrayCopy(noSuchElementArray, bestDownward);
428  vtkm::cont::ArrayCopy(noSuchElementArray, whichBranch);
429  auto bestUpwardPortal = bestUpward.WritePortal();
430  auto bestDownwardPortal = bestDownward.WritePortal();
431 
432  // STAGE II: Pick the best (largest volume) edge upwards and downwards
433  // II A. Pick the best upwards weight by sorting on lower vertex then processing by segments
434  // II A 1. Sort the superarcs by lower vertex
435  // II A 2. Per segment, best superarc writes to the best upwards array
438  superarcList);
439  auto superarcListPortal = superarcList.WritePortal();
440  vtkm::Id totalVolume = contourTree.Nodes.GetNumberOfValues();
441 #ifdef DEBUG_PRINT
442  std::cout << "Total Volume: " << totalVolume << std::endl;
443 #endif
444  auto superarcsPortal = contourTree.Superarcs.ReadPortal();
445 
446  // NB: Last element in array is guaranteed to be root superarc to infinity,
447  // so we can easily skip it by not indexing to the full size
448  for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++)
449  { // per superarc
450  if (IsAscending(superarcsPortal.Get(superarc)))
451  { // ascending superarc
452  superarcListPortal.Set(superarc,
453  EdgePair(superarc, MaskedIndex(superarcsPortal.Get(superarc))));
454  upWeightPortal.Set(superarc, superarcDependentWeightPortal.Get(superarc));
455  // at the inner end, dependent weight is the total in the subtree. Then there are vertices along the edge itself (intrinsic weight), including the supernode at the outer end
456  // So, to get the "dependent" weight in the other direction, we start with totalVolume - dependent, then subtract (intrinsic - 1)
457  downWeightPortal.Set(superarc,
458  (totalVolume - superarcDependentWeightPortal.Get(superarc)) +
459  (superarcIntrinsicWeightPortal.Get(superarc) - 1));
460  } // ascending superarc
461  else
462  { // descending superarc
463  superarcListPortal.Set(superarc,
464  EdgePair(MaskedIndex(superarcsPortal.Get(superarc)), superarc));
465  downWeightPortal.Set(superarc, superarcDependentWeightPortal.Get(superarc));
466  // at the inner end, dependent weight is the total in the subtree. Then there are vertices along the edge itself (intrinsic weight), including the supernode at the outer end
467  // So, to get the "dependent" weight in the other direction, we start with totalVolume - dependent, then subtract (intrinsic - 1)
468  upWeightPortal.Set(superarc,
469  (totalVolume - superarcDependentWeightPortal.Get(superarc)) +
470  (superarcIntrinsicWeightPortal.Get(superarc) - 1));
471  } // descending superarc
472  } // per superarc
473 
474 #ifdef DEBUG_PRINT
475  std::cout << "II A. Weights Computed" << std::endl;
476  PrintHeader(upWeight.GetNumberOfValues());
477  //PrintIndices("Intrinsic Weight", superarcIntrinsicWeight);
478  //PrintIndices("Dependent Weight", superarcDependentWeight);
479  PrintIndices("Upwards Weight", upWeight);
480  PrintIndices("Downwards Weight", downWeight);
481  std::cout << std::endl;
482 #endif
483 
484  // II B. Pick the best downwards weight by sorting on upper vertex then processing by segments
485  // II B 1. Sort the superarcs by upper vertex
486  IdArrayType superarcSorter;
487  superarcSorter.Allocate(nSuperarcs);
488  auto superarcSorterPortal = superarcSorter.WritePortal();
489  for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++)
490  superarcSorterPortal.Set(superarc, superarc);
491 
493  superarcSorter,
494  process_contourtree_inc_ns::SuperArcVolumetricComparator(upWeight, superarcList, false));
495 
496  // II B 2. Per segment, best superarc writes to the best upward array
497  for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++)
498  { // per superarc
499  vtkm::Id superarcID = superarcSorterPortal.Get(superarc);
500  const EdgePair& edge = superarcListPortal.Get(superarcID);
501  // if it's the last one
502  if (superarc == nSuperarcs - 1)
503  bestDownwardPortal.Set(edge.second, edge.first);
504  else
505  { // not the last one
506  const EdgePair& nextEdge = superarcListPortal.Get(superarcSorterPortal.Get(superarc + 1));
507  // if the next edge belongs to another, we're the highest
508  if (nextEdge.second != edge.second)
509  bestDownwardPortal.Set(edge.second, edge.first);
510  } // not the last one
511  } // per superarc
512 
513  // II B 3. Repeat for lower vertex
515  superarcSorter,
516  process_contourtree_inc_ns::SuperArcVolumetricComparator(downWeight, superarcList, true));
517 
518  // II B 2. Per segment, best superarc writes to the best upward array
519  for (vtkm::Id superarc = 0; superarc < nSuperarcs; superarc++)
520  { // per superarc
521  vtkm::Id superarcID = superarcSorterPortal.Get(superarc);
522  const EdgePair& edge = superarcListPortal.Get(superarcID);
523  // if it's the last one
524  if (superarc == nSuperarcs - 1)
525  bestUpwardPortal.Set(edge.first, edge.second);
526  else
527  { // not the last one
528  const EdgePair& nextEdge = superarcListPortal.Get(superarcSorterPortal.Get(superarc + 1));
529  // if the next edge belongs to another, we're the highest
530  if (nextEdge.first != edge.first)
531  bestUpwardPortal.Set(edge.first, edge.second);
532  } // not the last one
533  } // per superarc
534 
535 #ifdef DEBUG_PRINT
536  std::cout << "II. Best Edges Selected" << std::endl;
537  PrintHeader(bestUpward.GetNumberOfValues());
538  PrintIndices("Best Upwards", bestUpward);
539  PrintIndices("Best Downwards", bestDownward);
540  std::cout << std::endl;
541 #endif
542 
544  whichBranch,
545  branchMinimum,
546  branchMaximum,
547  branchSaddle,
548  branchParent,
549  bestUpward,
550  bestDownward);
551  }
552 
553  // routine to compute the branch decomposition by volume
554  void static ComputeBranchData(const ContourTree& contourTree,
555  IdArrayType& whichBranch,
556  IdArrayType& branchMinimum,
557  IdArrayType& branchMaximum,
558  IdArrayType& branchSaddle,
559  IdArrayType& branchParent,
560  IdArrayType& bestUpward,
561  IdArrayType& bestDownward)
562  { // ComputeBranchData()
563 
564  // Set up constants
565  vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
566  auto noSuchElementArray =
568  vtkm::cont::ArrayCopy(noSuchElementArray, whichBranch);
569 
570  // STAGE III: For each vertex, identify which neighbours are on same branch
571  // Let v = BestUp(u). Then if u = BestDown(v), copy BestUp(u) to whichBranch(u)
572  // Otherwise, let whichBranch(u) = BestUp(u) | TERMINAL to mark the end of the side branch
573  // NB 1: Leaves already have the flag set, but it's redundant so its safe
574  // NB 2: We don't need to do it downwards because it's symmetric
575  vtkm::cont::Invoker invoke;
577  propagateBestUpDownWorklet;
578  invoke(propagateBestUpDownWorklet, bestUpward, bestDownward, whichBranch);
579 
580 #ifdef DEBUG_PRINT
581  std::cout << "III. Branch Neighbours Identified" << std::endl;
582  PrintHeader(whichBranch.GetNumberOfValues());
583  PrintIndices("Which Branch", whichBranch);
584  std::cout << std::endl;
585 #endif
586 
587  // STAGE IV: Use pointer-doubling on whichBranch to propagate branches
588  // Compute the number of log steps required in this pass
589  vtkm::Id numLogSteps = 1;
590  for (vtkm::Id shifter = nSupernodes; shifter != 0; shifter >>= 1)
591  numLogSteps++;
592 
594  nSupernodes);
595 
596  // use pointer-doubling to build the branches
597  for (vtkm::Id iteration = 0; iteration < numLogSteps; iteration++)
598  { // per iteration
599  invoke(pointerDoubling, whichBranch);
600  } // per iteration
601 
602 
603 #ifdef DEBUG_PRINT
604  std::cout << "IV. Branch Chains Propagated" << std::endl;
605  PrintHeader(whichBranch.GetNumberOfValues());
606  PrintIndices("Which Branch", whichBranch);
607  std::cout << std::endl;
608 #endif
609 
610  // Initialise
611  IdArrayType chainToBranch;
613 
614  // Set 1 to every relevant
616  prepareChainToBranchWorklet;
617  invoke(prepareChainToBranchWorklet, whichBranch, chainToBranch);
618 
619  // Prefix scanto get IDs
620  vtkm::Id nBranches = vtkm::cont::Algorithm::ScanInclusive(chainToBranch, chainToBranch);
621 
623  finaliseChainToBranchWorklet;
624  invoke(finaliseChainToBranchWorklet, whichBranch, chainToBranch);
625 
626  // V B. Create the arrays for the branches
627  auto noSuchElementArrayNBranches =
629  vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchMinimum);
630  vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchMaximum);
631  vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchSaddle);
632  vtkm::cont::ArrayCopy(noSuchElementArrayNBranches, branchParent);
633 
634 #ifdef DEBUG_PRINT
635  std::cout << "V. Branch Arrays Created" << std::endl;
636  PrintHeader(chainToBranch.GetNumberOfValues());
637  PrintIndices("Chain To Branch", chainToBranch);
638  PrintHeader(nBranches);
639  PrintIndices("Branch Minimum", branchMinimum);
640  PrintIndices("Branch Maximum", branchMaximum);
641  PrintIndices("Branch Saddle", branchSaddle);
642  PrintIndices("Branch Parent", branchParent);
643 #endif
644 
645  IdArrayType supernodeSorter;
646  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleIndex(nSupernodes), supernodeSorter);
647 
649  supernodeSorter,
651 
652  IdArrayType permutedBranches;
653  permutedBranches.Allocate(nSupernodes);
654  PermuteArray<vtkm::Id>(whichBranch, supernodeSorter, permutedBranches);
655 
656  IdArrayType permutedRegularID;
657  permutedRegularID.Allocate(nSupernodes);
658  PermuteArray<vtkm::Id>(contourTree.Supernodes, supernodeSorter, permutedRegularID);
659 
660 #ifdef DEBUG_PRINT
661  std::cout << "VI A. Sorted into Branches" << std::endl;
662  PrintHeader(nSupernodes);
663  PrintIndices("Supernode IDs", supernodeSorter);
664  PrintIndices("Branch", permutedBranches);
665  PrintIndices("Regular ID", permutedRegularID);
666 #endif
667 
669  whichBranchNewIdWorklet;
670  invoke(whichBranchNewIdWorklet, chainToBranch, whichBranch);
671 
673  branchMinMaxSetWorklet(nSupernodes);
674  invoke(branchMinMaxSetWorklet, supernodeSorter, whichBranch, branchMinimum, branchMaximum);
675 
676 #ifdef DEBUG_PRINT
677  std::cout << "VI. Branches Set" << std::endl;
678  PrintHeader(nBranches);
679  PrintIndices("Branch Maximum", branchMaximum);
680  PrintIndices("Branch Minimum", branchMinimum);
681  PrintIndices("Branch Saddle", branchSaddle);
682  PrintIndices("Branch Parent", branchParent);
683 #endif
684 
686  branchSaddleParentSetWorklet;
687  invoke(branchSaddleParentSetWorklet,
688  whichBranch,
689  branchMinimum,
690  branchMaximum,
691  bestDownward,
692  bestUpward,
693  branchSaddle,
694  branchParent);
695 
696 #ifdef DEBUG_PRINT
697  std::cout << "VII. Branches Constructed" << std::endl;
698  PrintHeader(nBranches);
699  PrintIndices("Branch Maximum", branchMaximum);
700  PrintIndices("Branch Minimum", branchMinimum);
701  PrintIndices("Branch Saddle", branchSaddle);
702  PrintIndices("Branch Parent", branchParent);
703 #endif
704 
705  } // ComputeBranchData()
706 
707  // Create branch decomposition from contour tree
708  template <typename T, typename StorageType>
710  const IdArrayType& contourTreeSuperparents,
711  const IdArrayType& contourTreeSupernodes,
712  const IdArrayType& whichBranch,
713  const IdArrayType& branchMinimum,
714  const IdArrayType& branchMaximum,
715  const IdArrayType& branchSaddle,
716  const IdArrayType& branchParent,
717  const IdArrayType& sortOrder,
719  bool dataFieldIsSorted)
720  {
722  contourTreeSuperparents,
723  contourTreeSupernodes,
724  whichBranch,
725  branchMinimum,
726  branchMaximum,
727  branchSaddle,
728  branchParent,
729  sortOrder,
730  dataField,
731  dataFieldIsSorted);
732  }
733 
734  // routine to compute the branch decomposition by volume
735  void static ComputeVolumeBranchDecomposition(const ContourTree& contourTree,
736  const vtkm::Id nIterations,
737  IdArrayType& whichBranch,
738  IdArrayType& branchMinimum,
739  IdArrayType& branchMaximum,
740  IdArrayType& branchSaddle,
741  IdArrayType& branchParent)
742  { // ComputeHeightBranchDecomposition()
743 
744  vtkm::cont::Invoker Invoke;
745 
746  // STEP 1. Compute the number of nodes in every superarc, that's the intrinsic weight
747  IdArrayType superarcIntrinsicWeight;
748  superarcIntrinsicWeight.Allocate(contourTree.Superarcs.GetNumberOfValues());
749 
750  IdArrayType firstVertexForSuperparent;
751  firstVertexForSuperparent.Allocate(contourTree.Superarcs.GetNumberOfValues());
752 
753  // Compute the number of regular nodes on every superarcs (the intrinsic weight)
755  setFirstVertexForSuperparent;
756  Invoke(setFirstVertexForSuperparent,
757  contourTree.Nodes,
758  contourTree.Superparents,
759  firstVertexForSuperparent);
760 
762  computeIntrinsicWeight;
763  Invoke(computeIntrinsicWeight,
764  contourTree.Arcs,
765  contourTree.Superarcs,
766  firstVertexForSuperparent,
767  superarcIntrinsicWeight);
768 
769 
770  // Cache the number of non-root supernodes & superarcs
771  vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
772  auto noSuchElementArray =
774 
775  // Set up bestUpward and bestDownward array, these are the things we want to compute in this routine.
776  IdArrayType bestUpward, bestDownward;
777  vtkm::cont::ArrayCopy(noSuchElementArray, bestUpward);
778  vtkm::cont::ArrayCopy(noSuchElementArray, bestDownward);
779 
780  // We initiale with the weight of the superarcs, once we sum those up we'll get the hypersweep weight
781  IdArrayType sumValues;
782  vtkm::cont::ArrayCopy(superarcIntrinsicWeight, sumValues);
783 
784  // This should be 0 here, because we're not changing the root
788  howManyUsed);
789 
790  // Perform a sum hypersweep
791  hyperarcScan<decltype(vtkm::Sum())>(contourTree.Supernodes,
792  contourTree.Hypernodes,
793  contourTree.Hyperarcs,
794  contourTree.Hyperparents,
795  contourTree.Hyperparents,
796  contourTree.WhenTransferred,
797  howManyUsed,
798  nIterations,
799  vtkm::Sum(),
800  sumValues);
801 
802  // For every directed arc store the volume of it's associate subtree
804  arcs.Allocate(contourTree.Superarcs.GetNumberOfValues() * 2 - 2);
805 
806  vtkm::Id totalVolume = contourTree.Nodes.GetNumberOfValues();
808  totalVolume);
809  Invoke(initArcs, sumValues, superarcIntrinsicWeight, contourTree.Superarcs, arcs);
810 
811  // Sort arcs to obtain the best up and down
813 
815  Invoke(setBestUpDown, bestUpward, bestDownward, arcs);
816 
818  whichBranch,
819  branchMinimum,
820  branchMaximum,
821  branchSaddle,
822  branchParent,
823  bestUpward,
824  bestDownward);
825 
826  } // ComputeHeightBranchDecomposition()
827 
828  // routine to compute the branch decomposition by volume
829  void static ComputeHeightBranchDecomposition(const ContourTree& contourTree,
830  const cont::ArrayHandle<Float64> fieldValues,
831  const IdArrayType& ctSortOrder,
832  const vtkm::Id nIterations,
833  IdArrayType& whichBranch,
834  IdArrayType& branchMinimum,
835  IdArrayType& branchMaximum,
836  IdArrayType& branchSaddle,
837  IdArrayType& branchParent)
838  { // ComputeHeightBranchDecomposition()
839 
840  // Cache the number of non-root supernodes & superarcs
841  vtkm::Id nSupernodes = contourTree.Supernodes.GetNumberOfValues();
842  auto noSuchElementArray =
844 
845  // Set up bestUpward and bestDownward array, these are the things we want to compute in this routine.
846  IdArrayType bestUpward, bestDownward;
847  vtkm::cont::ArrayCopy(noSuchElementArray, bestUpward);
848  vtkm::cont::ArrayCopy(noSuchElementArray, bestDownward);
849 
850  // maxValues and minValues store the values from the max and min hypersweep respectively.
851  IdArrayType minValues, maxValues;
852  vtkm::cont::ArrayCopy(contourTree.Supernodes, maxValues);
853  vtkm::cont::ArrayCopy(contourTree.Supernodes, minValues);
854 
855  // Store the direction of the superarcs in the min and max hypersweep (find a way to get rid of these, the only differing direction is on the path from the root to the min/max).
856  IdArrayType minParents, maxParents;
857  vtkm::cont::ArrayCopy(contourTree.Superarcs, minParents);
858  vtkm::cont::ArrayCopy(contourTree.Superarcs, maxParents);
859 
860  auto minParentsPortal = minParents.WritePortal();
861  auto maxParentsPortal = maxParents.WritePortal();
862 
863  // Cache the glonal minimum and global maximum (these will be the roots in the min and max hypersweep)
864  Id minSuperNode = MaskedIndex(contourTree.Superparents.ReadPortal().Get(0));
865  Id maxSuperNode = MaskedIndex(
866  contourTree.Superparents.ReadPortal().Get(contourTree.Nodes.GetNumberOfValues() - 1));
867 
868  // Find the path from the global minimum to the root, not parallelisable (but it's fast, no need to parallelise)
869  auto minPath = findSuperPathToRoot(contourTree.Superarcs.ReadPortal(), minSuperNode);
870 
871  // Find the path from the global minimum to the root, not parallelisable (but it's fast, no need to parallelise)
872  auto maxPath = findSuperPathToRoot(contourTree.Superarcs.ReadPortal(), maxSuperNode);
873 
874  // Reserve the direction of the superarcs on the min path.
875  for (std::size_t i = 1; i < minPath.size(); i++)
876  {
877  minParentsPortal.Set(minPath[i], minPath[i - 1]);
878  }
879  minParentsPortal.Set(minPath[0], 0);
880 
881  // Reserve the direction of the superarcs on the max path.
882  for (std::size_t i = 1; i < maxPath.size(); i++)
883  {
884  maxParentsPortal.Set(maxPath[i], maxPath[i - 1]);
885  }
886  maxParentsPortal.Set(maxPath[0], 0);
887 
888  vtkm::cont::Invoker Invoke;
890  Invoke(unmaskArrayWorklet, minValues);
891  Invoke(unmaskArrayWorklet, maxValues);
892 
893  // Thse arrays hold the changes hyperarcs in the min and max hypersweep respectively
894  vtkm::cont::ArrayHandle<vtkm::Id> minHyperarcs, maxHyperarcs;
895  vtkm::cont::ArrayCopy(contourTree.Hyperarcs, minHyperarcs);
896  vtkm::cont::ArrayCopy(contourTree.Hyperarcs, maxHyperarcs);
897 
898  // These arrays hold the changed hyperarcs for the min and max hypersweep
899  vtkm::cont::ArrayHandle<vtkm::Id> minHyperparents, maxHyperparents;
900  vtkm::cont::ArrayCopy(contourTree.Hyperparents, minHyperparents);
901  vtkm::cont::ArrayCopy(contourTree.Hyperparents, maxHyperparents);
902 
903  auto minHyperparentsPortal = minHyperparents.WritePortal();
904  auto maxHyperparentsPortal = maxHyperparents.WritePortal();
905 
906  for (std::size_t i = 0; i < minPath.size(); i++)
907  {
908  // Set a unique dummy Id (something that the prefix scan by key will leave alone)
909  minHyperparentsPortal.Set(minPath[i],
910  contourTree.Hypernodes.GetNumberOfValues() + minPath[i]);
911  }
912 
913  for (std::size_t i = 0; i < maxPath.size(); i++)
914  {
915  // Set a unique dummy Id (something that the prefix scan by key will leave alone)
916  maxHyperparentsPortal.Set(maxPath[i],
917  contourTree.Hypernodes.GetNumberOfValues() + maxPath[i]);
918  }
919 
920  // These arrays hold the number of nodes in each hypearcs that are on the min or max path for the min and max hypersweep respectively.
921  vtkm::cont::ArrayHandle<vtkm::Id> minHowManyUsed, maxHowManyUsed;
924  minHowManyUsed);
927  maxHowManyUsed);
928 
929  // Min Hypersweep
930  const auto minOperator = vtkm::Minimum();
931 
932  // Cut hyperarcs at the first node on the path from the max to the root
933  editHyperarcs(contourTree.Hyperparents.ReadPortal(),
934  minPath,
935  minHyperarcs.WritePortal(),
936  minHowManyUsed.WritePortal());
937 
938  // Perform an ordinary hypersweep on those new hyperarcs
939  hyperarcScan<decltype(vtkm::Minimum())>(contourTree.Supernodes,
940  contourTree.Hypernodes,
941  minHyperarcs,
942  contourTree.Hyperparents,
943  minHyperparents,
944  contourTree.WhenTransferred,
945  minHowManyUsed,
946  nIterations,
947  vtkm::Minimum(),
948  minValues);
949 
950  // Prefix sum along the path from the min to the root
951  fixPath(vtkm::Minimum(), minPath, minValues.WritePortal());
952 
953  // Max Hypersweep
954  const auto maxOperator = vtkm::Maximum();
955 
956  // Cut hyperarcs at the first node on the path from the max to the root
957  editHyperarcs(contourTree.Hyperparents.ReadPortal(),
958  maxPath,
959  maxHyperarcs.WritePortal(),
960  maxHowManyUsed.WritePortal());
961 
962  // Perform an ordinary hypersweep on those new hyperarcs
963  hyperarcScan<decltype(vtkm::Maximum())>(contourTree.Supernodes,
964  contourTree.Hypernodes,
965  maxHyperarcs,
966  contourTree.Hyperparents,
967  maxHyperparents,
968  contourTree.WhenTransferred,
969  maxHowManyUsed,
970  nIterations,
971  vtkm::Maximum(),
972  maxValues);
973 
974  // Prefix sum along the path from the max to the root
975  fixPath(vtkm::Maximum(), maxPath, maxValues.WritePortal());
976 
977  // For every directed edge (a, b) consider that subtree who's root is b and does not contain a.
978  // We have so far found the min and max in all sub subtrees, now we compare those to a and incorporate a into that.
980  vtkm::Minimum())>
981  incorporateParentMinimumWorklet(minOperator);
982  Invoke(incorporateParentMinimumWorklet, minParents, contourTree.Supernodes, minValues);
983 
985  vtkm::Maximum())>
986  incorporateParentMaximumWorklet(maxOperator);
987  Invoke(incorporateParentMaximumWorklet, maxParents, contourTree.Supernodes, maxValues);
988 
989  // Initialise all directed superarcs in the contour tree. Those will correspond to subtrees whos height we need for the branch decomposition.
991  arcs.Allocate(contourTree.Superarcs.GetNumberOfValues() * 2 - 2);
992 
994  0, contourTree.Arcs.GetNumberOfValues() - 1, minPath[minPath.size() - 1]);
995 
996  Invoke(initArcs, minParents, maxParents, minValues, maxValues, contourTree.Superarcs, arcs);
997 
998  // Use the min & max to compute the height of all subtrees
1000  computeSubtreeHeight;
1001  Invoke(computeSubtreeHeight, fieldValues, ctSortOrder, contourTree.Supernodes, arcs);
1002 
1003  // Sort all directed edges based on the height of their subtree
1005 
1006  // Select a best up and best down neighbour for every vertex in the contour tree using heights of all subtrees
1008  Invoke(setBestUpDown, bestUpward, bestDownward, arcs);
1009 
1010  // Having computed the bestUp/Down we can propagte those to obtain the branches of the branch decomposition
1012  whichBranch,
1013  branchMinimum,
1014  branchMaximum,
1015  branchSaddle,
1016  branchParent,
1017  bestUpward,
1018  bestDownward);
1019 
1020  } // ComputeHeightBranchDecomposition()
1021 
1022  std::vector<Id> static findSuperPathToRoot(
1024  vtkm::Id vertex)
1025  {
1026  // Initialise the empty path and starting vertex
1027  std::vector<vtkm::Id> path;
1028  vtkm::Id current = vertex;
1029 
1030  // Go up the parent list until we reach the root
1031  while (MaskedIndex(parentsPortal.Get(current)) != 0)
1032  {
1033  path.push_back(current);
1034  current = MaskedIndex(parentsPortal.Get(current));
1035  }
1036  path.push_back(current);
1037 
1038  return path;
1039  }
1040 
1041  // Given a path from a leaf (the global min/max) to the root of the contour tree and a hypersweep (where all hyperarcs are cut at the path)
1042  // This function performs a prefix scan along that path to obtain the correct hypersweep values (as in the global min/max is the root of the hypersweep)
1043  void static fixPath(const std::function<vtkm::Id(vtkm::Id, vtkm::Id)> operation,
1044  const std::vector<vtkm::Id> path,
1046  {
1048 
1049  // Fix path from the old root to the new root. Parallelisble with a prefix scan, but sufficiently fast for now.
1050  for (auto i = path.size() - 2; i > 0; i--)
1051  {
1052  const auto vertex = path[i + 1];
1053  const auto parent = path[i];
1054 
1055  const auto vertexValue = minMaxIndex.Get(vertex);
1056  const auto parentValue = minMaxIndex.Get(parent);
1057 
1058  minMaxIndex.Set(parent, operation(vertexValue, parentValue));
1059  }
1060  }
1061 
1062  // This function edits all the hyperarcs which contain vertices which are on the supplied path. This path is usually the path between the global min/max to the root of the tree.
1063  // This function effectively cuts hyperarcs at the first node they encounter along that path.
1064  // In addition to this it computed the number of supernodes every hyperarc has on that path. This helps in the function hyperarcScan for choosing the new target of the cut hyperarcs.
1065  // NOTE: It is assumed that the supplied path starts at a leaf and ends at the root of the contour tree.
1066  void static editHyperarcs(
1067  const vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType hyperparentsPortal,
1068  const std::vector<vtkm::Id> path,
1071  {
1073 
1074  std::size_t i = 0;
1075  while (i < path.size())
1076  {
1077  // Cut the hyperacs at the first point
1078  hyperarcsPortal.Set(MaskedIndex(hyperparentsPortal.Get(path[i])), path[i]);
1079 
1080  Id currentHyperparent = MaskedIndex(hyperparentsPortal.Get(path[i]));
1081 
1082  // Skip the rest of the supernodes which are on the same hyperarc
1083  while (i < path.size() && MaskedIndex(hyperparentsPortal.Get(path[i])) == currentHyperparent)
1084  {
1085  const auto value = howManyUsedPortal.Get(MaskedIndex(hyperparentsPortal.Get(path[i])));
1086  howManyUsedPortal.Set(MaskedIndex(hyperparentsPortal.Get(path[i])), value + 1);
1087  i++;
1088  }
1089  }
1090  }
1091 
1092  template <class BinaryFunctor>
1093  void static hyperarcScan(const vtkm::cont::ArrayHandle<vtkm::Id> supernodes,
1094  const vtkm::cont::ArrayHandle<vtkm::Id> hypernodes,
1095  const vtkm::cont::ArrayHandle<vtkm::Id> hyperarcs,
1096  const vtkm::cont::ArrayHandle<vtkm::Id> hyperparents,
1097  const vtkm::cont::ArrayHandle<vtkm::Id> hyperparentKeys,
1098  const vtkm::cont::ArrayHandle<vtkm::Id> whenTransferred,
1099  const vtkm::cont::ArrayHandle<vtkm::Id> howManyUsed,
1100  const vtkm::Id nIterations,
1101  const BinaryFunctor operation,
1103  {
1105 
1106  vtkm::cont::Invoker invoke;
1107 
1108  auto supernodesPortal = supernodes.ReadPortal();
1109  auto hypernodesPortal = hypernodes.ReadPortal();
1110  auto hyperparentsPortal = hyperparents.ReadPortal();
1111 
1112  // Set the first supernode per iteration
1113  vtkm::cont::ArrayHandle<vtkm::Id> firstSupernodePerIteration;
1115  firstSupernodePerIteration);
1116 
1117  // The first different from the previous is the first in the iteration
1119  setFirstSupernodePerIteration;
1120  invoke(setFirstSupernodePerIteration, whenTransferred, firstSupernodePerIteration);
1121 
1122  auto firstSupernodePerIterationPortal = firstSupernodePerIteration.WritePortal();
1123  for (vtkm::Id iteration = 1; iteration < nIterations; ++iteration)
1124  {
1125  if (firstSupernodePerIterationPortal.Get(iteration) == 0)
1126  {
1127  firstSupernodePerIterationPortal.Set(iteration,
1128  firstSupernodePerIterationPortal.Get(iteration + 1));
1129  }
1130  }
1131 
1132  // set the sentinel at the end of the array
1133  firstSupernodePerIterationPortal.Set(nIterations, supernodesPortal.GetNumberOfValues());
1134 
1135  // Set the first hypernode per iteration
1136  vtkm::cont::ArrayHandle<vtkm::Id> firstHypernodePerIteration;
1138  firstHypernodePerIteration);
1139  auto firstHypernodePerIterationPortal = firstHypernodePerIteration.WritePortal();
1140 
1141  for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
1142  {
1143  firstHypernodePerIterationPortal.Set(
1144  iteration, hyperparentsPortal.Get(firstSupernodePerIterationPortal.Get(iteration)));
1145  }
1146 
1147  // Set the sentinel at the end of the array
1148  firstHypernodePerIterationPortal.Set(nIterations, hypernodesPortal.GetNumberOfValues());
1149 
1150  // This workled is used in every iteration of the following loop, so it's initialised outside.
1152  BinaryFunctor>
1153  addDependentWeightHypersweepWorklet(operation);
1154 
1155  // For every iteration do a prefix scan on all hyperarcs in that iteration and then transfer the scanned value to every hyperarc's target supernode
1156  for (vtkm::Id iteration = 0; iteration < nIterations; iteration++)
1157  {
1158  // Determine the first and last hypernode in the current iteration (all hypernodes between them are also in the current iteration)
1159  vtkm::Id firstHypernode = firstHypernodePerIterationPortal.Get(iteration);
1160  vtkm::Id lastHypernode = firstHypernodePerIterationPortal.Get(iteration + 1);
1161  lastHypernode = vtkm::Minimum()(lastHypernode, hypernodes.GetNumberOfValues() - 1);
1162 
1163  // Determine the first and last supernode in the current iteration (all supernode between them are also in the current iteration)
1164  vtkm::Id firstSupernode = MaskedIndex(hypernodesPortal.Get(firstHypernode));
1165  vtkm::Id lastSupernode = MaskedIndex(hypernodesPortal.Get(lastHypernode));
1166  lastSupernode = vtkm::Minimum()(lastSupernode, hyperparents.GetNumberOfValues() - 1);
1167 
1168  // Prefix scan along all hyperarcs in the current iteration
1169  auto subarrayValues = vtkm::cont::make_ArrayHandleView(
1170  minMaxIndex, firstSupernode, lastSupernode - firstSupernode);
1171  auto subarrayKeys = vtkm::cont::make_ArrayHandleView(
1172  hyperparentKeys, firstSupernode, lastSupernode - firstSupernode);
1174  subarrayKeys, subarrayValues, subarrayValues, operation);
1175 
1176  // Array containing the Ids of the hyperarcs in the current iteration
1177  vtkm::cont::ArrayHandleCounting<vtkm::Id> iterationHyperarcs(
1178  firstHypernode, 1, lastHypernode - firstHypernode);
1179 
1180  // Transfer the value accumulated in the last entry of the prefix scan to the hypernode's targe supernode
1181  invoke(addDependentWeightHypersweepWorklet,
1182  iterationHyperarcs,
1183  hypernodes,
1184  hyperarcs,
1185  howManyUsed,
1186  minMaxIndex);
1187  }
1188  }
1189 }; // class ProcessContourTree
1190 } // namespace contourtree_augmented
1191 } // worklet
1192 } // vtkm
1193 
1194 #endif
vtkm::cont::ArrayHandle::GetNumberOfValues
VTKM_CONT vtkm::Id GetNumberOfValues() const
Returns the number of entries in the array.
Definition: ArrayHandle.h:448
vtkm::worklet::contourtree_augmented::ContourTree::Hyperparents
IdArrayType Hyperparents
Definition: augmented/ContourTree.h:136
vtkm::worklet::contourtree_augmented::ProcessContourTree::ProcessContourTree
ProcessContourTree()
Definition: ProcessContourTree.h:105
vtkm::cont::make_ArrayHandle
VTKM_CONT vtkm::cont::ArrayHandleBasic< T > make_ArrayHandle(const T *array, vtkm::Id numberOfValues, vtkm::CopyFlag copy)
A convenience function for creating an ArrayHandle from a standard C array.
Definition: ArrayHandleBasic.h:217
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeBranchDecomposition
static process_contourtree_inc_ns::Branch< T > * ComputeBranchDecomposition(const IdArrayType &contourTreeSuperparents, const IdArrayType &contourTreeSupernodes, const IdArrayType &whichBranch, const IdArrayType &branchMinimum, const IdArrayType &branchMaximum, const IdArrayType &branchSaddle, const IdArrayType &branchParent, const IdArrayType &sortOrder, const vtkm::cont::ArrayHandle< T, StorageType > &dataField, bool dataFieldIsSorted)
Definition: ProcessContourTree.h:709
vtkm::worklet::contourtree_augmented::process_contourtree_inc::WhichBranchNewId
Definition: HypersweepWorklets.h:497
vtkm::cont::ArrayHandle< vtkm::Id >
vtkm::worklet::contourtree_augmented::IsAscending
VTKM_EXEC_CONT bool IsAscending(vtkm::Id flaggedIndex)
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:121
vtkm::worklet::contourtree_augmented::ProcessContourTree::fixPath
static void fixPath(const std::function< vtkm::Id(vtkm::Id, vtkm::Id)> operation, const std::vector< vtkm::Id > path, vtkm::cont::ArrayHandle< vtkm::Id >::WritePortalType minMaxIndex)
Definition: ProcessContourTree.h:1043
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetFirstSupernodePerIteration
Definition: HypersweepWorklets.h:181
vtkm::worklet::contourtree_augmented::ContourTree::Hyperarcs
IdArrayType Hyperarcs
Definition: augmented/ContourTree.h:149
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeHeightBranchDecomposition
static void ComputeHeightBranchDecomposition(const ContourTree &contourTree, const cont::ArrayHandle< Float64 > fieldValues, const IdArrayType &ctSortOrder, const vtkm::Id nIterations, IdArrayType &whichBranch, IdArrayType &branchMinimum, IdArrayType &branchMaximum, IdArrayType &branchSaddle, IdArrayType &branchParent)
Definition: ProcessContourTree.h:829
vtkm::cont::Algorithm::ScanInclusive
static VTKM_CONT T ScanInclusive(vtkm::cont::DeviceAdapterId devId, const vtkm::cont::ArrayHandle< T, CIn > &input, vtkm::cont::ArrayHandle< T, COut > &output)
Definition: Algorithm.h:731
vtkm::worklet::contourtree_augmented::process_contourtree_inc::InitialiseArcsVolume
Definition: HypersweepWorklets.h:72
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::cont::make_ArrayHandleView
ArrayHandleView< ArrayHandleType > make_ArrayHandleView(const ArrayHandleType &array, vtkm::Id startIndex, vtkm::Id numValues)
Definition: ArrayHandleView.h:222
vtkm::worklet::contourtree_augmented::ContourTree::Supernodes
IdArrayType Supernodes
Definition: augmented/ContourTree.h:125
Types.h
vtkm::worklet::contourtree_augmented::ContourTree
Definition: augmented/ContourTree.h:106
vtkm::cont::ArrayHandle::Allocate
VTKM_CONT void Allocate(vtkm::Id numberOfValues, vtkm::CopyFlag preserve, vtkm::cont::Token &token) const
Allocates an array large enough to hold the given number of values.
Definition: ArrayHandle.h:465
Pair.h
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeBranchData
static void ComputeBranchData(const ContourTree &contourTree, IdArrayType &whichBranch, IdArrayType &branchMinimum, IdArrayType &branchMaximum, IdArrayType &branchSaddle, IdArrayType &branchParent, IdArrayType &bestUpward, IdArrayType &bestDownward)
Definition: ProcessContourTree.h:554
vtkm::worklet::contourtree_augmented::ContourTree::Superparents
IdArrayType Superparents
Definition: augmented/ContourTree.h:118
vtkm::worklet::contourtree_augmented::ContourTree::Hypernodes
IdArrayType Hypernodes
Definition: augmented/ContourTree.h:144
vtkm::cont::Algorithm::ScanInclusiveByKey
static VTKM_CONT void ScanInclusiveByKey(vtkm::cont::DeviceAdapterId devId, const vtkm::cont::ArrayHandle< T, KIn > &keys, const vtkm::cont::ArrayHandle< U, VIn > &values, vtkm::cont::ArrayHandle< U, VOut > &values_output, BinaryFunctor binary_functor)
Definition: Algorithm.h:772
vtkm::worklet::contourtree_augmented::ProcessContourTree::findSuperPathToRoot
static std::vector< Id > findSuperPathToRoot(vtkm::cont::ArrayHandle< vtkm::Id >::ReadPortalType parentsPortal, vtkm::Id vertex)
Definition: ProcessContourTree.h:1022
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetBestUpDown
Definition: HypersweepWorklets.h:400
vtkm::worklet::contourtree_augmented::ContourTree::Nodes
IdArrayType Nodes
Definition: augmented/ContourTree.h:112
vtkm::cont::Algorithm::Sort
static VTKM_CONT void Sort(vtkm::cont::DeviceAdapterId devId, vtkm::cont::ArrayHandle< T, Storage > &values)
Definition: Algorithm.h:965
ArrayHandleConstant.h
vtkm::Maximum
Binary Predicate that takes two arguments argument x, and y and returns the x if x > y otherwise retu...
Definition: BinaryOperators.h:85
SuperNodeBranchComparator.h
PrintVectors.h
Invoker.h
ArrayTransforms.h
vtkm::worklet::contourtree_augmented::process_contourtree_inc::BranchSaddleParentSet
Definition: HypersweepWorklets.h:567
ArrayHandleView.h
vtkm::worklet::contourtree_augmented::MaskedIndex
VTKM_EXEC_CONT vtkm::Id MaskedIndex(vtkm::Id flaggedIndex)
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:127
vtkm::cont::ArrayHandle::ReadPortalType
typename StorageType::ReadPortalType ReadPortalType
Definition: ArrayHandle.h:294
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeIntrinsicWeight
Definition: HypersweepWorklets.h:150
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::worklet::contourtree_augmented::ProcessContourTree
Definition: ProcessContourTree.h:101
ArrayCopy.h
PointerDoubling.h
vtkm::worklet::contourtree_augmented::ProcessContourTree::hyperarcScan
static void hyperarcScan(const vtkm::cont::ArrayHandle< vtkm::Id > supernodes, const vtkm::cont::ArrayHandle< vtkm::Id > hypernodes, const vtkm::cont::ArrayHandle< vtkm::Id > hyperarcs, const vtkm::cont::ArrayHandle< vtkm::Id > hyperparents, const vtkm::cont::ArrayHandle< vtkm::Id > hyperparentKeys, const vtkm::cont::ArrayHandle< vtkm::Id > whenTransferred, const vtkm::cont::ArrayHandle< vtkm::Id > howManyUsed, const vtkm::Id nIterations, const BinaryFunctor operation, vtkm::cont::ArrayHandle< vtkm::Id > minMaxIndex)
Definition: ProcessContourTree.h:1093
vtkm::worklet::contourtree_augmented::ContourTree::Superarcs
IdArrayType Superarcs
Definition: augmented/ContourTree.h:129
vtkm::Sum
Binary Predicate that takes two arguments argument x, and y and returns sum (addition) of the two val...
Definition: BinaryOperators.h:33
vtkm::SortLess
Binary Predicate that takes two arguments argument x, and y and returns True if and only if x is less...
Definition: BinaryPredicates.h:45
vtkm::worklet::contourtree_augmented::NoSuchElement
VTKM_EXEC_CONT bool NoSuchElement(vtkm::Id flaggedIndex)
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:97
vtkm::worklet::contourtree_augmented::process_contourtree_inc::FinaliseChainToBranch
Definition: HypersweepWorklets.h:642
vtkm::worklet::contourtree_augmented::process_contourtree_inc::AddDependentWeightHypersweep
Definition: HypersweepWorklets.h:212
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PointerDoubling
Definition: processcontourtree/PointerDoubling.h:75
Timer.h
Algorithm.h
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeVolumeWeightsSerial
static void ComputeVolumeWeightsSerial(const ContourTree &contourTree, const vtkm::Id nIterations, IdArrayType &superarcIntrinsicWeight, IdArrayType &superarcDependentWeight, IdArrayType &supernodeTransferWeight, IdArrayType &hyperarcDependentWeight)
Definition: ProcessContourTree.h:211
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SuperArcVolumetricComparator
Definition: SuperArcVolumetricComparator.h:153
vtkm::cont::ArrayHandleCounting< vtkm::Id >
vtkm::worklet::contourtree_augmented::ProcessContourTree::CollectSortedSuperarcs
static void CollectSortedSuperarcs(const ContourTree &contourTree, const IdArrayType &sortOrder, EdgePairArray &saddlePeak)
Definition: ProcessContourTree.h:154
vtkm::cont::Invoker
Allows launching any worklet without a dispatcher.
Definition: Invoker.h:41
vtkm::Pair::first
FirstType first
The pair's first object.
Definition: Pair.h:50
vtkm::worklet::contourtree_augmented::process_contourtree_inc
Definition: Branch.h:69
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeVolumeBranchDecompositionSerial
static void ComputeVolumeBranchDecompositionSerial(const ContourTree &contourTree, const IdArrayType &superarcDependentWeight, const IdArrayType &superarcIntrinsicWeight, IdArrayType &whichBranch, IdArrayType &branchMinimum, IdArrayType &branchMaximum, IdArrayType &branchSaddle, IdArrayType &branchParent)
Definition: ProcessContourTree.h:399
vtkm::worklet::contourtree_augmented::process_contourtree_inc::Branch::ComputeBranchDecomposition
static Branch< T > * ComputeBranchDecomposition(const IdArrayType &contourTreeSuperparents, const IdArrayType &contourTreeSupernodes, const IdArrayType &whichBranch, const IdArrayType &branchMinimum, const IdArrayType &branchMaximum, const IdArrayType &branchSaddle, const IdArrayType &branchParent, const IdArrayType &sortOrder, const vtkm::cont::ArrayHandle< T, StorageType > &dataField, bool dataFieldIsSorted)
Definition: Branch.h:151
vtkm::cont::ArrayCopy
void ArrayCopy(const SourceArrayType &source, DestArrayType &destination)
Does a deep copy from one array to another array.
Definition: ArrayCopy.h:142
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SetFirstVertexForSuperparent
Definition: HypersweepWorklets.h:121
HypersweepWorklets.h
Types.h
vtkm::worklet::contourtree_augmented::process_contourtree_inc::BranchMinMaxSet
Definition: HypersweepWorklets.h:519
vtkm::cont::ArrayHandle::WritePortal
VTKM_CONT WritePortalType WritePortal() const
Get an array portal that can be used in the control environment.
Definition: ArrayHandle.h:435
vtkm::cont::ArrayHandle::WritePortalType
typename StorageType::WritePortalType WritePortalType
Definition: ArrayHandle.h:295
vtkm::worklet::contourtree_augmented::process_contourtree_inc::ComputeSubtreeHeight
Definition: HypersweepWorklets.h:361
vtkm::worklet::contourtree_augmented::process_contourtree_inc::SuperNodeBranchComparator
Definition: SuperNodeBranchComparator.h:116
vtkm::worklet::contourtree_augmented::ProcessContourTree::ComputeVolumeBranchDecomposition
static void ComputeVolumeBranchDecomposition(const ContourTree &contourTree, const vtkm::Id nIterations, IdArrayType &whichBranch, IdArrayType &branchMinimum, IdArrayType &branchMaximum, IdArrayType &branchSaddle, IdArrayType &branchParent)
Definition: ProcessContourTree.h:735
vtkm::cont::ArrayHandleConstant
An array handle with a constant value.
Definition: ArrayHandleConstant.h:63
vtkm::CopyFlag::On
@ On
Branch.h
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PropagateBestUpDown
Definition: HypersweepWorklets.h:466
vtkm::worklet::contourtree_augmented::ProcessContourTree::editHyperarcs
static void editHyperarcs(const vtkm::cont::ArrayHandle< vtkm::Id >::ReadPortalType hyperparentsPortal, const std::vector< vtkm::Id > path, vtkm::cont::ArrayHandle< vtkm::Id >::WritePortalType hyperarcsPortal, vtkm::cont::ArrayHandle< vtkm::Id >::WritePortalType howManyUsedPortal)
Definition: ProcessContourTree.h:1066
vtkm::worklet::contourtree_augmented::process_contourtree_inc::InitialiseArcs
Definition: HypersweepWorklets.h:275
SuperArcVolumetricComparator.h
vtkm::cont::ArrayHandle::ReadPortal
VTKM_CONT ReadPortalType ReadPortal() const
Get an array portal that can be used in the control environment.
Definition: ArrayHandle.h:414
vtkm::worklet::contourtree_augmented::ContourTree::Arcs
IdArrayType Arcs
Definition: augmented/ContourTree.h:115
vtkm::worklet::contourtree_augmented::process_contourtree_inc::IncorporateParent
Definition: HypersweepWorklets.h:672
vtkm::worklet::contourtree_augmented::ContourTree::FirstHypernodePerIteration
IdArrayType FirstHypernodePerIteration
Definition: augmented/ContourTree.h:157
vtkm::worklet::contourtree_augmented::EdgePair
vtkm::Pair< vtkm::Id, vtkm::Id > EdgePair
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:92
vtkm::worklet::contourtree_augmented::PrintIndices
void PrintIndices(std::string label, const vtkm::cont::ArrayHandle< T > &iVec, vtkm::Id nIndices=-1, std::ostream &outStream=std::cout)
Definition: augmented/PrintVectors.h:253
vtkm::worklet::contourtree_augmented::PrintHeader
void PrintHeader(vtkm::Id howMany, std::ostream &outStream=std::cout)
Definition: augmented/PrintVectors.h:151
vtkm::worklet::contourtree_augmented::process_contourtree_inc::Branch
Definition: Branch.h:74
vtkm::worklet::contourtree_augmented::ProcessContourTree::CollectSortedArcs
static void CollectSortedArcs(const ContourTree &contourTree, const IdArrayType &sortOrder, EdgePairArray &sortedArcs)
Definition: ProcessContourTree.h:110
vtkm::worklet::contourtree_augmented::SaddlePeakSort
Definition: augmented/ContourTree.h:87
ContourTree.h
vtkm::Pair
A vtkm::Pair is essentially the same as an STL pair object except that the methods (constructors and ...
Definition: Pair.h:29
vtkm::worklet::contourtree_augmented::ContourTree::FirstSupernodePerIteration
IdArrayType FirstSupernodePerIteration
Definition: augmented/ContourTree.h:156
vtkm::worklet::contourtree_augmented::process_contourtree_inc::UnmaskArray
Definition: HypersweepWorklets.h:445
vtkm::worklet::contourtree_augmented::NO_SUCH_ELEMENT
constexpr vtkm::Id NO_SUCH_ELEMENT
Definition: filter/scalar_topology/worklet/contourtree_augmented/Types.h:73
vtkm::Pair::second
SecondType second
The pair's second object.
Definition: Pair.h:55
vtkm::worklet::contourtree_augmented::process_contourtree_inc::PrepareChainToBranch
Definition: HypersweepWorklets.h:616
vtkm::cont::ArrayHandleIndex
An implicit array handle containing the its own indices.
Definition: ArrayHandleIndex.h:54
vtkm::Minimum
Binary Predicate that takes two arguments argument x, and y and returns the x if x < y otherwise retu...
Definition: BinaryOperators.h:99
vtkm::worklet::contourtree_augmented::ContourTree::WhenTransferred
IdArrayType WhenTransferred
Definition: augmented/ContourTree.h:139