// Returns a vector of the weight of the first sample, the location and weight of the middle sample, and weight of the last sample.
vec4 Cubic_GetWeights(float x)
{
  float x2 = x * x;
  float x3 = x2 * x;
  return vec4(
    -0.5 * x3 + x2 - x * 0.5,
    (-3.0 * x3 + 4.0 * x2 + x) / (x - x2 + 2.0),
    0.5 * x3 - x2 + 0.5 * (x2 - x3) + 0.5 * x + 1.0,
    (x3 - x2) * 0.5
  );
}

float SampleBicubic2D(sampler2D s, vec2 cPos)
{
  vec2 cSize = vec2(textureSize(s, 0));
  vec2 cSizeInv = 1.0 / cSize;
  vec2 cTmp = floor(cPos * cSize - vec2(0.5)) + vec2(0.5);
  vec2 t = fract(cPos * cSize - 0.5);
  cPos = cTmp * cSizeInv;

  vec4 cFastWeightsX = Cubic_GetWeights(t.x),
       cFastWeightsY = Cubic_GetWeights(t.y);
  cFastWeightsX.y *= cSizeInv.x;
  vec3 cOffsetsY = vec3(-cSizeInv.y, cFastWeightsY.y * cSizeInv.y, 2.0 * cSizeInv.y);
  vec3 cWeightsY = cFastWeightsY.xzw;
  float rValue = 0.0;
  for (int i = 0; i < 3; ++i)
  {
     rValue += (texture2D(s, cPos + vec2(-cSizeInv.x,      cOffsetsY.x)).r * cFastWeightsX.x +
                texture2D(s, cPos + vec2(cFastWeightsX.y,  cOffsetsY.x)).r * cFastWeightsX.z +
                texture2D(s, cPos + vec2(2.0 * cSizeInv.x, cOffsetsY.x)).r * cFastWeightsX.w) * cWeightsY.x;
     cWeightsY = cWeightsY.yzx;
     cOffsetsY = cOffsetsY.yzx;
  }
  return rValue;
}


float SampleTricubic3D(sampler3D s, vec3 cPos)
{
  vec3 cSize = vec3(textureSize(s, 0));
  vec3 cSizeInv = 1.0 / cSize;
  vec3 cTmp = floor(cPos * cSize - vec3(0.5)) + vec3(0.5);
  vec3 t = fract(cPos * cSize - vec3(0.5));
  cPos = cTmp * cSizeInv;

  vec4 cFastWeightsX = Cubic_GetWeights(t.x),
       cFastWeightsY = Cubic_GetWeights(t.y),
       cFastWeightsZ = Cubic_GetWeights(t.z);
  cFastWeightsX.y *= cSizeInv.x;
  cFastWeightsY.y *= cSizeInv.y;
  cFastWeightsZ.y *= cSizeInv.z;

  vec3 cOffsetsZ = vec3(-cSizeInv.z, cFastWeightsZ.y, 2.0 * cSizeInv.z);
  vec3 cWeightsZ = cFastWeightsZ.xzw;

  float rValue = 0.0;
  for (int i = 0; i < 3; ++i)
  {
    vec3 cOffsetsY = vec3(-cSizeInv.y, cFastWeightsY.y, 2.0 * cSizeInv.y);
    vec3 cWeightsY = cFastWeightsY.xzw * cWeightsZ.x;
    for (int j = 0; j < 3; ++j)
    {
      rValue += (texture3D(s, cPos + vec3(-cSizeInv.x,      cOffsetsY.x, cOffsetsZ.x)).r * cFastWeightsX.x +
                 texture3D(s, cPos + vec3(cFastWeightsX.y,  cOffsetsY.x, cOffsetsZ.x)).r * cFastWeightsX.z +
                 texture3D(s, cPos + vec3(2.0 * cSizeInv.x, cOffsetsY.x, cOffsetsZ.x)).r * cFastWeightsX.w) * cWeightsY.x;
      cWeightsY = cWeightsY.yzx;
      cOffsetsY = cOffsetsY.yzx;
    }
    cWeightsZ = cWeightsZ.yzx;
    cOffsetsZ = cOffsetsZ.yzx;
  }
  return rValue;
}

// Returns weights for four samples at texel centers.
vec4 GetCubicWeights(float rT)
{
  float rT2 = rT * rT,
        rT3 = rT2 * rT;

  return vec4(-0.5 * rT3 +       rT2 - 0.5 * rT,
               1.5 * rT3 - 2.5 * rT2 + 1.0,
              -1.5 * rT3 + 2.0 * rT2 + 0.5 * rT,
               0.5 * (rT3 - rT2));
}

float SampleBicubic2DWithNoValue(sampler2D s, vec2 cPos, float rNoValue)
{
  ivec2 cSize = textureSize(s, 0);
  vec2 cUV = cPos * vec2(cSize);
  ivec2 cNearestUV = ivec2(floor(cUV));
  ivec2 cFirstUV = ivec2(floor(cUV - 1.5)); // First of four sample points in each dimension.

  if (texelFetch(s, cNearestUV, 0).x == rNoValue) return rNoValue;

  vec2 cFracUV = cUV - 1.5 - vec2(cFirstUV);

  vec4 cWeightsU = GetCubicWeights(cFracUV.x);
  vec4 cWeightsV = GetCubicWeights(cFracUV.y);

  float rDataSum = 0.0;
  float rWeightSum = 0.0;

  for (int j = 0; j < 4; ++j)
  {
    vec4 cData = vec4(texelFetch(s, clamp(cFirstUV + ivec2(0, j), ivec2(0), cSize - 1), 0).x,
                      texelFetch(s, clamp(cFirstUV + ivec2(1, j), ivec2(0), cSize - 1), 0).x,
                      texelFetch(s, clamp(cFirstUV + ivec2(2, j), ivec2(0), cSize - 1), 0).x,
                      texelFetch(s, clamp(cFirstUV + ivec2(3, j), ivec2(0), cSize - 1), 0).x);
    vec4 cMask = vec4(notEqual(cData, vec4(rNoValue)));
    vec4 cWeights = cWeightsV.x * cWeightsU;
    rDataSum += dot(cData, cWeights * cMask);
    rWeightSum += dot(cWeights, cMask);

    cWeightsV = cWeightsV.yzwx;
  }

  return rDataSum / rWeightSum;
}

float SampleTricubic3DWithNoValue(sampler3D s, vec3 cPos, float rNoValue)
{
  ivec3 cSize = textureSize(s, 0);
  vec3 cUVW = cPos * vec3(cSize);
  ivec3 cNearestUVW = ivec3(floor(cUVW));
  ivec3 cFirstUVW = ivec3(floor(cUVW - 1.5)); // First of four sample points in each dimension.

  if (texelFetch(s, cNearestUVW, 0).x == rNoValue) return rNoValue;

  vec3 cFracUVW = cUVW - 1.5 - vec3(cFirstUVW);

  vec4 cWeightsU = GetCubicWeights(cFracUVW.x);
  vec4 cWeightsV = GetCubicWeights(cFracUVW.y);
  vec4 cWeightsW = GetCubicWeights(cFracUVW.z);

  float rDataSum = 0.0;
  float rWeightSum = 0.0;

  for (int k = 0; k < 4; ++k)
  {
    for (int j = 0; j < 4; ++j)
    {
      vec4 cData = vec4(texelFetch(s, clamp(cFirstUVW + ivec3(0, j, k), ivec3(0), cSize - 1), 0).x,
                        texelFetch(s, clamp(cFirstUVW + ivec3(1, j, k), ivec3(0), cSize - 1), 0).x,
                        texelFetch(s, clamp(cFirstUVW + ivec3(2, j, k), ivec3(0), cSize - 1), 0).x,
                        texelFetch(s, clamp(cFirstUVW + ivec3(3, j, k), ivec3(0), cSize - 1), 0).x);
      vec4 cMask = vec4(notEqual(cData, vec4(rNoValue)));
      vec4 cWeights = cWeightsW.x * cWeightsV.x * cWeightsU;
      rDataSum += dot(cData, cWeights * cMask);
      rWeightSum += dot(cWeights, cMask);

      cWeightsV = cWeightsV.yzwx;
    }
    cWeightsW = cWeightsW.yzwx;
  }

  return rDataSum / rWeightSum;
}
