Add Shrink filter

This filter shrinks the cells of a DataSet towards their centroid, computed as the average position of the cell points

docs/changelog/2.0/shrink-filter

@@ -0,0 +1,3 @@
+# New Shrink filter
+
+The Shrink filter shrinks the cells of a DataSet towards their centroid, computed as the average position of the cell points. This filter disconnects the cells, duplicating the points connected to multiple cells. The resulting CellSet is always an `ExplicitCellSet`.

vtkm/filter/geometry_refinement/CMakeLists.txt

@@ -8,6 +8,7 @@
 ##  PURPOSE.  See the above copyright notice for more information.
 ##============================================================================
 set(geometry_refinement_headers
+  Shrink.h
   SplitSharpEdges.h
   Tetrahedralize.h
   Triangulate.h
@@ -16,6 +17,7 @@ set(geometry_refinement_headers
   )
 
 set(geometry_refinement_sources
+  Shrink.cxx
   SplitSharpEdges.cxx
   Tetrahedralize.cxx
   Triangulate.cxx

vtkm/filter/geometry_refinement/Shrink.cxx

@@ -0,0 +1,65 @@
+//============================================================================
+//  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/UncertainCellSet.h>
+#include <vtkm/filter/MapFieldPermutation.h>
+#include <vtkm/filter/geometry_refinement/Shrink.h>
+#include <vtkm/filter/geometry_refinement/worklet/Shrink.h>
+
+namespace
+{
+VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result,
+                          const vtkm::cont::Field& inputField,
+                          const vtkm::cont::ArrayHandle<vtkm::Id>& outputToInputCellMap)
+{
+  if (inputField.IsCellField() || inputField.IsWholeDataSetField())
+  {
+    result.AddField(inputField); // pass through
+    return true;
+  }
+  else if (inputField.IsPointField())
+  {
+    return vtkm::filter::MapFieldPermutation(inputField, outputToInputCellMap, result);
+  }
+  else
+  {
+    return false;
+  }
+}
+} // anonymous namespace
+
+namespace vtkm
+{
+namespace filter
+{
+namespace geometry_refinement
+{
+VTKM_CONT vtkm::cont::DataSet Shrink::DoExecute(const vtkm::cont::DataSet& input)
+{
+  const vtkm::cont::UnknownCellSet& inCellSet = input.GetCellSet();
+  const auto& oldCoords = input.GetCoordinateSystem().GetDataAsMultiplexer();
+
+  vtkm::cont::ArrayHandle<vtkm::Vec3f> newCoords;
+  vtkm::cont::ArrayHandle<vtkm::Id> oldPointsMapping;
+  vtkm::cont::CellSetExplicit<> newCellset;
+  vtkm::worklet::Shrink worklet;
+
+  worklet.Run(inCellSet, this->ShrinkFactor, oldCoords, newCoords, oldPointsMapping, newCellset);
+
+  auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, oldPointsMapping); };
+
+  vtkm::cont::CoordinateSystem activeCoordSystem = input.GetCoordinateSystem();
+  activeCoordSystem = vtkm::cont::CoordinateSystem(activeCoordSystem.GetName(), newCoords);
+
+  return this->CreateResultCoordinateSystem(input, newCellset, activeCoordSystem, mapper);
+}
+} // namespace geometry_refinement
+} // namespace filter
+} // namespace vtkm

vtkm/filter/geometry_refinement/Shrink.h

@@ -0,0 +1,48 @@
+//============================================================================
+//  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_geometry_refinement_Shrink_h
+#define vtk_m_filter_geometry_refinement_Shrink_h
+
+#include <vtkm/filter/FilterField.h>
+#include <vtkm/filter/geometry_refinement/vtkm_filter_geometry_refinement_export.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace geometry_refinement
+{
+/// \brief Shrink cells of an arbitrary dataset by a constant factor
+/// The Shrink filter shrinks the cells of a DataSet towards their centroid,
+/// computed as the average position of the cell points.
+/// This filter disconnects the cells, duplicating the points connected to multiple cells.
+/// The resulting CellSet is always an `ExplicitCellSet`.
+class VTKM_FILTER_GEOMETRY_REFINEMENT_EXPORT Shrink : public vtkm::filter::FilterField
+{
+public:
+  VTKM_CONT
+  void SetShrinkFactor(const vtkm::FloatDefault& factor)
+  {
+    this->ShrinkFactor = vtkm::Min(vtkm::Max(0, factor), 1); // Clamp shrink factor value
+  }
+
+  VTKM_CONT
+  const vtkm::FloatDefault& GetShrinkFactor() const { return this->ShrinkFactor; }
+
+private:
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
+  vtkm::FloatDefault ShrinkFactor = 0.5f;
+};
+} // namespace geometry_refinement
+} // namespace filter
+} // namespace vtkm
+
+#endif // vtk_m_filter_geometry_refinement_Shrink_h

vtkm/filter/geometry_refinement/testing/CMakeLists.txt

@@ -9,6 +9,7 @@
 ##============================================================================
 
 set(unit_tests
+  UnitTestShrinkFilter.cxx
   UnitTestSplitSharpEdgesFilter.cxx
   UnitTestTetrahedralizeFilter.cxx
   UnitTestTriangulateFilter.cxx

vtkm/filter/geometry_refinement/testing/UnitTestShrinkFilter.cxx

@@ -0,0 +1,169 @@
+//============================================================================
+//  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/filter/geometry_refinement/Shrink.h>
+
+#include <vtkm/cont/DataSetBuilderExplicit.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/testing/MakeTestDataSet.h>
+#include <vtkm/cont/testing/Testing.h>
+
+using vtkm::cont::testing::MakeTestDataSet;
+
+namespace
+{
+
+const std::array<vtkm::FloatDefault, 7> expectedPointVar{
+  { 10.1f, 20.1f, 30.2f, 30.2f, 20.1f, 40.2f, 50.3f }
+};
+
+const std::array<vtkm::Id, 7> expectedConnectivityArray{ { 0, 1, 2, 3, 4, 5, 6 } };
+
+const vtkm::Vec3f expectedCoords[7]{ { 0.333333f, 0.166666f, 0.0f }, { 0.833333f, 0.166666f, 0.0f },
+                                     { 0.833333f, 0.666666f, 0.0f }, { 1.25f, 1.0f, 0.0f },
+                                     { 1.25f, 0.5f, 0.0f },          { 1.75f, 1.0f, 0.0f },
+                                     { 1.75f, 1.5f, 0.0f } };
+
+const vtkm::FloatDefault expectedPointValueCube1[8]{ 10.1f, 20.1f, 50.2f,  40.1f,
+                                                     70.2f, 80.2f, 110.3f, 100.3f };
+const vtkm::Vec3f expectedCoordsCell1[8]{ { 0.4f, 0.4f, 0.4f }, { 0.6f, 0.4f, 0.4f },
+                                          { 0.6f, 0.6f, 0.4f }, { 0.4f, 0.6f, 0.4f },
+                                          { 0.4f, 0.4f, 0.6f }, { 0.6f, 0.4f, 0.6f },
+                                          { 0.6f, 0.6f, 0.6f }, { 0.4f, 0.6f, 0.6f } };
+
+
+void TestWithExplicitData()
+{
+  vtkm::cont::DataSet dataSet = MakeTestDataSet().Make3DExplicitDataSet0();
+
+  vtkm::filter::geometry_refinement::Shrink shrink;
+  shrink.SetFieldsToPass({ "pointvar", "cellvar" });
+
+  VTKM_TEST_ASSERT(test_equal(shrink.GetShrinkFactor(), 0.5f), "Wrong shrink factor default value");
+
+  // Test shrink factor clamping
+  shrink.SetShrinkFactor(1.5f);
+  VTKM_TEST_ASSERT(test_equal(shrink.GetShrinkFactor(), 1.0f), "Shrink factor not limited to 1");
+
+  shrink.SetShrinkFactor(-0.5f);
+  VTKM_TEST_ASSERT(test_equal(shrink.GetShrinkFactor(), 0.0f),
+                   "Shrink factor is not always positive");
+
+  shrink.SetShrinkFactor(0.5f);
+
+  vtkm::cont::DataSet output = shrink.Execute(dataSet);
+  VTKM_TEST_ASSERT(test_equal(output.GetNumberOfCells(), dataSet.GetNumberOfCells()),
+                   "Wrong number of cells for Shrink filter");
+  VTKM_TEST_ASSERT(test_equal(output.GetNumberOfPoints(), 7), "Wrong number of points for Shrink");
+
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> outCellData =
+    output.GetField("cellvar").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>();
+
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(0), 100.1f), "Wrong cell field data");
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(1), 100.2f), "Wrong cell field data");
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> outPointData =
+    output.GetField("pointvar").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>();
+
+  for (vtkm::IdComponent i = 0; i < outPointData.GetNumberOfValues(); i++)
+  {
+    VTKM_TEST_ASSERT(test_equal(outPointData.ReadPortal().Get(i), expectedPointVar[i]),
+                     "Wrong point field data");
+  }
+
+  const auto& connectivityArray =
+    output.GetCellSet().AsCellSet<vtkm::cont::CellSetExplicit<>>().GetConnectivityArray(
+      vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
+  auto connectivityArrayPortal = connectivityArray.ReadPortal();
+
+  for (vtkm::IdComponent i = 0; i < connectivityArray.GetNumberOfValues(); i++)
+  {
+    VTKM_TEST_ASSERT(test_equal(connectivityArrayPortal.Get(i), expectedConnectivityArray[i]),
+                     "Wrong connectivity array value");
+  }
+
+  auto newCoords = output.GetCoordinateSystem().GetDataAsMultiplexer();
+  auto newCoordsP = newCoords.ReadPortal();
+
+  for (vtkm::IdComponent i = 0; i < newCoords.GetNumberOfValues(); i++)
+  {
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoords[i][0]),
+                     "Wrong point coordinates");
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoords[i][1]),
+                     "Wrong point coordinates");
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoords[i][2]),
+                     "Wrong point coordinates");
+  }
+}
+
+
+void TestWithUniformData()
+{
+  vtkm::cont::DataSet dataSet = MakeTestDataSet().Make3DUniformDataSet0();
+
+  vtkm::filter::geometry_refinement::Shrink shrink;
+  shrink.SetFieldsToPass({ "pointvar", "cellvar" });
+
+  shrink.SetShrinkFactor(0.2f);
+
+  vtkm::cont::DataSet output = shrink.Execute(dataSet);
+
+  VTKM_TEST_ASSERT(test_equal(output.GetNumberOfCells(), dataSet.GetNumberOfCells()),
+                   "Number of cells changed after filtering");
+  VTKM_TEST_ASSERT(test_equal(output.GetNumberOfPoints(), 4 * 8), "Wrong number of points");
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> outCellData =
+    output.GetField("cellvar").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>();
+
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(0), 100.1), "Wrong cell field data");
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(1), 100.2), "Wrong cell field data");
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(2), 100.3), "Wrong cell field data");
+  VTKM_TEST_ASSERT(test_equal(outCellData.ReadPortal().Get(3), 100.4), "Wrong cell field data");
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> outPointData =
+    output.GetField("pointvar").GetData().AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>();
+
+  for (vtkm::IdComponent i = 0; i < 8; i++) // Test for the first cell only
+  {
+    VTKM_TEST_ASSERT(test_equal(outPointData.ReadPortal().Get(i), expectedPointValueCube1[i]),
+                     "Wrong cell field data");
+  }
+
+  auto newCoords = output.GetCoordinateSystem().GetDataAsMultiplexer();
+  auto newCoordsP = newCoords.ReadPortal();
+
+  for (vtkm::IdComponent i = 0; i < 8; i++)
+  {
+    std::cout << newCoordsP.Get(i)[0] << " " << expectedCoordsCell1[i][0] << std::endl;
+    std::cout << newCoordsP.Get(i)[1] << " " << expectedCoordsCell1[i][1] << std::endl;
+    std::cout << newCoordsP.Get(i)[2] << " " << expectedCoordsCell1[i][2] << std::endl;
+
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoordsCell1[i][0]),
+                     "Wrong point coordinates");
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoordsCell1[i][1]),
+                     "Wrong point coordinates");
+    VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoordsCell1[i][2]),
+                     "Wrong point coordinates");
+  }
+}
+
+
+void TestShrinkFilter()
+{
+  TestWithExplicitData();
+  TestWithUniformData();
+}
+
+} // anonymous namespace
+
+int UnitTestShrinkFilter(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestShrinkFilter, argc, argv);
+}

vtkm/filter/geometry_refinement/worklet/CMakeLists.txt

@@ -9,6 +9,7 @@
 ##============================================================================
 
 set(headers
+  Shrink.h
   SplitSharpEdges.h
   Tetrahedralize.h
   Triangulate.h

vtkm/filter/geometry_refinement/worklet/Shrink.h

@@ -0,0 +1,150 @@
+//============================================================================
+//  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_Shrink_h
+#define vtk_m_worklet_Shrink_h
+
+#include <vtkm/worklet/CellDeepCopy.h>
+#include <vtkm/worklet/ScatterCounting.h>
+
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/Invoker.h>
+#include <vtkm/exec/ParametricCoordinates.h>
+
+
+namespace vtkm
+{
+namespace worklet
+{
+class Shrink
+{
+public:
+  struct PrepareCellsForShrink : vtkm::worklet::WorkletVisitCellsWithPoints
+  {
+    using ControlSignature = void(CellSetIn,
+                                  FieldOutCell numPoints,
+                                  FieldOutCell centroids,
+                                  FieldOutCell shapes,
+                                  FieldInPoint coords);
+    using ExecutionSignature =
+      void(PointCount, _2 numPoints, _3 centroids, _4 shapes, _5 coords, CellShape);
+
+    using InputDomain = _1;
+
+    template <typename CoordsArrayType, typename ShapeIdType, typename ShapeTagType>
+    VTKM_EXEC void operator()(vtkm::IdComponent numPointsInCell,
+                              vtkm::IdComponent& numPoints,
+                              vtkm::Vec3f& centroids,
+                              ShapeIdType& shapes,
+                              const CoordsArrayType& coords,
+                              ShapeTagType cellShape) const
+    {
+      numPoints = numPointsInCell;
+      shapes = cellShape.Id;
+
+      vtkm::Vec3f cellCenter;
+      vtkm::exec::ParametricCoordinatesCenter(numPoints, cellShape, cellCenter);
+      vtkm::exec::CellInterpolate(coords, cellCenter, cellShape, centroids);
+    }
+  };
+
+  struct ComputeNewPoints : vtkm::worklet::WorkletVisitCellsWithPoints
+  {
+    ComputeNewPoints(vtkm::FloatDefault shrinkFactor)
+      : ShrinkFactor(shrinkFactor)
+    {
+    }
+    using ControlSignature = void(CellSetIn,
+                                  FieldInCell offsets,
+                                  FieldInCell centroids,
+                                  FieldOutCell oldPointsMapping,
+                                  FieldOutCell newPoints,
+                                  FieldOutCell newCoords,
+                                  FieldInPoint coords);
+    using ExecutionSignature = void(_2 offsets,
+                                    _3 centroids,
+                                    _4 oldPointsMapping,
+                                    _5 newPoints,
+                                    _6 newCoords,
+                                    _7 coords,
+                                    VisitIndex localPointNum,
+                                    PointIndices globalPointIndex);
+    using InputDomain = _1;
+
+    using ScatterType = vtkm::worklet::ScatterCounting;
+
+    template <typename PointIndicesVecType, typename CoordsArrayTypeIn, typename CoordsArrayTypeOut>
+    VTKM_EXEC void operator()(const vtkm::Id& offsets,
+                              const vtkm::Vec3f& centroids,
+                              vtkm::Id& oldPointsMapping,
+                              vtkm::Id& newPoints,
+                              CoordsArrayTypeOut& newCoords,
+                              const CoordsArrayTypeIn& coords,
+                              vtkm::IdComponent localPtIndex,
+                              const PointIndicesVecType& globalPointIndex) const
+    {
+      newPoints = offsets + localPtIndex;
+      oldPointsMapping = globalPointIndex[localPtIndex];
+      newCoords = centroids + this->ShrinkFactor * (coords[localPtIndex] - centroids);
+    }
+
+  private:
+    vtkm::FloatDefault ShrinkFactor;
+  };
+
+  template <typename CellSetType,
+            typename CoordsComType,
+            typename CoordsInStorageType,
+            typename CoordsOutStorageType,
+            typename OldPointsMappingType,
+            typename NewCellSetType>
+  void Run(
+    const CellSetType& oldCellset,
+    const vtkm::FloatDefault shinkFactor,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsInStorageType>& oldCoords,
+    vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsOutStorageType>& newCoords,
+    vtkm::cont::ArrayHandle<vtkm::Id, OldPointsMappingType>& oldPointsMapping,
+    NewCellSetType& newCellset)
+  {
+    vtkm::cont::Invoker invoke;
+
+    // First pass : count the new number of points per cell, shapes and compute centroids
+    vtkm::cont::ArrayHandle<vtkm::IdComponent> CellPointCount;
+    vtkm::cont::ArrayHandle<vtkm::Vec3f> Centroids;
+    vtkm::cont::CellSetExplicit<>::ShapesArrayType shapeArray;
+    invoke(PrepareCellsForShrink{}, oldCellset, CellPointCount, Centroids, shapeArray, oldCoords);
+
+
+    // Second pass : compute new point positions and mappings to input points
+    vtkm::cont::ArrayHandle<vtkm::Id> newPoints;
+    vtkm::worklet::ScatterCounting scatter(CellPointCount, true);
+    vtkm::cont::ArrayHandle<vtkm::Id> offsets = scatter.GetInputToOutputMap();
+    vtkm::Id totalPoints = scatter.GetOutputRange(CellPointCount.GetNumberOfValues());
+
+    ComputeNewPoints worklet = ComputeNewPoints(shinkFactor);
+    invoke(worklet,
+           scatter,
+           oldCellset,
+           offsets,
+           Centroids,
+           oldPointsMapping,
+           newPoints,
+           newCoords,
+           oldCoords);
+
+    newCellset.Fill(totalPoints,
+                    shapeArray,
+                    newPoints,
+                    vtkm::cont::ConvertNumComponentsToOffsets(CellPointCount));
+  }
+};
+} // namespace vtkm::worklet
+} // namespace vtkm
+
+#endif // vtk_m_worklet_Shrink_h