Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-07-30 18:50:59 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Enable setting invalid value in probe filter

Browse Source 
Initially, the probe filter would simply not set a value if a sample was
outside the input `DataSet`. This is not great as the memory could be
left uninitalized and lead to unpredictable results. The testing
compared these invalid results to 0, which seemed to work but is
probably unstable.

This was partially fixed by a previous change that consolidated to
mapping of cell data with a general routine that permuted data. However,
the fix did not extend to point data in the input, and it was not
possible to specify a particular invalid value.

This change specifically updates the probe filter so that invalid values
are set to a user-specified value.
This commit is contained in:
Kenneth Moreland 2020-06-03 15:24:36 -06:00
parent c1c8e54940
commit 7be44c847f
9 changed files with 166 additions and 57 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
docs/changelog/probe-sets-invalid-value.md Normal file
View File

@ -0,0 +1,15 @@
# Enable setting invalid value in probe filter
Initially, the probe filter would simply not set a value if a sample was
outside the input `DataSet`. This is not great as the memory could be
left uninitalized and lead to unpredictable results. The testing
compared these invalid results to 0, which seemed to work but is
probably unstable.
This was partially fixed by a previous change that consolidated to
mapping of cell data with a general routine that permuted data. However,
the fix did not extend to point data in the input, and it was not
possible to specify a particular invalid value.
This change specifically updates the probe filter so that invalid values
are set to a user-specified value.

1
vtkm/cont/internal/CMakeLists.txt
View File

@ -20,6 +20,7 @@ set(headers
ArrayTransfer.h
AtomicInterfaceControl.h
AtomicInterfaceExecution.h
CastInvalidValue.h
ConnectivityExplicitInternals.h
DeviceAdapterAlgorithmGeneral.h
DeviceAdapterListHelpers.h

67
vtkm/cont/internal/CastInvalidValue.h Normal file
View File

@ -0,0 +1,67 @@
//============================================================================
// 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_cont_internal_CastInvalidValue_h
#define vtk_m_cont_internal_CastInvalidValue_h
#include <vtkm/Math.h>
#include <vtkm/Types.h>
#include <vtkm/VecTraits.h>
namespace vtkm
{
namespace cont
{
namespace internal
{
/// \brief Convert an invalid value to something type-appropriate.
///
/// There are algorithms in VTK-m that require a placeholder for invalid values in an array
/// or field. For example, when probing something, a probe location outside of the source data
/// has to be set to something.
///
/// Often we want to set this to something like NaN to make it clear that this is invalid.
/// However, integer types cannot represent these non-finite numbers.
///
/// For convenience, it is easiest to allow the user to specify the invalid value as a
/// vtkm::Float64 and use this function to convert it to something type-appropriate.
///
template <typename T>
T CastInvalidValue(vtkm::Float64 invalidValue)
{
using ComponentType = typename vtkm::VecTraits<T>::BaseComponentType;
if (std::is_same<vtkm::TypeTraitsIntegerTag, typename vtkm::TypeTraits<T>::NumericTag>::value)
{
// Casting to integer types
if (vtkm::IsFinite(invalidValue))
{
return T(static_cast<ComponentType>(invalidValue));
}
else if (vtkm::IsInf(invalidValue) && (invalidValue > 0))
{
return T(std::numeric_limits<ComponentType>::max());
}
else
{
return T(std::numeric_limits<ComponentType>::min());
}
}
else
{
// Not an integer type. Assume can be directly cast
return T(static_cast<ComponentType>(invalidValue));
}
}
}
}
} // namespace vtkm::cont::internal
#endif //vtk_m_cont_internal_CastInvalidValue_h

35
vtkm/filter/MapFieldPermutation.cxx
View File

@ -16,6 +16,8 @@
#include <vtkm/cont/Logging.h>
#include <vtkm/cont/internal/CastInvalidValue.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/filter/PolicyDefault.h>
@ -49,37 +51,6 @@ struct MapPermutationWorklet : vtkm::worklet::WorkletMapField
}
};
// For simplicity, the invalid value is specified as a single type (vtkm::Float64), and this is
// often a non-finite value, which is not well represented by integer types. This function does its
// best to find a reasonable cast for the value.
template <typename T>
T CastInvalidValue(vtkm::Float64 invalidValue)
{
using ComponentType = typename vtkm::VecTraits<T>::BaseComponentType;
if (std::is_same<vtkm::TypeTraitsIntegerTag, typename vtkm::TypeTraits<T>::NumericTag>::value)
{
// Casting to integer types
if (vtkm::IsFinite(invalidValue))
{
return T(static_cast<ComponentType>(invalidValue));
}
else if (vtkm::IsInf(invalidValue) && (invalidValue > 0))
{
return T(std::numeric_limits<ComponentType>::max());
}
else
{
return T(std::numeric_limits<ComponentType>::min());
}
}
else
{
// Not an integer type. Assume can be directly cast
return T(static_cast<ComponentType>(invalidValue));
}
}
struct DoMapFieldPermutation
{
bool CalledMap = false;
@ -91,7 +62,7 @@ struct DoMapFieldPermutation
vtkm::Float64 invalidValue)
{
vtkm::cont::ArrayHandle<T> outputArray;
MapPermutationWorklet<T> worklet(CastInvalidValue<T>(invalidValue));
MapPermutationWorklet<T> worklet(vtkm::cont::internal::CastInvalidValue<T>(invalidValue));
vtkm::cont::Invoker invoke;
invoke(worklet, permutation, inputArray, outputArray);
output = vtkm::cont::VariantArrayHandle(outputArray);

7
vtkm/filter/Probe.h
View File

@ -24,6 +24,12 @@ public:
VTKM_CONT
void SetGeometry(const vtkm::cont::DataSet& geometry);
VTKM_CONT
const vtkm::cont::DataSet& GetGeometry() const;
VTKM_CONT void SetInvalidValue(vtkm::Float64 invalidValue) { this->InvalidValue = invalidValue; }
VTKM_CONT vtkm::Float64 GetInvalidValue() const { return this->InvalidValue; }
template <typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
@ -44,6 +50,7 @@ public:
private:
vtkm::cont::DataSet Geometry;
vtkm::worklet::Probe Worklet;
vtkm::Float64 InvalidValue = vtkm::Nan64();
};
}
} // vtkm::filter

15
vtkm/filter/Probe.hxx
View File

@ -14,6 +14,8 @@
#include <vtkm/filter/MapFieldPermutation.h>
#include <vtkm/cont/internal/CastInvalidValue.h>
namespace vtkm
{
namespace filter
@ -27,6 +29,12 @@ inline void Probe::SetGeometry(const vtkm::cont::DataSet& geometry)
this->Geometry.AddCoordinateSystem(geometry.GetCoordinateSystem());
}
VTKM_CONT
inline const vtkm::cont::DataSet& Probe::GetGeometry() const
{
return this->Geometry;
}
template <typename DerivedPolicy>
VTKM_CONT inline vtkm::cont::DataSet Probe::DoExecute(
const vtkm::cont::DataSet& input,
@ -61,7 +69,8 @@ VTKM_CONT inline bool Probe::MapFieldOntoOutput(vtkm::cont::DataSet& result,
}
else if (field.IsFieldCell())
{
return vtkm::filter::MapFieldPermutation(field, this->Worklet.GetCellIds(), result);
return vtkm::filter::MapFieldPermutation(
field, this->Worklet.GetCellIds(), result, this->InvalidValue);
}
else if (field.IsFieldGlobal())
{
@ -82,7 +91,9 @@ VTKM_CONT inline bool Probe::DoMapField(vtkm::cont::DataSet& result,
{
VTKM_ASSERT(fieldMeta.IsPointField());
auto fieldArray =
this->Worklet.ProcessPointField(input, typename DerivedPolicy::AllCellSetList());
this->Worklet.ProcessPointField(input,
vtkm::cont::internal::CastInvalidValue<T>(this->InvalidValue),
typename DerivedPolicy::AllCellSetList());
result.AddField(fieldMeta.AsField(fieldArray));
return true;
}

21
vtkm/filter/testing/UnitTestProbe.cxx
View File

@ -67,13 +67,20 @@ vtkm::cont::DataSet ConvertDataSetUniformToExplicit(const vtkm::cont::DataSet& u
const std::vector<vtkm::Float32>& GetExpectedPointData()
{
static std::vector<vtkm::Float32> expected = {
1.05f, 1.155f, 1.26f, 1.365f, 1.47f, 1.575f, 1.68f, 0.0f, 0.0f, 1.47f, 1.575f, 1.68f,
1.785f, 1.89f, 1.995f, 2.1f, 0.0f, 0.0f, 1.89f, 1.995f, 2.1f, 2.205f, 2.31f, 2.415f,
2.52f, 0.0f, 0.0f, 2.31f, 2.415f, 2.52f, 2.625f, 2.73f, 2.835f, 2.94f, 0.0f, 0.0f,
2.73f, 2.835f, 2.94f, 3.045f, 3.15f, 3.255f, 3.36f, 0.0f, 0.0f, 3.15f, 3.255f, 3.36f,
3.465f, 3.57f, 3.675f, 3.78f, 0.0f, 0.0f, 3.57f, 3.675f, 3.78f, 3.885f, 3.99f, 4.095f,
4.2f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f,
0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f
1.05f, 1.155f, 1.26f, 1.365f, 1.47f, 1.575f,
1.68f, vtkm::Nan32(), vtkm::Nan32(), 1.47f, 1.575f, 1.68f,
1.785f, 1.89f, 1.995f, 2.1f, vtkm::Nan32(), vtkm::Nan32(),
1.89f, 1.995f, 2.1f, 2.205f, 2.31f, 2.415f,
2.52f, vtkm::Nan32(), vtkm::Nan32(), 2.31f, 2.415f, 2.52f,
2.625f, 2.73f, 2.835f, 2.94f, vtkm::Nan32(), vtkm::Nan32(),
2.73f, 2.835f, 2.94f, 3.045f, 3.15f, 3.255f,
3.36f, vtkm::Nan32(), vtkm::Nan32(), 3.15f, 3.255f, 3.36f,
3.465f, 3.57f, 3.675f, 3.78f, vtkm::Nan32(), vtkm::Nan32(),
3.57f, 3.675f, 3.78f, 3.885f, 3.99f, 4.095f,
4.2f, vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(),
vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(),
vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32(),
vtkm::Nan32(), vtkm::Nan32(), vtkm::Nan32()
};
return expected;
}

44
vtkm/worklet/Probe.h
View File

@ -178,9 +178,16 @@ public:
}
//============================================================================
template <typename T>
class InterpolatePointField : public vtkm::worklet::WorkletMapField
{
public:
T InvalidValue;
InterpolatePointField(const T& invalidValue)
: InvalidValue(invalidValue)
{
}
using ControlSignature = void(FieldIn cellIds,
FieldIn parametricCoords,
WholeCellSetIn<> inputCells,
@ -201,6 +208,10 @@ public:
auto pointVals = vtkm::make_VecFromPortalPermute(&indices, in);
vtkm::exec::CellInterpolate(pointVals, pc, cells.GetCellShape(cellId), out);
}
else
{
out = this->InvalidValue;
}
}
};
@ -210,23 +221,32 @@ public:
typename InputCellSetTypeList = VTKM_DEFAULT_CELL_SET_LIST>
vtkm::cont::ArrayHandle<T> ProcessPointField(
const vtkm::cont::ArrayHandle<T, Storage>& field,
const T& invalidValue,
InputCellSetTypeList icsTypes = InputCellSetTypeList()) const
{
vtkm::cont::ArrayHandle<T> result;
vtkm::worklet::DispatcherMapField<InterpolatePointField> dispatcher;
dispatcher.Invoke(this->CellIds,
this->ParametricCoordinates,
this->InputCellSet.ResetCellSetList(icsTypes),
field,
result);
vtkm::cont::Invoker invoke;
invoke(InterpolatePointField<T>(invalidValue),
this->CellIds,
this->ParametricCoordinates,
this->InputCellSet.ResetCellSetList(icsTypes),
field,
result);
return result;
}
//============================================================================
template <typename T>
class MapCellField : public vtkm::worklet::WorkletMapField
{
public:
T InvalidValue;
MapCellField(const T& invalidValue)
: InvalidValue(invalidValue)
{
}
using ControlSignature = void(FieldIn cellIds, WholeArrayIn inputField, FieldOut result);
using ExecutionSignature = void(_1, _2, _3);
@ -239,6 +259,10 @@ public:
{
out = in.Get(cellId);
}
else
{
out = this->InvalidValue;
}
}
};
@ -247,12 +271,12 @@ public:
/// cell is chosen arbitrarily.
///
template <typename T, typename Storage>
vtkm::cont::ArrayHandle<T> ProcessCellField(
const vtkm::cont::ArrayHandle<T, Storage>& field) const
vtkm::cont::ArrayHandle<T> ProcessCellField(const vtkm::cont::ArrayHandle<T, Storage>& field,
const T& invalidValue) const
{
vtkm::cont::ArrayHandle<T> result;
vtkm::worklet::DispatcherMapField<MapCellField> dispatcher;
dispatcher.Invoke(this->CellIds, field, result);
vtkm::cont::Invoker invoke;
invoke(MapCellField<T>(invalidValue), this->CellIds, field, result);
return result;
}

18
vtkm/worklet/testing/UnitTestProbe.cxx
View File

@ -133,8 +133,10 @@ private:
vtkm::worklet::Probe probe;
probe.Run(input.GetCellSet(), input.GetCoordinateSystem(), geometry.GetCoordinateSystem());
auto pf = probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>());
auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>());
auto pf =
probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>(), 0.0f);
auto cf =
probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>(), 0.0f);
auto hp = probe.GetHiddenPointsField();
auto hc = probe.GetHiddenCellsField(geometry.GetCellSet());
@ -154,8 +156,10 @@ private:
vtkm::worklet::Probe probe;
probe.Run(input.GetCellSet(), input.GetCoordinateSystem(), geometry.GetCoordinateSystem());
auto pf = probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>());
auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>());
auto pf =
probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>(), 0.0f);
auto cf =
probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>(), 0.0f);
auto hp = probe.GetHiddenPointsField();
auto hc = probe.GetHiddenCellsField(geometry.GetCellSet());
@ -176,8 +180,10 @@ private:
vtkm::worklet::Probe probe;
probe.Run(input.GetCellSet(), input.GetCoordinateSystem(), geometry.GetCoordinateSystem());
auto pf = probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>());
auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>());
auto pf =
probe.ProcessPointField(input.GetField("pointdata").GetData().Cast<FieldArrayType>(), 0.0f);
auto cf =
probe.ProcessCellField(input.GetField("celldata").GetData().Cast<FieldArrayType>(), 0.0f);
auto hp = probe.GetHiddenPointsField();
auto hc = probe.GetHiddenCellsField(geometry.GetCellSet());