Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Update clip filter's field map to work on any field type

Browse Source 
The previous implementation of the map field in the clip filters
(`ClipWithField` and `ClipWithImplicitFunction`) checked for common field
types and interpolated those. If the field value type did not match, it
would either convert the field to floats (which is at odds with what VTK
does) or fail outright if the `Vec` length is not supported.

The map field function for clip has been changed to support all possible
types. It does this by using the extract component functionality to get
data from any type of array.
This commit is contained in:
Kenneth Moreland 2023-01-13 09:45:25 -07:00
parent 0513eaf5fc
commit 1889447d82
4 changed files with 109 additions and 185 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
docs/changelog/clip-interp-any-field.md Normal file
View File

@ -0,0 +1,11 @@
# Update clip filter's field map to work on any field type
The previous implementation of the map field in the clip filters
(`ClipWithField` and `ClipWithImplicitFunction`) checked for common field
types and interpolated those. If the field value type did not match, it
would either convert the field to floats (which is at odds with what VTK
does) or fail outright if the `Vec` length is not supported.
The map field function for clip has been changed to support all possible
types. It does this by using the extract component functionality to get
data from any type of array.

42
vtkm/filter/contour/ClipWithField.cxx
View File

@ -23,19 +23,6 @@ namespace contour
{
namespace
{
struct ClipWithFieldProcessCoords
{
template <typename T, typename Storage>
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<T, Storage>& inCoords,
const std::string& coordsName,
const vtkm::worklet::Clip& worklet,
vtkm::cont::DataSet& output) const
{
vtkm::cont::ArrayHandle<T> outArray = worklet.ProcessPointField(inCoords);
vtkm::cont::CoordinateSystem outCoords(coordsName, outArray);
output.AddCoordinateSystem(outCoords);
}
};
bool DoMapField(vtkm::cont::DataSet& result,
const vtkm::cont::Field& field,
@ -43,17 +30,18 @@ bool DoMapField(vtkm::cont::DataSet& result,
{
if (field.IsPointField())
{
auto resolve = [&](const auto& concrete) {
vtkm::cont::UnknownArrayHandle inputArray = field.GetData();
vtkm::cont::UnknownArrayHandle outputArray = inputArray.NewInstanceBasic();
auto resolve = [&](const auto& concreteIn) {
// use std::decay to remove const ref from the decltype of concrete.
using T = typename std::decay_t<decltype(concrete)>::ValueType;
vtkm::cont::ArrayHandle<T> outputArray;
outputArray = worklet.ProcessPointField(concrete);
result.AddPointField(field.GetName(), outputArray);
using BaseT = typename std::decay_t<decltype(concreteIn)>::ValueType::ComponentType;
auto concreteOut = outputArray.ExtractArrayFromComponents<BaseT>();
worklet.ProcessPointField(concreteIn, concreteOut);
};
field.GetData()
.CastAndCallForTypesWithFloatFallback<VTKM_DEFAULT_TYPE_LIST, VTKM_DEFAULT_STORAGE_LIST>(
resolve);
inputArray.CastAndCallWithExtractedArray(resolve);
result.AddPointField(field.GetName(), outputArray);
return true;
}
else if (field.IsCellField())
@ -94,17 +82,7 @@ vtkm::cont::DataSet ClipWithField::DoExecute(const vtkm::cont::DataSet& input)
this->CastAndCallScalarField(this->GetFieldFromDataSet(input).GetData(), resolveFieldType);
auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); };
vtkm::cont::DataSet output = this->CreateResult(input, outputCellSet, mapper);
// Compute the new boundary points and add them to the output:
for (vtkm::IdComponent coordSystemId = 0; coordSystemId < input.GetNumberOfCoordinateSystems();
++coordSystemId)
{
const vtkm::cont::CoordinateSystem& coords = input.GetCoordinateSystem(coordSystemId);
coords.GetData().CastAndCall(ClipWithFieldProcessCoords{}, coords.GetName(), worklet, output);
}
return output;
return this->CreateResult(input, outputCellSet, mapper);
}
} // namespace contour
} // namespace filter

33
vtkm/filter/contour/ClipWithImplicitFunction.cxx
View File

@ -23,23 +23,9 @@ namespace contour
namespace
{
struct ClipWithImplicitFunctionProcessCoords
{
template <typename T, typename Storage>
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<T, Storage>& inCoords,
const std::string& coordsName,
const vtkm::worklet::Clip& worklet,
vtkm::cont::DataSet& output) const
{
vtkm::cont::ArrayHandle<T> outArray = worklet.ProcessPointField(inCoords);
vtkm::cont::CoordinateSystem outCoords(coordsName, outArray);
output.AddCoordinateSystem(outCoords);
}
};
bool DoMapField(vtkm::cont::DataSet& result,
const vtkm::cont::Field& field,
const vtkm::worklet::Clip& Worklet)
vtkm::worklet::Clip& worklet)
{
if (field.IsPointField())
{
@ -47,7 +33,7 @@ bool DoMapField(vtkm::cont::DataSet& result,
// use std::decay to remove const ref from the decltype of concrete.
using T = typename std::decay_t<decltype(concrete)>::ValueType;
vtkm::cont::ArrayHandle<T> outputArray;
outputArray = Worklet.ProcessPointField(concrete);
worklet.ProcessPointField(concrete, outputArray);
result.AddPointField(field.GetName(), outputArray);
};
@ -59,7 +45,7 @@ bool DoMapField(vtkm::cont::DataSet& result,
else if (field.IsCellField())
{
// Use the precompiled field permutation function.
vtkm::cont::ArrayHandle<vtkm::Id> permutation = Worklet.GetCellMapOutputToInput();
vtkm::cont::ArrayHandle<vtkm::Id> permutation = worklet.GetCellMapOutputToInput();
return vtkm::filter::MapFieldPermutation(field, permutation, result);
}
else if (field.IsWholeDataSetField())
@ -87,18 +73,7 @@ vtkm::cont::DataSet ClipWithImplicitFunction::DoExecute(const vtkm::cont::DataSe
worklet.Run(inputCellSet, this->Function, this->Offset, inputCoords, this->Invert);
auto mapper = [&](auto& result, const auto& f) { DoMapField(result, f, worklet); };
vtkm::cont::DataSet output = this->CreateResult(input, outputCellSet, mapper);
// compute output coordinates
for (vtkm::IdComponent coordSystemId = 0; coordSystemId < input.GetNumberOfCoordinateSystems();
++coordSystemId)
{
const vtkm::cont::CoordinateSystem& coords = input.GetCoordinateSystem(coordSystemId);
coords.GetData().CastAndCall(
ClipWithImplicitFunctionProcessCoords{}, coords.GetName(), worklet, output);
}
return output;
return this->CreateResult(input, outputCellSet, mapper);
}
} // namespace contour
} // namespace filter

208
vtkm/filter/contour/worklet/Clip.h
View File

@ -645,6 +645,7 @@ public:
this->InCellInterpolationKeys,
this->InCellInterpolationInfo,
this->CellMapOutputToInput);
this->InterpolationKeysBuilt = false;
// Get unique EdgeInterpolation : unique edge points.
// LowerBound for edgeInterpolation : get index into new edge points array.
@ -763,140 +764,97 @@ public:
return this->Run(cellSet, clipFunction, 0.0, coords, invert);
}
template <typename ArrayHandleType>
class InterpolateField
struct PerformEdgeInterpolations : public vtkm::worklet::WorkletMapField
{
public:
using ValueType = typename ArrayHandleType::ValueType;
using TypeMappedValue = vtkm::List<ValueType>;
using ControlSignature = void(FieldIn edgeInterpolations,
WholeArrayIn originalField,
FieldOut outputField);
using ExecutionSignature = void(_1, _2, _3);
InterpolateField(vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolationArray,
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationKeys,
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationInfo,
vtkm::Id edgePointsOffset,
vtkm::Id inCellPointsOffset,
ArrayHandleType* output)
: EdgeInterpolationArray(edgeInterpolationArray)
, InCellInterpolationKeys(inCellInterpolationKeys)
, InCellInterpolationInfo(inCellInterpolationInfo)
, EdgePointsOffset(edgePointsOffset)
, InCellPointsOffset(inCellPointsOffset)
, Output(output)
template <typename FieldPortal, typename T>
VTKM_EXEC void operator()(const EdgeInterpolation& edgeInterp,
const FieldPortal& originalField,
T& output) const
{
T v1 = originalField.Get(edgeInterp.Vertex1);
T v2 = originalField.Get(edgeInterp.Vertex2);
// Interpolate per-vertex because some vec-like objects do not allow intermediate variables
using VTraits = vtkm::VecTraits<T>;
using CType = typename VTraits::ComponentType;
VTKM_ASSERT(VTraits::GetNumberOfComponents(v1) == VTraits::GetNumberOfComponents(output));
VTKM_ASSERT(VTraits::GetNumberOfComponents(v2) == VTraits::GetNumberOfComponents(output));
for (vtkm::IdComponent component = 0; component < VTraits::GetNumberOfComponents(output);
++component)
{
CType c1 = VTraits::GetComponent(v1, component);
CType c2 = VTraits::GetComponent(v2, component);
CType o = static_cast<CType>(((c1 - c2) * edgeInterp.Weight) + c1);
VTraits::SetComponent(output, component, o);
}
}
class PerformEdgeInterpolations : public vtkm::worklet::WorkletMapField
{
public:
PerformEdgeInterpolations(vtkm::Id edgePointsOffset)
: EdgePointsOffset(edgePointsOffset)
{
}
using ControlSignature = void(FieldIn edgeInterpolations, WholeArrayInOut outputField);
using ExecutionSignature = void(_1, _2, WorkIndex);
template <typename EdgeInterp, typename OutputFieldPortal>
VTKM_EXEC void operator()(const EdgeInterp& ei,
OutputFieldPortal& field,
vtkm::Id workIndex) const
{
using T = typename OutputFieldPortal::ValueType;
T v1 = field.Get(ei.Vertex1);
T v2 = field.Get(ei.Vertex2);
field.Set(this->EdgePointsOffset + workIndex,
static_cast<T>(internal::Scale(T(v1 - v2), ei.Weight) + v1));
}
private:
vtkm::Id EdgePointsOffset;
};
class PerformInCellInterpolations : public vtkm::worklet::WorkletReduceByKey
{
public:
using ControlSignature = void(KeysIn keys, ValuesIn toReduce, ReducedValuesOut centroid);
using ExecutionSignature = void(_2, _3);
template <typename MappedValueVecType, typename MappedValueType>
VTKM_EXEC void operator()(const MappedValueVecType& toReduce, MappedValueType& centroid) const
{
vtkm::IdComponent numValues = toReduce.GetNumberOfComponents();
MappedValueType sum = toReduce[0];
for (vtkm::IdComponent i = 1; i < numValues; i++)
{
MappedValueType value = toReduce[i];
// static_cast is for when MappedValueType is a small int that gets promoted to int32.
sum = static_cast<MappedValueType>(sum + value);
}
centroid = internal::Scale(sum, 1. / static_cast<vtkm::Float64>(numValues));
}
};
template <typename Storage>
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<ValueType, Storage>& field) const
{
vtkm::worklet::Keys<vtkm::Id> interpolationKeys(InCellInterpolationKeys);
vtkm::Id numberOfOriginalValues = field.GetNumberOfValues();
vtkm::Id numberOfEdgePoints = EdgeInterpolationArray.GetNumberOfValues();
vtkm::Id numberOfInCellPoints = interpolationKeys.GetUniqueKeys().GetNumberOfValues();
ArrayHandleType result;
result.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
vtkm::cont::Algorithm::CopySubRange(field, 0, numberOfOriginalValues, result);
PerformEdgeInterpolations edgeInterpWorklet(numberOfOriginalValues);
vtkm::worklet::DispatcherMapField<PerformEdgeInterpolations> edgeInterpDispatcher(
edgeInterpWorklet);
edgeInterpDispatcher.Invoke(this->EdgeInterpolationArray, result);
// Perform a gather on output to get all required values for calculation of
// centroids using the interpolation info array.
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
using ValueHandle = vtkm::cont::ArrayHandle<ValueType>;
vtkm::cont::ArrayHandlePermutation<IdHandle, ValueHandle> toReduceValues(
InCellInterpolationInfo, result);
vtkm::cont::ArrayHandle<ValueType> reducedValues;
vtkm::worklet::DispatcherReduceByKey<PerformInCellInterpolations>
inCellInterpolationDispatcher;
inCellInterpolationDispatcher.Invoke(interpolationKeys, toReduceValues, reducedValues);
vtkm::Id inCellPointsOffset = numberOfOriginalValues + numberOfEdgePoints;
vtkm::cont::Algorithm::CopySubRange(
reducedValues, 0, reducedValues.GetNumberOfValues(), result, inCellPointsOffset);
*(this->Output) = result;
}
private:
vtkm::cont::ArrayHandle<EdgeInterpolation> EdgeInterpolationArray;
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
vtkm::Id EdgePointsOffset;
vtkm::Id InCellPointsOffset;
ArrayHandleType* Output;
};
template <typename ValueType, typename StorageType>
vtkm::cont::ArrayHandle<ValueType> ProcessPointField(
const vtkm::cont::ArrayHandle<ValueType, StorageType>& fieldData) const
struct PerformInCellInterpolations : public vtkm::worklet::WorkletReduceByKey
{
using ResultType = vtkm::cont::ArrayHandle<ValueType>;
using Worker = InterpolateField<ResultType>;
using ControlSignature = void(KeysIn keys, ValuesIn toReduce, ReducedValuesOut centroids);
using ExecutionSignature = void(_2, _3);
ResultType output;
template <typename MappedValueVecType, typename MappedValueType>
VTKM_EXEC void operator()(const MappedValueVecType& toReduce, MappedValueType& centroid) const
{
const vtkm::IdComponent numValues = toReduce.GetNumberOfComponents();
Worker worker = Worker(this->EdgePointsInterpolation,
this->InCellInterpolationKeys,
this->InCellInterpolationInfo,
this->EdgePointsOffset,
this->InCellPointsOffset,
&output);
worker(fieldData);
// Interpolate per-vertex because some vec-like objects do not allow intermediate variables
using VTraits = vtkm::VecTraits<MappedValueType>;
using CType = typename VTraits::ComponentType;
for (vtkm::IdComponent component = 0; component < VTraits::GetNumberOfComponents(centroid);
++component)
{
CType sum = VTraits::GetComponent(toReduce[0], component);
for (vtkm::IdComponent reduceI = 1; reduceI < numValues; ++reduceI)
{
// static_cast is for when MappedValueType is a small int that gets promoted to int32.
sum = static_cast<CType>(sum + VTraits::GetComponent(toReduce[reduceI], component));
}
VTraits::SetComponent(centroid, component, static_cast<CType>(sum / numValues));
}
}
};
return output;
template <typename InputType, typename OutputType>
void ProcessPointField(const InputType& input, OutputType& output)
{
if (!this->InterpolationKeysBuilt)
{
this->InterpolationKeys.BuildArrays(this->InCellInterpolationKeys, KeysSortType::Unstable);
}
vtkm::Id numberOfOriginalValues = input.GetNumberOfValues();
vtkm::Id numberOfEdgePoints = this->EdgePointsInterpolation.GetNumberOfValues();
vtkm::Id numberOfInCellPoints = this->InterpolationKeys.GetUniqueKeys().GetNumberOfValues();
// Copy over the original values. They are still part of the output. (Unused points are
// not culled. Use CleanGrid for that.)
output.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
vtkm::cont::Algorithm::CopySubRange(input, 0, numberOfOriginalValues, output);
// Interpolate all new points that lie on edges of the input mesh.
vtkm::cont::Invoker invoke;
invoke(PerformEdgeInterpolations{},
this->EdgePointsInterpolation,
input,
vtkm::cont::make_ArrayHandleView(output, numberOfOriginalValues, numberOfEdgePoints));
// Perform a gather on the output to get all the required values for calculation of centroids
// using the interpolation info array.
auto toReduceValues =
vtkm::cont::make_ArrayHandlePermutation(this->InCellInterpolationInfo, output);
invoke(PerformInCellInterpolations{},
this->InterpolationKeys,
toReduceValues,
vtkm::cont::make_ArrayHandleView(
output, numberOfOriginalValues + numberOfEdgePoints, numberOfInCellPoints));
}
vtkm::cont::ArrayHandle<vtkm::Id> GetCellMapOutputToInput() const
@ -912,6 +870,8 @@ private:
vtkm::cont::ArrayHandle<vtkm::Id> CellMapOutputToInput;
vtkm::Id EdgePointsOffset;
vtkm::Id InCellPointsOffset;
vtkm::worklet::Keys<vtkm::Id> InterpolationKeys;
bool InterpolationKeysBuilt = false;
};
}
} // namespace vtkm::worklet