diff --git a/docs/changelog/add-slice-multiple-filter.md b/docs/changelog/add-slice-multiple-filter.md new file mode 100644 index 000000000..d11dab22e --- /dev/null +++ b/docs/changelog/add-slice-multiple-filter.md @@ -0,0 +1,3 @@ +# Adding SliceMultiple filter + +The SliceMultiple filter can accept multiple implicit functions and output a merged dataset. This filter is added to the filter/contour. The code of this filter is adapted from vtk-h. The mechanism that merges multiple datasets in this filter is supposed to support more general datasets and work a separate filter in the future. diff --git a/vtkm/filter/contour/CMakeLists.txt b/vtkm/filter/contour/CMakeLists.txt index 9f45e5cb9..c7a4105b8 100644 --- a/vtkm/filter/contour/CMakeLists.txt +++ b/vtkm/filter/contour/CMakeLists.txt @@ -17,6 +17,7 @@ set(contour_headers ContourMarchingCells.h MIRFilter.h Slice.h + SliceMultiple.h ) set(contour_sources_device @@ -26,6 +27,7 @@ set(contour_sources_device ContourMarchingCells.cxx MIRFilter.cxx Slice.cxx + SliceMultiple.cxx ) set(contour_sources diff --git a/vtkm/filter/contour/SliceMultiple.cxx b/vtkm/filter/contour/SliceMultiple.cxx new file mode 100644 index 000000000..9be415797 --- /dev/null +++ b/vtkm/filter/contour/SliceMultiple.cxx @@ -0,0 +1,176 @@ +//============================================================================ +// 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 +{ +namespace contour +{ +class OffsetWorklet : public vtkm::worklet::WorkletMapField +{ +protected: + vtkm::Id OffsetValue; + +public: + VTKM_CONT + OffsetWorklet(const vtkm::Id offset) + : OffsetValue(offset) + { + } + typedef void ControlSignature(FieldInOut); + typedef void ExecutionSignature(_1); + VTKM_EXEC void operator()(vtkm::Id& value) const { value += this->OffsetValue; } +}; +// Original MergeContours code come from +// https://github.com/Alpine-DAV/ascent/blob/develop/src/libs/vtkh/filters/Slice.cpp#L517 +class MergeContours +{ + vtkm::cont::PartitionedDataSet DataSets; + +public: + MergeContours(vtkm::cont::PartitionedDataSet& dataSets) + : DataSets(dataSets) + { + } + vtkm::cont::DataSet MergeDataSets() + { + vtkm::cont::DataSet res; + vtkm::Id numOfDataSet = this->DataSets.GetNumberOfPartitions(); + vtkm::Id numCells = 0; + vtkm::Id numPoints = 0; + std::vector cellOffsets(numOfDataSet); + std::vector pointOffsets(numOfDataSet); + for (vtkm::Id i = 0; i < numOfDataSet; i++) + { + auto cellSet = this->DataSets.GetPartition(i).GetCellSet(); + //We assume all cells are triangles here + if (!cellSet.IsType>()) + { + throw vtkm::cont::ErrorFilterExecution( + "Expected singletype cell set as the result of contour."); + } + cellOffsets[i] = numCells; + numCells += cellSet.GetNumberOfCells(); + pointOffsets[i] = numPoints; + numPoints += this->DataSets.GetPartition(i).GetNumberOfPoints(); + } + const vtkm::Id connSize = numCells * 3; + // Calculating merged offsets for all domains + vtkm::cont::ArrayHandle conn; + conn.Allocate(connSize); + for (vtkm::Id i = 0; i < numOfDataSet; i++) + { + auto cellSet = this->DataSets.GetPartition(i).GetCellSet(); + if (!cellSet.IsType>()) + { + throw vtkm::cont::ErrorFilterExecution( + "Expected singletype cell set as the result of contour."); + } + // Grabing the connectivity and copy it into the larger array + vtkm::cont::CellSetSingleType<> singleType = + cellSet.AsCellSet>(); + const vtkm::cont::ArrayHandle connPerDataSet = singleType.GetConnectivityArray( + vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()); + vtkm::Id copySize = connPerDataSet.GetNumberOfValues(); + vtkm::Id start = 0; + vtkm::cont::Algorithm::CopySubRange( + connPerDataSet, start, copySize, conn, cellOffsets[i] * 3); + // We offset the connectiviy we just copied in so we references the + // correct points + if (cellOffsets[i] != 0) + { + vtkm::cont::Invoker invoker; + invoker(OffsetWorklet{ pointOffsets[i] }, + vtkm::cont::make_ArrayHandleView(conn, cellOffsets[i] * 3, copySize)); + } + } + vtkm::cont::CellSetSingleType<> cellSet; + cellSet.Fill(numPoints, vtkm::CELL_SHAPE_TRIANGLE, 3, conn); + res.SetCellSet(cellSet); + // Merging selected fields and coordinates + vtkm::IdComponent numFields = this->DataSets.GetPartition(0).GetNumberOfFields(); + for (vtkm::IdComponent i = 0; i < numFields; i++) + { + const vtkm::cont::Field& field = this->DataSets.GetPartition(0).GetField(i); + bool isSupported = (field.GetAssociation() == vtkm::cont::Field::Association::Points || + field.GetAssociation() == vtkm::cont::Field::Association::Cells); + if (!isSupported) + { + continue; + } + vtkm::cont::UnknownArrayHandle outFieldArray = field.GetData().NewInstanceBasic(); + bool assocPoints = field.GetAssociation() == vtkm::cont::Field::Association::Points; + if (assocPoints) + { + outFieldArray.Allocate(numPoints); + } + else + { + outFieldArray.Allocate(numCells); + } + auto resolveType = [&](auto& concreteOut) { + vtkm::Id offset = 0; + for (vtkm::Id partitionId = 0; partitionId < this->DataSets.GetNumberOfPartitions(); + ++partitionId) + { + vtkm::cont::UnknownArrayHandle in = this->DataSets.GetPartition(partitionId) + .GetField(field.GetName(), field.GetAssociation()) + .GetData(); + vtkm::Id copySize = in.GetNumberOfValues(); + auto viewOut = vtkm::cont::make_ArrayHandleView(concreteOut, offset, copySize); + vtkm::cont::ArrayCopy(in, viewOut); + offset += copySize; + } + }; + outFieldArray.CastAndCallWithExtractedArray(resolveType); + res.AddField(vtkm::cont::Field(field.GetName(), field.GetAssociation(), outFieldArray)); + } + //Labeling fields that belong to the coordinate system. + //There might be multiple coordinates systems, assuming all partitions have the same name of the coordinates system + vtkm::IdComponent numCoordsNames = + this->DataSets.GetPartition(0).GetNumberOfCoordinateSystems(); + for (vtkm::IdComponent i = 0; i < numCoordsNames; i++) + { + res.AddCoordinateSystem(this->DataSets.GetPartition(0).GetCoordinateSystemName(i)); + } + return res; + } +}; +vtkm::cont::DataSet SliceMultiple::DoExecute(const vtkm::cont::DataSet& input) +{ + vtkm::cont::PartitionedDataSet slices; + //Executing Slice filter several times and merge results together + for (vtkm::IdComponent i = 0; i < static_cast(this->FunctionList.size()); i++) + { + vtkm::filter::contour::Slice slice; + slice.SetImplicitFunction(this->GetImplicitFunction(i)); + slice.SetFieldsToPass(this->GetFieldsToPass()); + auto result = slice.Execute(input); + slices.AppendPartition(result); + } + if (slices.GetNumberOfPartitions() > 1) + { + //Since the slice filter have already selected fields + //the mergeCountours will copy all existing fields + MergeContours merger(slices); + vtkm::cont::DataSet mergedResults = merger.MergeDataSets(); + return mergedResults; + } + return slices.GetPartition(0); +} +} // namespace contour +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/contour/SliceMultiple.h b/vtkm/filter/contour/SliceMultiple.h new file mode 100644 index 000000000..03a680543 --- /dev/null +++ b/vtkm/filter/contour/SliceMultiple.h @@ -0,0 +1,52 @@ +//============================================================================ +// 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_Slice_Multi_h +#define vtk_m_filter_contour_Slice_Multi_h + +#include +#include + +#include + +namespace vtkm +{ +namespace filter +{ +namespace contour +{ +/// \brief This filter can accept multiple implicit functions used by the slice filter. +/// It returns a merged data set that contains multiple results returned by the slice filter. +class VTKM_FILTER_CONTOUR_EXPORT SliceMultiple : public vtkm::filter::contour::Contour +{ +public: + /// Set/Get the implicit function that is used to perform the slicing. + /// + VTKM_CONT + void AddImplicitFunction(const vtkm::ImplicitFunctionGeneral& func) + { + FunctionList.push_back(func); + } + VTKM_CONT + const vtkm::ImplicitFunctionGeneral& GetImplicitFunction(vtkm::Id id) const + { + VTKM_ASSERT(id < static_cast(FunctionList.size())); + return FunctionList[id]; + } + +private: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override; + + std::vector FunctionList; +}; +} // namespace contour +} // namespace filter +} // namespace vtkm + +#endif // vtk_m_filter_contour_Slice_Multi_h diff --git a/vtkm/filter/contour/testing/CMakeLists.txt b/vtkm/filter/contour/testing/CMakeLists.txt index 843140bed..617de0469 100644 --- a/vtkm/filter/contour/testing/CMakeLists.txt +++ b/vtkm/filter/contour/testing/CMakeLists.txt @@ -16,6 +16,7 @@ set(unit_tests set(unit_tests_device UnitTestContourFilter.cxx # Algorithm used, needs device compiler UnitTestMIRFilter.cxx # Algorithm used, needs device compiler + UnitTestSliceMultipleFilter.cxx # Algorithm used, needs device compiler ) if (VTKm_ENABLE_RENDERING) diff --git a/vtkm/filter/contour/testing/UnitTestSliceMultipleFilter.cxx b/vtkm/filter/contour/testing/UnitTestSliceMultipleFilter.cxx new file mode 100644 index 000000000..07c166737 --- /dev/null +++ b/vtkm/filter/contour/testing/UnitTestSliceMultipleFilter.cxx @@ -0,0 +1,134 @@ +//============================================================================ +// 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 +#include +#include +namespace +{ +class SetPointValuesWorklet : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldIn, FieldOut, FieldOut, FieldOut); + using ExecutionSignature = void(_1, _2, _3, _4); + template + VTKM_EXEC void operator()(const CoordinatesType& coordinates, + ScalarType& scalar, + V3Type& vec3, + V4Type& vec4) const + { + scalar = + static_cast((coordinates[2] * 3 * 3 + coordinates[1] * 3 + coordinates[0]) * 0.1); + vec3 = { coordinates[0] * 0.1, coordinates[1] * 0.1, coordinates[2] * 0.1 }; + vec4 = { + coordinates[0] * 0.1, coordinates[1] * 0.1, coordinates[2] * 0.1, coordinates[0] * 0.1 + }; + return; + } +}; + +class SetCellValuesWorklet : public vtkm::worklet::WorkletVisitCellsWithPoints +{ +public: + using ControlSignature = void(CellSetIn, FieldInPoint, FieldOutCell, FieldOutCell, FieldOutCell); + using ExecutionSignature = void(_2, _3, _4, _5); + using InputDomain = _1; + template + VTKM_EXEC void operator()(const PointFieldVecType& pointFieldVec, + ScalarType& scalar, + V3Type& vec3, + V4Type& vec4) const + { + //pointFieldVec has 8 values + scalar = static_cast(pointFieldVec[0]); + vec3 = { pointFieldVec[0] * 0.1, pointFieldVec[1] * 0.1, pointFieldVec[2] * 0.1 }; + vec4 = { + pointFieldVec[0] * 0.1, pointFieldVec[1] * 0.1, pointFieldVec[2] * 0.1, pointFieldVec[3] * 0.1 + }; + return; + } +}; + +vtkm::cont::DataSet MakeTestDatasetStructured3D() +{ + static constexpr vtkm::Id xdim = 3, ydim = 3, zdim = 3; + static const vtkm::Id3 dim(xdim, ydim, zdim); + vtkm::cont::DataSet ds; + ds = vtkm::cont::DataSetBuilderUniform::Create( + dim, vtkm::Vec3f(-1.0f, -1.0f, -1.0f), vtkm::Vec3f(1, 1, 1)); + vtkm::cont::ArrayHandle pointScalars; + vtkm::cont::ArrayHandle pointV3; + //a customized vector which is not in vtkm::TypeListCommon{} + vtkm::cont::ArrayHandle> pointV4; + pointScalars.Allocate(xdim * ydim * zdim); + pointV3.Allocate(xdim * ydim * zdim); + pointV4.Allocate(xdim * ydim * zdim); + vtkm::cont::Invoker invoker; + invoker( + SetPointValuesWorklet{}, ds.GetCoordinateSystem().GetData(), pointScalars, pointV3, pointV4); + ds.AddPointField("pointScalars", pointScalars); + ds.AddPointField("pointV3", pointV3); + ds.AddPointField("pointV4", pointV4); + //adding cell data + vtkm::cont::ArrayHandle cellScalars; + vtkm::cont::ArrayHandle cellV3; + vtkm::cont::ArrayHandle> cellV4; + vtkm::Id NumCells = ds.GetNumberOfCells(); + cellScalars.Allocate(NumCells); + cellV3.Allocate(NumCells); + cellV4.Allocate(NumCells); + invoker(SetCellValuesWorklet{}, ds.GetCellSet(), pointScalars, cellScalars, cellV3, cellV4); + ds.AddCellField("cellScalars", cellScalars); + ds.AddCellField("cellV3", cellV3); + ds.AddCellField("cellV4", cellV4); + return ds; +} + +void TestSliceMultipleFilter() +{ + auto ds = MakeTestDatasetStructured3D(); + vtkm::Plane plane1({ 0, 0, 0 }, { 0, 0, 1 }); + vtkm::Plane plane2({ 0, 0, 0 }, { 0, 1, 0 }); + vtkm::Plane plane3({ 0, 0, 0 }, { 1, 0, 0 }); + vtkm::filter::contour::SliceMultiple sliceMultiple; + sliceMultiple.AddImplicitFunction(plane1); + sliceMultiple.AddImplicitFunction(plane2); + sliceMultiple.AddImplicitFunction(plane3); + auto result = sliceMultiple.Execute(ds); + VTKM_TEST_ASSERT(result.GetNumberOfPoints() == 27, "wrong number of points in merged data set"); + VTKM_TEST_ASSERT(result.GetCoordinateSystem().GetData().GetNumberOfValues() == 27, + "wrong number of scalars in merged data set"); + vtkm::cont::ArrayHandle CheckingScalars; + vtkm::cont::ArrayHandle CheckingV3; + vtkm::cont::ArrayHandle> CheckingV4; + vtkm::cont::Invoker invoker; + invoker(SetPointValuesWorklet{}, + result.GetCoordinateSystem().GetData(), + CheckingScalars, + CheckingV3, + CheckingV4); + VTKM_TEST_ASSERT( + test_equal_ArrayHandles(CheckingScalars, result.GetField("pointScalars").GetData()), + "wrong scalar values"); + VTKM_TEST_ASSERT(test_equal_ArrayHandles(CheckingV3, result.GetField("pointV3").GetData()), + "wrong pointV3 values"); + VTKM_TEST_ASSERT(test_equal_ArrayHandles(CheckingV4, result.GetField("pointV4").GetData()), + "wrong pointV4 values"); + VTKM_TEST_ASSERT(result.GetNumberOfCells() == 24, "wrong number of cells in merged data set"); +} +} // anonymous namespace +int UnitTestSliceMultipleFilter(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestSliceMultipleFilter, argc, argv); +} diff --git a/vtkm/filter/contour/vtkm.module b/vtkm/filter/contour/vtkm.module index 56fe5b405..e1fbd164a 100644 --- a/vtkm/filter/contour/vtkm.module +++ b/vtkm/filter/contour/vtkm.module @@ -7,6 +7,7 @@ DEPENDS vtkm_filter_core vtkm_filter_vector_analysis vtkm_filter_mesh_info + vtkm_filter_multi_block PRIVATE_DEPENDS vtkm_worklet TEST_DEPENDS