diff --git a/docs/changelog/probe-sets-invalid-value.md b/docs/changelog/probe-sets-invalid-value.md new file mode 100644 index 000000000..1e179fd45 --- /dev/null +++ b/docs/changelog/probe-sets-invalid-value.md @@ -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. diff --git a/vtkm/cont/internal/CMakeLists.txt b/vtkm/cont/internal/CMakeLists.txt index c7f1c0247..e322a90ea 100644 --- a/vtkm/cont/internal/CMakeLists.txt +++ b/vtkm/cont/internal/CMakeLists.txt @@ -20,6 +20,7 @@ set(headers ArrayTransfer.h AtomicInterfaceControl.h AtomicInterfaceExecution.h + CastInvalidValue.h ConnectivityExplicitInternals.h DeviceAdapterAlgorithmGeneral.h DeviceAdapterListHelpers.h diff --git a/vtkm/cont/internal/CastInvalidValue.h b/vtkm/cont/internal/CastInvalidValue.h new file mode 100644 index 000000000..398639e97 --- /dev/null +++ b/vtkm/cont/internal/CastInvalidValue.h @@ -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 +#include +#include + +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 +T CastInvalidValue(vtkm::Float64 invalidValue) +{ + using ComponentType = typename vtkm::VecTraits::BaseComponentType; + + if (std::is_same::NumericTag>::value) + { + // Casting to integer types + if (vtkm::IsFinite(invalidValue)) + { + return T(static_cast(invalidValue)); + } + else if (vtkm::IsInf(invalidValue) && (invalidValue > 0)) + { + return T(std::numeric_limits::max()); + } + else + { + return T(std::numeric_limits::min()); + } + } + else + { + // Not an integer type. Assume can be directly cast + return T(static_cast(invalidValue)); + } +} +} +} +} // namespace vtkm::cont::internal + +#endif //vtk_m_cont_internal_CastInvalidValue_h diff --git a/vtkm/filter/MapFieldPermutation.cxx b/vtkm/filter/MapFieldPermutation.cxx index d7948fae6..76a92f61d 100644 --- a/vtkm/filter/MapFieldPermutation.cxx +++ b/vtkm/filter/MapFieldPermutation.cxx @@ -16,6 +16,8 @@ #include +#include + #include #include @@ -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 -T CastInvalidValue(vtkm::Float64 invalidValue) -{ - using ComponentType = typename vtkm::VecTraits::BaseComponentType; - - if (std::is_same::NumericTag>::value) - { - // Casting to integer types - if (vtkm::IsFinite(invalidValue)) - { - return T(static_cast(invalidValue)); - } - else if (vtkm::IsInf(invalidValue) && (invalidValue > 0)) - { - return T(std::numeric_limits::max()); - } - else - { - return T(std::numeric_limits::min()); - } - } - else - { - // Not an integer type. Assume can be directly cast - return T(static_cast(invalidValue)); - } -} - struct DoMapFieldPermutation { bool CalledMap = false; @@ -91,7 +62,7 @@ struct DoMapFieldPermutation vtkm::Float64 invalidValue) { vtkm::cont::ArrayHandle outputArray; - MapPermutationWorklet worklet(CastInvalidValue(invalidValue)); + MapPermutationWorklet worklet(vtkm::cont::internal::CastInvalidValue(invalidValue)); vtkm::cont::Invoker invoke; invoke(worklet, permutation, inputArray, outputArray); output = vtkm::cont::VariantArrayHandle(outputArray); diff --git a/vtkm/filter/Probe.h b/vtkm/filter/Probe.h index db8e90728..55ef54d80 100644 --- a/vtkm/filter/Probe.h +++ b/vtkm/filter/Probe.h @@ -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 VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input, vtkm::filter::PolicyBase policy); @@ -44,6 +50,7 @@ public: private: vtkm::cont::DataSet Geometry; vtkm::worklet::Probe Worklet; + vtkm::Float64 InvalidValue = vtkm::Nan64(); }; } } // vtkm::filter diff --git a/vtkm/filter/Probe.hxx b/vtkm/filter/Probe.hxx index e7266f3e5..8e164bfa3 100644 --- a/vtkm/filter/Probe.hxx +++ b/vtkm/filter/Probe.hxx @@ -14,6 +14,8 @@ #include +#include + 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 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(this->InvalidValue), + typename DerivedPolicy::AllCellSetList()); result.AddField(fieldMeta.AsField(fieldArray)); return true; } diff --git a/vtkm/filter/testing/UnitTestProbe.cxx b/vtkm/filter/testing/UnitTestProbe.cxx index 30402231c..2928b10ec 100644 --- a/vtkm/filter/testing/UnitTestProbe.cxx +++ b/vtkm/filter/testing/UnitTestProbe.cxx @@ -67,13 +67,20 @@ vtkm::cont::DataSet ConvertDataSetUniformToExplicit(const vtkm::cont::DataSet& u const std::vector& GetExpectedPointData() { static std::vector 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; } diff --git a/vtkm/worklet/Probe.h b/vtkm/worklet/Probe.h index 47c8a16b0..14603bca9 100644 --- a/vtkm/worklet/Probe.h +++ b/vtkm/worklet/Probe.h @@ -178,9 +178,16 @@ public: } //============================================================================ + template 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 ProcessPointField( const vtkm::cont::ArrayHandle& field, + const T& invalidValue, InputCellSetTypeList icsTypes = InputCellSetTypeList()) const { vtkm::cont::ArrayHandle result; - vtkm::worklet::DispatcherMapField dispatcher; - dispatcher.Invoke(this->CellIds, - this->ParametricCoordinates, - this->InputCellSet.ResetCellSetList(icsTypes), - field, - result); + vtkm::cont::Invoker invoke; + invoke(InterpolatePointField(invalidValue), + this->CellIds, + this->ParametricCoordinates, + this->InputCellSet.ResetCellSetList(icsTypes), + field, + result); return result; } //============================================================================ + template 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 - vtkm::cont::ArrayHandle ProcessCellField( - const vtkm::cont::ArrayHandle& field) const + vtkm::cont::ArrayHandle ProcessCellField(const vtkm::cont::ArrayHandle& field, + const T& invalidValue) const { vtkm::cont::ArrayHandle result; - vtkm::worklet::DispatcherMapField dispatcher; - dispatcher.Invoke(this->CellIds, field, result); + vtkm::cont::Invoker invoke; + invoke(MapCellField(invalidValue), this->CellIds, field, result); return result; } diff --git a/vtkm/worklet/testing/UnitTestProbe.cxx b/vtkm/worklet/testing/UnitTestProbe.cxx index 6be8323b1..b6006a5f1 100644 --- a/vtkm/worklet/testing/UnitTestProbe.cxx +++ b/vtkm/worklet/testing/UnitTestProbe.cxx @@ -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()); - auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast()); + auto pf = + probe.ProcessPointField(input.GetField("pointdata").GetData().Cast(), 0.0f); + auto cf = + probe.ProcessCellField(input.GetField("celldata").GetData().Cast(), 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()); - auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast()); + auto pf = + probe.ProcessPointField(input.GetField("pointdata").GetData().Cast(), 0.0f); + auto cf = + probe.ProcessCellField(input.GetField("celldata").GetData().Cast(), 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()); - auto cf = probe.ProcessCellField(input.GetField("celldata").GetData().Cast()); + auto pf = + probe.ProcessPointField(input.GetField("pointdata").GetData().Cast(), 0.0f); + auto cf = + probe.ProcessCellField(input.GetField("celldata").GetData().Cast(), 0.0f); auto hp = probe.GetHiddenPointsField(); auto hc = probe.GetHiddenCellsField(geometry.GetCellSet());