/*
 * Ashkihmin-Shirley BRDF
 * Peter Stuart
 * $Id: ash-shir.sl,v 1.9 2003/05/05 05:22:41 ophi Exp $
 * description: Ashikmin and Shirley BRDF lighting model as described
 * in "An Anisotropic Phong Light Reflection Model." A local illumination
 * version and an environment map version are both included. The environment
 * map version uses the importance sampling as described in the paper to reduce
 * variance.
 * 
 * ---- local illumination surface parameters ----------
 * N: normal direction
 * V: view direction
 * xdir: direction of anisotropy - can be primary or secondary
 * rd: diffuse reflectance
 * rs: specular reflectance
 * nu: specular exponent along xdir
 * nv: specular exponent perpendicular to xdir
 *
 * ---- environment illumination surface parameters ----------
 * N: normal direction
 * V: view direction
 * xdir: direction of anisotropy - can be primary or secondary
 * rd: diffuse reflectance
 * rs: specular reflectance
 * nu: specular exponent along xdir
 * nv: specular exponent perpendicular to xdir
 * samples: number of samples to take of the environment map
 * envmap: name of environment map
*/

/* ---- fast fresnel --------
 * H: vector halfway between V and L ((V+L) / 2)
 * V: view vector (-I/|I|)
 * f0: reflectance at normal incidence
 *
 * Schlick's approximation to the fresnel function
 */
float fastFresnel(vector H, V; float f0;)
{
  return f0 + (1-f0)*pow(1 - (H.V), 5);
}

/******** local illumination model **********/

color LocIllumAshShir (normal N; vector k1; vector xdir;
		       float rd, rs; float nu, nv;)
{
  float sqr(float x) { return x*x; }
  
  vector u_dir = xdir;
  vector v_dir = N ^ xdir;

  vector k2;
  vector h;
  float rho_d, rho_s;

  /* calculations that can be done one time ... */
  color diff = color(0,0,0), spec = color(0,0,0);
  float diffcont = ((28*rd) / (23)) * (1-rs)*(1-pow(1-(N.k1)/2, 5));
  float speccont = (sqrt((nu + 1)*(nv + 1)) / (8));

  color C = 0;
  extern point P;
  extern color Cs;
  illuminance(P, N, PI/2) {
    extern vector L;
    extern color Cl;

    k2 = normalize(L);
    h = normalize(k1 + k2);

    rho_d = diffcont * (1-pow(1-(N.k1)/2, 5));
    rho_s = speccont * (pow(N.h, (nu*sqr(h.u_dir) + nv*sqr(h.v_dir)) /
			    (1-sqr(h.N)))) / ((h.k1)*max(N.k1, N.k2)) *
      fastFresnel(h, k1, rs);

    C += (rho_d*Cs + rho_s) * Cl * max(0, (k2.N));
  }

  return C;
}

/******* "global" illumination model ********/

/* ---- function to sample environment ----
 * P: location of sample
 * R: reflection vector
 * envname: name of environment map - use raytracing if ""
 */
color SampleEnvironment (point P; vector R; uniform string envname;)
{
  color C = color(0,0,0);
  if (envname != "") {
    C = color environment(envname, R);
  } else {
    C = trace(P, R);
  }
  return C;
}

/* ---- helper function -----
 * nu, nv: anisotropic specular exponents
 * e1: random variable [0,1]
 */
float calcPhi(float nu, nv; float e1)
{
  return atan(sqrt((nu+1)/(nv+1)) * tan(PI*e1/2));
}

color EnvIllumAshShir (normal N; vector V; vector xdir;
		       float rd, rs; float nu, nv;
		       uniform float samples;
		       uniform string envmap;)
{
  uniform float i;
  color C = color(0,0,0);

  for (i=0; i<samples; i+=1) {
    float e1=random(), e2=random();
    float ph;

    /* Split up the hemisphere into quarters and remap e1 to one quarter.
       This helps to stratify the samples, hopefully decreasing variance... */
    /* FIXME this can be done better... */
    if (e1 >= 0 && e1 < 0.25)
      ph = calcPhi(nu, nv, 1 - 4*(0.25 - e1));
    else if (e1 >= 0.25 && e1 < 0.5)
      ph = calcPhi(nu, nv, 1 - 4*(0.5 - e1)) + PI/2;
    else if (e1 >= 0.5 && e1 < 0.75)
      ph = calcPhi(nu, nv, 1 - 4*(0.75 - e1)) + PI;
    else
      ph = calcPhi(nu, nv, 1 - 4*(1 - e1)) - PI/2;
    
    float cosph = cos(ph);
    float sinph = sin(ph);
    float th = acos(pow(1-e2, 1/(nu*cosph*cosph+nv*sinph*sinph+1)));

    if (th <= PI/2) { /* we don't want thetas below hemisphere */
    
      /* transform H (halfway vector) into space where frame is
	 <xdir, ydir, N> */
      vector ydir = N ^ xdir;
      vector H = cosph*sin(th)*xdir + sinph*sin(th)*ydir + cos(th)*N;
    
      extern point P;
      C += rs * SampleEnvironment(P, 2*H - V, envmap) * fastFresnel(H,V,rs);

      /* now generate the diffuse part using cos-weighted distribution. */
      e1 = random();
      e2 = random();
      H = cos(2*PI*e1)*sqrt(1-e2)*xdir +
	sin(2*PI*e1)*sqrt(1-e2)*ydir +
	sqrt(e2)*N;
      
      C += rd * SampleEnvironment(P, H, envmap);

      /* FIXME - instead of sampling over the hemisphere, get the cos-weighted
	 average value, which can be done with a cos-weighted filtered call to
	 environment. The problem is that such a filter doesn't seem to
	 exist :( */
/*     C += rd * color environment(envmap, xdir, ydir, -xdir, -ydir, */
/* 				"filter", "gaussian"); */
    }
  }

  if (samples > 0) C /= samples;
  return C;
}