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

vec4 WrapToClosestAngle(vec4 cData, float rFixed, vec2 cValueRange)
{
  float rRangeSize = cValueRange.y - cValueRange.x;
  return cData - floor((cData - rFixed) / rRangeSize + 0.5) * rRangeSize;
}

float WrapToRange(float rData, vec2 cValueRange)
{
  float rRangeSize = cValueRange.y - cValueRange.x;
  return rData - floor((rData - cValueRange.x) / rRangeSize) * rRangeSize;
}

float SampleAngular2D(sampler2D s, vec2 cPos, vec2 cValueRange)
{
  vec2 cSize = vec2(textureSize(s, 0));
  vec2 cUV = clamp(cPos * cSize, vec2(0.5), cSize - 0.5);
  ivec2 cNearestUV = ivec2(floor(cUV));
  ivec2 cFirstUV = ivec2(floor(cUV - 0.5)); // First of two sample points in each dimension.

  float rNearest = texelFetch(s, cNearestUV, 0).x;

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

  vec2 cWeightsU = GetLinearWeights(cFracUV.x);
  vec2 cWeightsV = GetLinearWeights(cFracUV.y);

  vec4 cData = vec4(texelFetchOffset(s, cFirstUV, 0, ivec2(0, 0)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(1, 0)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(0, 1)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(1, 1)).x);
  vec4 cWeights = cWeightsV.xxyy * cWeightsU.xyxy;
  cData = WrapToClosestAngle(cData, rNearest, cValueRange);
  float rDataSum = dot(cData, cWeights);

  return WrapToRange(rDataSum, cValueRange);
}

float SampleAngular3D(sampler3D s, vec3 cPos, vec2 cValueRange)
{
  vec3 cSize = vec3(textureSize(s, 0));
  vec3 cUVW = clamp(cPos * cSize, vec3(0.5), cSize - 0.5);
  ivec3 cNearestUVW = ivec3(floor(cUVW));
  ivec3 cFirstUVW = ivec3(floor(cUVW - 0.5)); // First of two sample points in each dimension.

  float rNearest = texelFetch(s, cNearestUVW, 0).x;

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

  vec2 cWeightsU = GetLinearWeights(cFracUVW.x);
  vec2 cWeightsV = GetLinearWeights(cFracUVW.y);
  vec2 cWeightsW = GetLinearWeights(cFracUVW.z);

  float rDataSum = 0.0;

  for (int k = 0; k < 2; ++k)
  {
    vec4 cData = vec4(texelFetch(s, cFirstUVW + ivec3(0, 0, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(1, 0, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(0, 1, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(1, 1, k), 0).x);
    vec4 cWeights = cWeightsW.x * cWeightsV.xxyy * cWeightsU.xyxy;
    cData = WrapToClosestAngle(cData, rNearest, cValueRange);
    rDataSum += dot(cData, cWeights);

    cWeightsW = cWeightsW.yx;
  }

  return WrapToRange(rDataSum, cValueRange);
}

float SampleAngular2DWithNoValue(sampler2D s, vec2 cPos, float rNoValue, vec2 cValueRange)
{
  vec2 cSize = vec2(textureSize(s, 0));
  vec2 cUV = clamp(cPos * cSize, vec2(0.5), cSize - 0.5);
  ivec2 cNearestUV = ivec2(floor(cUV));
  ivec2 cFirstUV = ivec2(floor(cUV - 0.5)); // First of two sample points in each dimension.

  float rNearest = texelFetch(s, cNearestUV, 0).x;

  if (rNearest == rNoValue) return rNoValue;

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

  vec2 cWeightsU = GetLinearWeights(cFracUV.x);
  vec2 cWeightsV = GetLinearWeights(cFracUV.y);

  vec4 cData = vec4(texelFetchOffset(s, cFirstUV, 0, ivec2(0, 0)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(1, 0)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(0, 1)).x,
                    texelFetchOffset(s, cFirstUV, 0, ivec2(1, 1)).x);
  vec4 cMask = vec4(notEqual(cData, vec4(rNoValue)));
  vec4 cWeights = cWeightsV.xxyy * cWeightsU.xyxy;
  cData = WrapToClosestAngle(cData, rNearest, cValueRange);
  float rDataSum = dot(cData, cWeights * cMask);
  float rWeightSum = dot(cWeights, cMask);

  return WrapToRange(rDataSum / rWeightSum, cValueRange);
}

float SampleAngular3DWithNoValue(sampler3D s, vec3 cPos, float rNoValue, vec2 cValueRange)
{
  vec3 cSize = vec3(textureSize(s, 0));
  vec3 cUVW = clamp(cPos * cSize, vec3(0.5), cSize - 0.5);
  ivec3 cNearestUVW = ivec3(floor(cUVW));
  ivec3 cFirstUVW = ivec3(floor(cUVW - 0.5)); // First of two sample points in each dimension.

  float rNearest = texelFetch(s, cNearestUVW, 0).x;

  if (rNearest == rNoValue) return rNoValue;

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

  vec2 cWeightsU = GetLinearWeights(cFracUVW.x);
  vec2 cWeightsV = GetLinearWeights(cFracUVW.y);
  vec2 cWeightsW = GetLinearWeights(cFracUVW.z);

  float rDataSum = 0.0;
  float rWeightSum = 0.0;

  for (int k = 0; k < 2; ++k)
  {
    vec4 cData = vec4(texelFetch(s, cFirstUVW + ivec3(0, 0, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(1, 0, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(0, 1, k), 0).x,
                      texelFetch(s, cFirstUVW + ivec3(1, 1, k), 0).x);
    vec4 cMask = vec4(notEqual(cData, vec4(rNoValue)));
    vec4 cWeights = cWeightsW.x * cWeightsV.xxyy * cWeightsU.xyxy;
    cData = WrapToClosestAngle(cData, rNearest, cValueRange);
    rDataSum += dot(cData, cWeights * cMask);
    rWeightSum += dot(cWeights, cMask);

    cWeightsW = cWeightsW.yx;
  }

  return WrapToRange(rDataSum / rWeightSum, cValueRange);
}
