// ===================================================================
// $Id$
//
// ward.h
//    (Hemi)spherical harmonics representation of BRDFs.
//
// Class:
//
// Initial coding by Jaroslav Krivanek, January 2004.

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

__BEGIN_BRDFLIB_SOURCE


/**********************************************************************
 *
 *             CSchlickSurvey94 methods
 *
 **********************************************************************/
// describes the surface to a given stream
//========================================================================
void
CSchlickSurvey94::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";
    }

    if (_anisotropic) {
      indent(app, ind+1);
      app << "  anisotropy    = true\n";
    }
    else
    {
      indent(app, ind+1);
      app << "  anisotropy    = false\n";
    }

    int fresnelSize = _termsFresnel.size();
    assert(fresnelSize > 0);

    indent(app, ind+1);
    app << "  Fresnel terms:\n";

    for(int i = 0; i < fresnelSize; i++)
    {
      indent(app, ind+1);
      app << "  term[" << i << "]:\n";
      _termsFresnel[i]->Describe(app, ind);
      app << "\n";
    }
    app << "\n";
  }

  app << flush;

  return;
}


/*
  Integral of this function 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
CSchlickSurvey94::GetAlbedo(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
  //CVector3D N       = hitInfo->GetNormal();
  //CVector3D dirInit = hitInfo->GetRay()->GetDir();

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

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

  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;
    //cout << "albedo from GetAlbedo = " << albedo << "\n";
    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;
//  GetAlbedoMax(albedo, hitInfo, distType);
}


//#ifndef WTINY
//#define WTINY 1e-06
//#endif
//========================================================================
/*
                / w_i * Distrib_i * Fresnel_i(Light) * GeometricFactor_i  \
   BRDF  = Sum | -------------------------------------------------------- |
               \            4 * Pi * lightCos * cosEye                   /

  w0 + w1 + ... + wN = 1,  delta = angle(H, N), phi = angle(H_prj,X)

  Distrib_i =                 1
              -------------------------------,  m, n - roughness in X,Y dir
	      m*n*(p + cos(delta)^2*(1-p))^2

	      p = cos(phi)^2*(1/m^2-1/n^2) + (1-cos(phi)^2)/n^2

  (see "A Survey of Shading and Reflectance Models", Schlick, 1994)
*/
bool
CSchlickSurvey94::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
  //cout << "evaluateBSDF 0 Iso\n";
  bool retval = false;
  float d = Magnitude(dirIn);
  CVector3D N       =  hitInfo.GetNormal();
  CVector3D dirInit =  (hitInfo.GetRay())->GetDir();
  float lightCos    =  DotProd(dirIn, N);
   
  float incEyeCos   =  -(hitInfo.GetDotProd());
  CVector3D anisAxis;

  if(_anisotropic)
   anisAxis =  hitInfo.GetAnisAxis();

  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(dirIn);

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

    if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) {
      int slopes_size  = _slope->GetSlopesSize();
      int weights_size = _slope->GetWeightsSize();
      int fresnel_size = _termsFresnel.size();

      if(fresnel_size <= 0)
      {
		  WARNING_INFO("CSchlickSurvey94::EvaluateBSDF() : \n\
Schlick BRDF model must include one or number of slopes Fresnel functions.");
		  return false;
      }


      //    Schlick model supports anisotropy and multi-sloppiness.
      // But for the sake of clarity and simplicity we decided to use
      // random multi-sloppiness instead of support of several axes of
      // anisotropy.
      //    Therefore, if 2 slopes are specified and parameter _anisotropic
      // is true, model is considered as anisotropic and slopes are oriented
      // along the tangent and binormal directions. Otherwise, model is
      // considered isotropic and slopes are oriented in random directions.
      if(slopes_size > 0 && (slopes_size == weights_size))
      {
        float *a2       = new float[slopes_size];
	float *weights  = new float[slopes_size];
	float sumf      = 0;
	CValues sumcol(0,0,0);
        sumcol.SetRepresentation(CValues::EE_Weigth3);
        vector<CValues> termsF;

	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;

	  if(fresnel_size > 1)
	  {
	    CValues termF;
            hitInfo.GetTextureValue(termF, *_termsFresnel[i], CMappingBase::EE_KFM);
	    termsF.push_back(termF);
	  }
	}

	if(fresnel_size == 1)
	{
	  CValues termF;
          hitInfo.GetTextureValue(termF, *_termsFresnel[0], CMappingBase::EE_KFM);
	  termsF.push_back(termF);
	}

	CValues rho_s;
        hitInfo.GetTextureValue(rho_s, *_specular, CMappingBase::EE_KS);

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

        if(cos_HN < WTINY)  cos_HN = WTINY;
	float cos2_HN = cos_HN * cos_HN;

        if(_anisotropic)
	{
	  // Now let us get phi polar angle for half-vector
	  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          CVector3D lx(1,0,0), ly(0,1,0), lz(0,0,1), hl(1,0,0);
	  MakeLocalCS(N, anisAxis, lx, ly, lz);

	  hl[0] = DotProd(H, lx);
	  hl[1] = DotProd(H, ly);
	  hl[2] = DotProd(H, lz);
	  assert(Magnitude(hl) > 1.0-1e-06 && Magnitude(hl) < 1.0+1e-06);
          //cout << "local H = " << hl;

          float hx2 = hl[0]*hl[0], hy2 = hl[1]*hl[1];
	  float len_hprj2 = hx2 + hy2;
	  float cos2_phi   = hx2 / len_hprj2;
	  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

          float i_aX2 = 1.0/a2[0];
          float i_aY2 = 1.0/a2[1];
	  float aXaY   = sqrt(a2[0]*a2[1]);
	  float coef1  = sqrt(aXaY * M_2_PI);
	  float coef2  = 1.0 - coef1;

          //cout << "~~~~~~~~~~~~ Anisotropic!!!!!!!!!!!!!!!!!!\n";
	  //cout << "~~~~~~~~~~~~ cos2_phi  = " << cos2_phi << "\n";
	  //cout << "~~~~~~~~~~~~ aXaY      = " << aXaY << "\n";
	  //cout << "~~~~~~~~~~~~ termsF[0] = " << termsF[0] << "\n";

	  float anisoMemb    = cos2_phi*(i_aX2 - i_aY2) + i_aY2;
          float denominator  = anisoMemb + cos2_HN*(1.0 - anisoMemb);
	  float distribution = 1.0/(aXaY*denominator*denominator);

	  sumf   = distribution / ((lightCos *coef2 + coef1)*
	                           (incEyeCos*coef2 + coef1)*
				    4.0*M_PI);
	  sumcol = sumf * termsF[0];
        }
        else
	{
          //cout << "~~~~~~~~~~~~ Isotropic!!!!!!!!!!!!!!!!!!\n";
	  for(int i = 0; i < slopes_size; i++)
	  {
	    float coef1  = sqrt(a2[i] * M_2_PI);
	    float coef2  = 1.0 - coef1;

            float denominator  = 1.0 + cos2_HN*(a2[i] - 1.0);
	    float distribution = a2[i]/(denominator*denominator);
	    //cout << "~~~~~~~~~~~~ weights["<< i << "] = " << weights[i] << "\n";
	    //cout << "~~~~~~~~~~~~ a2[" << i << "] = " << a2[i] << "\n";

            if(fresnel_size > 1)
	    {
	      sumcol += termsF[i] * (weights[i] * distribution /
	                ((lightCos *coef2 + coef1)*(incEyeCos*coef2 + coef1)));
	      //cout << "~~~~~~~~~~~~ termsF["<< i << "] = " << termsF[i] << "\n";
	    }
	    else
	    {
	      sumf   += weights[i] * distribution /
	                ((lightCos *coef2 + coef1)*(incEyeCos*coef2 + coef1));
	    }

	  }

	  if(fresnel_size == 1)
	  {
            sumcol  = termsF[0] * (sumf / (4.0*M_PI));
            //cout << "~~~~~~~~~~~~ 1->termsF[0] = " << termsF[0] << "\n";
	  }
	  else
            sumcol *= 1.0/(4.0*M_PI);
	}

        //cout << "m2 = " << a2 << ", tan(H,N)^2 = " << tan2_HN << ", exp(-tan^2/m^2) = " << ex << "\n";
        //cout << "coefficient = " << k << ", ";
        result += rho_s * sumcol * lightCos;
        retval = true;
        delete [] a2;
	delete [] weights;
      }
    }

    /*if(retval)
    {
      //cout << "result = " << result << ", after cos: ";
      result *= lightCos;
      //cout << result << "\n";
    }*/
  }
  else
  {
    retval = false;
  }
  //cout << "evaluateBSDF 1\n";
  return retval;
}

//========================================================================
int
CSchlickSurvey94::SpawnStochasticDiffuseRays(
  const CValues         &kd_local,
  const CValues         &ks_local,
  const vector<CValues> &termsF,
  CContributionArray    &carray,
  const CHitPointInfo   &hitInfo,
  CSamplingInfo         &samplingInfo,
  int                   maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  //cout << "normal = " << N << "\n";
  float incEyeCos = hitInfo.GetDotProd();
  if (incEyeCos > 0.0)
    N = -N;

  incEyeCos = fabs(incEyeCos);
  // It is considered as X axis for local coordinate system
  //CVector3D anisAxis = hitInfo.GetAnisAxis();
  //assert(Magnitude(anisAxis) >= 1.0 - 1e-06 &&
  //       Magnitude(anisAxis) <= 1.0 + 1e+06);

  CVector3D lx(1,0,0), ly(0,1,0), lz(0,0,1), hl(1,0,0);
  CVector3D anisAxis(1,0,0);

  if(_anisotropic)
  {
    anisAxis =  hitInfo.GetAnisAxis();
    MakeLocalCS(N, anisAxis, lx, ly, lz);
  }
  else
  {
    lx = ArbitraryNormal(N); // perpendicular to N
    ly = CrossProd(N, lx);   // orthonormal base with N and lx
    lz = N;
  }

  // number of rays created
  int num = 0;
  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;
  int fresnel_size = termsF.size();
  float i_aX2 = 1, i_aY2 = 1 , aXaY = 1;
  //vector<CValues> termsF;

  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;
      if(!_anisotropic)
        avg_w   += weight;

      /*if(fresnel_size > 1)
      {
        CValues termF;
        hitInfo.GetTextureValue(termF, *_termsFresnel[i], CMappingBase::EE_KFM);
        termsF.push_back(termF);
      }*/
    }
    if(!_anisotropic)
    {
      avg_w     /= slopes_size;
    }
    else
    {
      i_aX2 = 1.0/a2[0];
      i_aY2 = 1.0/a2[1];
      aXaY   = sqrt(a2[0]*a2[1]);
    }
  }

 /* if(fresnel_size == 1)
  {
    CValues termF;
    hitInfo.GetTextureValue(termF, *_termsFresnel[0], CMappingBase::EE_KFM);
    termsF.push_back(termF);
  }*/

  for(int i = maxRays; i--; num++) {
    //cout << "3pd " << i << "\n";
    //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 to real coordinates
    CVector3D dirRes = dirN.x * lx + dirN.y * ly + dirN.z * lz ;
    // 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);
    // Cosine between the incoming angle and normal is already taken into
    // account!
    contrib->weigth = kd_local; // * M_1_PI;

    if(slopes_size > 0 && (slopes_size == weights_size)) {
      float sumf = 0;
      CValues sumcol(0,0,0);
      sumcol.SetRepresentation(CValues::EE_Weigth3);
      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 cos2_HN =  cos_HN * cos_HN;

      if(avg_w > 1e-06)
      {
        if(ks_local.Luminance()/avg_w > 1e-06)
	{
	  for(int i = 0; i < slopes_size; i++)
	  {
	    float coef1  = sqrt(a2[i] * M_2_PI);
	    float coef2  = 1.0 - coef1;
            //cout << "~~~~~~~~~~~~ a2[" << i << "] = " << a2[i] << "\n";

            float denominator  = 1.0 + cos2_HN*(a2[i] - 1.0);
	    float distribution = a2[i]/(denominator*denominator);

            if(fresnel_size > 1)
	    {
	      sumcol += termsF[i] * (weights[i] * distribution /
	                ((lightCos *coef2 + coef1)*(incEyeCos*coef2 + coef1)));
              //cout << "~~~~~~~~~~~~ 1-> termsF[" << i << "] = " << termsF[i] << "\n";
	    }
	    else
	      sumf   += weights[i] * distribution /
	                ((lightCos *coef2 + coef1)*(incEyeCos*coef2 + coef1));
	  }

	  if(fresnel_size == 1)
	  {
            sumcol  = termsF[0] * (sumf / (4.0*M_PI));
            //cout << "~~~~~~~~~~~~ 1-> termsF[0] = " << termsF[0] << "\n";
          }
	  else
            sumcol *= (0.25/M_PI);

        //cout << "after diffuse -> weight = " << contrib->weigth << ", probability = " << contrib->prob << ", pure = " << contrib->weigth * (contrib->prob/lightCos) << "\n";
        }
      }
      else
      {
        if(_anisotropic)
	{
	  // Now let us get phi polar angle for half-vector
	  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          float coef1  = sqrt(aXaY * M_2_PI);
          float coef2  = 1.0 - coef1;

	  hl[0] = DotProd(H, lx);
	  hl[1] = DotProd(H, ly);
	  hl[2] = DotProd(H, lz);
	  assert(Magnitude(hl) > 1.0-1e-06 && Magnitude(hl) < 1.0+1e-06);
          //cout << "local H = " << hl;

          float hx2       = hl[0]*hl[0], hy2 = hl[1]*hl[1];
	  float len_hprj2 = hx2 + hy2;
	  float cos2_phi  = hx2 / len_hprj2;
	  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	  float anisoMemb    = cos2_phi*(i_aX2 - i_aY2) + i_aY2;
          float denominator  = anisoMemb + cos2_HN*(1.0 - anisoMemb);
	  float distribution = 1.0/(aXaY*denominator*denominator);

	  sumcol  = termsF[0] * (distribution / ((lightCos *coef2 + coef1)*
	                                         (incEyeCos*coef2 + coef1)*
					          4.0*M_PI));
          //cout << "~~~~~~~~~~~~ Anisotropic Diffuse!\n";
	  //cout << "~~~~~~~~~~~~ 1-> termsF[0] = " << termsF[0] << "\n";
          //cout << "~~~~~~~~~~~~ cos2_phi = " << cos2_phi << "\n";
          //cout << "~~~~~~~~~~~~ i_aX2 = " << i_aX2 << ", i_aY2" << i_aY2 << "\n";
	}
      }

      //cout << "~~~~~~~~~~~~ sumcol = " << sumcol << "\n";
      // specular = spec_part * cosLight / PDF =
      //            spec_part * cosLight / (cosLight/PI) = spec_part * PI
      contrib->weigth += ks_local * sumcol * M_PI;
    }

    //if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    //  hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
    //  contrib->weigth += ks_local;
    //}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

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

int
CSchlickSurvey94::SpawnStochasticSpecularRays(
  const CValues         &ks_local,
  const CValues         &kd_local,
  const vector<CValues> &termsF,
  CContributionArray    &carray,
  const CHitPointInfo   &hitInfo,
  CSamplingInfo         &samplingInfo,
  int                   maxRays) const
{
  //cout << "==================================\n";
  // the number of samples generated at the end of this function
  int num = 0;
  //int countSamples = 0;

  float incEyeCos = hitInfo.GetDotProd();
  //cout << "incEysCos = " << incEyeCos;
  CVector3D N = hitInfo.GetNormal();

  // now holds DotProd(N, incomingRay) < 0.0
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  CVector3D anisAxis(1,0,0);
  CVector3D lx(1,0,0), ly(0,1,0), lz(0,0,1), hl(1,0,0);

  if (incEyeCos > 0.0)
  {
    N = -N;
    //cout << "cos > 0!\n";
  }

  if(_anisotropic)
  {
    anisAxis =  hitInfo.GetAnisAxis();
    MakeLocalCS(N, anisAxis, lx, ly, lz);
  }
  else
  {
    // Make new basis with R as Z axis
    lx = ArbitraryNormal(N);  // perpendicular to N        - X-axis
    ly = CrossProd(N, lx);    // perpendicular to N and lx - Y-axis
    lz = N;
  }

  // Choose one of the distribution fucntions:
  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("CSchlickSurvey94::SpawnStochasticSpecularRays() :\n\
For Schlick BRDF model non-empty slopes structure\n\
with non-zero weights must be specified");
    return 0;
  }
  CTexture *chSlope = NULL;

  // Make array of alphas for weights for choosing of distribution function
  // according to its weight.
  //vector<float> alphas;
  //alphas.push_back(0);
  int fresnel_size = termsF.size();
  float alpha      = 0;
  float *a2        = new float[slopes_size];
  float *weights   = new float[slopes_size];
  float i_aX2   = 1, i_aY2  = 1 , aXaY = 1, aXaY2 = 1, aX_aY = 1, aX_aY2 = 1;
  float i_aXaY2 = 1, i_aXaY = 1;
  float coef1 = 1, coef2 = 1;

  bool choosen;
  if(!_anisotropic) choosen = false;
  else              choosen = true;

  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;

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

  //cout << "~~~~~~~~~~ chosen distrib = " << choose_dist << "\n";
  //cout << "~~~~~~~~~~ chWeight = " << chWeight << "\n";

  CValues slopeC;
  if(_anisotropic)
  {
    i_aX2  = 1.0/a2[0];
    i_aY2  = 1.0/a2[1];
    aXaY2  = a2[0]*a2[1];
    aXaY   = sqrt(aXaY2);
    aX_aY2 = a2[0]/a2[1];
    aX_aY  = sqrt(aX_aY2);
  }
  else
  {
    hitInfo.GetTextureValue(slopeC, *chSlope, CMappingBase::EE_KSI);
    aXaY   = (slopeC.Get(0)) * (slopeC.Get(0));
    i_aXaY = 1.0/aXaY;
  }
  //cout << "~~~~~~~~~~ aXaY = " << aXaY << "\n";
  coef1  = sqrt(aXaY * M_2_PI);
  coef2  = 1.0 - coef1;

  // diffuseness of the bounce
  //!!!!!!!!!!!!!!!!!!!!!!    it is incorrect     !!!!!!!!!!!!!!!!!!!!!!
  //!!!!!!!!!!!!!!!!!!!!!! it is just temporarily !!!!!!!!!!!!!!!!!!!!!!
  float diffuseness = aXaY;//_nToDiffuseness(shininess.Sum() / shininess.GetCount());

  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(ks_local.Luminance()/aXaY > 1e-06)
    {
      float cos_HI = -1;
      CVector3D hvector(0,0,1);
      float cos2_HN = 0, cos_HN = 0, sin_HN, cos2_ph = 0, anisPart = 1.0;
      bool rejected = false;
      //int ooo = 0;

      while(outCos < 0) {
      //while(outCos < 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));

	float cos_ph = 0, sin_ph = 0;
        float phi_h = 2.0 * M_PI * h.yy;

	// See Compendium, p.31, original Ward'92 paper
        if(!_anisotropic)
	{
	  cos_ph   = cos(phi_h);
	  sin_ph   = sin(phi_h);
	  //cos2_ph  = cos_phi_h * cos_phi_h;
	  cos2_HN  = (1.0-h.xx)/((aXaY-1.0)*h.xx+1.0);
        }
	else
	{
          float cos_r  = cos(phi_h);
          float sin_r  = sin(phi_h);
	  float len    = (sqrt(cos_r*cos_r*(aX_aY2-1.0) + 1.0));
          cos_ph       = cos_r * (aX_aY) / len; //cos(phi_h);
          sin_ph       = sin_r / len;
          cos2_ph      = cos_ph * cos_ph;
          anisPart     = cos2_ph*(i_aX2 - i_aY2) + i_aY2;
	  float aPart2 = anisPart * (1.0-h.xx);
	  cos2_HN      = aPart2/(aPart2 + h.xx);
	}
        cos_HN   = sqrt(cos2_HN);
        sin_HN   = sqrt(1.0 - cos2_HN);

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

        CVector3D localDir(sin_HN * cos_ph, sin_HN * sin_ph, 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 * lx + localDir.y * ly + localDir.z * lz;

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

	cos_HI = -DotProd(hvector,dirInit);

	//if(cos_HI < 0 && fabs(cos_HI) < 1e-06) cos_HI = 0;
	if(cos_HI < 0 && fabs(cos_HI) < 1e-06) cos_HI = WTINY;

        if(cos_HI >= 0)
        {
          dirRes = (2.0 * cos_HI) * hvector + dirInit;

          //cout << "local = " << localDir << ", H = " << hvector << ", dirRes = " << dirRes << "\n";
          //cout << "dirInit = " << dirInit << "|localDir| = " << Magnitude(localDir) << ", |H| = " << Magnitude(hvector) << ",  |dirRes| = " << Magnitude(dirRes) << "\n";

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

	  if(outCos < 0)
          {
	    rejected = true;
            //cos_HI = -1;
	  }
        }
        else //(outCos < 0)
        {
          rejected = true;
	  //cout << "HI < 0!\n";
        }
        //ooo++;
        //assert(ooo < 10);
      }

      /************************************************************************
      1. Anisotropic case
      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         BRDF =                     Fresnel
	       -------------------------------------------------------------- *
	      4*PI*cos(Lit)*cos(Eye)*m*n*[anisPart+cos(delta)^2*(1-anisPart)]^2

                        cos(Lit)                 cos(Eye)
	       ------------------------ * ------------------------
	       [cos(Lit)*(1-coef)+coef]   [cos(Eye)*(1-coef)+coef]

	 where anisPart = cos(phi)^2*(1/m^2-1/n^2)+1/n^2
	       coef     = sqrt(2*m*n/PI)

         PDF =                      cos(delta)
	       ------------------------------------------------
	       PI*m*n*[anisPart + cos(delta)^2*(1-anisPart)]^2

	 => weight = BRDF*cos(Lit)/PDF =>

	 weight =          Fresnel * cos(Lit)/cos(delta)
	       ----------------------------------------------------
	       4*[cos(Lit)*(1-coef)+coef]*[cos(Eye)*(1-coef)+coef]


      2. Isotropic case
      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         BRDF =     1                        /   w_i * m_i^2 * Fresnel_i
	       ---------------------- * Sum | ---------------------------- *
	       4*PI*cos(Lit)*cos(Eye)       \  [cos(delta_i)^2*(m_i^2-1)+1]^2

                        cos(Lit)                        cos(Eye)            \
	       ---------------------------- * ----------------------------  |
	       [cos(Lit)*(1-coef_i)+coef_i]   [cos(Eye)*(1-coef_i)+coef_i] /

	 where  coef_i = sqrt(2*m_i^2/PI)

         PDF =      cos(delta_c)*m_c^2
	       ------------------------------------------------
	       PI*[cos(delta_c)^2*(m_c^2-1) + 1]^2

	 => weight = BRDF*cos(Lit)/PDF =>

	 wei_c =         w_c * Fresnel_c * cos(Lit)/cos(delta_c)
	       ----------------------------------------------------------- +
	       4*[cos(Lit)*(1-coef_c)+coef_c]*[cos(Eye)*(1-coef_c)+coef_c]

                    1           /            \
	       --------- * Sum | all i != c  |
	       4*PI*PDF        \            /

      ************************************************************************/
      //cout << "~~~~~~~~~~~~ cos2_ph = " << cos2_ph << ", anisPart = " << anisPart << "\n";
      //cout << "~~~~~~~~~~~~ cos_HN  = " << cos_HN << ", sin_HN = " << sin_HN << "\n";


      incEyeCos = fabs(incEyeCos);
      if(incEyeCos < WTINY) incEyeCos = WTINY;
      assert(outCos >= 0);
      float sumf = 0;
      CValues sumcol(0,0,0);
      sumcol.SetRepresentation(CValues::EE_Weigth3);

      if(_anisotropic)
      {
        //cout << "~~~~~~~~~~~~ Anisotropic imps!\n";
        float denom   = (anisPart + cos2_HN*(1.0-anisPart));
	contrib->prob = cos_HN / (M_PI * aXaY * denom * denom);

        float denominator = 4.0*(outCos*coef2+coef1)*(incEyeCos*coef2+coef1);

        sumcol  = termsF[0] * (1.0/ (denominator * cos_HN));
        //cout << "~~~~~~~~~~~~ 1-> termsF[0] = " << termsF[0] << "\n";
        //cout << "~~~~~~~~~~~~ PDF = " << contrib->prob << "\n";
      }
      else
      {
        //cout << "~~~~~~~~~~~~ Isotropic imps!\n";
	//float coef_c1   = sqrt(aXaY * M_2_PI);
	//float coef_c2   = 1.0 - coef_c1;
        float denom_c1  = (outCos*coef2 + coef1)*(incEyeCos*coef2 + coef1);
        float denom_c2  = (1.0 + cos2_HN*(aXaY-1.0));

	contrib->prob = cos_HN * aXaY / (M_PI * denom_c2 * denom_c2);
        //cout << "~~~~~~~~~~~~ PDF = " << contrib->prob << "\n";

        if(slopes_size > 1)
	{
	  for(int i = 0; i < slopes_size; i++)
	  {
	    if(i != choose_dist)
	    {
	      float coef_t1   = sqrt(a2[i] * M_2_PI);
	      float coef_t2   = 1.0 - coef_t1;

	      float denom1  = (outCos*coef_t2 + coef_t1)*
	                      (incEyeCos*coef_t2 + coef_t1);
	      float denom2  = 1.0 + cos2_HN*(a2[i] - 1.0);

              float distrib = a2[i]/(denom2*denom2);

	      if(fresnel_size > 1)
	      {
	        sumcol += termsF[i] * weights[i] * (distrib / denom1);
                //cout << "~~~~~~~~~~~~ 1-> termsF[" << i << "] = " << termsF[i] << "\n";
                //cout << "~~~~~~~~~~~~ a2[" << i << "] = " << a2[i] << "\n";
	      }
	      else
	        sumf += weights[i] * distrib / denom1;
	    }
	  }

	  if(fresnel_size == 1)
	  {
           //cout << "~~~~~~~~~~~~ 1-> termsF[0] = " << termsF[0] << "\n";
           sumf /= (M_PI*contrib->prob);
	  }
	  else
            sumcol *= (M_1_PI/contrib->prob);
	}
	if(fresnel_size == 1)
	{
          sumf  += chWeight / (denom_c1 * cos_HN);
	  //cout << "~~~~~~~~~~~~ choosen part = " << chWeight / (denom_c1 * cos_HN) << "\n";
	  sumcol = termsF[0] * (sumf/ 4.0);
	}
	else
	{
          sumcol += termsF[choose_dist] * (chWeight / (denom_c1 * cos_HN));
	  sumcol *= 0.25;
	}
      }
      //cout << "~~~~~~~~~~~~ sumcol = " << sumcol << "\n";
      //cout << "probability = " << contrib->prob << "\n";

      // Set the value of BSDF type
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      contrib->weigth = ks_local * sumcol;

      //cout << "coef = " << multi*contrib->prob/outCos << ", result = " << contrib->weigth*contrib->prob*(1.0/outCos) << "\n";
      //cout << "before specular -> weight = " << contrib->weigth << ", probability = " << contrib->prob << "\n"; // << "outCos = " << outCos << ", incEyeCos = " << incEyeCos << "\n";
     ////", incoming = " << dirInit << ", outgoing = " <<  dirRes << "\n";
      //// 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 kd_local(CValues::EE_Weigth3);
    //if (_diffuse) {
    //  hitInfo.GetTextureValue(kd_local, *_diffuse, CMappingBase::EE_KD);
    //}
    if(kd_local.Luminance() > 1e-06)
    {
      //cout << "difffffffffffffffffffffff = " << kd_local << ", multed = " << kd_local * (outCos/(M_PI*contrib->prob)) << "\n";
    //if(kd_local.Luminance() > 1e-05) {
       //kd_local *= 1.0/(M_PI*contrib->prob);
       contrib->weigth += kd_local * (M_1_PI/contrib->prob);
    }

    contrib->weigth *= outCos;  // (*)
    //cout << "after diffuse -> weight = " << contrib->weigth << ", without PDF,cos: " << contrib->weigth * contrib->prob*(1.0/outCos) << "\n";

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

      // the number of valid samples taken
      //countSamples++;
  } // for i
  //cout << "4p\n";

  delete [] a2;
  delete [] weights;

 // return countSamples;
  assert(num == maxRays);
  return num;
}
//========================================================================




//-----------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// Ward'92 original paper "Measuring and Modeling Anisotropic Reflection".
// PDF is not explicitly shown in this paper, fortunately, paper of
// Kelemen'2001 "A Microfacet Based Coupled Specular-Matte BRDF Model
// with Importance Sampling" contains several PDFs for microfacet
// based BRDF model and for Ward model among them.
//
//  PDF = exp(-tan(delta)^2/slopeX^2) / (slopeX^2 * Pi * cos(delta)^3)
//
//  where slopeX - is the standard deviation of the surface slope, and
//  delta - is the angle between the half-vector and normal
//
// See Documentation for mathematical derivation of sampling scheme
// according to Ward's probability density.
//-----------------------------------------------------------------------------
int
CSchlickSurvey94::SpawnStochasticRays(CContributionArray  &carray,
		                    const CHitPointInfo &hitInfo,
                                    CSamplingInfo       &samplingInfo,
				    EDistributionType distType) const
{
  //cout << "spawnStohasticRays 0\n";
  int maxRays = samplingInfo.GetMaxSamples();
  assert(maxRays > 0);

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

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

  // diffuse component
  CValues kd_local(CValues::EE_Weigth3);
  // specular component
  CValues ks_local(CValues::EE_Weigth3);
  vector<CValues> termsF;

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

  if (_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
    hitInfo.GetTextureValue(kd_local, *_diffuse, CMappingBase::EE_KD);
    akd = kd_local.Luminance();
  }

  if (_specular && (distType & EE_REFLECTED_SPECULAR)) {
    int fresnel_size = _termsFresnel.size();

    if(fresnel_size <= 0)
    {
      WARNING_INFO("CSchlickSurvey94::SpawnStochasticRays() :\n\
Schlick BRDF model must include one or number of slopes Fresnel functions.");
		  return 0;
    }

    for(int i = 0; i < fresnel_size; i++)
    {
      CValues termF(CValues::EE_Weigth3);
      hitInfo.GetTextureValue(termF, *_termsFresnel[i], CMappingBase::EE_KFM);
      termsF.push_back(termF);
    }

    hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
    aks = ks_local.Luminance();//*(n+2.0)/(2.0*M_PI))*

    if(!_slope)
    {
       WARNING_INFO("CSchlickSurvey94::SpawnStochasticRays() :\n\
slope structure should be set to compute specular part of\n\
Schlick model");
      return 0;
    }

    CTexture *slope = _slope->GetSlope();
    if(!slope)
    {
	  WARNING_INFO("CSchlickSurvey94::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope' which consists of at least one 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);*/
  }

  float albedo_whole = akd + aks;
  //cout << "999999999999999999999999999999\n";
  //cout << "akd = " << akd << ", aks = " << aks << ", alb_whole = " << albedo_whole << "\n";

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

  // Compute the number of rays to reflect diffusely / specularly
  int numDiffuse = 0;
  int numSpecular = 0;

  int count = maxRays;
  while(count--)
  {
     float rnd_num = samplingInfo.Sample(0);
     //cout << rnd_num << "\n";
     //cout << "6p\n";

     // This RND number is from the User's program. Actually, it should be
     // fit into the range 0..albedo_whole if some event besides absorbtion
     // was choosen. Therefore, we shoul make RND fall into the [0..1)
     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) {
       // 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));
       //cout << "here\n";
       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++;
     }
     //cout << "rnd = " << rnd_num << ", albedo_whole = " << albedo_whole << "\n";
     //cout << "alpha1 = " << alpha1 << ", alpha2 = " << alpha2 << "\n";
     // For a now it does not work for Halton sequence.
     samplingInfo.IncDimension(2);
     samplingInfo.RecoverSeed();
  }
  //cout << "numDiffuse:  " << numDiffuse  << ", numSpecular: " << numSpecular << endl;

  int numRays = 0;
  //cout << "=============================\n";

  if (numDiffuse) {
    numRays += SpawnStochasticDiffuseRays(kd_local, ks_local, termsF,
    					  carray, hitInfo, samplingInfo,
					  numDiffuse);
  }

  // The number of rays for diffuse interaction
  int indexDiffuse = carray.GetUsedItems();

  if (numSpecular) {
    numRays += SpawnStochasticSpecularRays(ks_local, kd_local, termsF,
					   carray, hitInfo, samplingInfo,
					   numSpecular);
  }

  int indexMax = carray.GetUsedItems();
  assert(indexMax > 0);

  //cout << "indexMax = " << indexMax << ", interaction = " << item->interaction << ", weight = " << item->weigth << "\n";
  //cout << "indexMax = " << indexMax << ", numRefract = " << numRefract << ", numDiffuse = " << numDiffuse << ", numSpecular = " << numSpecular << "\n";

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

  //  DEBUG << "returning " << numRays << endl;
  samplingInfo.IncDimension(2);

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

  if(indexMax < maxRays)
  {
    //cout << "indexMax = " << indexMax << ", numDiffuse = " << numDiffuse << ", numSpecular = " << numSpecular << "\n";
  }
  assert(indexMax == maxRays);
  //cout << "spawnStohasticRays 1\n";
  return numRays;
  //return maxRays;
}

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

const char* 
CSchlickSurvey94::GetHelp() const
{
	return BrdfStrings::KSchlickSurveyDesc;
}
//========================================================================

int 
CSchlickSurvey94::GetNumParameters() const
{
	// for this time, model is considered as isotropic, thus parameter _anisotropic is not used
	return 3 + _slopesCount; // diffuse + specular + fresnel term + number of slopes (1 at least is default)
}
//========================================================================

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

	switch (parameterIndex)
	{
	case 0 :
		_diffuse->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 1 :
		_specular->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	/*case 2 :
		c.Set(0,static_cast<float>(_anisotropic));
		break;*/
	case 2 :
		CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_termsFresnel[0]);
		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-3))
	{
		float weight;
		_slope->GetSlope(parameterIndex-3,weight)->GetValue(c,CHitPointInfo::constMapInfo);
		c.Set(3,weight);
	}

	return c;
}

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

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

	switch (parameterIndex)
	{
	case 0 :
	case 1 :
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 2 :
		_setBasicConstraints(EE_Fresnel,c);
		c.Set(14,5.0f); // default ext coef
		break;
	default :
		_setBasicConstraints(EE_Weighted_Float,c);
		c.Set(0,0); // min value
		c.Set(1,M_PI/2.f); // max value
		c.Set(2,0.5f); // default value
		break;
	}

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

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

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

const char*
CSchlickSurvey94::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseTitle;
	case 1 :
		return BrdfStrings::KSpecularTitle;
	case 2 :
		return BrdfStrings::KFresnelTitle;
	}
	//	if (parameterIndex >= 3)
	return BrdfStrings::KSlopeTitle;
}
//========================================================================

void
CSchlickSurvey94::SetParameter(int parameterIndex, const CValues& value)
{
	switch (parameterIndex)
	{
	case 0 :{
				CConstantTexture* tex = new CConstantTexture(value);
				SetDiffuse(tex);
				break;
			}
	case 1 :{
				CConstantTexture* tex = new CConstantTexture(value);
				SetSpecular(tex);
				break;
			}
	case 2 :{
				CFresnelTexture* tex;
				if (_termsFresnel.size() == 0)
				{
					tex = new CFresnelTexture(CValues(CValues::EE_Weigth3));
					_termsFresnel.push_back(tex);
				}
				else
				{
					tex = dynamic_cast<CFresnelTexture*>(_termsFresnel[0]);
				}
				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 >= 3)
	{
		if (value.GetType() != CValues::EE_Weigth2)
		{
			return;
		}
		float w = value.Get(1);
		int slopeIndex = parameterIndex - 3;
		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
