/*
 * levoy.sl -- Implement Levoy volume rendering in Shading Language
 *
 * DESCRIPTION:
 *   Believe it or not, this surface renders volume data, like from a
 *   CT scan.  It's a "magic" surface which actually invokes ray marching
 *   to determine how the volumetric data _behind_ the surface should
 *   look.  Yowza!
 *
 * INSTRUCTIONS:
 *   First, convert your volume data into one BIG texture (TIFF files for
 *   rendrib, or use txmake for prman).  Tile the images into a big square
 *   like this:
 *
 *     +---------+---------+---------+
 *     |         |         |         |
 *     | slice 0 | slice 1 | slice 2 |
 *     |         |         |         |
 *     +---------+---------+---------+
 *     |         |         |         |
 *     | slice 3 | slice 4 | slice 5 |
 *     |         |         |         |
 *     +---------+---------+---------+
 *     |         |         |         |
 *     | slice 6 | slice 7 | slice 8 |
 *     |         |         |         |
 *     +---------+---------+---------+
 *
 *   Obviously, you aren't limited to 9 slices.  If N = your actual
 *   number of CT slices, then you need to tile n x n, where
 *   n = ceil(sqrt(N)).  Just leave the last ones blank.
 *
 *   In your scene, make a box -- a unit cube centered at the origin
 *   of its shader space.  So in shader space, it spans [-.5,.5] in
 *   in x, y, and z.  Scale the box appropriately to get the slice
 *   spacing correct, or to otherwise scale or warp the volume.
 *
 *   For the six sides of the box, use this surface shader, "levoy".
 *   Be sure to set xpanels and ypanels appropriately.  The stepsize
 *   parameter may also need fiddling.
 *
 *   There are two ways to shade.  If you want something that looks
 *   like an isosurface, set the "isovalue" parameter to the value of
 *   the isosurface, "isovaluealpha" to the opacity, and "isovaluethickness"
 *   to the apparent thickness of the surface.  These parameters are
 *   explained in Levoy's paper.
 *
 *   If you want a more "jello-like" appearance rather than a single
 *   isosurface, leave "isosurface" alone (or set it to something <= 0),
 *   and set the various mat?_val, mat?_alpha, mat?_color for the various
 *   materials you expect in the volume.  Again, refer to Levoy's paper
 *   for more details on what this all means.
 *
 *   Then just render the image.  Wherever the renderer sees this surface
 *   (i.e. your box), it will shade the pixel appropriately by ray marching
 *   in the direction of I, looking up volume values as it goes. Shading
 *   Language has no 3-D volumetric texture mapping lookup, which is why
 *   we tiled the volume onto one big texture image -- that lets us use SL's
 *   2D texture mapping functions.
 *
 *   Since there is nothing special about this volume object, you can put
 *   it in a scene with any other (volumetric or non-volumetric) objects,
 *   and just let the renderer go at it!
 *
 *
 *
 * PARAMETERS
 *   Ka, Kd, Ks, roughness, specularcolor - control shading like "plastic"
 *   stepsize - step size (in shader space) for ray marching.  Smaller is
 *        more accurate, but also slower.
 *   filename - the name of the file containing your volume data.
 *   xpanels, ypanels - give the configuration of the slice data in your
 *        texture file.
 *   isovalue, isovaluealpha, isovaluethickness - control appearance in the
 *        case where you want to do an isovalue (only if isovalue > 0).
 *   mat?_val, mat?_alpha, mat?_color - control the materials when the
 *        multiple material option is being used (if isovalue < 0).
 *   nmaterials - total number of materials.
 *   debug - when nonzero, prints lots of useless debugging information.
 *
 *
 *
 * AUTHOR: Larry Gritz, 1995
 *         gritz@seas.gwu.edu
 *         The George Washington University
 *         801 22nd St. Rm. T-624, Washington DC 20052
 *
 *
 * REFERENCES:
 *   Levoy, Marc.  "Display of Surfaces From Volume Data", IEEE CG&A,
 *       volume 8, number 3, pp. 29-37, May 1988.
 *
 * WARNINGS:
 *   I suspect this will not work properly under the current version of
 *   PRMan.
 *
 * HINTS:
 *   For best results, turn off light caching in BMRT by putting this in
 *      front of any LightSource declarations:
 *              Attribute "light" "cache" "off"
 *
 *
 * HISTORY:
 *   16 March 95 (lg) - finally got everything working
 *
 *   last modified 16 March 95
 */





#define DEB if (debug > 0)


#define linterp(val,val0,val1,r0,r1) (r1 * (val - val0) / (val1 - val0) + \
          r0 * (val1 - val) / (val1 - val0))


#define lookup_data(P,val) \
  Ptemp=(P); \
  z = floor((zcomp(Ptemp))*(xpanels*ypanels)); \
  yp = floor(z/xpanels)/ypanels; \
  xp = (z - xpanels*floor(z/xpanels))/xpanels; \
  x = xp + (xcomp(Ptemp))/xpanels; \
  y = yp + (ycomp(Ptemp))/ypanels; \
  val = texture(filename,x,y,x,y,x,y,x,y)



#define lookup_data_sameslice(P,val) \
  Ptemp=(P); \
  x = xp + (xcomp(Ptemp))/xpanels; \
  y = yp + (ycomp(Ptemp))/ypanels; \
  val = texture(filename,x,y,x,y,x,y,x,y)


#define calc_interpolate_gradient(P,file,G,stepsize,here)  \
  interpolate_data_sameslice(P+point(stepsize,0,0),v0); \
  interpolate_data_sameslice(P+point(0,stepsize,0),v1); \
  interpolate_data(P+point(0,0,stepsize),v2); \
  G = (0.25 / stepsize) * (point (v0 - here, v1-here, v2-here))


#define calc_gradient(P,file,G,stepsize,here)  \
  lookup_data_sameslice(P+point(stepsize,0,0),v0); \
  lookup_data_sameslice(P+point(0,stepsize,0),v1); \
  lookup_data(P+point(0,0,stepsize),v2); \
  G = (0.25 / stepsize) * (point (v0 - here, v1-here, v2-here))



/* #define access_volume(P,val) lookup_data(P,val) */



surface
levoy3 (float  Ka = 1, Kd = 0.5, Ks = 0.5;
       float  roughness = 0.02;
       color  specularcolor = 1;
       float  stepsize = 0.005;
       float  gradientstepsize = 0.01;
       string filename = "";
       float  xpanels = 8, ypanels = 8;
       float  isovalue = -1, isovaluealpha = 10, isovaluethickness = 0.1;
       float  nmaterials = 0;
       float  mat0_val = 0, mat0_alpha = 0;
       color  mat0_color = 0;
       float  mat1_val = 0.4, mat1_alpha = 0;
       color  mat1_color = 1;
       float  mat2_val = 0.7, mat2_alpha = 5;
       color  mat2_color = 1;
       float  mat3_val = 1, mat3_alpha = 5;
       color  mat3_color = 0;
       float  mat4_val = 1, mat4_alpha = 0;
       color  mat4_color = 0;
       float  filter_it = 0;
       float  skip_below = 0.2;
       float  skip_step = 0.075;
       float  debug = 0;
    )
{
  point PP, dir, Nf, Ptemp, V = normalize(-I);
  float falloff, val, alpha, strength, here;
  color C;
  float done = 0, skipped_last = 0;
  float localthickness;
  color volcolor;
  color amb = Ka * ambient();
  point gradient, step;
  float v0, v1, v2, v3, v4, v5, v6, x, y, z, xp, yp;
  /* interpolation values */
  float x0, x1, y0, y1, z0, z1, xp0, xp1, yp0, yp1, val0, val1, nrpanels;


void interpolate_data(point P; output float val)
{
  extern float x0, y0, z0, x1, y1, z1, xp0, yp0, xp1, yp1, val0, val1, xpanels, ypanels, nrpanels;
  extern point Ptemp;

  Ptemp=(P);
  z0 = floor((zcomp(Ptemp))*nrpanels);
  yp0 = floor(z0/xpanels)/ypanels;
  xp0 = (z0 - xpanels*floor(z0/xpanels))/xpanels;
  x0 = xp0 + (xcomp(Ptemp))/xpanels;
  y0 = yp0 + (ycomp(Ptemp))/ypanels;
  val0 = texture(filename,x0,y0,x0,y0,x0,y0,x0,y0);
  z1 = ceil((zcomp(Ptemp))*nrpanels);
  if (z0 != z1)
  {
    yp1 = floor(z1/xpanels)/ypanels;
    xp1 = (z1 - xpanels*floor(z1/xpanels))/xpanels;
    x1 = xp1 + (xcomp(Ptemp))/xpanels;
    y1 = yp1 + (ycomp(Ptemp))/ypanels;
    val1 = texture(filename,x1,y1,x1,y1,x1,y1,x1,y1);
    val = linterp(zcomp(Ptemp), z0/nrpanels, z1/nrpanels, val0, val1);
  }
  else val = val0;
}

void interpolate_data_sameslice(point P; output float val)
{
  extern float x0, y0, z0, x1, y1, z1, xp0, yp0, xp1, yp1, val0, val1, xpanels, ypanels, nrpanels;
  extern point Ptemp;
  Ptemp=(P);
  x0 = xp0 + (xcomp(Ptemp))/xpanels;
  y0 = yp0 + (ycomp(Ptemp))/ypanels;
  val0 = texture(filename,x0,y0,x0,y0,x0,y0,x0,y0);
  if (z0 != z1)
  {
    x1 = xp1 + (xcomp(Ptemp))/xpanels;
    y1 = yp1 + (ycomp(Ptemp))/ypanels;
    val1 = texture(filename,x1,y1,x1,y1,x1,y1,x1,y1);
    val = linterp(zcomp(Ptemp), z0/nrpanels, z1/nrpanels, val0, val1);
  }
  else val = val0;
}

  /* Transform the current sampling point and ray direction to shader space */
  PP = transform ("shader", P) + 0.5;
  dir = normalize (vtransform ("shader", I));
  step = stepsize * dir;

  Oi = 0;  Ci = 0;

  nrpanels = xpanels*ypanels;
/*
  DEB  printf ("Starting levoy: stepsize = %f\n", stepsize);
  */

  while (done == 0) {

      /* Request color, alpha, gradient from volume */
/*      lookup_data(PP,here); */
      interpolate_data(PP,here);
      if (here < skip_below) {
          PP += skip_step * dir;
          skipped_last = 1;
          /* We're done if we're outside of the unit cube */
          x = xcomp(PP);  y = ycomp(PP);  z = zcomp(PP);
          if (x < 0 || x > 1 || y < 0 || y > 1 || z < 0 || z > 1)
              done = 1;
          continue;
        }
      else if (skipped_last > 0) {
          PP -= dir * (skip_step - stepsize);
/*          lookup_data (PP,here); */
          interpolate_data (PP,here);
        }
      skipped_last = 0;

/*      calc_gradient (PP, filename, gradient, gradientstepsize, here); */
      calc_interpolate_gradient (PP, filename, gradient, gradientstepsize, here);
      strength = length (gradient);
      val = here;
/*      access_volume (PP,val); */

      if (isovalue >= 0.0) {   /* we want an isovalue surface */
          localthickness = strength * isovaluethickness;
          if ((val - localthickness <= isovalue) &&
              (isovalue <= val + localthickness))
            {
              if (strength > 0.001) {
                alpha = isovaluealpha * (1 - abs ((isovalue - val) / localthickness));
                volcolor = 1;
              }
            }
          else {
              alpha = 0;
            }
        }
      else {     /* We want a region boundary surface */
          if (val >= mat0_val && val < mat1_val) {
              alpha = strength * linterp (val, mat0_val, mat1_val,
                                          mat0_alpha, mat1_alpha);
              volcolor = linterp (val, mat0_val, mat1_val,
                                  mat0_color, mat1_color);
            }
          else if (val >= mat1_val && val < mat2_val) {
              alpha = strength * linterp (val, mat1_val, mat2_val,
                                          mat1_alpha, mat2_alpha);
              volcolor = linterp (val, mat1_val, mat2_val,
                                  mat1_color, mat2_color);
            }
          else if (val >= mat2_val && val < mat3_val) {
              alpha = strength * linterp (val, mat2_val, mat3_val,
                                          mat2_alpha, mat3_alpha);
              volcolor = linterp (val, mat2_val, mat3_val,
                                  mat2_color, mat3_color);
            }
          else if (val >= mat3_val && val < mat4_val) {
              alpha = strength * linterp (val, mat3_val, mat4_val,
                                          mat3_alpha, mat4_alpha);
              volcolor = linterp (val, mat3_val, mat4_val,
                                  mat3_color, mat4_color);
            }
          else alpha = 0;
        }

      if (alpha > 0.001  &&  strength > 0.0001) {
          /* Transform gradient to world space and normalize for
           * surface normal
           */
          Nf = normalize (faceforward (vtransform ("shader", "current", gradient), I));

          C = 0;
          illuminance (PP, Nf, PI/2) {
              C += Cl * (Nf . normalize(L));
            }
          C = volcolor * (amb + Kd * C);

          /* Scale alpha by 1/stepsize */
          alpha *= stepsize;
          if (alpha > 1)
              alpha = 1;

          /* composite current step behind previous steps */
          Ci += (1-Oi) * C * alpha;
          Oi = 1 - (1-Oi)*(1-alpha);
          if (comp(Oi,1)+comp(Oi,2)+comp(Oi,3) > 3)
              done = 1;
        }
/*
DEB      printf ("    Ci = %p, Oi = %p\n", Ci, Oi);
*/

      /* Increment sampling point through the volume */
      PP += step;

      /* We're done if we're outside of the unit cube */
      x = xcomp(PP);  y = ycomp(PP);  z = zcomp(PP);
      if (x < 0 || x > 1 || y < 0 || y > 1 || z < 0 || z > 1)
        {
          done = 1;
        }
    }
}