/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageBSplineInterpolator.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm 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 "vtkImageBSplineInterpolator.h"
#include "vtkDataArray.h"
#include "vtkImageBSplineInternals.h"
#include "vtkImageData.h"
#include "vtkImageInterpolatorInternals.h"
#include "vtkObjectFactory.h"

#include "vtkTemplateAliasMacro.h"
// turn off 64-bit ints when templating over all types, because
// they cannot be faithfully represented by doubles
#undef VTK_USE_INT64
#define VTK_USE_INT64 0
#undef VTK_USE_UINT64
#define VTK_USE_UINT64 0

#define VTK_BSPLINE_KERNEL_SIZE_MAX (VTK_IMAGE_BSPLINE_DEGREE_MAX + 1)
#define VTK_BSPLINE_INT_INITIALIZER                                                                \
  {                                                                                                \
    0, 0, 0, 0, 0, 0, 0, 0, 0, 0,                                                                  \
  }

// kernel lookup table size must be 256*n where n is kernel half-width
// in order to provide sufficient precision for 16-bit images
#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
#define VTK_BSPLINE_KERNEL_TABLE_DIVISIONS 256
#endif

vtkStandardNewMacro(vtkImageBSplineInterpolator);

//------------------------------------------------------------------------------
vtkImageBSplineInterpolator::vtkImageBSplineInterpolator()
{
  this->SplineDegree = 3;
  this->KernelLookupTable = nullptr;
}

//------------------------------------------------------------------------------
vtkImageBSplineInterpolator::~vtkImageBSplineInterpolator()
{
  if (this->KernelLookupTable)
  {
    this->FreeKernelLookupTable();
  }
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "SplineDegree: " << this->SplineDegree << "\n";
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::ComputeSupportSize(
  const double vtkNotUsed(matrix)[16], int size[3])
{
  int n = this->SplineDegree + 1;
  size[0] = n;
  size[1] = n;
  size[2] = n;
}

//------------------------------------------------------------------------------
bool vtkImageBSplineInterpolator::IsSeparable()
{
  return true;
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::SetSplineDegree(int degree)
{
  degree = vtkMath::ClampValue(degree, 0, VTK_IMAGE_BSPLINE_DEGREE_MAX);
  if (this->SplineDegree != degree)
  {
    this->SplineDegree = degree;
    this->Modified();
  }
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::InternalDeepCopy(vtkAbstractImageInterpolator* a)
{
  vtkImageBSplineInterpolator* obj = vtkImageBSplineInterpolator::SafeDownCast(a);

  if (obj)
  {
    this->SetSplineDegree(obj->SplineDegree);
  }

  if (this->KernelLookupTable)
  {
    this->FreeKernelLookupTable();
  }
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::InternalUpdate()
{
  int mode = this->SplineDegree;

  if (this->InterpolationInfo->InterpolationMode != mode || this->KernelLookupTable == nullptr)
  {
    this->BuildKernelLookupTable();
  }

  this->InterpolationInfo->InterpolationMode = mode;
  this->InterpolationInfo->ExtraInfo = this->KernelLookupTable;
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//  Interpolation subroutines and associated code
//------------------------------------------------------------------------------

namespace
{

#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
//------------------------------------------------------------------------------
// B-spline kernel computation: compute half of the symmetric kernel.
// In the table, x=0.0 corresponds to index position zero, and each
// "m" bins corresponds to a unit spacing.  The full size of the table
// will be m*(n + 1)/2.

struct vtkBSplineKernel
{
  template <class F>
  static void BSpline(F* kernel, int size, int n);
};

template <class F>
void vtkBSplineKernel::BSpline(F* kernel, int m, int n)
{
  long order = n;
  int kn = (n + 2) / 2;
  int km = n - kn + 1;
  double offset = 0;
  double delta = 1.0 / m;
  double weights[VTK_BSPLINE_KERNEL_SIZE_MAX];

  // special case for order zero
  if (n == 0)
  {
    int k = m * (n + 1) / 2;
    do
    {
      *kernel++ = 1.0;
    } while (--k);
    return;
  }

  // for offset of zero, weights are symmetrical
  vtkImageBSplineInternals::GetInterpolationWeights(weights, offset, order);
  double* weights2 = weights + kn;
  F* kernel2 = kernel;
  int k2 = kn;
  do
  {
    --weights2;
    *kernel2 = *weights2;
    kernel2 += m;
  } while (--k2);

  // need the opposite end of the kernel array
  kernel2 = kernel + km * m;

  int j = m / 2;
  if (j)
    do
    {
      kernel++;
      kernel2--;
      offset += delta;

      vtkImageBSplineInternals::GetInterpolationWeights(weights, offset, order);

      double* weights1 = weights + n + 1;
      F* kernel1 = kernel2;
      int k = km;
      do
      {
        --weights1;
        *kernel1 = *weights1;
        kernel1 -= m;
      } while (--k);
      kernel1 = kernel;
      k = kn;
      do
      {
        --weights1;
        *kernel1 = *weights1;
        kernel1 += m;
      } while (--k);
    } while (--j);
}

//------------------------------------------------------------------------------
// Use kernel lookup table, might improve performance for large kernels
template <class T, class F>
void vtkBSplineInterpWeights(T* kernel, F* fX, F fx, int m)
{
  if (m == 0)
  {
    fX[0] = 1;
    return;
  }

  // table bins per unit
  int p = VTK_BSPLINE_KERNEL_TABLE_DIVISIONS;

  // compute table interpolation info
  F f = fx * p;
  int offset = static_cast<int>(f);
  f -= offset;
  F r = 1 - f;

  // interpolate the table
  int n = m + 1;
  F s = 0;
  int i = (1 - ((n + 1) >> 1)) * p - offset;
  do
  {
    int i0 = i;
    int i1 = i + 1;
    int ni = -i0;
    i0 = ((i0 >= 0) ? i0 : ni);
    ni = -i1;
    i1 = ((i1 >= 0) ? i1 : ni);
    F y = r * kernel[i0] + f * kernel[i1];
    *fX = y;
    s += y;
    i += p;
    fX++;
  } while (--n);
}
#endif

//------------------------------------------------------------------------------
template <class F, class T>
struct vtkImageBSplineInterpolate
{
  static void BSpline(vtkInterpolationInfo* info, const F point[3], F* outPtr);
};

//------------------------------------------------------------------------------
template <class F, class T>
void vtkImageBSplineInterpolate<F, T>::BSpline(
  vtkInterpolationInfo* info, const F point[3], F* outPtr)
{
  const T* inPtr = static_cast<const T*>(info->Pointer);
  int* inExt = info->Extent;
  vtkIdType* inInc = info->Increments;
  int numscalars = info->NumberOfComponents;

#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
  // kernel lookup table
  float* kernel = static_cast<float*>(info->ExtraInfo);
#endif

  // size of kernel is degree of spline plus one
  int n = info->InterpolationMode;
  int m = n + 1;

  // index to kernel midpoint position
  int m2 = (n >> 1);

  // offset for odd-size kernels
  F offset = 0.5 * (m & 1);

  F fx, fy, fz;
  int inIdX0 = vtkInterpolationMath::Floor(point[0] + offset, fx);
  int inIdY0 = vtkInterpolationMath::Floor(point[1] + offset, fy);
  int inIdZ0 = vtkInterpolationMath::Floor(point[2] + offset, fz);

  fx -= offset;
  fy -= offset;
  fz -= offset;

  // change arrays into locals
  vtkIdType inIncX = inInc[0];
  vtkIdType inIncY = inInc[1];
  vtkIdType inIncZ = inInc[2];

  int minX = inExt[0];
  int maxX = inExt[1];
  int minY = inExt[2];
  int maxY = inExt[3];
  int minZ = inExt[4];
  int maxZ = inExt[5];

  // the memory offsets
  vtkIdType factX[VTK_BSPLINE_KERNEL_SIZE_MAX + 4];
  vtkIdType factY[VTK_BSPLINE_KERNEL_SIZE_MAX];
  vtkIdType factZ[VTK_BSPLINE_KERNEL_SIZE_MAX];

  // handle the borders
  int xi = inIdX0 - m2;
  int yi = inIdY0 - m2;
  int zi = inIdZ0 - m2;

  switch (info->BorderMode)
  {
    case VTK_IMAGE_BORDER_REPEAT:
    {
      int l = 0;
      int mm = m;
      do
      {
        factX[l] = vtkInterpolationMath::Wrap(xi, minX, maxX) * inIncX;
        factY[l] = vtkInterpolationMath::Wrap(yi, minY, maxY) * inIncY;
        factZ[l] = vtkInterpolationMath::Wrap(zi, minZ, maxZ) * inIncZ;
        l++;
        xi++;
        yi++;
        zi++;
      } while (--mm);
    }
    break;

    case VTK_IMAGE_BORDER_MIRROR:
    {
      int l = 0;
      int mm = m;
      do
      {
        factX[l] = vtkInterpolationMath::Mirror(xi, minX, maxX) * inIncX;
        factY[l] = vtkInterpolationMath::Mirror(yi, minY, maxY) * inIncY;
        factZ[l] = vtkInterpolationMath::Mirror(zi, minZ, maxZ) * inIncZ;
        l++;
        xi++;
        yi++;
        zi++;
      } while (--mm);
    }
    break;

    default:
    {
      int l = 0;
      int mm = m;
      do
      {
        factX[l] = vtkInterpolationMath::Clamp(xi, minX, maxX) * inIncX;
        factY[l] = vtkInterpolationMath::Clamp(yi, minY, maxY) * inIncY;
        factZ[l] = vtkInterpolationMath::Clamp(zi, minZ, maxZ) * inIncZ;
        l++;
        xi++;
        yi++;
        zi++;
      } while (--mm);
    }
    break;
  }

  // compute the kernel weights (pad X for loop unrolling)
  F fX[VTK_BSPLINE_KERNEL_SIZE_MAX + 4];
  F fY[VTK_BSPLINE_KERNEL_SIZE_MAX];
  F fZ[VTK_BSPLINE_KERNEL_SIZE_MAX];

  // check if only one slice in a particular direction
  int nx = n * (minX != maxX);
  int ny = n * (minY != maxY);
  int nz = n * (minZ != maxZ);

#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
  vtkBSplineInterpWeights(kernel, fX, fx, nx);
  vtkBSplineInterpWeights(kernel, fY, fy, ny);
  vtkBSplineInterpWeights(kernel, fZ, fz, nz);
#else
  vtkImageBSplineInternals::GetInterpolationWeights(fX, fx, nx);
  vtkImageBSplineInternals::GetInterpolationWeights(fY, fy, ny);
  vtkImageBSplineInternals::GetInterpolationWeights(fZ, fz, nz);
#endif

  // pad coeffs to allow unrolling of inner loop
  int lm = ((nx + 4) >> 2);
  vtkIdType factXl = factX[nx];
  vtkIdType* tmpfactXl = &factX[nx + 1];
  F* tmpfXl = &fX[nx + 1];
  tmpfactXl[0] = factXl;
  tmpfactXl[1] = factXl;
  tmpfactXl[2] = factXl;
  tmpfXl[0] = 0;
  tmpfXl[1] = 0;
  tmpfXl[2] = 0;

  do // loop over components
  {
    F val = 0;
    int k = 0;
    do // loop over z
    {
      F ifz = fZ[k];
      vtkIdType factz = factZ[k];
      int j = 0;
      do // loop over y
      {
        F ify = fY[j];
        F fzy = ifz * ify;
        vtkIdType factzy = factz + factY[j];
        // loop over x
        const T* tmpPtr = inPtr + factzy;
        const F* tmpfX = fX;
        const vtkIdType* tmpfactX = factX;
        F tmpval = 0;
        // unroll inner loop for efficiency
        int l = lm;
        do
        {
          tmpval += tmpfX[0] * tmpPtr[tmpfactX[0]];
          tmpval += tmpfX[1] * tmpPtr[tmpfactX[1]];
          tmpval += tmpfX[2] * tmpPtr[tmpfactX[2]];
          tmpval += tmpfX[3] * tmpPtr[tmpfactX[3]];
          tmpfX += 4;
          tmpfactX += 4;
        } while (--l);
        val += fzy * tmpval;
      } while (++j <= ny);
    } while (++k <= nz);

    *outPtr++ = val;
    inPtr++;
  } while (--numscalars);
}

//------------------------------------------------------------------------------
// Get the interpolation function for the specified data types
template <class F>
void vtkImageBSplineInterpolatorGetInterpolationFunc(
  void (**interpolate)(vtkInterpolationInfo*, const F[3], F*), int dataType,
  int vtkNotUsed(interpolationMode))
{
  switch (dataType)
  {
    vtkTemplateAliasMacro(*interpolate = &(vtkImageBSplineInterpolate<F, VTK_TT>::BSpline));
    default:
      *interpolate = nullptr;
  }
}

//------------------------------------------------------------------------------
// Interpolation for precomputed weights

template <class F, class T>
struct vtkImageBSplineRowInterpolate
{
  static void BSpline(
    vtkInterpolationWeights* weights, int idX, int idY, int idZ, F* outPtr, int n);
};

//------------------------------------------------------------------------------
// helper function for high-order interpolation
template <class F, class T>
void vtkImageBSplineRowInterpolate<F, T>::BSpline(
  vtkInterpolationWeights* weights, int idX, int idY, int idZ, F* outPtr, int n)
{
  int stepX = weights->KernelSize[0];
  int stepY = weights->KernelSize[1];
  int stepZ = weights->KernelSize[2];
  idX *= stepX;
  idY *= stepY;
  idZ *= stepZ;
  const F* fX = static_cast<F*>(weights->Weights[0]) + idX;
  const F* fY = static_cast<F*>(weights->Weights[1]) + idY;
  const F* fZ = static_cast<F*>(weights->Weights[2]) + idZ;
  const vtkIdType* factX = weights->Positions[0] + idX;
  const vtkIdType* factY = weights->Positions[1] + idY;
  const vtkIdType* factZ = weights->Positions[2] + idZ;
  const T* inPtr = static_cast<const T*>(weights->Pointer);

  // part of the loop unrolling scheme
  int lm = ((stepX + 3) >> 2);

  int numscalars = weights->NumberOfComponents;
  for (int i = n; i > 0; --i)
  {
    // allow unrolling of inner loop
    F fX1[VTK_BSPLINE_KERNEL_SIZE_MAX + 4];
    vtkIdType factX1[VTK_BSPLINE_KERNEL_SIZE_MAX + 4];
    F* tmpfX = fX1;
    vtkIdType* tmpfactX = factX1;
    int ii = stepX;
    do
    {
      *tmpfX++ = *fX++;
      *tmpfactX++ = *factX++;
    } while (--ii);
    vtkIdType lfactX = *(tmpfactX - 1);
    tmpfX[0] = 0.0;
    tmpfactX[0] = lfactX;
    tmpfX[1] = 0.0;
    tmpfactX[1] = lfactX;
    tmpfX[2] = 0.0;
    tmpfactX[2] = lfactX;

    const T* inPtr0 = inPtr;
    int c = numscalars;
    do // loop over components
    {
      F val = 0;
      int k = 0;
      do // loop over z
      {
        F ifz = fZ[k];
        vtkIdType factz = factZ[k];
        int j = 0;
        do // loop over y
        {
          F ify = fY[j];
          F fzy = ifz * ify;
          vtkIdType factzy = factz + factY[j];
          // loop over x
          const T* tmpPtr = inPtr0 + factzy;
          tmpfX = fX1;
          tmpfactX = factX1;
          F tmpval = 0;
          // unroll inner loop for efficiency
          int l = lm;
          do
          {
            tmpval += tmpfX[0] * tmpPtr[tmpfactX[0]];
            tmpval += tmpfX[1] * tmpPtr[tmpfactX[1]];
            tmpval += tmpfX[2] * tmpPtr[tmpfactX[2]];
            tmpval += tmpfX[3] * tmpPtr[tmpfactX[3]];
            tmpfX += 4;
            tmpfactX += 4;
          } while (--l);
          val += fzy * tmpval;
        } while (++j < stepY);
      } while (++k < stepZ);

      *outPtr++ = val;
      inPtr0++;
    } while (--c);
  }
}

//------------------------------------------------------------------------------
// get row interpolation function for different interpolation modes
// and different scalar types
template <class F>
void vtkImageBSplineInterpolatorGetRowInterpolationFunc(
  void (**summation)(vtkInterpolationWeights* weights, int idX, int idY, int idZ, F* outPtr, int n),
  int scalarType, int vtkNotUsed(interpolationMode))
{
  switch (scalarType)
  {
    vtkTemplateAliasMacro(*summation = &(vtkImageBSplineRowInterpolate<F, VTK_TT>::BSpline));
    default:
      *summation = nullptr;
  }
}

//------------------------------------------------------------------------------
template <class F>
void vtkImageBSplineInterpolatorPrecomputeWeights(const F newmat[16], const int outExt[6],
  int clipExt[6], const F bounds[6], vtkInterpolationWeights* weights)
{
#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
  float* kernel = static_cast<float*>(weights->ExtraInfo);
#endif
  weights->WeightType = vtkTypeTraits<F>::VTKTypeID();
  int degree = weights->InterpolationMode;
  int m = degree + 1;

  // set up input positions table for interpolation
  bool validClip = true;
  for (int j = 0; j < 3; j++)
  {
    // set k to the row for which the element in column j is nonzero,
    // and set matrow to the elements of that row
    int k = 0;
    const F* matrow = newmat;
    while (k < 3 && matrow[j] == 0)
    {
      k++;
      matrow += 4;
    }

    // get the extents
    clipExt[2 * j] = outExt[2 * j];
    clipExt[2 * j + 1] = outExt[2 * j + 1];
    int minExt = weights->Extent[2 * k];
    int maxExt = weights->Extent[2 * k + 1];
    F minBounds = bounds[2 * k];
    F maxBounds = bounds[2 * k + 1];

    // the kernel size should not exceed the input dimension
    int m2 = ((m - 1) >> 1);
    int step = m;
    int inCount = maxExt - minExt + 1;
    step = ((step < inCount) ? step : inCount);

    // allocate space for the weights
    vtkIdType size = step * (outExt[2 * j + 1] - outExt[2 * j] + 1);
    vtkIdType* positions = new vtkIdType[size];
    positions -= step * outExt[2 * j];
    F* constants = new F[size];
    constants -= step * outExt[2 * j];

    weights->KernelSize[j] = step;
    weights->Positions[j] = positions;
    weights->Weights[j] = constants;
    weights->WeightExtent[2 * j] = outExt[2 * j];
    weights->WeightExtent[2 * j + 1] = outExt[2 * j + 1];

    int region = 0;
    for (int i = outExt[2 * j]; i <= outExt[2 * j + 1]; i++)
    {
      F point = matrow[3] + i * matrow[j];

      // offset for odd-size kernels
      F offset = 0.5 * (m & 1);
      F f = 0;
      int idx = vtkInterpolationMath::Floor(point + offset, f);
      f -= offset;
      if (step > 1)
      {
        idx -= m2;
      }

      // compute the weights and offsets
      vtkIdType inInc = weights->Increments[k];
      if (inCount == 1)
      {
        positions[step * i] = 0;
        constants[step * i] = static_cast<F>(1);
      }
      else
      {
        // initialization is needed to avoid a warning for gcc 4.9.2, but
        // not for other compilers or for valgrind
        int inId[VTK_BSPLINE_KERNEL_SIZE_MAX] = VTK_BSPLINE_INT_INITIALIZER;

        int l = 0;
        switch (weights->BorderMode)
        {
          case VTK_IMAGE_BORDER_REPEAT:
            do
            {
              inId[l] = vtkInterpolationMath::Wrap(idx++, minExt, maxExt);
            } while (++l < m);
            break;

          case VTK_IMAGE_BORDER_MIRROR:
            do
            {
              inId[l] = vtkInterpolationMath::Mirror(idx++, minExt, maxExt);
            } while (++l < m);
            break;

          default:
            do
            {
              inId[l] = vtkInterpolationMath::Clamp(idx++, minExt, maxExt);
            } while (++l < m);
            break;
        }

        F g[VTK_BSPLINE_KERNEL_SIZE_MAX];
#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
        vtkBSplineInterpWeights(kernel, g, f, m - 1);
#else
        vtkImageBSplineInternals::GetInterpolationWeights(g, f, m - 1);
#endif
        if (step == m)
        {
          int ll = 0;
          do
          {
            positions[step * i + ll] = inId[ll] * inInc;
            constants[step * i + ll] = g[ll];
          } while (++ll < step);
        }
        else
        {
          // it gets tricky if the data is thinner than the kernel
          F gg[VTK_BSPLINE_KERNEL_SIZE_MAX];
          int ll = 0;
          do
          {
            gg[ll] = 0;
          } while (++ll < m);
          ll = 0;
          do
          {
            int rIdx = inId[ll];
            gg[rIdx] += g[ll];
          } while (++ll < m);
          ll = 0;
          do
          {
            positions[step * i + ll] = ll * inInc;
            constants[step * i + ll] = gg[ll];
          } while (++ll < step);
        }
      }

      if (point >= minBounds && point <= maxBounds)
      {
        if (region == 0)
        { // entering the input extent
          region = 1;
          clipExt[2 * j] = i;
        }
      }
      else
      {
        if (region == 1)
        { // leaving the input extent
          region = 2;
          clipExt[2 * j + 1] = i - 1;
        }
      }
    }

    if (region == 0 || clipExt[2 * j] > clipExt[2 * j + 1])
    { // never entered input extent!
      validClip = false;
    }
  }

  if (!validClip)
  {
    // output extent doesn't itersect input extent
    for (int j = 0; j < 3; j++)
    {
      clipExt[2 * j] = outExt[2 * j];
      clipExt[2 * j + 1] = outExt[2 * j] - 1;
    }
  }
}

//------------------------------------------------------------------------------
} // ends anonymous namespace

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::GetInterpolationFunc(
  void (**func)(vtkInterpolationInfo*, const double[3], double*))
{
  vtkImageBSplineInterpolatorGetInterpolationFunc(
    func, this->InterpolationInfo->ScalarType, this->SplineDegree);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::GetInterpolationFunc(
  void (**func)(vtkInterpolationInfo*, const float[3], float*))
{
  vtkImageBSplineInterpolatorGetInterpolationFunc(
    func, this->InterpolationInfo->ScalarType, this->SplineDegree);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::GetRowInterpolationFunc(
  void (**func)(vtkInterpolationWeights*, int, int, int, double*, int))
{
  vtkImageBSplineInterpolatorGetRowInterpolationFunc(
    func, this->InterpolationInfo->ScalarType, this->SplineDegree);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::GetRowInterpolationFunc(
  void (**func)(vtkInterpolationWeights*, int, int, int, float*, int))
{
  vtkImageBSplineInterpolatorGetRowInterpolationFunc(
    func, this->InterpolationInfo->ScalarType, this->SplineDegree);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::PrecomputeWeightsForExtent(
  const double matrix[16], const int extent[6], int newExtent[6], vtkInterpolationWeights*& weights)
{
  weights = new vtkInterpolationWeights(*this->InterpolationInfo);

  vtkImageBSplineInterpolatorPrecomputeWeights(
    matrix, extent, newExtent, this->StructuredBoundsDouble, weights);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::PrecomputeWeightsForExtent(
  const float matrix[16], const int extent[6], int newExtent[6], vtkInterpolationWeights*& weights)
{
  weights = new vtkInterpolationWeights(*this->InterpolationInfo);

  vtkImageBSplineInterpolatorPrecomputeWeights(
    matrix, extent, newExtent, this->StructuredBoundsFloat, weights);
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::FreePrecomputedWeights(vtkInterpolationWeights*& weights)
{
  this->Superclass::FreePrecomputedWeights(weights);
}

//------------------------------------------------------------------------------
// build any tables required for the interpolation
void vtkImageBSplineInterpolator::BuildKernelLookupTable()
{
#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
  if (this->KernelLookupTable)
  {
    this->FreeKernelLookupTable();
  }

  float* kernel = 0;

  // kernel parameters
  int m = this->SplineDegree + 1;

  // compute lookup table size and step size
  int n = VTK_BSPLINE_KERNEL_TABLE_DIVISIONS;
  int size = n * m / 2;

  // allocate a little extra space
  kernel = new float[size + 4];

  // compute the table
  vtkBSplineKernel::BSpline(kernel, n, m - 1);

  // pad with a few zeros
  for (int j = 0; j < 4; j++)
  {
    kernel[size + j] = 0.0;
  }

  this->KernelLookupTable = kernel;
#endif
}

//------------------------------------------------------------------------------
void vtkImageBSplineInterpolator::FreeKernelLookupTable()
{
#ifdef VTK_BSPLINE_USE_KERNEL_TABLE
  float* kernel = this->KernelLookupTable;
  delete[] kernel;
#endif
}