// testint.C
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <iostream.h>
#include "nied31.h"
#include "vector.H"
#include "matrix.H"

namespace xrml {

  static unsigned idx;

  Vec3 uniformdir(double xi_1, double xi_2)
  {
    double phi, cos_phi, sin_phi, cos_theta, sin_theta;

    phi = 2. * M_PI * xi_1;
    cos_phi = cos(phi);
    sin_phi = sin(phi);
    cos_theta = 1.0-xi_2;
    sin_theta = sqrt(1.-cos_theta*cos_theta);

    return Vec3( cos_phi * sin_theta,
		 sin_phi * sin_theta,
		 cos_theta );
  }

  static Vec3 IdealReflectedDirection(const Vec3& inDir)
  {
    return Vec3(-inDir.x, -inDir.y, inDir.z);
  }

  class testfunc {
  public:
    virtual double f(const Vec3& inDir, const Vec3& outDir) { return 0.; }
    virtual double integral(const Vec3& inDir) { return 0.; }
  };

  class testf1: public testfunc {
  public:
    // constant function
    double f(const Vec3& inDir, const Vec3& outDir)
    {
      return outDir.z > 0.? 1. : 0.;
    }

    double integral(const Vec3& inDir)
    {
      return 2. * M_PI;
    }
  };

  class testf2: public testfunc {
  public:
    // cos(outDir, normal)
    double f(const Vec3& inDir, const Vec3& outDir)
    {
      return outDir.z > 0. ? outDir.z : 0.;
    }

    double integral(const Vec3& inDir)
    {
      return M_PI;
    }
  };

  class testf3: public testfunc {
  public:
    double sharpness;

    // Phong lobe with sharpness given below
    testf3(double s = 10.)
    {
      sharpness = s;
    }

    double f(const Vec3& inDir, const Vec3& outDir)
    {
      double cosa = IdealReflectedDirection(inDir) & outDir;
      return (cosa <= 0.) ? 0. /* angle >= 90deg */ : pow(cosa, sharpness);
    }

    double integral(const Vec3& inDir)
    {
      return (2 * M_PI) / (sharpness + 1);
    }
  };

  class bsdf {
  public:
    virtual Vec3 sample(const Vec3& inDir, double xi_1, double xi_2) { abort(); return Vec3(0,0,0); }
    virtual float eval(const Vec3& inDir, const Vec3& outDir, float *pdf) { abort(); return 0.; }
    virtual float albedo(const Vec3& inDir) {abort(); return 0.; }
  };

  class diffuse: public bsdf {
  public:
    Vec3 sample(const Vec3& inDir, double xi1, double xi2)
    {
      double phi = 2. * M_PI * xi1;
      double cos_phi = cos(phi);
      double sin_phi = sin(phi);
      double cos_theta = sqrt(xi2);
      double sin_theta = sqrt(1.-xi2);

      // generate outDir on the same side of the surface
      // as inDir.
      return Vec3(cos_phi * sin_theta,
		  sin_phi * sin_theta,
		  cos_theta * (inDir.z >= 0. ? 1. : -1.));
    }

    float eval(const Vec3& inDir, const Vec3& outDir, float *pdf)
    {
      if (pdf) *pdf = 0.;
      if ((inDir.z >= 0.) != (outDir.z >= 0))
	return 0.;    // refracted ray. Also pdf is 0.

      if (pdf) *pdf = outDir.z / M_PI;
      return 1. / M_PI;
    }

    float albedo(const Vec3& inDir)
    {
      return 1.;
    }
  };

  class phong: public bsdf {
  public:
    double sharpness;

    phong(double sharp = 10.)
    {
      sharpness = sharp;
    }

    Vec3 sample(const Vec3& inDir, double xi1, double xi2)
    {
      // rotation transforming the Z axis into
      // the ideal reflected direction.
      // The axis X is chosen to lay in the surface
      // tangent plane.
      Vec3 Z = IdealReflectedDirection(inDir);
      double zz = sqrt(1 - Z.z*Z.z);
      Vec3 X = (zz < EPSILON)	// normal incidence
	? Vec3(1.,0.,0.)
	: Vec3(Z.y/zz, -Z.x/zz, 0.);
      Vec3 Y = Z ^ X;
      Mat3 rot(X,Y,Z);

      double phi = 2. * M_PI * xi1;
      double cos_phi = cos(phi);
      double sin_phi = sin(phi);
      double cos_theta = pow(xi2, 1./(sharpness+1.));
      double sin_theta = sqrt(1. - cos_theta*cos_theta);
      return Vec3(cos_phi * sin_theta,
		  sin_phi * sin_theta,
		  cos_theta) * rot;
    }

    float eval(const Vec3& inDir, const Vec3& outDir, float *pdf)
    {
      if (pdf) *pdf = 0.;

      double cosa = IdealReflectedDirection(inDir) & outDir;
      if (cosa <= 0.)	// angle >= 90 degrees
	// the sampling algorithm above cannot yield such
	// directions, so also the pdf is 0.
	return 0.;

      double val = pow(cosa, sharpness);
      if (pdf) *pdf = val * (sharpness+1.) / (2. * M_PI);

      if ((outDir.z > 0.) != (inDir.z > 0.))
	// refracted direction
	// pdf is not 0: the sampling method above can yield rays 
	// leaving on the wrong side of the surface.
	return 0.;

      return val * (sharpness+2.) / (2. * M_PI);
    }

    // returns approximate albedo of modified Phong BRDF: we assume that
    // the cosine of the angle between outgoing direction and surface normal
    // remains constant over the Phong lobe. This is a good approximation for
    // sharp distributions.
    float albedo(const Vec3& inDir)
    {
      return fabs(inDir.z);
    }
  };

  class ward: public bsdf {
    // G. Ward, "Modeling and measuring anisotropic reflectance", SIGGRAPH'92
  public:
    double standardDeviationX, standardDeviationY;

    ward(double stddevX = 0.1, double stddevY = 0.1)
    {
      standardDeviationX = stddevX;
      standardDeviationY = stddevY;
    }

    // Samples outgoing direction for given incident directon and random
    // numbers in the range [0,1].
    // pdf is proportional to Wards eliptic Gaussian
    Vec3 sample(const Vec3& inDir, double xi_1, double xi_2)
    {
      double r = (xi_1 > EPSILON*EPSILON) ? sqrt(-log(xi_1)) : HUGE;
      xi_2 *= 2.*M_PI;
      double x = standardDeviationX * r * cos(xi_2);
      double y = standardDeviationY * r * sin(xi_2);
      double z = 1. / sqrt(x*x + y*y + 1.);	// choose z so h will be normalised
      x *= z; y *= z;
      Vec3 h = Vec3(x, y, z);			// halfway vector (microfacet normal)
      return h * (2. * (h & inDir)) - inDir;	// mirror inDir w.r.t. h
    }

    // evaluates Wards SIGGRAPH'92 scattering funtion for given incident and
    // outgoing direction. If 'pdf' is not nil, also computes the probability that
    // the routine sample() above would pick outDir.
    float eval(const Vec3& inDir, const Vec3& outDir, float *pdf)
    {
      if (pdf) *pdf = 0.;				// initialise

      Vec3 h = (Vec3)inDir + (Vec3)outDir;	// halfway vector (unnormalised).
      double xfac = (h.x / (h.z * standardDeviationX));	// normalisation of h doesn't matter.
      double yfac = (h.y / (h.z * standardDeviationY));
      double val = exp( - (xfac*xfac + yfac*yfac) );	// Wards eliptic Gaussian

      // For normalised h, one division more, by the square length of the unnormalised
      // halfway vector, would be necessary.
      double N = M_PI * standardDeviationX * standardDeviationY;
      if (pdf) *pdf = val * (h & outDir) / (h.z*h.z*h.z) / N;

      if ((inDir.z >= 0) != (outDir.z >= 0.))
	// outDir and inDir are on different sides of the surface
	// pdf is not 0: the sampling method above can yield rays leaving 
	// on the wrong side of the surface!
	return 0.;

      double normalisation = sqrt( inDir.z * outDir.z ) * 4 * N;
      if (normalisation < EPSILON)
	// vector lies in plane (inDir.z == 0 || outDir.z == 0) or illegal parameter values (standardDeviationX == 0 || standardDeviationY == 0)
	return 0.;

      return val / normalisation;
    }

    // Computes the integral over all outgoing direction of Wards scattering
    // function times the outgoing cosine. This is the fraction of energy incident
    // from 'inDir' that gets reflected.
    // This approximation for the albedo is good for small standardDeviationX, 
    // standardDeviationY and for normal incidence. For grazing angles, the error can be 
    // large because part of the lobe will be under the surface.
    float albedo(const Vec3& inDir)
    {
      return fabs(inDir.z);
    }
  };

  class schlick: public bsdf {
    // C. Schlick, "An inexpensive BRDF model for physically based rendering",
    // Computer Graphics Forum Vol 13, nr 3 (Proceedings of EuroGraphics 1994), p 233
    // Modified so that the albedo will not depend on roughness and with proper
    // pdf.
  public:
    double roughness;
    double isotropy;

  protected:
    bool parmschecked;

    inline void checkparms(void)
    {
      if (roughness < EPSILON || isotropy < EPSILON) {
	cerr << "Invalid parameters\n";
      }
      parmschecked = true;
    }

  public:
    schlick(double roughn =1., double isot =1.)
    {
      roughness = roughn;
      isotropy = isot;
      parmschecked = false;
    }

    Vec3 sample(const Vec3& inDir, double xi_1, double xi_2)
    {
      if (!parmschecked) checkparms();

      // zenith angle
      double c2 = xi_1 / (roughness + (1.-roughness) * xi_1);
      double s = sqrt(1. - c2);

      // determine in which quadrant the azimuthal angle must be
      if (xi_2 > 1.-EPSILON) xi_2 = 1.-EPSILON;
      xi_2 *= 4.;
      double fxi2 = floor(xi_2);
      int quadrant = (int)fxi2;

      // make random nr in the rand [0,1) again
      // mirrored in quadrant 1 and 3, so result will change
      // continuously as a function of xi2.
      xi_2 = xi_2 - fxi2;
      if (quadrant&1) xi_2 = 1. - xi_2;

      // azimuth angle in first quadrant
      double p2 = isotropy * isotropy;
      double b2 = xi_2 * xi_2;
      double psi = 0.5*M_PI * sqrt( p2 * b2 / (1. - (1.-p2) * b2) );

      Vec3 h(cos(psi)*s, sin(psi)*s, sqrt(c2));

      switch (quadrant) {
      case 0: break;
      case 1: h.x = -h.x; break;
      case 2: h.x = -h.x; h.y = -h.y; break;
      case 3: h.y = -h.y; break;
      }
    
      return h * (2. * (h & inDir)) - inDir;	// mirror inDir w.r.t. h
    }

    float eval(const Vec3& inDir, const Vec3& outDir, float *pdf)
    {
      if (!parmschecked) checkparms();

      // halfway vector (microfacet normal)
      Vec3 h = (Vec3)inDir + (Vec3)outDir;
      double H = (h & h);
      if (H < EPSILON) {
	// inDir == -outDir, rarely results from the sampling algorithm above.
	if (pdf) *pdf = 1.;	// avoid divisions by zero.
	return 0.;
      }
      H = sqrt(H);
      h *= (1./H);

      // zenith angle dependence
      double t2 = h.z*h.z; 					// t^2
      double Z = (1-t2) + roughness * t2;
      Z = roughness / (Z*Z);

      // azimuth angle dependence
      double w2 = (t2 > 1.-EPSILON) ? 0. : h.x*h.x / (1. - t2) ; 	// w^2
      double A = sqrt( isotropy / (isotropy*isotropy * (1-w2) + w2) );

      if (pdf) {
	// The pdf corresponding to Schlicks sampling formula (formula 29-right
	// on page C-244) is NOT A(w) dw / sqrt(1-w^2) but something completely 
	// different !!! It is not even particularly suited to integrate Schlicks
	// own BRDF when anisotropic !!!
	double x = (2./M_PI) * ((w2<EPSILON*EPSILON) ? (M_PI/2.) : (w2>=1.) ? 0. : acos(sqrt(w2)));
	double p = isotropy;
	double p2 = p*p;
	double A1 = 1 + (1-p2)/p2 * x*x;
	*pdf = fabs(((Vec3)h & (Vec3)outDir) * h.z) / (M_PI * H*H) * Z / (p * A1 * sqrt(A1));
      }

      double vi = inDir.z;
      double vo = outDir.z;

      if (fabs(vi) < EPSILON || fabs(vo) < EPSILON || (vi>0.) != (vo>0.))
	// incoming or outgoing direction in surface tangent plane or
	// they are directions on different sides of the surface
	// pdf is not zero: the sampling algorithm above can yield such
	// directions.
	return 0.;
    
      double Gi = vi / (roughness - roughness * vi + vi);
      double Go = vo / (roughness - roughness * vo + vo);

      return (Gi * Go * A * Z + (1. - Gi * Go)) / (4.*M_PI * vi * vo);
    }

    float albedo(const Vec3& inDir)
    {
      // The value below is an approximation only.
      return fabs(inDir.z);
    }
  };


  class ashikmin: public bsdf {

  protected:
    float nu, nv;

  public:
    ashikmin(double _nu =1., double _nv =1.)
    {
      nu = _nu;
      nv = _nv;
    }
  
    float albedo(const Vec3& inDir)
    {
      return fabs(inDir.z);
    }

    Vec3 sample(const Vec3& inDir, float xi1, float xi2)
    {
      double theta = atan( sqrt( (nu + 1) / (nv + 1) ) * tan( M_PI * xi1 / 2 ) );
      double sinTheta = sin(theta);
      double cosTheta = cos(theta);

      double cosDelta = pow( (1 - xi2), 1. / (nu * cosTheta * cosTheta + nv * sinTheta * sinTheta + 1) );
      double sinDelta = sqrt(1 - cosDelta * cosDelta);

      Vec3 h = Vec3( sinTheta * cosDelta, sinTheta * sinDelta, cosTheta );
      
      return Vec3(-inDir + h * 2 * (inDir & h));
    }

    float eval(const Vec3& inDir,const Vec3& outDir, float *pdf)
    {
      Vec3 h = Vec3(inDir + outDir).normalize();
      double hu = h.x;
      double hv = h.y;
      double hn = h.z;
      double hk1 = h & inDir;
      double hk2 = h & outDir;
      double nk1 = inDir.z;
      double nk2 = outDir.z;

      if(hn < EPSILON) return(1.);
      if(hk2 < EPSILON) return(1.);
      if(hk1 < EPSILON) return(1.);

      double c = sqrt( (nu + 1) * (nv + 1) ) / (8 * M_PI);
      double max = (nk1 > nk2) ? nk1 : nk2;
      double p = pow(hn, (nu*hu*hu + nv*hv*hv) / (1 - hn*hn) );
      
      if(pdf) *pdf = (float)(c * p / hk1);
  
      return( (float)((c * p) / (hk2 * max)) );
    }
  };

  class neumann : public bsdf
  {
  protected:
    float sharpness;

  public:
    neumann(double s=0.)
    {
      sharpness = s;
    }

    static Vec3 Mirror(Vec3 dir)
    {
      return Vec3(-dir.x, -dir.y, dir.z);
    }

    float albedo(const Vec3& inDir)
    {
    return fabs(inDir.z);
    }

    Vec3 sample(const Vec3& inDir, float xi1, float xi2)
    {
      // uses the sample sampling as (and pdf) as in:
      //
      // "Using the Modified Phong Reflectance Model for Physically Based Rendering"
      // Eric P. Lafortune, Yves D. Willems
      // Technical Report CW 197, November 1994
      // (Department of Computing Science, K.U.Leuven)
      //
      float u = pow(xi1, 2. / (sharpness + 1));
      float v = sqrt(1 - u*u);
      float phi = 2 * M_PI * xi2;
      
      Vec3 dir = Vec3(v * cos(phi), v * sin(phi), u);
      
      return Mirror(dir);
    }

    float eval(const Vec3& inDir, const Vec3& outDir, float *pdf)
    {
      if (pdf) *pdf = 0.;

      double cosa = Mirror(inDir) & outDir;
      if (cosa <= 0.)	// angle >= 90 degrees
	// pdf is 0: such outgoing rays cannot result from the sampling method
	// above.
	return 0.;

      double val = pow(cosa, sharpness);
      double cos_max = (inDir.z > outDir.z) ? inDir.z : outDir.z;
      
      if (pdf) *pdf = val * (sharpness+1.) / (2. * M_PI);

      if ((outDir.z > 0.) != (inDir.z > 0.))
	// refracted direction
	// pdf is not 0: the sampling method above can yield rays
	// leaving the wrong side of the surface!
	return 0.;

      return val * (sharpness+2.) / (2. * M_PI * cos_max);
    }
  };

  static int nrsamples = 100000000;

  // computes albedo by means of uniform random sampling and compares the result
  // with what we think it should be.
  // Monte Carlo integration: generate random samples with certain pdf, compute
  // average ratio of function value (here: bsdf x cos) over sample pdf.
  double TestAlbedo(class bsdf* sampling, class bsdf *s, Vec3 inDir)
  {
    int n = 1;
    double sum = 0., ratio = 1.;
    while (n <= nrsamples) {
      unsigned *zeta = Nied31(idx++);
      Vec3 outDir = sampling->sample(inDir, (double)zeta[0]*RECIP, (double)zeta[1]*RECIP);
      float pdf;
      sampling->eval(inDir, outDir, &pdf);

      float bsdf = s->eval(inDir, outDir, 0);
      sum += bsdf * outDir.z / pdf;   	// outDir.z = cos(outDir, normal=Z)

      if ((n>0 && floor(log10((double)n)) != floor(log10((double)n-1))) ||
	  n%1000000 == 0)
	fprintf(stderr, "%8d dirs: albedo = %f (verhouding = %f)\n", n, 
		sum / (double)n ,
		ratio = (sum / (double)n) / s->albedo(inDir) );

      n++;
    }
    return ratio;
  }

  double TestSampling(class bsdf *s, Vec3 inDir, class testfunc *test)
  {
    int n = 1;
    double sum = 0., norm = 1.;
    while (n <= nrsamples) {
      unsigned *zeta = Nied31(idx++);
      Vec3 outDir = s->sample(inDir, (double)zeta[0]*RECIP, (double)zeta[1]*RECIP);
      float pdf;
      s->eval(inDir, outDir, &pdf);
      sum += test->f(inDir, outDir) / pdf;

      if ((n>0 && floor(log10((double)n)) != floor(log10((double)n-1))) ||
	  n%1000000 == 0)
	fprintf(stderr, "%8d dirs: integral = %f (verhouding = %f)\n", n, 
		sum / (double)n , 
		norm = (sum / (double)n) / test->integral(inDir) );

      n++;
    }
    return norm;
  }

  // evaluates bsdf and pdf for directoins returns by s->sample() and writes out the
  // result in gnuplot format
  void PlotBsdf(class bsdf *s, Vec3 inDir)
  {
    int divs;
    fprintf(stderr, "Geef fijnheid (geheel getal, bvb 10)\n");
    scanf("%d", &divs);

    int i, j;
    printf("# theta phi bsdf inDir = (%g,%g,%g)\n", inDir.x, inDir.y, inDir.z);
    for (i=0; i<=divs; i++) {
      for (j=0; j<=divs; j++) {
	Vec3 outDir = s->sample(inDir, (double)j/(double)divs, (double)i/(double)divs);
	float pdf, bsdf = s->eval(inDir, outDir, &pdf);
	printf("%f %f %f   %f %f %f     %f %f %f  %f %f %f\n",
	       outDir.x * bsdf, outDir.y * bsdf, outDir.z * bsdf,
	       outDir.x * pdf, outDir.y * pdf, outDir.z * pdf,
	       outDir.x * bsdf * outDir.z, outDir.y * bsdf * outDir.z, outDir.z * bsdf * outDir.z,
	       outDir.x * pdf / outDir.z, outDir.y * pdf / outDir.z, pdf);
      }
      printf("\n");
    }
  }

  Vec3 ChooseIndir(void)
  {
    Vec3 inDir;
    fprintf(stderr, "Geef inDir (x,y en z, z is volgens de normaal):\n");
    scanf("%f %f %f", &inDir.x, &inDir.y, &inDir.z);
    inDir.normalize();

    return inDir;
  }

  class testfunc* ChooseTestf1(void)
  {
    return new testf1;
  }

  class testfunc* ChooseTestf2(void)
  {
    return new testf2;
  }

  class testfunc* ChooseTestf3(void)
  {
    double sharpness;
    fprintf(stderr, "Enter sharpness:\n");
    scanf("%lf", &sharpness);
    return new testf3(sharpness);
  }

  class testfunc* ChooseTestfunc(void)
  {
    fprintf(stderr, "Welke testfunctie?\n");
    fprintf(stderr, "1 = constante functie\n");
    fprintf(stderr, "2 = cosinus tussen uitgaande hoek en normaal (= diffuus model)\n");
    fprintf(stderr, "3 = Phong lobe\n");
    int which;
    scanf("%d", &which);

    switch (which) {
    case 1: return ChooseTestf1(); 
    case 2: return ChooseTestf2(); 
    case 3: return ChooseTestf3(); 
    default:
      fprintf(stderr, "No such test function.\n");
      exit(0);
    }
    return 0;
  }

  class bsdf* ChooseDiffuse(void)
  {
    return new diffuse;
  }

  class bsdf* ChoosePhong(void)
  {
    double sharpness;
    fprintf(stderr, "Enter sharpness:\n");
    scanf("%lf", &sharpness);
    return new phong(sharpness);
  }

  class bsdf* ChooseWard(void)
  {
    fprintf(stderr, "Enter standard deviation X and Y:\n");
    double stddevX, stddevY;
    scanf("%lf %lf", &stddevX, &stddevY);
    return new ward(stddevX, stddevY);
  }

  class bsdf* ChooseSchlick(void)
  {
    fprintf(stderr, "Enter roughness and isotropy:\n");
    double roughness, isotropy;
    scanf("%lf %lf", &roughness, &isotropy);
    return new schlick(roughness, isotropy);
  }

  class bsdf* ChooseAshikmin(void)
  {
    fprintf(stderr, "Enter nu and nv:\n");
    double nu, nv;
    scanf("%lf %lf", &nu, &nv);
    return new ashikmin(nu, nv);
  }

  class bsdf* ChooseNeumann(void)
  {
    fprintf(stderr, "Enter sharpness:\n");
    double s;
    scanf("%lf", &s);
    return new neumann(s);
  }

  class bsdf *ChooseBsdf(char *question)
  {
    fprintf(stderr, "%s\n", question);
    fprintf(stderr, "1 = diffuse\n");
    fprintf(stderr, "2 = Phong\n");
    fprintf(stderr, "3 = Ward\n");
    fprintf(stderr, "4 = Schlick\n");
    fprintf(stderr, "5 = Ashikmin\n");
    fprintf(stderr, "6 = Neumann\n");

    int which;
    scanf("%d", &which);
    switch (which) {
    case 1: return ChooseDiffuse(); break;
    case 2: return ChoosePhong(); break;
    case 3: return ChooseWard(); break;
    case 4: return ChooseSchlick(); break;
    case 5: return ChooseAshikmin(); break;
    case 6: return ChooseNeumann(); break;
    default:
      fprintf(stderr, "Ongeldige keuze.\n"); exit(-1);
    }

    return 0;
  }

  void do_bla(void)
  {
  }

}

using namespace xrml;

int main(int argc, char **argv)
{
  char *progname = argv[0];
  fprintf(stderr, "progname = '%s'\n", progname);
  
  idx = 1;
  fprintf(stderr, "Seeding random nr generator with %d\n", idx);

  if (strcmp(argv[0], "testalbedo") == 0) {
    TestAlbedo(ChooseBsdf("Welke sampling methode"), ChooseBsdf("Welke BSDF"), ChooseIndir());
  }
  
  else if (strcmp(argv[0], "plotbsdf") == 0) {
    PlotBsdf(ChooseBsdf("Welke BSDF"), ChooseIndir());
  }
  
  else if (strcmp(argv[0], "testsampling") == 0) {
    TestSampling(ChooseBsdf("Welke BSDF"), ChooseIndir(), ChooseTestfunc());
  }
  
  else {
    do_bla();
  }
    
  return 0;
}