#version 130
#extension GL_ARB_compatibility : enable
#extension GL_EXT_gpu_shader4 : enable

vec4 WrapToClosestAngle(vec4 value, float rFixed, vec2 range)
{
  float rangeSize = range.y - range.x;
  return value - floor((value - rFixed) / rangeSize + 0.5) * rangeSize;
}

float WrapToRange(float value, vec2 range)
{
  float rangeSize = range.y - range.x;
  return value - floor((value - range.x) / rangeSize) * rangeSize;
}

// Returns weights for two samples at texel centers.
vec2 GetLinearWeights(float t)
{
  return vec2(1.0 - t, t);
}

// Returns weights for four samples at texel centers.
vec4 GetCubicWeights(float t)
{
  float t2 = t * t,
        t3 = t2 * t;

  return vec4(-0.5 * t3 +       t2 - 0.5 * t,
               1.5 * t3 - 2.5 * t2 + 1.0,
              -1.5 * t3 + 2.0 * t2 + 0.5 * t,
               0.5 * (t3 - t2));
}

// 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 Get3PointCubicWeights(float t)
{
  float t2 = t * t;
  float t3 = t2 * t;
  return vec4(-0.5 * t3 + t2 - t * 0.5,
              (-3.0 * t3 + 4.0 * t2 + t) / (t - t2 + 2.0),
              0.5 * t3 - t2 + 0.5 * (t2 - t3) + 0.5 * t + 1.0,
              (t3 - t2) * 0.5);
}

float Sample1Bit(sampler3D s, vec3 coord, sampler2D lookup1Bit)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  float value = texelFetch(s, ivec3(uvw), 0).x;
  return texture(lookup1Bit, vec2(value, fract(uvw.x))).w;
}

float SampleNearest(sampler3D s, vec3 coord)
{
  return texelFetch(s, ivec3(coord * vec3(textureSize(s, 0))), 0).x;
}

int SampleNearest(isampler3D s, vec3 coord)
{
  return texelFetch(s, ivec3(coord * vec3(textureSize(s, 0))), 0).x;
}

float SampleLinear(sampler3D s, vec3 coord)
{
  return texture(s, coord).x;
}

float SampleLinear2DWithNoValue(sampler3D s, vec3 coord, float noValue)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw.xy - 0.5), nearestUVW.z); // First of two sample points in each dimension.

  if (texelFetch(s, nearestUVW, 0).x == noValue) return noValue;

  vec2 fracUV = uvw.xy - 0.5 - vec2(firstUVW.xy);

  vec2 weightsU = GetLinearWeights(fracUV.x);
  vec2 weightsV = GetLinearWeights(fracUV.y);

  vec4 value = vec4(texelFetchOffset(s, firstUVW, 0, ivec3(0, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(0, 1, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 1, 0)).x);
  vec4 mask = vec4(notEqual(value, vec4(noValue)));
  vec4 weights = weightsV.xxyy * weightsU.xyxy;
  float valueSum = dot(value, weights * mask);
  float weightSum = dot(weights, mask);

  return valueSum / weightSum;
}

float SampleLinear3DWithNoValue(sampler3D s, vec3 coord, float noValue)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw - 0.5)); // First of two sample points in each dimension.

  if (texelFetch(s, nearestUVW, 0).x == noValue) return noValue;

  vec3 fracUVW = uvw - 0.5 - vec3(firstUVW);

  vec2 weightsU = GetLinearWeights(fracUVW.x);
  vec2 weightsV = GetLinearWeights(fracUVW.y);
  vec2 weightsW = GetLinearWeights(fracUVW.z);

  float valueSum = 0.0;
  float weightSum = 0.0;

  for (int k = 0; k < 2; ++k)
  {
    vec4 value = vec4(texelFetch(s, firstUVW + ivec3(0, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(0, 1, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 1, k), 0).x);
    vec4 mask = vec4(notEqual(value, vec4(noValue)));
    vec4 weights = weightsW.x * weightsV.xxyy * weightsU.xyxy;
    valueSum += dot(value, weights * mask);
    weightSum += dot(weights, mask);

    weightsW = weightsW.yx;
  }

  return valueSum / weightSum;
}

float SampleAngular2D(sampler3D s, vec3 coord, vec2 angularRange)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw.xy - 0.5), nearestUVW.z); // First of two sample points in each dimension.

  float nearestValue = texelFetch(s, nearestUVW, 0).x;

  vec2 fracUV = uvw.xy - 0.5 - vec2(firstUVW.xy);

  vec2 weightsU = GetLinearWeights(fracUV.x);
  vec2 weightsV = GetLinearWeights(fracUV.y);

  vec4 value = vec4(texelFetchOffset(s, firstUVW, 0, ivec3(0, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(0, 1, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 1, 0)).x);
  vec4 weights = weightsV.xxyy * weightsU.xyxy;
  value = WrapToClosestAngle(value, nearestValue, angularRange);
  float valueSum = dot(value, weights);

  return WrapToRange(valueSum, angularRange);
}

float SampleAngular3D(sampler3D s, vec3 coord, vec2 angularRange)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw - 0.5)); // First of two sample points in each dimension.

  float nearestValue = texelFetch(s, nearestUVW, 0).x;

  vec3 fracUVW = uvw - 0.5 - vec3(firstUVW);

  vec2 weightsU = GetLinearWeights(fracUVW.x);
  vec2 weightsV = GetLinearWeights(fracUVW.y);
  vec2 weightsW = GetLinearWeights(fracUVW.z);

  float valueSum = 0.0;

  for (int k = 0; k < 2; ++k)
  {
    vec4 value = vec4(texelFetch(s, firstUVW + ivec3(0, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(0, 1, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 1, k), 0).x);
    vec4 weights = weightsW.x * weightsV.xxyy * weightsU.xyxy;
    value = WrapToClosestAngle(value, nearestValue, angularRange);
    valueSum += dot(value, weights);

    weightsW = weightsW.yx;
  }

  return WrapToRange(valueSum, angularRange);
}

float SampleAngular2DWithNoValue(sampler3D s, vec3 coord, float noValue, vec2 angularRange)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw.xy - 0.5), nearestUVW.z); // First of two sample points in each dimension.

  float nearestValue = texelFetch(s, nearestUVW, 0).x;

  if (nearestValue == noValue) return noValue;

  vec2 fracUV = uvw.xy - 0.5 - vec2(firstUVW.xy);

  vec2 weightsU = GetLinearWeights(fracUV.x);
  vec2 weightsV = GetLinearWeights(fracUV.y);

  vec4 value = vec4(texelFetchOffset(s, firstUVW, 0, ivec3(0, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 0, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(0, 1, 0)).x,
                    texelFetchOffset(s, firstUVW, 0, ivec3(1, 1, 0)).x);
  vec4 mask = vec4(notEqual(value, vec4(noValue)));
  vec4 weights = weightsV.xxyy * weightsU.xyxy;
  value = WrapToClosestAngle(value, nearestValue, angularRange);
  float valueSum = dot(value, weights * mask);
  float weightSum = dot(weights, mask);

  return WrapToRange(valueSum / weightSum, angularRange);
}

float SampleAngular3DWithNoValue(sampler3D s, vec3 coord, float noValue, vec2 angularRange)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw - 0.5)); // First of two sample points in each dimension.

  float nearestValue = texelFetch(s, nearestUVW, 0).x;

  if (nearestValue == noValue) return noValue;

  vec3 fracUVW = uvw - 0.5 - vec3(firstUVW);

  vec2 weightsU = GetLinearWeights(fracUVW.x);
  vec2 weightsV = GetLinearWeights(fracUVW.y);
  vec2 weightsW = GetLinearWeights(fracUVW.z);

  float valueSum = 0.0;
  float weightSum = 0.0;

  for (int k = 0; k < 2; ++k)
  {
    vec4 value = vec4(texelFetch(s, firstUVW + ivec3(0, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 0, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(0, 1, k), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, 1, k), 0).x);
    vec4 mask = vec4(notEqual(value, vec4(noValue)));
    vec4 weights = weightsW.x * weightsV.xxyy * weightsU.xyxy;
    value = WrapToClosestAngle(value, nearestValue, angularRange);
    valueSum += dot(value, weights * mask);
    weightSum += dot(weights, mask);

    weightsW = weightsW.yx;
  }

  return WrapToRange(valueSum / weightSum, angularRange);
}

float SampleCubic2D(sampler3D s, vec3 coord)
{
  vec3 size = vec3(textureSize(s, 0));
  vec3 sizeInv = 1.0 / size;
  vec3 uvw = coord * size;
  vec3 secondUVW = vec3(floor(uvw.xy - 0.5) + 0.5, uvw.z);
  vec2 fracUV = uvw.xy - secondUVW.xy;
  coord = secondUVW * sizeInv;

  vec4 threePointWeightsX = Get3PointCubicWeights(fracUV.x);
  vec4 threePointWeightsY = Get3PointCubicWeights(fracUV.y);

  threePointWeightsX.y *= sizeInv.x;
  threePointWeightsY.y *= sizeInv.y;

  float valueSum = 0.0;

  vec3 offsetsY = vec3(-sizeInv.y, threePointWeightsY.y, 2.0 * sizeInv.y);
  vec3 weightsY = threePointWeightsY.xzw;

  for (int j = 0; j < 3; ++j)
  {
    valueSum += dot(vec3(texture(s, coord + vec3(-sizeInv.x, offsetsY.x, 0.0)).x,
                         texture(s, coord + vec3(threePointWeightsX.y, offsetsY.x, 0.0)).x,
                         texture(s, coord + vec3(2.0 * sizeInv.x, offsetsY.x, 0.0)).x), threePointWeightsX.xzw) * weightsY.x;

    weightsY = weightsY.yzx;
    offsetsY = offsetsY.yzx;
  }

  return valueSum;
}

float SampleCubic3D(sampler3D s, vec3 coord)
{
  vec3 size = vec3(textureSize(s, 0));
  vec3 sizeInv = 1.0 / size;
  vec3 uvw = coord * size;
  vec3 secondUVW = floor(uvw - 0.5) + 0.5;
  vec3 fracUVW = uvw - secondUVW;
  coord = secondUVW * sizeInv;

  vec4 threePointWeightsX = Get3PointCubicWeights(fracUVW.x);
  vec4 threePointWeightsY = Get3PointCubicWeights(fracUVW.y);
  vec4 threePointWeightsZ = Get3PointCubicWeights(fracUVW.z);

  threePointWeightsX.y *= sizeInv.x;
  threePointWeightsY.y *= sizeInv.y;
  threePointWeightsZ.y *= sizeInv.z;

  vec3 offsetsZ = vec3(-sizeInv.z, threePointWeightsZ.y, 2.0 * sizeInv.z);
  vec3 weightsZ = threePointWeightsZ.xzw;

  float valueSum = 0.0;

  for (int k = 0; k < 3; ++k)
  {
    vec3 offsetsY = vec3(-sizeInv.y, threePointWeightsY.y, 2.0 * sizeInv.y);
    vec3 weightsY = threePointWeightsY.xzw * weightsZ.x;

    for (int j = 0; j < 3; ++j)
    {
      valueSum += dot(vec3(texture(s, coord + vec3(-sizeInv.x, offsetsY.x, offsetsZ.x)).x,
                           texture(s, coord + vec3(threePointWeightsX.y, offsetsY.x, offsetsZ.x)).x,
                           texture(s, coord + vec3(2.0 * sizeInv.x, offsetsY.x, offsetsZ.x)).x), threePointWeightsX.xzw) * weightsY.x;

      weightsY = weightsY.yzx;
      offsetsY = offsetsY.yzx;
    }

    weightsZ = weightsZ.yzx;
    offsetsZ = offsetsZ.yzx;
  }

  return valueSum;
}

float SampleCubic2DWithNoValue(sampler3D s, vec3 coord, float noValue)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw.xy - 1.5), nearestUVW.z); // First of four sample points in each dimension.

  if (texelFetch(s, nearestUVW, 0).x == noValue) return noValue;

  vec2 fracUV = uvw.xy - 1.5 - vec2(firstUVW.xy);

  vec4 weightsU = GetCubicWeights(fracUV.x);
  vec4 weightsV = GetCubicWeights(fracUV.y);

  float valueSum = 0.0;
  float weightSum = 0.0;

  for (int j = 0; j < 4; ++j)
  {
    vec4 value = vec4(texelFetch(s, firstUVW + ivec3(0, j, 0), 0).x,
                      texelFetch(s, firstUVW + ivec3(1, j, 0), 0).x,
                      texelFetch(s, firstUVW + ivec3(2, j, 0), 0).x,
                      texelFetch(s, firstUVW + ivec3(3, j, 0), 0).x);
    vec4 mask = vec4(notEqual(value, vec4(noValue)));
    vec4 weights = weightsV.x * weightsU;
    valueSum += dot(value, weights * mask);
    weightSum += dot(weights, mask);

    weightsV = weightsV.yzwx;
  }

  return valueSum / weightSum;
}

float SampleCubic3DWithNoValue(sampler3D s, vec3 coord, float noValue)
{
  vec3 uvw = coord * vec3(textureSize(s, 0));
  ivec3 nearestUVW = ivec3(floor(uvw));
  ivec3 firstUVW = ivec3(floor(uvw - 1.5)); // First of four sample points in each dimension.

  if (texelFetch(s, nearestUVW, 0).x == noValue) return noValue;

  vec3 fracUVW = uvw - 1.5 - vec3(firstUVW);

  vec4 weightsU = GetCubicWeights(fracUVW.x);
  vec4 weightsV = GetCubicWeights(fracUVW.y);
  vec4 weightsW = GetCubicWeights(fracUVW.z);

  float valueSum = 0.0;
  float weightSum = 0.0;

  for (int k = 0; k < 4; ++k)
  {
    for (int j = 0; j < 4; ++j)
    {
      vec4 value = vec4(texelFetch(s, firstUVW + ivec3(0, j, k), 0).x,
                        texelFetch(s, firstUVW + ivec3(1, j, k), 0).x,
                        texelFetch(s, firstUVW + ivec3(2, j, k), 0).x,
                        texelFetch(s, firstUVW + ivec3(3, j, k), 0).x);
      vec4 mask = vec4(notEqual(value, vec4(noValue)));
      vec4 weights = weightsW.x * weightsV.x * weightsU;
      valueSum += dot(value, weights * mask);
      weightSum += dot(weights, mask);

      weightsV = weightsV.yzwx;
    }
    weightsW = weightsW.yzwx;
  }

  return valueSum / weightSum;
}
