
#define EPSILON 0.00001

uniform sampler2D solid;
uniform sampler2D back;
uniform sampler2D front;
uniform sampler2D front_color;
uniform sampler2D occlusion;
uniform vec3 origin;
uniform vec3 imageSize;
uniform ivec3 texSize;
uniform float slices;
uniform ivec3 has_clipped;

varying vec3 objectPos;

vec4 correctColor(vec4 color, float thickness)
{
  float eqa = (color.a < 1.0) ? ((pow(color.a, thickness)-1.0)/(color.a-1.0)) : 1.0;
  color.rgb *= eqa;
  color.a = pow(color.a, thickness);
  return color;
}

vec4 blendColor(vec4 src, vec4 dst)
{
  return src.a*src + (1-src.a)*dst;
}

vec4 backBlend(vec4 src, vec4 dst)
{
  vec4 ret;
  ret.rgb = src.a*src.rgb*(1-dst.a);
  ret.a = (1-dst.a)*(1-src.a);
  return ret;
}

vec4 blendPremulColor(vec4 src, vec4 tgt)
{
  vec4 ret;
  ret.a = src.a*tgt.a;
  ret.rgb = src.rgb + src.a * tgt.rgb;
  return ret;
}

float equivAlpha(float a0, float a1)
{
  return a0+a1-a0*a1;
}

vec4 equivColor(vec4 far, vec4 near)
{
  float ea = equivAlpha(far.a, near.a);
  vec4 ret;
  ret.rgb = far.rgb + (near.a*near.rgb - near.a*near.rgb)/ea;
  ret.a = ea;
  return ret;
}

void main()
{
  float alpha = 1.0;

  vec3 bbox_max = origin + imageSize;
  vec3 bbox_min = origin;
  vec3 pos = objectPos; // A point on the ray

  vec3 viewDir = vec3(gl_ModelViewMatrix[0][2], gl_ModelViewMatrix[1][2], gl_ModelViewMatrix[2][2]);
  viewDir /= -1.0*length(viewDir); // Orientation of the ray

  // Find intersection of the current ray with the cube, in image coordinate system
  vec2 lspan = vec2(-1e20, 1e20);

  for(int i = 0 ; i < 3 ; ++i)
  {
    if(fabs(viewDir[i]) > 1e-8)
    {
      bool inc = viewDir[i] > 0;
      float bb_min = inc ? bbox_min[i] : bbox_max[i];
      float bb_max = inc ? bbox_max[i] : bbox_min[i];
      float min_intersection = (bb_min - pos[i]) / viewDir[i];
      float max_intersection = (bb_max - pos[i]) / viewDir[i];
      if(min_intersection > lspan.x) lspan.x = min_intersection;
      if(max_intersection < lspan.y) lspan.y = max_intersection;
    }
    else
    {
      if(pos[i] > bbox_max[i] ||
         pos[i] < bbox_min[i])
      {
        discard;
        return;
      }
    }
  }

// Compute intersection with clipping planes
  if(lspan.x > lspan.y) { discard; return; }
  else
  {
    bool to_draw = true;
    vec3 startPos = pos + lspan.x * viewDir;
    vec3 endPos = pos + lspan.y * viewDir;

    vec4 realStart = gl_ModelViewMatrix * vec4(startPos, 1);
    vec4 realEnd = gl_ModelViewMatrix * vec4(endPos, 1);

    realStart /= realStart.w;
    realEnd /= realEnd.w;

    for(int i = 0 ; i < 3 ; ++i)
    {
      if(has_clipped[i] != 0)
      {
        vec4 c0 = gl_ClipPlane[2*i];
        vec4 c1 = gl_ClipPlane[2*i+1];

        float testS0 = dot(c0, realStart);
        float testS1 = dot(c1, realStart);
        float testE0 = dot(c0, realEnd);
        float testE1 = dot(c1, realEnd);

        if((testS0 < 0 && testE0 < 0) ||
           (testS1 < 0 && testE1 < 0))
        {
          discard;
          return;
        }
        else if((testS0 * testE0 < 0) ||
                (testS1 * testE1 < 0))
        {
          float lambda0 = testS0 / (testS0 - testE0);
          float lambda1 = testS1 / (testS1 - testE1);
          if(lambda0 > lambda1)
          {
            float a = lambda0;
            lambda0 = lambda1;
            lambda1 = a;
          }
          vec4 newStart = (lambda0 > 0) ? (1-lambda0)*realStart + lambda0 * realEnd : realStart;
          vec4 newEnd = (lambda1 < 1) ? (1-lambda1)*realStart + lambda1 * realEnd : realEnd;
          realStart = newStart;
          realEnd = newEnd;
        }
      }
    }

    vec4 startPos4 = gl_ModelViewMatrixInverse * realStart;
    vec4 endPos4 = gl_ModelViewMatrixInverse * realEnd;
    startPos = startPos4.xyz / startPos4.w;
    endPos = endPos4.xyz / endPos4.w;

    vec4 screenPos = gl_ModelViewProjectionMatrix * startPos4;
    screenPos.xy = (screenPos.xy / screenPos.w + 1.0) / 2.0;

    float front_alpha = texture2D(front_color, screenPos.xy).a;

    {
      bool occluded = texture2D(occlusion, screenPos.xy).a == 0;

      if(occluded || front_alpha > 0.99)
      {
        discard;
        return;
      }
    }

    {
      float depth_back = texture2D(back, screenPos.xy).a;
      float depth_solid = texture2D(solid, screenPos.xy).a;
      float max_depth = min(depth_back, depth_solid);

      vec4 realPos = screenPos;
      realPos.z = max_depth;
      realPos.xyz = 2.0*realPos.xyz - 1.0;
      realPos.w = 1.0;
      vec4 limitPos4 = gl_ModelViewProjectionMatrixInverse * realPos;

      vec3 limitPos = limitPos4.xyz / limitPos4.w;
      vec3 dv = endPos - startPos;

      float dv2 = dot(dv, dv);
      float dist = dot(limitPos - startPos, dv);

      if(dist < 0)
      {
        discard;
        return;
      }
      else if(dist < dv2)
      {
        endPos = limitPos;
      }
    }

    {
      vec4 realPos = screenPos;
      realPos.z = texture2D(front, screenPos.xy).a;
      realPos.xyz = 2.0*realPos.xyz - 1.0;
      realPos.w = 1.0;
      vec4 limitPos4 = gl_ModelViewProjectionMatrixInverse * realPos;

      vec3 limitPos = limitPos4.xyz / limitPos4.w;
      vec3 dv = startPos - endPos;

      float dv2 = dot(dv, dv);
      float dist = dot(limitPos - endPos, dv);

      if(dist < 0)
      {
        discard;
        return;
      }
      else if(dist < dv2)
      {
        startPos = limitPos;
      }
    }

    {
      vec3 startTex = (startPos - origin) / imageSize;
      vec3 endTex = (endPos - origin) / imageSize;

      vec3 dTex = (endTex - startTex);
      vec3 dPos = (endPos - startPos);

      /*float squares = dot(vabs(dTex), texSize);
       */
      vec3 texStep = imageSize / vec3(texSize);
      float unitTexSize = min(texStep.x, min(texStep.y, texStep.z));
      vec3 cs = dPos / unitTexSize;
      float squares = max(abs(cs.x), max(abs(cs.y), abs(cs.z)));

      int nb_slices = int(8*ceil(squares*slices));
      if(nb_slices < 1) nb_slices = 1;

      vec3 curTex = startTex;// + start*dTex;
      vec4 c = vec4(0,0,0,1);// color_premul(curTex);
      vec3 curPos = startPos;// + start*dPos;
      dTex /= float(nb_slices);
      dPos /= float(nb_slices);

      curPos += dPos/2;

      float thickness = length(dPos) / unitTexSize; // squares/float(nb_slices);

      front_alpha = 1-front_alpha;

      for(int n = 0 ; n < nb_slices ; ++n)
      {
        vec4 nc = color(curTex);
        vec4 rc = correctColor(nc, thickness);
        c = blendPremulColor(c, rc);
        /*c = maxColor(c, nc);*/
        /*c = equivColor(c, nc);*/
        /*c = blendColor(nc, c);*/
        /*c = backBlend(nc, c);*/
        /*if(c.a > .95) break;*/
        front_alpha *= rc.a;
        if(front_alpha < 0.01) { break; }
        curPos += dPos;
        curTex += dTex;
      }
      c.a = 1-c.a;
      if(c.a > 0)
        c.rgb /= c.a;
      gl_FragColor = c;
    }
  }
}
