diff --git a/vtkm/ErrorCode.h b/vtkm/ErrorCode.h index 7ef94a0fcf89f92c8fabe15f7a50dd198fe975c8..ba35d7240f3158df79ac0e68ded2e74982c14e88 100644 --- a/vtkm/ErrorCode.h +++ b/vtkm/ErrorCode.h @@ -42,6 +42,7 @@ enum class ErrorCode MalformedCellDetected, OperationOnEmptyCell, CellNotFound, + ValueOutOfRange, UnknownError }; @@ -160,6 +161,8 @@ inline const char* ErrorString(vtkm::ErrorCode code) noexcept return "Operation on empty cell"; case vtkm::ErrorCode::CellNotFound: return "Cell not found"; + case vtkm::ErrorCode::ValueOutOfRange: + return "Value out of range"; case vtkm::ErrorCode::UnknownError: return "Unknown error"; } diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index ca47f3cbb673e02b74e6e188fe27d2e0ce40a0b4..c686c386daed19746bd24632a3a0a621f17c20c7 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -80,6 +80,7 @@ set(headers ColorTableSamples.h ConvertNumComponentsToOffsets.h CoordinateSystem.h + CubicHermiteSpline.h DataSet.h DataSetBuilderCurvilinear.h DataSetBuilderExplicit.h @@ -199,6 +200,7 @@ set(device_sources CellSetExtrude.cxx ColorTable.cxx ConvertNumComponentsToOffsets.cxx + CubicHermiteSpline.cxx Field.cxx internal/ArrayCopyUnknown.cxx internal/ArrayRangeComputeUtils.cxx diff --git a/vtkm/cont/CubicHermiteSpline.cxx b/vtkm/cont/CubicHermiteSpline.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d8e739ab73c732c7559061e1bc08752304a83268 --- /dev/null +++ b/vtkm/cont/CubicHermiteSpline.cxx @@ -0,0 +1,148 @@ +//============================================================================ +// 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 +#include +#include +#include + +namespace vtkm +{ +namespace cont +{ +namespace +{ +struct CalcNeighborDistanceWorklet : public vtkm::worklet::WorkletMapField +{ + using ControlSignature = void(FieldOut, WholeArrayIn); + using ExecutionSignature = void(InputIndex, _1, _2); + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, + vtkm::FloatDefault& val, + const ArrayType& data) const + { + if (idx == 0) + val = 0.0; + else + val = vtkm::Magnitude(data.Get(idx) - data.Get(idx - 1)); + } +}; + +struct CalcTangentsWorklet : public vtkm::worklet::WorkletMapField +{ + CalcTangentsWorklet(const vtkm::Id& numPoints) + : NumPoints(numPoints) + { + } + + using ControlSignature = void(FieldOut, WholeArrayIn, WholeArrayIn); + using ExecutionSignature = void(InputIndex, _1, _2, _3); + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, + TangentArrayType& tangent, + const PointArrayType& points, + const KnotArrayType& knots) const + { + vtkm::Id idx0, idx1; + if (idx == 0) // Forward difference + { + idx0 = 0; + idx1 = 1; + } + else if (idx == NumPoints - 1) // Backward difference + { + idx0 = NumPoints - 2; + idx1 = NumPoints - 1; + } + else // central difference + { + idx0 = idx - 1; + idx1 = idx + 1; + } + + auto dX = points.Get(idx1) - points.Get(idx0); + auto dT = knots.Get(idx1) - knots.Get(idx0); + + tangent = dX / dT; + } + + vtkm::Id NumPoints; +}; +} //anonymous namespace + +VTKM_CONT vtkm::Range CubicHermiteSpline::GetParametricRange() +{ + auto n = this->Knots.GetNumberOfValues(); + if (n == 0) + { + this->ComputeKnots(); + n = this->Knots.GetNumberOfValues(); + } + const auto ids = vtkm::cont::make_ArrayHandle({ 0, n - 1 }); + const std::vector output = vtkm::cont::ArrayGetValues(ids, this->Knots); + return { output[0], output[1] }; +} + +VTKM_CONT void CubicHermiteSpline::ComputeKnots() +{ + vtkm::Id n = this->Data.GetNumberOfValues(); + this->Knots.Allocate(n); + + vtkm::cont::Invoker invoker; + + //uses chord length parameterization. + invoker(CalcNeighborDistanceWorklet{}, this->Knots, this->Data); + vtkm::FloatDefault sum = vtkm::cont::Algorithm::ScanInclusive(this->Knots, this->Knots); + + if (sum == 0.0) + throw std::invalid_argument("Error: accumulated distance between data is zero."); + + auto divideBy = vtkm::cont::make_ArrayHandleConstant(1.0 / sum, this->Knots.GetNumberOfValues()); + vtkm::cont::Algorithm::Transform(this->Knots, divideBy, this->Knots, vtkm::Product{}); +} + +VTKM_CONT void CubicHermiteSpline::ComputeTangents() +{ + vtkm::Id n = this->Data.GetNumberOfValues(); + this->Tangents.Allocate(n); + + vtkm::cont::Invoker invoker; + + invoker(CalcTangentsWorklet{ n }, this->Tangents, this->Data, this->Knots); +} + +vtkm::exec::CubicHermiteSpline CubicHermiteSpline::PrepareForExecution( + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) +{ + vtkm::Id n = this->Data.GetNumberOfValues(); + if (n < 2) + throw vtkm::cont::ErrorBadValue("At least two points are required for spline interpolation."); + if (this->Knots.GetNumberOfValues() == 0) + this->ComputeKnots(); + if (this->Tangents.GetNumberOfValues() == 0) + this->ComputeTangents(); + + if (n != this->Knots.GetNumberOfValues()) + throw vtkm::cont::ErrorBadValue("Number of data points must match the number of knots."); + if (n != this->Tangents.GetNumberOfValues()) + throw vtkm::cont::ErrorBadValue("Number of data points must match the number of tangents."); + + using ExecObjType = vtkm::exec::CubicHermiteSpline; + + return ExecObjType(this->Data, this->Knots, this->Tangents, device, token); +} +} +} // namespace vtkm::cont diff --git a/vtkm/cont/CubicHermiteSpline.h b/vtkm/cont/CubicHermiteSpline.h new file mode 100644 index 0000000000000000000000000000000000000000..92f3cb114c6e114d324d008773e4fef4a38f20c9 --- /dev/null +++ b/vtkm/cont/CubicHermiteSpline.h @@ -0,0 +1,90 @@ +//============================================================================ +// 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_CubicHermiteSpline_h +#define vtk_m_cont_CubicHermiteSpline_h + +#include + +#include +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace cont +{ + +class VTKM_CONT_EXPORT CubicHermiteSpline : public vtkm::cont::ExecutionObjectBase +{ +public: + CubicHermiteSpline() = default; + virtual ~CubicHermiteSpline() = default; + + VTKM_CONT void SetData(const vtkm::cont::ArrayHandle& data) { this->Data = data; } + VTKM_CONT void SetData(const std::vector& data, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + { + this->Data = vtkm::cont::make_ArrayHandle(data, copy); + } + + VTKM_CONT void SetKnots(const vtkm::cont::ArrayHandle& knots) + { + this->Knots = knots; + } + VTKM_CONT void SetKnots(const std::vector& knots, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + { + this->Knots = vtkm::cont::make_ArrayHandle(knots, copy); + } + + VTKM_CONT void SetTangents(const vtkm::cont::ArrayHandle& tangents) + { + this->Tangents = tangents; + } + VTKM_CONT void SetTangents(const std::vector& tangents, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + { + this->Tangents = vtkm::cont::make_ArrayHandle(tangents, copy); + } + + VTKM_CONT vtkm::Range GetParametricRange(); + + VTKM_CONT const vtkm::cont::ArrayHandle& GetData() const { return this->Data; } + VTKM_CONT const vtkm::cont::ArrayHandle& GetTangents() + { + if (this->Tangents.GetNumberOfValues() == 0) + this->ComputeTangents(); + return this->Tangents; + } + VTKM_CONT const vtkm::cont::ArrayHandle& GetKnots() + { + if (this->Knots.GetNumberOfValues() == 0) + this->ComputeKnots(); + return this->Knots; + } + + VTKM_CONT vtkm::exec::CubicHermiteSpline PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token); + +private: + VTKM_CONT void ComputeKnots(); + VTKM_CONT void ComputeTangents(); + + vtkm::cont::ArrayHandle Data; + vtkm::cont::ArrayHandle Knots; + vtkm::cont::ArrayHandle Tangents; +}; +} +} // vtkm::cont + +#endif //vtk_m_cont_CubicHermiteSpline_h diff --git a/vtkm/cont/testing/CMakeLists.txt b/vtkm/cont/testing/CMakeLists.txt index 250a5f787a7ba92acc5a7a42ed06a98565d4c793..5041d4d405c6a243ad3ee6e01f221fec5300c4cf 100644 --- a/vtkm/cont/testing/CMakeLists.txt +++ b/vtkm/cont/testing/CMakeLists.txt @@ -105,6 +105,7 @@ set(unit_tests_device UnitTestCellSetExplicit.cxx UnitTestCellSetPermutation.cxx UnitTestColorTable.cxx + UnitTestCubicHermiteSpline.cxx UnitTestDataSetPermutation.cxx UnitTestDataSetSingleType.cxx UnitTestDeviceAdapterAlgorithmDependency.cxx diff --git a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4d145c6a21c40136df27ee9a0cc210453f2333c0 --- /dev/null +++ b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx @@ -0,0 +1,142 @@ +//============================================================================ +// 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 +#include +#include +#include + +namespace +{ + +class SplineEvalWorklet : public vtkm::worklet::WorkletMapField +{ +public: + SplineEvalWorklet() {} + + using ControlSignature = void(FieldIn param, ExecObject cubicSpline, FieldOut value); + using ExecutionSignature = void(_1, _2, _3); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const ParamType& param, + const CubicSplineType& spline, + ResultType& value) const + { + auto res = spline.Evaluate(param, value); + if (res != vtkm::ErrorCode::Success) + this->RaiseError("Spline evaluation failed."); + } +}; + +void CheckEvaluation(const vtkm::cont::CubicHermiteSpline& spline, + const vtkm::cont::ArrayHandle& params, + const std::vector& answer) +{ + vtkm::cont::Invoker invoke; + vtkm::cont::ArrayHandle result; + invoke(SplineEvalWorklet{}, params, spline, result); + + VTKM_TEST_ASSERT( + test_equal_ArrayHandles(result, vtkm::cont::make_ArrayHandle(answer, vtkm::CopyFlag::Off))); +} + +void CheckEvaluation(const vtkm::cont::CubicHermiteSpline& spline, + const std::vector& params, + const std::vector& answer) +{ + return CheckEvaluation(spline, vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off), answer); +} + +void CubicHermiteSplineTest() +{ + std::vector pts = { { 0, 0, 0 }, { 1, 1, 1 }, { 2, 1, 0 }, { 3, -.5, -1 }, + { 4, -1, 0 }, { 5, -1, 1 }, { 6, 0, 0 } }; + + vtkm::cont::CubicHermiteSpline spline; + spline.SetData(pts); + //Evaluation at knots gives the sample pts. + CheckEvaluation(spline, spline.GetKnots(), pts); + + //Evaluation at non-knot values. + std::vector params = { 0.21, 0.465, 0.501, 0.99832 }; + std::vector result = { { 1.23261, 1.08861, 0.891725 }, + { 2.68524, -0.0560059, -0.855685 }, + { 2.85574, -0.32766, -0.970523 }, + { 5.99045, -0.00959875, 0.00964856 } }; + CheckEvaluation(spline, params, result); + + //Explicitly set knots and check. + std::vector knots = { 0, 1, 2, 3, 4, 5, 6 }; + spline = vtkm::cont::CubicHermiteSpline(); + spline.SetData(pts); + spline.SetKnots(knots); + CheckEvaluation(spline, knots, pts); + + //Evaluation at non-knot values. + params = { 0.84, 1.399, 2.838, 4.930, 5.001, 5.993 }; + result = { { 0.84, 0.896448, 0.952896 }, { 1.399, 1.14382, 0.745119 }, + { 2.838, -0.297388, -0.951764 }, { 4.93, -1.03141, 0.990543 }, + { 5.001, -0.999499, 0.999998 }, { 5.993, -0.00702441, 0.00704873 } }; + CheckEvaluation(spline, params, result); + + //Non-uniform knots. + knots = { 0, 1, 2, 2.1, 2.2, 2.3, 3 }; + spline = vtkm::cont::CubicHermiteSpline(); + spline.SetData(pts); + spline.SetKnots(knots); + CheckEvaluation(spline, knots, pts); + + params = { 1.5, 2.05, 2.11, 2.299, 2.8 }; + result = { { 1.39773, 1.23295, 0.727273 }, + { 2.39773, 0.357954, -0.522727 }, + { 3.1, -0.59275, -0.981 }, + { 4.99735, -1.00125, 0.999801 }, + { 5.75802, -0.293003, 0.344023 } }; + CheckEvaluation(spline, params, result); + + //Create a more complex spline from analytical functions. + vtkm::Id n = 500; + vtkm::FloatDefault t = 0.0, dt = vtkm::TwoPi() / static_cast(n); + + pts.clear(); + knots.clear(); + while (t <= vtkm::TwoPi()) + { + vtkm::FloatDefault x = vtkm::Cos(t); + vtkm::FloatDefault y = vtkm::Sin(t); + vtkm::FloatDefault z = x * y; + pts.push_back({ x, y, z }); + knots.push_back(t); + t += dt; + } + spline = vtkm::cont::CubicHermiteSpline(); + spline.SetData(pts); + spline.SetKnots(knots); + CheckEvaluation(spline, knots, pts); + + //Evaluate at a few points and check against analytical results. + params = { 0.15, 1.83, 2.38, 3.0291, 3.8829, 4.92, 6.2 }; + result.clear(); + for (const auto& p : params) + result.push_back({ vtkm::Cos(p), vtkm::Sin(p), vtkm::Cos(p) * vtkm::Sin(p) }); + CheckEvaluation(spline, params, result); +} + +} // anonymous namespace + +int UnitTestCubicHermiteSpline(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(CubicHermiteSplineTest, argc, argv); +} diff --git a/vtkm/exec/CMakeLists.txt b/vtkm/exec/CMakeLists.txt index 2f27cb90771b6a45ee2e6b33ac0463ae0c8b2b64..621e6cc863a84583abaa40d8cacca9dadb8bf11d 100644 --- a/vtkm/exec/CMakeLists.txt +++ b/vtkm/exec/CMakeLists.txt @@ -29,6 +29,8 @@ set(headers ConnectivityExtrude.h ConnectivityPermuted.h ConnectivityStructured.h + CubicHermiteSpline.h + CubicSpline.h FieldNeighborhood.h FunctorBase.h ParametricCoordinates.h diff --git a/vtkm/exec/CubicHermiteSpline.h b/vtkm/exec/CubicHermiteSpline.h new file mode 100644 index 0000000000000000000000000000000000000000..0b49686ebe961b5cf95151ef79c7b0d34d7c3dad --- /dev/null +++ b/vtkm/exec/CubicHermiteSpline.h @@ -0,0 +1,112 @@ +//============================================================================ +// 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_exec_CubicHermiteSpline_h +#define vtk_m_exec_CubicHermiteSpline_h + +#include +#include +#include + +namespace vtkm +{ +namespace exec +{ + +class VTKM_ALWAYS_EXPORT CubicHermiteSpline +{ +private: + using DataArrayHandleType = vtkm::cont::ArrayHandle; + using KnotsArrayHandleType = vtkm::cont::ArrayHandle; + using TangentArrayHandleType = vtkm::cont::ArrayHandle; + using DataPortalType = typename DataArrayHandleType::ReadPortalType; + using KnotsPortalType = typename KnotsArrayHandleType::ReadPortalType; + using TangentPortalType = typename TangentArrayHandleType::ReadPortalType; + +public: + VTKM_CONT CubicHermiteSpline(const vtkm::cont::ArrayHandle& data, + const vtkm::cont::ArrayHandle& knots, + const vtkm::cont::ArrayHandle& tangents, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) + { + this->Data = data.PrepareForInput(device, token); + this->Knots = knots.PrepareForInput(device, token); + this->Tangents = tangents.PrepareForInput(device, token); + } + + VTKM_EXEC + vtkm::ErrorCode Evaluate(const vtkm::FloatDefault& tVal, vtkm::Vec3f& val) const + { + vtkm::Id idx = this->FindInterval(tVal); + if (idx < 0) + { + idx = this->FindInterval(tVal); + return vtkm::ErrorCode::ValueOutOfRange; + } + + auto m0 = this->Tangents.Get(idx); + auto m1 = this->Tangents.Get(idx + 1); + + vtkm::FloatDefault t0 = this->Knots.Get(idx); + vtkm::FloatDefault t1 = this->Knots.Get(idx + 1); + vtkm::FloatDefault t = (tVal - t0) / (t1 - t0); + vtkm::FloatDefault t2 = t * t, t3 = t2 * t; + + // Hermite basis functions. + vtkm::FloatDefault h00 = 2 * t3 - 3 * t2 + 1; + vtkm::FloatDefault h10 = t3 - 2 * t2 + t; + vtkm::FloatDefault h01 = -2 * t3 + 3 * t2; + vtkm::FloatDefault h11 = t3 - t2; + + const auto& d0 = this->Data.Get(idx); + const auto& d1 = this->Data.Get(idx + 1); + for (vtkm::Id i = 0; i < 3; ++i) + val[i] = h00 * d0[i] + h10 * (t1 - t0) * m0[i] + h01 * d1[i] + h11 * (t1 - t0) * m1[i]; + + return vtkm::ErrorCode::Success; + } + +private: + vtkm::Id FindInterval(const vtkm::FloatDefault& t) const + { + vtkm::Id n = this->Knots.GetNumberOfValues(); + + if (t < this->Knots.Get(0) || t > this->Knots.Get(n - 1)) + return -1; + + //Binary search for the interval + vtkm::Id left = 0; + vtkm::Id right = n - 1; + + while (left < right) + { + vtkm::Id mid = left + (right - left) / 2; + + if (t >= this->Knots.Get(mid) && t <= this->Knots.Get(mid + 1)) + return mid; + else if (t < this->Knots.Get(mid)) + right = mid; + else + left = mid; + } + + // t not within the interval. We should not get here. + return -1; + } + + DataPortalType Data; + KnotsPortalType Knots; + TangentPortalType Tangents; +}; +} //namespace exec +} //namespace vtkm + +#endif //vtk_m_exec_CubicHermiteSpline_h diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h new file mode 100644 index 0000000000000000000000000000000000000000..6a5c861a04d5ddff20d42456593da8729dee2100 --- /dev/null +++ b/vtkm/exec/CubicSpline.h @@ -0,0 +1,99 @@ +//============================================================================ +// 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 vtkm_exec_cubicspline_h +#define vtkm_exec_cubicspline_h + +#include +#include +#include + +namespace vtkm +{ + +namespace exec +{ +class VTKM_ALWAYS_EXPORT CubicSpline +{ +private: + using ArrayHandleType = vtkm::cont::ArrayHandle; + using PortalType = typename ArrayHandleType::ReadPortalType; + +public: + VTKM_CONT CubicSpline(const vtkm::cont::ArrayHandle& controlPoints, + const vtkm::cont::ArrayHandle& values, + const vtkm::cont::ArrayHandle& coefficientsB, + const vtkm::cont::ArrayHandle& coefficientsC, + const vtkm::cont::ArrayHandle& coefficientsD, + vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) + { + this->ControlPointsPortal = controlPoints.PrepareForInput(device, token); + this->ValuesPortal = values.PrepareForInput(device, token); + this->CoefficientsBPortal = coefficientsB.PrepareForInput(device, token); + this->CoefficientsCPortal = coefficientsC.PrepareForInput(device, token); + this->CoefficientsDPortal = coefficientsD.PrepareForInput(device, token); + } + + VTKM_EXEC + vtkm::ErrorCode Evaluate(const vtkm::FloatDefault& param, vtkm::FloatDefault& val) const + { + val = 0; + + vtkm::Id idx = this->FindInterval(param); + if (idx < 0) + return vtkm::ErrorCode::ValueOutOfRange; + + vtkm::FloatDefault dx = param - this->ControlPointsPortal.Get(idx); + auto B = this->CoefficientsBPortal.Get(idx); + auto C = this->CoefficientsCPortal.Get(idx); + auto D = this->CoefficientsDPortal.Get(idx); + val = this->ValuesPortal.Get(idx) + B * dx + C * dx * dx + D * dx * dx * dx; + + return vtkm::ErrorCode::Success; + } + +private: + vtkm::Id FindInterval(const vtkm::FloatDefault& x) const + { + vtkm::Id numPoints = this->ControlPointsPortal.GetNumberOfValues(); + + if (x < this->ControlPointsPortal.Get(0) || x > this->ControlPointsPortal.Get(numPoints - 1)) + return -1; + + //Binary search for the interval + vtkm::Id left = 0; + vtkm::Id right = numPoints - 1; + + while (left < right) + { + vtkm::Id mid = left + (right - left) / 2; + + if (x >= this->ControlPointsPortal.Get(mid) && x <= this->ControlPointsPortal.Get(mid + 1)) + return mid; + else if (x < this->ControlPointsPortal.Get(mid)) + right = mid; + else + left = mid; + } + + // x not within the interval. We should not get here. + return -1; + } + + PortalType ControlPointsPortal; + PortalType ValuesPortal; + PortalType CoefficientsBPortal; + PortalType CoefficientsCPortal; + PortalType CoefficientsDPortal; +}; +} //namespace exec +} //namespace vtkm + +#endif //vtkm_exec_cubicspline_h