add multislice filter

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.

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

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 <vtkm/cont/ArrayCopyDevice.h>
+#include <vtkm/cont/ArrayHandleTransform.h>
+#include <vtkm/filter/contour/Slice.h>
+#include <vtkm/filter/contour/SliceMultiple.h>
+#include <vtkm/filter/multi_block/AmrArrays.h>
+#include <vtkm/worklet/WorkletMapField.h>
+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<vtkm::Id> cellOffsets(numOfDataSet);
+    std::vector<vtkm::Id> 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<vtkm::cont::CellSetSingleType<>>())
+      {
+        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<vtkm::Id> conn;
+    conn.Allocate(connSize);
+    for (vtkm::Id i = 0; i < numOfDataSet; i++)
+    {
+      auto cellSet = this->DataSets.GetPartition(i).GetCellSet();
+      if (!cellSet.IsType<vtkm::cont::CellSetSingleType<>>())
+      {
+        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<vtkm::cont::CellSetSingleType<>>();
+      const vtkm::cont::ArrayHandle<vtkm::Id> 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<vtkm::IdComponent>(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

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 <vtkm/filter/contour/Contour.h>
+#include <vtkm/filter/contour/vtkm_filter_contour_export.h>
+
+#include <vtkm/ImplicitFunction.h>
+
+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<vtkm::Id>(FunctionList.size()));
+    return FunctionList[id];
+  }
+
+private:
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
+
+  std::vector<vtkm::ImplicitFunctionGeneral> FunctionList;
+};
+} // namespace contour
+} // namespace filter
+} // namespace vtkm
+
+#endif // vtk_m_filter_contour_Slice_Multi_h

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)

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 <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/Invoker.h>
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/filter/contour/SliceMultiple.h>
+#include <vtkm/io/VTKDataSetWriter.h>
+#include <vtkm/worklet/WorkletMapField.h>
+#include <vtkm/worklet/WorkletMapTopology.h>
+namespace
+{
+class SetPointValuesWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn, FieldOut, FieldOut, FieldOut);
+  using ExecutionSignature = void(_1, _2, _3, _4);
+  template <typename CoordinatesType, typename ScalarType, typename V3Type, typename V4Type>
+  VTKM_EXEC void operator()(const CoordinatesType& coordinates,
+                            ScalarType& scalar,
+                            V3Type& vec3,
+                            V4Type& vec4) const
+  {
+    scalar =
+      static_cast<ScalarType>((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 <typename PointFieldVecType, typename ScalarType, typename V3Type, typename V4Type>
+  VTKM_EXEC void operator()(const PointFieldVecType& pointFieldVec,
+                            ScalarType& scalar,
+                            V3Type& vec3,
+                            V4Type& vec4) const
+  {
+    //pointFieldVec has 8 values
+    scalar = static_cast<ScalarType>(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<vtkm::Float64> pointScalars;
+  vtkm::cont::ArrayHandle<vtkm::Vec3f_64> pointV3;
+  //a customized vector which is not in vtkm::TypeListCommon{}
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 4>> 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<vtkm::Float64> cellScalars;
+  vtkm::cont::ArrayHandle<vtkm::Vec3f_64> cellV3;
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 4>> 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<vtkm::Float64> CheckingScalars;
+  vtkm::cont::ArrayHandle<vtkm::Vec3f_64> CheckingV3;
+  vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 4>> 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);
+}

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