   // this initialization corresponds to mirror boundaries
   // accelerated loop
   float zn = IU_BSPLINE_POLE;
-  floatN Sum = *c;
+  PixelType Sum = *c;
   for (unsigned int n = 1; n < Horizon; n++) {
     c += step;
     Sum += zn * *c;

  return(Sum);
}

template<class floatN> __device__ floatN InitialAntiCausalCoefficient(
  floatN* c,        // last coefficient
  unsigned int DataLength,  // number of samples or coefficients
  int step)
{
  // this initialization corresponds to mirror boundaries
  return((IU_BSPLINE_POLE / (IU_BSPLINE_POLE * IU_BSPLINE_POLE - 1.0f)) * (IU_BSPLINE_POLE * c[-step] + *c));
}

template<class floatN> __device__ void ConvertToInterpolationCoefficients(
  floatN* coeffs,   // input samples --> output coefficients
  unsigned int DataLength,  // number of samples or coefficients
  int step)
{
  // compute the overall gain
  const float Lambda = (1.0f - IU_BSPLINE_POLE) * (1.0f - 1.0f / IU_BSPLINE_POLE);

  // causal initialization
  floatN* c = coeffs;
  floatN previous_c;  //cache the previously calculated c rather than look it up again (faster!)
  *c = previous_c = Lambda * InitialCausalCoefficient(c, DataLength, step);
  // causal recursion
  for (unsigned int n = 1; n < DataLength; ++n) {
    c += step;
    *c = previous_c = Lambda * (*c) + IU_BSPLINE_POLE * previous_c;
  }

  // anticausal initialization

#undef IU_BSPLINE_POLE

//-----------------------------------------------------------------------------
template<class PixelType> __global__ void cuSamplesToCoefficients2DX(
  PixelType* image, unsigned int width, unsigned int height, size_t stride)
{
  const unsigned int y = blockIdx.x * blockDim.x + threadIdx.x;
  if (y < height)
  {
    // process lines in x-direction
    PixelType* line = image + y * stride;  //direct access

    ConvertToInterpolationCoefficients(line, width, 1);
  }
}

//-----------------------------------------------------------------------------
template<class PixelType> __global__ void cuSamplesToCoefficients2DY(
  PixelType* image, unsigned int width, unsigned int height, size_t stride)
{
  const unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
  if (x < width)
  {
    // process lines in y-direction
    PixelType* line = image + x;  //direct access

    ConvertToInterpolationCoefficients(line, height, stride);
  }
