// ===================================================================
// $Id$
//
// cook_torrance.cpp
//    Implementation of Cook-Torrance BRDF.
//
// Class: CCookTorrance
//
// Initial coding by Polina Kondratieva (PK), January 2004.
// ===================================================================

// GOLEM headers
#include "config.h"
#include "cook_torrance.h"
#include <iostream>

#define	M_4PI 4.0 * M_PI

__BEGIN_BRDFLIB_SOURCE


/**********************************************************************
 *
 *             CCookTorrance methods
 *
 **********************************************************************/
// describes the surface to a given stream
//========================================================================
void
CCookTorrance::Describe(ostream &app, int ind) const
{
  indent(app, ind);
  app << GetName() << ":\n";

  if (_diffuse) {
    indent(app, ind+1);
    app << "  diffuse  = \t";
    _diffuse->Describe(app, ind); app << "\n";
  }

  if (_specular) {
    indent(app, ind+1);
    app << "  specular = \t";
    _specular->Describe(app, ind); app << "\n";

    if (_slope) {
      indent(app, ind+1);
      app << "  slope    = \t";
      _slope->Describe(app, ind); app << "\n";
    }
  }

  app << flush;

  return;
}

// Deep copy of the class entity
//========================================================================
CSurface *CCookTorrance::Dup() const
{
  CCookTorrance *ret = CreateNewObject();

  if (_diffuse)
    ret->SetDiffuse(_diffuse);
  if (_specular)
  {
     ret->_specular = _specular;

     if(_slope)
        ret->_slope = _slope;
  }

  ret->_supportedType = _supportedType;
  return ret;
}
//========================================================================

/*
  Albedo is the integral of BRDF multiplied by cos(Eye,N) over the hemisphere.
  Integral of Cook-Torrance BRDF cannot be evaluated analitically. However, it
  can be evaluated using Pade approximation. Unfortunately, it requres a lot
  of computations. So, for a now we leave this albedo function in this way.

  Possibly later on, we'll provide this integral evaluation in the way
  similar to the Arvo's procedures (Arvo'1995 "Applications of Irradiance
  Tensors to the Simulation of Non-Lambertian Phenomena"). */
//========================================================================
bool
CCookTorrance::GetAlbedo(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
  // Aldebo is a function of the angle bewtween the incoming direction
  // and normal, but in most cases only the maximum of this function
  // computed below is intresting
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

  // Diffuse part of albedo
  if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
    hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
    c += d;
    retval = true;
  }

  // Specular part of albedo
  if (((distType & _microfacetReflectionType()) && _specular) && _slope) {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
    c += d;
    retval = true;
  }

//#define _DEBUG_PHYSICAL_PLAUSIBILITY
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    cout << "Current parameter settings for Cook-Torrance model make it\n";
    cout << "physically unplausible\n";
    exit(0);
  }
#endif // _DEBUG_PHYSICAL_PLAUSIBILITY

  if (retval) {
    albedo = c;
    return true;
  }

  // unable to compute the value for some reason
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  if ((!_diffuse) && (!_specular)) {
    WARNING << "This is the black body - albedo cannot be computed" << endl;
  }
#endif // _DEBUG

  return false;
}

/*=============================================================================
   BRDF  =   Distribution * Fresnel(Light) * GeometricFactor
            ---------------------------------------------------------
                4 * Pi * lightCos * cosEye

  Distribution     = sum(i=0..N, wi*Distribution(ai))
                     (w0 + w1 + ... + wN = 1)

  Distribution(ai) = exp(-tan(H,N)^2/ai^2) / (ai^2 * cos(H,N)^4)

  (see "A Reflectance Model for Computer Graphics", Cook, Torrance, 1981)
=============================================================================*/
// Direct lighting computation
//=============================================================================
bool
CCookTorrance::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirLit,
			 EDistributionType   distType) const
{
  bool retval 	    = false;
  float d           = Magnitude(dirLit);
  CVector3D N       = hitInfo.GetNormal();
  CVector3D dirEye  = (hitInfo.GetRay())->GetDir();
  float lightCos    = DotProd(dirLit, N);
  float incEyeCos   = -(hitInfo.GetDotProd());

  // Provide 'zero' result if Light source or viewer is under the surface
  if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
  {
    if(lightCos < WTINY) lightCos       = WTINY;
    if(incEyeCos < WTINY) incEyeCos  = WTINY;

    if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
      Normalize(dirLit);

    // Diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
      CValues rho_d;
      hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
      result = rho_d * M_1_PI * lightCos;
      retval = true;
    }

   // Specular part
   //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) {
      int slopes_size  = _slope->GetSlopesSize();
      int weights_size = _slope->GetWeightsSize();

      // Compute specular part only if slope structure is defined
      if(slopes_size > 0 && (slopes_size == weights_size))
      {
        float *a2      = new float[slopes_size];
	float *weights = new float[slopes_size];

	// Get values for slopes and their weights
	for(int i = 0; i < slopes_size; i++)
	{
          float weight = 1.0;
	  CTexture *slope = _slope->GetSlope(i, weight);
	  assert(slope);
          CValues s_tmp(CValues::EE_Weigth3);
	  hitInfo.GetTextureValue(s_tmp, *slope, CMappingBase::EE_KSI);

	  a2[i]      = (s_tmp.Get(0)) * (s_tmp.Get(0));
          if(a2[i] < WTINY) a2[i] = WTINY;

	  weights[i] = weight;
	}

	CValues rho_s;
        bool use_fresnel = _specular->IsFresnel();

    // Specular texture for Cook-Torrance model must be defined
	// as Fresnel texture
	if(use_fresnel)
       hitInfo.GetTextureValue(rho_s, *_specular, CMappingBase::EE_KFM);
    else
	{
		WARNING_INFO("CSurfaceCookTorrance()::EvaluateBsdf() : Fresnel texture is not defined!");
		return false;
	}

	CVector3D  H  =  Normalize(dirLit - (hitInfo.GetRay()->GetDir()));
        float cos_HN  =  DotProd(H, N);
	float cos_HV  =  DotProd(H, -dirEye);
        assert(cos_HN >= 0);

        if(cos_HN < WTINY)  cos_HN = WTINY;
	float cos_HN4 = pow(cos_HN, 4.0f);

	float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
        if(tan2_HN < 0) tan2_HN = 0;

	// compute term Distribution
	//--------------------------------------------------------------
	float distribution = 0;
	for(int i = 0; i < slopes_size; i++)
	{
          float memb    = exp( -tan2_HN/a2[i] )/(a2[i] * cos_HN4);
	  distribution += weights[i] * memb;
	}

	/* k =    Distribution * lightCos      =   Distribution
	       ------------------------------     ----------------
	         4 * Pi * lightCos * cosEye       4 * Pi * cosEye */
	//--------------------------------------------------------------
	float k = distribution / (incEyeCos * M_4PI);

	// Compute geometrical attenuation factor (accounting for
	// shadowing and masking)
        //--------------------------------------------------------------
	float part1 = 2.0 * cos_HN * incEyeCos/cos_HV;
	float part2 = 2.0 * cos_HN * lightCos/cos_HV;
	float G = 2.0;
	if(part1 < part2)
	  G = ((part1 < 1.0) ? part1 : 1.0);
	else
	  G = ((part2 < 1.0) ? part2 : 1.0);

        // Add all together
	result += rho_s * k * G;
        retval = true;

	// delete temporary structures
	delete [] a2;
	delete [] weights;
      }
    }
    /* returned result = BRDF * lightCos
       But BRDF model contains lightCos in denominator => result
       should not be multiplied on it */
  }
  else
  {
    retval = false;
  }
  return retval;
}

// Diffuse part for stochastic sampling
//========================================================================
int
CCookTorrance::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  const CValues       &ks_local,
  CVector2D           &slopeXY,
  CVector2D           &slopeXY_sq,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N      = hitInfo.GetNormal();
  CVector3D dirEye = (hitInfo.GetRay())->GetDir();
  float incEyeCos  = hitInfo.GetDotProd();

  // If light source/viewer is under the surface
  if (incEyeCos > 0.0)
    N = -N;

  incEyeCos = fabs(incEyeCos);

  CVector3D U = ArbitraryNormal(N);    // U is perpendicular to N
  CVector3D V = CrossProd(N, U);       // orthonormal base with N and AnisAxis

  int num = 0;			       // number of rays created

  // Make some precomutations for additional parts (for vizialization)
  //----------------------------------------------------------------------
  int slopes_size  = _slope->GetSlopesSize();
  int weights_size = _slope->GetWeightsSize();
  float *a2        = new float[slopes_size];
  float *weights   = new float[slopes_size];
  float avg_w      = 0;

  if(slopes_size > 0 && (slopes_size == weights_size))
  {
    // Gather all the distributions specified for this function
    for(int i = 0; i < slopes_size; i++)
    {
      float weight = 1.0;
      CTexture *slope = _slope->GetSlope(i, weight);
      assert(slope);
      CValues s_tmp(CValues::EE_Weigth3);
      hitInfo.GetTextureValue(s_tmp, *slope, CMappingBase::EE_KSI);
       //slopeC[i]  = s_tmp;
      a2[i]      = (s_tmp.Get(0)) * (s_tmp.Get(0));
      if(a2[i] < WTINY) a2[i] = WTINY;

      weights[i] = weight;
      avg_w     += weight;
    }
    avg_w       /= slopes_size;
  }

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
   // TODO
    // For Halton sequence random numbers generation
    //if (num) {
      // increase the seed not the generate the same sample again
      //samplingInfo.IncSeed();
      // generate a new sample for this time
      //aux_rnd = samplingInfo.Sample(0);
    //}
     // allocate next contribution
    CContributionItem *contrib = carray.AllocateNext();
    // generate the sample in the cannonical space
    CVector3D dirN;
    // Generate a sample in 2D, range [0;1)x[0;1)
    CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

    // probability is set to value cos(theta)/M_PI
    contrib->prob = RND::hemisph_cos(dirN, vv);
    // transform computed sample dirN to real coordinates
    CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
    // set the reflected ray
    contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
    // set the flags concerning ray type
    contrib->SetFlags(CContributionItem::EE_ReflectedRay);
    // type of interaction
    contrib->interaction = EE_REFLECTED_DIFFUSE;
    // diffuseness of the bounce
    contrib->diffuseness = 1.0;
    // Set the value of BSDF type
    // compute = BRDF * cos(normal, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // This means that we should not multiply the weight by the cosine
    // of the angle between the incoming light direction and normal
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = kd_local;

    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
    if(slopes_size > 0 && (slopes_size == weights_size) && avg_w > 0) {
      if(ks_local.Luminance()/avg_w > 1e-06)
      {
        float lightCos = DotProd(dirRes, N);
        assert(lightCos >= -WTINY);

        if(incEyeCos < WTINY) incEyeCos = WTINY;
        if(lightCos  < WTINY) lightCos  = WTINY;

        CVector3D  H  =  Normalize(dirRes - (hitInfo.GetRay()->GetDir()));
        float cos_HN  =  DotProd(H, N);
        float cos_HV  =  DotProd(H, -dirEye);
        if(cos_HN < WTINY)  cos_HN = WTINY;

        float cos_HN4 =  pow(cos_HN, 4.0f);
	float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
        if(tan2_HN < 0) tan2_HN = 0;

	// compute distribution
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	float distribution = 0;
	for(int j = 0; j < slopes_size; j++)
	{
          float memb    = exp( -tan2_HN/a2[j] )/(a2[j] * cos_HN4);
	  distribution += weights[i] * memb;
	}

	float k = distribution / (incEyeCos * 4.0 * M_PI);
 	k      *= 1.0/contrib->prob;

	// Compute geometrical attenuation factor (accounting for
	// shadowing and masking)
        //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	float part1 = 2.0 * cos_HN * incEyeCos/cos_HV;
	float part2 = 2.0 * cos_HN * lightCos/cos_HV;
	float G = 2.0;
	if(part1 < part2)
	   G = ((part1 < 1.0) ? part1 : 1.0);
	else
	   G = ((part2 < 1.0) ? part2 : 1.0);

	// add specular part to computed weight
	contrib->weigth += ks_local * k * G;
      }
    }
    //------------------------------------------------------------------------

#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

  } // for i
  // delete temporary arrays
  delete [] a2;
  delete [] weights;

  return maxRays;
}
//========================================================================

/* Specular part for stochastic sampling which is implemented in a similar
   to importance sampling scheme for Ward model way. (see paper Kelemen'2001)

   Distribution function contains several members (in other words, surface
   consists of several types of microfacets). In this way, to provide the
   importance sampling only one type of microfacet should be choosen.

   For vizualisation not only the diffuse part should be added to
   the computed sample weight, but the `rest` of the specular part represented
   by unchoosen distribution members. */
//========================================================================
int
CCookTorrance::SpawnStochasticSpecularRays(
  const CValues       &ks_local,
  CVector2D           &slopeXY,
  CVector2D           &slopeXY_sq,
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  int num = 0;		// the number of samples generated

  float incEyeCos  = hitInfo.GetDotProd();
  CVector3D N 	   = hitInfo.GetNormal();
  CVector3D dirEye = (hitInfo.GetRay())->GetDir();

  if (incEyeCos > 0.0)
    N = -N;

  // Choose one of the members of the distribution fucntion for
  // importance sampling scheme
  //----------------------------------------------------------------------
  float rnd_dist    = specular_rnd[0];
  int   choose_dist = 0;
  float chWeight    = 1.0;
  int   slopes_size = _slope->GetSlopesSize();
  int  weights_size = _slope->GetWeightsSize();
  if(slopes_size == 0 || (slopes_size != weights_size))
  {
	WARNING_INFO("For Cook-Torrance BRDF model non-empty slopes structure\nwith non-zero weights must be specified");
    return 0;
  }
  CTexture *chSlope = NULL;

  float alpha    = 0;
  float *a2      = new float[slopes_size];
  float *weights = new float[slopes_size];
  bool choosen   = false;

  // choose the distribution member according to its weight
  for(int i = 0; i < slopes_size; i++)
  {
    float weight = 1.0;
    CTexture *slope = _slope->GetSlope(i, weight);
    assert(slope);
    CValues s_tmp(CValues::EE_Weigth3);
    hitInfo.GetTextureValue(s_tmp, *slope, CMappingBase::EE_KSI);

    a2[i] = (s_tmp.Get(0)) * (s_tmp.Get(0));
    if(a2[i] < WTINY) a2[i] = WTINY;

    weights[i] = weight;

    if(!choosen)
    {
      alpha += weight;
      if(rnd_dist <= alpha)
      {
        choose_dist = i;
        chSlope  = slope;
        chWeight = weight;
        choosen = true;
      }
    }
  }
  //----------------------------------------------------------------------

  // Get value of the choosen distribution member
  CValues slopeC;
  hitInfo.GetTextureValue(slopeC, *chSlope, CMappingBase::EE_KSI);
  float slopeC_sq = (slopeC.Get(0)) * (slopeC.Get(0));

  // Make new basis with R as Z axis
  CVector3D U = ArbitraryNormal(N); // U is perpendicular to R - X-axis
  CVector3D V = CrossProd(N, U);    // V is perpendicular to R and V - Y-axis

  // TODO diffuseness of the bounce
  float diffuseness = slopeXY.xx;

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    // the ray is sampled in this direction
    CContributionItem *contrib = carray.AllocateNext();
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

    // set the flags concerning ray type
    contrib->SetFlags(CContributionItem::EE_ReflectedRay);

    float outCos = -1;
    CVector3D dirRes(0,0,1);
    // If slope is not very small
    //-------------------------------------------------------------------------
    if(slopeC_sq > 1e-06)
    {
      float cos_HI = -1;
      CVector3D hvector(0,0,1);
      float cos_HN = 0, cos2_HN = 0, sin_HN = 1, tan2_HN = 0;
      bool rejected = false;

      // rejection sampling
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      while(outCos < 0) {
        // TODO: Halton sequence
        //if (num) {
           // increase the seed not the generate the same sample again
           //samplingInfo.IncSeed();
           // generate a new sample for this time
           //aux_rnd = samplingInfo.Sample(0);
         //}
        CVector2D h;
        if( !rejected )
          h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
        else
          h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

        tan2_HN  = slopeC_sq * (-log(1.0-h.xx));
        assert(tan2_HN >= 0);
        cos2_HN  = (1.0/(1.0 + tan2_HN));
        cos_HN   = sqrt(cos2_HN);
        sin_HN   = sqrt(1.0 - cos2_HN);

        assert(sin_HN >= 0 && cos_HN >= 0);

        float phi_h   = 2.0 * M_PI * h.yy;

        // Outgoing direction (half-vector H) in local coordinate system
        CVector3D localDir(sin_HN * cos(phi_h), sin_HN * sin(phi_h), cos_HN);

        // Change the basis from the VUN coordinate system (local coordinate
        // system where normal R is Z-axis) to the World coordinate system
        hvector = localDir.x * U + localDir.y * V + localDir.z * N;

        assert(Magnitude(hvector) >= 1.0-1e-06 &&
	       Magnitude(hvector) <= 1.0+1e-06);

	cos_HI = -DotProd(hvector,dirEye);

	if(cos_HI < 0 && fabs(cos_HI) < 1e-06) cos_HI = WTINY;
        if(cos_HI >= 0)
        {
          // Compute Light ray direction
          dirRes = (2.0 * cos_HI) * hvector + dirEye;

          assert(Magnitude(dirRes) >= 1.0-1e-06 &&
	         Magnitude(dirRes) <= 1.0+1e-06);
          outCos = DotProd(N, dirRes);  // (angle btw the light ray and normal)

	  // reject unappropriate samples
	  if(outCos < 0)
	    rejected = true;
        }
        else //(outCos < 0)
        {
          rejected = true;
        }
      }
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      float cos3_HN = pow(cos_HN, 3.0f);

      // set PDF value
      contrib->prob = exp(-tan2_HN/slopeC_sq)/(slopeC_sq * M_PI * cos3_HN);

      // compute distribution function for the rest distribution members
      float distribution = 0;
      if(slopes_size > 1)
      {
        float cos4_HN = cos3_HN * cos_HN;
        for(int i = 0; i < slopes_size; i++)
        {
          if(i != choose_dist)
	  {
	    float memb    = exp( -tan2_HN/a2[i] )/(a2[i]);
 	    distribution += weights[i] * memb;
	  }
        }
	distribution /= cos4_HN;
      }
      incEyeCos = fabs(incEyeCos);
      if(incEyeCos < WTINY) incEyeCos = WTINY;

      // Compute geometrical attenuation factor (accounting for
      // shadowing and masking)
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      float part1 = 2.0 * cos_HN * incEyeCos/cos_HI;
      float part2 = 2.0 * cos_HN * outCos/cos_HI;
      float G = 2.0;
      if(part1 < part2)
	G = ((part1 < 1.0) ? part1 : 1.0);
      else
        G = ((part2 < 1.0) ? part2 : 1.0);
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      float m1     = G  / (4.0 * incEyeCos);
      // choosen distribution member
      float multi  = chWeight/ cos_HN;
      // the `rest` of the unchoosen distribution members
      float multi2 = distribution/(M_PI * contrib->prob);

      // Set the value of BSDF weight for choosen distribution member
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      /* weight = BRDF * cos(Inc,N)/ probability_density_of_generated
         direction or more precisely:
         weight = BRDF * cos(Inc,N)/ pdf =
                = [Fresnel * exp(-tan(H,N)^2/m^2) * G * cos(Inc,N)/
	              {4*Pi*m^2*cos(Out,N)*cos(Inc,N)*cos(H,N)^4}]/
                   [exp(-tan(H,N)^2/m^2) / {cos(H,N)^3*Pi*m^2}] =
                = ps * G / {4 * cos(N,V) * cos(H,N)}
         (m   - is standard deviation of surface slope,
	  Out - eye ray direction, Inc - light direction)      */
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      // just avoid additional multiplication on ks_local * m1
      contrib->weigth = ks_local * m1 * (multi + multi2);

       //// set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
      // type of interaction
      contrib->interaction = EE_REFLECTED_SPECULAR;
      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;
    }
    else
    { // If slope is very small, surface behaves like a mirror. However,
      // we should keep normalization parameter 1/(4*Pi*alpha^2) and
      // only approximate outgoing direction as mirror-reflected
      // direction and omit the exponent:
      /*
          alpha != 0, slopeXY -> 0  => exp(-tan(alpha)/slopeXY)/slopeXY -> 0
	  alpha -> 0, slopeXY -> 0  => exp(-tan(alpha)/slopeXY)/slopeXY ->
	                               1/slopeXY

	  In other words, exp(-tan(alpha)/slopeXY) behaves like delta function
      */
      // set the reflected ray, direction is normalized here
      dirRes = hitInfo.GetReflectedDir();
      outCos = DotProd(dirRes, N);
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
      // type of interaction
      contrib->interaction = EE_REFLECTED_MIRROR;
      // diffuseness of the bounce
      contrib->diffuseness = 0.0;

      // exp(-tan(0)/slopeXY) = 1, cos(H,N) = 1 =>
      contrib->prob = 1.0 / (slopeC_sq * M_PI);

      // For mirror cosine of incident equal to the cosine of
      // reflected ray => sqrt(outCos/incEyeCos) =>
      contrib->weigth = ks_local * 0.25;
    }
    //-------------------------------------------------------------------------

    // Additional (diffuse) components in order to provide nice and simple
    // vizualisation
    //------------------------------------------------------------------------
    if(kd_local.Luminance() > 1e-06)
       contrib->weigth += kd_local * (outCos/(M_PI*contrib->prob));

#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
  } // for i
  // delete temporary arrays
  delete [] a2;
  delete [] weights;

 // return number of generated samples
  assert(num == maxRays);
  return num;
}
//========================================================================


// Gather all parts for stochastic sampling together. Scaling approach
// of random values for stochastic methods was adopted from [3]
//========================================================================
int
CCookTorrance::SpawnStochasticRays(CContributionArray  &carray,
		                    const CHitPointInfo &hitInfo,
                                    CSamplingInfo       &samplingInfo,
				    EDistributionType distType) const
{
  int maxRays = samplingInfo.GetMaxSamples();
  assert(maxRays > 0);

  // TODO
  // It should be changed for Halton sequence. But for simple random
  // numbers it is in this way.
  diffuse_rnd.clear();
  specular_rnd.clear();

  // store the seed to remain the same if possible, depending
  // on the random generator.
  samplingInfo.StoreSeed();

  CValues kd_local(CValues::EE_Weigth3);	// diffuse component
  CValues ks_local(CValues::EE_Weigth3);	// specular component
  CVector2D slopeXY(0,0), slopeXY_sq(0,0);

  // compute diffuse and specular
  float akd = 0.0f, aks = 0.0f, tmp_ratio = 0;

  // compute albedo for diffuse part
  //----------------------------------------------------------------------
  if (_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
    hitInfo.GetTextureValue(kd_local, *_diffuse, CMappingBase::EE_KD);
    akd = kd_local.Luminance();
  }

  // compute albedo for specular part
  //----------------------------------------------------------------------
  if (_specular && (distType & EE_REFLECTED_SPECULAR)) {
    // Check whether fresnel parameters are set
    bool use_fresnel = _specular->IsFresnel();

    if(use_fresnel)
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
    else{
	  WARNING_INFO("CCookTorrance::SpawnStochasticRays() : Specular term of cook torrance should be fresnel type! 0 rays has been spawned.");
	  return 0;
	}

    aks = ks_local.Luminance();

    if(!_slope)
    {
	  WARNING_INFO("CCookTorrance::SpawnStochasticRays() : slope structure should be set to compute specular part of Cook Torrance model.");
      return 0;
    }

    // Check whether slope structure (distribution members) is set
    CTexture *slope = _slope->GetSlope();
    if(!slope)
    {
      WARNING_INFO("CCookTorrance::SpawnStochasticRays() : Glossy term requires structure 'slope' which consists of 'slopeX' term.");
      return 0;
    }
    CValues slopeXc(0,0,0);
    hitInfo.GetTextureValue(slopeXc, *slope, CMappingBase::EE_KSI);
    float slopeXf = slopeXc.Get(0);
    (slopeXf < WTINY) ? slopeXf = WTINY : slopeXf = slopeXf;

    slopeXY.xx    = slopeXf;
    slopeXY_sq.xx = slopeXf * slopeXf;
    tmp_ratio = aks/(slopeXY_sq.xx * 4.0 * M_PI);
  }

  // Complete albedo of Beard-Maxwell model
  //----------------------------------------------------------------------
  float albedo_whole = akd + aks;

  if(albedo_whole < 1e-05) {
    WARNING_INFO("CCookTorrance::SpawnStochasticRays() : Albedo of this surface is almost zero => it is the black body!");
    return 0;
  }
  
  if( albedo_whole >= 1.0 ) {
    WARNING_INFO("CCookTorrance::SpawnStochasticRays() : \nAlbedo > 1.0\n\
To keep the energy conservation requirement the total albedo\n\
(albedo_specularity + albedo_diffuse) should be less than 1.\n\
(If it is equal to 1, termination of a random walk with\nRussian Roulette will be impossible)");
  }

  // 'Distribute' (precalculate random numbers for each type of interaction:
  // diffuse, specular) maxRays rays of different type according to the
  // corresponding albedos
  //----------------------------------------------------------------------
  int numDiffuse = 0, numSpecular = 0;
  int count = maxRays;
  while(count--)
  {
     float rnd_num = samplingInfo.Sample(0);

     // This RND number comes from the User's program. In the case some event
     // besides absorption was choosen by the user application this RND number
     // fits into the range 0..albedo_whole. In order to reuse this RND number
     // later, it should be scaled (to fit into the range (0..1]) as following
     // (See for details [3], and phong.cpp file)
     if(count >= maxRays-1)
     {
       if(rnd_num <= albedo_whole)
         rnd_num /= albedo_whole;  // make first RND fall into the [0..1)
     }
     assert(rnd_num <= 1.0);

     float alpha1 = akd/albedo_whole;
     float alpha2 = aks/albedo_whole;

     if((rnd_num <= (alpha1) && alpha1 >= 1e-6) || (tmp_ratio <= 1e-06)) {
       // Reuse random number for generation of one component of the random
       // direction. Scale it to fall into the range [0..1)
       diffuse_rnd.push_back(rnd_num/alpha1);
       // Choose another component of direction randomly
       diffuse_rnd.push_back(samplingInfo.Sample(1));
       numDiffuse++;
     }
     else {
       // Reuse random number for generation of one component of the random
       // direction. Scale it to fall into the range [0..1)
       specular_rnd.push_back((rnd_num-alpha1)/alpha2);
       // Choose another component of direction randomly
       specular_rnd.push_back(samplingInfo.Sample(1));
       numSpecular++;
     }
     // TODO For a now it does not work for Halton sequence.
     samplingInfo.IncDimension(2);
     samplingInfo.RecoverSeed();
  }
  int numRays = 0;

  // The number of rays for diffuse interaction
  //----------------------------------------------------------------------
  if (numDiffuse) {
    numRays += SpawnStochasticDiffuseRays(kd_local, ks_local, slopeXY,
    		slopeXY_sq, carray, hitInfo, samplingInfo, numDiffuse);
  }

  //int indexDiffuse = carray.GetUsedItems();

  // The number of rays for specular interaction
  //----------------------------------------------------------------------
  if (numSpecular) {
    numRays += SpawnStochasticSpecularRays(ks_local, slopeXY, slopeXY_sq,
    		kd_local, carray, hitInfo, samplingInfo, numSpecular);
  }

  samplingInfo.IncDimension(2);
  // recover the seed for the next dimension to be sampled
  samplingInfo.RecoverSeed();

  //assert(indexMax == maxRays);
  return numRays;
}
//========================================================================

//========================================================================

const char* 
CCookTorrance::GetHelp() const
{
	return BrdfStrings::KCookTorranceDesc;
}
//========================================================================

int 
CCookTorrance::GetNumParameters() const
{
	return 2 + _slopesCount; // diffuse + fresnel term + number of slopes (1 at least is default)
}
//========================================================================

CValues
CCookTorrance::GetParameter(int parameterIndex) const
{
	CValues c(4);
	//float weight;

	switch (parameterIndex)
	{
	case 0 :
		_diffuse->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 1 :
		CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_specular);
		c.SetRepresentation(CValues::EE_WeigthN);
		CValues aux(CValues::EE_Weigth3);
		tex->GetColor(aux);
		c.Set(0,aux.Get(0));
		c.Set(1,aux.Get(1));
		c.Set(2,aux.Get(2));
		tex->GetIndexR(aux);
		c.Set(3,aux.Get(0));
		tex->GetExtinction(aux);
		c.Set(4,aux.Get(0));
		bool schlick;
		tex->GetApproxSchlick(schlick);
		c.Set(5,schlick);
		break;
	}
	
	if (_slope->GetSlopesSize() < (parameterIndex-2))
	{
		float weight;
		_slope->GetSlope(parameterIndex-2,weight)->GetValue(c,CHitPointInfo::constMapInfo);
		c.Set(3,weight);
	}

	return c;
}

//========================================================================

CValues
CCookTorrance::GetParameterConstraints(int parameterIndex) const
{
	CValues c(CValues::EE_WeigthN);

	switch (parameterIndex)
	{
	case 0 :
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 1 :
		_setBasicConstraints(EE_Fresnel,c);
		break;
	}

	if (_slopesCount <= (parameterIndex-1))
	{
		// set as default weighted float constraint
		_setBasicConstraints(EE_Weighted_Float,c);
		// edit for slope constraints
		c.Set(0,0); // min value
		c.Set(1,M_PI/2.f); // max value
		c.Set(2,0.5f); // default value
	}

	return c;
}
//========================================================================

const char* 
CCookTorrance::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseDesc;
	case 1 :
		return BrdfStrings::KFresnelDesc;
	}
	//	if (parameterIndex >= 2)
	return BrdfStrings::KSlopeDesc;
}
//========================================================================

CSurface::EParameterUIType
CCookTorrance::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return EE_ColorRGB;
	case 1 :
		return EE_Fresnel;
	}
	//	if (parameterIndex >= 2)
	return EE_Weighted_Float;
}
//========================================================================

const char*
CCookTorrance::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseTitle;
	case 1 :
		return BrdfStrings::KFresnelTitle;
	}

	//	if (parameterIndex >= 2)
	return BrdfStrings::KSlopeTitle;
}
//========================================================================

void
CCookTorrance::SetParameter(int parameterIndex, const CValues& value)
{
	switch (parameterIndex)
	{
	case 0 :{
				CConstantTexture* tex = new CConstantTexture(value);
				SetDiffuse(tex);
				break;
			}
	case 1 :{
				if (!_specular)
				{
					_specular = new CFresnelTexture(CValues(CValues::EE_Weigth3));
				}
				CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_specular);
				CValues auxVal(CValues::EE_Weigth3);
				auxVal.Set(0,value.Get(0));
				auxVal.Set(1,value.Get(1));
				auxVal.Set(2,value.Get(2));
				tex->SetColor(auxVal);
				auxVal.SetRepresentation(CValues::EE_Weigth1);
				auxVal.Set(value.Get(3));
				tex->SetIndexR(auxVal);
				auxVal.Set(value.Get(4));
				tex->SetExtinction(auxVal);
				bool useSchlick = value.Get(5);
				tex->SetApproxSchlick(useSchlick);
				break;
			}
				
	}

	if (parameterIndex >= 2)
	{
		if (value.GetType() != CValues::EE_Weigth2)
		{
			return;
		}
		float w = value.Get(1);
		int slopeIndex = parameterIndex - 2;
		CValues slope_value(value[0],value[0],value[0]);
		CConstantTexture* newTex = new CConstantTexture(slope_value);
		if (!_slope)
		{
			_slope = new CSlopeMicrofacets();
			_slopesCount = 1;
		}
		if (slopeIndex >= _slope->GetSlopesSize())
		{
			_slope->AddSlope(newTex,w);
			_slopesCount++;
		}
		else
		{
			_slope->Set(newTex,slopeIndex);
			_slope->SetWeight(w,slopeIndex);
		}
	}
}

__END_BRDFLIB_SOURCE
