From 67d6780fba29b0070da16ea685895adc212ba6a9 Mon Sep 17 00:00:00 2001 From: Louis Gombert Date: Thu, 2 Mar 2023 22:34:59 +0100 Subject: [PATCH] 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 | 3 + .../filter/geometry_refinement/CMakeLists.txt | 2 + vtkm/filter/geometry_refinement/Shrink.cxx | 65 +++++++ vtkm/filter/geometry_refinement/Shrink.h | 48 +++++ .../testing/CMakeLists.txt | 1 + .../testing/UnitTestShrinkFilter.cxx | 169 ++++++++++++++++++ .../worklet/CMakeLists.txt | 1 + .../geometry_refinement/worklet/Shrink.h | 150 ++++++++++++++++ 8 files changed, 439 insertions(+) create mode 100644 docs/changelog/2.0/shrink-filter create mode 100644 vtkm/filter/geometry_refinement/Shrink.cxx create mode 100644 vtkm/filter/geometry_refinement/Shrink.h create mode 100644 vtkm/filter/geometry_refinement/testing/UnitTestShrinkFilter.cxx create mode 100644 vtkm/filter/geometry_refinement/worklet/Shrink.h diff --git a/docs/changelog/2.0/shrink-filter b/docs/changelog/2.0/shrink-filter new file mode 100644 index 000000000..36e523c31 --- /dev/null +++ b/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`. diff --git a/vtkm/filter/geometry_refinement/CMakeLists.txt b/vtkm/filter/geometry_refinement/CMakeLists.txt index ff52259c9..b146ff5c3 100644 --- a/vtkm/filter/geometry_refinement/CMakeLists.txt +++ b/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 diff --git a/vtkm/filter/geometry_refinement/Shrink.cxx b/vtkm/filter/geometry_refinement/Shrink.cxx new file mode 100644 index 000000000..ed59e53b4 --- /dev/null +++ b/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 +#include +#include +#include + +namespace +{ +VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result, + const vtkm::cont::Field& inputField, + const vtkm::cont::ArrayHandle& 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 newCoords; + vtkm::cont::ArrayHandle 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 diff --git a/vtkm/filter/geometry_refinement/Shrink.h b/vtkm/filter/geometry_refinement/Shrink.h new file mode 100644 index 000000000..173865f1c --- /dev/null +++ b/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 +#include + +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 diff --git a/vtkm/filter/geometry_refinement/testing/CMakeLists.txt b/vtkm/filter/geometry_refinement/testing/CMakeLists.txt index 98f75c805..3e563d003 100644 --- a/vtkm/filter/geometry_refinement/testing/CMakeLists.txt +++ b/vtkm/filter/geometry_refinement/testing/CMakeLists.txt @@ -9,6 +9,7 @@ ##============================================================================ set(unit_tests + UnitTestShrinkFilter.cxx UnitTestSplitSharpEdgesFilter.cxx UnitTestTetrahedralizeFilter.cxx UnitTestTriangulateFilter.cxx diff --git a/vtkm/filter/geometry_refinement/testing/UnitTestShrinkFilter.cxx b/vtkm/filter/geometry_refinement/testing/UnitTestShrinkFilter.cxx new file mode 100644 index 000000000..ff88c443d --- /dev/null +++ b/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 + +#include +#include +#include +#include + +using vtkm::cont::testing::MakeTestDataSet; + +namespace +{ + +const std::array expectedPointVar{ + { 10.1f, 20.1f, 30.2f, 30.2f, 20.1f, 40.2f, 50.3f } +}; + +const std::array 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 outCellData = + output.GetField("cellvar").GetData().AsArrayHandle>(); + + 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 outPointData = + output.GetField("pointvar").GetData().AsArrayHandle>(); + + 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>().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 outCellData = + output.GetField("cellvar").GetData().AsArrayHandle>(); + + 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 outPointData = + output.GetField("pointvar").GetData().AsArrayHandle>(); + + 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); +} diff --git a/vtkm/filter/geometry_refinement/worklet/CMakeLists.txt b/vtkm/filter/geometry_refinement/worklet/CMakeLists.txt index 7124cb6f2..a81485a76 100644 --- a/vtkm/filter/geometry_refinement/worklet/CMakeLists.txt +++ b/vtkm/filter/geometry_refinement/worklet/CMakeLists.txt @@ -9,6 +9,7 @@ ##============================================================================ set(headers + Shrink.h SplitSharpEdges.h Tetrahedralize.h Triangulate.h diff --git a/vtkm/filter/geometry_refinement/worklet/Shrink.h b/vtkm/filter/geometry_refinement/worklet/Shrink.h new file mode 100644 index 000000000..5d2a6c5bb --- /dev/null +++ b/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 +#include + +#include +#include +#include + + +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 + 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 + 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 + void Run( + const CellSetType& oldCellset, + const vtkm::FloatDefault shinkFactor, + const vtkm::cont::ArrayHandle, CoordsInStorageType>& oldCoords, + vtkm::cont::ArrayHandle, CoordsOutStorageType>& newCoords, + vtkm::cont::ArrayHandle& oldPointsMapping, + NewCellSetType& newCellset) + { + vtkm::cont::Invoker invoke; + + // First pass : count the new number of points per cell, shapes and compute centroids + vtkm::cont::ArrayHandle CellPointCount; + vtkm::cont::ArrayHandle 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 newPoints; + vtkm::worklet::ScatterCounting scatter(CellPointCount, true); + vtkm::cont::ArrayHandle 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