/*
 * lafortune.sl -- Implement the BRDF representation of Lafortune et al.
 *
 * DESCRIPTION:
 *	Implement the BRDF representation described in "Non-Linear
 *	Approximation of Reflectance Functions", by Eric P. F. Lafortune,
 *	Sing-Choong Foo, Kenneth E. Torrance, and Donald P. Greenberg, in
 *	SIGGRAPH 97 Proceedings.
 *
 * PARAMETERS:
 *   coeff:	Coefficients for each lobe. There will be
 *		3 coefficients per lobe per wavelength.
 *		for each lobe of Nlobes:
 *		  for each channel of N_WAVES:
 *		    cxy, cz, n
 *		where cxy, cz are the directional and scale components
 *		n is the exponent for the cosine
 *
 *   Cs:	Diffuse color of surface. This provides the diffuse
 *		term that is an integral part of the model. This
 *		will be the diffuse component from the Lafortune
 *		model multiplied by pi, as the RenderMan definition
 *		of reflectance is different.
 * colormatrix:	3x3 color transformation matrix, if needed.
 *		This is here because measured data may be in the
 *		RGB space of some sensor (e.g. a digital camera)
 *		and may need correction for display. Since the three
 *		channels are represented independently (i.e. each
 *		channel will have different directionality), this
 *		correction must be done on a per-pixel basis.
 *		Since the correction may be different for each
 *		shader, we must do it here. Default is the identity
 *		matrix.
 *
 * ANTIALIASING: should antialias itself fairly well
 *
 * AUTHOR: written by Stephen H. Westin, Cornell University Program of
 *	Computer Graphics
 *
 * HISTORY:
 *
 * initial version 22 March 1999 S. H. Westin
 *
 * 22 March 1999 S. H. Westin
 *		Flipped sign on Cx to agree with description in SIGGRAPH paper
 * 25 October 2000 S. H. Westin
 *		Multiply lobes by incident cosine.
 * 27 October 2000 S. H. Westin
 *		Multiply lobes by incident cosine. Before, I was trying to
 *		use diffuse(), which doesn't work.
 *
 */

/* Number of coefficients per lobe: 3 for an isotropic surface */
#define LOBESIZE 3
/* Number of wavelengths (wired in) */
#define N_WAVES 3
/* Number of lobes (wired in) */
#define Nlobes 3

#define COEFFLEN (LOBESIZE*N_WAVES*Nlobes)

/* Default BRDF: flat blue paint */
/* Use "Color  [ .3094 .39667 .70837 ]" for consistency with default */

surface
lafortune ( uniform float coeff[27] =	 /* Change 27 to equal COEFFLEN */
	    { -0.870567, 0.803624, 21.820103,
	      -0.857255, 0.774290, 18.597755,
	      -0.670982, 0.586674, 7.472717,

	      0.451218, 0.023123, 2.774499,
	      0.406681, 0.017625, 2.581499,
	      0.477976, 0.227295, 3.677653,

	      1.031545, 0.706734, 66.899060,
	      1.029426, 0.696530, 63.767912,
	      1.026588, 0.687715, 57.489181 };
	    uniform float colormatrix[9] =
	    { 1, 0, 0, 0, 1, 0, 0, 0, 1 };
            uniform string texturename = ""; )
{

  varying vector local_z;
  varying vector local_x, local_y;	 /* Unit vector in "u" and "v" directions */
  varying vector V;                      /* Normalized eye vector */
  varying vector Ln;                     /* Normalized vector to light */
  varying float x, z, f;	 /* subterms */
  varying float fr, fg, fb;	 /* RGB components of the non-Lambertian term */
  uniform float basepointer, j;	 /* loop counters */

  /* Get unit vector in "u" parameter direction */
  local_x = normalize ( dPdu );

  local_z = faceforward (normalize(N), I);
  V = normalize (I);

  /* Get a local coordinate system. */
  local_y = local_z^local_x;
  
  /* The first term is the diffuse component. This should be the
     diffuse component in the Lafortune model multiplied by pi. */

  if ( texturename == "" )
    Ci = (diffuse(local_z) + ambient()) * Cs;
  else
    Ci = (diffuse(local_z) + ambient()) * color texture (texturename);

  illuminance ( P, local_z, 1.57079632679489661923 /* Hemisphere */ ) {
    Ln = normalize ( L );
    /* Compute the terms
       x = x_in * x_view  +  y_in * y_view
       z = z_in * z_view
    */
    x = local_x.V * local_x.Ln + local_y.V * local_y.Ln;
    z = - ( local_z.V * local_z.Ln );

    /* Coefficient structure:
       for each lobe of Nlobes:
          for each channel of N_WAVES:
            cxy, cz, n
       where cxy, cz are the directional and scale components
             n is the exponent for the cosine
    */

    for ( basepointer=0;
  	  basepointer<COEFFLEN;
  	  basepointer += LOBESIZE * N_WAVES ) {
      uniform float rexponent = coeff[basepointer+2],
	gexponent = coeff[basepointer+LOBESIZE+2],
	bexponent = coeff[basepointer+2*LOBESIZE+2];
      fr = fg = fb = 0.0;
      f = -x * coeff[basepointer] + z * coeff[basepointer+1];
      if ( f > 0.001*rexponent )
  	fr = pow ( f, rexponent ) ;
      f = -x * coeff[basepointer+LOBESIZE]
          + z * coeff[basepointer+LOBESIZE+1];
      if ( f > 0.001*gexponent )
	fg = pow ( f, gexponent ) ;
      f = -x * coeff[basepointer+2*LOBESIZE]
          + z * coeff[basepointer+2*LOBESIZE+1];
      if ( f > 0.001*bexponent )
	fb = pow ( f, bexponent ) ;
      Ci += ( Cl * color "rgb" ( fr, fg, fb ) ) *  (local_z.Ln);
    }
  }

  /* Color correction from camera space */
  /* use dot product with rows of matrix */
  fr = colormatrix[0]*comp(Ci,0)
    + colormatrix[1]*comp(Ci,1) + colormatrix[2]*comp(Ci,2);
  fg = colormatrix[3]*comp(Ci,0)
    + colormatrix[4]*comp(Ci,1) + colormatrix[5]*comp(Ci,2);
  fb = colormatrix[6]*comp(Ci,0)
    + colormatrix[7]*comp(Ci,1) + colormatrix[8]*comp(Ci,2);
  Ci = color "rgb" ( fr, fg, fb );
  Ci *= Os;
}