Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'compute-moments-any-vec-size'

Browse Source 
4064a3bd0 Allow ComputeMoments to operate on any scalar field

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !3065
This commit is contained in:
Kenneth Moreland 2023-06-05 18:42:03 +00:00 committed by Kitware Robot
commit 7488ad731d
3 changed files with 65 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
docs/changelog/compute-moments-any-vec-size.md Normal file
View File

@ -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.

15
vtkm/filter/image_processing/ComputeMoments.cxx
View File

@ -21,18 +21,6 @@ namespace filter
{
namespace image_processing
{
using SupportedTypes = vtkm::List<vtkm::Float32,
vtkm::Float64,
vtkm::Vec<vtkm::Float32, 2>,
vtkm::Vec<vtkm::Float64, 2>,
vtkm::Vec<vtkm::Float32, 3>,
vtkm::Vec<vtkm::Float64, 3>,
vtkm::Vec<vtkm::Float32, 4>,
vtkm::Vec<vtkm::Float64, 4>,
vtkm::Vec<vtkm::Float32, 6>,
vtkm::Vec<vtkm::Float64, 6>,
vtkm::Vec<vtkm::Float32, 9>,
vtkm::Vec<vtkm::Float64, 9>>;
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<SupportedTypes, VTKM_DEFAULT_STORAGE_LIST>(
resolveType);
this->CastAndCallVariableVecField(field, resolveType);
return output;
}

79
vtkm/filter/image_processing/worklet/ComputeMoments.h
View File

@ -15,6 +15,9 @@
#include <vtkm/Math.h>
#include <vtkm/worklet/WorkletPointNeighborhood.h>
#include <vtkm/cont/ArrayHandleRecombineVec.h>
#include <vtkm/cont/ArrayHandleRuntimeVec.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/UncertainArrayHandle.h>
#include <vtkm/cont/UncertainCellSet.h>
@ -54,13 +57,19 @@ public:
using ExecutionSignature = void(_2, Boundary, _3);
template <typename NeighIn, typename T>
template <typename NeighIn, typename TOut>
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<T>::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<ComponentType>::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<T>(vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q) * image.Get(i, j, 0));
ComponentType multiplier =
static_cast<ComponentType>(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<ComponentType>(this->SpacingProduct);
}
}
private:
@ -126,13 +145,19 @@ public:
using ExecutionSignature = void(_2, Boundary, _3);
template <typename NeighIn, typename T>
template <typename NeighIn, typename TOut>
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<T>::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<ComponentType>::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<T>(vtkm::Pow(radius[0], p) * vtkm::Pow(radius[1], q) *
vtkm::Pow(radius[2], r) * image.Get(i, j, k));
ComponentType multiplier = static_cast<ComponentType>(
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<ComponentType>(this->SpacingProduct);
}
}
private:
@ -195,9 +230,9 @@ public:
class ResolveUnknownCellSet
{
public:
template <typename T, typename S>
template <typename T>
void operator()(const vtkm::cont::CellSetStructured<2>& input,
const vtkm::cont::ArrayHandle<T, S>& pixels,
const vtkm::cont::ArrayHandleRecombineVec<T>& pixels,
vtkm::Vec3f spacing,
vtkm::Float64 radius,
int maxOrder,
@ -212,7 +247,7 @@ public:
{
const int q = order - p;
vtkm::cont::ArrayHandle<T> moments;
vtkm::cont::ArrayHandleRuntimeVec<T> moments{ pixels.GetNumberOfComponents() };
DispatcherType dispatcher(WorkletType{ spacing, radius, p, q });
dispatcher.Invoke(input, pixels, moments);
@ -226,9 +261,9 @@ public:
}
}
template <typename T, typename S>
template <typename T>
void operator()(const vtkm::cont::CellSetStructured<3>& input,
const vtkm::cont::ArrayHandle<T, S>& pixels,
const vtkm::cont::ArrayHandleRecombineVec<T>& pixels,
vtkm::Vec3f spacing,
vtkm::Float64 radius,
int maxOrder,
@ -246,7 +281,7 @@ public:
{
const int p = order - r - q;
vtkm::cont::ArrayHandle<T> moments;
vtkm::cont::ArrayHandleRuntimeVec<T> moments{ pixels.GetNumberOfComponents() };
DispatcherType dispatcher(WorkletType{ spacing, radius, p, q, r });
dispatcher.Invoke(input, pixels, moments);
@ -263,9 +298,9 @@ public:
}
};
template <typename T, typename S>
template <typename T>
void Run(const vtkm::cont::UnknownCellSet& input,
const vtkm::cont::ArrayHandle<T, S>& pixels,
const vtkm::cont::ArrayHandleRecombineVec<T>& pixels,
int maxOrder,
vtkm::cont::DataSet& output) const
{