From 8999ba0926cac27c22fd55fb4cd010dee85ac0d8 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 18 Mar 2025 08:48:03 -0400 Subject: [PATCH 01/13] 1d cubic splines. --- vtkm/cont/CMakeLists.txt | 2 + vtkm/cont/CubicSpline.cxx | 144 ++++++++++++++++++++++ vtkm/cont/CubicSpline.h | 72 +++++++++++ vtkm/cont/testing/CMakeLists.txt | 1 + vtkm/cont/testing/UnitTestCubicSpline.cxx | 103 ++++++++++++++++ vtkm/exec/CMakeLists.txt | 1 + vtkm/exec/CubicSpline.h | 87 +++++++++++++ 7 files changed, 410 insertions(+) create mode 100644 vtkm/cont/CubicSpline.cxx create mode 100644 vtkm/cont/CubicSpline.h create mode 100644 vtkm/cont/testing/UnitTestCubicSpline.cxx create mode 100644 vtkm/exec/CubicSpline.h diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index ca47f3cbb6..381abad219 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -80,6 +80,7 @@ set(headers ColorTableSamples.h ConvertNumComponentsToOffsets.h CoordinateSystem.h + CubicSpline.h DataSet.h DataSetBuilderCurvilinear.h DataSetBuilderExplicit.h @@ -153,6 +154,7 @@ set(sources CellSetStructured.cxx ColorTablePresets.cxx CoordinateSystem.cxx + CubicSpline.cxx DataSet.cxx DataSetBuilderCurvilinear.cxx DataSetBuilderExplicit.cxx diff --git a/vtkm/cont/CubicSpline.cxx b/vtkm/cont/CubicSpline.cxx new file mode 100644 index 0000000000..fd4566b040 --- /dev/null +++ b/vtkm/cont/CubicSpline.cxx @@ -0,0 +1,144 @@ +//============================================================================ +// 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 + +namespace vtkm +{ +namespace cont +{ + +void CubicSpline::Update() const +{ + if (this->Modified) + { + // Although the data of the derived class may change, the logical state + // of the class should not. Thus, we will instruct the compiler to relax + // the const constraint. + const_cast(this)->Build(); + this->Modified = false; + } +} + +void CubicSpline::Build() +{ + //Calculate the spline coeficients. + vtkm::Id n = this->ControlPoints.GetNumberOfValues(); + if (n < 2) + throw std::invalid_argument("At least two points are required for spline interpolation."); + + // Allocate memory for spline coefficients + this->CoefficientsB.Allocate(n - 1); + this->CoefficientsC.Allocate(n); + this->CoefficientsD.Allocate(n - 1); + //a1D = y1D; + //b1D.resize(n - 1); + //c1D.resize(n); + //d1D.resize(n - 1); + + //std::vector h(n - 1), alpha(n - 2); // Fix: correct alpha size + vtkm::cont::ArrayHandle h, alpha; + h.Allocate(n - 1); + alpha.Allocate(n - 2); + + + for (vtkm::Id i = 0; i < n - 1; i++) + { + //h[i] = x1D[i + 1] - x1D[i]; + h.WritePortal().Set( + i, this->ControlPoints.ReadPortal().Get(i + 1) - this->ControlPoints.ReadPortal().Get(i)); + } + + for (vtkm::Id i = 1; i < n - 1; ++i) + { + //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); + alpha.WritePortal().Set( + i - 1, + (3.0 * (this->Values.ReadPortal().Get(i + 1) - this->Values.ReadPortal().Get(i)) / + h.ReadPortal().Get(i)) - + (3.0 * (this->Values.ReadPortal().Get(i) - this->Values.ReadPortal().Get(i - 1)) / + h.ReadPortal().Get(i - 1))); + } + + //std::vector l(n), mu(n), z(n); + //l[0] = 1.0; + //mu[0] = z[0] = 0.0; + vtkm::cont::ArrayHandle l, mu, z; + l.Allocate(n); + mu.Allocate(n); + z.Allocate(n); + l.WritePortal().Set(0, 1.0); + mu.WritePortal().Set(0, 0.0); + z.WritePortal().Set(0, 0.0); + + + for (vtkm::Id i = 1; i < n - 1; ++i) + { + //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; + l.WritePortal().Set(i, + 2.0 * + (this->ControlPoints.ReadPortal().Get(i + 1) - + this->ControlPoints.ReadPortal().Get(i - 1)) - + h.ReadPortal().Get(i - 1) * mu.ReadPortal().Get(i - 1)); + //mu[i] = h[i] / l[i]; + mu.WritePortal().Set(i, h.ReadPortal().Get(i) / l.ReadPortal().Get(i)); + //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; + z.WritePortal().Set( + i, + (alpha.ReadPortal().Get(i - 1) - h.ReadPortal().Get(i - 1) * z.ReadPortal().Get(i - 1)) / + l.ReadPortal().Get(i)); + } + + //l[n - 1] = 1.0; + //z[n - 1] = c1D[n - 1] = 0.0; + l.WritePortal().Set(n - 1, 1.0); + z.WritePortal().Set(n - 1, 0.0); + this->CoefficientsC.WritePortal().Set(n - 1, 0.0); + + for (vtkm::Id j = n - 2; j >= 0; j--) + { + //c1D[j] = z[j] - mu[j] * c1D[j + 1]; + this->CoefficientsC.WritePortal().Set( + j, + z.ReadPortal().Get(j) - mu.ReadPortal().Get(j) * this->CoefficientsC.ReadPortal().Get(j + 1)); + //b1D[j] = (a1D[j + 1] - a1D[j]) / h[j] - h[j] * (c1D[j + 1] + 2.0 * c1D[j]) / 3.0; + auto val0 = (this->Values.ReadPortal().Get(j + 1) - this->Values.ReadPortal().Get(j)) / + h.ReadPortal().Get(j); + auto val1 = h.ReadPortal().Get(j) * + (this->CoefficientsC.ReadPortal().Get(j + 1) + + 2.0 * this->CoefficientsC.ReadPortal().Get(j)) / + 3.0; + this->CoefficientsB.WritePortal().Set(j, val0 - val1); + //d1D[j] = (c1D[j + 1] - c1D[j]) / (3.0 * h[j]); + this->CoefficientsD.WritePortal().Set( + j, + (this->CoefficientsC.ReadPortal().Get(j + 1) - this->CoefficientsC.ReadPortal().Get(j)) / + (3.0 * h.ReadPortal().Get(j))); + } +} + +vtkm::exec::CubicSpline CubicSpline::PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const +{ + this->Update(); + + using ExecObjType = vtkm::exec::CubicSpline; + + return ExecObjType(this->GetControlPoints(), + this->GetValues(), + this->CoefficientsB, + this->CoefficientsC, + this->CoefficientsD, + device, + token); +} + +} +} // namespace vtkm::cont diff --git a/vtkm/cont/CubicSpline.h b/vtkm/cont/CubicSpline.h new file mode 100644 index 0000000000..80a517ec84 --- /dev/null +++ b/vtkm/cont/CubicSpline.h @@ -0,0 +1,72 @@ +//============================================================================ +// 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_CubicSpline_h +#define vtk_m_cont_CubicSpline_h + +#include + +#include +#include +#include +#include + +namespace vtkm +{ +namespace cont +{ + +class VTKM_CONT_EXPORT CubicSpline : public vtkm::cont::ExecutionObjectBase +{ +public: + CubicSpline() = default; + virtual ~CubicSpline() = default; + + VTKM_CONT void SetControlPoints(const vtkm::cont::ArrayHandle& controlPoints) + { + this->ControlPoints = controlPoints; + this->SetModified(); + } + VTKM_CONT void SetValues(const vtkm::cont::ArrayHandle& values) + { + this->Values = values; + this->SetModified(); + } + + VTKM_CONT vtkm::cont::ArrayHandle GetControlPoints() const + { + return this->ControlPoints; + } + + VTKM_CONT vtkm::cont::ArrayHandle GetValues() const { return this->Values; } + + VTKM_CONT void Update() const; + + VTKM_CONT vtkm::exec::CubicSpline PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const; + +private: + void SetModified() { this->Modified = true; } + bool GetModified() const { return this->Modified; } + + void Build(); + + vtkm::cont::ArrayHandle ControlPoints; + vtkm::cont::ArrayHandle Values; + vtkm::cont::ArrayHandle CoefficientsA; + vtkm::cont::ArrayHandle CoefficientsB; + vtkm::cont::ArrayHandle CoefficientsC; + vtkm::cont::ArrayHandle CoefficientsD; + mutable bool Modified = true; +}; + +} +} // vtkm::cont + +#endif //vtk_m_cont_CubicSpline_h diff --git a/vtkm/cont/testing/CMakeLists.txt b/vtkm/cont/testing/CMakeLists.txt index 250a5f787a..10431a2e3e 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 + UnitTestCubicSpline.cxx UnitTestDataSetPermutation.cxx UnitTestDataSetSingleType.cxx UnitTestDeviceAdapterAlgorithmDependency.cxx diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx new file mode 100644 index 0000000000..fcbc2848d5 --- /dev/null +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -0,0 +1,103 @@ +//============================================================================ +// 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, + FieldOut valid); + using ExecutionSignature = void(_1, _2, _3, _4); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const ParamType& param, + const CubicSplineType& spline, + ResultType& value, + bool& valid) const + { + valid = spline.Evaluate(param, value); + } + +private: +}; + +void CheckEvaluation(const vtkm::cont::CubicSpline& spline, + const std::vector& params, + const std::vector& answer) +{ + auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); + + vtkm::cont::Invoker invoke; + SplineEvalWorklet worklet; + vtkm::cont::ArrayHandle result; + vtkm::cont::ArrayHandle valid; + invoke(worklet, paramsAH, spline, result, valid); + + VTKM_TEST_ASSERT(result.GetNumberOfValues() == static_cast(answer.size()), + "Result wrong length."); + return; + + for (std::size_t i = 0; i < answer.size(); i++) + { + VTKM_TEST_ASSERT(valid.ReadPortal().Get(static_cast(i)), "Evaluation failed."); + auto val = result.ReadPortal().Get(static_cast(i)); + std::cout << params[i] << " = " << val << " " << answer[i] + << " :: diff: " << vtkm::Abs(val - answer[i]) << std::endl; + VTKM_TEST_ASSERT(vtkm::Abs(val - answer[i]) < 1e-6, "Result has wrong value."); + } +} + +void CubicSplineTest() +{ + + std::vector xVals = { 0.0f, 1.0f, 1.5f, 2.0f, 3.0f, 4.0f }; + std::vector yVals = { 0.0f, -1.0f, 2.0f, 1.0f, -1.0f, 0.0f }; + + auto xValsAH = vtkm::cont::make_ArrayHandle(xVals, vtkm::CopyFlag::On); + auto yValsAH = vtkm::cont::make_ArrayHandle(yVals, vtkm::CopyFlag::On); + + vtkm::cont::CubicSpline cubicSpline; + cubicSpline.SetControlPoints(xValsAH); + cubicSpline.SetValues(yValsAH); + cubicSpline.Update(); + + //Ensure that the values at control points are properly interpolated. + auto tVals = xVals; + auto res = yVals; + CheckEvaluation(cubicSpline, tVals, res); + + //Evaluate between control points. + tVals = { 0.5, 1.3, 2.1, 3.7 }; + res = { -1.68413, 1.00957, 3.02636, -0.461989 }; + CheckEvaluation(cubicSpline, tVals, res); +} + +} // anonymous namespace + +int UnitTestCubicSpline(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(CubicSplineTest, argc, argv); +} diff --git a/vtkm/exec/CMakeLists.txt b/vtkm/exec/CMakeLists.txt index 2f27cb9077..a3283634d5 100644 --- a/vtkm/exec/CMakeLists.txt +++ b/vtkm/exec/CMakeLists.txt @@ -29,6 +29,7 @@ set(headers ConnectivityExtrude.h ConnectivityPermuted.h ConnectivityStructured.h + CubicSpline.h FieldNeighborhood.h FunctorBase.h ParametricCoordinates.h diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h new file mode 100644 index 0000000000..0aff21a942 --- /dev/null +++ b/vtkm/exec/CubicSpline.h @@ -0,0 +1,87 @@ +//============================================================================ +// 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 + +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 + bool Evaluate(const vtkm::FloatDefault& param, vtkm::FloatDefault& val) const + { + val = 0; + + vtkm::Id idx = this->FindInterval(param); + if (idx < 0) + return false; + + 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 true; + } + +private: + vtkm::Id FindInterval(const vtkm::FloatDefault& x) const + { + vtkm::Id numPoints = this->ControlPointsPortal.GetNumberOfValues(); + auto p0 = this->ControlPointsPortal.Get(0); + auto p1 = this->ControlPointsPortal.Get(numPoints - 1); + if (x < this->ControlPointsPortal.Get(0) || x > this->ControlPointsPortal.Get(numPoints - 1)) + return -1; + + for (vtkm::Id i = 0; i < numPoints - 1; ++i) + { + if (x >= this->ControlPointsPortal.Get(i) && x <= this->ControlPointsPortal.Get(i + 1)) + return i; + } + return -1; + } + + PortalType ControlPointsPortal; + PortalType ValuesPortal; + PortalType CoefficientsBPortal; + PortalType CoefficientsCPortal; + PortalType CoefficientsDPortal; +}; +} //namespace exec +} //namespace vtkm + +#endif //vtkm_exec_cubicspline_h -- GitLab From 09f2cb1bf5fa752a55ebbfc46503cd67381dff18 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 18 Mar 2025 09:06:36 -0400 Subject: [PATCH 02/13] Add comment to worklet-ize the coeficient building. --- vtkm/cont/CubicSpline.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/vtkm/cont/CubicSpline.cxx b/vtkm/cont/CubicSpline.cxx index fd4566b040..8821c93ac0 100644 --- a/vtkm/cont/CubicSpline.cxx +++ b/vtkm/cont/CubicSpline.cxx @@ -29,6 +29,9 @@ void CubicSpline::Update() const void CubicSpline::Build() { + //TODO: convert this a worklet. + + //Calculate the spline coeficients. vtkm::Id n = this->ControlPoints.GetNumberOfValues(); if (n < 2) -- GitLab From 148ce759a76d7fd03c9a5a930f4452c28a528c17 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Thu, 20 Mar 2025 10:52:18 -0400 Subject: [PATCH 03/13] precalculate the coefficients. --- vtkm/cont/CubicSpline.cxx | 259 +++++++++++++++------- vtkm/cont/CubicSpline.h | 9 + vtkm/cont/testing/UnitTestCubicSpline.cxx | 1 - vtkm/exec/CubicSpline.h | 8 +- 4 files changed, 187 insertions(+), 90 deletions(-) diff --git a/vtkm/cont/CubicSpline.cxx b/vtkm/cont/CubicSpline.cxx index 8821c93ac0..f793232562 100644 --- a/vtkm/cont/CubicSpline.cxx +++ b/vtkm/cont/CubicSpline.cxx @@ -9,11 +9,172 @@ //============================================================================ #include +#include namespace vtkm { namespace cont { +namespace +{ +class CalculateH : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldInOut h, WholeArrayIn xVals); + using ExecutionSignature = void(InputIndex, _1, _2); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, ResultType& h, const InputType& xVals) const + { + //h[i] = x1D[i + 1] - x1D[i]; + h = xVals.Get(idx + 1) - xVals.Get(idx); + } +}; + +class CalculateAlpha : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(FieldInOut alpha, WholeArrayIn h, WholeArrayIn yVals); + using ExecutionSignature = void(InputIndex, _1, _2, _3); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, + ResultType& alpha, + const HArrayType& h, + const YArrayType& y) const + { + //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); + // alpha: i-1 --> i. Add +1 to each i. + //alpha[i] = (3.0 * (a1D[i + 2] - a1D[i+1]) / h[i+1]) - (3.0 * (a1D[i+1] - a1D[i]) / h[i]); + + alpha = (3.0 * (y.Get(idx + 2) - y.Get(idx + 1)) / h.Get(idx + 1)) - + (3.0 * (y.Get(idx + 1) - y.Get(idx)) / h.Get(idx)); + } +}; + +class CalculateLMuZ : public vtkm::worklet::WorkletMapField +{ +public: + using ControlSignature = void(WholeArrayOut L, + WholeArrayOut Mu, + WholeArrayOut Z, + WholeArrayIn xVals, + WholeArrayIn alpha, + WholeArrayIn h); + using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6); + using InputDomain = _1; + + CalculateLMuZ(const vtkm::Id& size) + : Size(size) + { + } + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, + LType& L, + MuType& Mu, + ZType& Z, + const XType& xVal, + const AlphaType& alpha, + const HType& H) const + { + if (idx == 0) + { + L.Set(idx, 1.0); + Mu.Set(idx, 0.0); + Z.Set(idx, 0.0); + } + else if (idx == this->Size - 1) + { + L.Set(idx, 1.0); + Z.Set(idx, 1.0); + } + else + { + //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; + L.Set(idx, 2.0 * (xVal.Get(idx + 1) - xVal.Get(idx - 1)) - H.Get(idx - 1) * Mu.Get(idx - 1)); + + //mu[i] = h[i] / l[i]; + Mu.Set(idx, H.Get(idx) / L.Get(idx)); + + //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; + Z.Set(idx, alpha.Get(idx - 1) - H.Get(idx - 1) * Z.Get(idx - 1) / L.Get(idx)); + } + } + +private: + vtkm::Id Size; +}; + +} //namespace + +template +void CubicSpline::CalculateLMuZ(ArrayType& mu, + ArrayType& z, + const ArrayType& alpha, + const ArrayType& h) const +{ + auto Mu = mu.WritePortal(); + auto Z = z.WritePortal(); + auto Alpha = alpha.ReadPortal(); + const auto H = h.ReadPortal(); + const auto xVals = this->ControlPoints.ReadPortal(); + vtkm::Id size = h.GetNumberOfValues(); + + // initial values. + //L.Set(0, 1.0); + Mu.Set(0, 0.0); + Z.Set(0, 0.0); + + for (vtkm::Id i = 1; i < size - 1; i++) + { + auto L = 2.0 * (xVals.Get(i + 1) - xVals.Get(i - 1)) - H.Get(i - 1) * Mu.Get(i - 1); + Mu.Set(i, H.Get(i) / L); + Z.Set(i, (Alpha.Get(i - 1) - H.Get(i - 1) * Z.Get(i - 1)) / L); + } + //L.Set(size - 1, 1.0); + Z.Set(size - 1, 0.0); +} + +template +void CubicSpline::CalculateCoefficients(vtkm::Id n, + const ArrayPortalType& H, + const ArrayPortalType& Mu, + const ArrayPortalType& Z) +{ + this->CoefficientsB.Allocate(n - 1); + this->CoefficientsC.Allocate(n); + this->CoefficientsD.Allocate(n - 1); + auto B = this->CoefficientsB.WritePortal(); + auto C = this->CoefficientsC.WritePortal(); + auto D = this->CoefficientsD.WritePortal(); + const auto yVals = this->Values.ReadPortal(); + + C.Set(n - 1, 0.0); + + for (vtkm::Id i = n - 2; i >= 0; i--) + { + //c1D[j] = z[j] - mu[j] * c1D[j + 1]; + auto val = Z.Get(i) - Mu.Get(i) * C.Get(i + 1); + C.Set(i, val); + + //b1D[j] = (a1D[j + 1] - a1D[j]) / h[j] - h[j] * (c1D[j + 1] + 2.0 * c1D[j]) / 3.0; + val = (yVals.Get(i + 1) - yVals.Get(i)) / H.Get(i) - + H.Get(i) * (C.Get(i + 1) + 2.0 * C.Get(i)) / 3.0; + B.Set(i, val); + + //d1D[j] = (c1D[j + 1] - c1D[j]) / (3.0 * h[j]); + val = (C.Get(i + 1) - C.Get(i)) / (3.0 * H.Get(i)); + D.Set(i, val); + } +} void CubicSpline::Update() const { @@ -31,100 +192,30 @@ void CubicSpline::Build() { //TODO: convert this a worklet. - //Calculate the spline coeficients. vtkm::Id n = this->ControlPoints.GetNumberOfValues(); if (n < 2) throw std::invalid_argument("At least two points are required for spline interpolation."); - // Allocate memory for spline coefficients - this->CoefficientsB.Allocate(n - 1); - this->CoefficientsC.Allocate(n); - this->CoefficientsD.Allocate(n - 1); - //a1D = y1D; - //b1D.resize(n - 1); - //c1D.resize(n); - //d1D.resize(n - 1); - - //std::vector h(n - 1), alpha(n - 2); // Fix: correct alpha size - vtkm::cont::ArrayHandle h, alpha; + //calculate h + vtkm::cont::ArrayHandle h; h.Allocate(n - 1); - alpha.Allocate(n - 2); + vtkm::cont::Invoker invoker; + invoker(CalculateH{}, h, this->ControlPoints); + //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); + // for loop from 0 to n-2 + // alpha[i] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i + 1]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i]); + vtkm::cont::ArrayHandle alpha; + alpha.Allocate(n - 2); + invoker(CalculateAlpha{}, alpha, h, this->Values); - for (vtkm::Id i = 0; i < n - 1; i++) - { - //h[i] = x1D[i + 1] - x1D[i]; - h.WritePortal().Set( - i, this->ControlPoints.ReadPortal().Get(i + 1) - this->ControlPoints.ReadPortal().Get(i)); - } - - for (vtkm::Id i = 1; i < n - 1; ++i) - { - //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); - alpha.WritePortal().Set( - i - 1, - (3.0 * (this->Values.ReadPortal().Get(i + 1) - this->Values.ReadPortal().Get(i)) / - h.ReadPortal().Get(i)) - - (3.0 * (this->Values.ReadPortal().Get(i) - this->Values.ReadPortal().Get(i - 1)) / - h.ReadPortal().Get(i - 1))); - } - - //std::vector l(n), mu(n), z(n); - //l[0] = 1.0; - //mu[0] = z[0] = 0.0; - vtkm::cont::ArrayHandle l, mu, z; - l.Allocate(n); + //calculate mu, z. + vtkm::cont::ArrayHandle mu, z; mu.Allocate(n); z.Allocate(n); - l.WritePortal().Set(0, 1.0); - mu.WritePortal().Set(0, 0.0); - z.WritePortal().Set(0, 0.0); - - - for (vtkm::Id i = 1; i < n - 1; ++i) - { - //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; - l.WritePortal().Set(i, - 2.0 * - (this->ControlPoints.ReadPortal().Get(i + 1) - - this->ControlPoints.ReadPortal().Get(i - 1)) - - h.ReadPortal().Get(i - 1) * mu.ReadPortal().Get(i - 1)); - //mu[i] = h[i] / l[i]; - mu.WritePortal().Set(i, h.ReadPortal().Get(i) / l.ReadPortal().Get(i)); - //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; - z.WritePortal().Set( - i, - (alpha.ReadPortal().Get(i - 1) - h.ReadPortal().Get(i - 1) * z.ReadPortal().Get(i - 1)) / - l.ReadPortal().Get(i)); - } - - //l[n - 1] = 1.0; - //z[n - 1] = c1D[n - 1] = 0.0; - l.WritePortal().Set(n - 1, 1.0); - z.WritePortal().Set(n - 1, 0.0); - this->CoefficientsC.WritePortal().Set(n - 1, 0.0); - - for (vtkm::Id j = n - 2; j >= 0; j--) - { - //c1D[j] = z[j] - mu[j] * c1D[j + 1]; - this->CoefficientsC.WritePortal().Set( - j, - z.ReadPortal().Get(j) - mu.ReadPortal().Get(j) * this->CoefficientsC.ReadPortal().Get(j + 1)); - //b1D[j] = (a1D[j + 1] - a1D[j]) / h[j] - h[j] * (c1D[j + 1] + 2.0 * c1D[j]) / 3.0; - auto val0 = (this->Values.ReadPortal().Get(j + 1) - this->Values.ReadPortal().Get(j)) / - h.ReadPortal().Get(j); - auto val1 = h.ReadPortal().Get(j) * - (this->CoefficientsC.ReadPortal().Get(j + 1) + - 2.0 * this->CoefficientsC.ReadPortal().Get(j)) / - 3.0; - this->CoefficientsB.WritePortal().Set(j, val0 - val1); - //d1D[j] = (c1D[j + 1] - c1D[j]) / (3.0 * h[j]); - this->CoefficientsD.WritePortal().Set( - j, - (this->CoefficientsC.ReadPortal().Get(j + 1) - this->CoefficientsC.ReadPortal().Get(j)) / - (3.0 * h.ReadPortal().Get(j))); - } + this->CalculateLMuZ(mu, z, alpha, h); + this->CalculateCoefficients(n, h.ReadPortal(), mu.ReadPortal(), z.ReadPortal()); } vtkm::exec::CubicSpline CubicSpline::PrepareForExecution(vtkm::cont::DeviceAdapterId device, diff --git a/vtkm/cont/CubicSpline.h b/vtkm/cont/CubicSpline.h index 80a517ec84..9fe988680b 100644 --- a/vtkm/cont/CubicSpline.h +++ b/vtkm/cont/CubicSpline.h @@ -55,6 +55,15 @@ private: void SetModified() { this->Modified = true; } bool GetModified() const { return this->Modified; } + template + void CalculateLMuZ(ArrayType& Mu, ArrayType& Z, const ArrayType& alpha, const ArrayType& h) const; + + template + void CalculateCoefficients(vtkm::Id n, + const ArrayPortalType& H, + const ArrayPortalType& Mu, + const ArrayPortalType& Z); + void Build(); vtkm::cont::ArrayHandle ControlPoints; diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx index fcbc2848d5..5572f379f3 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -72,7 +72,6 @@ void CheckEvaluation(const vtkm::cont::CubicSpline& spline, void CubicSplineTest() { - std::vector xVals = { 0.0f, 1.0f, 1.5f, 2.0f, 3.0f, 4.0f }; std::vector yVals = { 0.0f, -1.0f, 2.0f, 1.0f, -1.0f, 0.0f }; diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h index 0aff21a942..f21ed0118f 100644 --- a/vtkm/exec/CubicSpline.h +++ b/vtkm/exec/CubicSpline.h @@ -50,9 +50,9 @@ public: return false; 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); + const auto B = this->CoefficientsBPortal.Get(idx); + const auto C = this->CoefficientsCPortal.Get(idx); + const auto D = this->CoefficientsDPortal.Get(idx); val = this->ValuesPortal.Get(idx) + B * dx + C * dx * dx + D * dx * dx * dx; return true; @@ -62,8 +62,6 @@ private: vtkm::Id FindInterval(const vtkm::FloatDefault& x) const { vtkm::Id numPoints = this->ControlPointsPortal.GetNumberOfValues(); - auto p0 = this->ControlPointsPortal.Get(0); - auto p1 = this->ControlPointsPortal.Get(numPoints - 1); if (x < this->ControlPointsPortal.Get(0) || x > this->ControlPointsPortal.Get(numPoints - 1)) return -1; -- GitLab From 99e147744619026ba7340a3572532b953410eccb Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Mon, 24 Mar 2025 11:30:01 -0400 Subject: [PATCH 04/13] Code cleanup. --- vtkm/cont/CubicSpline.cxx | 237 +++++++------------ vtkm/cont/CubicSpline.h | 15 +- vtkm/cont/testing/UnitTestCubicSpline.cxx | 271 +++++++++++++++++++++- vtkm/exec/CubicSpline.h | 7 +- 4 files changed, 356 insertions(+), 174 deletions(-) diff --git a/vtkm/cont/CubicSpline.cxx b/vtkm/cont/CubicSpline.cxx index f793232562..77226d4f89 100644 --- a/vtkm/cont/CubicSpline.cxx +++ b/vtkm/cont/CubicSpline.cxx @@ -15,207 +15,128 @@ namespace vtkm { namespace cont { -namespace -{ -class CalculateH : public vtkm::worklet::WorkletMapField -{ -public: - using ControlSignature = void(FieldInOut h, WholeArrayIn xVals); - using ExecutionSignature = void(InputIndex, _1, _2); - using InputDomain = _1; - template - VTKM_EXEC void operator()(const vtkm::Id& idx, ResultType& h, const InputType& xVals) const +void CubicSpline::Update() const +{ + if (this->Modified) { - //h[i] = x1D[i + 1] - x1D[i]; - h = xVals.Get(idx + 1) - xVals.Get(idx); + // Although the data of the derived class may change, the logical state + // of the class should not. Thus, we will instruct the compiler to relax + // the const constraint. + const_cast(this)->Build(); + this->Modified = false; } -}; +} -class CalculateAlpha : public vtkm::worklet::WorkletMapField +vtkm::cont::ArrayHandle CubicSpline::CalcH() const { -public: - using ControlSignature = void(FieldInOut alpha, WholeArrayIn h, WholeArrayIn yVals); - using ExecutionSignature = void(InputIndex, _1, _2, _3); - using InputDomain = _1; - - template - VTKM_EXEC void operator()(const vtkm::Id& idx, - ResultType& alpha, - const HArrayType& h, - const YArrayType& y) const - { - //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); - // alpha: i-1 --> i. Add +1 to each i. - //alpha[i] = (3.0 * (a1D[i + 2] - a1D[i+1]) / h[i+1]) - (3.0 * (a1D[i+1] - a1D[i]) / h[i]); + vtkm::Id n = this->ControlPoints.GetNumberOfValues(); + vtkm::cont::ArrayHandle h; + h.Allocate(n - 1); - alpha = (3.0 * (y.Get(idx + 2) - y.Get(idx + 1)) / h.Get(idx + 1)) - - (3.0 * (y.Get(idx + 1) - y.Get(idx)) / h.Get(idx)); - } -}; + auto hPortal = h.WritePortal(); + auto cPortal = this->ControlPoints.ReadPortal(); + for (vtkm::Id i = 0; i < n - 1; i++) + hPortal.Set(i, cPortal.Get(i + 1) - cPortal.Get(i)); -class CalculateLMuZ : public vtkm::worklet::WorkletMapField + return h; +} + +vtkm::cont::ArrayHandle CubicSpline::CalcAlpha( + const vtkm::cont::ArrayHandle& h) const { -public: - using ControlSignature = void(WholeArrayOut L, - WholeArrayOut Mu, - WholeArrayOut Z, - WholeArrayIn xVals, - WholeArrayIn alpha, - WholeArrayIn h); - using ExecutionSignature = void(InputIndex, _1, _2, _3, _4, _5, _6); - using InputDomain = _1; - - CalculateLMuZ(const vtkm::Id& size) - : Size(size) - { - } + vtkm::cont::ArrayHandle alpha; + vtkm::Id n = this->ControlPoints.GetNumberOfValues(); + alpha.Allocate(n - 2); + + auto Alpha = alpha.WritePortal(); + auto H = h.ReadPortal(); + auto Y = this->Values.ReadPortal(); - template - VTKM_EXEC void operator()(const vtkm::Id& idx, - LType& L, - MuType& Mu, - ZType& Z, - const XType& xVal, - const AlphaType& alpha, - const HType& H) const + for (vtkm::Id i = 1; i < n - 1; ++i) { - if (idx == 0) - { - L.Set(idx, 1.0); - Mu.Set(idx, 0.0); - Z.Set(idx, 0.0); - } - else if (idx == this->Size - 1) - { - L.Set(idx, 1.0); - Z.Set(idx, 1.0); - } - else - { - //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; - L.Set(idx, 2.0 * (xVal.Get(idx + 1) - xVal.Get(idx - 1)) - H.Get(idx - 1) * Mu.Get(idx - 1)); - - //mu[i] = h[i] / l[i]; - Mu.Set(idx, H.Get(idx) / L.Get(idx)); - - //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; - Z.Set(idx, alpha.Get(idx - 1) - H.Get(idx - 1) * Z.Get(idx - 1) / L.Get(idx)); - } + //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); + Alpha.Set(i - 1, + (3.0 * (Y.Get(i + 1) - Y.Get(i)) / H.Get(i)) - + (3.0 * (Y.Get(i) - Y.Get(i - 1)) / H.Get(i - 1))); } -private: - vtkm::Id Size; -}; - -} //namespace + return alpha; +} -template -void CubicSpline::CalculateLMuZ(ArrayType& mu, - ArrayType& z, - const ArrayType& alpha, - const ArrayType& h) const +void CubicSpline::CalcCoefficients(const vtkm::cont::ArrayHandle& H, + const vtkm::cont::ArrayHandle& Alpha) { - auto Mu = mu.WritePortal(); - auto Z = z.WritePortal(); - auto Alpha = alpha.ReadPortal(); - const auto H = h.ReadPortal(); - const auto xVals = this->ControlPoints.ReadPortal(); - vtkm::Id size = h.GetNumberOfValues(); - - // initial values. - //L.Set(0, 1.0); - Mu.Set(0, 0.0); - Z.Set(0, 0.0); - - for (vtkm::Id i = 1; i < size - 1; i++) + vtkm::Id n = this->ControlPoints.GetNumberOfValues(); + vtkm::cont::ArrayHandle L, Mu, Z; + + L.Allocate(n); + Mu.Allocate(n); + Z.Allocate(n); + + auto l = L.WritePortal(); + auto mu = Mu.WritePortal(); + auto z = Z.WritePortal(); + auto h = H.ReadPortal(); + auto alpha = Alpha.ReadPortal(); + + l.Set(0, 1.0); + mu.Set(0, 0.0); + z.Set(0, 0.0); + + auto xVals = this->ControlPoints.ReadPortal(); + + for (vtkm::Id i = 1; i < n - 1; ++i) { - auto L = 2.0 * (xVals.Get(i + 1) - xVals.Get(i - 1)) - H.Get(i - 1) * Mu.Get(i - 1); - Mu.Set(i, H.Get(i) / L); - Z.Set(i, (Alpha.Get(i - 1) - H.Get(i - 1) * Z.Get(i - 1)) / L); + //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; + l.Set(i, 2.0 * (xVals.Get(i + 1) - xVals.Get(i - 1)) - h.Get(i - 1) * mu.Get(i - 1)); + //mu[i] = h[i] / l[i]; + mu.Set(i, h.Get(i) / l.Get(i)); + //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; + z.Set(i, (alpha.Get(i - 1) - h.Get(i - 1) * z.Get(i - 1)) / l.Get(i)); } - //L.Set(size - 1, 1.0); - Z.Set(size - 1, 0.0); -} -template -void CubicSpline::CalculateCoefficients(vtkm::Id n, - const ArrayPortalType& H, - const ArrayPortalType& Mu, - const ArrayPortalType& Z) -{ + //l[n - 1] = 1.0; + //z[n - 1] = c1D[n - 1] = 0.0; + l.Set(n - 1, 1.0); + z.Set(n - 1, 0.0); + this->CoefficientsB.Allocate(n - 1); this->CoefficientsC.Allocate(n); this->CoefficientsD.Allocate(n - 1); + auto B = this->CoefficientsB.WritePortal(); auto C = this->CoefficientsC.WritePortal(); auto D = this->CoefficientsD.WritePortal(); - const auto yVals = this->Values.ReadPortal(); C.Set(n - 1, 0.0); + auto yVals = this->Values.ReadPortal(); for (vtkm::Id i = n - 2; i >= 0; i--) { //c1D[j] = z[j] - mu[j] * c1D[j + 1]; - auto val = Z.Get(i) - Mu.Get(i) * C.Get(i + 1); - C.Set(i, val); - + C.Set(i, z.Get(i) - mu.Get(i) * C.Get(i + 1)); //b1D[j] = (a1D[j + 1] - a1D[j]) / h[j] - h[j] * (c1D[j + 1] + 2.0 * c1D[j]) / 3.0; - val = (yVals.Get(i + 1) - yVals.Get(i)) / H.Get(i) - - H.Get(i) * (C.Get(i + 1) + 2.0 * C.Get(i)) / 3.0; - B.Set(i, val); - + auto val0 = (yVals.Get(i + 1) - yVals.Get(i)) / h.Get(i); + auto val1 = h.Get(i) * (C.Get(i + 1) + 2.0 * C.Get(i)) / 3.0; + B.Set(i, val0 - val1); //d1D[j] = (c1D[j + 1] - c1D[j]) / (3.0 * h[j]); - val = (C.Get(i + 1) - C.Get(i)) / (3.0 * H.Get(i)); - D.Set(i, val); + D.Set(i, (C.Get(i + 1) - C.Get(i)) / (3.0 * h.Get(i))); } } -void CubicSpline::Update() const -{ - if (this->Modified) - { - // Although the data of the derived class may change, the logical state - // of the class should not. Thus, we will instruct the compiler to relax - // the const constraint. - const_cast(this)->Build(); - this->Modified = false; - } -} void CubicSpline::Build() { - //TODO: convert this a worklet. - //Calculate the spline coeficients. vtkm::Id n = this->ControlPoints.GetNumberOfValues(); if (n < 2) throw std::invalid_argument("At least two points are required for spline interpolation."); - //calculate h - vtkm::cont::ArrayHandle h; - h.Allocate(n - 1); - vtkm::cont::Invoker invoker; - invoker(CalculateH{}, h, this->ControlPoints); + auto h = this->CalcH(); + auto alpha = this->CalcAlpha(h); - //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); - // for loop from 0 to n-2 - // alpha[i] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i + 1]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i]); - vtkm::cont::ArrayHandle alpha; - alpha.Allocate(n - 2); - invoker(CalculateAlpha{}, alpha, h, this->Values); - - //calculate mu, z. - vtkm::cont::ArrayHandle mu, z; - mu.Allocate(n); - z.Allocate(n); - this->CalculateLMuZ(mu, z, alpha, h); - this->CalculateCoefficients(n, h.ReadPortal(), mu.ReadPortal(), z.ReadPortal()); + this->CalcCoefficients(h, alpha); } vtkm::exec::CubicSpline CubicSpline::PrepareForExecution(vtkm::cont::DeviceAdapterId device, diff --git a/vtkm/cont/CubicSpline.h b/vtkm/cont/CubicSpline.h index 9fe988680b..c3976e40df 100644 --- a/vtkm/cont/CubicSpline.h +++ b/vtkm/cont/CubicSpline.h @@ -55,17 +55,14 @@ private: void SetModified() { this->Modified = true; } bool GetModified() const { return this->Modified; } - template - void CalculateLMuZ(ArrayType& Mu, ArrayType& Z, const ArrayType& alpha, const ArrayType& h) const; - - template - void CalculateCoefficients(vtkm::Id n, - const ArrayPortalType& H, - const ArrayPortalType& Mu, - const ArrayPortalType& Z); - void Build(); + vtkm::cont::ArrayHandle CalcH() const; + vtkm::cont::ArrayHandle CalcAlpha( + const vtkm::cont::ArrayHandle& h) const; + void CalcCoefficients(const vtkm::cont::ArrayHandle& H, + const vtkm::cont::ArrayHandle& Alpha); + vtkm::cont::ArrayHandle ControlPoints; vtkm::cont::ArrayHandle Values; vtkm::cont::ArrayHandle CoefficientsA; diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx index 5572f379f3..bb357b8099 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -17,6 +17,215 @@ #include #include +class Spline +{ +public: + // Constructor: expects x and y to have the same size and at least 2 points. + Spline(const std::vector& x, const std::vector& y) + : x_(x) + , y_(y) + , n_(x.size()) + { + if (n_ < 2) + { + throw std::invalid_argument("At least two points are required for spline interpolation."); + } + if (x_.size() != y_.size()) + { + throw std::invalid_argument("x and y must have the same length."); + } + computeCoefficients(); + } + + // Evaluate the spline at a given x value. + vtkm::FloatDefault evaluate(vtkm::FloatDefault x_val) const + { + // Ensure x_val is within the interpolation range. + if (x_val < x_.front() || x_val > x_.back()) + { + throw std::out_of_range("x_val is outside the interpolation range."); + } + // Find the interval index i such that x_[i] <= x_val <= x_[i+1]. + size_t i = std::upper_bound(x_.begin(), x_.end(), x_val) - x_.begin() - 1; + vtkm::FloatDefault dx = x_val - x_[i]; + return a_[i] + b_[i] * dx + c_[i] * dx * dx + d_[i] * dx * dx * dx; + } + +private: + std::vector x_; // x coordinates + std::vector y_; // y coordinates + size_t n_; // number of points + + // Spline coefficients for each interval i = 0, ..., n-2. + // For each segment, the cubic polynomial is: + // S_i(x) = a_[i] + b_[i]*(x-x_[i]) + c_[i]*(x-x_[i])^2 + d_[i]*(x-x_[i])^3 + std::vector a_, b_, c_, d_; + + // Compute the spline coefficients using the standard algorithm. + void computeCoefficients() + { + // Natural cubic spline: second derivatives at endpoints are zero. + // Step 1: Compute h[i] = x[i+1] - x[i] for i = 0 ... n-2. + std::vector h(n_ - 1); + for (size_t i = 0; i < n_ - 1; ++i) + { + h[i] = x_[i + 1] - x_[i]; + } + + // Step 2: Compute the array 'alpha' for i = 1 ... n-2: + // alpha[i] = 3/h[i]*(y[i+1]-y[i]) - 3/h[i-1]*(y[i]-y[i-1]) + std::vector alpha(n_, 0.0); + for (size_t i = 1; i < n_ - 1; ++i) + { + alpha[i] = (3.0 / h[i]) * (y_[i + 1] - y_[i]) - (3.0 / h[i - 1]) * (y_[i] - y_[i - 1]); + } + + // Step 3: Set up the tridiagonal system and solve it. + std::vector l(n_), mu(n_), z(n_); + l[0] = 1.0; + mu[0] = z[0] = 0.0; + for (size_t i = 1; i < n_ - 1; ++i) + { + l[i] = 2.0 * (x_[i + 1] - x_[i - 1]) - h[i - 1] * mu[i - 1]; + mu[i] = h[i] / l[i]; + z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]; + } + l[n_ - 1] = 1.0; + z[n_ - 1] = 0.0; + c_.resize(n_); + c_[n_ - 1] = 0.0; + + // Step 4: Back substitution to solve for c[i]. + for (int j = n_ - 2; j >= 0; --j) + { + c_[j] = z[j] - mu[j] * c_[j + 1]; + } + + // Step 5: Compute the coefficients b and d for each interval. + a_.resize(n_ - 1); + b_.resize(n_ - 1); + d_.resize(n_ - 1); + for (size_t i = 0; i < n_ - 1; ++i) + { + a_[i] = y_[i]; // a[i] is simply the y value at x[i] + b_[i] = (y_[i + 1] - y_[i]) / h[i] - h[i] * (c_[i + 1] + 2.0 * c_[i]) / 3.0; + d_[i] = (c_[i + 1] - c_[i]) / (3.0 * h[i]); + } + } +}; + +class Spline2 +{ +public: + Spline2(const std::vector& x, const std::vector& y) + : x1D(x) + , y1D(y) + , n(x.size()) + { + if (x1D.size() != y1D.size()) + throw std::invalid_argument("x and y must have the same length."); + if (n < 2) + throw std::invalid_argument("At least two points are required for spline interpolation."); + } + + // Evaluate the spline at a given x value. + vtkm::FloatDefault evaluate(vtkm::FloatDefault x_val) const + { + // Compute the spline coefficients on the fly. + Coefs coefs = computeCoefficients(); + + // Ensure x_val is within the interpolation range. + if (x_val < x1D.front() || x_val > x1D.back()) + throw std::out_of_range("x_val is outside the interpolation range."); + + // Find the interval index i such that x1D[i] <= x_val <= x1D[i+1] + size_t i = findInterval(x1D, x_val); + + vtkm::FloatDefault dx = x_val - x1D[i]; + // Evaluate S_i(x) = a[i] + b[i]*(x - x_i) + c[i]*(x - x_i)^2 + d[i]*(x - x_i)^3 + return coefs.a[i] + coefs.b[i] * dx + coefs.c[i] * dx * dx + coefs.d[i] * dx * dx * dx; + } + +private: + std::vector x1D, y1D; + size_t n; + + // Structure to hold the spline coefficients for each interval. + struct Coefs + { + std::vector a; + std::vector b; + std::vector c; + std::vector d; + }; + + // Compute the cubic spline coefficients using the natural cubic spline algorithm. + // This function computes the full set of coefficients for all segments. + Coefs computeCoefficients() const + { + Coefs coefs; + coefs.a = y1D; // a[i] = y1D[i] + + // Allocate coefficient arrays for each segment (there are n-1 segments) + coefs.b.resize(n - 1); + coefs.d.resize(n - 1); + coefs.c.resize(n); // second derivatives for each point + + // Step 1: Compute h[i] = x1D[i+1] - x1D[i] + std::vector h(n - 1); + for (size_t i = 0; i < n - 1; ++i) + { + h[i] = x1D[i + 1] - x1D[i]; + } + + // Step 2: Compute alpha[i] for i=1..n-2 + std::vector alpha(n, 0.0); + for (size_t i = 1; i < n - 1; ++i) + { + alpha[i] = (3.0 / h[i]) * (y1D[i + 1] - y1D[i]) - (3.0 / h[i - 1]) * (y1D[i] - y1D[i - 1]); + } + + // Step 3: Set up and solve the tridiagonal system + std::vector l(n, 0.0), mu(n, 0.0), z(n, 0.0); + l[0] = 1.0; + mu[0] = 0.0; + z[0] = 0.0; + for (size_t i = 1; i < n - 1; ++i) + { + l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; + mu[i] = h[i] / l[i]; + z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]; + } + l[n - 1] = 1.0; + z[n - 1] = 0.0; + coefs.c[n - 1] = 0.0; + + // Step 4: Back substitution to solve for c[i] + for (int j = n - 2; j >= 0; --j) + { + coefs.c[j] = z[j] - mu[j] * coefs.c[j + 1]; + } + + // Step 5: Compute b[i] and d[i] for i=0..n-2 + for (size_t i = 0; i < n - 1; ++i) + { + coefs.b[i] = (y1D[i + 1] - y1D[i]) / h[i] - h[i] * (coefs.c[i + 1] + 2.0 * coefs.c[i]) / 3.0; + coefs.d[i] = (coefs.c[i + 1] - coefs.c[i]) / (3.0 * h[i]); + } + + return coefs; + } + + // Find the interval index for a given value. + size_t findInterval(const std::vector& x, vtkm::FloatDefault val) const + { + if (val < x.front() || val > x.back()) + throw std::out_of_range("Value is outside the interpolation range."); + size_t idx = std::upper_bound(x.begin(), x.end(), val) - x.begin() - 1; + return (idx >= x.size() - 1) ? x.size() - 2 : idx; + } +}; + namespace { @@ -40,8 +249,6 @@ public: { valid = spline.Evaluate(param, value); } - -private: }; void CheckEvaluation(const vtkm::cont::CubicSpline& spline, @@ -70,10 +277,61 @@ void CheckEvaluation(const vtkm::cont::CubicSpline& spline, } } +void SaveSamples(const vtkm::cont::CubicSpline& spline_) +{ + vtkm::Id n = spline_.GetControlPoints().GetNumberOfValues(); + vtkm::FloatDefault t = spline_.GetControlPoints().ReadPortal().Get(0); + vtkm::FloatDefault tEnd = spline_.GetControlPoints().ReadPortal().Get(n - 1); + + std::vector x1D, y1D; + for (vtkm::Id i = 0; i < n; i++) + { + x1D.push_back(spline_.GetControlPoints().ReadPortal().Get(i)); + y1D.push_back(spline_.GetValues().ReadPortal().Get(i)); + } + + //Spline spline(x1D, y1D); + Spline2 spline(x1D, y1D); + + vtkm::FloatDefault dt = (tEnd - t) / (n * 100); + std::vector params; + while (t < tEnd) + { + params.push_back(t); + t += dt; + } + + /* + auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); + + vtkm::cont::Invoker invoke; + SplineEvalWorklet worklet; + vtkm::cont::ArrayHandle result; + vtkm::cont::ArrayHandle valid; + invoke(worklet, paramsAH, spline, result, valid); +*/ + + std::cout << "SaveSamples" << std::endl; + std::ofstream fout("/Users/dpn/output.txt"), ptsOut("/Users/dpn/pts.txt"); + fout << "X,Y" << std::endl; + for (std::size_t i = 0; i < params.size(); i++) + { + fout << params[i] << "," << spline.evaluate(params[i]) << std::endl; + //fout << params[i] << "," << result.ReadPortal().Get(i) << std::endl; + } + + ptsOut << "X,Y" << std::endl; + for (std::size_t i = 0; i < x1D.size(); i++) + ptsOut << x1D[i] << ", " << y1D[i] << std::endl; +} + + void CubicSplineTest() { - std::vector xVals = { 0.0f, 1.0f, 1.5f, 2.0f, 3.0f, 4.0f }; - std::vector yVals = { 0.0f, -1.0f, 2.0f, 1.0f, -1.0f, 0.0f }; + //std::vector xVals = { 0.0f, 1.0f, 1.5f, 2.0f, 3.0f, 4.0f }; + //std::vector yVals = { 0.0f, -1.0f, 2.0f, 1.0f, -1.0f, 0.0f }; + std::vector xVals = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f }; + std::vector yVals = { 0.0f, 1.0f, 4.4f, 1.0f, 0.0f }; auto xValsAH = vtkm::cont::make_ArrayHandle(xVals, vtkm::CopyFlag::On); auto yValsAH = vtkm::cont::make_ArrayHandle(yVals, vtkm::CopyFlag::On); @@ -83,15 +341,20 @@ void CubicSplineTest() cubicSpline.SetValues(yValsAH); cubicSpline.Update(); + SaveSamples(cubicSpline); + + //Ensure that the values at control points are properly interpolated. auto tVals = xVals; auto res = yVals; CheckEvaluation(cubicSpline, tVals, res); + /* //Evaluate between control points. tVals = { 0.5, 1.3, 2.1, 3.7 }; res = { -1.68413, 1.00957, 3.02636, -0.461989 }; CheckEvaluation(cubicSpline, tVals, res); + */ } } // anonymous namespace diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h index f21ed0118f..29ea7d8d1d 100644 --- a/vtkm/exec/CubicSpline.h +++ b/vtkm/exec/CubicSpline.h @@ -50,9 +50,9 @@ public: return false; vtkm::FloatDefault dx = param - this->ControlPointsPortal.Get(idx); - const auto B = this->CoefficientsBPortal.Get(idx); - const auto C = this->CoefficientsCPortal.Get(idx); - const auto D = this->CoefficientsDPortal.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 true; @@ -62,6 +62,7 @@ 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; -- GitLab From 6abef8691e0d02c889d6889e1ce767ba2d499d3c Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Mon, 24 Mar 2025 11:47:00 -0400 Subject: [PATCH 05/13] Code cleanup. --- vtkm/cont/testing/UnitTestCubicSpline.cxx | 261 ++-------------------- 1 file changed, 18 insertions(+), 243 deletions(-) diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx index bb357b8099..2e9abdd320 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -17,215 +17,6 @@ #include #include -class Spline -{ -public: - // Constructor: expects x and y to have the same size and at least 2 points. - Spline(const std::vector& x, const std::vector& y) - : x_(x) - , y_(y) - , n_(x.size()) - { - if (n_ < 2) - { - throw std::invalid_argument("At least two points are required for spline interpolation."); - } - if (x_.size() != y_.size()) - { - throw std::invalid_argument("x and y must have the same length."); - } - computeCoefficients(); - } - - // Evaluate the spline at a given x value. - vtkm::FloatDefault evaluate(vtkm::FloatDefault x_val) const - { - // Ensure x_val is within the interpolation range. - if (x_val < x_.front() || x_val > x_.back()) - { - throw std::out_of_range("x_val is outside the interpolation range."); - } - // Find the interval index i such that x_[i] <= x_val <= x_[i+1]. - size_t i = std::upper_bound(x_.begin(), x_.end(), x_val) - x_.begin() - 1; - vtkm::FloatDefault dx = x_val - x_[i]; - return a_[i] + b_[i] * dx + c_[i] * dx * dx + d_[i] * dx * dx * dx; - } - -private: - std::vector x_; // x coordinates - std::vector y_; // y coordinates - size_t n_; // number of points - - // Spline coefficients for each interval i = 0, ..., n-2. - // For each segment, the cubic polynomial is: - // S_i(x) = a_[i] + b_[i]*(x-x_[i]) + c_[i]*(x-x_[i])^2 + d_[i]*(x-x_[i])^3 - std::vector a_, b_, c_, d_; - - // Compute the spline coefficients using the standard algorithm. - void computeCoefficients() - { - // Natural cubic spline: second derivatives at endpoints are zero. - // Step 1: Compute h[i] = x[i+1] - x[i] for i = 0 ... n-2. - std::vector h(n_ - 1); - for (size_t i = 0; i < n_ - 1; ++i) - { - h[i] = x_[i + 1] - x_[i]; - } - - // Step 2: Compute the array 'alpha' for i = 1 ... n-2: - // alpha[i] = 3/h[i]*(y[i+1]-y[i]) - 3/h[i-1]*(y[i]-y[i-1]) - std::vector alpha(n_, 0.0); - for (size_t i = 1; i < n_ - 1; ++i) - { - alpha[i] = (3.0 / h[i]) * (y_[i + 1] - y_[i]) - (3.0 / h[i - 1]) * (y_[i] - y_[i - 1]); - } - - // Step 3: Set up the tridiagonal system and solve it. - std::vector l(n_), mu(n_), z(n_); - l[0] = 1.0; - mu[0] = z[0] = 0.0; - for (size_t i = 1; i < n_ - 1; ++i) - { - l[i] = 2.0 * (x_[i + 1] - x_[i - 1]) - h[i - 1] * mu[i - 1]; - mu[i] = h[i] / l[i]; - z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]; - } - l[n_ - 1] = 1.0; - z[n_ - 1] = 0.0; - c_.resize(n_); - c_[n_ - 1] = 0.0; - - // Step 4: Back substitution to solve for c[i]. - for (int j = n_ - 2; j >= 0; --j) - { - c_[j] = z[j] - mu[j] * c_[j + 1]; - } - - // Step 5: Compute the coefficients b and d for each interval. - a_.resize(n_ - 1); - b_.resize(n_ - 1); - d_.resize(n_ - 1); - for (size_t i = 0; i < n_ - 1; ++i) - { - a_[i] = y_[i]; // a[i] is simply the y value at x[i] - b_[i] = (y_[i + 1] - y_[i]) / h[i] - h[i] * (c_[i + 1] + 2.0 * c_[i]) / 3.0; - d_[i] = (c_[i + 1] - c_[i]) / (3.0 * h[i]); - } - } -}; - -class Spline2 -{ -public: - Spline2(const std::vector& x, const std::vector& y) - : x1D(x) - , y1D(y) - , n(x.size()) - { - if (x1D.size() != y1D.size()) - throw std::invalid_argument("x and y must have the same length."); - if (n < 2) - throw std::invalid_argument("At least two points are required for spline interpolation."); - } - - // Evaluate the spline at a given x value. - vtkm::FloatDefault evaluate(vtkm::FloatDefault x_val) const - { - // Compute the spline coefficients on the fly. - Coefs coefs = computeCoefficients(); - - // Ensure x_val is within the interpolation range. - if (x_val < x1D.front() || x_val > x1D.back()) - throw std::out_of_range("x_val is outside the interpolation range."); - - // Find the interval index i such that x1D[i] <= x_val <= x1D[i+1] - size_t i = findInterval(x1D, x_val); - - vtkm::FloatDefault dx = x_val - x1D[i]; - // Evaluate S_i(x) = a[i] + b[i]*(x - x_i) + c[i]*(x - x_i)^2 + d[i]*(x - x_i)^3 - return coefs.a[i] + coefs.b[i] * dx + coefs.c[i] * dx * dx + coefs.d[i] * dx * dx * dx; - } - -private: - std::vector x1D, y1D; - size_t n; - - // Structure to hold the spline coefficients for each interval. - struct Coefs - { - std::vector a; - std::vector b; - std::vector c; - std::vector d; - }; - - // Compute the cubic spline coefficients using the natural cubic spline algorithm. - // This function computes the full set of coefficients for all segments. - Coefs computeCoefficients() const - { - Coefs coefs; - coefs.a = y1D; // a[i] = y1D[i] - - // Allocate coefficient arrays for each segment (there are n-1 segments) - coefs.b.resize(n - 1); - coefs.d.resize(n - 1); - coefs.c.resize(n); // second derivatives for each point - - // Step 1: Compute h[i] = x1D[i+1] - x1D[i] - std::vector h(n - 1); - for (size_t i = 0; i < n - 1; ++i) - { - h[i] = x1D[i + 1] - x1D[i]; - } - - // Step 2: Compute alpha[i] for i=1..n-2 - std::vector alpha(n, 0.0); - for (size_t i = 1; i < n - 1; ++i) - { - alpha[i] = (3.0 / h[i]) * (y1D[i + 1] - y1D[i]) - (3.0 / h[i - 1]) * (y1D[i] - y1D[i - 1]); - } - - // Step 3: Set up and solve the tridiagonal system - std::vector l(n, 0.0), mu(n, 0.0), z(n, 0.0); - l[0] = 1.0; - mu[0] = 0.0; - z[0] = 0.0; - for (size_t i = 1; i < n - 1; ++i) - { - l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; - mu[i] = h[i] / l[i]; - z[i] = (alpha[i] - h[i - 1] * z[i - 1]) / l[i]; - } - l[n - 1] = 1.0; - z[n - 1] = 0.0; - coefs.c[n - 1] = 0.0; - - // Step 4: Back substitution to solve for c[i] - for (int j = n - 2; j >= 0; --j) - { - coefs.c[j] = z[j] - mu[j] * coefs.c[j + 1]; - } - - // Step 5: Compute b[i] and d[i] for i=0..n-2 - for (size_t i = 0; i < n - 1; ++i) - { - coefs.b[i] = (y1D[i + 1] - y1D[i]) / h[i] - h[i] * (coefs.c[i + 1] + 2.0 * coefs.c[i]) / 3.0; - coefs.d[i] = (coefs.c[i + 1] - coefs.c[i]) / (3.0 * h[i]); - } - - return coefs; - } - - // Find the interval index for a given value. - size_t findInterval(const std::vector& x, vtkm::FloatDefault val) const - { - if (val < x.front() || val > x.back()) - throw std::out_of_range("Value is outside the interpolation range."); - size_t idx = std::upper_bound(x.begin(), x.end(), val) - x.begin() - 1; - return (idx >= x.size() - 1) ? x.size() - 2 : idx; - } -}; - namespace { @@ -277,21 +68,12 @@ void CheckEvaluation(const vtkm::cont::CubicSpline& spline, } } -void SaveSamples(const vtkm::cont::CubicSpline& spline_) +//Convenience function to save a bunch of samples to a file. +void SaveSamples(const vtkm::cont::CubicSpline& spline) { - vtkm::Id n = spline_.GetControlPoints().GetNumberOfValues(); - vtkm::FloatDefault t = spline_.GetControlPoints().ReadPortal().Get(0); - vtkm::FloatDefault tEnd = spline_.GetControlPoints().ReadPortal().Get(n - 1); - - std::vector x1D, y1D; - for (vtkm::Id i = 0; i < n; i++) - { - x1D.push_back(spline_.GetControlPoints().ReadPortal().Get(i)); - y1D.push_back(spline_.GetValues().ReadPortal().Get(i)); - } - - //Spline spline(x1D, y1D); - Spline2 spline(x1D, y1D); + vtkm::Id n = spline.GetControlPoints().GetNumberOfValues(); + vtkm::FloatDefault t = spline.GetControlPoints().ReadPortal().Get(0); + vtkm::FloatDefault tEnd = spline.GetControlPoints().ReadPortal().Get(n - 1); vtkm::FloatDefault dt = (tEnd - t) / (n * 100); std::vector params; @@ -300,38 +82,33 @@ void SaveSamples(const vtkm::cont::CubicSpline& spline_) params.push_back(t); t += dt; } - - /* auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); vtkm::cont::Invoker invoke; - SplineEvalWorklet worklet; vtkm::cont::ArrayHandle result; vtkm::cont::ArrayHandle valid; - invoke(worklet, paramsAH, spline, result, valid); -*/ + invoke(SplineEvalWorklet{}, paramsAH, spline, result, valid); std::cout << "SaveSamples" << std::endl; - std::ofstream fout("/Users/dpn/output.txt"), ptsOut("/Users/dpn/pts.txt"); + std::ofstream fout("output.txt"), ptsOut("pts.txt"); fout << "X,Y" << std::endl; - for (std::size_t i = 0; i < params.size(); i++) + for (vtkm::Id i = 0; i < paramsAH.GetNumberOfValues(); i++) { - fout << params[i] << "," << spline.evaluate(params[i]) << std::endl; - //fout << params[i] << "," << result.ReadPortal().Get(i) << std::endl; + fout << paramsAH.ReadPortal().Get(i) << "," << result.ReadPortal().Get(i) << std::endl; } ptsOut << "X,Y" << std::endl; - for (std::size_t i = 0; i < x1D.size(); i++) - ptsOut << x1D[i] << ", " << y1D[i] << std::endl; + auto pts = spline.GetControlPoints().ReadPortal(); + auto vals = spline.GetValues().ReadPortal(); + for (vtkm::Id i = 0; i < pts.GetNumberOfValues(); i++) + ptsOut << pts.Get(i) << ", " << vals.Get(i) << std::endl; } void CubicSplineTest() { - //std::vector xVals = { 0.0f, 1.0f, 1.5f, 2.0f, 3.0f, 4.0f }; - //std::vector yVals = { 0.0f, -1.0f, 2.0f, 1.0f, -1.0f, 0.0f }; std::vector xVals = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f }; - std::vector yVals = { 0.0f, 1.0f, 4.4f, 1.0f, 0.0f }; + std::vector yVals = { 0.0f, 1.0f, -1.0f, 1.0f, 0.0f }; auto xValsAH = vtkm::cont::make_ArrayHandle(xVals, vtkm::CopyFlag::On); auto yValsAH = vtkm::cont::make_ArrayHandle(yVals, vtkm::CopyFlag::On); @@ -341,20 +118,18 @@ void CubicSplineTest() cubicSpline.SetValues(yValsAH); cubicSpline.Update(); - SaveSamples(cubicSpline); - + //SaveSamples(cubicSpline); //Ensure that the values at control points are properly interpolated. auto tVals = xVals; auto res = yVals; CheckEvaluation(cubicSpline, tVals, res); - /* //Evaluate between control points. - tVals = { 0.5, 1.3, 2.1, 3.7 }; - res = { -1.68413, 1.00957, 3.02636, -0.461989 }; + tVals = { 0.6, 1.4, 2.16, 3.51198 }; + res = { 1.03886, 0.110853, -0.890431, 0.91292 }; + CheckEvaluation(cubicSpline, tVals, res); - */ } } // anonymous namespace -- GitLab From 5f005fa7f311d375716b438bdefd4a14ff3103af Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 1 Apr 2025 07:37:53 -0400 Subject: [PATCH 06/13] Add cubic hermite splines. --- vtkm/cont/ArrayHandleIsMonotonic.h | 36 ++++ vtkm/cont/ArrayHandleIsMonotonic.hxx | 106 ++++++++++ vtkm/cont/CMakeLists.txt | 6 +- vtkm/cont/CubicHermiteSpline.cxx | 159 +++++++++++++++ vtkm/cont/CubicHermiteSpline.h | 95 +++++++++ vtkm/cont/testing/CMakeLists.txt | 1 + .../UnitTestArrayHandleIsMonotonic.cxx | 86 ++++++++ vtkm/cont/testing/UnitTestCubicSpline.cxx | 184 +++++++++++++----- vtkm/exec/CMakeLists.txt | 1 + vtkm/exec/CubicHermiteSpline.h | 133 +++++++++++++ 10 files changed, 755 insertions(+), 52 deletions(-) create mode 100644 vtkm/cont/ArrayHandleIsMonotonic.h create mode 100644 vtkm/cont/ArrayHandleIsMonotonic.hxx create mode 100644 vtkm/cont/CubicHermiteSpline.cxx create mode 100644 vtkm/cont/CubicHermiteSpline.h create mode 100644 vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx create mode 100644 vtkm/exec/CubicHermiteSpline.h diff --git a/vtkm/cont/ArrayHandleIsMonotonic.h b/vtkm/cont/ArrayHandleIsMonotonic.h new file mode 100644 index 0000000000..ccfb64880f --- /dev/null +++ b/vtkm/cont/ArrayHandleIsMonotonic.h @@ -0,0 +1,36 @@ +//============================================================================ +// 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_ArrayHandleIsMonotonic_h +#define vtk_m_cont_ArrayHandleIsMonotonic_h + +#include +#include + +namespace vtkm +{ +namespace cont +{ + +template +VTKM_ALWAYS_EXPORT inline bool IsMonotonicIncreasing( + const vtkm::cont::ArrayHandle& input); + +template +VTKM_ALWAYS_EXPORT inline bool IsMonotonicDecreasing( + const vtkm::cont::ArrayHandle& input); + +} +} // namespace vtkm::cont + +#ifndef vtk_m_cont_ArrayHandleIsMonotonic_hxx +#include +#endif //vtk_m_cont_ArrayHandleIsMonotonic_hxx + +#endif //vtk_m_cont_ArrayHandleIsMonotonic_h diff --git a/vtkm/cont/ArrayHandleIsMonotonic.hxx b/vtkm/cont/ArrayHandleIsMonotonic.hxx new file mode 100644 index 0000000000..1f1900db80 --- /dev/null +++ b/vtkm/cont/ArrayHandleIsMonotonic.hxx @@ -0,0 +1,106 @@ +//============================================================================ +// 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_ArrayHandleIsMonotonic_hxx +#define vtk_m_cont_ArrayHandleIsMonotonic_hxx + +#include +#include +#include + +namespace vtkm +{ +namespace cont +{ + +namespace +{ +struct MonotonicIncreasing : public vtkm::worklet::WorkletMapField +{ + MonotonicIncreasing(const vtkm::Id numPoints) + : NumPoints(numPoints) + { + } + + using ControlSignature = void(WholeArrayIn, FieldOut); + using ExecutionSignature = void(InputIndex, _1, _2); + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, const ArrayType& input, bool& result) const + { + if (idx == 0) + result = true; + else + result = input.Get(idx) >= input.Get(idx - 1); + } + + vtkm::Id NumPoints; +}; + +struct MonotonicDecreasing : public vtkm::worklet::WorkletMapField +{ + MonotonicDecreasing(const vtkm::Id numPoints) + : NumPoints(numPoints) + { + } + + using ControlSignature = void(WholeArrayIn, FieldOut); + using ExecutionSignature = void(InputIndex, _1, _2); + + template + VTKM_EXEC void operator()(const vtkm::Id& idx, const ArrayType& input, bool& result) const + { + if (idx == 0) + result = true; + else + result = input.Get(idx) <= input.Get(idx - 1); + } + + vtkm::Id NumPoints; +}; +} //anonymous namespace + + +template +VTKM_ALWAYS_EXPORT inline bool IsMonotonicIncreasing( + const vtkm::cont::ArrayHandle& input) +{ + // TODO or false for empty array? + vtkm::Id numValues = input.GetNumberOfValues(); + if (numValues < 2) + return true; + + vtkm::cont::Invoker invoke; + + vtkm::cont::ArrayHandle result; + invoke(MonotonicIncreasing(input.GetNumberOfValues()), input, result); + return vtkm::cont::Algorithm::Reduce(result, true, vtkm::LogicalAnd()); +} + +template +VTKM_ALWAYS_EXPORT inline bool IsMonotonicDecreasing( + const vtkm::cont::ArrayHandle& input) +{ + // TODO or false for empty array? + vtkm::Id numValues = input.GetNumberOfValues(); + if (numValues < 2) + return true; + + vtkm::cont::Invoker invoke; + + vtkm::cont::ArrayHandle result; + invoke(MonotonicDecreasing(input.GetNumberOfValues()), input, result); + return vtkm::cont::Algorithm::Reduce(result, true, vtkm::LogicalAnd()); +} + +} +} // namespace vtkm::cont + +#endif //vtk_m_cont_ArrayHandleIsMonotonic_hxx diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index 381abad219..371a4ad87f 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -30,6 +30,7 @@ set(headers ArrayHandleGroupVecVariable.h ArrayHandleImplicit.h ArrayHandleIndex.h + ArrayHandleIsMonotonic.h ArrayHandleMultiplexer.h ArrayHandleOffsetsToNumComponents.h ArrayHandlePermutation.h @@ -80,6 +81,7 @@ set(headers ColorTableSamples.h ConvertNumComponentsToOffsets.h CoordinateSystem.h + CubicHermiteSpline.h CubicSpline.h DataSet.h DataSetBuilderCurvilinear.h @@ -131,6 +133,7 @@ set(headers ) set(template_sources + ArrayHandleIsMonotonic.hxx CellSetExplicit.hxx ParticleArrayCopy.hxx ) @@ -154,7 +157,6 @@ set(sources CellSetStructured.cxx ColorTablePresets.cxx CoordinateSystem.cxx - CubicSpline.cxx DataSet.cxx DataSetBuilderCurvilinear.cxx DataSetBuilderExplicit.cxx @@ -201,6 +203,8 @@ set(device_sources CellSetExtrude.cxx ColorTable.cxx ConvertNumComponentsToOffsets.cxx + CubicHermiteSpline.cxx + CubicSpline.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 0000000000..9e9ec78a2b --- /dev/null +++ b/vtkm/cont/CubicHermiteSpline.cxx @@ -0,0 +1,159 @@ +//============================================================================ +// 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 + +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 DivideBy : public vtkm::worklet::WorkletMapField +{ + DivideBy(const vtkm::FloatDefault& divisor) + : Divisor(divisor) + { + } + + using ControlSignature = void(FieldInOut); + using ExecutionSignature = void(_1); + + VTKM_EXEC void operator()(vtkm::FloatDefault& val) const { val = val / this->Divisor; } + + vtkm::FloatDefault Divisor; +}; + +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; + std::cout << "Tan: " << idx << " = " << tangent << " dx: " << dX << " " << dT << std::endl; + } + + vtkm::Id NumPoints; +}; +} //anonymous namespace + +VTKM_CONT vtkm::Range CubicHermiteSpline::GetParametricRange() const +{ + auto portal = this->Knots.ReadPortal(); + auto n = portal.GetNumberOfValues(); + return { portal.Get(0), portal.Get(n - 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."); + + invoker(DivideBy{ sum }, this->Knots); + std::cout << ":: params= "; + vtkm::cont::printSummary_ArrayHandle(this->Knots, std::cout); +} + +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) const +{ + vtkm::Id n = this->Data.GetNumberOfValues(); + if (n < 2) + throw std::invalid_argument("At least two points are required for spline interpolation."); + + bool isMonotonic = vtkm::cont::IsMonotonicIncreasing(this->Knots); + if (!isMonotonic) + throw std::invalid_argument("Error. Knots must be monotonic increasing."); + + if (n != this->Knots.GetNumberOfValues()) + throw std::invalid_argument("Number of data points must match the number of knots."); + if (n != this->Tangents.GetNumberOfValues()) + throw std::invalid_argument("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 0000000000..a45febfef5 --- /dev/null +++ b/vtkm/cont/CubicHermiteSpline.h @@ -0,0 +1,95 @@ +//============================================================================ +// 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; + + CubicHermiteSpline(const std::vector& data, vtkm::CopyFlag copy = vtkm::CopyFlag::On) + : Data(vtkm::cont::make_ArrayHandle(data, copy)) + { + this->ComputeKnots(); + this->ComputeTangents(); + } + + CubicHermiteSpline(const std::vector& data, + const std::vector& tangents, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + : Data(vtkm::cont::make_ArrayHandle(data, copy)) + , Tangents(vtkm::cont::make_ArrayHandle(tangents, copy)) + { + this->ComputeKnots(); + } + + CubicHermiteSpline(const std::vector& data, + const std::vector& knots, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + : Data(vtkm::cont::make_ArrayHandle(data, copy)) + , Knots(vtkm::cont::make_ArrayHandle(knots, copy)) + { + this->ComputeTangents(); + } + + CubicHermiteSpline(const std::vector& data, + const std::vector& tangents, + const std::vector& knots, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) + : Data(vtkm::cont::make_ArrayHandle(data, copy)) + , Knots(vtkm::cont::make_ArrayHandle(knots, copy)) + , Tangents(vtkm::cont::make_ArrayHandle(tangents, copy)) + { + } + + VTKM_CONT vtkm::Range GetParametricRange() const; + + VTKM_CONT const vtkm::cont::ArrayHandle& GetData() const { return this->Data; } + VTKM_CONT const vtkm::cont::ArrayHandle& GetTangents() const + { + return this->Tangents; + } + VTKM_CONT const vtkm::cont::ArrayHandle& GetKnots() const + { + return this->Knots; + } + + VTKM_CONT vtkm::exec::CubicHermiteSpline PrepareForExecution(vtkm::cont::DeviceAdapterId device, + vtkm::cont::Token& token) const; + +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 10431a2e3e..8c300b87b8 100644 --- a/vtkm/cont/testing/CMakeLists.txt +++ b/vtkm/cont/testing/CMakeLists.txt @@ -22,6 +22,7 @@ set(unit_tests UnitTestArrayHandleCounting.cxx UnitTestArrayHandleDiscard.cxx UnitTestArrayHandleIndex.cxx + UnitTestArrayHandleIsMonotonic.cxx UnitTestArrayHandleOffsetsToNumComponents.cxx UnitTestArrayHandleRandomUniformBits.cxx UnitTestArrayHandleReverse.cxx diff --git a/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx b/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx new file mode 100644 index 0000000000..abe893431c --- /dev/null +++ b/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx @@ -0,0 +1,86 @@ +//============================================================================ +// 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 + +namespace +{ +template +void CheckArray(const std::vector& input) +{ + auto array = vtkm::cont::make_ArrayHandle(input, vtkm::CopyFlag::Off); + bool isInc = vtkm::cont::IsMonotonicIncreasing(array); + bool isDec = vtkm::cont::IsMonotonicDecreasing(array); + + if (input.empty() || input.size() == 1) + { + VTKM_TEST_ASSERT(isInc && isDec, + "Array with zero or 1 should be both monotonic increasing and decreasing"); + return; + } + + VTKM_TEST_ASSERT(isInc, "Array is not monotonic increasing"); + VTKM_TEST_ASSERT(!isDec, "Array should not be monotonic decreasing"); + + //Check the reverse of the array + auto copy = input; + std::reverse(copy.begin(), copy.end()); + array = vtkm::cont::make_ArrayHandle(copy, vtkm::CopyFlag::Off); + isInc = vtkm::cont::IsMonotonicIncreasing(array); + isDec = vtkm::cont::IsMonotonicDecreasing(array); + + VTKM_TEST_ASSERT(!isInc, "Reversed array is not monotonic decreasing"); + VTKM_TEST_ASSERT(isDec, "Reversed array is not monotonic increasing"); +} + +template +std::vector ConvertVec(const std::vector& input) +{ + std::vector output; + output.reserve(input.size()); + + for (const auto& value : input) + output.push_back(static_cast(value)); + return output; +} + +void CheckTypes(const std::vector& input) +{ + CheckArray(input); + CheckArray(ConvertVec(input)); + CheckArray(ConvertVec(input)); + CheckArray(ConvertVec(input)); +} + +} //namespace anonymous + +void TestArrayHandleIsMonotonic() +{ + CheckTypes({ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }); + CheckTypes({ -5, -4, -3, -2, -1, 0, 1, 2, 3, 4 }); + + //check duplicate values + CheckTypes({ 0, 1, 1, 2, 3, 4, 4, 5, 6 }); + CheckTypes({ -3, -2, -2, -1, 0, 0, 1, 2, 3 }); + + //check empty and single element arrays + CheckTypes({}); + CheckTypes({ 0 }); +} + +int UnitTestArrayHandleIsMonotonic(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestArrayHandleIsMonotonic, argc, argv); +} diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx index 2e9abdd320..3fc6b4e401 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -11,7 +11,9 @@ #include #include +#include #include +#include #include #include #include @@ -25,55 +27,48 @@ class SplineEvalWorklet : public vtkm::worklet::WorkletMapField public: SplineEvalWorklet() {} - using ControlSignature = void(FieldIn param, - ExecObject cubicSpline, - FieldOut value, - FieldOut valid); - using ExecutionSignature = void(_1, _2, _3, _4); + 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, - bool& valid) const + ResultType& value) const { - valid = spline.Evaluate(param, value); + bool valid = spline.Evaluate(param, value); + if (!valid) + throw vtkm::cont::ErrorBadValue("Spline evaluation failed."); } }; -void CheckEvaluation(const vtkm::cont::CubicSpline& spline, +template +void CheckEvaluation(const SplineType& spline, const std::vector& params, - const std::vector& answer) + const std::vector& answer) { auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); vtkm::cont::Invoker invoke; - SplineEvalWorklet worklet; - vtkm::cont::ArrayHandle result; - vtkm::cont::ArrayHandle valid; - invoke(worklet, paramsAH, spline, result, valid); + vtkm::cont::ArrayHandle result; - VTKM_TEST_ASSERT(result.GetNumberOfValues() == static_cast(answer.size()), - "Result wrong length."); - return; + invoke(SplineEvalWorklet{}, paramsAH, spline, result); - for (std::size_t i = 0; i < answer.size(); i++) - { - VTKM_TEST_ASSERT(valid.ReadPortal().Get(static_cast(i)), "Evaluation failed."); - auto val = result.ReadPortal().Get(static_cast(i)); - std::cout << params[i] << " = " << val << " " << answer[i] - << " :: diff: " << vtkm::Abs(val - answer[i]) << std::endl; - VTKM_TEST_ASSERT(vtkm::Abs(val - answer[i]) < 1e-6, "Result has wrong value."); - } + VTKM_TEST_ASSERT( + test_equal_ArrayHandles(result, vtkm::cont::make_ArrayHandle(answer, vtkm::CopyFlag::Off))); } //Convenience function to save a bunch of samples to a file. -void SaveSamples(const vtkm::cont::CubicSpline& spline) +void SaveSamples(const std::vector& x, const std::vector& y) { - vtkm::Id n = spline.GetControlPoints().GetNumberOfValues(); - vtkm::FloatDefault t = spline.GetControlPoints().ReadPortal().Get(0); - vtkm::FloatDefault tEnd = spline.GetControlPoints().ReadPortal().Get(n - 1); + vtkm::cont::CubicSpline cubicSpline; + cubicSpline.SetControlPoints(vtkm::cont::make_ArrayHandle(x, vtkm::CopyFlag::Off)); + cubicSpline.SetValues(vtkm::cont::make_ArrayHandle(y, vtkm::CopyFlag::Off)); + cubicSpline.Update(); + + vtkm::Id n = cubicSpline.GetControlPoints().GetNumberOfValues(); + vtkm::FloatDefault t = x[0]; + vtkm::FloatDefault tEnd = x[x.size() - 1]; vtkm::FloatDefault dt = (tEnd - t) / (n * 100); std::vector params; @@ -82,12 +77,12 @@ void SaveSamples(const vtkm::cont::CubicSpline& spline) params.push_back(t); t += dt; } + auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); vtkm::cont::Invoker invoke; vtkm::cont::ArrayHandle result; - vtkm::cont::ArrayHandle valid; - invoke(SplineEvalWorklet{}, paramsAH, spline, result, valid); + invoke(SplineEvalWorklet{}, paramsAH, cubicSpline, result); std::cout << "SaveSamples" << std::endl; std::ofstream fout("output.txt"), ptsOut("pts.txt"); @@ -98,38 +93,125 @@ void SaveSamples(const vtkm::cont::CubicSpline& spline) } ptsOut << "X,Y" << std::endl; - auto pts = spline.GetControlPoints().ReadPortal(); - auto vals = spline.GetValues().ReadPortal(); + auto pts = cubicSpline.GetControlPoints().ReadPortal(); + auto vals = cubicSpline.GetValues().ReadPortal(); for (vtkm::Id i = 0; i < pts.GetNumberOfValues(); i++) ptsOut << pts.Get(i) << ", " << vals.Get(i) << std::endl; } +void DoTest(const std::vector& x, + const std::vector& y, + const std::vector& xSamples, + const std::vector& vals) +{ + /* + vtkm::cont::CubicSpline cubicSpline; + cubicSpline.SetControlPoints(vtkm::cont::make_ArrayHandle(x, vtkm::CopyFlag::Off)); + cubicSpline.SetValues(vtkm::cont::make_ArrayHandle(y, vtkm::CopyFlag::Off)); + cubicSpline.Update(); + + //Make sure we the spline interpolates the control points. + CheckEvaluation(cubicSpline, x, y); -void CubicSplineTest() + //Make sure the spline evaluates to the correct values at the given points. + CheckEvaluation(cubicSpline, xSamples, vals); + */ +} + +void CubicHermiteSplineTest() { - std::vector xVals = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f }; - std::vector yVals = { 0.0f, 1.0f, -1.0f, 1.0f, 0.0f }; + std::vector pts = { { 0, 0, 0 }, { 1, 1, 0 }, { 2, 1, 0 }, { 3, -.5, 0 }, + { 4, -1, 0 }, { 5, -1, 0 }, { 6, 0, 0 } }; - auto xValsAH = vtkm::cont::make_ArrayHandle(xVals, vtkm::CopyFlag::On); - auto yValsAH = vtkm::cont::make_ArrayHandle(yVals, vtkm::CopyFlag::On); + std::vector knots; + pts.clear(); - vtkm::cont::CubicSpline cubicSpline; - cubicSpline.SetControlPoints(xValsAH); - cubicSpline.SetValues(yValsAH); - cubicSpline.Update(); + vtkm::Id n = 100; + vtkm::FloatDefault t = 0.0, dt = vtkm::TwoPi() / static_cast(n); + + while (t <= vtkm::TwoPi()) + { + vtkm::FloatDefault x = vtkm::Cos(t); + vtkm::FloatDefault y = vtkm::Sin(t); + vtkm::FloatDefault z = vtkm::Cos(t) * vtkm::Sin(t); + pts.push_back({ x, y, z }); + knots.push_back(t); + t += dt; + } + vtkm::cont::CubicHermiteSpline spline0(pts, knots); + auto parametricRange = spline0.GetParametricRange(); + + t = parametricRange.Min; + dt = parametricRange.Length() / 1500; + std::vector vals; + while (t <= parametricRange.Max) + { + vals.push_back(t); + t += dt; + } + std::cout << "KNOTS= "; + vtkm::cont::printSummary_ArrayHandle(spline0.GetKnots(), std::cout); + + vtkm::cont::Invoker invoke; + vtkm::cont::ArrayHandle result; + invoke( + SplineEvalWorklet{}, vtkm::cont::make_ArrayHandle(vals, vtkm::CopyFlag::Off), spline0, result); + + /* + vtkm::cont::CubicHermiteSpline spline(pts); - //SaveSamples(cubicSpline); + vtkm::Id n = spline.GetData().GetNumberOfValues(); + auto t0 = spline.GetKnots().ReadPortal().Get(0); + auto t1 = spline.GetKnots().ReadPortal().Get(n-1); - //Ensure that the values at control points are properly interpolated. - auto tVals = xVals; - auto res = yVals; - CheckEvaluation(cubicSpline, tVals, res); + vtkm::FloatDefault dt = (t1 - t0) / (n * 100); + std::vector params; + vtkm::FloatDefault t = t0; + + while (t < t1) + { + params.push_back(t); + t += dt; + } + + vtkm::cont::Invoker invoke; + vtkm::cont::ArrayHandle result; + invoke(SplineEvalWorklet{}, vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off), spline, result); + */ + + std::cout << "SaveSamples" << std::endl; + std::ofstream fout("/Users/dpn/output.txt"), ptsOut("/Users/dpn/pts.txt"); + fout << "T, X, Y, Z" << std::endl; + ptsOut << "T, X, Y, Z" << std::endl; + auto knots_ = spline0.GetKnots().ReadPortal(); + for (vtkm::Id i = 0; i < pts.size(); i++) + ptsOut << knots_.Get(i) << ", " << pts[i][0] << ", " << pts[i][1] << ", " << pts[i][2] + << std::endl; + + auto res = result.ReadPortal(); + for (vtkm::Id i = 0; i < res.GetNumberOfValues(); i++) + fout << vals[i] << ", " << res.Get(i)[0] << ", " << res.Get(i)[1] << ", " << res.Get(i)[2] + << std::endl; +} + +void CubicSplineTest() +{ + CubicHermiteSplineTest(); + + //uniform spacing wavy function + std::vector xVals = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f }; + std::vector yVals = { 0.0f, 1.0f, -1.0f, 1.0f, 0.0f }; + std::vector samples = { 0.6, 1.4, 2.16, 3.51198 }; + std::vector res = { 1.03886, 0.110853, -0.890431, 0.91292 }; - //Evaluate between control points. - tVals = { 0.6, 1.4, 2.16, 3.51198 }; - res = { 1.03886, 0.110853, -0.890431, 0.91292 }; + //DoTest(xVals, yVals, samples, res); - CheckEvaluation(cubicSpline, tVals, res); + //non-uniform spacing wavy function + xVals = { 0.0f, 1.0f, 1.1f, 3.9f, 4.0f, 5.0f }; + yVals = { 0.0f, 1.0f, 1.2f, 1.0f, 0.0f, 0.0f }; + samples = { 0.7f, 1.45f, 3.65f, 4.38f }; + res = { 0.548318, 2.15815f, 3.13104f, -1.77143 }; + //DoTest(xVals, yVals, samples, res); } } // anonymous namespace diff --git a/vtkm/exec/CMakeLists.txt b/vtkm/exec/CMakeLists.txt index a3283634d5..621e6cc863 100644 --- a/vtkm/exec/CMakeLists.txt +++ b/vtkm/exec/CMakeLists.txt @@ -29,6 +29,7 @@ set(headers ConnectivityExtrude.h ConnectivityPermuted.h ConnectivityStructured.h + CubicHermiteSpline.h CubicSpline.h FieldNeighborhood.h FunctorBase.h diff --git a/vtkm/exec/CubicHermiteSpline.h b/vtkm/exec/CubicHermiteSpline.h new file mode 100644 index 0000000000..e2bfa4cac9 --- /dev/null +++ b/vtkm/exec/CubicHermiteSpline.h @@ -0,0 +1,133 @@ +//============================================================================ +// 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 + +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 + bool Evaluate(const vtkm::FloatDefault& tVal, vtkm::Vec3f& val) const + { + vtkm::Id idx = this->FindInterval(tVal); + if (idx < 0) + { + idx = this->FindInterval(tVal); + return false; + } + + 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 true; + } + +private: + vtkm::Vec3f ComputeTangent(vtkm::Id idx) const + { + const vtkm::Id n = this->Data.GetNumberOfValues(); + vtkm::Id idx0, idx1; + if (idx == 0) // Forward difference + { + idx0 = 0; + idx1 = 1; + } + else if (idx == n - 1) // Backward difference + { + idx0 = n - 2; + idx1 = n - 1; + } + else // central difference + { + idx0 = idx - 1; + idx1 = idx + 1; + } + + auto dX = this->Data.Get(idx1) - this->Data.Get(idx0); + auto dT = this->Knots.Get(idx1) - this->Knots.Get(idx0); + std::cout << "dX= " << dX << std::endl; + std::cout << "dT= " << dT << std::endl; + + vtkm::Vec3f tanVec = dX / dT; + std::cout << "tanVec= " << tanVec << std::endl; + return tanVec; + } + + vtkm::Id FindInterval(const vtkm::FloatDefault& t) const + { + vtkm::Id n = this->Data.GetNumberOfValues(); + auto k0 = this->Knots.Get(0); + auto k1 = this->Knots.Get(n - 1); + + if (t < this->Knots.Get(0) || t > this->Knots.Get(n - 1)) + return -1; + if (t == this->Knots.Get(n - 1)) + return n - 2; + + for (vtkm::Id i = 0; i < n - 1; ++i) + { + if (t >= this->Knots.Get(i) && t < this->Knots.Get(i + 1)) + return i; + } + return -1; + } + + DataPortalType Data; + KnotsPortalType Knots; + TangentPortalType Tangents; +}; +} //namespace exec +} //namespace vtkm + +#endif //vtk_m_exec_CubicHermiteSpline_h -- GitLab From 9abe6b91fd634cf6e5ba979722da1a6ea81a7f06 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 1 Apr 2025 07:50:23 -0400 Subject: [PATCH 07/13] Evaluate returns vtkm::ErrorCode --- vtkm/ErrorCode.h | 3 +++ vtkm/cont/testing/UnitTestCubicSpline.cxx | 4 ++-- vtkm/exec/CubicHermiteSpline.h | 9 ++++----- vtkm/exec/CubicSpline.h | 7 ++++--- 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/vtkm/ErrorCode.h b/vtkm/ErrorCode.h index 7ef94a0fcf..ba35d7240f 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/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicSpline.cxx index 3fc6b4e401..776e9e259b 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicSpline.cxx @@ -36,8 +36,8 @@ public: const CubicSplineType& spline, ResultType& value) const { - bool valid = spline.Evaluate(param, value); - if (!valid) + auto res = spline.Evaluate(param, value); + if (res != vtkm::ErrorCode::Success) throw vtkm::cont::ErrorBadValue("Spline evaluation failed."); } }; diff --git a/vtkm/exec/CubicHermiteSpline.h b/vtkm/exec/CubicHermiteSpline.h index e2bfa4cac9..c9e00ef04a 100644 --- a/vtkm/exec/CubicHermiteSpline.h +++ b/vtkm/exec/CubicHermiteSpline.h @@ -11,6 +11,7 @@ #ifndef vtk_m_exec_CubicHermiteSpline_h #define vtk_m_exec_CubicHermiteSpline_h +#include #include #include @@ -42,13 +43,13 @@ public: } VTKM_EXEC - bool Evaluate(const vtkm::FloatDefault& tVal, vtkm::Vec3f& val) const + 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 false; + return vtkm::ErrorCode::ValueOutOfRange; } auto m0 = this->Tangents.Get(idx); @@ -70,7 +71,7 @@ public: 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 true; + return vtkm::ErrorCode::Success; } private: @@ -107,8 +108,6 @@ private: vtkm::Id FindInterval(const vtkm::FloatDefault& t) const { vtkm::Id n = this->Data.GetNumberOfValues(); - auto k0 = this->Knots.Get(0); - auto k1 = this->Knots.Get(n - 1); if (t < this->Knots.Get(0) || t > this->Knots.Get(n - 1)) return -1; diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h index 29ea7d8d1d..f3e84d6ab1 100644 --- a/vtkm/exec/CubicSpline.h +++ b/vtkm/exec/CubicSpline.h @@ -10,6 +10,7 @@ #ifndef vtkm_exec_cubicspline_h #define vtkm_exec_cubicspline_h +#include #include #include @@ -41,13 +42,13 @@ public: } VTKM_EXEC - bool Evaluate(const vtkm::FloatDefault& param, vtkm::FloatDefault& val) const + vtkm::ErrorCode Evaluate(const vtkm::FloatDefault& param, vtkm::FloatDefault& val) const { val = 0; vtkm::Id idx = this->FindInterval(param); if (idx < 0) - return false; + return vtkm::ErrorCode::ValueOutOfRange; vtkm::FloatDefault dx = param - this->ControlPointsPortal.Get(idx); auto B = this->CoefficientsBPortal.Get(idx); @@ -55,7 +56,7 @@ public: auto D = this->CoefficientsDPortal.Get(idx); val = this->ValuesPortal.Get(idx) + B * dx + C * dx * dx + D * dx * dx * dx; - return true; + return vtkm::ErrorCode::Success; } private: -- GitLab From 69ec3549edcfb2accf091ae571661c7968d52678 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 1 Apr 2025 08:26:33 -0400 Subject: [PATCH 08/13] remo --- vtkm/cont/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index 371a4ad87f..accb96f275 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -203,7 +203,7 @@ set(device_sources CellSetExtrude.cxx ColorTable.cxx ConvertNumComponentsToOffsets.cxx - CubicHermiteSpline.cxx + CubicHermiteSpline.cxx CubicSpline.cxx Field.cxx internal/ArrayCopyUnknown.cxx -- GitLab From 57d8a53136d183d4bad76bbf11b7b5630aba5746 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Tue, 1 Apr 2025 09:08:58 -0400 Subject: [PATCH 09/13] Use binary search in FindInterval --- vtkm/exec/CubicHermiteSpline.h | 52 +++++++++++----------------------- vtkm/exec/CubicSpline.h | 18 ++++++++++-- 2 files changed, 31 insertions(+), 39 deletions(-) diff --git a/vtkm/exec/CubicHermiteSpline.h b/vtkm/exec/CubicHermiteSpline.h index c9e00ef04a..0b49686ebe 100644 --- a/vtkm/exec/CubicHermiteSpline.h +++ b/vtkm/exec/CubicHermiteSpline.h @@ -75,50 +75,30 @@ public: } private: - vtkm::Vec3f ComputeTangent(vtkm::Id idx) const - { - const vtkm::Id n = this->Data.GetNumberOfValues(); - vtkm::Id idx0, idx1; - if (idx == 0) // Forward difference - { - idx0 = 0; - idx1 = 1; - } - else if (idx == n - 1) // Backward difference - { - idx0 = n - 2; - idx1 = n - 1; - } - else // central difference - { - idx0 = idx - 1; - idx1 = idx + 1; - } - - auto dX = this->Data.Get(idx1) - this->Data.Get(idx0); - auto dT = this->Knots.Get(idx1) - this->Knots.Get(idx0); - std::cout << "dX= " << dX << std::endl; - std::cout << "dT= " << dT << std::endl; - - vtkm::Vec3f tanVec = dX / dT; - std::cout << "tanVec= " << tanVec << std::endl; - return tanVec; - } - vtkm::Id FindInterval(const vtkm::FloatDefault& t) const { - vtkm::Id n = this->Data.GetNumberOfValues(); + vtkm::Id n = this->Knots.GetNumberOfValues(); if (t < this->Knots.Get(0) || t > this->Knots.Get(n - 1)) return -1; - if (t == this->Knots.Get(n - 1)) - return n - 2; - for (vtkm::Id i = 0; i < n - 1; ++i) + //Binary search for the interval + vtkm::Id left = 0; + vtkm::Id right = n - 1; + + while (left < right) { - if (t >= this->Knots.Get(i) && t < this->Knots.Get(i + 1)) - return i; + 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; } diff --git a/vtkm/exec/CubicSpline.h b/vtkm/exec/CubicSpline.h index f3e84d6ab1..6a5c861a04 100644 --- a/vtkm/exec/CubicSpline.h +++ b/vtkm/exec/CubicSpline.h @@ -67,11 +67,23 @@ private: if (x < this->ControlPointsPortal.Get(0) || x > this->ControlPointsPortal.Get(numPoints - 1)) return -1; - for (vtkm::Id i = 0; i < numPoints - 1; ++i) + //Binary search for the interval + vtkm::Id left = 0; + vtkm::Id right = numPoints - 1; + + while (left < right) { - if (x >= this->ControlPointsPortal.Get(i) && x <= this->ControlPointsPortal.Get(i + 1)) - return i; + 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; } -- GitLab From 2bf68fad099a3774e15e7b34b6eff1747b094417 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Fri, 4 Apr 2025 12:43:35 -0400 Subject: [PATCH 10/13] Remove the natural cubic spline and montonic algo. --- vtkm/cont/ArrayHandleIsMonotonic.h | 36 ---- vtkm/cont/ArrayHandleIsMonotonic.hxx | 106 ------------ vtkm/cont/CMakeLists.txt | 4 - vtkm/cont/CubicHermiteSpline.cxx | 5 - vtkm/cont/CubicSpline.cxx | 159 ------------------ vtkm/cont/CubicSpline.h | 78 --------- vtkm/cont/testing/CMakeLists.txt | 3 +- .../UnitTestArrayHandleIsMonotonic.cxx | 86 ---------- ...ine.cxx => UnitTestCubicHermiteSpline.cxx} | 108 +----------- 9 files changed, 3 insertions(+), 582 deletions(-) delete mode 100644 vtkm/cont/ArrayHandleIsMonotonic.h delete mode 100644 vtkm/cont/ArrayHandleIsMonotonic.hxx delete mode 100644 vtkm/cont/CubicSpline.cxx delete mode 100644 vtkm/cont/CubicSpline.h delete mode 100644 vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx rename vtkm/cont/testing/{UnitTestCubicSpline.cxx => UnitTestCubicHermiteSpline.cxx} (50%) diff --git a/vtkm/cont/ArrayHandleIsMonotonic.h b/vtkm/cont/ArrayHandleIsMonotonic.h deleted file mode 100644 index ccfb64880f..0000000000 --- a/vtkm/cont/ArrayHandleIsMonotonic.h +++ /dev/null @@ -1,36 +0,0 @@ -//============================================================================ -// 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_ArrayHandleIsMonotonic_h -#define vtk_m_cont_ArrayHandleIsMonotonic_h - -#include -#include - -namespace vtkm -{ -namespace cont -{ - -template -VTKM_ALWAYS_EXPORT inline bool IsMonotonicIncreasing( - const vtkm::cont::ArrayHandle& input); - -template -VTKM_ALWAYS_EXPORT inline bool IsMonotonicDecreasing( - const vtkm::cont::ArrayHandle& input); - -} -} // namespace vtkm::cont - -#ifndef vtk_m_cont_ArrayHandleIsMonotonic_hxx -#include -#endif //vtk_m_cont_ArrayHandleIsMonotonic_hxx - -#endif //vtk_m_cont_ArrayHandleIsMonotonic_h diff --git a/vtkm/cont/ArrayHandleIsMonotonic.hxx b/vtkm/cont/ArrayHandleIsMonotonic.hxx deleted file mode 100644 index 1f1900db80..0000000000 --- a/vtkm/cont/ArrayHandleIsMonotonic.hxx +++ /dev/null @@ -1,106 +0,0 @@ -//============================================================================ -// 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_ArrayHandleIsMonotonic_hxx -#define vtk_m_cont_ArrayHandleIsMonotonic_hxx - -#include -#include -#include - -namespace vtkm -{ -namespace cont -{ - -namespace -{ -struct MonotonicIncreasing : public vtkm::worklet::WorkletMapField -{ - MonotonicIncreasing(const vtkm::Id numPoints) - : NumPoints(numPoints) - { - } - - using ControlSignature = void(WholeArrayIn, FieldOut); - using ExecutionSignature = void(InputIndex, _1, _2); - - template - VTKM_EXEC void operator()(const vtkm::Id& idx, const ArrayType& input, bool& result) const - { - if (idx == 0) - result = true; - else - result = input.Get(idx) >= input.Get(idx - 1); - } - - vtkm::Id NumPoints; -}; - -struct MonotonicDecreasing : public vtkm::worklet::WorkletMapField -{ - MonotonicDecreasing(const vtkm::Id numPoints) - : NumPoints(numPoints) - { - } - - using ControlSignature = void(WholeArrayIn, FieldOut); - using ExecutionSignature = void(InputIndex, _1, _2); - - template - VTKM_EXEC void operator()(const vtkm::Id& idx, const ArrayType& input, bool& result) const - { - if (idx == 0) - result = true; - else - result = input.Get(idx) <= input.Get(idx - 1); - } - - vtkm::Id NumPoints; -}; -} //anonymous namespace - - -template -VTKM_ALWAYS_EXPORT inline bool IsMonotonicIncreasing( - const vtkm::cont::ArrayHandle& input) -{ - // TODO or false for empty array? - vtkm::Id numValues = input.GetNumberOfValues(); - if (numValues < 2) - return true; - - vtkm::cont::Invoker invoke; - - vtkm::cont::ArrayHandle result; - invoke(MonotonicIncreasing(input.GetNumberOfValues()), input, result); - return vtkm::cont::Algorithm::Reduce(result, true, vtkm::LogicalAnd()); -} - -template -VTKM_ALWAYS_EXPORT inline bool IsMonotonicDecreasing( - const vtkm::cont::ArrayHandle& input) -{ - // TODO or false for empty array? - vtkm::Id numValues = input.GetNumberOfValues(); - if (numValues < 2) - return true; - - vtkm::cont::Invoker invoke; - - vtkm::cont::ArrayHandle result; - invoke(MonotonicDecreasing(input.GetNumberOfValues()), input, result); - return vtkm::cont::Algorithm::Reduce(result, true, vtkm::LogicalAnd()); -} - -} -} // namespace vtkm::cont - -#endif //vtk_m_cont_ArrayHandleIsMonotonic_hxx diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index accb96f275..c686c386da 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -30,7 +30,6 @@ set(headers ArrayHandleGroupVecVariable.h ArrayHandleImplicit.h ArrayHandleIndex.h - ArrayHandleIsMonotonic.h ArrayHandleMultiplexer.h ArrayHandleOffsetsToNumComponents.h ArrayHandlePermutation.h @@ -82,7 +81,6 @@ set(headers ConvertNumComponentsToOffsets.h CoordinateSystem.h CubicHermiteSpline.h - CubicSpline.h DataSet.h DataSetBuilderCurvilinear.h DataSetBuilderExplicit.h @@ -133,7 +131,6 @@ set(headers ) set(template_sources - ArrayHandleIsMonotonic.hxx CellSetExplicit.hxx ParticleArrayCopy.hxx ) @@ -204,7 +201,6 @@ set(device_sources ColorTable.cxx ConvertNumComponentsToOffsets.cxx CubicHermiteSpline.cxx - CubicSpline.cxx Field.cxx internal/ArrayCopyUnknown.cxx internal/ArrayRangeComputeUtils.cxx diff --git a/vtkm/cont/CubicHermiteSpline.cxx b/vtkm/cont/CubicHermiteSpline.cxx index 9e9ec78a2b..b3ad06be28 100644 --- a/vtkm/cont/CubicHermiteSpline.cxx +++ b/vtkm/cont/CubicHermiteSpline.cxx @@ -10,7 +10,6 @@ #include #include -#include #include #include #include @@ -141,10 +140,6 @@ vtkm::exec::CubicHermiteSpline CubicHermiteSpline::PrepareForExecution( if (n < 2) throw std::invalid_argument("At least two points are required for spline interpolation."); - bool isMonotonic = vtkm::cont::IsMonotonicIncreasing(this->Knots); - if (!isMonotonic) - throw std::invalid_argument("Error. Knots must be monotonic increasing."); - if (n != this->Knots.GetNumberOfValues()) throw std::invalid_argument("Number of data points must match the number of knots."); if (n != this->Tangents.GetNumberOfValues()) diff --git a/vtkm/cont/CubicSpline.cxx b/vtkm/cont/CubicSpline.cxx deleted file mode 100644 index 77226d4f89..0000000000 --- a/vtkm/cont/CubicSpline.cxx +++ /dev/null @@ -1,159 +0,0 @@ -//============================================================================ -// 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 - -namespace vtkm -{ -namespace cont -{ - -void CubicSpline::Update() const -{ - if (this->Modified) - { - // Although the data of the derived class may change, the logical state - // of the class should not. Thus, we will instruct the compiler to relax - // the const constraint. - const_cast(this)->Build(); - this->Modified = false; - } -} - -vtkm::cont::ArrayHandle CubicSpline::CalcH() const -{ - vtkm::Id n = this->ControlPoints.GetNumberOfValues(); - vtkm::cont::ArrayHandle h; - h.Allocate(n - 1); - - auto hPortal = h.WritePortal(); - auto cPortal = this->ControlPoints.ReadPortal(); - for (vtkm::Id i = 0; i < n - 1; i++) - hPortal.Set(i, cPortal.Get(i + 1) - cPortal.Get(i)); - - return h; -} - -vtkm::cont::ArrayHandle CubicSpline::CalcAlpha( - const vtkm::cont::ArrayHandle& h) const -{ - vtkm::cont::ArrayHandle alpha; - vtkm::Id n = this->ControlPoints.GetNumberOfValues(); - alpha.Allocate(n - 2); - - auto Alpha = alpha.WritePortal(); - auto H = h.ReadPortal(); - auto Y = this->Values.ReadPortal(); - - for (vtkm::Id i = 1; i < n - 1; ++i) - { - //alpha[i - 1] = (3.0 * (a1D[i + 1] - a1D[i]) / h[i]) - (3.0 * (a1D[i] - a1D[i - 1]) / h[i - 1]); - Alpha.Set(i - 1, - (3.0 * (Y.Get(i + 1) - Y.Get(i)) / H.Get(i)) - - (3.0 * (Y.Get(i) - Y.Get(i - 1)) / H.Get(i - 1))); - } - - return alpha; -} - -void CubicSpline::CalcCoefficients(const vtkm::cont::ArrayHandle& H, - const vtkm::cont::ArrayHandle& Alpha) -{ - vtkm::Id n = this->ControlPoints.GetNumberOfValues(); - vtkm::cont::ArrayHandle L, Mu, Z; - - L.Allocate(n); - Mu.Allocate(n); - Z.Allocate(n); - - auto l = L.WritePortal(); - auto mu = Mu.WritePortal(); - auto z = Z.WritePortal(); - auto h = H.ReadPortal(); - auto alpha = Alpha.ReadPortal(); - - l.Set(0, 1.0); - mu.Set(0, 0.0); - z.Set(0, 0.0); - - auto xVals = this->ControlPoints.ReadPortal(); - - for (vtkm::Id i = 1; i < n - 1; ++i) - { - //l[i] = 2.0 * (x1D[i + 1] - x1D[i - 1]) - h[i - 1] * mu[i - 1]; - l.Set(i, 2.0 * (xVals.Get(i + 1) - xVals.Get(i - 1)) - h.Get(i - 1) * mu.Get(i - 1)); - //mu[i] = h[i] / l[i]; - mu.Set(i, h.Get(i) / l.Get(i)); - //z[i] = (alpha[i - 1] - h[i - 1] * z[i - 1]) / l[i]; - z.Set(i, (alpha.Get(i - 1) - h.Get(i - 1) * z.Get(i - 1)) / l.Get(i)); - } - - //l[n - 1] = 1.0; - //z[n - 1] = c1D[n - 1] = 0.0; - l.Set(n - 1, 1.0); - z.Set(n - 1, 0.0); - - this->CoefficientsB.Allocate(n - 1); - this->CoefficientsC.Allocate(n); - this->CoefficientsD.Allocate(n - 1); - - auto B = this->CoefficientsB.WritePortal(); - auto C = this->CoefficientsC.WritePortal(); - auto D = this->CoefficientsD.WritePortal(); - - C.Set(n - 1, 0.0); - auto yVals = this->Values.ReadPortal(); - - for (vtkm::Id i = n - 2; i >= 0; i--) - { - //c1D[j] = z[j] - mu[j] * c1D[j + 1]; - C.Set(i, z.Get(i) - mu.Get(i) * C.Get(i + 1)); - //b1D[j] = (a1D[j + 1] - a1D[j]) / h[j] - h[j] * (c1D[j + 1] + 2.0 * c1D[j]) / 3.0; - auto val0 = (yVals.Get(i + 1) - yVals.Get(i)) / h.Get(i); - auto val1 = h.Get(i) * (C.Get(i + 1) + 2.0 * C.Get(i)) / 3.0; - B.Set(i, val0 - val1); - //d1D[j] = (c1D[j + 1] - c1D[j]) / (3.0 * h[j]); - D.Set(i, (C.Get(i + 1) - C.Get(i)) / (3.0 * h.Get(i))); - } -} - - -void CubicSpline::Build() -{ - //Calculate the spline coeficients. - vtkm::Id n = this->ControlPoints.GetNumberOfValues(); - if (n < 2) - throw std::invalid_argument("At least two points are required for spline interpolation."); - - auto h = this->CalcH(); - auto alpha = this->CalcAlpha(h); - - this->CalcCoefficients(h, alpha); -} - -vtkm::exec::CubicSpline CubicSpline::PrepareForExecution(vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const -{ - this->Update(); - - using ExecObjType = vtkm::exec::CubicSpline; - - return ExecObjType(this->GetControlPoints(), - this->GetValues(), - this->CoefficientsB, - this->CoefficientsC, - this->CoefficientsD, - device, - token); -} - -} -} // namespace vtkm::cont diff --git a/vtkm/cont/CubicSpline.h b/vtkm/cont/CubicSpline.h deleted file mode 100644 index c3976e40df..0000000000 --- a/vtkm/cont/CubicSpline.h +++ /dev/null @@ -1,78 +0,0 @@ -//============================================================================ -// 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_CubicSpline_h -#define vtk_m_cont_CubicSpline_h - -#include - -#include -#include -#include -#include - -namespace vtkm -{ -namespace cont -{ - -class VTKM_CONT_EXPORT CubicSpline : public vtkm::cont::ExecutionObjectBase -{ -public: - CubicSpline() = default; - virtual ~CubicSpline() = default; - - VTKM_CONT void SetControlPoints(const vtkm::cont::ArrayHandle& controlPoints) - { - this->ControlPoints = controlPoints; - this->SetModified(); - } - VTKM_CONT void SetValues(const vtkm::cont::ArrayHandle& values) - { - this->Values = values; - this->SetModified(); - } - - VTKM_CONT vtkm::cont::ArrayHandle GetControlPoints() const - { - return this->ControlPoints; - } - - VTKM_CONT vtkm::cont::ArrayHandle GetValues() const { return this->Values; } - - VTKM_CONT void Update() const; - - VTKM_CONT vtkm::exec::CubicSpline PrepareForExecution(vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const; - -private: - void SetModified() { this->Modified = true; } - bool GetModified() const { return this->Modified; } - - void Build(); - - vtkm::cont::ArrayHandle CalcH() const; - vtkm::cont::ArrayHandle CalcAlpha( - const vtkm::cont::ArrayHandle& h) const; - void CalcCoefficients(const vtkm::cont::ArrayHandle& H, - const vtkm::cont::ArrayHandle& Alpha); - - vtkm::cont::ArrayHandle ControlPoints; - vtkm::cont::ArrayHandle Values; - vtkm::cont::ArrayHandle CoefficientsA; - vtkm::cont::ArrayHandle CoefficientsB; - vtkm::cont::ArrayHandle CoefficientsC; - vtkm::cont::ArrayHandle CoefficientsD; - mutable bool Modified = true; -}; - -} -} // vtkm::cont - -#endif //vtk_m_cont_CubicSpline_h diff --git a/vtkm/cont/testing/CMakeLists.txt b/vtkm/cont/testing/CMakeLists.txt index 8c300b87b8..5041d4d405 100644 --- a/vtkm/cont/testing/CMakeLists.txt +++ b/vtkm/cont/testing/CMakeLists.txt @@ -22,7 +22,6 @@ set(unit_tests UnitTestArrayHandleCounting.cxx UnitTestArrayHandleDiscard.cxx UnitTestArrayHandleIndex.cxx - UnitTestArrayHandleIsMonotonic.cxx UnitTestArrayHandleOffsetsToNumComponents.cxx UnitTestArrayHandleRandomUniformBits.cxx UnitTestArrayHandleReverse.cxx @@ -106,7 +105,7 @@ set(unit_tests_device UnitTestCellSetExplicit.cxx UnitTestCellSetPermutation.cxx UnitTestColorTable.cxx - UnitTestCubicSpline.cxx + UnitTestCubicHermiteSpline.cxx UnitTestDataSetPermutation.cxx UnitTestDataSetSingleType.cxx UnitTestDeviceAdapterAlgorithmDependency.cxx diff --git a/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx b/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx deleted file mode 100644 index abe893431c..0000000000 --- a/vtkm/cont/testing/UnitTestArrayHandleIsMonotonic.cxx +++ /dev/null @@ -1,86 +0,0 @@ -//============================================================================ -// 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 - -namespace -{ -template -void CheckArray(const std::vector& input) -{ - auto array = vtkm::cont::make_ArrayHandle(input, vtkm::CopyFlag::Off); - bool isInc = vtkm::cont::IsMonotonicIncreasing(array); - bool isDec = vtkm::cont::IsMonotonicDecreasing(array); - - if (input.empty() || input.size() == 1) - { - VTKM_TEST_ASSERT(isInc && isDec, - "Array with zero or 1 should be both monotonic increasing and decreasing"); - return; - } - - VTKM_TEST_ASSERT(isInc, "Array is not monotonic increasing"); - VTKM_TEST_ASSERT(!isDec, "Array should not be monotonic decreasing"); - - //Check the reverse of the array - auto copy = input; - std::reverse(copy.begin(), copy.end()); - array = vtkm::cont::make_ArrayHandle(copy, vtkm::CopyFlag::Off); - isInc = vtkm::cont::IsMonotonicIncreasing(array); - isDec = vtkm::cont::IsMonotonicDecreasing(array); - - VTKM_TEST_ASSERT(!isInc, "Reversed array is not monotonic decreasing"); - VTKM_TEST_ASSERT(isDec, "Reversed array is not monotonic increasing"); -} - -template -std::vector ConvertVec(const std::vector& input) -{ - std::vector output; - output.reserve(input.size()); - - for (const auto& value : input) - output.push_back(static_cast(value)); - return output; -} - -void CheckTypes(const std::vector& input) -{ - CheckArray(input); - CheckArray(ConvertVec(input)); - CheckArray(ConvertVec(input)); - CheckArray(ConvertVec(input)); -} - -} //namespace anonymous - -void TestArrayHandleIsMonotonic() -{ - CheckTypes({ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }); - CheckTypes({ -5, -4, -3, -2, -1, 0, 1, 2, 3, 4 }); - - //check duplicate values - CheckTypes({ 0, 1, 1, 2, 3, 4, 4, 5, 6 }); - CheckTypes({ -3, -2, -2, -1, 0, 0, 1, 2, 3 }); - - //check empty and single element arrays - CheckTypes({}); - CheckTypes({ 0 }); -} - -int UnitTestArrayHandleIsMonotonic(int argc, char* argv[]) -{ - return vtkm::cont::testing::Testing::Run(TestArrayHandleIsMonotonic, argc, argv); -} diff --git a/vtkm/cont/testing/UnitTestCubicSpline.cxx b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx similarity index 50% rename from vtkm/cont/testing/UnitTestCubicSpline.cxx rename to vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx index 776e9e259b..32ec21723a 100644 --- a/vtkm/cont/testing/UnitTestCubicSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx @@ -12,11 +12,9 @@ #include #include -#include #include #include #include -#include #include namespace @@ -58,66 +56,6 @@ void CheckEvaluation(const SplineType& spline, test_equal_ArrayHandles(result, vtkm::cont::make_ArrayHandle(answer, vtkm::CopyFlag::Off))); } -//Convenience function to save a bunch of samples to a file. -void SaveSamples(const std::vector& x, const std::vector& y) -{ - vtkm::cont::CubicSpline cubicSpline; - cubicSpline.SetControlPoints(vtkm::cont::make_ArrayHandle(x, vtkm::CopyFlag::Off)); - cubicSpline.SetValues(vtkm::cont::make_ArrayHandle(y, vtkm::CopyFlag::Off)); - cubicSpline.Update(); - - vtkm::Id n = cubicSpline.GetControlPoints().GetNumberOfValues(); - vtkm::FloatDefault t = x[0]; - vtkm::FloatDefault tEnd = x[x.size() - 1]; - - vtkm::FloatDefault dt = (tEnd - t) / (n * 100); - std::vector params; - while (t < tEnd) - { - params.push_back(t); - t += dt; - } - - auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); - - vtkm::cont::Invoker invoke; - vtkm::cont::ArrayHandle result; - invoke(SplineEvalWorklet{}, paramsAH, cubicSpline, result); - - std::cout << "SaveSamples" << std::endl; - std::ofstream fout("output.txt"), ptsOut("pts.txt"); - fout << "X,Y" << std::endl; - for (vtkm::Id i = 0; i < paramsAH.GetNumberOfValues(); i++) - { - fout << paramsAH.ReadPortal().Get(i) << "," << result.ReadPortal().Get(i) << std::endl; - } - - ptsOut << "X,Y" << std::endl; - auto pts = cubicSpline.GetControlPoints().ReadPortal(); - auto vals = cubicSpline.GetValues().ReadPortal(); - for (vtkm::Id i = 0; i < pts.GetNumberOfValues(); i++) - ptsOut << pts.Get(i) << ", " << vals.Get(i) << std::endl; -} - -void DoTest(const std::vector& x, - const std::vector& y, - const std::vector& xSamples, - const std::vector& vals) -{ - /* - vtkm::cont::CubicSpline cubicSpline; - cubicSpline.SetControlPoints(vtkm::cont::make_ArrayHandle(x, vtkm::CopyFlag::Off)); - cubicSpline.SetValues(vtkm::cont::make_ArrayHandle(y, vtkm::CopyFlag::Off)); - cubicSpline.Update(); - - //Make sure we the spline interpolates the control points. - CheckEvaluation(cubicSpline, x, y); - - //Make sure the spline evaluates to the correct values at the given points. - CheckEvaluation(cubicSpline, xSamples, vals); - */ -} - void CubicHermiteSplineTest() { std::vector pts = { { 0, 0, 0 }, { 1, 1, 0 }, { 2, 1, 0 }, { 3, -.5, 0 }, @@ -157,28 +95,6 @@ void CubicHermiteSplineTest() invoke( SplineEvalWorklet{}, vtkm::cont::make_ArrayHandle(vals, vtkm::CopyFlag::Off), spline0, result); - /* - vtkm::cont::CubicHermiteSpline spline(pts); - - vtkm::Id n = spline.GetData().GetNumberOfValues(); - auto t0 = spline.GetKnots().ReadPortal().Get(0); - auto t1 = spline.GetKnots().ReadPortal().Get(n-1); - - vtkm::FloatDefault dt = (t1 - t0) / (n * 100); - std::vector params; - vtkm::FloatDefault t = t0; - - while (t < t1) - { - params.push_back(t); - t += dt; - } - - vtkm::cont::Invoker invoke; - vtkm::cont::ArrayHandle result; - invoke(SplineEvalWorklet{}, vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off), spline, result); - */ - std::cout << "SaveSamples" << std::endl; std::ofstream fout("/Users/dpn/output.txt"), ptsOut("/Users/dpn/pts.txt"); fout << "T, X, Y, Z" << std::endl; @@ -194,29 +110,9 @@ void CubicHermiteSplineTest() << std::endl; } -void CubicSplineTest() -{ - CubicHermiteSplineTest(); - - //uniform spacing wavy function - std::vector xVals = { 0.0f, 1.0f, 2.0f, 3.0f, 4.0f }; - std::vector yVals = { 0.0f, 1.0f, -1.0f, 1.0f, 0.0f }; - std::vector samples = { 0.6, 1.4, 2.16, 3.51198 }; - std::vector res = { 1.03886, 0.110853, -0.890431, 0.91292 }; - - //DoTest(xVals, yVals, samples, res); - - //non-uniform spacing wavy function - xVals = { 0.0f, 1.0f, 1.1f, 3.9f, 4.0f, 5.0f }; - yVals = { 0.0f, 1.0f, 1.2f, 1.0f, 0.0f, 0.0f }; - samples = { 0.7f, 1.45f, 3.65f, 4.38f }; - res = { 0.548318, 2.15815f, 3.13104f, -1.77143 }; - //DoTest(xVals, yVals, samples, res); -} - } // anonymous namespace -int UnitTestCubicSpline(int argc, char* argv[]) +int UnitTestCubicHermiteSpline(int argc, char* argv[]) { - return vtkm::cont::testing::Testing::Run(CubicSplineTest, argc, argv); + return vtkm::cont::testing::Testing::Run(CubicHermiteSplineTest, argc, argv); } -- GitLab From 6bdfa46df7f16774701f74b88f888ad1e1809cec Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Fri, 4 Apr 2025 16:00:58 -0400 Subject: [PATCH 11/13] Code cleanup..... --- vtkm/cont/CubicHermiteSpline.cxx | 24 +---- .../testing/UnitTestCubicHermiteSpline.cxx | 99 ++++++++++--------- 2 files changed, 58 insertions(+), 65 deletions(-) diff --git a/vtkm/cont/CubicHermiteSpline.cxx b/vtkm/cont/CubicHermiteSpline.cxx index b3ad06be28..3750f335b5 100644 --- a/vtkm/cont/CubicHermiteSpline.cxx +++ b/vtkm/cont/CubicHermiteSpline.cxx @@ -8,8 +8,10 @@ // PURPOSE. See the above copyright notice for more information. //============================================================================ +#include #include #include +#include #include #include #include @@ -37,21 +39,6 @@ struct CalcNeighborDistanceWorklet : public vtkm::worklet::WorkletMapField } }; -struct DivideBy : public vtkm::worklet::WorkletMapField -{ - DivideBy(const vtkm::FloatDefault& divisor) - : Divisor(divisor) - { - } - - using ControlSignature = void(FieldInOut); - using ExecutionSignature = void(_1); - - VTKM_EXEC void operator()(vtkm::FloatDefault& val) const { val = val / this->Divisor; } - - vtkm::FloatDefault Divisor; -}; - struct CalcTangentsWorklet : public vtkm::worklet::WorkletMapField { CalcTangentsWorklet(const vtkm::Id& numPoints) @@ -89,7 +76,6 @@ struct CalcTangentsWorklet : public vtkm::worklet::WorkletMapField auto dT = knots.Get(idx1) - knots.Get(idx0); tangent = dX / dT; - std::cout << "Tan: " << idx << " = " << tangent << " dx: " << dX << " " << dT << std::endl; } vtkm::Id NumPoints; @@ -117,9 +103,8 @@ VTKM_CONT void CubicHermiteSpline::ComputeKnots() if (sum == 0.0) throw std::invalid_argument("Error: accumulated distance between data is zero."); - invoker(DivideBy{ sum }, this->Knots); - std::cout << ":: params= "; - vtkm::cont::printSummary_ArrayHandle(this->Knots, std::cout); + 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() @@ -149,6 +134,5 @@ vtkm::exec::CubicHermiteSpline CubicHermiteSpline::PrepareForExecution( return ExecObjType(this->Data, this->Knots, this->Tangents, device, token); } - } } // namespace vtkm::cont diff --git a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx index 32ec21723a..86c15d7010 100644 --- a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx @@ -40,33 +40,66 @@ public: } }; -template -void CheckEvaluation(const SplineType& spline, - const std::vector& params, - const std::vector& answer) +void CheckEvaluation(const vtkm::cont::CubicHermiteSpline& spline, + const vtkm::cont::ArrayHandle& params, + const std::vector& answer) { - auto paramsAH = vtkm::cont::make_ArrayHandle(params, vtkm::CopyFlag::Off); - vtkm::cont::Invoker invoke; - vtkm::cont::ArrayHandle result; - - invoke(SplineEvalWorklet{}, paramsAH, spline, result); + 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, 0 }, { 2, 1, 0 }, { 3, -.5, 0 }, { 4, -1, 0 }, { 5, -1, 0 }, { 6, 0, 0 } }; - std::vector knots; - pts.clear(); - - vtkm::Id n = 100; + vtkm::cont::CubicHermiteSpline spline(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.20072, 1.07379, 0 }, + { 2.64681, 0.012108, 0 }, + { 2.79356, -0.239254, 0 }, + { 5.99081, -0.00924034, 0 } }; + CheckEvaluation(spline, params, result); + + //Explicitly set knots and check. + std::vector knots = { 0, 1, 2, 3, 4, 5, 6 }; + spline = vtkm::cont::CubicHermiteSpline(pts, knots); + CheckEvaluation(spline, knots, pts); + + //Non-uniform knots. + knots = { 0, 1, 2, 2.1, 2.2, 2.3, 3 }; + spline = vtkm::cont::CubicHermiteSpline(pts, knots); + CheckEvaluation(spline, knots, pts); + + params = { 1.5, 2.05, 2.11, 2.299, 2.8 }; + result = { { 1.39773, 1.23295, 0 }, + { 2.39773, 0.357954, 0 }, + { 3.1, -0.59275, 0 }, + { 4.99735, -1.00125, 0 }, + { 5.75802, -0.293003, 0 } }; + 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); @@ -76,38 +109,14 @@ void CubicHermiteSplineTest() knots.push_back(t); t += dt; } - vtkm::cont::CubicHermiteSpline spline0(pts, knots); - auto parametricRange = spline0.GetParametricRange(); - - t = parametricRange.Min; - dt = parametricRange.Length() / 1500; - std::vector vals; - while (t <= parametricRange.Max) - { - vals.push_back(t); - t += dt; - } - std::cout << "KNOTS= "; - vtkm::cont::printSummary_ArrayHandle(spline0.GetKnots(), std::cout); - - vtkm::cont::Invoker invoke; - vtkm::cont::ArrayHandle result; - invoke( - SplineEvalWorklet{}, vtkm::cont::make_ArrayHandle(vals, vtkm::CopyFlag::Off), spline0, result); - - std::cout << "SaveSamples" << std::endl; - std::ofstream fout("/Users/dpn/output.txt"), ptsOut("/Users/dpn/pts.txt"); - fout << "T, X, Y, Z" << std::endl; - ptsOut << "T, X, Y, Z" << std::endl; - auto knots_ = spline0.GetKnots().ReadPortal(); - for (vtkm::Id i = 0; i < pts.size(); i++) - ptsOut << knots_.Get(i) << ", " << pts[i][0] << ", " << pts[i][1] << ", " << pts[i][2] - << std::endl; - - auto res = result.ReadPortal(); - for (vtkm::Id i = 0; i < res.GetNumberOfValues(); i++) - fout << vals[i] << ", " << res.Get(i)[0] << ", " << res.Get(i)[1] << ", " << res.Get(i)[2] - << std::endl; + spline = vtkm::cont::CubicHermiteSpline(pts, knots); + CheckEvaluation(spline, knots, pts); + + params = { 0.15, 2.38, 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 -- GitLab From f84ea4a64948526262aeb51d651a786dc02534d4 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Mon, 7 Apr 2025 07:58:56 -0400 Subject: [PATCH 12/13] Increase the coverage of the unit tests. --- .../testing/UnitTestCubicHermiteSpline.cxx | 34 ++++++++++++------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx index 86c15d7010..d9ffe90b72 100644 --- a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx @@ -61,8 +61,8 @@ void CheckEvaluation(const vtkm::cont::CubicHermiteSpline& spline, void CubicHermiteSplineTest() { - std::vector pts = { { 0, 0, 0 }, { 1, 1, 0 }, { 2, 1, 0 }, { 3, -.5, 0 }, - { 4, -1, 0 }, { 5, -1, 0 }, { 6, 0, 0 } }; + 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(pts); //Evaluation at knots gives the sample pts. @@ -70,10 +70,10 @@ void CubicHermiteSplineTest() //Evaluation at non-knot values. std::vector params = { 0.21, 0.465, 0.501, 0.99832 }; - std::vector result = { { 1.20072, 1.07379, 0 }, - { 2.64681, 0.012108, 0 }, - { 2.79356, -0.239254, 0 }, - { 5.99081, -0.00924034, 0 } }; + 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. @@ -81,17 +81,24 @@ void CubicHermiteSplineTest() spline = vtkm::cont::CubicHermiteSpline(pts, 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(pts, knots); CheckEvaluation(spline, knots, pts); params = { 1.5, 2.05, 2.11, 2.299, 2.8 }; - result = { { 1.39773, 1.23295, 0 }, - { 2.39773, 0.357954, 0 }, - { 3.1, -0.59275, 0 }, - { 4.99735, -1.00125, 0 }, - { 5.75802, -0.293003, 0 } }; + 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. @@ -104,7 +111,7 @@ void CubicHermiteSplineTest() { vtkm::FloatDefault x = vtkm::Cos(t); vtkm::FloatDefault y = vtkm::Sin(t); - vtkm::FloatDefault z = vtkm::Cos(t) * vtkm::Sin(t); + vtkm::FloatDefault z = x * y; pts.push_back({ x, y, z }); knots.push_back(t); t += dt; @@ -112,7 +119,8 @@ void CubicHermiteSplineTest() spline = vtkm::cont::CubicHermiteSpline(pts, knots); CheckEvaluation(spline, knots, pts); - params = { 0.15, 2.38, 4.92, 6.2 }; + //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) }); -- GitLab From 4ad59d1e6f514617685bbca64e977faf1c832964 Mon Sep 17 00:00:00 2001 From: David Pugmire Date: Mon, 7 Apr 2025 10:59:01 -0400 Subject: [PATCH 13/13] Update the interface for using. --- vtkm/cont/CubicHermiteSpline.cxx | 26 ++++++--- vtkm/cont/CubicHermiteSpline.h | 53 +++++++++---------- .../testing/UnitTestCubicHermiteSpline.cxx | 17 ++++-- 3 files changed, 54 insertions(+), 42 deletions(-) diff --git a/vtkm/cont/CubicHermiteSpline.cxx b/vtkm/cont/CubicHermiteSpline.cxx index 3750f335b5..d8e739ab73 100644 --- a/vtkm/cont/CubicHermiteSpline.cxx +++ b/vtkm/cont/CubicHermiteSpline.cxx @@ -82,11 +82,17 @@ struct CalcTangentsWorklet : public vtkm::worklet::WorkletMapField }; } //anonymous namespace -VTKM_CONT vtkm::Range CubicHermiteSpline::GetParametricRange() const +VTKM_CONT vtkm::Range CubicHermiteSpline::GetParametricRange() { - auto portal = this->Knots.ReadPortal(); - auto n = portal.GetNumberOfValues(); - return { portal.Get(0), portal.Get(n - 1) }; + 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() @@ -119,16 +125,20 @@ VTKM_CONT void CubicHermiteSpline::ComputeTangents() vtkm::exec::CubicHermiteSpline CubicHermiteSpline::PrepareForExecution( vtkm::cont::DeviceAdapterId device, - vtkm::cont::Token& token) const + vtkm::cont::Token& token) { vtkm::Id n = this->Data.GetNumberOfValues(); if (n < 2) - throw std::invalid_argument("At least two points are required for spline interpolation."); + 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 std::invalid_argument("Number of data points must match the number of knots."); + throw vtkm::cont::ErrorBadValue("Number of data points must match the number of knots."); if (n != this->Tangents.GetNumberOfValues()) - throw std::invalid_argument("Number of data points must match the number of tangents."); + throw vtkm::cont::ErrorBadValue("Number of data points must match the number of tangents."); using ExecObjType = vtkm::exec::CubicHermiteSpline; diff --git a/vtkm/cont/CubicHermiteSpline.h b/vtkm/cont/CubicHermiteSpline.h index a45febfef5..92f3cb114c 100644 --- a/vtkm/cont/CubicHermiteSpline.h +++ b/vtkm/cont/CubicHermiteSpline.h @@ -30,55 +30,51 @@ public: CubicHermiteSpline() = default; virtual ~CubicHermiteSpline() = default; - CubicHermiteSpline(const std::vector& data, vtkm::CopyFlag copy = vtkm::CopyFlag::On) - : Data(vtkm::cont::make_ArrayHandle(data, copy)) + 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->ComputeKnots(); - this->ComputeTangents(); + this->Data = vtkm::cont::make_ArrayHandle(data, copy); } - CubicHermiteSpline(const std::vector& data, - const std::vector& tangents, - vtkm::CopyFlag copy = vtkm::CopyFlag::On) - : Data(vtkm::cont::make_ArrayHandle(data, copy)) - , Tangents(vtkm::cont::make_ArrayHandle(tangents, copy)) + VTKM_CONT void SetKnots(const vtkm::cont::ArrayHandle& knots) { - this->ComputeKnots(); + this->Knots = knots; } - - CubicHermiteSpline(const std::vector& data, - const std::vector& knots, - vtkm::CopyFlag copy = vtkm::CopyFlag::On) - : Data(vtkm::cont::make_ArrayHandle(data, copy)) - , Knots(vtkm::cont::make_ArrayHandle(knots, copy)) + VTKM_CONT void SetKnots(const std::vector& knots, + vtkm::CopyFlag copy = vtkm::CopyFlag::On) { - this->ComputeTangents(); + this->Knots = vtkm::cont::make_ArrayHandle(knots, copy); } - CubicHermiteSpline(const std::vector& data, - const std::vector& tangents, - const std::vector& knots, - vtkm::CopyFlag copy = vtkm::CopyFlag::On) - : Data(vtkm::cont::make_ArrayHandle(data, copy)) - , Knots(vtkm::cont::make_ArrayHandle(knots, copy)) - , Tangents(vtkm::cont::make_ArrayHandle(tangents, 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() const; + VTKM_CONT vtkm::Range GetParametricRange(); VTKM_CONT const vtkm::cont::ArrayHandle& GetData() const { return this->Data; } - VTKM_CONT const vtkm::cont::ArrayHandle& GetTangents() const + VTKM_CONT const vtkm::cont::ArrayHandle& GetTangents() { + if (this->Tangents.GetNumberOfValues() == 0) + this->ComputeTangents(); return this->Tangents; } - VTKM_CONT const vtkm::cont::ArrayHandle& GetKnots() const + 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) const; + vtkm::cont::Token& token); private: VTKM_CONT void ComputeKnots(); @@ -88,7 +84,6 @@ private: vtkm::cont::ArrayHandle Knots; vtkm::cont::ArrayHandle Tangents; }; - } } // vtkm::cont diff --git a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx index d9ffe90b72..4d145c6a21 100644 --- a/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx +++ b/vtkm/cont/testing/UnitTestCubicHermiteSpline.cxx @@ -36,7 +36,7 @@ public: { auto res = spline.Evaluate(param, value); if (res != vtkm::ErrorCode::Success) - throw vtkm::cont::ErrorBadValue("Spline evaluation failed."); + this->RaiseError("Spline evaluation failed."); } }; @@ -64,7 +64,8 @@ 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(pts); + vtkm::cont::CubicHermiteSpline spline; + spline.SetData(pts); //Evaluation at knots gives the sample pts. CheckEvaluation(spline, spline.GetKnots(), pts); @@ -78,7 +79,9 @@ void CubicHermiteSplineTest() //Explicitly set knots and check. std::vector knots = { 0, 1, 2, 3, 4, 5, 6 }; - spline = vtkm::cont::CubicHermiteSpline(pts, knots); + spline = vtkm::cont::CubicHermiteSpline(); + spline.SetData(pts); + spline.SetKnots(knots); CheckEvaluation(spline, knots, pts); //Evaluation at non-knot values. @@ -90,7 +93,9 @@ void CubicHermiteSplineTest() //Non-uniform knots. knots = { 0, 1, 2, 2.1, 2.2, 2.3, 3 }; - spline = vtkm::cont::CubicHermiteSpline(pts, knots); + 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 }; @@ -116,7 +121,9 @@ void CubicHermiteSplineTest() knots.push_back(t); t += dt; } - spline = vtkm::cont::CubicHermiteSpline(pts, knots); + spline = vtkm::cont::CubicHermiteSpline(); + spline.SetData(pts); + spline.SetKnots(knots); CheckEvaluation(spline, knots, pts); //Evaluate at a few points and check against analytical results. -- GitLab