diff --git a/benchmarking/BenchmarkFilters.cxx b/benchmarking/BenchmarkFilters.cxx index dcc11f58d..fc40aa08d 100644 --- a/benchmarking/BenchmarkFilters.cxx +++ b/benchmarking/BenchmarkFilters.cxx @@ -438,8 +438,7 @@ void BenchContour(::benchmark::State& state) filter.SetMergeDuplicatePoints(mergePoints); filter.SetGenerateNormals(normals); - filter.SetComputeFastNormalsForStructured(fastNormals); - filter.SetComputeFastNormalsForUnstructured(fastNormals); + filter.SetComputeFastNormals(fastNormals); vtkm::cont::Timer timer{ device }; diff --git a/benchmarking/BenchmarkInSitu.cxx b/benchmarking/BenchmarkInSitu.cxx index d3bc990d8..3b629ffa9 100644 --- a/benchmarking/BenchmarkInSitu.cxx +++ b/benchmarking/BenchmarkInSitu.cxx @@ -340,8 +340,7 @@ void BenchContour(::benchmark::State& state) filter.SetActiveField(PointScalarsName, vtkm::cont::Field::Association::Points); filter.SetMergeDuplicatePoints(true); filter.SetGenerateNormals(true); - filter.SetComputeFastNormalsForStructured(true); - filter.SetComputeFastNormalsForUnstructured(true); + filter.SetComputeFastNormals(true); state.ResumeTiming(); // And resume timers. filterTimer.Start(); diff --git a/docs/changelog/split-contour-filters.md b/docs/changelog/split-contour-filters.md new file mode 100644 index 000000000..3dd58bfab --- /dev/null +++ b/docs/changelog/split-contour-filters.md @@ -0,0 +1,11 @@ +# Split flying edges and marching cells into separate filters + +The contour filter contains 2 separate implementations, Marching Cells and Flying Edges, the latter only available if the input has a `CellSetStructured<3>` and `ArrayHandleUniformPointCoordinates` for point coordinates. The compilation of this filter was lenghty and resource-heavy, because both algorithms were part of the same translation unit. + +Now, this filter is separated into two new filters, `ContourFlyingEdges` and `ContourMarchingCells`, compiling more efficiently into two translation units. The `Contour` API is left unchanged. All 3 filters `Contour`, `ContourFlyingEdges` and `ContourMarchingCells` rely on a new abstract class `AbstractContour` to provide configuration and common utility functions. + +Although `Contour` is still the preferred option for most cases because it selects the best implementation according to the input, `ContourMarchingCells` is usable on any kind of 3D Dataset. For now, `ContourFlyingEdges` operates only on structured uniform datasets. + +Deprecate functions `GetComputeFastNormalsForStructured`, `SetComputeFastNormalsForStructured`, `GetComputeFastNormalsForUnstructured` and `GetComputeFastNormalsForUnstructured`, to use the more general `GetComputeFastNormals` and `SetComputeFastNormals` instead. + + By default, for the `Contour` filter, `GenerateNormals` is now `true`, and `ComputeFastNormals` is `false`. diff --git a/vtkm/filter/contour/AbstractContour.h b/vtkm/filter/contour/AbstractContour.h new file mode 100644 index 000000000..fad11af86 --- /dev/null +++ b/vtkm/filter/contour/AbstractContour.h @@ -0,0 +1,181 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_filter_contour_AbstractContour_h +#define vtk_m_filter_contour_AbstractContour_h + +#include +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace contour +{ +/// \brief Contour filter interface +/// +/// Provides common configuration & execution methods for contour filters +/// Only the method \c DoExecute executing the contour algorithm needs to be implemented +class VTKM_FILTER_CONTOUR_EXPORT AbstractContour : public vtkm::filter::FilterField +{ +public: + void SetNumberOfIsoValues(vtkm::Id num) + { + if (num >= 0) + { + this->IsoValues.resize(static_cast(num)); + } + } + + vtkm::Id GetNumberOfIsoValues() const { return static_cast(this->IsoValues.size()); } + + void SetIsoValue(vtkm::Float64 v) { this->SetIsoValue(0, v); } + + void SetIsoValue(vtkm::Id index, vtkm::Float64 v) + { + std::size_t i = static_cast(index); + if (i >= this->IsoValues.size()) + { + this->IsoValues.resize(i + 1); + } + this->IsoValues[i] = v; + } + + void SetIsoValues(const std::vector& values) { this->IsoValues = values; } + + vtkm::Float64 GetIsoValue(vtkm::Id index) const + { + return this->IsoValues[static_cast(index)]; + } + + /// Set/Get whether normals should be generated. Off by default. + VTKM_CONT + void SetGenerateNormals(bool on) { this->GenerateNormals = on; } + VTKM_CONT + bool GetGenerateNormals() const { return this->GenerateNormals; } + + /// Set/Get whether to append the ids of the intersected edges to the vertices of the isosurface triangles. Off by default. + VTKM_CONT + void SetAddInterpolationEdgeIds(bool on) { this->AddInterpolationEdgeIds = on; } + VTKM_CONT + bool GetAddInterpolationEdgeIds() const { return this->AddInterpolationEdgeIds; } + + /// Set/Get whether the fast path should be used for normals computation. Off by default. + VTKM_CONT + void SetComputeFastNormals(bool on) { this->ComputeFastNormals = on; } + VTKM_CONT + bool GetComputeFastNormals() const { return this->ComputeFastNormals; } + + VTKM_CONT + void SetNormalArrayName(const std::string& name) { this->NormalArrayName = name; } + + VTKM_CONT + const std::string& GetNormalArrayName() const { return this->NormalArrayName; } + + /// Set/Get whether the points generated should be unique for every triangle + /// or will duplicate points be merged together. Duplicate points are identified + /// by the unique edge it was generated from. + /// + VTKM_CONT + void SetMergeDuplicatePoints(bool on) { this->MergeDuplicatedPoints = on; } + + VTKM_CONT + bool GetMergeDuplicatePoints() { return this->MergeDuplicatedPoints; } + +protected: + /// \brief Map a given field to the output \c DataSet , depending on its type. + /// + /// The worklet needs to implement \c ProcessPointField to process point fields as arrays + /// and \c GetCellIdMap function giving the cell id mapping from input to output + template + VTKM_CONT static bool DoMapField(vtkm::cont::DataSet& result, + const vtkm::cont::Field& field, + WorkletType& worklet) + { + if (field.IsPointField()) + { + vtkm::cont::UnknownArrayHandle inputArray = field.GetData(); + vtkm::cont::UnknownArrayHandle outputArray = inputArray.NewInstanceBasic(); + + auto functor = [&](const auto& concrete) { + using ComponentType = typename std::decay_t::ValueType::ComponentType; + auto fieldArray = outputArray.ExtractArrayFromComponents(); + worklet.ProcessPointField(concrete, fieldArray); + }; + inputArray.CastAndCallWithExtractedArray(functor); + result.AddPointField(field.GetName(), outputArray); + return true; + } + else if (field.IsCellField()) + { + // Use the precompiled field permutation function. + vtkm::cont::ArrayHandle permutation = worklet.GetCellIdMap(); + return vtkm::filter::MapFieldPermutation(field, permutation, result); + } + else if (field.IsWholeDataSetField()) + { + result.AddField(field); + return true; + } + return false; + } + + VTKM_CONT void ExecuteGenerateNormals(vtkm::cont::DataSet& output, + const vtkm::cont::ArrayHandle& normals) + { + if (this->GenerateNormals) + { + if (this->GetComputeFastNormals()) + { + vtkm::filter::vector_analysis::SurfaceNormals surfaceNormals; + surfaceNormals.SetPointNormalsName(this->NormalArrayName); + surfaceNormals.SetGeneratePointNormals(true); + output = surfaceNormals.Execute(output); + } + else + { + output.AddField(vtkm::cont::make_FieldPoint(this->NormalArrayName, normals)); + } + } + } + + template + VTKM_CONT void ExecuteAddInterpolationEdgeIds(vtkm::cont::DataSet& output, WorkletType& worklet) + { + if (this->AddInterpolationEdgeIds) + { + vtkm::cont::Field interpolationEdgeIdsField(this->InterpolationEdgeIdsArrayName, + vtkm::cont::Field::Association::Points, + worklet.GetInterpolationEdgeIds()); + output.AddField(interpolationEdgeIdsField); + } + } + + VTKM_CONT + virtual vtkm::cont::DataSet DoExecute( + const vtkm::cont::DataSet& result) = 0; // Needs to be overridden by contour implementations + + std::vector IsoValues; + bool GenerateNormals = true; + bool ComputeFastNormals = false; + + bool AddInterpolationEdgeIds = false; + bool MergeDuplicatedPoints = true; + std::string NormalArrayName = "normals"; + std::string InterpolationEdgeIdsArrayName = "edgeIds"; +}; +} // namespace contour +} // namespace filter +} // namespace vtkm + +#endif // vtk_m_filter_contour_AbstractContour_h diff --git a/vtkm/filter/contour/CMakeLists.txt b/vtkm/filter/contour/CMakeLists.txt index 351d55546..e9a85badd 100644 --- a/vtkm/filter/contour/CMakeLists.txt +++ b/vtkm/filter/contour/CMakeLists.txt @@ -9,9 +9,12 @@ ##============================================================================ set(contour_headers + AbstractContour.h ClipWithField.h ClipWithImplicitFunction.h Contour.h + ContourFlyingEdges.h + ContourMarchingCells.h MIRFilter.h Slice.h ) @@ -19,15 +22,23 @@ set(contour_headers set(contour_sources_device ClipWithField.cxx ClipWithImplicitFunction.cxx - Contour.cxx + ContourFlyingEdges.cxx + ContourMarchingCells.cxx MIRFilter.cxx Slice.cxx ) +set(contour_sources + # Contour defers worklet compilation to other filters, + # so it does not need to be compiled with a device adapter. + Contour.cxx +) + set_source_files_properties(Contour.cxx PROPERTIES SKIP_UNITY_BUILD_INCLUSION ON) vtkm_library( NAME vtkm_filter_contour + SOURCES ${contour_sources} HEADERS ${contour_headers} DEVICE_SOURCES ${contour_sources_device} USE_VTKM_JOB_POOL diff --git a/vtkm/filter/contour/Contour.cxx b/vtkm/filter/contour/Contour.cxx index 6b75a9c2c..ebbd896c0 100644 --- a/vtkm/filter/contour/Contour.cxx +++ b/vtkm/filter/contour/Contour.cxx @@ -7,16 +7,13 @@ // 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 #include -#include -#include +#include +#include namespace vtkm { @@ -25,155 +22,48 @@ namespace filter using SupportedTypes = vtkm::List; -namespace -{ - -inline bool IsCellSetStructured(const vtkm::cont::UnknownCellSet& cellset) -{ - if (cellset.template IsType>() || - cellset.template IsType>() || - cellset.template IsType>()) - { - return true; - } - return false; -} - -VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result, - const vtkm::cont::Field& field, - vtkm::worklet::Contour& worklet) -{ - if (field.IsPointField()) - { - vtkm::cont::UnknownArrayHandle inputArray = field.GetData(); - vtkm::cont::UnknownArrayHandle outputArray = inputArray.NewInstanceBasic(); - - auto functor = [&](const auto& concrete) { - using ComponentType = typename std::decay_t::ValueType::ComponentType; - auto fieldArray = outputArray.ExtractArrayFromComponents(); - worklet.ProcessPointField(concrete, fieldArray); - }; - inputArray.CastAndCallWithExtractedArray(functor); - result.AddPointField(field.GetName(), outputArray); - return true; - } - else if (field.IsCellField()) - { - // Use the precompiled field permutation function. - vtkm::cont::ArrayHandle permutation = worklet.GetCellIdMap(); - return vtkm::filter::MapFieldPermutation(field, permutation, result); - } - else if (field.IsWholeDataSetField()) - { - result.AddField(field); - return true; - } - else - { - return false; - } -} - -} // anonymous namespace - namespace contour { -//----------------------------------------------------------------------------- -void Contour::SetMergeDuplicatePoints(bool on) -{ - this->MergeDuplicatedPoints = on; -} - -VTKM_CONT -bool Contour::GetMergeDuplicatePoints() const -{ - return MergeDuplicatedPoints; -} - //----------------------------------------------------------------------------- vtkm::cont::DataSet Contour::DoExecute(const vtkm::cont::DataSet& inDataSet) { - vtkm::worklet::Contour worklet; - worklet.SetMergeDuplicatePoints(this->GetMergeDuplicatePoints()); + // Switch between Marching Cubes and Flying Edges implementation of contour, + // depending on the type of CellSet we are processing - if (!this->GetFieldFromDataSet(inDataSet).IsPointField()) + vtkm::cont::UnknownCellSet inCellSet = inDataSet.GetCellSet(); + auto inCoords = inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()).GetData(); + std::unique_ptr implementation; + + // For now, Flying Edges is only used for 3D Structured CellSets, + // using Uniform coordinates. + if (inCellSet.template IsType>() && + inCoords.template IsType< + vtkm::cont::ArrayHandle>()) { - throw vtkm::cont::ErrorFilterExecution("Point field expected."); + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Using flying edges"); + implementation.reset(new vtkm::filter::contour::ContourFlyingEdges); + implementation->SetComputeFastNormals(this->GetComputeFastNormals()); + } + else + { + VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Using marching cells"); + implementation.reset(new vtkm::filter::contour::ContourMarchingCells); + implementation->SetComputeFastNormals(this->GetComputeFastNormals()); } - if (this->IsoValues.empty()) + implementation->SetMergeDuplicatePoints(this->GetMergeDuplicatePoints()); + implementation->SetGenerateNormals(this->GetGenerateNormals()); + implementation->SetAddInterpolationEdgeIds(this->GetAddInterpolationEdgeIds()); + implementation->SetNormalArrayName(this->GetNormalArrayName()); + implementation->SetActiveField(this->GetActiveFieldName()); + implementation->SetFieldsToPass(this->GetFieldsToPass()); + implementation->SetNumberOfIsoValues(this->GetNumberOfIsoValues()); + for (int i = 0; i < this->GetNumberOfIsoValues(); i++) { - throw vtkm::cont::ErrorFilterExecution("No iso-values provided."); + implementation->SetIsoValue(i, this->GetIsoValue(i)); } - //get the inputCells and coordinates of the dataset - const vtkm::cont::UnknownCellSet& inputCells = inDataSet.GetCellSet(); - const vtkm::cont::CoordinateSystem& inputCoords = - inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()); - - using Vec3HandleType = vtkm::cont::ArrayHandle; - Vec3HandleType vertices; - Vec3HandleType normals; - - vtkm::cont::CellSetSingleType<> outputCells; - - bool generateHighQualityNormals = IsCellSetStructured(inputCells) - ? !this->ComputeFastNormalsForStructured - : !this->ComputeFastNormalsForUnstructured; - - auto resolveFieldType = [&](const auto& concrete) { - // use std::decay to remove const ref from the decltype of concrete. - using T = typename std::decay_t::ValueType; - std::vector ivalues(this->IsoValues.size()); - for (std::size_t i = 0; i < ivalues.size(); ++i) - { - ivalues[i] = static_cast(this->IsoValues[i]); - } - - if (this->GenerateNormals && generateHighQualityNormals) - { - outputCells = - worklet.Run(ivalues, inputCells, inputCoords.GetData(), concrete, vertices, normals); - } - else - { - outputCells = worklet.Run(ivalues, inputCells, inputCoords.GetData(), concrete, vertices); - } - }; - - this->GetFieldFromDataSet(inDataSet) - .GetData() - .CastAndCallForTypesWithFloatFallback( - resolveFieldType); - - auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); }; - vtkm::cont::DataSet output = this->CreateResultCoordinateSystem( - inDataSet, outputCells, inputCoords.GetName(), vertices, mapper); - - if (this->GenerateNormals) - { - if (!generateHighQualityNormals) - { - vtkm::filter::vector_analysis::SurfaceNormals surfaceNormals; - surfaceNormals.SetPointNormalsName(this->NormalArrayName); - surfaceNormals.SetGeneratePointNormals(true); - output = surfaceNormals.Execute(output); - } - else - { - output.AddField(vtkm::cont::make_FieldPoint(this->NormalArrayName, normals)); - } - } - - if (this->AddInterpolationEdgeIds) - { - vtkm::cont::Field interpolationEdgeIdsField(InterpolationEdgeIdsArrayName, - vtkm::cont::Field::Association::Points, - worklet.GetInterpolationEdgeIds()); - output.AddField(interpolationEdgeIdsField); - } - - return output; + return implementation->Execute(inDataSet); } } // namespace contour } // namespace filter diff --git a/vtkm/filter/contour/Contour.h b/vtkm/filter/contour/Contour.h index 6ef69d9d2..615f10cf0 100644 --- a/vtkm/filter/contour/Contour.h +++ b/vtkm/filter/contour/Contour.h @@ -11,7 +11,7 @@ #ifndef vtk_m_filter_contour_Contour_h #define vtk_m_filter_contour_Contour_h -#include +#include #include namespace vtkm @@ -25,103 +25,33 @@ namespace contour /// Takes as input a volume (e.g., 3D structured point set) and generates on /// output one or more isosurfaces. /// Multiple contour values must be specified to generate the isosurfaces. +/// This filter automatically selects the right implmentation between Marching Cells +/// and Flying Edges algorithms depending on the type of input \c DataSet : Flying Edges +/// is only available for 3-dimensional datasets using uniform point coordinates. /// @warning /// This filter is currently only supports 3D volumes. -class VTKM_FILTER_CONTOUR_EXPORT Contour : public vtkm::filter::FilterField +class VTKM_FILTER_CONTOUR_EXPORT Contour : public vtkm::filter::contour::AbstractContour { public: - void SetNumberOfIsoValues(vtkm::Id num) - { - if (num >= 0) - { - this->IsoValues.resize(static_cast(num)); - } - } - - vtkm::Id GetNumberOfIsoValues() const { return static_cast(this->IsoValues.size()); } - - void SetIsoValue(vtkm::Float64 v) { this->SetIsoValue(0, v); } - - void SetIsoValue(vtkm::Id index, vtkm::Float64 v) - { - std::size_t i = static_cast(index); - if (i >= this->IsoValues.size()) - { - this->IsoValues.resize(i + 1); - } - this->IsoValues[i] = v; - } - - void SetIsoValues(const std::vector& values) { this->IsoValues = values; } - - vtkm::Float64 GetIsoValue(vtkm::Id index) const - { - return this->IsoValues[static_cast(index)]; - } - - /// Set/Get whether the points generated should be unique for every triangle - /// or will duplicate points be merged together. Duplicate points are identified - /// by the unique edge it was generated from. - /// - VTKM_CONT - void SetMergeDuplicatePoints(bool on); - - VTKM_CONT - bool GetMergeDuplicatePoints() const; - - /// Set/Get whether normals should be generated. Off by default. If enabled, - /// the default behaviour is to generate high quality normals for structured - /// datasets, using gradients, and generate fast normals for unstructured - /// datasets based on the result triangle mesh. - /// - VTKM_CONT - void SetGenerateNormals(bool on) { this->GenerateNormals = on; } - VTKM_CONT - bool GetGenerateNormals() const { return this->GenerateNormals; } - - /// Set/Get whether to append the ids of the intersected edges to the vertices of the isosurface triangles. Off by default. - VTKM_CONT - void SetAddInterpolationEdgeIds(bool on) { this->AddInterpolationEdgeIds = on; } - VTKM_CONT - bool GetAddInterpolationEdgeIds() const { return this->AddInterpolationEdgeIds; } - /// Set/Get whether the fast path should be used for normals computation for /// structured datasets. Off by default. - VTKM_CONT - void SetComputeFastNormalsForStructured(bool on) { this->ComputeFastNormalsForStructured = on; } - VTKM_CONT - bool GetComputeFastNormalsForStructured() const { return this->ComputeFastNormalsForStructured; } + VTKM_DEPRECATED(2.1, "Use SetComputeFastNormals.") + VTKM_CONT void SetComputeFastNormalsForStructured(bool on) { this->SetComputeFastNormals(on); } + VTKM_DEPRECATED(2.1, "Use GetComputeFastNormals.") + VTKM_CONT bool GetComputeFastNormalsForStructured() const + { + return this->GetComputeFastNormals(); + } /// Set/Get whether the fast path should be used for normals computation for /// unstructured datasets. On by default. - VTKM_CONT - void SetComputeFastNormalsForUnstructured(bool on) + VTKM_DEPRECATED(2.1, "Use SetComputeFastNormals.") + VTKM_CONT void SetComputeFastNormalsForUnstructured(bool on) { this->SetComputeFastNormals(on); } + VTKM_DEPRECATED(2.1, "Use GetComputeFastNormals.") + VTKM_CONT bool GetComputeFastNormalsForUnstructured() const { - this->ComputeFastNormalsForUnstructured = on; + return this->GetComputeFastNormals(); } - VTKM_CONT - bool GetComputeFastNormalsForUnstructured() const - { - return this->ComputeFastNormalsForUnstructured; - } - - VTKM_CONT - void SetNormalArrayName(const std::string& name) { this->NormalArrayName = name; } - - VTKM_CONT - const std::string& GetNormalArrayName() const { return this->NormalArrayName; } - -private: - VTKM_CONT - - std::vector IsoValues; - bool GenerateNormals = false; - bool AddInterpolationEdgeIds = false; - bool ComputeFastNormalsForStructured = false; - bool ComputeFastNormalsForUnstructured = true; - bool MergeDuplicatedPoints = true; - std::string NormalArrayName = "normals"; - std::string InterpolationEdgeIdsArrayName = "edgeIds"; protected: // Needed by the subclass Slice diff --git a/vtkm/filter/contour/ContourFlyingEdges.cxx b/vtkm/filter/contour/ContourFlyingEdges.cxx new file mode 100644 index 000000000..f457a0be8 --- /dev/null +++ b/vtkm/filter/contour/ContourFlyingEdges.cxx @@ -0,0 +1,112 @@ +//============================================================================ +// 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 + +namespace vtkm +{ +namespace filter +{ + +using SupportedTypes = vtkm::List; + +namespace contour +{ +//----------------------------------------------------------------------------- +vtkm::cont::DataSet ContourFlyingEdges::DoExecute(const vtkm::cont::DataSet& inDataSet) +{ + vtkm::worklet::ContourFlyingEdges worklet; + worklet.SetMergeDuplicatePoints(this->GetMergeDuplicatePoints()); + + if (!this->GetFieldFromDataSet(inDataSet).IsPointField()) + { + throw vtkm::cont::ErrorFilterExecution("Point field expected."); + } + + if (this->IsoValues.empty()) + { + throw vtkm::cont::ErrorFilterExecution("No iso-values provided."); + } + + vtkm::cont::UnknownCellSet inCellSet = inDataSet.GetCellSet(); + const vtkm::cont::CoordinateSystem& inCoords = + inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()); + + if (!inCellSet.template IsType>() || + !inCoords.GetData() + .template IsType< + vtkm::cont::ArrayHandle>()) + { + throw vtkm::cont::ErrorFilterExecution("This filter is only available for 3-Dimensional " + "Structured Cell Sets using uniform point coordinates."); + } + + // Get the CellSet's known dynamic type + const vtkm::cont::CellSetStructured<3>& inputCells = + inDataSet.GetCellSet().AsCellSet>(); + + using Vec3HandleType = vtkm::cont::ArrayHandle; + Vec3HandleType vertices; + Vec3HandleType normals; + + vtkm::cont::CellSetSingleType<> outputCells; + + auto resolveFieldType = [&](const auto& concrete) { + // use std::decay to remove const ref from the decltype of concrete. + using T = typename std::decay_t::ValueType; + std::vector ivalues(this->IsoValues.size()); + for (std::size_t i = 0; i < ivalues.size(); ++i) + { + ivalues[i] = static_cast(this->IsoValues[i]); + } + + if (this->GenerateNormals && !this->GetComputeFastNormals()) + { + outputCells = worklet.Run( + ivalues, + inputCells, + inCoords.GetData().AsArrayHandle(), + concrete, + vertices, + normals); + } + else + { + outputCells = worklet.Run( + ivalues, + inputCells, + inCoords.GetData().AsArrayHandle(), + concrete, + vertices); + } + }; + + this->GetFieldFromDataSet(inDataSet) + .GetData() + .CastAndCallForTypesWithFloatFallback( + resolveFieldType); + + auto mapper = [&](auto& result, const auto& f) { this->DoMapField(result, f, worklet); }; + vtkm::cont::DataSet output = this->CreateResultCoordinateSystem( + inDataSet, outputCells, inCoords.GetName(), vertices, mapper); + + this->ExecuteGenerateNormals(output, normals); + this->ExecuteAddInterpolationEdgeIds(output, worklet); + + return output; +} +} // namespace contour +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/contour/ContourFlyingEdges.h b/vtkm/filter/contour/ContourFlyingEdges.h new file mode 100644 index 000000000..cd26d3300 --- /dev/null +++ b/vtkm/filter/contour/ContourFlyingEdges.h @@ -0,0 +1,41 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_filter_contour_ContourFlyingEdges_h +#define vtk_m_filter_contour_ContourFlyingEdges_h + +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace contour +{ +/// \brief generate isosurface(s) from a 3-dimensional structured mesh + +/// Takes as input a 3D structured mesh and generates on +/// output one or more isosurfaces using the Flying Edges algorithm. +/// Multiple contour values must be specified to generate the isosurfaces. +/// +/// This implementation only accepts \c CellSetStructured<3> inputs using +/// \c ArrayHandleUniformPointCoordinates for point coordinates, +/// and is only used as part of the more general \c Contour filter +class VTKM_FILTER_CONTOUR_EXPORT ContourFlyingEdges : public vtkm::filter::contour::AbstractContour +{ +protected: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& result) override; +}; +} // namespace contour +} // namespace filter +} // namespace vtkm + +#endif // vtk_m_filter_contour_ContourFlyingEdges_h diff --git a/vtkm/filter/contour/ContourMarchingCells.cxx b/vtkm/filter/contour/ContourMarchingCells.cxx new file mode 100644 index 000000000..855e5b154 --- /dev/null +++ b/vtkm/filter/contour/ContourMarchingCells.cxx @@ -0,0 +1,85 @@ +//============================================================================ +// 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 + + +namespace vtkm +{ +namespace filter +{ +namespace contour +{ +//----------------------------------------------------------------------------- +vtkm::cont::DataSet ContourMarchingCells::DoExecute(const vtkm::cont::DataSet& inDataSet) +{ + vtkm::worklet::ContourMarchingCells worklet; + worklet.SetMergeDuplicatePoints(this->GetMergeDuplicatePoints()); + + if (!this->GetFieldFromDataSet(inDataSet).IsPointField()) + { + throw vtkm::cont::ErrorFilterExecution("Point field expected."); + } + + if (this->IsoValues.empty()) + { + throw vtkm::cont::ErrorFilterExecution("No iso-values provided."); + } + + //get the inputCells and coordinates of the dataset + const vtkm::cont::UnknownCellSet& inputCells = inDataSet.GetCellSet(); + const vtkm::cont::CoordinateSystem& inputCoords = + inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()); + + using Vec3HandleType = vtkm::cont::ArrayHandle; + Vec3HandleType vertices; + Vec3HandleType normals; + + vtkm::cont::CellSetSingleType<> outputCells; + + auto resolveFieldType = [&](const auto& concrete) { + // use std::decay to remove const ref from the decltype of concrete. + using T = typename std::decay_t::ValueType; + std::vector ivalues(this->IsoValues.size()); + for (std::size_t i = 0; i < ivalues.size(); ++i) + { + ivalues[i] = static_cast(this->IsoValues[i]); + } + + if (this->GenerateNormals && !this->GetComputeFastNormals()) + { + outputCells = + worklet.Run(ivalues, inputCells, inputCoords.GetData(), concrete, vertices, normals); + } + else + { + outputCells = worklet.Run(ivalues, inputCells, inputCoords.GetData(), concrete, vertices); + } + }; + + this->CastAndCallScalarField(this->GetFieldFromDataSet(inDataSet), resolveFieldType); + + auto mapper = [&](auto& result, const auto& f) { this->DoMapField(result, f, worklet); }; + vtkm::cont::DataSet output = this->CreateResultCoordinateSystem( + inDataSet, outputCells, inputCoords.GetName(), vertices, mapper); + + this->ExecuteGenerateNormals(output, normals); + this->ExecuteAddInterpolationEdgeIds(output, worklet); + + + return output; +} +} // namespace contour +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/contour/ContourMarchingCells.h b/vtkm/filter/contour/ContourMarchingCells.h new file mode 100644 index 000000000..653a1eaad --- /dev/null +++ b/vtkm/filter/contour/ContourMarchingCells.h @@ -0,0 +1,43 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_filter_contour_ContourMarchingCells_h +#define vtk_m_filter_contour_ContourMarchingCells_h + +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace contour +{ +/// \brief generate isosurface(s) from a Volume using the Marching Cells algorithm +/// +/// Takes as input a volume (e.g., 3D structured point set) and generates on +/// output one or more isosurfaces. +/// Multiple contour values must be specified to generate the isosurfaces. +/// +/// This implementation is not optimized for all use cases, it is used by +/// the more general \c Contour filter which selects the best implementation +/// for all types of \c DataSet . . +class VTKM_FILTER_CONTOUR_EXPORT ContourMarchingCells + : public vtkm::filter::contour::AbstractContour +{ +protected: + VTKM_CONT + vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& result) override; +}; +} // namespace contour +} // namespace filter +} // namespace vtkm + +#endif // vtk_m_filter_contour_ContourMarchingCells_h diff --git a/vtkm/filter/contour/testing/UnitTestContourFilter.cxx b/vtkm/filter/contour/testing/UnitTestContourFilter.cxx index 61e41e5c6..d1e159ce5 100644 --- a/vtkm/filter/contour/testing/UnitTestContourFilter.cxx +++ b/vtkm/filter/contour/testing/UnitTestContourFilter.cxx @@ -11,10 +11,13 @@ #include #include #include +#include #include #include #include +#include +#include #include #include @@ -26,7 +29,8 @@ namespace class TestContourFilter { public: - void TestContourUniformGrid() const + template + void TestContourUniformGrid(vtkm::IdComponent numPointsNoMergeDuplicate) const { std::cout << "Testing Contour filter on a uniform grid" << std::endl; @@ -38,14 +42,14 @@ public: genIds.SetCellFieldName("cellvar"); vtkm::cont::DataSet dataSet = genIds.Execute(tangle.Execute()); - vtkm::filter::contour::Contour mc; + ContourFilterType filter; - mc.SetGenerateNormals(true); - mc.SetIsoValue(0, 0.5); - mc.SetActiveField("tangle"); - mc.SetFieldsToPass(vtkm::filter::FieldSelection::Mode::None); + filter.SetGenerateNormals(true); + filter.SetIsoValue(0, 0.5); + filter.SetActiveField("tangle"); + filter.SetFieldsToPass(vtkm::filter::FieldSelection::Mode::None); - auto result = mc.Execute(dataSet); + auto result = filter.Execute(dataSet); { VTKM_TEST_ASSERT(result.GetNumberOfCoordinateSystems() == 1, "Wrong number of coordinate systems in the output dataset"); @@ -55,8 +59,8 @@ public: } // let's execute with mapping fields. - mc.SetFieldsToPass({ "tangle", "cellvar" }); - result = mc.Execute(dataSet); + filter.SetFieldsToPass({ "tangle", "cellvar" }); + result = filter.Execute(dataSet); { const bool isMapped = result.HasField("tangle"); VTKM_TEST_ASSERT(isMapped, "mapping should pass"); @@ -99,16 +103,13 @@ public: VTKM_TEST_ASSERT(cells.GetNumberOfCells() == 160, ""); } - //Now try with vertex merging disabled. Since this - //we use FlyingEdges we now which does point merging for free - //so we should see the number of points not change - mc.SetMergeDuplicatePoints(false); - mc.SetFieldsToPass(vtkm::filter::FieldSelection::Mode::All); - result = mc.Execute(dataSet); + //Now try with vertex merging disabled. + filter.SetMergeDuplicatePoints(false); + filter.SetFieldsToPass(vtkm::filter::FieldSelection::Mode::All); + result = filter.Execute(dataSet); { vtkm::cont::CoordinateSystem coords = result.GetCoordinateSystem(); - - VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 72, + VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == numPointsNoMergeDuplicate, "Shouldn't have less coordinates than the unmerged version"); //verify that the number of cells is correct (160) @@ -120,20 +121,24 @@ public: } } + template void Test3DUniformDataSet0() const { vtkm::cont::testing::MakeTestDataSet maker; vtkm::cont::DataSet inputData = maker.Make3DUniformDataSet0(); std::string fieldName = "pointvar"; + // Defend the test against changes to Make3DUniformDataSet0(): VTKM_TEST_ASSERT(inputData.HasField(fieldName)); vtkm::cont::Field pointField = inputData.GetField(fieldName); + vtkm::Range range; pointField.GetRange(&range); vtkm::FloatDefault isovalue = 100.0; // Range = [10.1, 180.5] VTKM_TEST_ASSERT(range.Contains(isovalue)); - vtkm::filter::contour::Contour filter; + + ContourFilterType filter; filter.SetGenerateNormals(false); filter.SetMergeDuplicatePoints(true); filter.SetIsoValue(isovalue); @@ -143,6 +148,7 @@ public: VTKM_TEST_ASSERT(outputData.GetNumberOfPoints() == 9); } + template void TestContourWedges() const { std::cout << "Testing Contour filter on wedge cells" << std::endl; @@ -158,7 +164,7 @@ public: vtkm::cont::ArrayHandle fieldArray; dataSet.GetPointField("gyroid").GetData().AsArrayHandle(fieldArray); - vtkm::filter::contour::Contour isosurfaceFilter; + ContourFilterType isosurfaceFilter; isosurfaceFilter.SetActiveField("gyroid"); isosurfaceFilter.SetMergeDuplicatePoints(false); isosurfaceFilter.SetIsoValue(0.0); @@ -167,11 +173,54 @@ public: VTKM_TEST_ASSERT(result.GetNumberOfCells() == 52); } + void TestUnsupportedFlyingEdges() const + { + vtkm::cont::testing::MakeTestDataSet maker; + + vtkm::cont::DataSet rectilinearDataset = maker.Make3DRectilinearDataSet0(); + vtkm::cont::DataSet explicitDataSet = maker.Make3DExplicitDataSet0(); + + vtkm::filter::contour::ContourFlyingEdges filter; + filter.SetIsoValue(2.0); + filter.SetActiveField("pointvar"); + + try + { + filter.Execute(rectilinearDataset); + VTKM_TEST_FAIL("Flying Edges filter should not run on datasets with rectilinear coordinates"); + } + catch (vtkm::cont::ErrorFilterExecution&) + { + std::cout << "Execution successfully aborted" << std::endl; + } + + try + { + filter.Execute(explicitDataSet); + VTKM_TEST_FAIL("Flying Edges filter should not run on explicit datasets"); + } + catch (vtkm::cont::ErrorFilterExecution&) + { + std::cout << "Execution successfully aborted" << std::endl; + } + } + void operator()() const { - this->Test3DUniformDataSet0(); - this->TestContourUniformGrid(); - this->TestContourWedges(); + this->TestContourUniformGrid(72); + this->TestContourUniformGrid(72); + // Unlike flying edges, marching cells does not have point merging for free, + // So the number of points should increase when disabling duplicate point merging. + this->TestContourUniformGrid(480); + + this->Test3DUniformDataSet0(); + this->Test3DUniformDataSet0(); + this->Test3DUniformDataSet0(); + + this->TestContourWedges(); + this->TestContourWedges(); + + this->TestUnsupportedFlyingEdges(); } }; // class TestContourFilter diff --git a/vtkm/filter/contour/testing/UnitTestContourFilterNormals.cxx b/vtkm/filter/contour/testing/UnitTestContourFilterNormals.cxx index d9159b152..be37c3057 100644 --- a/vtkm/filter/contour/testing/UnitTestContourFilterNormals.cxx +++ b/vtkm/filter/contour/testing/UnitTestContourFilterNormals.cxx @@ -105,6 +105,14 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) vtkm::filter::contour::Contour mc; mc.SetIsoValue(0, 200); mc.SetGenerateNormals(true); + if (structured) + { + mc.SetComputeFastNormals(false); + } + else + { + mc.SetComputeFastNormals(true); + } // Test default normals generation: high quality for structured, fast for unstructured. auto expected = structured ? hq_sg : fast; @@ -136,7 +144,7 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) // Test the other normals generation method if (structured) { - mc.SetComputeFastNormalsForStructured(true); + mc.SetComputeFastNormals(true); expected = fast; if (using_fe_y_alg_ordering) { @@ -145,7 +153,7 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured) } else { - mc.SetComputeFastNormalsForUnstructured(false); + mc.SetComputeFastNormals(false); expected = hq_ug; } diff --git a/vtkm/filter/contour/worklet/CMakeLists.txt b/vtkm/filter/contour/worklet/CMakeLists.txt index 4c97cf5b1..c7f0be5bf 100644 --- a/vtkm/filter/contour/worklet/CMakeLists.txt +++ b/vtkm/filter/contour/worklet/CMakeLists.txt @@ -10,7 +10,8 @@ set(headers Clip.h - Contour.h + ContourFlyingEdges.h + ContourMarchingCells.h MIR.h ) diff --git a/vtkm/filter/contour/worklet/ContourFlyingEdges.h b/vtkm/filter/contour/worklet/ContourFlyingEdges.h new file mode 100644 index 000000000..85dabcf78 --- /dev/null +++ b/vtkm/filter/contour/worklet/ContourFlyingEdges.h @@ -0,0 +1,114 @@ +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_ContourFlyingEdges_h +#define vtk_m_worklet_ContourFlyingEdges_h + + +#include +#include +#include +#include + +namespace vtkm +{ +namespace worklet +{ + +/// \brief Compute the isosurface of a given \c CellSetStructured<3> input with +/// \c ArrayHandleUniformPointCoordinates for point coordinates using the Flying Edges algorithm. +class ContourFlyingEdges +{ +public: + //---------------------------------------------------------------------------- + ContourFlyingEdges(bool mergeDuplicates = true) + : SharedState(mergeDuplicates) + { + } + + //---------------------------------------------------------------------------- + vtkm::cont::ArrayHandle GetInterpolationEdgeIds() const + { + return this->SharedState.InterpolationEdgeIds; + } + + //---------------------------------------------------------------------------- + void SetMergeDuplicatePoints(bool merge) { this->SharedState.MergeDuplicatePoints = merge; } + + //---------------------------------------------------------------------------- + bool GetMergeDuplicatePoints() const { return this->SharedState.MergeDuplicatePoints; } + + //---------------------------------------------------------------------------- + vtkm::cont::ArrayHandle GetCellIdMap() const { return this->SharedState.CellIdMap; } + + //---------------------------------------------------------------------------- + template + void ProcessPointField(const InArrayType& input, const OutArrayType& output) const + { + + using vtkm::worklet::contour::MapPointField; + vtkm::worklet::DispatcherMapField applyFieldDispatcher; + + applyFieldDispatcher.Invoke(this->SharedState.InterpolationEdgeIds, + this->SharedState.InterpolationWeights, + input, + output); + } + + //---------------------------------------------------------------------------- + void ReleaseCellMapArrays() { this->SharedState.CellIdMap.ReleaseResources(); } + + // Filter called without normals generation + template + vtkm::cont::CellSetSingleType<> Run( + const std::vector& isovalues, + const vtkm::cont::CellSetStructured<3>& cells, + const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle, StorageTagVertices>& vertices) + { + this->SharedState.GenerateNormals = false; + vtkm::cont::ArrayHandle> normals; + + vtkm::cont::CellSetSingleType<> outputCells; + return flying_edges::execute( + cells, coordinateSystem, isovalues, input, vertices, normals, this->SharedState); + } + + // Filter called with normals generation + template + vtkm::cont::CellSetSingleType<> Run( + const std::vector& isovalues, + const vtkm::cont::CellSetStructured<3>& cells, + const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle, StorageTagVertices>& vertices, + vtkm::cont::ArrayHandle, StorageTagNormals>& normals) + { + this->SharedState.GenerateNormals = true; + vtkm::cont::CellSetSingleType<> outputCells; + return flying_edges::execute( + cells, coordinateSystem, isovalues, input, vertices, normals, this->SharedState); + } + +private: + vtkm::worklet::contour::CommonState SharedState; +}; +} +} // namespace vtkm::worklet + +#endif // vtk_m_worklet_ContourFlyingEdges_h diff --git a/vtkm/filter/contour/worklet/Contour.h b/vtkm/filter/contour/worklet/ContourMarchingCells.h similarity index 85% rename from vtkm/filter/contour/worklet/Contour.h rename to vtkm/filter/contour/worklet/ContourMarchingCells.h index e60189b63..1a27de5dc 100644 --- a/vtkm/filter/contour/worklet/Contour.h +++ b/vtkm/filter/contour/worklet/ContourMarchingCells.h @@ -8,16 +8,11 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ -#ifndef vtk_m_worklet_Contour_h -#define vtk_m_worklet_Contour_h - -#include -#include -#include +#ifndef vtk_m_worklet_ContourMarchingCells_h +#define vtk_m_worklet_ContourMarchingCells_h #include #include -#include #include @@ -25,8 +20,6 @@ namespace vtkm { namespace worklet { - - namespace contour { struct DeduceCoordType @@ -39,16 +32,6 @@ struct DeduceCoordType { result = marching_cells::execute(cells, coords, std::forward(args)...); } - - template - void operator()( - const vtkm::cont::ArrayHandle& coords, - const vtkm::cont::CellSetStructured<3>& cells, - vtkm::cont::CellSetSingleType<>& result, - Args&&... args) const - { - result = flying_edges::execute(cells, coords, std::forward(args)...); - } }; struct DeduceCellType @@ -62,13 +45,12 @@ struct DeduceCellType }; } -/// \brief Compute the isosurface of a given 3D data set, supports all -/// linear cell types -class Contour +/// \brief Compute the isosurface of a given 3D data set, supports all linear cell types +class ContourMarchingCells { public: //---------------------------------------------------------------------------- - Contour(bool mergeDuplicates = true) + ContourMarchingCells(bool mergeDuplicates = true) : SharedState(mergeDuplicates) { } @@ -89,6 +71,23 @@ public: vtkm::cont::ArrayHandle GetCellIdMap() const { return this->SharedState.CellIdMap; } //---------------------------------------------------------------------------- + template + void ProcessPointField(const InArrayType& input, const OutArrayType& output) const + { + + using vtkm::worklet::contour::MapPointField; + vtkm::worklet::DispatcherMapField applyFieldDispatcher; + + applyFieldDispatcher.Invoke(this->SharedState.InterpolationEdgeIds, + this->SharedState.InterpolationWeights, + input, + output); + } + + //---------------------------------------------------------------------------- + void ReleaseCellMapArrays() { this->SharedState.CellIdMap.ReleaseResources(); } + + // Filter called without normals generation template - void ProcessPointField(const InArrayType& input, const OutArrayType& output) const - { - - using vtkm::worklet::contour::MapPointField; - vtkm::worklet::DispatcherMapField applyFieldDispatcher; - - applyFieldDispatcher.Invoke(this->SharedState.InterpolationEdgeIds, - this->SharedState.InterpolationWeights, - input, - output); - } - - //---------------------------------------------------------------------------- - void ReleaseCellMapArrays() { this->SharedState.CellIdMap.ReleaseResources(); } private: vtkm::worklet::contour::CommonState SharedState; @@ -172,4 +155,4 @@ private: } } // namespace vtkm::worklet -#endif // vtk_m_worklet_Contour_h +#endif // vtk_m_worklet_ContourMarchingCells_h diff --git a/vtkm/filter/geometry_refinement/testing/UnitTestSplitSharpEdgesFilter.cxx b/vtkm/filter/geometry_refinement/testing/UnitTestSplitSharpEdgesFilter.cxx index 76f0e1ff2..260c46dc0 100644 --- a/vtkm/filter/geometry_refinement/testing/UnitTestSplitSharpEdgesFilter.cxx +++ b/vtkm/filter/geometry_refinement/testing/UnitTestSplitSharpEdgesFilter.cxx @@ -234,7 +234,7 @@ void TestWithStructuredData() contour.SetIsoValue(192); contour.SetMergeDuplicatePoints(true); contour.SetGenerateNormals(true); - contour.SetComputeFastNormalsForStructured(true); + contour.SetComputeFastNormals(true); contour.SetNormalArrayName("normals"); dataSet = contour.Execute(dataSet);