// ===================================================================
// $Id$
//
// bear_maxwell.cpp
//
// Class: CBeardMaxwell
//
// Initial coding by Polina Kondratieva (PK), Mai 2004.
// ===================================================================

// GOLEM headers
#include "config.h"
#include "beard_maxwell.h"
#include "texture.h"
#include "functions.h"
#include <iostream>

__BEGIN_BRDFLIB_SOURCE


typedef vector<double> rnd_array;
rnd_array CBeardMaxwell::volume_rnd;

/**********************************************************************
 *
 *             CBeardMaxwell methods
 *
 **********************************************************************/
// describes the surface to a given stream
//========================================================================
void
CBeardMaxwell::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 (_volume) {
    indent(app, ind+1);
    app << "  volume   = \t";
    _volume->Describe(app, ind); app << "\n";
  }

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

    if (_fresnel) {
      indent(app, ind+1);
      app << "  fresnel  = \n";
      _fresnel->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 *CBeardMaxwell::Dup() const
{
  CBeardMaxwell *ret = CreateNewObject();

  if (_diffuse)
    ret->SetDiffuse(_diffuse);
  if (_specular)
  {
     ret->_specular = _specular;
    
     if(_slope)
        ret->_slope = _slope;
     if(_fresnel)
        ret->_fresnel = _fresnel;
  }
  if (_volume)
    ret->_volume = _volume;

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

// Computes maximal albedo
//========================================================================
bool
CBeardMaxwell::GetAlbedoMax(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
	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 && _fresnel))  {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
    c += d;
    retval = true;
  }

  // Volumetric part of albedo
  if ((distType & _microfacetReflectionType()) && _volume) {
    hitInfo.GetTextureValue(d, *_volume, CMappingBase::EE_KV);
    c += d;
    retval = true;
  }

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


// Computes albedo
//========================================================================
bool
CBeardMaxwell::GetAlbedo(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
  // As default for all the Microfacet-based models. Actually, for each
  // model aldebo is a function of the angle bewtween the incoming
  // direction and normal.
  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 && _fresnel))  {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
    c += d;
    retval = true;
  }

  // Volumetric part of albedo
  if ((distType & _microfacetReflectionType()) && _volume) {
    hitInfo.GetTextureValue(d, *_volume, CMappingBase::EE_KV);
    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";
    }
    if (_volume) {
      hitInfo.GetTextureValue(d, *_volume, CMappingBase::EE_KV);
      DEBUG << "volume = " << d << "\n";
    }
  }
#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;
}


/*=============================================================================
 Beard-Maxwell model:

 fr = pho_d/Pi + spec_part + volume_part

 dif_part  = pho_d/Pi
 
 spec_part = Fresnel(B)   pho_s * cos(H,N)^2    / 1+ Th *exp(-2*B/tau)/omega \
             ---------- * ------------------ * | --------------------------  |
             Fresnel(0)   cos(L,N)*cos(R,N)    \        1+ Th/omega         /

 volume_part = 2*pho_v/(cos(L,N)+cos(R,N))/(4*Pi)=  /// 18-10-2011 : Nekola - what is 4*Pi term? its not included in original paper
             =   pho_v/(cos(L,N)+cos(R,N))/(2*Pi)
=============================================================================*/
// Direct lighting computation
//=============================================================================
bool
CBeardMaxwell::EvaluateBSDF(CValues &result,
		const CHitPointInfo &hitInfo,
		const CVector3D     &dirLit,
		EDistributionType   distType) const
{
  bool retval      = false;
  float d          = Magnitude(dirLit);
  CVector3D N      = hitInfo.GetNormal();
  float lightCos   = DotProd(dirLit, N);
  float incEyeCos  = -(hitInfo.GetDotProd());

  CHitPointInfo hitInfoTmp = hitInfo;

  // Provide 'zero' result if Light source or viewer is under the surface 
  if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
  {
     //if(lightCos < COS_MIN_BM) lightCos    = COS_MIN_BM;
     //if(incEyeCos < COS_MIN_BM) incEyeCos  = COS_MIN_BM;
 
     if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
       Normalize(dirLit);

	 result = CValues(CValues::EE_Weigth3);

     // Diffuse part
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
       CValues rho_d;
       hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
       //result = rho_d * lightCos * M_1_PI;  // NEKOLA 18-10-2011 : why *lightCos? it is multiplied at the end of this method..
	   result += rho_d * M_1_PI;
       retval = true;
     }

     // Specular part
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if((_specular && (distType & EE_REFLECTED_GLOSSY)) &&
       (_slope && _fresnel))
     {
       CTexture *omegaT = _slope->GetSlopeX();
       CTexture *tauT   = _slope->GetSlopeY();
       if(omegaT && tauT)
       {
         CValues rho_s;
         hitInfo.GetTextureValue(rho_s, *_specular, CMappingBase::EE_KS);

         if(rho_s.Luminance() > 1e-06)
         {
           CValues fresnel, fresnel0, omegaCol, tauCol;
           CVector3D dirEye =  hitInfo.GetRay()->GetDir();
           CVector3D  H     =  Normalize(dirLit - dirEye);
           float cos_HN     =  DotProd(H, N);
           assert(cos_HN >= 0);
	  
	   float cos2_HN =  cos_HN * cos_HN;
	   float theta_h =  acos(cos_HN);
       float cos_2B  =  DotProd(dirLit,-dirEye);
	   float beta_2  =  acos(cos_2B);

           // According to the papers we should compute Fresnel reflection for
  	   // bistatic angle Beta. In other words, we consider Fresnel
	   // reflection for a single microfactet whose normal is oriented 
	   // along the H-vector.
	   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
           hitInfoTmp.SetNormal(H);
           hitInfoTmp.GetTextureValue(fresnel,*_fresnel, CMappingBase::EE_KFM);
           hitInfoTmp.GetTextureValue(fresnel0,*_fresnel, CMappingBase::EE_KM);
	   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
          
	   hitInfo.GetTextureValue(omegaCol, *omegaT, CMappingBase::EE_KSI);
       hitInfo.GetTextureValue(tauCol,   *tauT,   CMappingBase::EE_KSI);
       float omega = omegaCol.Get(0);
	   float tau   = tauCol.Get(0);

	   CValues scaled_F;
	   for(int i = 0; i < fresnel.GetCount(); i++)
	      scaled_F[i] = fresnel[i]/fresnel0[i];

       // Specular part includes lightCos in denominator => we can
	   // take it into account and do not divide/multiply on it.
	   float k = (cos2_HN / (incEyeCos * lightCos)) * 
				 ((1.0f + (theta_h/omega) * exp(-beta_2/tau))/(1.0f + theta_h/omega)) ;
	  
	   result += rho_s * scaled_F * k; 
	  
	   /* The model is little bit strange. For angles > 80 degrees, value 
	    of BRDF is extremely high. In order to avoid this, it is possible
	    to do the following. The expression  cos2_HN/(incEyeCos * lightCos)
	    has its local maximum at  theta_i = theta_r (when eye vector
	    coincides with mirror reflected vector). The value of local
	    maximum can be used as maximum value for entire function.
	    
	    float k = (1.0f + theta_h*exp(-beta_2/tau)/omega)/
	              ((1.0f + theta_h/omega));
	
            float cosineTerm =  cos2_HN/(incEyeCos * lightCos);
	    
	    if(cosineTerm > 1.0/(lightCos * lightCos)) 
	    	cosineTerm = 1.0/(lightCos * lightCos);
	    
	    result += rho_s * (scaled_F * cosineTerm * lightCos * k);
	  	    
	    If it is needed use the code above instead of current version.
	    In current version the limit of cos(85°) for incEyeCos and
	    lightCos is used. (In NEFDS all the measurments are accomplished
	    for angles between 10° and 80°)	  */
          
	   retval = true;
         }
       }
     }

     // Volumetric part
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     if(_volume) {
       CValues rho_v;
       hitInfo.GetTextureValue(rho_v, *_volume, CMappingBase::EE_KV);
       //result += rho_v * (0.5*lightCos/(M_PI*(lightCos + incEyeCos)));
	   result += rho_v * Max((2.0/(lightCos + incEyeCos)),WTINY);
       retval = true;
     }
  }
  else // Light source or viewer is 'under' the surface
  {
    retval = false;
  }
  result *= lightCos;
  return retval;
}

// Diffuse part for stochastic sampling
//========================================================================
int
CBeardMaxwell::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CVector2D           &omega_tau,
  const CValues       &ks_local,
  const CValues       &kv_local,
  const CValues       &kf0_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CHitPointInfo hitInfoTmp = hitInfo;
  CVector3D N              = hitInfo.GetNormal();
  CVector3D iEye	   = -((hitInfo.GetRay())->GetDir()); //inverse Eye ray
  float incEyeCos          = hitInfo.GetDotProd();
  
  // If light source/viewer is under the surface
  if (incEyeCos > 0.0)
    N = -N;

  incEyeCos = fabs(incEyeCos);
  if(incEyeCos < COS_MIN_BM) incEyeCos = COS_MIN_BM;

  // Make new basis with R as Z axis
  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

  // Compute data 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 (local coordinate system)
    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 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 weight
    //------------------------------------------------------------------------
    // compute = BRDF * cos(normal, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // It 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; 

    float lightCos = DotProd(dirRes, N);
    assert(lightCos >= -WTINY);
    if(lightCos < COS_MIN_BM)  lightCos  = COS_MIN_BM;

    // Additional components in order to provide nice and simple vizualisation
    //------------------------------------------------------------------------
    CValues addWeight(0,0,0);
    addWeight.SetRepresentation(CValues::EE_Weigth3);

    // Add specular part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if(ks_local.Luminance() > 1e-06) {
      CVector3D  H  =  Normalize(dirRes + iEye);
      float cos_HN  =  DotProd(H, N);
      assert(cos_HN >= 0 && cos_HN <= 1.0f);
      float cos2_HN =  cos_HN * cos_HN;
      if(cos2_HN < WTINY)  cos2_HN = WTINY;

      float cos_2B  =  DotProd(dirRes, iEye);
      float theta_h = acos(cos_HN);
      float beta_2  = acos(cos_2B);

      CValues fresnel, fresnel0, scaled_F;
      hitInfoTmp.SetNormal(H);
      hitInfoTmp.GetTextureValue(fresnel,  *_fresnel,  CMappingBase::EE_KFM);
      for(int i = 0; i < fresnel.GetCount(); i++)
         scaled_F[i] = fresnel[i]/kf0_local[i];

      /* spec_part =  cos^2(H,N) * SO                 1
                     ----------------- * lightCos * -----  = 
		            incEyeCos*lightCos               pdf

		  =  cos^2(H,N) * SO * Pi/(incEyeCos*lightCos)    */ 	  
      //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  
      float k = cos2_HN *
                (1.0f + theta_h*exp(-beta_2/omega_tau.yy)/omega_tau.xx)/
	       ((1.0f + theta_h) * incEyeCos);

      addWeight += ks_local * scaled_F * (k/contrib->prob);
    }

    // Add volume part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if(kv_local.Luminance() > 1e-06) {
      // This is already divided over the PDF = cos(L,N)/Pi
      //addWeight += kv_local * (0.5/(lightCos + incEyeCos));
	  addWeight += kv_local * (2.0*M_PI/(lightCos + incEyeCos)); /// 02-11-11 Nekola
    }

    // Add computed parts to diffuse component
    if (addWeight.Max() > 1e-06) {
      contrib->weigth += addWeight;
    }
    //-----------------------------------------------------------------------
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

  } // for i
  return maxRays;
}
//========================================================================


// Volumetric part for stochastic sampling
//========================================================================
int
CBeardMaxwell::SpawnStochasticVolumeRays(
  const CValues       &kv_local,
  CVector2D           &omega_tau,
  const CValues       &kd_local,
  const CValues       &ks_local,
  const CValues       &kf0_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CHitPointInfo hitInfoTmp = hitInfo;
  CVector3D N              = hitInfo.GetNormal();
  CVector3D iEye           = -((hitInfo.GetRay())->GetDir());
  float incEyeCos          = hitInfo.GetDotProd();
  
  if (incEyeCos > 0.0)
    N = -N;

  incEyeCos = fabs(incEyeCos);
  if(incEyeCos < COS_MIN_BM) incEyeCos = COS_MIN_BM;

  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

  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 (local coordinate system)
    CVector3D dirN;
    // Generate a sample in 2D, range [0;1)x[0;1)
    CVector2D vv(volume_rnd[2*num], volume_rnd[2*num+1]);

    // probability is set to value cos(theta)/M_PI
    contrib->prob = RND::hemisph_cos(dirN, vv);
    // transform computed sample 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_GLOSSY;
    // diffuseness of the bounce
    contrib->diffuseness = 1.0; // TODO

    float lightCos = DotProd(dirRes, N);
    assert(lightCos >= 0);

    // Set the value of BSDF weight
    //------------------------------------------------------------------------
    // weight = BRDF * cos(normal, ray.dir) / pdf =
    //        = [pv * cos(N, L)/(2 * Pi * (cos(N,L)+cos(N,R)))]/[cos(N,L)/Pi]
    //        =  pv / (2 * (cos(N,L)+cos(N,R)))
    // It 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 = kv_local * (0.5 /(lightCos + incEyeCos));
	contrib->weigth = kv_local * (2.0*M_PI /(lightCos + incEyeCos)); // 02-11-11 Nekola
    
    // Additional components in order to provide nice and simple vizualisation
    //------------------------------------------------------------------------
    CValues addWeight(0,0,0);
    addWeight.SetRepresentation(CValues::EE_Weigth3);
    
    // Add specular part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if(ks_local.Luminance() > 1e-06) {
      CVector3D  H  =  Normalize(dirRes + iEye);
      float cos_HN  =  DotProd(H, N);
      assert(cos_HN >= 0 && cos_HN <= 1.0f);
      float cos2_HN = cos_HN * cos_HN;
      if(cos2_HN < WTINY)  cos2_HN = WTINY;

      float cos_2B  =  DotProd(dirRes, iEye);
      float theta_h = acos(cos_HN);
      float beta_2  = acos(cos_2B);

      CValues fresnel, fresnel0, scaled_F;
      hitInfoTmp.SetNormal(H);
      hitInfoTmp.GetTextureValue(fresnel,  *_fresnel,  CMappingBase::EE_KFM);
      
      for(int i = 0; i < fresnel.GetCount(); i++)
        scaled_F[i] = fresnel[i]/kf0_local[i];

      float k = cos2_HN *
      		(1.0f + theta_h*exp(-beta_2/omega_tau.yy)/omega_tau.xx) /
	       ((1.0f + theta_h/omega_tau.xx) * incEyeCos);

      addWeight += ks_local * scaled_F * (k/contrib->prob);
    }

    // Add diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if(kd_local.Luminance() > 1e-06) {
      addWeight += kd_local;
    }

    if (addWeight.Max() > 1e-06) {
      contrib->weigth += addWeight;
    }
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

  } // for i
  return maxRays;
}
//========================================================================


// Specular part for stochastic sampling
//========================================================================
int
CBeardMaxwell::SpawnStochasticSpecularRays(
  const CValues       &ks_local,
  CVector2D           &omega_tau,
  const CValues       &kd_local,
  const CValues       &kv_local,
  const CValues       &kf0_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  int num = 0;		// the number of rays created

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

  if (incEyeCos > 0.0)
    N = -N;
  
  incEyeCos = fabs(incEyeCos);
  if(incEyeCos < COS_MIN_BM) incEyeCos = COS_MIN_BM;

  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

  float diffuseness = omega_tau.xx * omega_tau.yy;   // TODO

  // Compute outgoing directions
  //---------------------------------------------------------------------------
  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);

    CVector3D dirRes(0,0,1), hvector(0,0,1);

    float lightCos = -1;
    float cos_HI   = -1;
    float cos_HN   = 0, cos2_HN = 0, sin_HN = 1;
    bool rejected  = false;

    // rejection sampling
    //--------------------------------------------------------------------
    while(lightCos < 0) {
      // TODO: Halton sequence
      //while(lightCos < 0 || cos_HI < 0) {
      //while(cos_HI < 0) {
        //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(aux_rnd, samplingInfo.Sample(1));
      CVector2D h;
      if( !rejected )
        h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
      else
        h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

      cos_HN   = pow((1.0-h.xx),(1.0/3.0));
      cos2_HN  = cos_HN * cos_HN;
      sin_HN   = sqrt(1.0f - 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 = 0;
      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);
        lightCos = DotProd(N, dirRes);  // (angle btw the light ray and normal)

	// reject unappropriate samples
	if(lightCos < 0)
	  rejected = true;
      }
      else //(outCos < 0)
      {
        rejected = true;
      }
    }
    //--------------------------------------------------------------------
    float theta_h  = acos(cos_HN);
    float cos_2B   = DotProd(-dirEye, dirRes);
    float beta_2   = acos(cos_2B);
    
    if(cos2_HN < WTINY) cos2_HN = WTINY;
    
    contrib->prob = cos2_HN * 1.5/M_PI;    // set PDF value

    // Set the value of BSDF weight
    /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    weight = BRDF/PDF * cos(L,N);

    weight = Fresnel(B)   pho_s * cos(H,N)^2   /(1+ Th *exp(-2*B/tau)/omega)\
            ---------- * ------------------ * | --------------------------  | *
             Fresnel(0)   cos(L,N)*cos(R,N)    \       (1+ Th/omega)        /
	    ---------------------------------------------------------------
	                    cos(H,N)^2 * 3/(2*Pi)

	      * cos(L,N)

           = Fresnel(B)    pho_s      / (1+ Th *exp(-2*B/tau)/omega)*2*Pi/3 \
            ---------- * --------- * | -----------------------------------  |
             Fresnel(0)   cos(R,N)    \          (1+ Th/omega)             /
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

    // According to the papers we should compute Fresnel reflection for
    // bistatic angle Beta. In other words, we consider Fresnel
    // reflection for a single microfactet whose normal is oriented along
    // the H-vector.
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CValues fresnel, scaled_F;
    hitInfoTmp.SetNormal(hvector);
    hitInfoTmp.GetTextureValue(fresnel,  *_fresnel,  CMappingBase::EE_KFM);
    
    for(int i = 0; i < fresnel.GetCount(); i++)
      scaled_F[i] = fresnel[i]/kf0_local[i];

    assert(lightCos >= 0);
  
    float multi = (1.0f+theta_h*exp(-beta_2/omega_tau.yy)/omega_tau.xx)*M_PI/
                 ((1.0f+theta_h/omega_tau.xx)* incEyeCos * 1.5);
     
    // set the weight of the sample
    contrib->weigth = ks_local * scaled_F * multi; 

    // 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;

    CValues addWeight(0,0,0);
    addWeight.SetRepresentation(CValues::EE_Weigth3);

    // Additional components in order to provide nice and simple vizualisation
    //------------------------------------------------------------------------
    if(kd_local.Luminance() > 1e-06)
    {
      addWeight += kd_local;
    }

    if(kv_local.Luminance() > 1e-06) {
      if(lightCos < COS_MIN_BM)
         lightCos = COS_MIN_BM;
      addWeight += kv_local * (2.0*M_PI / (lightCos + incEyeCos));
    }

    if(addWeight.Luminance() > 1e-06)
      contrib->weigth += addWeight * (lightCos/(M_PI*contrib->prob));
    //------------------------------------------------------------------------

#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
  } // for i
  assert(num == maxRays);
  return num;
}
//========================================================================


// Gather all parts for stochastic sampling together. Scaling approach
// of random values for stochastic methods was adopted from [5]
//========================================================================
int
CBeardMaxwell::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 works in this way.
  diffuse_rnd.clear();
  specular_rnd.clear();
  volume_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
  CValues kv_local(CValues::EE_Weigth3);    // volume component
  CValues kf0_local(CValues::EE_Weigth3);   // fresnel term
  CVector2D omega_tau(0,0);

  float akd = 0.0f, aks = 0.0f, akv = 0.0f, af = 0.0f;

  // 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)) {
    hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
    aks = ks_local.Luminance();

    // Check whether slope structure (coeffitients omega, tau) is set
    if(!_slope)
    {
	  WARNING_INFO("CBeardMaxwell::SpawnStochasticRays() :\n\
Slope structure should be set to compute specular part of Ward model.");
      return 0;
    }

    // Check whether fresnel parameters are set
    if(_fresnel)
    {
      hitInfo.GetTextureValue(kf0_local, *_fresnel, CMappingBase::EE_KM);
      assert(kf0_local[0] > 0 && kf0_local[1] > 0 && kf0_local[2] > 0);
    }
    else
	{
	  WARNING_INFO("CBeardMaxwell::SpawnStochasticRays() :\n\
Glossy term requires parameter for Fresnel reflection.");
      return 0;
    }
    af = kf0_local.Luminance();
    assert(af >=0 && af <= 1.0f);

    CTexture *omegaT = _slope->GetSlopeX();  // get omega texture
    CTexture *tauT   = _slope->GetSlopeY();  // get tau texture
    if(!(omegaT && tauT))
    {
	  WARNING_INFO("CBeardMaxwell::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope'\n\
which consists of 'omega' and 'tau' terms.");
      return 0;
    }
    // get omega and tau color values for particular point on the surface 
    CValues omegaC(0,0,0), tauC(0,0,0);
    hitInfo.GetTextureValue(omegaC, *omegaT, CMappingBase::EE_KSI);
    hitInfo.GetTextureValue(tauC, *tauT, CMappingBase::EE_KSI);
    
    // for a now Beard-Maxwell model does not support difraction =>
    // omega, tau just floating point numbers 
    float omega = omegaC.Get(0);
    float tau   = tauC.Get(0);

    // In order to avoid division by 0
    (omega < WTINY) ? omega = WTINY : omega = omega;
    (tau < WTINY) ? tau = WTINY : tau = omega;
    omega_tau.xx = omega;
    omega_tau.yy = tau;
  }

  // compute albedo for volumetric part
  //----------------------------------------------------------------------
  if (_volume && (distType & EE_REFLECTED_SPECULAR)) {
    hitInfo.GetTextureValue(kv_local, *_volume, CMappingBase::EE_KV);
    akv = kv_local.Luminance();
  }

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

  if(albedo_whole < 1e-05) {
	WARNING_INFO("CBeardMaxwell::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
    return 0;
  }
  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CBeardMaxwell::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, volumetric) maxRays rays of different type according
  // to the corresponding albedos
  //----------------------------------------------------------------------
  int numDiffuse = 0, numSpecular = 0, numVolume = 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 [5], 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;
     float alpha3 = akv/albedo_whole;

     if((rnd_num <= alpha1) && (alpha1 >= 1e-6)) {
       // 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 {
       if((rnd_num <= (alpha1+alpha2)) && (alpha2 >= 1e-6))
       {
         // 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++;
       }
       else
       { // Reuse random number for generation of one component of the
	 // random direction. Scale it to fall into the range [0..1)
         volume_rnd.push_back((rnd_num-alpha1-alpha2)/alpha3);
         // Choose another component of direction randomly
         volume_rnd.push_back(samplingInfo.Sample(1));
         numVolume++;
       }
     }
     // 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, omega_tau, ks_local,
                        kv_local, kf0_local, carray, hitInfo, samplingInfo,
		        numDiffuse);
  }

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

  // The number of rays for volumetric interaction
   //----------------------------------------------------------------------
 if (numVolume) {
    numRays += SpawnStochasticVolumeRays(kv_local, omega_tau, kd_local,
			ks_local, kf0_local, carray, hitInfo, samplingInfo,
			numVolume);
  }
  int indexMax = carray.GetUsedItems();
  assert(indexMax > 0);

  // TODO distributed ray tracing
/*  for(int i = 0; i < indexMax; i++) {
    CContributionItem* ci = carray.GetItem(i);
    float mult;
    if(i < indexDiffuse) {
      mult = akd;
    }
    else {
      mult = aks;
      assert(numSpecular > 0);
      ci->weigth *= (1.0 / ((float)numSpecular));
    }
    mult /= albedo_whole;
    //    ci->weigth *= (1.0 / (float)maxRays);
    //ci->weigth *= (1.0/mult);
    //ci->prob *= mult;
  } // for*/
  samplingInfo.IncDimension(2);
  // recover the seed for the next dimension to be sampled
  samplingInfo.RecoverSeed();

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

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

const char* 
CBeardMaxwell::GetHelp() const
{
	return BrdfStrings::KBeardMaxwellDesc;
}
//========================================================================

int 
CBeardMaxwell::GetNumParameters() const
{
	return 6; // diffuse,specular,slope X,slope Y,volume,fresnel
}
//========================================================================

CValues
CBeardMaxwell::GetParameter(int parameterIndex) const
{
	CValues c(0);

	float w;
	switch (parameterIndex)
	{
		case 0 : //diffuse
			_diffuse->GetValue(c,CHitPointInfo::constMapInfo);
			break;
		case 1 : //specular
			_specular->GetValue(c,CHitPointInfo::constMapInfo);
			break;
		case 2 :
			c.SetRepresentation(CValues::EE_WeigthN);
			_slope->GetSlope(0,w)->GetValue(c,CHitPointInfo::constMapInfo);
			c.Set(3,w);
			break;
		case 3 :
			c.SetRepresentation(CValues::EE_WeigthN);
			_slope->GetSlope(1,w)->GetValue(c,CHitPointInfo::constMapInfo);
			c.Set(3,w);
			break;
		case 4 :
			_volume->GetValue(c,CHitPointInfo::constMapInfo);
			break;
		case 5 :
			CFresnelTexture* tex = static_cast<CFresnelTexture*>(_fresnel);
			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;
	}
	return c;
}
//========================================================================

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

	switch (parameterIndex)
	{
		case 0 : //diffuse
		case 1 : //specular
			_setBasicConstraints(EE_ColorRGB,c);
			break;
		case 2 :
		case 3 : // slopes
			c = CValues(3);
			c.Set(0,0); // min deviation
			c.Set(1,M_PI/2.f); // max deviation - Pi/2 radians 
			c.Set(2,0.5f); // default 
			break;
		case 4 : // volume
			_setBasicConstraints(EE_ColorRGB,c);
			break;
		case 5 : // fresnel
			_setBasicConstraints(EE_Fresnel,c);
	}
	return c;
}
//========================================================================

const char* 
CBeardMaxwell::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KDiffuseDesc;
		case 1 :
			return BrdfStrings::KSpecularDesc;
		case 2 :
			return BrdfStrings::KSlopeXDesc;
		case 3 :
			return BrdfStrings::KSlopeYDesc;
		case 4 :
			return BrdfStrings::KVolumeDesc;
		case 5 :
			return BrdfStrings::KFresnelDesc;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CBeardMaxwell::GetParameterSyntax(int parameterIndex) const
{
	
	switch (parameterIndex)
	{
		case 0 : //diffuse
		case 1 : //specular
			return EE_ColorRGB;
		case 2 :
		case 3 :
			return EE_Float;
		case 4 : // volume
			return EE_ColorRGB;
		case 5 : // fresnel
			return EE_Fresnel;
		default :
			return EE_Unknown;
	}
}
//========================================================================

const char*
CBeardMaxwell::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KDiffuseTitle;
		case 1 :
			return BrdfStrings::KSpecularTitle;
		case 2 :
			return BrdfStrings::KSlopeXTitle;
		case 3 :
			return BrdfStrings::KSlopeYTitle;
		case 4 :
			return BrdfStrings::KVolumeTitle;
		case 5 :
			return BrdfStrings::KFresnelTitle;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CBeardMaxwell::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex<2) // parent parameter
	{
		CSurfaceMicrofacet::SetParameter(parameterIndex, value);
		return;
	}

	switch (parameterIndex)
	{
		case 2 :
			{
				CValues slope_value(value[0],value[0],value[0]);
				CConstantTexture* newTex = new CConstantTexture(slope_value);
				if (!_slope)
				{
					_slope = new CSlopeMicrofacets();
				}
				_slope->Set(newTex,0);
				_slope->SetWeight(0.5f,0); // doesnt care in this model
			}
			break;
		case 3 :
			{
				CValues slope_value(value[0],value[0],value[0]);
				CConstantTexture* newTex = new CConstantTexture(slope_value);
				if (!_slope)
				{
					_slope = new CSlopeMicrofacets();
				}
				_slope->Set(newTex,1);
				_slope->SetWeight(0.5f,1);
				break;
			}
		case 4 :
			{
				CConstantTexture* newTex = new CConstantTexture(value);
				SetVolume(newTex);
				break;
			}
			
		case 5 :
			if (!_fresnel)
			{
				_fresnel = new CFresnelTexture(CValues(CValues::EE_Weigth3));
			}
			CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_fresnel);
			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;
	}
}
//========================================================================


__END_BRDFLIB_SOURCE
