diff --git a/examples/contour_tree_distributed/ContourTreeApp.cxx b/examples/contour_tree_distributed/ContourTreeApp.cxx index b00059ba5..b3ff96568 100644 --- a/examples/contour_tree_distributed/ContourTreeApp.cxx +++ b/examples/contour_tree_distributed/ContourTreeApp.cxx @@ -74,6 +74,7 @@ #include #include +#include #include #include #include @@ -213,8 +214,8 @@ int main(int argc, char* argv[]) if (!augmentHierarchicalTree) { VTKM_LOG_S(vtkm::cont::LogLevel::Warn, - "Warning: --computeVolumeBranchDecomposition only " - "allowed augmentation. Enabling --augmentHierarchicalTree option."); + "Warning: --computeVolumeBranchDecomposition requires augmentation. " + "Enabling --augmentHierarchicalTree option."); augmentHierarchicalTree = true; } } @@ -265,6 +266,31 @@ int main(int argc, char* argv[]) } } + vtkm::Id numBranches = 0; + if (parser.hasOption("--numBranches")) + { + numBranches = std::stoi(parser.getOption("--numBranches")); + if (!computeHierarchicalVolumetricBranchDecomposition) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Warn, + "Warning: --numBranches requires computing branch decomposition. " + "Enabling --computeHierarchicalVolumetricBranchDecomposition option."); + computeHierarchicalVolumetricBranchDecomposition = true; + } + if (!augmentHierarchicalTree) + { + VTKM_LOG_S(vtkm::cont::LogLevel::Warn, + "Warning: --numBranches requires augmentation. " + "Enabling --augmentHierarchicalTree option."); + augmentHierarchicalTree = true; + } + } + + ValueType eps = 0.00001f; + + if (parser.hasOption("--eps")) + eps = std::stof(parser.getOption("--eps")); + #ifdef ENABLE_HDFIO std::string dataset_name = "data"; if (parser.hasOption("--dataset")) @@ -336,7 +362,12 @@ int main(int argc, char* argv[]) std::cout << "--augmentHierarchicalTree Augment the hierarchical tree." << std::endl; std::cout << "--computeVolumeBranchDecomposition Compute the volume branch decomposition. " << std::endl; - std::cout << " Requries --augmentHierarchicalTree to be set." << std::endl; + std::cout << " Requires --augmentHierarchicalTree to be set." << std::endl; + std::cout << "--numBranches Number of top volume branches to select." << std::endl; + std::cout << " Requires --computeVolumeBranchDecomposition." << std::endl; + std::cout + << "--eps= Floating point offset awary from the critical point. (default=0.00001)" + << std::endl; std::cout << "--preSplitFiles Input data is already pre-split into blocks." << std::endl; std::cout << "--saveDot Save DOT files of the distributed contour tree " << std::endl << " computation (Default=False). " << std::endl; @@ -383,6 +414,8 @@ int main(int argc, char* argv[]) << " saveOutputData=" << saveOutputData << std::endl << " forwardSummary=" << forwardSummary << std::endl << " nblocks=" << numBlocks << std::endl + << " nbranches=" << numBranches << std::endl + << " eps=" << eps << std::endl #ifdef ENABLE_HDFIO << " dataset=" << dataset_name << " (HDF5 only)" << std::endl << " blocksPerDim=" << blocksPerDimIn[0] << "," << blocksPerDimIn[1] << "," @@ -642,6 +675,15 @@ int main(int argc, char* argv[]) vtkm::Float64 branchDecompTime = currTime - prevTime; prevTime = currTime; + // Compute SelectTopVolumeContours if needed + vtkm::cont::PartitionedDataSet tp_result; + if (numBranches > 0) + { + vtkm::filter::scalar_topology::SelectTopVolumeContoursFilter tp_filter; + tp_filter.SetSavedBranches(numBranches); + tp_result = tp_filter.Execute(bd_result); + } + // Save output if (saveOutputData) { @@ -684,6 +726,66 @@ int main(int argc, char* argv[]) << lowerEndGRId.Get(branch) << std::endl; } } + + if (numBranches > 0) + { +#ifndef DEBUG_PRINT + bool print_to_files = (rank == 0); + vtkm::Id max_blocks_to_print = 1; +#else + bool print_to_files = true; + vtkm::Id max_blocks_to_print = result.GetNumberOfPartitions(); +#endif + + for (vtkm::Id ds_no = 0; print_to_files && ds_no < max_blocks_to_print; ++ds_no) + { + auto ds = tp_result.GetPartition(ds_no); + std::string topVolumeBranchFileName = std::string("TopVolumeBranch_Rank_") + + std::to_string(static_cast(rank)) + std::string("_Block_") + + std::to_string(static_cast(ds_no)) + std::string(".txt"); + std::ofstream topVolumeBranchStream(topVolumeBranchFileName.c_str()); + auto topVolBranchGRId = ds.GetField("TopVolumeBranchGlobalRegularIds") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchVolume = ds.GetField("TopVolumeBranchVolume") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchSaddleEpsilon = ds.GetField("TopVolumeBranchSaddleEpsilon") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchSaddleIsoValue = ds.GetField("TopVolumeBranchSaddleIsoValue") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + + vtkm::Id nSelectedBranches = topVolBranchGRId.GetNumberOfValues(); + for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch) + { + topVolumeBranchStream << std::setw(12) << topVolBranchGRId.Get(branch) + << std::setw(14) << topVolBranchVolume.Get(branch) + << std::setw(5) << topVolBranchSaddleEpsilon.Get(branch) + << std::setw(14) << topVolBranchSaddleIsoValue.Get(branch) + << std::endl; + } + + std::string isoValuesFileName = std::string("IsoValues_Rank_") + + std::to_string(static_cast(rank)) + std::string("_Block_") + + std::to_string(static_cast(ds_no)) + std::string(".txt"); + std::ofstream isoValuesStream(isoValuesFileName.c_str()); + + for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch) + { + isoValuesStream << (topVolBranchSaddleIsoValue.Get(branch) + + (eps * topVolBranchSaddleEpsilon.Get(branch))) + << " "; + } + + isoValuesStream << std::endl; + } + } } else { diff --git a/vtkm/filter/scalar_topology/CMakeLists.txt b/vtkm/filter/scalar_topology/CMakeLists.txt index d3eec065a..5603ebaa1 100644 --- a/vtkm/filter/scalar_topology/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/CMakeLists.txt @@ -12,17 +12,21 @@ set(scalar_topology_headers ContourTreeUniformAugmented.h ContourTreeUniformDistributed.h DistributedBranchDecompositionFilter.h + SelectTopVolumeContoursFilter.h ) set(scalar_topology_sources internal/BranchDecompositionBlock.cxx + internal/SelectTopVolumeContoursBlock.cxx internal/ComputeBlockIndices.cxx internal/ComputeDistributedBranchDecompositionFunctor.cxx + internal/SelectTopVolumeContoursFunctor.cxx internal/ExchangeBranchEndsFunctor.cxx ContourTreeUniform.cxx ContourTreeUniformAugmented.cxx ContourTreeUniformDistributed.cxx DistributedBranchDecompositionFilter.cxx + SelectTopVolumeContoursFilter.cxx ) vtkm_library( diff --git a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx index 03714050a..9992fd64a 100644 --- a/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx +++ b/vtkm/filter/scalar_topology/ContourTreeUniformDistributed.cxx @@ -832,6 +832,15 @@ VTKM_CONT void ContourTreeUniformDistributed::DoPostExecute( "from information in CellSetStructured."); diyBounds = vtkm::filter::scalar_topology::internal::ComputeBlockIndices( input, diyDivisions, vtkmdiyLocalBlockGids); + + // Set BlocksPerDimension fromn diyDivisions result as add them + // as information to the output data set for use in subsequent + // filters + this->BlocksPerDimension = vtkm::Id3{ 1, 1, 1 }; + for (unsigned int d = 0; d < diyDivisions.size(); ++d) + { + this->BlocksPerDimension[d] = diyDivisions[d]; + } } else { diff --git a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx index a7000420f..0f522424b 100644 --- a/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx +++ b/vtkm/filter/scalar_topology/DistributedBranchDecompositionFilter.cxx @@ -331,7 +331,6 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D outputDataSets[ds_no] = input.GetPartition(ds_no); } - branch_decomposition_master.foreach ( [&](BranchDecompositionBlock* b, const vtkmdiy::Master::ProxyWithLink&) { vtkm::cont::Field branchRootField( @@ -347,6 +346,51 @@ VTKM_CONT vtkm::cont::PartitionedDataSet DistributedBranchDecompositionFilter::D vtkm::cont::Field::Association::WholeDataSet, b->VolumetricBranchDecomposer.LowerEndGRId); outputDataSets[b->LocalBlockNo].AddField(LowerEndGRIdField); + + vtkm::cont::Field UpperEndIntrinsicVolume( + "UpperEndIntrinsicVolume", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.UpperEndIntrinsicVolume); + outputDataSets[b->LocalBlockNo].AddField(UpperEndIntrinsicVolume); + vtkm::cont::Field UpperEndDependentVolume( + "UpperEndDependentVolume", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.UpperEndDependentVolume); + outputDataSets[b->LocalBlockNo].AddField(UpperEndDependentVolume); + vtkm::cont::Field LowerEndIntrinsicVolume( + "LowerEndIntrinsicVolume", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.LowerEndIntrinsicVolume); + outputDataSets[b->LocalBlockNo].AddField(LowerEndIntrinsicVolume); + vtkm::cont::Field LowerEndDependentVolume( + "LowerEndDependentVolume", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.LowerEndDependentVolume); + outputDataSets[b->LocalBlockNo].AddField(LowerEndDependentVolume); + vtkm::cont::Field LowerEndSuperarcId("LowerEndSuperarcId", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.LowerEndSuperarcId); + outputDataSets[b->LocalBlockNo].AddField(LowerEndSuperarcId); + vtkm::cont::Field UpperEndSuperarcId("UpperEndSuperarcId", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.UpperEndSuperarcId); + outputDataSets[b->LocalBlockNo].AddField(UpperEndSuperarcId); + vtkm::cont::Field LowerEndValue("LowerEndValue", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.LowerEndValue); + outputDataSets[b->LocalBlockNo].AddField(LowerEndValue); + vtkm::cont::Field UpperEndValue("UpperEndValue", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.UpperEndValue); + outputDataSets[b->LocalBlockNo].AddField(UpperEndValue); + vtkm::cont::Field BranchRoot("BranchRoot", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.BranchRoot); + outputDataSets[b->LocalBlockNo].AddField(BranchRoot); + vtkm::cont::Field BranchRootGRId("BranchRootGRId", + vtkm::cont::Field::Association::WholeDataSet, + b->VolumetricBranchDecomposer.BranchRootGRId); + outputDataSets[b->LocalBlockNo].AddField(BranchRootGRId); }); timingsStream << " " << std::setw(38) << std::left diff --git a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx new file mode 100644 index 000000000..c7d0b8f31 --- /dev/null +++ b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.cxx @@ -0,0 +1,183 @@ +//============================================================================ +// 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. +//============================================================================ + +#include +#include +#include +#include +#include +#include + + +// vtkm includes +#include + +// DIY includes +// clang-format off +VTKM_THIRDPARTY_PRE_INCLUDE +#include +#include +VTKM_THIRDPARTY_POST_INCLUDE +// clang-format on + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ + +VTKM_CONT vtkm::cont::DataSet SelectTopVolumeContoursFilter::DoExecute(const vtkm::cont::DataSet&) +{ + throw vtkm::cont::ErrorFilterExecution( + "SelectTopVolumeContoursFilter expects PartitionedDataSet as input."); +} + +VTKM_CONT vtkm::cont::PartitionedDataSet SelectTopVolumeContoursFilter::DoExecutePartitions( + const vtkm::cont::PartitionedDataSet& input) +{ + auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator(); + int rank = comm.rank(); + int size = comm.size(); + + using SelectTopVolumeContoursBlock = + vtkm::filter::scalar_topology::internal::SelectTopVolumeContoursBlock; + vtkmdiy::Master branch_top_volume_master(comm, + 1, // Use 1 thread, VTK-M will do the treading + -1, // All blocks in memory + 0, // No create function + SelectTopVolumeContoursBlock::Destroy); + + auto firstDS = input.GetPartition(0); + vtkm::Id3 firstPointDimensions, firstGlobalPointDimensions, firstGlobalPointIndexStart; + firstDS.GetCellSet().CastAndCallForTypes( + vtkm::worklet::contourtree_augmented::GetLocalAndGlobalPointDimensions(), + firstPointDimensions, + firstGlobalPointDimensions, + firstGlobalPointIndexStart); + int numDims = firstGlobalPointDimensions[2] > 1 ? 3 : 2; + auto vtkmBlocksPerDimensionRP = input.GetPartition(0) + .GetField("vtkmBlocksPerDimension") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + + int globalNumberOfBlocks = 1; + + for (vtkm::IdComponent d = 0; d < static_cast(numDims); ++d) + { + globalNumberOfBlocks *= static_cast(vtkmBlocksPerDimensionRP.Get(d)); + } + + vtkmdiy::DynamicAssigner assigner(comm, size, globalNumberOfBlocks); + for (vtkm::Id localBlockIndex = 0; localBlockIndex < input.GetNumberOfPartitions(); + ++localBlockIndex) + { + const vtkm::cont::DataSet& ds = input.GetPartition(localBlockIndex); + int globalBlockId = static_cast( + vtkm::cont::ArrayGetValue(0, + ds.GetField("vtkmGlobalBlockId") + .GetData() + .AsArrayHandle>())); + + SelectTopVolumeContoursBlock* b = + new SelectTopVolumeContoursBlock(localBlockIndex, globalBlockId); + + branch_top_volume_master.add(globalBlockId, b, new vtkmdiy::Link()); + assigner.set_rank(rank, globalBlockId); + } + + vtkmdiy::fix_links(branch_top_volume_master, assigner); + + branch_top_volume_master.foreach ( + [&](SelectTopVolumeContoursBlock* b, const vtkmdiy::Master::ProxyWithLink&) { + const auto& globalSize = firstGlobalPointDimensions; + vtkm::Id totalVolume = globalSize[0] * globalSize[1] * globalSize[2]; + const vtkm::cont::DataSet& ds = input.GetPartition(b->LocalBlockNo); + + b->SortBranchByVolume(ds, totalVolume); + + // copy the top volume branches into a smaller array + // we skip index 0 because it must be the main branch (which has the highest volume) + vtkm::Id nActualSavedBranches = + std::min(this->nSavedBranches, b->SortedBranchByVolume.GetNumberOfValues() - 1); + + vtkm::worklet::contourtree_augmented::IdArrayType topVolumeBranch; + vtkm::cont::Algorithm::CopySubRange( + b->SortedBranchByVolume, 1, nActualSavedBranches, topVolumeBranch); + + auto branchRootGRId = + ds.GetField("BranchRootGRId").GetData().AsArrayHandle>(); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchRootGRId, topVolumeBranch, b->TopVolumeBranchRootGRId); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + branchRootGRId, topVolumeBranch, b->TopVolumeBranchRootGRId); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + b->BranchVolume, topVolumeBranch, b->TopVolumeBranchVolume); + + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + b->BranchSaddleEpsilon, topVolumeBranch, b->TopVolumeBranchSaddleEpsilon); + + auto resolveArray = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + InArrayHandleType topVolBranchSaddleIsoValue; + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + inArray, topVolumeBranch, topVolBranchSaddleIsoValue); + b->TopVolumeBranchSaddleIsoValue = topVolBranchSaddleIsoValue; + }; + + b->BranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); + }); + + // We apply all-to-all broadcast to collect the top nSavedBranches branches by volume + vtkmdiy::all_to_all( + branch_top_volume_master, + assigner, + vtkm::filter::scalar_topology::internal::SelectTopVolumeContoursFunctor(this->nSavedBranches)); + + // For each block, we compute the get the extracted isosurface for every selected branch + // storing format: key (branch ID) - Value (list of meshes in the isosurface) + + std::vector outputDataSets(input.GetNumberOfPartitions()); + + branch_top_volume_master.foreach ( + [&](SelectTopVolumeContoursBlock* b, const vtkmdiy::Master::ProxyWithLink&) { + vtkm::cont::Field TopVolBranchGRIdField("TopVolumeBranchGlobalRegularIds", + vtkm::cont::Field::Association::WholeDataSet, + b->TopVolumeBranchRootGRId); + outputDataSets[b->LocalBlockNo].AddField(TopVolBranchGRIdField); + vtkm::cont::Field TopVolBranchVolumeField("TopVolumeBranchVolume", + vtkm::cont::Field::Association::WholeDataSet, + b->TopVolumeBranchVolume); + outputDataSets[b->LocalBlockNo].AddField(TopVolBranchVolumeField); + vtkm::cont::Field TopVolBranchSaddleEpsilonField("TopVolumeBranchSaddleEpsilon", + vtkm::cont::Field::Association::WholeDataSet, + b->TopVolumeBranchSaddleEpsilon); + outputDataSets[b->LocalBlockNo].AddField(TopVolBranchSaddleEpsilonField); + + auto resolveArray = [&](const auto& inArray) { + vtkm::cont::Field TopVolBranchSaddleIsoValueField( + "TopVolumeBranchSaddleIsoValue", vtkm::cont::Field::Association::WholeDataSet, inArray); + outputDataSets[b->LocalBlockNo].AddField(TopVolBranchSaddleIsoValueField); + }; + b->TopVolumeBranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); + }); + + return vtkm::cont::PartitionedDataSet{ outputDataSets }; +} + +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h new file mode 100644 index 000000000..cdc8c2f3f --- /dev/null +++ b/vtkm/filter/scalar_topology/SelectTopVolumeContoursFilter.h @@ -0,0 +1,83 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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. +// +//============================================================================= + +#ifndef vtk_m_filter_scalar_topology_SelectTopVolumeContoursFilter_h +#define vtk_m_filter_scalar_topology_SelectTopVolumeContoursFilter_h + +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +/// \brief Compute branch decompostion from distributed contour tree + +class VTKM_FILTER_SCALAR_TOPOLOGY_EXPORT SelectTopVolumeContoursFilter : public vtkm::filter::Filter +{ +public: + VTKM_CONT SelectTopVolumeContoursFilter() = default; + + VTKM_CONT void SetSavedBranches(const vtkm::Id& numBranches) + { + this->nSavedBranches = numBranches; + } + +private: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet&) override; + VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions( + const vtkm::cont::PartitionedDataSet& inData) override; + + vtkm::Id nSavedBranches; + vtkm::worklet::contourtree_augmented::IdArrayType BranchVolume; + vtkm::worklet::contourtree_augmented::IdArrayType BranchSaddleEpsilon; + vtkm::worklet::contourtree_augmented::IdArrayType SortedBranchByVolume; + vtkm::cont::UnknownArrayHandle BranchSaddleIsoValue; +}; + +} // namespace scalar_topology +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/internal/CMakeLists.txt b/vtkm/filter/scalar_topology/internal/CMakeLists.txt index 687c753ae..3280d72b8 100644 --- a/vtkm/filter/scalar_topology/internal/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/internal/CMakeLists.txt @@ -10,8 +10,10 @@ set(headers BranchDecompositionBlock.h + SelectTopVolumeContoursBlock.h ComputeBlockIndices.h ComputeDistributedBranchDecompositionFunctor.h + SelectTopVolumeContoursFunctor.h ExchangeBranchEndsFunctor.h ) #----------------------------------------------------------------------------- diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx new file mode 100644 index 000000000..9a3e9d2ca --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.cxx @@ -0,0 +1,229 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#include +#include +#include +#include +#include + +#ifdef DEBUG_PRINT +#include +#endif + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +namespace internal +{ + +SelectTopVolumeContoursBlock::SelectTopVolumeContoursBlock(vtkm::Id localBlockNo, int globalBlockId) + : LocalBlockNo(localBlockNo) + , GlobalBlockId(globalBlockId) +{ +} + +void SelectTopVolumeContoursBlock::SortBranchByVolume( + const vtkm::cont::DataSet& hierarchicalTreeDataSet, + const vtkm::Id totalVolume) +{ + /// Pipeline to compute the branch volume + /// 1. check both ends of the branch. If both leaves, then main branch, volume = totalVolume + /// 2. for other branches, check the direction of the inner superarc + /// branch volume = (inner superarc points to the senior-most node) ? + /// dependentVolume[innerSuperarc] : + /// reverseVolume[innerSuperarc] + /// NOTE: reverseVolume = totalVolume - dependentVolume + intrinsicVolume + + // Generally, if ending superarc has intrinsicVol == dependentVol, then it is a leaf node + vtkm::cont::ArrayHandle isLowerLeaf; + vtkm::cont::ArrayHandle isUpperLeaf; + + auto upperEndIntrinsicVolume = hierarchicalTreeDataSet.GetField("UpperEndIntrinsicVolume") + .GetData() + .AsArrayHandle>(); + auto upperEndDependentVolume = hierarchicalTreeDataSet.GetField("UpperEndDependentVolume") + .GetData() + .AsArrayHandle>(); + auto lowerEndIntrinsicVolume = hierarchicalTreeDataSet.GetField("LowerEndIntrinsicVolume") + .GetData() + .AsArrayHandle>(); + auto lowerEndDependentVolume = hierarchicalTreeDataSet.GetField("LowerEndDependentVolume") + .GetData() + .AsArrayHandle>(); + auto lowerEndSuperarcId = hierarchicalTreeDataSet.GetField("LowerEndSuperarcId") + .GetData() + .AsArrayHandle>(); + auto upperEndSuperarcId = hierarchicalTreeDataSet.GetField("UpperEndSuperarcId") + .GetData() + .AsArrayHandle>(); + auto branchRoot = hierarchicalTreeDataSet.GetField("BranchRoot") + .GetData() + .AsArrayHandle>(); + + vtkm::cont::Algorithm::Transform( + upperEndIntrinsicVolume, upperEndDependentVolume, isUpperLeaf, vtkm::Equal()); + vtkm::cont::Algorithm::Transform( + lowerEndIntrinsicVolume, lowerEndDependentVolume, isLowerLeaf, vtkm::Equal()); + + + // NOTE: special cases (one-superarc branches) exist + // if the upper end superarc == lower end superarc == branch root superarc + // then it's probably not a leaf-leaf branch (Both equality has to be satisfied!) + // exception: the entire domain has only one superarc (intrinsic == dependent == total - 1) + // then it is a leaf-leaf branch + vtkm::cont::Invoker invoke; + + vtkm::worklet::scalar_topology::select_top_volume_contours::ClarifyBranchEndSupernodeTypeWorklet + clarifyNodeTypeWorklet(totalVolume); + + invoke(clarifyNodeTypeWorklet, + lowerEndSuperarcId, + lowerEndIntrinsicVolume, + upperEndSuperarcId, + upperEndIntrinsicVolume, + branchRoot, + isLowerLeaf, + isUpperLeaf); + + vtkm::cont::UnknownArrayHandle upperEndValue = + hierarchicalTreeDataSet.GetField("UpperEndValue").GetData(); + + // Based on the direction info of the branch, store epsilon direction and isovalue of the saddle + auto resolveArray = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + using ValueType = typename InArrayHandleType::ValueType; + + vtkm::cont::ArrayHandle branchSaddleIsoValue; + branchSaddleIsoValue.Allocate(isLowerLeaf.GetNumberOfValues()); + this->BranchSaddleEpsilon.Allocate(isLowerLeaf.GetNumberOfValues()); + + vtkm::worklet::scalar_topology::select_top_volume_contours::UpdateInfoByBranchDirectionWorklet< + ValueType> + updateInfoWorklet; + auto lowerEndValue = hierarchicalTreeDataSet.GetField("LowerEndValue") + .GetData() + .AsArrayHandle>(); + + invoke(updateInfoWorklet, + isLowerLeaf, + isUpperLeaf, + inArray, + lowerEndValue, + this->BranchSaddleEpsilon, + branchSaddleIsoValue); + this->BranchSaddleIsoValue = branchSaddleIsoValue; + }; + + upperEndValue.CastAndCallForTypes( + resolveArray); + + vtkm::worklet::contourtree_augmented::IdArrayType branchVolume; + vtkm::worklet::scalar_topology::select_top_volume_contours::GetBranchVolumeWorklet + getBranchVolumeWorklet(totalVolume); + + invoke(getBranchVolumeWorklet, // worklet + lowerEndSuperarcId, // input + lowerEndIntrinsicVolume, // input + lowerEndDependentVolume, // input + upperEndSuperarcId, // input + upperEndIntrinsicVolume, // input + upperEndDependentVolume, // input + isLowerLeaf, + isUpperLeaf, + branchVolume); // output + +#ifdef DEBUG_PRINT + std::stringstream resultStream; + resultStream << "Branch Volume In The Block" << std::endl; + const vtkm::Id nVolume = branchVolume.GetNumberOfValues(); + vtkm::worklet::contourtree_augmented::PrintHeader(nVolume, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "BranchVolume", branchVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices("isLowerLeaf", isLowerLeaf, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices("isUpperLeaf", isUpperLeaf, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "LowerEndIntrinsicVol", lowerEndIntrinsicVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "LowerEndDependentVol", lowerEndDependentVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "UpperEndIntrinsicVol", upperEndIntrinsicVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "UpperEndDependentVol", upperEndDependentVolume, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "LowerEndSuperarc", lowerEndSuperarcId, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices( + "UpperEndSuperarc", upperEndSuperarcId, -1, resultStream); + vtkm::worklet::contourtree_augmented::PrintIndices("BranchRoot", branchRoot, -1, resultStream); + resultStream << std::endl; + VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); +#endif + + vtkm::cont::Algorithm::Copy(branchVolume, this->BranchVolume); + + const vtkm::Id nBranches = lowerEndSuperarcId.GetNumberOfValues(); + vtkm::cont::ArrayHandleIndex branchesIdx(nBranches); + vtkm::worklet::contourtree_augmented::IdArrayType sortedBranches; + vtkm::cont::Algorithm::Copy(branchesIdx, sortedBranches); + + // sort the branch volume + vtkm::cont::Algorithm::SortByKey(branchVolume, sortedBranches, vtkm::SortGreater()); + vtkm::cont::Algorithm::Copy(sortedBranches, this->SortedBranchByVolume); +} + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h new file mode 100644 index 000000000..7d08433ed --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursBlock.h @@ -0,0 +1,98 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_internal_SelectTopVolumeContoursBlock_h +#define vtk_m_filter_scalar_topology_internal_SelectTopVolumeContoursBlock_h + +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +namespace internal +{ + +struct SelectTopVolumeContoursBlock +{ + SelectTopVolumeContoursBlock(vtkm::Id localBlockNo, int globalBlockId); + + void SortBranchByVolume(const vtkm::cont::DataSet& hierarchicalTreeDataSet, + const vtkm::Id totalVolume); + + // Block metadata + vtkm::Id LocalBlockNo; + int GlobalBlockId; // TODO/FIXME: Check whether really needed. Possibly only during debugging + + vtkm::worklet::contourtree_augmented::IdArrayType BranchVolume; + vtkm::worklet::contourtree_augmented::IdArrayType BranchSaddleEpsilon; + vtkm::worklet::contourtree_augmented::IdArrayType SortedBranchByVolume; + vtkm::cont::UnknownArrayHandle BranchSaddleIsoValue; + + // Output Datasets. + vtkm::worklet::contourtree_augmented::IdArrayType TopVolumeBranchRootGRId; + vtkm::worklet::contourtree_augmented::IdArrayType TopVolumeBranchVolume; + vtkm::cont::UnknownArrayHandle TopVolumeBranchSaddleIsoValue; + vtkm::worklet::contourtree_augmented::IdArrayType TopVolumeBranchSaddleEpsilon; + + // Destroy function allowing DIY to own blocks and clean them up after use + static void Destroy(void* b) { delete static_cast(b); } +}; + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm +#endif diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.cxx b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.cxx new file mode 100644 index 000000000..485c14d2e --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.cxx @@ -0,0 +1,384 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#include +#include +#include +#include + +#include + +#ifdef DEBUG_PRINT +#define DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH +#include +#endif + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +namespace internal +{ + +void SelectTopVolumeContoursFunctor::operator()( + SelectTopVolumeContoursBlock* b, + const vtkmdiy::ReduceProxy& rp // communication proxy +) const +{ + if (this->nSavedBranches < 1) + return; + // Get our rank and DIY id + //const vtkm::Id rank = vtkm::cont::EnvironmentTracker::GetCommunicator().rank(); + const auto selfid = rp.gid(); + + // Aliases to reduce verbosity + using IdArrayType = vtkm::worklet::contourtree_augmented::IdArrayType; + + vtkm::cont::Invoker invoke; + + if (rp.in_link().size() == 0) + { + for (int cc = 0; cc < rp.out_link().size(); ++cc) + { + auto target = rp.out_link().target(cc); + if (target.gid != selfid) + { +#ifdef DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH + rp.enqueue(target, b->GlobalBlockId); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + "Block " << b->GlobalBlockId << " enqueue to Block " << target.gid); +#endif + auto topVolBranchRootGRIdPortal = b->TopVolumeBranchRootGRId.ReadPortal(); + auto topVolBranchVolumePortal = b->TopVolumeBranchVolume.ReadPortal(); + auto topVolBranchSaddleEpsilonPortal = b->TopVolumeBranchSaddleEpsilon.ReadPortal(); + + vtkm::Id nBranches = topVolBranchRootGRIdPortal.GetNumberOfValues(); + + rp.enqueue(target, nBranches); + for (vtkm::Id branch = 0; branch < nBranches; ++branch) + rp.enqueue(target, topVolBranchRootGRIdPortal.Get(branch)); + for (vtkm::Id branch = 0; branch < nBranches; ++branch) + rp.enqueue(target, topVolBranchVolumePortal.Get(branch)); + for (vtkm::Id branch = 0; branch < nBranches; ++branch) + rp.enqueue(target, topVolBranchSaddleEpsilonPortal.Get(branch)); + + auto resolveArray = [&](const auto& inArray) { + using InArrayHandleType = std::decay_t; + using ValueType = typename InArrayHandleType::ValueType; + auto topVolBranchSaddleIsoValuePortal = inArray.ReadPortal(); + for (vtkm::Id branch = 0; branch < nBranches; ++branch) + rp.enqueue(target, topVolBranchSaddleIsoValuePortal.Get(branch)); + }; + b->TopVolumeBranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); + + // rp.enqueue(target, b->TopVolumeBranchRootGRId); + // rp.enqueue(target, b->TopVolumeBranchVolume); + } + } + } + else + { + for (int i = 0; i < rp.in_link().size(); ++i) + { + int ingid = rp.in_link().target(i).gid; + if (ingid == selfid) + continue; + + // copy incoming to the block +#ifdef DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH + int incomingGlobalBlockId; + rp.dequeue(ingid, incomingGlobalBlockId); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + "Combining local block " << b->GlobalBlockId << " with incoming block " + << incomingGlobalBlockId); +#endif + + // dequeue the data from other blocks. + // nIncomingBranches + // array of incoming branch global regular ID + // array of incoming branch volume + // array of branch epsilon direction + // array of branch saddle end value + // TODO/FIXME: This is a workaround for a bug in DIY/vtk-m. + // Replace with dequeuing ArrayHandles once bug is fixed. + vtkm::Id nIncoming; + rp.dequeue(ingid, nIncoming); + + IdArrayType incomingTopVolBranchGRId; + incomingTopVolBranchGRId.Allocate(nIncoming); + auto incomingTopVolBranchGRIdPortal = incomingTopVolBranchGRId.WritePortal(); + for (vtkm::Id branch = 0; branch < nIncoming; ++branch) + { + vtkm::Id incomingTmpBranchGRId; + rp.dequeue(ingid, incomingTmpBranchGRId); + incomingTopVolBranchGRIdPortal.Set(branch, incomingTmpBranchGRId); + } + + // TODO/FIXME: This is a workaround for a bug in DIY/vtk-m. + // Replace with dequeuing ArrayHandles once bug is fixed. + // rp.dequeue(ingid, incomingTopVolBranchGRId); + + IdArrayType incomingTopVolBranchVolume; + incomingTopVolBranchVolume.Allocate(nIncoming); + auto incomingTopVolBranchVolumePortal = incomingTopVolBranchVolume.WritePortal(); + for (vtkm::Id branch = 0; branch < nIncoming; ++branch) + { + vtkm::Id incomingTmpBranchVolume; + rp.dequeue(ingid, incomingTmpBranchVolume); + incomingTopVolBranchVolumePortal.Set(branch, incomingTmpBranchVolume); + } + + // TODO/FIXME: This is a workaround for a bug in DIY/vtk-m. + // Replace with dequeuing ArrayHandles once bug is fixed. + // rp.dequeue(ingid, incomingTopVolBranchVolume); + + IdArrayType incomingTopVolBranchSaddleEpsilon; + incomingTopVolBranchSaddleEpsilon.Allocate(nIncoming); + auto incomingTopVolBranchSaddleEpsilonPortal = + incomingTopVolBranchSaddleEpsilon.WritePortal(); + for (vtkm::Id branch = 0; branch < nIncoming; ++branch) + { + vtkm::Id incomingTmpBranchSaddleEpsilon; + rp.dequeue(ingid, incomingTmpBranchSaddleEpsilon); + incomingTopVolBranchSaddleEpsilonPortal.Set(branch, incomingTmpBranchSaddleEpsilon); + } + + // TODO/FIXME: This is a workaround for a bug in DIY/vtk-m. + // Replace with dequeuing ArrayHandles once bug is fixed. + // rp.dequeue(ingid, incomingTopVolBranchSaddleEpsilon); + + auto resolveArray = [&](auto& inArray) { + using InArrayHandleType = std::decay_t; + using ValueType = typename InArrayHandleType::ValueType; + InArrayHandleType incomingTopVolBranchSaddleIsoValue; + incomingTopVolBranchSaddleIsoValue.Allocate(nIncoming); + auto incomingTopVolBranchSaddleIsoValuePortal = + incomingTopVolBranchSaddleIsoValue.WritePortal(); + for (vtkm::Id branch = 0; branch < nIncoming; ++branch) + { + ValueType incomingSaddleValue; + rp.dequeue(ingid, incomingSaddleValue); + incomingTopVolBranchSaddleIsoValuePortal.Set(branch, incomingSaddleValue); + } + + // TODO/FIXME: This is a workaround for a bug in DIY/vtk-m. + // Replace with dequeuing ArrayHandles once bug is fixed. + // rp.dequeue(ingid, incomingTopVolBranchSaddleIsoValue); + + vtkm::Id nSelf = b->TopVolumeBranchRootGRId.GetNumberOfValues(); + +#ifdef DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH + VTKM_LOG_S(vtkm::cont::LogLevel::Info, + "nIncoming = " << nIncoming << ", nSelf = " << nSelf); + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader(nIncoming, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingTopBranchId", incomingTopVolBranchGRId, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "incomingTopBranchVol", incomingTopVolBranchVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintValues( + "incomingSaddleVal", incomingTopVolBranchSaddleIsoValue, -1, rs); + + vtkm::worklet::contourtree_augmented::PrintHeader(nSelf, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "selfTopBranchId", b->TopVolumeBranchRootGRId, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "selfTopBranchVol", b->TopVolumeBranchVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintValues( + "selfTopSaddleVal", inArray, -1, rs); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, rs.str()); + } +#endif + // merge incoming branches with self branches + IdArrayType mergedTopVolBranchGRId; + IdArrayType mergedTopVolBranchVolume; + IdArrayType mergedTopVolBranchSaddleEpsilon; + InArrayHandleType mergedTopVolBranchSaddleIsoValue; + mergedTopVolBranchGRId.Allocate(nIncoming + nSelf); + mergedTopVolBranchVolume.Allocate(nIncoming + nSelf); + mergedTopVolBranchSaddleEpsilon.Allocate(nIncoming + nSelf); + mergedTopVolBranchSaddleIsoValue.Allocate(nIncoming + nSelf); + + vtkm::cont::Algorithm::CopySubRange( + incomingTopVolBranchGRId, 0, nIncoming, mergedTopVolBranchGRId, 0); + vtkm::cont::Algorithm::CopySubRange( + incomingTopVolBranchVolume, 0, nIncoming, mergedTopVolBranchVolume, 0); + vtkm::cont::Algorithm::CopySubRange( + incomingTopVolBranchSaddleEpsilon, 0, nIncoming, mergedTopVolBranchSaddleEpsilon, 0); + vtkm::cont::Algorithm::CopySubRange( + incomingTopVolBranchSaddleIsoValue, 0, nIncoming, mergedTopVolBranchSaddleIsoValue, 0); + vtkm::cont::Algorithm::CopySubRange( + b->TopVolumeBranchRootGRId, 0, nSelf, mergedTopVolBranchGRId, nIncoming); + vtkm::cont::Algorithm::CopySubRange( + b->TopVolumeBranchVolume, 0, nSelf, mergedTopVolBranchVolume, nIncoming); + vtkm::cont::Algorithm::CopySubRange( + b->TopVolumeBranchSaddleEpsilon, 0, nSelf, mergedTopVolBranchSaddleEpsilon, nIncoming); + vtkm::cont::Algorithm::CopySubRange( + inArray, 0, nSelf, mergedTopVolBranchSaddleIsoValue, nIncoming); + + // Sort all branches (incoming + self) based on volume + // sorting key: (volume, branch global regular ID) + // the highest volume comes first, the lowest branch GR ID comes first + vtkm::cont::ArrayHandleIndex mergedBranchId(nIncoming + nSelf); + IdArrayType sortedBranchId; + vtkm::cont::Algorithm::Copy(mergedBranchId, sortedBranchId); + vtkm::worklet::scalar_topology::select_top_volume_contours::BranchVolumeComparator + branchVolumeComparator(mergedTopVolBranchGRId, mergedTopVolBranchVolume); + vtkm::cont::Algorithm::Sort(sortedBranchId, branchVolumeComparator); + + // permute the branch information based on sorting + IdArrayType permutedTopVolBranchGRId; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + mergedTopVolBranchGRId, sortedBranchId, permutedTopVolBranchGRId); + IdArrayType permutedTopVolBranchVolume; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + mergedTopVolBranchVolume, sortedBranchId, permutedTopVolBranchVolume); + IdArrayType permutedTopVolBranchSaddleEpsilon; + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( + mergedTopVolBranchSaddleEpsilon, sortedBranchId, permutedTopVolBranchSaddleEpsilon); + InArrayHandleType permutedTopVolBranchSaddleIsoValue; + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( + mergedTopVolBranchSaddleIsoValue, sortedBranchId, permutedTopVolBranchSaddleIsoValue); + +#ifdef DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader(nIncoming + nSelf, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "permutedTopBranchId", permutedTopVolBranchGRId, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "permutedTopBranchVol", permutedTopVolBranchVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintValues( + "permutedTopSaddleVal", permutedTopVolBranchSaddleIsoValue, -1, rs); + } +#endif + + // there may be duplicate branches. We remove duplicate branches based on global regular IDs + // We can reuse the filter from removing duplicate branches in the process of collecting branches + IdArrayType oneIfUniqueBranch; + oneIfUniqueBranch.Allocate(nIncoming + nSelf); + vtkm::worklet::scalar_topology::hierarchical_volumetric_branch_decomposer:: + OneIfBranchEndWorklet oneIfUniqueWorklet; + invoke(oneIfUniqueWorklet, mergedBranchId, permutedTopVolBranchGRId, oneIfUniqueBranch); + + // Remove duplicate + IdArrayType mergedUniqueBranchGRId; + IdArrayType mergedUniqueBranchVolume; + IdArrayType mergedUniqueBranchSaddleEpsilon; + InArrayHandleType mergedUniqueBranchSaddleIsoValue; + + vtkm::cont::Algorithm::CopyIf( + permutedTopVolBranchGRId, oneIfUniqueBranch, mergedUniqueBranchGRId); + vtkm::cont::Algorithm::CopyIf( + permutedTopVolBranchVolume, oneIfUniqueBranch, mergedUniqueBranchVolume); + vtkm::cont::Algorithm::CopyIf( + permutedTopVolBranchSaddleEpsilon, oneIfUniqueBranch, mergedUniqueBranchSaddleEpsilon); + vtkm::cont::Algorithm::CopyIf( + permutedTopVolBranchSaddleIsoValue, oneIfUniqueBranch, mergedUniqueBranchSaddleIsoValue); + + vtkm::Id nMergedUnique = mergedUniqueBranchGRId.GetNumberOfValues(); + +#ifdef DEBUG_PRINT_COMBINED_HIGH_VOLUME_BRANCH + { + std::stringstream rs; + vtkm::worklet::contourtree_augmented::PrintHeader(nMergedUnique, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "mergedUniqueBranchId", mergedUniqueBranchGRId, -1, rs); + vtkm::worklet::contourtree_augmented::PrintIndices( + "mergedUniqueBranchVol", mergedUniqueBranchVolume, -1, rs); + vtkm::worklet::contourtree_augmented::PrintValues( + "mergedUniqueSaddleVal", mergedUniqueBranchSaddleIsoValue, -1, rs); + } +#endif + + // After removing duplicate, if there are more branches than we need + // We only save the top nSavedBranches branches + if (nMergedUnique > this->nSavedBranches) + { + vtkm::cont::Algorithm::CopySubRange( + mergedUniqueBranchGRId, 0, this->nSavedBranches, b->TopVolumeBranchRootGRId); + vtkm::cont::Algorithm::CopySubRange( + mergedUniqueBranchVolume, 0, this->nSavedBranches, b->TopVolumeBranchVolume); + vtkm::cont::Algorithm::CopySubRange(mergedUniqueBranchSaddleEpsilon, + 0, + this->nSavedBranches, + b->TopVolumeBranchSaddleEpsilon); + // InArrayHandleType subRangeUniqueBranchSaddleIsoValue; + inArray.Allocate(this->nSavedBranches); + vtkm::cont::Algorithm::CopySubRange( + mergedUniqueBranchSaddleIsoValue, 0, this->nSavedBranches, inArray); + // inArray = subRangeUniqueBranchSaddleIsoValue; + } + else + { + vtkm::cont::Algorithm::Copy(mergedUniqueBranchGRId, b->TopVolumeBranchRootGRId); + vtkm::cont::Algorithm::Copy(mergedUniqueBranchVolume, b->TopVolumeBranchVolume); + vtkm::cont::Algorithm::Copy(mergedUniqueBranchSaddleEpsilon, + b->TopVolumeBranchSaddleEpsilon); + inArray.Allocate(nMergedUnique); + vtkm::cont::Algorithm::Copy(mergedUniqueBranchSaddleIsoValue, inArray); + } + }; + b->TopVolumeBranchSaddleIsoValue + .CastAndCallForTypes(resolveArray); + } + } +} + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.h b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.h new file mode 100644 index 000000000..023047def --- /dev/null +++ b/vtkm/filter/scalar_topology/internal/SelectTopVolumeContoursFunctor.h @@ -0,0 +1,93 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_internal_SelectTopVolumeContoursFunctor_h +#define vtk_m_filter_scalar_topology_internal_SelectTopVolumeContoursFunctor_h + +#include + +// clang-format off +VTKM_THIRDPARTY_PRE_INCLUDE +#include +VTKM_THIRDPARTY_POST_INCLUDE +// clang-format on + + +namespace vtkm +{ +namespace filter +{ +namespace scalar_topology +{ +namespace internal +{ + +struct SelectTopVolumeContoursFunctor +{ + SelectTopVolumeContoursFunctor(const vtkm::Id& nSB) + : nSavedBranches(nSB) + { + } + + void operator()(SelectTopVolumeContoursBlock* b, + const vtkmdiy::ReduceProxy& rp // communication proxy + ) const; + + const vtkm::Id nSavedBranches; +}; + +} // namespace internal +} // namespace scalar_topology +} // namespace filter +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h b/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h index 98419f64e..909cf7076 100644 --- a/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h +++ b/vtkm/filter/scalar_topology/testing/TestingContourTreeUniformDistributedFilter.h @@ -63,6 +63,7 @@ #include #include #include +#include #include #include #include @@ -466,6 +467,14 @@ inline void TestContourTreeUniformDistributedBranchDecomposition8x9(int nBlocks, augmentHierarchicalTree, computeHierarchicalVolumetricBranchDecomposition); + using vtkm::filter::scalar_topology::SelectTopVolumeContoursFilter; + + vtkm::Id numBranches = 2; + SelectTopVolumeContoursFilter tp_filter; + + tp_filter.SetSavedBranches(numBranches); + + auto tp_result = tp_filter.Execute(result); if (vtkm::cont::EnvironmentTracker::GetCommunicator().rank() == 0) { @@ -503,27 +512,99 @@ inline void TestContourTreeUniformDistributedBranchDecomposition8x9(int nBlocks, std::sort(computed.begin(), computed.end()); std::sort(expected.begin(), expected.end()); - if (computed == expected) + if (computed != expected) { - std::cout << "Branch Decomposition: Results Match!" << std::endl; - return; + std::cout << "Branch Decomposition Results:" << std::endl; + std::cout << "Computed Contour Tree" << std::endl; + for (std::size_t i = 0; i < computed.size(); i++) + { + std::cout << std::setw(12) << computed[i].low << std::setw(14) << computed[i].high + << std::endl; + } + + std::cout << "Expected Contour Tree" << std::endl; + for (std::size_t i = 0; i < expected.size(); i++) + { + std::cout << std::setw(12) << expected[i].low << std::setw(14) << expected[i].high + << std::endl; + } + VTKM_TEST_FAIL("Branch Decomposition Failed!"); } - std::cout << "Branch Decomposition Results:" << std::endl; - std::cout << "Computed Contour Tree" << std::endl; - for (std::size_t i = 0; i < computed.size(); i++) + std::cout << "Branch Decomposition: Results Match!" << std::endl; + + for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no) { - std::cout << std::setw(12) << computed[i].low << std::setw(14) << computed[i].high - << std::endl; + auto ds = tp_result.GetPartition(ds_no); + auto topVolBranchGRId = ds.GetField("TopVolumeBranchGlobalRegularIds") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchVolume = ds.GetField("TopVolumeBranchVolume") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchSaddleEpsilon = ds.GetField("TopVolumeBranchSaddleEpsilon") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + auto topVolBranchSaddleIsoValue = ds.GetField("TopVolumeBranchSaddleIsoValue") + .GetData() + .AsArrayHandle>() + .ReadPortal(); + + vtkm::Id nSelectedBranches = topVolBranchGRId.GetNumberOfValues(); + Edge expectedGRIdVolumeAtBranch0(38, 6); + Edge expectedEpsilonIsoAtBranch0(1, 50); + Edge expectedGRIdVolumeAtBranch1(50, 2); + Edge expectedEpsilonIsoAtBranch1(-1, 30); + + for (vtkm::Id branch = 0; branch < nSelectedBranches; ++branch) + { + bool failed = false; + Edge computedGRIdVolume(topVolBranchGRId.Get(branch), topVolBranchVolume.Get(branch)); + Edge computedEpsilonIso(topVolBranchSaddleEpsilon.Get(branch), + (vtkm::Id)topVolBranchSaddleIsoValue.Get(branch)); + + switch (branch) + { + case 0: + failed = !(computedGRIdVolume == expectedGRIdVolumeAtBranch0); + failed = (failed || !(computedEpsilonIso == expectedEpsilonIsoAtBranch0)); + break; + case 1: + failed = !(computedGRIdVolume == expectedGRIdVolumeAtBranch1); + failed = (failed || !(computedEpsilonIso == expectedEpsilonIsoAtBranch1)); + break; + default: + VTKM_TEST_ASSERT(false); + } + + if (failed) + { + std::vector expectedGRIdVolume{ expectedGRIdVolumeAtBranch0, + expectedGRIdVolumeAtBranch1 }; + + std::vector expectedEpsilonIso{ expectedEpsilonIsoAtBranch0, + expectedEpsilonIsoAtBranch1 }; + + std::cout << "Top Branch Volume Results:" << std::endl; + std::cout << "Computed Top Branch Volume:branch=" << branch << std::endl; + std::cout << computedGRIdVolume.low << std::setw(14) << computedGRIdVolume.high + << std::setw(5) << computedEpsilonIso.low << std::setw(14) + << computedEpsilonIso.high << std::endl; + + std::cout << "Expected Top Branch Volume:branch=" << branch << std::endl; + std::cout << expectedGRIdVolume[branch].low << std::setw(14) + << expectedGRIdVolume[branch].high << std::setw(5) + << expectedEpsilonIso[branch].low << std::setw(14) + << expectedEpsilonIso[branch].high << std::endl; + VTKM_TEST_FAIL("Top Branch Volume Computation Failed!"); + } + } } - std::cout << "Expected Contour Tree" << std::endl; - for (std::size_t i = 0; i < expected.size(); i++) - { - std::cout << std::setw(12) << expected[i].low << std::setw(14) << expected[i].high - << std::endl; - } - VTKM_TEST_FAIL("Branch Decomposition Failed!"); + std::cout << "Top Branch Volume: Results Match!" << std::endl; } } diff --git a/vtkm/filter/scalar_topology/worklet/CMakeLists.txt b/vtkm/filter/scalar_topology/worklet/CMakeLists.txt index ffcf92fa5..e0a604670 100644 --- a/vtkm/filter/scalar_topology/worklet/CMakeLists.txt +++ b/vtkm/filter/scalar_topology/worklet/CMakeLists.txt @@ -16,6 +16,7 @@ set(headers vtkm_declare_headers(${headers}) add_subdirectory(branch_decomposition) +add_subdirectory(select_top_volume_contours) add_subdirectory(contourtree) add_subdirectory(contourtree_augmented) add_subdirectory(contourtree_distributed) diff --git a/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h b/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h index bea7e8772..35a8b8609 100644 --- a/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h +++ b/vtkm/filter/scalar_topology/worklet/branch_decomposition/HierarchicalVolumetricBranchDecomposer.h @@ -600,7 +600,7 @@ inline void HierarchicalVolumetricBranchDecomposer::CollectEndsOfBranches( { IdArrayType superarcGRId; vtkm::cont::make_ArrayHandlePermutation(supernodes, globalRegularIds); - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( globalRegularIds, supernodes, superarcGRId); std::stringstream resultStream; @@ -615,7 +615,7 @@ inline void HierarchicalVolumetricBranchDecomposer::CollectEndsOfBranches( using InArrayHandleType = std::decay_t; using ValueType = typename InArrayHandleType::ValueType; vtkm::cont::ArrayHandle superarcValue; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithRawIndex( inArray, supernodes, superarcValue); vtkm::worklet::contourtree_augmented::PrintValues( "Data Values", superarcValue, -1, resultStream); @@ -627,7 +627,6 @@ inline void HierarchicalVolumetricBranchDecomposer::CollectEndsOfBranches( "Intrinsic Volumes", intrinsicVolumes, -1, resultStream); vtkm::worklet::contourtree_augmented::PrintIndices( "Dependent Volumes", dependentVolumes, -1, resultStream); - resultStream << std::endl; VTKM_LOG_S(vtkm::cont::LogLevel::Info, resultStream.str()); } @@ -690,27 +689,27 @@ inline void HierarchicalVolumetricBranchDecomposer::CollectEndsOfBranches( // actualOuterNodeValues[i] = dataValues[supernodes[outerNodes[actualSuperarcs[i]]]]; // actualOuterNodeRegularIds[i] = globalRegularIds[supernodes[outerNodes[actualSuperarcs[i]]]]; // } - // Solution: PermuteArray helps allocate the space so no need for explicit allocation - // PermuteArray also calls MaskedIndex + // Solution: PermuteArrayWithMaskedIndex helps allocate the space so no need for explicit allocation + // PermuteArrayWithMaskedIndex also calls MaskedIndex // IdArrayType, size: nActualSuperarcs IdArrayType actualBranchRoots; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( branchRoots, actualSuperarcs, actualBranchRoots); // IdArrayType, size: nActualSuperarcs IdArrayType actualOuterNodes; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( outerNodes, actualSuperarcs, actualOuterNodes); // IdArrayType, size: nActualSuperarcs IdArrayType actualOuterNodeLocalIds; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( supernodes, actualOuterNodes, actualOuterNodeLocalIds); // IdArrayType, size: nActualSuperarcs IdArrayType actualOuterNodeRegularIds; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( globalRegularIds, actualOuterNodeLocalIds, actualOuterNodeRegularIds); auto resolveArray = [&](const auto& inArray) { @@ -744,7 +743,7 @@ inline void HierarchicalVolumetricBranchDecomposer::CollectEndsOfBranches( // This is the real superarc local ID after permutation IdArrayType permutedActualSuperarcs; - vtkm::worklet::contourtree_augmented::PermuteArray( + vtkm::worklet::contourtree_augmented::PermuteArrayWithMaskedIndex( actualSuperarcs, sortedSuperarcs, permutedActualSuperarcs); #ifdef DEBUG_HIERARCHICAL_VOLUMETRIC_BRANCH_DECOMPOSER diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ActiveGraph.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ActiveGraph.h index c38efcbe7..e6b8d1ba8 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ActiveGraph.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ActiveGraph.h @@ -942,13 +942,13 @@ inline void ActiveGraph::DebugPrint(const char* message, const char* fileName, l // Active Vertex Arrays IdArrayType activeIndices; - PermuteArray(this->GlobalIndex, this->ActiveVertices, activeIndices); + PermuteArrayWithMaskedIndex(this->GlobalIndex, this->ActiveVertices, activeIndices); IdArrayType activeFirst; - PermuteArray(this->FirstEdge, this->ActiveVertices, activeFirst); + PermuteArrayWithMaskedIndex(this->FirstEdge, this->ActiveVertices, activeFirst); IdArrayType activeOutdegree; - PermuteArray(this->Outdegree, this->ActiveVertices, activeOutdegree); + PermuteArrayWithMaskedIndex(this->Outdegree, this->ActiveVertices, activeOutdegree); IdArrayType activeHyperarcs; - PermuteArray(this->Hyperarcs, this->ActiveVertices, activeHyperarcs); + PermuteArrayWithMaskedIndex(this->Hyperarcs, this->ActiveVertices, activeHyperarcs); std::cout << "Active Vertex Arrays - Size: " << this->ActiveVertices.GetNumberOfValues() << std::endl; PrintHeader(this->ActiveVertices.GetNumberOfValues()); @@ -961,9 +961,9 @@ inline void ActiveGraph::DebugPrint(const char* message, const char* fileName, l // Full Edge Arrays IdArrayType farIndices; - PermuteArray(this->GlobalIndex, this->EdgeFar, farIndices); + PermuteArrayWithMaskedIndex(this->GlobalIndex, this->EdgeFar, farIndices); IdArrayType nearIndices; - PermuteArray(this->GlobalIndex, this->EdgeNear, nearIndices); + PermuteArrayWithMaskedIndex(this->GlobalIndex, this->EdgeNear, nearIndices); std::cout << "Full Edge Arrays - Size: " << this->EdgeNear.GetNumberOfValues() << std::endl; PrintHeader(this->EdgeFar.GetNumberOfValues()); PrintIndices("Near", this->EdgeNear); @@ -974,9 +974,9 @@ inline void ActiveGraph::DebugPrint(const char* message, const char* fileName, l // Active Edge Arrays IdArrayType activeFarIndices; - PermuteArray(this->EdgeFar, this->ActiveEdges, activeFarIndices); + PermuteArrayWithMaskedIndex(this->EdgeFar, this->ActiveEdges, activeFarIndices); IdArrayType activeNearIndices; - PermuteArray(this->EdgeNear, this->ActiveEdges, activeNearIndices); + PermuteArrayWithMaskedIndex(this->EdgeNear, this->ActiveEdges, activeNearIndices); std::cout << "Active Edge Arrays - Size: " << this->ActiveEdges.GetNumberOfValues() << std::endl; PrintHeader(this->ActiveEdges.GetNumberOfValues()); @@ -987,9 +987,9 @@ inline void ActiveGraph::DebugPrint(const char* message, const char* fileName, l // Edge Sorter Array IdArrayType sortedFarIndices; - PermuteArray(this->EdgeFar, this->EdgeSorter, sortedFarIndices); + PermuteArrayWithMaskedIndex(this->EdgeFar, this->EdgeSorter, sortedFarIndices); IdArrayType sortedNearIndices; - PermuteArray(this->EdgeNear, this->EdgeSorter, sortedNearIndices); + PermuteArrayWithMaskedIndex(this->EdgeNear, this->EdgeSorter, sortedNearIndices); std::cout << "Edge Sorter - Size: " << this->EdgeSorter.GetNumberOfValues() << std::endl; PrintHeader(this->EdgeSorter.GetNumberOfValues()); PrintIndices("Edge Sorter", this->EdgeSorter); diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ArrayTransforms.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ArrayTransforms.h index 95874384a..7e93db9fd 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ArrayTransforms.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ArrayTransforms.h @@ -75,7 +75,9 @@ namespace contourtree_augmented // permute routines template -inline void PermuteArray(const ArrayType& input, IdArrayType& permute, ArrayType& output) +inline void PermuteArrayWithMaskedIndex(const ArrayType& input, + IdArrayType& permute, + ArrayType& output) { // permuteValues() using transform_type = vtkm::cont::ArrayHandleTransform>; @@ -110,6 +112,38 @@ inline void PermuteArray(const ArrayType& input, IdArrayType& permute, ArrayType } // permuteValues() +// permute value type arrays +template +inline void PermuteArrayWithRawIndex(const ArrayType& input, + IdArrayType& permute, + ArrayType& output) +{ // PermuteArrayWithRawIndex() + using permute_type = vtkm::cont::ArrayHandlePermutation; + + // resize the output, i.e., output.resize(permute.size()); + vtkm::Id permNumValues = permute.GetNumberOfValues(); + vtkm::Id outNumValues = output.GetNumberOfValues(); + if (permNumValues > outNumValues) + { + output.Allocate(permNumValues); + } + else if (permNumValues < outNumValues) + { + output.Allocate(permNumValues, vtkm::CopyFlag::On); + } // else the output has already the correct size + + // The following is equivalent to doing the following in serial + // + // for (vtkm::Id entry = 0; entry < permute.size(); entry++) + // output[entry] = input[permute[entry]]; + // + // fancy vtkm array so that we do not actually copy any data here + permute_type permutedInput(permute, input); + // Finally, copy the permuted values to the output array + vtkm::cont::ArrayCopyDevice(permutedInput, output); +} // PermuteArrayWithRawIndex() + + // transform functor used in ContourTreeMesh to flag indicies as other when using the CombinedVectorClass struct MarkOther { diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ContourTreeMaker.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ContourTreeMaker.h index 039a56b4d..dfcfa1f82 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ContourTreeMaker.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ContourTreeMaker.h @@ -268,13 +268,13 @@ inline void ContourTreeMaker::ComputeHyperAndSuperStructure() // we have to permute a bunch of arrays, so let's have some temporaries to store them IdArrayType permutedHyperparents; - PermuteArray( + PermuteArrayWithMaskedIndex( this->ContourTreeResult.Hyperparents, this->ContourTreeResult.Hypernodes, permutedHyperparents); IdArrayType permutedSupernodes; - PermuteArray( + PermuteArrayWithMaskedIndex( this->ContourTreeResult.Supernodes, this->ContourTreeResult.Hypernodes, permutedSupernodes); IdArrayType permutedSuperarcs; - PermuteArray( + PermuteArrayWithMaskedIndex( this->ContourTreeResult.Superarcs, this->ContourTreeResult.Hypernodes, permutedSuperarcs); // now we establish the reverse index array @@ -284,7 +284,7 @@ inline void ContourTreeMaker::ComputeHyperAndSuperStructure() // for (vtkm::Id supernode = 0; supernode < this->ContourTreeResult.Supernodes.size(); supernode++) // superSortIndex[this->ContourTreeResult.Hypernodes[supernode]] = supernode; - //typedef vtkm::cont::ArrayHandlePermutation PermuteArrayHandleIndex; + //typedef vtkm::cont::ArrayHandlePermutation PermuteArrayWithMaskedIndexHandleIndex; vtkm::cont::ArrayHandlePermutation permutedSuperSortIndex( this->ContourTreeResult.Hypernodes, // index array superSortIndex); // value array @@ -307,7 +307,7 @@ inline void ContourTreeMaker::ComputeHyperAndSuperStructure() // we will permute the hyperarcs & copy them back with the new supernode target IDs IdArrayType permutedHyperarcs; - PermuteArray( + PermuteArrayWithMaskedIndex( this->ContourTreeResult.Hyperarcs, this->ContourTreeResult.Hypernodes, permutedHyperarcs); contourtree_maker_inc_ns::ComputeHyperAndSuperStructure_PermuteArcs permuteHyperarcsWorklet; this->Invoke( @@ -315,9 +315,9 @@ inline void ContourTreeMaker::ComputeHyperAndSuperStructure() // now swizzle the WhenTransferred value IdArrayType permutedWhenTransferred; - PermuteArray(this->ContourTreeResult.WhenTransferred, - this->ContourTreeResult.Hypernodes, - permutedWhenTransferred); + PermuteArrayWithMaskedIndex(this->ContourTreeResult.WhenTransferred, + this->ContourTreeResult.Hypernodes, + permutedWhenTransferred); vtkm::cont::Algorithm::Copy(permutedWhenTransferred, this->ContourTreeResult.WhenTransferred); // now we compress both the hypernodes & Hyperarcs diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h index 9d1aaf7d4..289bc686b 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/PrintVectors.h @@ -72,7 +72,7 @@ namespace contourtree_augmented { // local constants to allow changing the spacing as needed -constexpr int PRINT_WIDTH = 12; +constexpr int PRINT_WIDTH = 14; constexpr int PREFIX_WIDTH = 30; template diff --git a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ProcessContourTree.h b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ProcessContourTree.h index cc36eed87..f620781e0 100644 --- a/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ProcessContourTree.h +++ b/vtkm/filter/scalar_topology/worklet/contourtree_augmented/ProcessContourTree.h @@ -659,11 +659,12 @@ public: IdArrayType permutedBranches; permutedBranches.Allocate(nSupernodes); - PermuteArray(whichBranch, supernodeSorter, permutedBranches); + PermuteArrayWithMaskedIndex(whichBranch, supernodeSorter, permutedBranches); IdArrayType permutedRegularID; permutedRegularID.Allocate(nSupernodes); - PermuteArray(contourTree.Supernodes, supernodeSorter, permutedRegularID); + PermuteArrayWithMaskedIndex( + contourTree.Supernodes, supernodeSorter, permutedRegularID); #ifdef DEBUG_PRINT std::cout << "VI A. Sorted into Branches" << std::endl; diff --git a/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/BranchVolumeComparator.h b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/BranchVolumeComparator.h new file mode 100644 index 000000000..ddb9edc8f --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/BranchVolumeComparator.h @@ -0,0 +1,146 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_worklet_branch_decomposition_select_top_volume_contours_BranchVolumeComparator_h +#define vtk_m_filter_scalar_topology_worklet_branch_decomposition_select_top_volume_contours_BranchVolumeComparator_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace scalar_topology +{ +namespace select_top_volume_contours +{ + +using IdArrayType = vtkm::worklet::contourtree_augmented::IdArrayType; + +// Implementation of BranchVolumeComparator +class BranchVolumeComparatorImpl +{ +public: + using IdPortalType = typename IdArrayType::ReadPortalType; + + // constructor + VTKM_CONT + BranchVolumeComparatorImpl(const IdArrayType& branchRoots, + const IdArrayType& branchVolume, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) + : branchRootsPortal(branchRoots.PrepareForInput(device, token)) + , branchVolumePortal(branchVolume.PrepareForInput(device, token)) + { // constructor + } // constructor + + // () operator - gets called to do comparison + VTKM_EXEC + bool operator()(const vtkm::Id& i, const vtkm::Id& j) const + { // operator() + vtkm::Id volumeI = this->branchVolumePortal.Get(i); + vtkm::Id volumeJ = this->branchVolumePortal.Get(j); + + // primary sort on branch volume + if (volumeI > volumeJ) + return true; + if (volumeI < volumeJ) + return false; + + vtkm::Id branchI = + vtkm::worklet::contourtree_augmented::MaskedIndex(this->branchRootsPortal.Get(i)); + vtkm::Id branchJ = + vtkm::worklet::contourtree_augmented::MaskedIndex(this->branchRootsPortal.Get(j)); + + // secondary sort on branch ID + return (branchI < branchJ); + } // operator() + +private: + IdPortalType branchRootsPortal; + IdPortalType branchVolumePortal; + +}; // BranchVolumeComparatorImpl + +/// +/// Comparator of branch volume. Higher volume comes first +/// +class BranchVolumeComparator : public vtkm::cont::ExecutionObjectBase +{ + +public: + // constructor + VTKM_CONT + BranchVolumeComparator(const IdArrayType& branchRoots, const IdArrayType& branchVolume) + : BranchRoots(branchRoots) + , BranchVolume(branchVolume) + { + } + + VTKM_CONT BranchVolumeComparatorImpl PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const + { + return BranchVolumeComparatorImpl(this->BranchRoots, this->BranchVolume, device, token); + } + +private: + IdArrayType BranchRoots; + IdArrayType BranchVolume; +}; // BranchVolumeComparator + + +} // namespace select_top_volume_contours +} // namespace scalar_topology +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/CMakeLists.txt b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/CMakeLists.txt new file mode 100644 index 000000000..234966323 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/CMakeLists.txt @@ -0,0 +1,19 @@ +##============================================================================ +## 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. +##============================================================================ + +set(headers + ClarifyBranchEndSupernodeTypeWorklet.h + UpdateInfoByBranchDirectionWorklet.h + GetBranchVolumeWorklet.h + BranchVolumeComparator.h + ) +#----------------------------------------------------------------------------- + +vtkm_declare_headers(${headers}) diff --git a/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/ClarifyBranchEndSupernodeTypeWorklet.h b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/ClarifyBranchEndSupernodeTypeWorklet.h new file mode 100644 index 000000000..39ecdbdf3 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/ClarifyBranchEndSupernodeTypeWorklet.h @@ -0,0 +1,133 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_ClarifyBranchEndSupernodeTypeWorklet_h +#define vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_ClarifyBranchEndSupernodeTypeWorklet_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace scalar_topology +{ +namespace select_top_volume_contours +{ + +// For special branches that only have one superarc +// Clarify which end is leaf and which is saddle +class ClarifyBranchEndSupernodeTypeWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void( + FieldIn lowerSuperarcId, // (input) lower end superarc ID + FieldIn lowerIntrinsic, // (input) lower end superarc intrisic volume + FieldIn upperSuperarcId, // (input) upper end superarc ID + FieldIn upperIntrinsic, // (input) upper end superarc intrisic volume + FieldIn branchRoot, // (input) branch root superarc ID + FieldInOut isLowerLeaf, // (input/output) bool, whether the lower end is a leaf + FieldInOut isUpperLeaf // (input/output) bool, whether the upper end is a leaf + ); + using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, _7); + using InputDomain = _1; + + /// Constructor + VTKM_EXEC_CONT + ClarifyBranchEndSupernodeTypeWorklet(const vtkm::Id tVol) + : totalVolume(tVol) + { + } + + /// The functor checks the direction of the branch + VTKM_EXEC void operator()(const vtkm::Id& lowerSuperarcId, + const vtkm::Id& lowerIntrinsic, + const vtkm::Id& upperSuperarcId, + const vtkm::Id& upperIntrinsic, + const vtkm::Id& branchRoot, + bool& isLowerLeaf, + bool& isUpperLeaf) const + { + // do nothing: not a "leaf-leaf" branch + if (!isLowerLeaf || !isUpperLeaf) + return; + + // do nothing: actual leaf-leaf branch + if (lowerIntrinsic == totalVolume - 1 && lowerIntrinsic == upperIntrinsic) + return; + + // do something: fake leaf-leaf branch + // we already exclude the case of only one superarc above + vtkm::Id maskedLowerId = vtkm::worklet::contourtree_augmented::MaskedIndex(lowerSuperarcId); + vtkm::Id maskedUpperId = vtkm::worklet::contourtree_augmented::MaskedIndex(upperSuperarcId); + vtkm::Id maskedRoot = vtkm::worklet::contourtree_augmented::MaskedIndex(branchRoot); + + if (maskedLowerId == maskedRoot && maskedUpperId == maskedRoot) + { + const bool isAscending = vtkm::worklet::contourtree_augmented::IsAscending(lowerSuperarcId); + if (isAscending) + isUpperLeaf = false; + else + isLowerLeaf = false; + } + } + +private: + const vtkm::Id totalVolume; +}; // ClarifyBranchEndSupernodeTypeWorklet + +} // namespace select_top_volume_contours +} // namespace scalar_topology +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/GetBranchVolumeWorklet.h b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/GetBranchVolumeWorklet.h new file mode 100644 index 000000000..4d808f8ed --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/GetBranchVolumeWorklet.h @@ -0,0 +1,137 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_GetBranchVolumeWorklet_h +#define vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_GetBranchVolumeWorklet_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace scalar_topology +{ +namespace select_top_volume_contours +{ + +using IdArrayType = vtkm::worklet::contourtree_augmented::IdArrayType; + +/// +/// worklet to check the direction of branch +/// return true if the branch inner superarc points to the senior-most node +/// return false if the branch inner superarc comes from the senior-most node +/// trick: if the branch is the main branch, we return true for computation later +/// +class GetBranchVolumeWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void( + FieldIn lowerDirection, // (input) lower end superarc ID with direction information + FieldIn lowerIntrinsic, // (input) lower end superarc intrisic volume + FieldIn lowerDependent, // (input) lower end superarc dependent volume + FieldIn upperDirection, // (input) upper end superarc ID with direction information + FieldIn upperIntrinsic, // (input) upper end superarc intrisic volume + FieldIn upperDependent, // (input) upper end superarc dependent volume + FieldIn isLowerLeaf, // (input) bool, whether the lower end is a leaf + FieldIn isUpperLeaf, // (input) bool, whether the upper end is a leaf + FieldOut branchVolume // (output) volume of the branch + ); + using ExecutionSignature = _9(_1, _2, _3, _4, _5, _6, _7, _8); + using InputDomain = _1; + + /// Constructor + VTKM_EXEC_CONT + GetBranchVolumeWorklet(const vtkm::Id tVol) + : totalVolume(tVol) + { + } + + /// The functor checks the direction of the branch + VTKM_EXEC vtkm::Id operator()(const vtkm::Id& lowerDirection, + const vtkm::Id& lowerIntrinsic, + const vtkm::Id& lowerDependent, + const vtkm::Id& upperDirection, + const vtkm::Id& upperIntrinsic, + const vtkm::Id& upperDependent, + const bool& isLowerLeaf, + const bool& isUpperLeaf) const + { + if (isLowerLeaf && isUpperLeaf) + return totalVolume; + // if the branch is a minimum-saddle branch + // if the upper end superarc direction is pointing up, then dependent; otherwise, reverse + if (isLowerLeaf) + return contourtree_augmented::IsAscending(upperDirection) + ? upperDependent + : totalVolume - upperDependent + upperIntrinsic; + // if the branch is a maximum-saddle branch + // if the lower end superarc direction is pointing down, then true; otherwise, false + if (isUpperLeaf) + return !contourtree_augmented::IsAscending(lowerDirection) + ? lowerDependent + : totalVolume - lowerDependent + lowerIntrinsic; + + // in case of fallout, should never reach + return 0; + } + +private: + const vtkm::Id totalVolume; +}; // GetBranchVolumeWorklet + +} // namespace select_top_volume_contours +} // namespace scalar_topology +} // namespace worklet +} // namespace vtkm + +#endif diff --git a/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/UpdateInfoByBranchDirectionWorklet.h b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/UpdateInfoByBranchDirectionWorklet.h new file mode 100644 index 000000000..39568d0b7 --- /dev/null +++ b/vtkm/filter/scalar_topology/worklet/select_top_volume_contours/UpdateInfoByBranchDirectionWorklet.h @@ -0,0 +1,124 @@ +//============================================================================ +// 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 (c) 2018, The Regents of the University of California, through +// Lawrence Berkeley National Laboratory (subject to receipt of any required approvals +// from the U.S. Dept. of Energy). All rights reserved. +// +// 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 the University of California, Lawrence Berkeley National +// Laboratory, U.S. Dept. of Energy 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 THE COPYRIGHT HOLDERS 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 THE COPYRIGHT OWNER 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 an extension of 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. +// +// The PPP2 algorithm and software were jointly developed by +// Hamish Carr (University of Leeds), Gunther H. Weber (LBNL), and +// Oliver Ruebel (LBNL) +//============================================================================== + +#ifndef vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_UpdateInfoByBranchDirectionWorklet_h +#define vtk_m_filter_scalar_topology_worklet_select_top_volume_contours_UpdateInfoByBranchDirectionWorklet_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace scalar_topology +{ +namespace select_top_volume_contours +{ + +/// +/// 1. Get the saddle end isovalue +/// 2. Get Epsilon direction near the branch saddle end +/// If main branch, epsilon is 0. +/// Otherwise, -1 if the branch is the lower leaf branch, or 1 if upper leaf branch. +/// +template +class UpdateInfoByBranchDirectionWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void( + FieldIn isLowerLeaf, // (input) bool, whether the lower end is a leaf + FieldIn isUpperLeaf, // (input) bool, whether the upper end is a leaf + FieldIn upperEndValue, // (input) value type, the branch upper end value + FieldIn lowerEndValue, // (input) value type, the branch lower end value + FieldOut saddleEpsilon, // (output) vtkm::Id, epsilon direction around the saddle isovalue + FieldOut saddleValue // (output) value type, the saddle isovalue + ); + using ExecutionSignature = void(_1, _2, _3, _4, _5, _6); + using InputDomain = _1; + + /// Constructor + VTKM_EXEC_CONT + UpdateInfoByBranchDirectionWorklet() + { // constructor + } // constructor + + /// The functor returns the isovalue and the epsilon direction around the saddle end of the branch + VTKM_EXEC void operator()(const bool isLowerLeaf, + const bool isUpperLeaf, + const ValueType upperEndValue, + const ValueType lowerEndValue, + vtkm::Id& saddleEpsilon, + ValueType& saddleValue) const + { + // NOTE: for the main branch, the saddle value is undefined, + // because both upper and lower ends are leaf nodes. + // Let's use upperEndValue here to make the output not random. + if (isLowerLeaf && isUpperLeaf) + { + saddleEpsilon = 0; + saddleValue = upperEndValue; + } + else + { + saddleEpsilon = isLowerLeaf ? -1 : 1; + saddleValue = isLowerLeaf ? upperEndValue : lowerEndValue; + } + } +}; // EpsilonFromBranchDirection + + +} // namespace select_top_volume_contours +} // namespace scalar_topology +} // namespace worklet +} // namespace vtkm + +#endif