// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.08.13
#pragma once
#include <Mathematics/Logger.h>
#include <Mathematics/Math.h>
#include <Mathematics/Polynomial1.h>
#include <array>
// IntpBSplineUniform is the class for B-spline interpolation of uniformly
// spaced N-dimensional data. The algorithm is described in
// https://www.geometrictools.com/Documentation/BSplineInterpolation.pdf
// A sample application for this topic is
// GeometricTools/GTEngine/Samples/Imagics/BSplineInterpolation
//
// The Controls adapter allows access to your control points without regard
// to how you organize your data. You can even defer the computation of a
// control point until it is needed via the operator()(...) calls that
// Controls must provide, and you can cache the points according to your own
// needs. The minimal interface for a Controls adapter is
//
// struct Controls
// {
// // The control_point_type is of your choosing. It must support
// // assignment, scalar multiplication and addition. Specifically, if
// // C0, C1 and C2 are control points and s is a scalar, the interpolator
// // needs to perform operations
// // C1 = C0;
// // C1 = C0 * s;
// // C2 = C0 + C1;
// typedef control_point_type Type;
//
// // The number of elements in the specified dimension.
// int GetSize(int dimension) const;
//
// // Get a control point based on an n-tuple lookup. The interpolator
// // does not need to know your organization; all it needs is the
// // desired control point. The 'tuple' input must have n elements.
// Type operator() (int const* tuple) const;
//
// // If you choose to use the specialized interpolators for dimensions
// // 1, 2 or 3, you must also provide the following accessor, where
// // the input is an n-tuple listed component-by-component (1, 2 or 3
// // components).
// Type operator() (int i0, int i1, ..., int inm1) const;
// }
namespace WwiseGTE
{
template <typename Real, typename Controls>
class IntpBSplineUniformShared
{
protected:
// Abstract base class construction. A virtual destructor is not
// provided because there are no required side effects in the base
// class when destroying objects from the derived classes.
IntpBSplineUniformShared(int numDimensions, int const* degrees,
Controls const& controls, typename Controls::Type ctZero, int cacheMode)
:
mNumDimensions(numDimensions),
mDegree(numDimensions),
mControls(&controls),
mCTZero(ctZero),
mCacheMode(cacheMode),
mNumLocalControls(0),
mDegreeP1(numDimensions),
mNumControls(numDimensions),
mTMin(numDimensions),
mTMax(numDimensions),
mBlender(numDimensions),
mDCoefficient(numDimensions),
mLMax(numDimensions),
mPowerDSDT(numDimensions),
mITuple(numDimensions),
mJTuple(numDimensions),
mKTuple(numDimensions),
mLTuple(numDimensions),
mSumIJTuple(numDimensions),
mUTuple(numDimensions),
mPTuple(numDimensions)
{
// The condition c+1 > d+1 is required so that when s = c+1-d, its
// maximum value, we have at least two s-knots (d and d + 1).
for (int dim = 0; dim < mNumDimensions; ++dim)
{
if (mControls->GetSize(dim) <= degrees[dim] + 1)
{
LogError("Incompatible degree and number of controls.");
}
}
mNumLocalControls = 1;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mDegree[dim] = degrees[dim];
mDegreeP1[dim] = degrees[dim] + 1;
mNumLocalControls *= mDegreeP1[dim];
mNumControls[dim] = controls.GetSize(dim);
mTMin[dim] = (Real)-0.5;
mTMax[dim] = static_cast<Real>(mNumControls[dim]) - (Real)0.5;
ComputeBlendingMatrix(mDegree[dim], mBlender[dim]);
ComputeDCoefficients(mDegree[dim], mDCoefficient[dim], mLMax[dim]);
ComputePowers(mDegree[dim], mNumControls[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim]);
}
if (mCacheMode == NO_CACHING)
{
mPhi.resize(mNumDimensions);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mPhi[dim].resize(mDegreeP1[dim]);
}
}
else
{
InitializeTensors();
}
}
#if !defined(GTE_INTP_BSPLINE_UNIFORM_NO_SPECIALIZATION)
IntpBSplineUniformShared()
:
mNumDimensions(0),
mControls(nullptr),
mCTZero(),
mCacheMode(0),
mNumLocalControls(0)
{
}
#endif
public:
// Support for caching the intermediate tensor product of control
// points with the blending matrices. A precached container has all
// elements precomputed before any Evaluate(...) calls. The 'bool'
// flags are all set to 'true'. A cached container fills the elements
// on demand. The 'bool' flags are initially 'false', indicating the
// EvalType component has not yet been computed. After it is computed
// and stored, the flag is set to 'true'.
enum
{
NO_CACHING,
PRE_CACHING,
ON_DEMAND_CACHING
};
// Member access.
inline int GetDegree(int dim) const
{
return mDegree[dim];
}
inline int GetNumControls(int dim) const
{
return mNumControls[dim];
}
inline Real GetTMin(int dim) const
{
return mTMin[dim];
}
inline Real GetTMax(int dim) const
{
return mTMax[dim];
}
inline int GetCacheMode() const
{
return mCacheMode;
}
protected:
// Disallow copying and moving.
IntpBSplineUniformShared(IntpBSplineUniformShared const&) = delete;
IntpBSplineUniformShared& operator=(IntpBSplineUniformShared const&) = delete;
IntpBSplineUniformShared(IntpBSplineUniformShared&&) = delete;
IntpBSplineUniformShared& operator=(IntpBSplineUniformShared&&) = delete;
// ComputeBlendingMatrix, ComputeDCoefficients, ComputePowers and
// GetKey are used by the general-dimension derived classes and by
// the specializations for dimensions 1, 2 and 3.
// Compute the blending matrix that combines the control points and
// the polynomial vector. The matrix A is stored in row-major order.
static void ComputeBlendingMatrix(int degree, std::vector<Real>& A)
{
int const degreeP1 = degree + 1;
A.resize(degreeP1 * degreeP1);
if (degree == 0)
{
A[0] = (Real)1;
return;
}
// P_{0,0}(s)
std::vector<Polynomial1<Real>> P(degreeP1);
P[0][0] = (Real)1;
// L0 = s/j
Polynomial1<Real> L0(1);
L0[0] = (Real)0;
// L1(s) = (j + 1 - s)/j
Polynomial1<Real> L1(1);
// s-1 is used in computing translated P_{j-1,k-1}(s-1)
Polynomial1<Real> sm1 = { (Real)-1, (Real)1 };
// Compute
// P_{j,k}(s) = L0(s)*P_{j-1,k}(s) + L1(s)*P_{j-1,k-1}(s-1)
// for 0 <= k <= j where 1 <= j <= degree. When k = 0,
// P_{j-1,-1}(s) = 0, so P_{j,0}(s) = L0(s)*P_{j-1,0}(s). When
// k = j, P_{j-1,j}(s) = 0, so P_{j,j}(s) = L1(s)*P_{j-1,j-1}(s).
// The polynomials at level j-1 are currently stored in P[0]
// through P[j-1]. The polynomials at level j are computed and
// stored in P[0] through P[j]; that is, they are computed in
// place to reduce memory usage and copying. This requires
// computing P[k] (level j) from P[k] (level j-1) and P[k-1]
// (level j-1), which means we have to process k = j down to
// k = 0.
for (int j = 1; j <= degree; ++j)
{
Real invJ = (Real)1 / (Real)j;
L0[1] = invJ;
L1[0] = (Real)1 + invJ;
L1[1] = -invJ;
for (int k = j; k >= 0; --k)
{
Polynomial1<Real> result = { (Real)0 };
if (k > 0)
{
result += L1 * P[k - 1].GetTranslation((Real)1);
}
if (k < j)
{
result += L0 * P[k];
}
P[k] = result;
}
}
// Compute Q_{d,k}(s) = P_{d,k}(s + k).
std::vector<Polynomial1<Real>> Q(degreeP1);
for (int k = 0; k <= degree; ++k)
{
Q[k] = P[k].GetTranslation(static_cast<Real>(-k));
}
// Extract the matrix A from the Q-polynomials. Row r of A
// contains the coefficients of Q_{d,d-r}(s).
for (int k = 0, row = degree; k <= degree; ++k, --row)
{
for (int col = 0; col <= degree; ++col)
{
A[col + degreeP1 * row] = Q[k][col];
}
}
}
// Compute the coefficients for the derivative polynomial terms.
static void ComputeDCoefficients(int degree, std::vector<Real>& dCoefficients,
std::vector<int>& ellMax)
{
int numDCoefficients = (degree + 1) * (degree + 2) / 2;
dCoefficients.resize(numDCoefficients);
for (int i = 0; i < numDCoefficients; ++i)
{
dCoefficients[i] = 1;
}
for (int order = 1, col0 = 0, col1 = degree + 1; order <= degree; ++order)
{
++col0;
for (int c = order, m = 1; c <= degree; ++c, ++m, ++col0, ++col1)
{
dCoefficients[col1] = dCoefficients[col0] * m;
}
}
ellMax.resize(degree + 1);
ellMax[0] = degree;
for (int i0 = 0, i1 = 1; i1 <= degree; i0 = i1++)
{
ellMax[i1] = ellMax[i0] + degree - i0;
}
}
// Compute powers of ds/dt.
void ComputePowers(int degree, int numControls, Real tmin, Real tmax,
std::vector<Real>& powerDSDT)
{
Real dsdt = static_cast<Real>(numControls - degree) / (tmax - tmin);
powerDSDT.resize(degree + 1);
powerDSDT[0] = (Real)1;
powerDSDT[1] = dsdt;
for (int i = 2; i <= degree; ++i)
{
powerDSDT[i] = powerDSDT[i - 1] * dsdt;
}
}
// Determine the interval [index,index+1) corresponding to the
// specified value of t and compute u in that interval.
static void GetKey(Real t, Real tmin, Real tmax, Real dsdt, int numControls,
int degree, int& index, Real& u)
{
// Compute s - d = ((c + 1 - d)/(c + 1))(t + 1/2), the index for
// which d + index <= s < d + index + 1. Let u = s - d - index so
// that 0 <= u < 1.
if (t > tmin)
{
if (t < tmax)
{
Real smd = dsdt * (t - tmin);
index = static_cast<uint32_t>(std::floor(smd));
u = smd - static_cast<float>(index);
}
else
{
// In the evaluation, s = c + 1 - d and i = c - d. This
// causes s-d-i to be 1 in G_c(c+1-d). Effectively, the
// selection of i extends the s-domain [d,c+1) to its
// support [d,c+1].
index = numControls - 1 - degree;
u = (Real)1;
}
}
else
{
index = 0;
u = (Real)0;
}
}
// The remaining functions are used only by the general-dimension
// derived classes when caching is enabled.
// For the multidimensional tensor Phi{iTuple, kTuple), compute the
// portion of the 1-dimensional index that corresponds to iTuple.
int GetRowIndex(std::vector<int> const& i) const
{
int rowIndex = i[mNumDimensions - 1];
int j1 = 2 * mNumDimensions - 2;
for (int j0 = mNumDimensions - 2; j0 >= 0; --j0, --j1)
{
rowIndex = mTBound[j1] * rowIndex + i[j0];
}
rowIndex = mTBound[j1] * rowIndex;
return rowIndex;
}
// For the multidimensional tensor Phi{iTuple, kTuple), combine the
// GetRowIndex(...) output with kTuple to produce the full
// 1-dimensional index.
int GetIndex(int rowIndex, std::vector<int> const& k) const
{
int index = rowIndex + k[mNumDimensions - 1];
for (int j = mNumDimensions - 2; j >= 0; --j)
{
index = mTBound[j] * index + k[j];
}
return index;
}
// Compute Phi(iTuple, kTuple). The 'index' value is an already
// computed 1-dimensional index for the tensor.
void ComputeTensor(int const* i, int const* k, int index)
{
auto element = mCTZero;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mComputeJTuple[dim] = 0;
}
for (int iterate = 0; iterate < mNumLocalControls; ++iterate)
{
Real blend(1);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
blend *= mBlender[dim][k[dim] + mDegreeP1[dim] * mComputeJTuple[dim]];
mComputeSumIJTuple[dim] = i[dim] + mComputeJTuple[dim];
}
element = element + (*mControls)(mComputeSumIJTuple.data()) * blend;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
if (++mComputeJTuple[dim] < mDegreeP1[dim])
{
break;
}
mComputeJTuple[dim] = 0;
}
}
mTensor[index] = element;
}
// Allocate the containers used for caching and fill in the tensor
// for precaching when that mode is selected.
void InitializeTensors()
{
mTBound.resize(2 * mNumDimensions);
mComputeJTuple.resize(mNumDimensions);
mComputeSumIJTuple.resize(mNumDimensions);
mDegreeMinusOrder.resize(mNumDimensions);
mTerm.resize(mNumDimensions);
int current = 0;
int numCached = 1;
for (int dim = 0; dim < mNumDimensions; ++dim, ++current)
{
mTBound[current] = mDegreeP1[dim];
numCached *= mTBound[current];
}
for (int dim = 0; dim < mNumDimensions; ++dim, ++current)
{
mTBound[current] = mNumControls[dim] - mDegree[dim];
numCached *= mTBound[current];
}
mTensor.resize(numCached);
mCached.resize(numCached);
if (mCacheMode == PRE_CACHING)
{
std::vector<int> tuple(2 * mNumDimensions, 0);
for (int index = 0; index < numCached; ++index)
{
ComputeTensor(&tuple[mNumDimensions], &tuple[0], index);
for (int i = 0; i < 2 * mNumDimensions; ++i)
{
if (++tuple[i] < mTBound[i])
{
break;
}
tuple[i] = 0;
}
}
std::fill(mCached.begin(), mCached.end(), true);
}
else
{
std::fill(mCached.begin(), mCached.end(), false);
}
}
// Evaluate the interpolator. Each element of 'order' indicates the
// order of the derivative you want to compute. For the function
// value itself, pass in 'order' that has all 0 elements.
typename Controls::Type EvaluateNoCaching(int const* order, Real const* t)
{
auto result = mCTZero;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
if (order[dim] < 0 || order[dim] > mDegree[dim])
{
return result;
}
}
for (int dim = 0; dim < mNumDimensions; ++dim)
{
GetKey(t[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim][1],
mNumControls[dim], mDegree[dim], mITuple[dim], mUTuple[dim]);
}
for (int dim = 0; dim < mNumDimensions; ++dim)
{
int jIndex = 0;
for (int j = 0; j <= mDegree[dim]; ++j)
{
int kjIndex = mDegree[dim] + jIndex;
int ell = mLMax[dim][order[dim]];
mPhi[dim][j] = (Real)0;
for (int k = mDegree[dim]; k >= order[dim]; --k)
{
mPhi[dim][j] = mPhi[dim][j] * mUTuple[dim] +
mBlender[dim][kjIndex--] * mDCoefficient[dim][ell--];
}
jIndex += mDegreeP1[dim];
}
}
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mJTuple[dim] = 0;
mSumIJTuple[dim] = mITuple[dim];
mPTuple[dim] = mPhi[dim][0];
}
for (int iterate = 0; iterate < mNumLocalControls; ++iterate)
{
Real product(1);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
product *= mPTuple[dim];
}
result = result + (*mControls)(mSumIJTuple.data()) * product;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
if (++mJTuple[dim] <= mDegree[dim])
{
mSumIJTuple[dim] = mITuple[dim] + mJTuple[dim];
mPTuple[dim] = mPhi[dim][mJTuple[dim]];
break;
}
mJTuple[dim] = 0;
mSumIJTuple[dim] = mITuple[dim];
mPTuple[dim] = mPhi[dim][0];
}
}
Real adjust(1);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
adjust *= mPowerDSDT[dim][order[dim]];
}
result = result * adjust;
return result;
}
typename Controls::Type EvaluateCaching(int const* order, Real const* t)
{
int numIterates = 1;
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mDegreeMinusOrder[dim] = mDegree[dim] - order[dim];
if (mDegreeMinusOrder[dim] < 0 || mDegreeMinusOrder[dim] > mDegree[dim])
{
return mCTZero;
}
numIterates *= mDegreeMinusOrder[dim] + 1;
}
for (int dim = 0; dim < mNumDimensions; ++dim)
{
GetKey(t[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim][1],
mNumControls[dim], mDegree[dim], mITuple[dim], mUTuple[dim]);
}
int rowIndex = GetRowIndex(mITuple);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
mJTuple[dim] = 0;
mKTuple[dim] = mDegree[dim];
mLTuple[dim] = mLMax[dim][order[dim]];
mTerm[dim] = mCTZero;
}
for (int iterate = 0; iterate < numIterates; ++iterate)
{
int index = GetIndex(rowIndex, mKTuple);
if (mCacheMode == ON_DEMAND_CACHING && !mCached[index])
{
ComputeTensor(mITuple.data(), mKTuple.data(), index);
mCached[index] = true;
}
mTerm[0] = mTerm[0] * mUTuple[0] + mTensor[index] * mDCoefficient[0][mLTuple[0]];
for (int dim = 0; dim < mNumDimensions; ++dim)
{
if (++mJTuple[dim] <= mDegreeMinusOrder[dim])
{
--mKTuple[dim];
--mLTuple[dim];
break;
}
int dimp1 = dim + 1;
if (dimp1 < mNumDimensions)
{
mTerm[dimp1] = mTerm[dimp1] * mUTuple[dimp1] + mTerm[dim] * mDCoefficient[dimp1][mLTuple[dimp1]];
mTerm[dim] = mCTZero;
mJTuple[dim] = 0;
mKTuple[dim] = mDegree[dim];
mLTuple[dim] = mLMax[dim][order[dim]];
}
}
}
auto result = mTerm[mNumDimensions - 1];
Real adjust(1);
for (int dim = 0; dim < mNumDimensions; ++dim)
{
adjust *= mPowerDSDT[dim][order[dim]];
}
result = result * adjust;
return result;
}
// Constructor inputs.
int const mNumDimensions; // N
std::vector<int> mDegree; // degree[N]
Controls const* mControls;
typename Controls::Type const mCTZero;
int const mCacheMode;
// Parameters for B-spline evaluation. All std::vector containers
// have N elements.
int mNumLocalControls; // product of (degree[]+1)
std::vector<int> mDegreeP1;
std::vector<int> mNumControls;
std::vector<Real> mTMin;
std::vector<Real> mTMax;
std::vector<std::vector<Real>> mBlender;
std::vector<std::vector<Real>> mDCoefficient;
std::vector<std::vector<int>> mLMax;
std::vector<std::vector<Real>> mPowerDSDT;
std::vector<int> mITuple;
std::vector<int> mJTuple;
std::vector<int> mKTuple;
std::vector<int> mLTuple;
std::vector<int> mSumIJTuple;
std::vector<Real> mUTuple;
std::vector<Real> mPTuple;
// Support for no-cached B-spline evaluation. The std::vector
// container has N elements.
std::vector<std::vector<Real>> mPhi;
// Support for cached B-spline evaluation.
std::vector<int> mTBound; // tbound[2*N]
std::vector<int> mComputeJTuple; // computejtuple[N]
std::vector<int> mComputeSumIJTuple; // computesumijtuple[N]
std::vector<int> mDegreeMinusOrder; // degreeminusorder[N]
std::vector<typename Controls::Type> mTerm; // mTerm[N]
std::vector<typename Controls::Type> mTensor; // depends on numcontrols
std::vector<bool> mCached; // same size as mTensor
};
}
// Implementation for B-spline interpolation whose dimension is known at
// compile time.
namespace WwiseGTE
{
template <typename Real, typename Controls, int N = 0>
class IntpBSplineUniform : public IntpBSplineUniformShared<Real, Controls>
{
public:
// The caller is responsible for ensuring that the IntpBSplineUniform
// object persists as long as the input 'controls' exists.
IntpBSplineUniform(std::array<int, N> const& degrees, Controls const& controls,
typename Controls::Type ctZero, int cacheMode)
:
IntpBSplineUniformShared<Real, Controls>(N, degrees.data(), controls,
ctZero, cacheMode)
{
}
// Disallow copying and moving.
IntpBSplineUniform(IntpBSplineUniform const&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform const&) = delete;
IntpBSplineUniform(IntpBSplineUniform&&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform&&) = delete;
// Evaluate the interpolator. Each element of 'order' indicates the
// order of the derivative you want to compute. For the function
// value itself, pass in 'order' that has all 0 elements.
typename Controls::Type Evaluate(std::array<int, N> const& order,
std::array<Real, N> const& t)
{
if (this->mCacheMode == this->NO_CACHING)
{
return this->EvaluateNoCaching(order.data(), t.data());
}
else
{
return this->EvaluateCaching(order.data(), t.data());
}
}
};
}
// Implementation for B-spline interpolation whose dimension is known only
// at run time.
namespace WwiseGTE
{
template <typename Real, typename Controls>
class IntpBSplineUniform<Real, Controls, 0> : public IntpBSplineUniformShared<Real, Controls>
{
public:
// The caller is responsible for ensuring that the IntpBSplineUniform
// object persists as long as the input 'controls' exists.
IntpBSplineUniform(std::vector<int> const& degrees, Controls const& controls,
typename Controls::Type ctZero, int cacheMode)
:
IntpBSplineUniformShared<Real, Controls>(static_cast<int>(degrees.size()),
degrees.data(), controls, ctZero, cacheMode)
{
}
// Disallow copying and moving.
IntpBSplineUniform(IntpBSplineUniform const&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform const&) = delete;
IntpBSplineUniform(IntpBSplineUniform&&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform&&) = delete;
// Evaluate the interpolator. Each element of 'order' indicates the
// order of the derivative you want to compute. For the function
// value itself, pass in 'order' that has all 0 elements.
typename Controls::Type Evaluate(std::vector<int> const& order,
std::vector<Real> const& t)
{
if (static_cast<int>(order.size()) >= this->mNumDimensions
&& static_cast<int>(t.size()) >= this->mNumDimensions)
{
if (this->mCacheMode == this->NO_CACHING)
{
return this->EvaluateNoCaching(order.data(), t.data());
}
else
{
return this->EvaluateCaching(order.data(), t.data());
}
}
else
{
return this->mCTZero;
}
}
};
}
// To use only the N-dimensional template code above, define the symbol
// GTE_INTP_BSPLINE_UNIFORM_NO_SPECIALIZATION to disable the specializations
// that occur below. Notice that the interfaces are different between the
// specializations and the general code.
#if !defined(GTE_INTP_BSPLINE_UNIFORM_NO_SPECIALIZATION)
// Specialization for 1-dimensional data.
namespace WwiseGTE
{
template <typename Real, typename Controls>
class IntpBSplineUniform<Real, Controls, 1> : public IntpBSplineUniformShared<Real, Controls>
{
public:
// The caller is responsible for ensuring that the IntpBSplineUniform
// object persists as long as the input 'controls' exists.
IntpBSplineUniform(int degree, Controls const& controls,
typename Controls::Type ctZero, int cacheMode)
:
IntpBSplineUniformShared<Real, Controls>(),
mDegree(degree),
mControls(&controls),
mCTZero(ctZero),
mCacheMode(cacheMode)
{
// The condition c+1 > d+1 is required so that when s = c+1-d, its
// maximum value, we have at least two s-knots (d and d + 1).
if (mControls->GetSize(0) <= mDegree + 1)
{
LogError("Incompatible degree and number of controls.");
}
mDegreeP1 = mDegree + 1;
mNumControls = mControls->GetSize(0);
mTMin = (Real)-0.5;
mTMax = static_cast<Real>(mNumControls) - (Real)0.5;
mNumTRows = 0;
mNumTCols = 0;
this->ComputeBlendingMatrix(mDegree, mBlender);
this->ComputeDCoefficients(mDegree, mDCoefficient, mLMax);
this->ComputePowers(mDegree, mNumControls, mTMin, mTMax, mPowerDSDT);
if (mCacheMode == this->NO_CACHING)
{
mPhi.resize(mDegreeP1);
}
else
{
InitializeTensors();
}
}
// Disallow copying and moving.
IntpBSplineUniform(IntpBSplineUniform const&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform const&) = delete;
IntpBSplineUniform(IntpBSplineUniform&&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform&&) = delete;
// Member access.
inline int GetDegree(int) const
{
return mDegree;
}
inline int GetNumControls(int) const
{
return mNumControls;
}
inline Real GetTMin(int) const
{
return mTMin;
}
inline Real GetTMax(int) const
{
return mTMax;
}
inline int GetCacheMode() const
{
return mCacheMode;
}
// Evaluate the interpolator. The order is 0 when you want the B-spline
// function value itself. The order is 1 for the first derivative of the
// function, and so on.
typename Controls::Type Evaluate(std::array<int, 1> const& order,
std::array<Real,1> const& t)
{
auto result = mCTZero;
if (0 <= order[0] && order[0] <= mDegree)
{
int i;
Real u;
this->GetKey(t[0], mTMin, mTMax, mPowerDSDT[1], mNumControls, mDegree, i, u);
if (mCacheMode == this->NO_CACHING)
{
int jIndex = 0;
for (int j = 0; j <= mDegree; ++j)
{
int kjIndex = mDegree + jIndex;
int ell = mLMax[order[0]];
mPhi[j] = (Real)0;
for (int k = mDegree; k >= order[0]; --k)
{
mPhi[j] = mPhi[j] * u + mBlender[kjIndex--] * mDCoefficient[ell--];
}
jIndex += mDegreeP1;
}
for (int j = 0; j <= mDegree; ++j)
{
result = result + (*mControls)(i + j) * mPhi[j];
}
}
else
{
int iIndex = mNumTCols * i;
int kiIndex = mDegree + iIndex;
int ell = mLMax[order[0]];
for (int k = mDegree; k >= order[0]; --k)
{
if (mCacheMode == this->ON_DEMAND_CACHING && !mCached[kiIndex])
{
ComputeTensor(i, k, kiIndex);
mCached[kiIndex] = true;
}
result = result * u + mTensor[kiIndex--] * mDCoefficient[ell--];
}
}
result = result * mPowerDSDT[order[0]];
}
return result;
}
protected:
void ComputeTensor(int r, int c, int index)
{
auto element = mCTZero;
for (int j = 0; j <= mDegree; ++j)
{
element = element + (*mControls)(r + j) * mBlender[c + mDegreeP1 * j];
}
mTensor[index] = element;
}
void InitializeTensors()
{
mNumTRows = mNumControls - mDegree;
mNumTCols = mDegreeP1;
int numCached = mNumTRows * mNumTCols;
mTensor.resize(numCached);
mCached.resize(numCached);
if (mCacheMode == this->PRE_CACHING)
{
for (int r = 0, index = 0; r < mNumTRows; ++r)
{
for (int c = 0; c < mNumTCols; ++c, ++index)
{
ComputeTensor(r, c, index);
}
}
std::fill(mCached.begin(), mCached.end(), true);
}
else
{
std::fill(mCached.begin(), mCached.end(), false);
}
}
// Constructor inputs.
int mDegree;
Controls const* mControls;
typename Controls::Type mCTZero;
int mCacheMode;
// Parameters for B-spline evaluation.
int mDegreeP1;
int mNumControls;
Real mTMin, mTMax;
std::vector<Real> mBlender;
std::vector<Real> mDCoefficient;
std::vector<int> mLMax;
std::vector<Real> mPowerDSDT;
// Support for no-cached B-spline evaluation.
std::vector<Real> mPhi;
// Support for cached B-spline evaluation.
int mNumTRows, mNumTCols;
std::vector<typename Controls::Type> mTensor;
std::vector<bool> mCached;
};
}
// Specialization for 2-dimensional data.
namespace WwiseGTE
{
template <typename Real, typename Controls>
class IntpBSplineUniform<Real, Controls, 2> : public IntpBSplineUniformShared<Real, Controls>
{
public:
// The caller is responsible for ensuring that the IntpBSplineUniform2
// object persists as long as the input 'controls' exists.
IntpBSplineUniform(std::array<int, 2> const& degrees, Controls const& controls,
typename Controls::Type ctZero, int cacheMode)
:
IntpBSplineUniformShared<Real, Controls>(),
mDegree(degrees),
mControls(&controls),
mCTZero(ctZero),
mCacheMode(cacheMode)
{
// The condition c+1 > d+1 is required so that when s = c+1-d, its
// maximum value, we have at least two s-knots (d and d + 1).
for (int dim = 0; dim < 2; ++dim)
{
if (mControls->GetSize(dim) <= mDegree[dim] + 1)
{
LogError("Incompatible degree and number of controls.");
}
}
for (int dim = 0; dim < 2; ++dim)
{
mDegreeP1[dim] = mDegree[dim] + 1;
mNumControls[dim] = mControls->GetSize(dim);
mTMin[dim] = (Real)-0.5;
mTMax[dim] = static_cast<Real>(mNumControls[dim]) - (Real)0.5;
mNumTRows[dim] = 0;
mNumTCols[dim] = 0;
this->ComputeBlendingMatrix(mDegree[dim], mBlender[dim]);
this->ComputeDCoefficients(mDegree[dim], mDCoefficient[dim], mLMax[dim]);
this->ComputePowers(mDegree[dim], mNumControls[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim]);
}
if (mCacheMode == this->NO_CACHING)
{
for (int dim = 0; dim < 2; ++dim)
{
mPhi[dim].resize(mDegreeP1[dim]);
}
}
else
{
InitializeTensors();
}
}
// Disallow copying and moving.
IntpBSplineUniform(IntpBSplineUniform const&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform const&) = delete;
IntpBSplineUniform(IntpBSplineUniform&&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform&&) = delete;
// Member access.
inline int GetDegree(int dim) const
{
return mDegree[dim];
}
inline int GetNumControls(int dim) const
{
return mNumControls[dim];
}
inline Real GetTMin(int dim) const
{
return mTMin[dim];
}
inline Real GetTMax(int dim) const
{
return mTMax[dim];
}
inline int GetCacheMode() const
{
return mCacheMode;
}
// Evaluate the interpolator. The order is (0,0) when you want the
// B-spline function value itself. The order0 is 1 for the first
// derivative with respect to t0 and the order1 is 1 for the first
// derivative with respect to t1. Higher-order derivatives in other
// t-inputs are computed similarly.
typename Controls::Type Evaluate(std::array<int, 2> const& order,
std::array<Real, 2> const& t)
{
auto result = mCTZero;
if (0 <= order[0] && order[0] <= mDegree[0]
&& 0 <= order[1] && order[1] <= mDegree[1])
{
std::array<int, 2> i;
std::array<Real, 2> u;
for (int dim = 0; dim < 2; ++dim)
{
this->GetKey(t[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim][1],
mNumControls[dim], mDegree[dim], i[dim], u[dim]);
}
if (mCacheMode == this->NO_CACHING)
{
for (int dim = 0; dim < 2; ++dim)
{
int jIndex = 0;
for (int j = 0; j <= mDegree[dim]; ++j)
{
int kjIndex = mDegree[dim] + jIndex;
int ell = mLMax[dim][order[dim]];
mPhi[dim][j] = (Real)0;
for (int k = mDegree[dim]; k >= order[dim]; --k)
{
mPhi[dim][j] = mPhi[dim][j] * u[dim] +
mBlender[dim][kjIndex--] * mDCoefficient[dim][ell--];
}
jIndex += mDegreeP1[dim];
}
}
for (int j1 = 0; j1 <= mDegree[1]; ++j1)
{
Real phi1 = mPhi[1][j1];
for (int j0 = 0; j0 <= mDegree[0]; ++j0)
{
Real phi0 = mPhi[0][j0];
Real phi01 = phi0 * phi1;
result = result + (*mControls)(i[0] + j0, i[1] + j1) * phi01;
}
}
}
else
{
int i0i1Index = mNumTCols[1] * (i[0] + mNumTRows[0] * i[1]);
int k1i0i1Index = mDegree[1] + i0i1Index;
int ell1 = mLMax[1][order[1]];
for (int k1 = mDegree[1]; k1 >= order[1]; --k1)
{
int k0k1i0i1Index = mDegree[0] + mNumTCols[0] * k1i0i1Index;
int ell0 = mLMax[0][order[0]];
auto term = mCTZero;
for (int k0 = mDegree[0]; k0 >= order[0]; --k0)
{
if (mCacheMode == this->ON_DEMAND_CACHING && !mCached[k0k1i0i1Index])
{
ComputeTensor(i[0], i[1], k0, k1, k0k1i0i1Index);
mCached[k0k1i0i1Index] = true;
}
term = term * u[0] + mTensor[k0k1i0i1Index--] * mDCoefficient[0][ell0--];
}
result = result * u[1] + term * mDCoefficient[1][ell1--];
--k1i0i1Index;
}
}
Real adjust(1);
for (int dim = 0; dim < 2; ++dim)
{
adjust *= mPowerDSDT[dim][order[dim]];
}
result = result * adjust;
}
return result;
}
void ComputeTensor(int r0, int r1, int c0, int c1, int index)
{
auto element = mCTZero;
for (int j1 = 0; j1 <= mDegree[1]; ++j1)
{
Real blend1 = mBlender[1][c1 + mDegreeP1[1] * j1];
for (int j0 = 0; j0 <= mDegree[0]; ++j0)
{
Real blend0 = mBlender[0][c0 + mDegreeP1[0] * j0];
Real blend01 = blend0 * blend1;
element = element + (*mControls)(r0 + j0, r1 + j1) * blend01;
}
}
mTensor[index] = element;
}
void InitializeTensors()
{
int numCached = 1;
for (int dim = 0; dim < 2; ++dim)
{
mNumTRows[dim] = mNumControls[dim] - mDegree[dim];
mNumTCols[dim] = mDegreeP1[dim];
numCached *= mNumTRows[dim] * mNumTCols[dim];
}
mTensor.resize(numCached);
mCached.resize(numCached);
if (mCacheMode == this->PRE_CACHING)
{
for (int r1 = 0, index = 0; r1 < mNumTRows[1]; ++r1)
{
for (int r0 = 0; r0 < mNumTRows[0]; ++r0)
{
for (int c1 = 0; c1 < mNumTCols[1]; ++c1)
{
for (int c0 = 0; c0 < mNumTCols[0]; ++c0, ++index)
{
ComputeTensor(r0, r1, c0, c1, index);
}
}
}
}
std::fill(mCached.begin(), mCached.end(), true);
}
else
{
std::fill(mCached.begin(), mCached.end(), false);
}
}
// Constructor inputs.
std::array<int, 2> mDegree;
Controls const* mControls;
typename Controls::Type mCTZero;
int mCacheMode;
// Parameters for B-spline evaluation.
std::array<int, 2> mDegreeP1;
std::array<int, 2> mNumControls;
std::array<Real, 2> mTMin, mTMax;
std::array<std::vector<Real>, 2> mBlender;
std::array<std::vector<Real>, 2> mDCoefficient;
std::array<std::vector<int>, 2> mLMax;
std::array<std::vector<Real>, 2> mPowerDSDT;
// Support for no-cached B-spline evaluation.
std::array<std::vector<Real>, 2> mPhi;
// Support for cached B-spline evaluation.
std::array<int, 2> mNumTRows, mNumTCols;
std::vector<typename Controls::Type> mTensor;
std::vector<bool> mCached;
};
}
// Specialization for 3-dimensional data.
namespace WwiseGTE
{
template <typename Real, typename Controls>
class IntpBSplineUniform<Real, Controls, 3> : public IntpBSplineUniformShared<Real, Controls>
{
public:
// The caller is responsible for ensuring that the IntpBSplineUniform3
// object persists as long as the input 'controls' exists.
IntpBSplineUniform(std::array<int, 3> const& degrees, Controls const& controls,
typename Controls::Type ctZero, int cacheMode)
:
IntpBSplineUniformShared<Real, Controls>(),
mDegree(degrees),
mControls(&controls),
mCTZero(ctZero),
mCacheMode(cacheMode)
{
// The condition c+1 > d+1 is required so that when s = c+1-d, its
// maximum value, we have at least two s-knots (d and d + 1).
for (int dim = 0; dim < 3; ++dim)
{
if (mControls->GetSize(dim) <= mDegree[dim] + 1)
{
LogError("Incompatible degree and number of controls.");
}
}
for (int dim = 0; dim < 3; ++dim)
{
mDegreeP1[dim] = mDegree[dim] + 1;
mNumControls[dim] = mControls->GetSize(dim);
mTMin[dim] = (Real)-0.5;
mTMax[dim] = static_cast<Real>(mNumControls[dim]) - (Real)0.5;
mNumTRows[dim] = 0;
mNumTCols[dim] = 0;
this->ComputeBlendingMatrix(mDegree[dim], mBlender[dim]);
this->ComputeDCoefficients(mDegree[dim], mDCoefficient[dim], mLMax[dim]);
this->ComputePowers(mDegree[dim], mNumControls[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim]);
}
if (mCacheMode == this->NO_CACHING)
{
for (int dim = 0; dim < 3; ++dim)
{
mPhi[dim].resize(mDegreeP1[dim]);
}
}
else
{
InitializeTensors();
}
}
// Disallow copying and moving.
IntpBSplineUniform(IntpBSplineUniform const&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform const&) = delete;
IntpBSplineUniform(IntpBSplineUniform&&) = delete;
IntpBSplineUniform& operator=(IntpBSplineUniform&&) = delete;
// Member access. The input i specifies the dimension (0, 1, 2).
inline int GetDegree(int dim) const
{
return mDegree[dim];
}
inline int GetNumControls(int dim) const
{
return mNumControls[dim];
}
inline Real GetTMin(int dim) const
{
return mTMin[dim];
}
inline Real GetTMax(int dim) const
{
return mTMax[dim];
}
// Evaluate the interpolator. The order is (0,0,0) when you want the
// B-spline function value itself. The order0 is 1 for the first
// derivative with respect to t0, the order1 is 1 for the first
// derivative with respect to t1 or the order2 is 1 for the first
// derivative with respect to t2. Higher-order derivatives in other
// t-inputs are computed similarly.
typename Controls::Type Evaluate(std::array<int, 3> const& order,
std::array<Real, 3> const& t)
{
auto result = mCTZero;
if (0 <= order[0] && order[0] <= mDegree[0]
&& 0 <= order[1] && order[1] <= mDegree[1]
&& 0 <= order[2] && order[2] <= mDegree[2])
{
std::array<int, 3> i;
std::array<Real, 3> u;
for (int dim = 0; dim < 3; ++dim)
{
this->GetKey(t[dim], mTMin[dim], mTMax[dim], mPowerDSDT[dim][1],
mNumControls[dim], mDegree[dim], i[dim], u[dim]);
}
if (mCacheMode == this->NO_CACHING)
{
for (int dim = 0; dim < 3; ++dim)
{
int jIndex = 0;
for (int j = 0; j <= mDegree[dim]; ++j)
{
int kjIndex = mDegree[dim] + jIndex;
int ell = mLMax[dim][order[dim]];
mPhi[dim][j] = (Real)0;
for (int k = mDegree[dim]; k >= order[dim]; --k)
{
mPhi[dim][j] = mPhi[dim][j] * u[dim] +
mBlender[dim][kjIndex--] * mDCoefficient[dim][ell--];
}
jIndex += mDegreeP1[dim];
}
}
for (int j2 = 0; j2 <= mDegree[2]; ++j2)
{
Real phi2 = mPhi[2][j2];
for (int j1 = 0; j1 <= mDegree[1]; ++j1)
{
Real phi1 = mPhi[1][j1];
Real phi12 = phi1 * phi2;
for (int j0 = 0; j0 <= mDegree[0]; ++j0)
{
Real phi0 = mPhi[0][j0];
Real phi012 = phi0 * phi12;
result = result + (*mControls)(i[0] + j0, i[1] + j1, i[2] + j2) * phi012;
}
}
}
}
else
{
int i0i1i2Index = mNumTCols[2] * (i[0] + mNumTRows[0] * (i[1] + mNumTRows[1] * i[2]));
int k2i0i1i2Index = mDegree[2] + i0i1i2Index;
int ell2 = mLMax[2][order[2]];
for (int k2 = mDegree[2]; k2 >= order[2]; --k2)
{
int k1k2i0i1i2Index = mDegree[1] + mNumTCols[1] * k2i0i1i2Index;
int ell1 = mLMax[1][order[1]];
auto term1 = mCTZero;
for (int k1 = mDegree[1]; k1 >= order[1]; --k1)
{
int k0k1k2i0i1i2Index = mDegree[0] + mNumTCols[0] * k1k2i0i1i2Index;
int ell0 = mLMax[0][order[0]];
auto term0 = mCTZero;
for (int k0 = mDegree[0]; k0 >= order[0]; --k0)
{
if (mCacheMode == this->ON_DEMAND_CACHING && !mCached[k0k1k2i0i1i2Index])
{
ComputeTensor(i[0], i[1], i[2], k0, k1, k2, k0k1k2i0i1i2Index);
mCached[k0k1k2i0i1i2Index] = true;
}
term0 = term0 * u[0] + mTensor[k0k1k2i0i1i2Index--] * mDCoefficient[0][ell0--];
}
term1 = term1 * u[1] + term0 * mDCoefficient[1][ell1--];
--k1k2i0i1i2Index;
}
result = result * u[2] + term1 * mDCoefficient[2][ell2--];
--k2i0i1i2Index;
}
}
Real adjust(1);
for (int dim = 0; dim < 3; ++dim)
{
adjust *= mPowerDSDT[dim][order[dim]];
}
result = result * adjust;
}
return result;
}
protected:
void ComputeTensor(int r0, int r1, int r2, int c0, int c1, int c2, int index)
{
auto element = mCTZero;
for (int j2 = 0; j2 <= mDegree[2]; ++j2)
{
Real blend2 = mBlender[2][c2 + mDegreeP1[2] * j2];
for (int j1 = 0; j1 <= mDegree[1]; ++j1)
{
Real blend1 = mBlender[1][c1 + mDegreeP1[1] * j1];
Real blend12 = blend1 * blend2;
for (int j0 = 0; j0 <= mDegree[0]; ++j0)
{
Real blend0 = mBlender[0][c0 + mDegreeP1[0] * j0];
Real blend012 = blend0 * blend12;
element = element + (*mControls)(r0 + j0, r1 + j1, r2 + j2) * blend012;
}
}
}
mTensor[index] = element;
}
void InitializeTensors()
{
int numCached = 1;
for (int dim = 0; dim < 3; ++dim)
{
mNumTRows[dim] = mNumControls[dim] - mDegree[dim];
mNumTCols[dim] = mDegreeP1[dim];
numCached *= mNumTRows[dim] * mNumTCols[dim];
}
mTensor.resize(numCached);
mCached.resize(numCached);
if (mCacheMode == this->PRE_CACHING)
{
for (int r2 = 0, index = 0; r2 < mNumTRows[2]; ++r2)
{
for (int r1 = 0; r1 < mNumTRows[1]; ++r1)
{
for (int r0 = 0; r0 < mNumTRows[0]; ++r0)
{
for (int c2 = 0; c2 < mNumTCols[2]; ++c2)
{
for (int c1 = 0; c1 < mNumTCols[1]; ++c1)
{
for (int c0 = 0; c0 < mNumTCols[0]; ++c0, ++index)
{
ComputeTensor(r0, r1, r2, c0, c1, c2, index);
}
}
}
}
}
}
std::fill(mCached.begin(), mCached.end(), true);
}
else
{
std::fill(mCached.begin(), mCached.end(), false);
}
}
// Constructor inputs.
std::array<int, 3> mDegree;
Controls const* mControls;
typename Controls::Type mCTZero;
int mCacheMode;
// Parameters for B-spline evaluation.
std::array<int, 3> mDegreeP1;
std::array<int, 3> mNumControls;
std::array<Real, 3> mTMin, mTMax;
std::array<std::vector<Real>, 3> mBlender;
std::array<std::vector<Real>, 3> mDCoefficient;
std::array<std::vector<int>, 3> mLMax;
std::array<std::vector<Real>, 3> mPowerDSDT;
// Support for no-cached B-spline evaluation.
std::array<std::vector<Real>, 3> mPhi;
// Support for cached B-spline evaluation.
std::array<int, 3> mNumTRows, mNumTCols;
std::vector<typename Controls::Type> mTensor;
std::vector<bool> mCached;
};
}
#endif