diff --git a/data/baseline/filter/moments.png b/data/baseline/filter/moments.png new file mode 100644 index 000000000..63ae037f7 --- /dev/null +++ b/data/baseline/filter/moments.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c077570d0acc5ee81707e1b763e5edbf285af840474d48ab5b051599cb809697 +size 81863 diff --git a/data/baseline/filter/moments0.png b/data/baseline/filter/moments0.png new file mode 100644 index 000000000..2c0f1ff2a --- /dev/null +++ b/data/baseline/filter/moments0.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d4979a7c0e46943bd7a2c94581023c7db12a30ee12bda735b0041c64cff52196 +size 65438 diff --git a/data/baseline/filter/moments12.png b/data/baseline/filter/moments12.png new file mode 100644 index 000000000..37e3f4c07 --- /dev/null +++ b/data/baseline/filter/moments12.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cbd7292ca4006866afc2339e39e14a3edf8c5983e04a5ceea80b5d698293fc86 +size 94505 diff --git a/docs/changelog/compute-moments-any-vec-size.md b/docs/changelog/compute-moments-any-vec-size.md new file mode 100644 index 000000000..75efc6c9f --- /dev/null +++ b/docs/changelog/compute-moments-any-vec-size.md @@ -0,0 +1,7 @@ +# ComputeMoments filter now operates on any scalar field + +Previously, the `ComputeMoments` filter only operated on a finite set of +array types as its input field. This included a prescribed list of `Vec` +sizes for the input. The filter has been updated to use more generic +interfaces to the field's array (and float fallback) to enable the +computation of moments on any type of scalar field. diff --git a/vtkm/cont/testing/UnitTestDataSetRectilinear.cxx b/vtkm/cont/testing/UnitTestDataSetRectilinear.cxx index f1aa27f25..a953e57f6 100644 --- a/vtkm/cont/testing/UnitTestDataSetRectilinear.cxx +++ b/vtkm/cont/testing/UnitTestDataSetRectilinear.cxx @@ -92,7 +92,7 @@ static void TwoDimRectilinearTest() vtkm::Id cells[2][4] = { { 0, 1, 4, 3 }, { 1, 2, 5, 4 } }; for (vtkm::Id cellIndex = 0; cellIndex < 2; cellIndex++) { - vtkm::Id4 pointIds = pointToCell.GetIndices(pointToCell.FlatToLogicalToIndex(cellIndex)); + vtkm::Id4 pointIds = pointToCell.GetIndices(pointToCell.FlatToLogicalVisitIndex(cellIndex)); for (vtkm::IdComponent localPointIndex = 0; localPointIndex < 4; localPointIndex++) { VTKM_TEST_ASSERT(pointIds[localPointIndex] == cells[cellIndex][localPointIndex], @@ -106,7 +106,7 @@ static void TwoDimRectilinearTest() for (vtkm::Id pointIndex = 0; pointIndex < 6; pointIndex++) { vtkm::VecVariable retrievedCellIds = - cellToPoint.GetIndices(cellToPoint.FlatToLogicalToIndex(pointIndex)); + cellToPoint.GetIndices(cellToPoint.FlatToLogicalVisitIndex(pointIndex)); VTKM_TEST_ASSERT(retrievedCellIds.GetNumberOfComponents() <= 4, "Got wrong number of cell ids."); for (vtkm::IdComponent cellIndex = 0; cellIndex < retrievedCellIds.GetNumberOfComponents(); diff --git a/vtkm/cont/testing/UnitTestDataSetUniform.cxx b/vtkm/cont/testing/UnitTestDataSetUniform.cxx index ff400c015..3f0c57929 100644 --- a/vtkm/cont/testing/UnitTestDataSetUniform.cxx +++ b/vtkm/cont/testing/UnitTestDataSetUniform.cxx @@ -94,7 +94,7 @@ static void TwoDimUniformTest() vtkm::Id cells[2][4] = { { 0, 1, 4, 3 }, { 1, 2, 5, 4 } }; for (vtkm::Id cellIndex = 0; cellIndex < 2; cellIndex++) { - vtkm::Id4 pointIds = pointToCell.GetIndices(pointToCell.FlatToLogicalToIndex(cellIndex)); + vtkm::Id4 pointIds = pointToCell.GetIndices(pointToCell.FlatToLogicalVisitIndex(cellIndex)); for (vtkm::IdComponent localPointIndex = 0; localPointIndex < 4; localPointIndex++) { VTKM_TEST_ASSERT(pointIds[localPointIndex] == cells[cellIndex][localPointIndex], @@ -108,7 +108,7 @@ static void TwoDimUniformTest() for (vtkm::Id pointIndex = 0; pointIndex < 6; pointIndex++) { vtkm::VecVariable retrievedCellIds = - cellToPoint.GetIndices(cellToPoint.FlatToLogicalToIndex(pointIndex)); + cellToPoint.GetIndices(cellToPoint.FlatToLogicalVisitIndex(pointIndex)); VTKM_TEST_ASSERT(retrievedCellIds.GetNumberOfComponents() <= 4, "Got wrong number of cell ids."); for (vtkm::IdComponent cellIndex = 0; cellIndex < retrievedCellIds.GetNumberOfComponents(); diff --git a/vtkm/exec/ConnectivityStructured.h b/vtkm/exec/ConnectivityStructured.h index 0f14923c3..fd3fce942 100644 --- a/vtkm/exec/ConnectivityStructured.h +++ b/vtkm/exec/ConnectivityStructured.h @@ -11,6 +11,7 @@ #ifndef vtk_m_exec_ConnectivityStructured_h #define vtk_m_exec_ConnectivityStructured_h +#include #include #include #include @@ -77,28 +78,54 @@ public: return Helper::GetIndices(this->Internals, index); } + VTKM_EXEC_CONT SchedulingRangeType FlatToLogicalVisitIndex(vtkm::Id flatVisitIndex) const + { + return Helper::FlatToLogicalVisitIndex(this->Internals, flatVisitIndex); + } + + VTKM_EXEC_CONT SchedulingRangeType FlatToLogicalIncidentIndex(vtkm::Id flatIncidentIndex) const + { + return Helper::FlatToLogicalIncidentIndex(this->Internals, flatIncidentIndex); + } + + VTKM_EXEC_CONT vtkm::Id LogicalToFlatVisitIndex( + const SchedulingRangeType& logicalVisitIndex) const + { + return Helper::LogicalToFlatVisitIndex(this->Internals, logicalVisitIndex); + } + + VTKM_EXEC_CONT vtkm::Id LogicalToFlatIncidentIndex( + const SchedulingRangeType& logicalIncidentIndex) const + { + return Helper::LogicalToFlatIncidentIndex(this->Internals, logicalIncidentIndex); + } + VTKM_EXEC_CONT + VTKM_DEPRECATED(2.1, "Use FlatToLogicalIncidentIndex.") SchedulingRangeType FlatToLogicalFromIndex(vtkm::Id flatFromIndex) const { - return Helper::FlatToLogicalFromIndex(this->Internals, flatFromIndex); + return this->FlatToLogicalIncidentIndex(flatFromIndex); } VTKM_EXEC_CONT + VTKM_DEPRECATED(2.1, "Use LogicalToFlatIncidentIndex.") vtkm::Id LogicalToFlatFromIndex(const SchedulingRangeType& logicalFromIndex) const { - return Helper::LogicalToFlatFromIndex(this->Internals, logicalFromIndex); + return this->LogicalToFlatIncidentIndex(logicalFromIndex); } VTKM_EXEC_CONT + VTKM_DEPRECATED(2.1, "Use FlatToLogicalVisitIndex.") SchedulingRangeType FlatToLogicalToIndex(vtkm::Id flatToIndex) const { - return Helper::FlatToLogicalToIndex(this->Internals, flatToIndex); + return this->FlatToLogicalVisitIndex(flatToIndex); } VTKM_EXEC_CONT + VTKM_DEPRECATED(2.1, "Use LogicalToFlatVisitIndex.") vtkm::Id LogicalToFlatToIndex(const SchedulingRangeType& logicalToIndex) const { - return Helper::LogicalToFlatToIndex(this->Internals, logicalToIndex); + return this->LogicalToFlatVisitIndex(logicalToIndex); } VTKM_EXEC_CONT diff --git a/vtkm/exec/arg/ThreadIndicesPointNeighborhood.h b/vtkm/exec/arg/ThreadIndicesPointNeighborhood.h index 00318e25a..4bcf7f915 100644 --- a/vtkm/exec/arg/ThreadIndicesPointNeighborhood.h +++ b/vtkm/exec/arg/ThreadIndicesPointNeighborhood.h @@ -72,7 +72,7 @@ public: inputIndex, visitIndex, outputIndex, - vtkm::exec::BoundaryState{ detail::To3D(connectivity.FlatToLogicalToIndex(inputIndex)), + vtkm::exec::BoundaryState{ detail::To3D(connectivity.FlatToLogicalVisitIndex(inputIndex)), detail::To3D(connectivity.GetPointDimensions()) }) { } diff --git a/vtkm/exec/arg/ThreadIndicesTopologyMap.h b/vtkm/exec/arg/ThreadIndicesTopologyMap.h index 95f2c9165..970afa0fa 100644 --- a/vtkm/exec/arg/ThreadIndicesTopologyMap.h +++ b/vtkm/exec/arg/ThreadIndicesTopologyMap.h @@ -181,7 +181,7 @@ public: this->InputIndex = inputIndex; this->VisitIndex = visitIndex; this->OutputIndex = outputIndex; - this->LogicalIndex = connectivity.FlatToLogicalToIndex(this->InputIndex); + this->LogicalIndex = connectivity.FlatToLogicalVisitIndex(this->InputIndex); this->IndicesIncident = connectivity.GetIndices(this->LogicalIndex); this->CellShape = connectivity.GetCellShape(this->InputIndex); } @@ -338,7 +338,7 @@ public: const ConnectivityType& connectivity) { this->ThreadIndex = threadIndex; - this->LogicalIndex = connectivity.FlatToLogicalToIndex(inputIndex); + this->LogicalIndex = connectivity.FlatToLogicalVisitIndex(inputIndex); this->IndicesIncident = connectivity.GetIndices(this->LogicalIndex); this->CellShape = connectivity.GetCellShape(inputIndex); } @@ -503,7 +503,7 @@ public: this->OutputIndex = outputIndex; const vtkm::Id permutedIndex = permutation.Portal.Get(this->InputIndex); - this->LogicalIndex = permutation.Connectivity.FlatToLogicalToIndex(permutedIndex); + this->LogicalIndex = permutation.Connectivity.FlatToLogicalVisitIndex(permutedIndex); this->IndicesIncident = permutation.Connectivity.GetIndices(this->LogicalIndex); this->CellShape = permutation.Connectivity.GetCellShape(permutedIndex); } diff --git a/vtkm/filter/image_processing/ComputeMoments.cxx b/vtkm/filter/image_processing/ComputeMoments.cxx index 78499739c..96f350694 100644 --- a/vtkm/filter/image_processing/ComputeMoments.cxx +++ b/vtkm/filter/image_processing/ComputeMoments.cxx @@ -21,18 +21,6 @@ namespace filter { namespace image_processing { -using SupportedTypes = vtkm::List, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec, - vtkm::Vec>; VTKM_CONT ComputeMoments::ComputeMoments() { @@ -53,8 +41,7 @@ VTKM_CONT vtkm::cont::DataSet ComputeMoments::DoExecute(const vtkm::cont::DataSe auto resolveType = [&](const auto& concrete) { worklet.Run(input.GetCellSet(), concrete, this->Order, output); }; - field.GetData().CastAndCallForTypesWithFloatFallback( - resolveType); + this->CastAndCallVariableVecField(field, resolveType); return output; } diff --git a/vtkm/filter/image_processing/testing/CMakeLists.txt b/vtkm/filter/image_processing/testing/CMakeLists.txt index e159e9ddc..5ac223e45 100644 --- a/vtkm/filter/image_processing/testing/CMakeLists.txt +++ b/vtkm/filter/image_processing/testing/CMakeLists.txt @@ -9,13 +9,17 @@ ##============================================================================ set(unit_tests + RenderTestComputeMoments.cxx UnitTestImageDifferenceFilter.cxx UnitTestImageMedianFilter.cxx ) set(libraries vtkm_filter_image_processing - vtkm_source) + vtkm_source + vtkm_rendering + vtkm_rendering_testing +) vtkm_unit_tests( SOURCES ${unit_tests} diff --git a/vtkm/filter/image_processing/testing/RenderTestComputeMoments.cxx b/vtkm/filter/image_processing/testing/RenderTestComputeMoments.cxx new file mode 100644 index 000000000..9514069c2 --- /dev/null +++ b/vtkm/filter/image_processing/testing/RenderTestComputeMoments.cxx @@ -0,0 +1,44 @@ +//============================================================================ +// 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 +{ + +void TestComputeMoments() +{ + vtkm::source::Wavelet source; + vtkm::cont::DataSet data = source.Execute(); + + vtkm::filter::image_processing::ComputeMoments filter; + filter.SetActiveField("RTData"); + filter.SetOrder(2); + filter.SetRadius(2); + vtkm::cont::DataSet result = filter.Execute(data); + + vtkm::rendering::testing::RenderTestOptions testOptions; + testOptions.AllowedPixelErrorRatio = 0.001f; + testOptions.ColorTable = vtkm::cont::ColorTable("inferno"); + testOptions.EnableAnnotations = false; + vtkm::rendering::testing::RenderTest(result, "index", "filter/moments.png", testOptions); + vtkm::rendering::testing::RenderTest(result, "index0", "filter/moments0.png", testOptions); + vtkm::rendering::testing::RenderTest(result, "index12", "filter/moments12.png", testOptions); +} +} // namespace + +int RenderTestComputeMoments(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestComputeMoments, argc, argv); +} diff --git a/vtkm/filter/image_processing/vtkm.module b/vtkm/filter/image_processing/vtkm.module index 5fdd1257c..129750e9e 100644 --- a/vtkm/filter/image_processing/vtkm.module +++ b/vtkm/filter/image_processing/vtkm.module @@ -9,3 +9,5 @@ PRIVATE_DEPENDS vtkm_worklet TEST_DEPENDS vtkm_source + vtkm_rendering + vtkm_rendering_testing diff --git a/vtkm/filter/image_processing/worklet/ComputeMoments.h b/vtkm/filter/image_processing/worklet/ComputeMoments.h index d99722496..283989abe 100644 --- a/vtkm/filter/image_processing/worklet/ComputeMoments.h +++ b/vtkm/filter/image_processing/worklet/ComputeMoments.h @@ -15,6 +15,9 @@ #include #include +#include +#include +#include #include #include #include @@ -54,13 +57,19 @@ public: using ExecutionSignature = void(_2, Boundary, _3); - template + template VTKM_EXEC void operator()(const NeighIn& image, const vtkm::exec::BoundaryState& boundary, - T& moment) const + TOut& moment) const { - // TODO: type safety and numerical precision - auto sum = vtkm::TypeTraits::ZeroInitialization(); + using ComponentType = typename TOut::ComponentType; + const vtkm::IdComponent numComponents = moment.GetNumberOfComponents(); + + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] = vtkm::TypeTraits::ZeroInitialization(); + } // Clamp the radius to the dataset bounds (discard out-of-bounds points). const auto minRadius = boundary.ClampNeighborIndex(-this->RadiusDiscrete); @@ -85,13 +94,23 @@ public: if (vtkm::Dot(radius, radius) <= 1) { - sum += - static_cast(vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q) * image.Get(i, j, 0)); + ComponentType multiplier = + static_cast(vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q)); + auto inputField = image.Get(i, j, 0); + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] += multiplier * inputField[componentI]; + } } } } - moment = T(sum * this->SpacingProduct); + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] *= static_cast(this->SpacingProduct); + } } private: @@ -126,13 +145,19 @@ public: using ExecutionSignature = void(_2, Boundary, _3); - template + template VTKM_EXEC void operator()(const NeighIn& image, const vtkm::exec::BoundaryState& boundary, - T& moment) const + TOut& moment) const { - // TODO: type safety and numerical precision - auto sum = vtkm::TypeTraits::ZeroInitialization(); + using ComponentType = typename TOut::ComponentType; + const vtkm::IdComponent numComponents = moment.GetNumberOfComponents(); + + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] = vtkm::TypeTraits::ZeroInitialization(); + } // Clamp the radius to the dataset bounds (discard out-of-bounds points). const auto minRadius = boundary.ClampNeighborIndex(-this->RadiusDiscrete); @@ -165,14 +190,24 @@ public: if (vtkm::Dot(radius, radius) <= 1) { - sum += static_cast(vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q) * - vtkm::Pow(radius[2], r) * image.Get(i, j, k)); + ComponentType multiplier = static_cast( + vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q) * vtkm::Pow(radius[2], r)); + auto inputField = image.Get(i, j, k); + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] += multiplier * inputField[componentI]; + } } } } } - moment = T(sum * this->SpacingProduct); + // For variable sized Vecs, need to iterate over each component. + for (vtkm::IdComponent componentI = 0; componentI < numComponents; ++componentI) + { + moment[componentI] *= static_cast(this->SpacingProduct); + } } private: @@ -195,9 +230,9 @@ public: class ResolveUnknownCellSet { public: - template + template void operator()(const vtkm::cont::CellSetStructured<2>& input, - const vtkm::cont::ArrayHandle& pixels, + const vtkm::cont::ArrayHandleRecombineVec& pixels, vtkm::Vec3f spacing, vtkm::Float64 radius, int maxOrder, @@ -212,7 +247,7 @@ public: { const int q = order - p; - vtkm::cont::ArrayHandle moments; + vtkm::cont::ArrayHandleRuntimeVec moments{ pixels.GetNumberOfComponents() }; DispatcherType dispatcher(WorkletType{ spacing, radius, p, q }); dispatcher.Invoke(input, pixels, moments); @@ -226,9 +261,9 @@ public: } } - template + template void operator()(const vtkm::cont::CellSetStructured<3>& input, - const vtkm::cont::ArrayHandle& pixels, + const vtkm::cont::ArrayHandleRecombineVec& pixels, vtkm::Vec3f spacing, vtkm::Float64 radius, int maxOrder, @@ -246,7 +281,7 @@ public: { const int p = order - r - q; - vtkm::cont::ArrayHandle moments; + vtkm::cont::ArrayHandleRuntimeVec moments{ pixels.GetNumberOfComponents() }; DispatcherType dispatcher(WorkletType{ spacing, radius, p, q, r }); dispatcher.Invoke(input, pixels, moments); @@ -263,9 +298,9 @@ public: } }; - template + template void Run(const vtkm::cont::UnknownCellSet& input, - const vtkm::cont::ArrayHandle& pixels, + const vtkm::cont::ArrayHandleRecombineVec& pixels, int maxOrder, vtkm::cont::DataSet& output) const { diff --git a/vtkm/internal/ConnectivityStructuredInternals.h b/vtkm/internal/ConnectivityStructuredInternals.h index 7e2e0717a..448faf6d6 100644 --- a/vtkm/internal/ConnectivityStructuredInternals.h +++ b/vtkm/internal/ConnectivityStructuredInternals.h @@ -607,29 +607,29 @@ struct ConnectivityStructuredIndexHelper cellIndices = Connectivity.GetIndices(cellId); // Look up the offset into the face list for each cell type diff --git a/vtkm/rendering/raytracing/VolumeRendererStructured.cxx b/vtkm/rendering/raytracing/VolumeRendererStructured.cxx index 86d084d2c..9fff4ef35 100644 --- a/vtkm/rendering/raytracing/VolumeRendererStructured.cxx +++ b/vtkm/rendering/raytracing/VolumeRendererStructured.cxx @@ -61,7 +61,7 @@ public: vtkm::Id cellId{}; auto self = static_cast(this); self->Locator.FindCell(point, cellId, parametric); - cell = self->Conn.FlatToLogicalToIndex(cellId); + cell = self->Conn.FlatToLogicalVisitIndex(cellId); self->ComputeInvSpacing(cell, point, invSpacing, parametric); } @@ -74,7 +74,7 @@ public: VTKM_EXEC inline vtkm::Id GetCellIndex(const vtkm::Id3& cell) const { - return static_cast(this)->Conn.LogicalToFlatToIndex(cell); + return static_cast(this)->Conn.LogicalToFlatVisitIndex(cell); } VTKM_EXEC @@ -88,7 +88,7 @@ public: inline void GetMinPoint(const vtkm::Id3& cell, vtkm::Vec3f_32& point) const { const vtkm::Id pointIndex = - static_cast(this)->Conn.LogicalToFlatFromIndex(cell); + static_cast(this)->Conn.LogicalToFlatIncidentIndex(cell); point = static_cast(this)->Coordinates.Get(pointIndex); } };