// ===================================================================
// $Id$
//
// lafortune.cpp
//     Implementation of Lafortune shading model.
//
// Initial coding by Polina Kondratieva (PK), March 2004.
// ===================================================================


// BRDFLIB headers
#include "config.h"
#include "lafortune.h"
#include "rnd.h"
#include "texture.h"
#include <iostream>

__BEGIN_BRDFLIB_SOURCE


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

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

    assert(((_cx.size() == _cy.size()) && (_cx.size() == _cy.size())) &&
            (_cx.size() == _specnums.size()));

    for(int i = 0; i < (int)_cx.size(); i++)
    {
      //if (_specnums[i]) {
        indent(app, ind+2);
        app << "  [" << i << "]: specnum = ";
        _specnums[i]->Describe(app, ind, 1.0/128.0);
        app << "  Cx = ";
        _cx[i]->Describe(app, ind, 1.0/128.0);
        app << "  Cy = ";
        _cy[i]->Describe(app, ind, 1.0/128.0);
        app << "  Cz = ";
        _cz[i]->Describe(app, ind, 1.0/128.0);
      //}
      app << "\n";
    }
  }

  app << flush;

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


// Deep copy of the class entity
//========================================================================
CSurface *
CSurfaceLafortune::Dup() const
{
  assert(((_cx.size() == _cy.size()) && (_cx.size() == _cz.size())) &&
          (_cz.size() == _specnums.size()));
  CSurfaceLafortune *ret = CreateNewObject();

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

  if((int)_cx.size() > 0)
  {
    ret->_specnums = _specnums;
    ret->_cx = _cx;
    ret->_cy = _cy;
    ret->_cz = _cz;
  }

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


/* Albedo for Lafortune reflectance model is different for different
   incoming light directions. Moreover, it differs of ordinary Phong model.*/
//========================================================================
bool
CSurfaceLafortune::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  assert(((_cx.size() == _cy.size()) && (_cx.size() == _cz.size())) &&
          (_cz.size() == _specnums.size()));

  bool retval = false;

  CValues d(CValues::EE_Weigth3);
  CVector3D N = hitInfo.GetNormal();

  if(hitInfo.GetDotProd() > 0.0)
  {
    N = -N;
  }

  float cosine_inc = hitInfo.GetDotProd();

  CValues c(CValues::EE_Weigth3); // result - zero by default
  vector<float> ccx, ccy, ccz, nn;
  for(int i = 0; i < (int)_cx.size(); i++)
  {
     CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
     CValues tmp3(CValues::EE_Weigth1);

     hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
     hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
     hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

     // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
	 // Nekola 14-11-2011 : this is not good..if this condition is not fullfilled
	 // (== lobe is "empty") tmp1/2/3 is not pushed to ccx/y/z vectors, but in
	 // the next iteration, indexing into these ccx/y/z vectors, while checking
	 // correctness (line marked with "MARK") is provided by "i" -> wrong indexing error
	 // proposed correction >> do not check empty lobes
     //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
     //    (tmp3.Get(0) > 1e-06))
      {
        ccx.push_back(tmp1.Get(0));
        ccy.push_back(tmp2.Get(0));
        ccz.push_back(tmp3.Get(0));

	CValues shin;
        hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
        nn.push_back(shin.Get(0));

		// MARK
        float check = pow(fabs(ccx[i]), nn[i]);
        float t1 = pow(fabs(ccy[i]), nn[i]);
        float t2 = pow(fabs(ccz[i]), nn[i]);

        if(check <= t1) check = t1;
        if(check <= t2) check = t2;
        if(check > (nn[i]+2.0)/(2.0*M_PI))
        {
          WARNING_INFO("CSurfaceLafortune::GetAlbedo() :\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		  return 0;
        }
      }
    }

    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      c += d;
      retval = true;
      //cout << "albedo for diffuse term\n";
    }

    if ((distType & _phongReflectionType()) &&
        (_specular && ((int)ccx.size() > 0))) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      CVector3D anisAxis = hitInfo.GetAnisAxis();
      CVector3D dirInit = (hitInfo.GetRay())->GetDir();
      CVector3D initDirInv = -dirInit;  // Incident eye direction in local CS
      //cout << "phong: anisDir = " << anisAxis << ", dirInit = " << dirInit << ", normal = " << N << "\n";

      // Transform incident/exitant directions to the local coordinate system
      CVector3D locAxisX, locAxisY, locAxisZ;
      MakeLocalCS(N, anisAxis, locAxisX, locAxisY, locAxisZ);
      //cout << "lx = " << locAxisX << ", ly = " << locAxisY << ", lz = " << locAxisZ << "\n";

      CVector3D inc_eye, norm_loc;
      // incoming direction in local coordinate system
      inc_eye[0] = DotProd(locAxisX, initDirInv);
      inc_eye[1] = DotProd(locAxisY, initDirInv);
      inc_eye[2] = DotProd(locAxisZ, initDirInv);

      // normal in local coordinate system
      norm_loc[0] = DotProd(locAxisX, N);
      norm_loc[1] = DotProd(locAxisY, N);
      norm_loc[2] = DotProd(locAxisZ, N);

      float cos_lobes = 0;

      for(int i = 0; i < (int)ccx.size(); i++)
      {
        if(ccx[i] > 1e-06 || ccy[i] > 1e-06 || ccz[i] > 1e-06)
	{
	  CVector3D transfIncDir(inc_eye[0]*ccx[i],
                                 inc_eye[1]*ccy[i],
		  	         inc_eye[2]*ccz[i]);
          float transfIncDirLength = Magnitude(transfIncDir);
          transfIncDir.Normalize();

	  float coefficient = 0;
          if(nn[i] > 1e-05) coefficient = pow(transfIncDirLength, nn[i]);
	  else coefficient = M_1_PI;

          float integral = IntegralCosNCos(transfIncDir, norm_loc, nn[i]);
          //float integral = IntegralCosNCos(N, N, nn[i]);
          //cout << "inc_eye = " << inc_eye << ",  transfIncDir = " << transfIncDir << ", length = " << transfIncDirLength <<  "\n";
	  //cout << "ccx = " << ccx[i] << ", ccy = " << ccy[i] << ", ccz = " << ccz[i] << ", nn = " << nn[i] << "\n";
          //cout << "integral = " << integral << ", coef = " << coefficient << ", tmp = " << integral * coefficient << ", lum = " << (d*integral * coefficient).Luminance() << "\n";
	  cos_lobes += integral * coefficient;
	}
      }
      c += d * (cos_lobes/ccx.size());
      //c += d * cos_lobes;
      retval = true;
    }
    //88 cout << "albedo for specular term\t";
    //cout << "incoming = " << initDirInv << "cos_lobes = " << cos_lobes << ", d * cos_lobes = " << d * cos_lobes << "\n";
    //88 cout << ", c = " << c << ", alb_lum = " << c.Luminance() << "\n";
    //else - incoming direction is under the surface
  //}


//#define _DEBUG_PHYSICAL_PLAUSIBILIT

#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";
      }
    }
#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) && (!_refract)) {
      WARNING << "This is the black body - albedo cannot be computed" << endl;
    }
#endif // _DEBUG
  //}
  //else
  //{
  //  cout << "cosine is not appropiate\n";
  //}

  return false;
}
//========================================================================

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

/* Albedo for Lafortune reflectance model is different for different
   incoming light directions. Moreover, it differs of ordinary Phong model.
*/
//========================================================================
bool
CSurfaceLafortune::GetAlbedoMax(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  assert(((_cx.size() == _cy.size()) && (_cx.size() == _cz.size())) &&
          (_cz.size() == _specnums.size()));

  CValues d(CValues::EE_Weigth3);
  CVector3D N = hitInfo.GetNormal();

  if(hitInfo.GetDotProd() > 0.0)
  {
    N = -N;
  }

  bool retval = false;

  // Clean the auxiliary arrays
  //_ccx.clear(); _ccy.clear(); _ccz.clear(); _nn.clear();

  //cout << "~~~~~~~~~~~~~~~~~~~~~~~\n";

  //float cosine_inc = -hitInfo.GetDotProd();
  //cout << "cosine_income = " << cosine_inc << ", inc = " << (hitInfo.GetRay())->GetDir() << ", norm = " << hitInfo.GetNormal() << "\n";
  //if (cosine_inc >= -1e-6)
  //{
    CValues c(CValues::EE_Weigth3); // result - zero by default
    vector<float> ccx, ccy, ccz, nn;
    for(int i = 0; i < (int)_cx.size(); i++)
    {
      CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
      CValues tmp3(CValues::EE_Weigth1);

      hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
      hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
      hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

      // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
      //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
      //    (tmp3.Get(0) > 1e-06))
      {
        ccx.push_back(tmp1.Get(0));
        ccy.push_back(tmp2.Get(0));
        ccz.push_back(tmp3.Get(0));

	CValues shin;
        hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
        nn.push_back(shin.Get(0));

        float check = pow(fabs(ccx[i]), nn[i]);
        float t1 = pow(fabs(ccy[i]), nn[i]);
        float t2 = pow(fabs(ccz[i]), nn[i]);

        if(check <= t1) check = t1;
        if(check <= t2) check = t2;
        if(check > (nn[i]+2.0)/(2.0*M_PI))
        {
			WARNING_INFO("CSurfaceLafortune::GetAlbedoMax():\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		  return 0;
        }
      }
    }

    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      c += d;
      retval = true;
      //cout << "albedo for diffuse term\n";
    }

    // It can happen, that user set ks > 0, but all the coefficients
    // Cx = Cy = Cz = 0. In this case specular term should be ignored
    if ((distType & _phongReflectionType()) &&
        (_specular && ((int)ccx.size() > 0))) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      /*CVector3D anisAxis = hitInfo.GetAnisAxis();
      CVector3D dirInit = (hitInfo.GetRay())->GetDir();
      CVector3D initDirInv = -dirInit;  // Incident eye direction in local CS
      //cout << "phong: anisDir = " << anisAxis << ", dirInit = " << dirInit << ", normal = " << N << "\n";

      // Transform incident/exitant directions to the local coordinate system
      CVector3D locAxisX, locAxisY, locAxisZ;
      MakeLocalCS(N, anisAxis, locAxisX, locAxisY, locAxisZ);
      //cout << "lx = " << locAxisX << ", ly = " << locAxisY << ", lz = " << locAxisZ << "\n";

      CVector3D inc_eye, norm_loc;
      // incoming direction in local coordinate system
      inc_eye[0] = DotProd(locAxisX, initDirInv);
      inc_eye[1] = DotProd(locAxisY, initDirInv);
      inc_eye[2] = DotProd(locAxisZ, initDirInv);

      // normal in local coordinate system
      norm_loc[0] = DotProd(locAxisX, N);
      norm_loc[1] = DotProd(locAxisY, N);
      norm_loc[2] = DotProd(locAxisZ, N);

      float cos_lobes = 0;

      for(int i = 0; i < (int)ccx.size(); i++)
      {
        if(ccx[i] > 1e-06 || ccy[i] > 1e-06 || ccz[i] > 1e-06)
	{
	  CVector3D transfIncDir(inc_eye[0]*ccx[i],
                                 inc_eye[1]*ccy[i],
		  	         inc_eye[2]*ccz[i]);
          float transfIncDirLength = Magnitude(transfIncDir);
          transfIncDir.Normalize();

	  float coefficient = 0;
          if(nn[i] > 1e-05) coefficient = pow(transfIncDirLength, nn[i]);
	  else coefficient = M_1_PI;

          float integral = IntegralCosNCos(transfIncDir, norm_loc, nn[i]);
          //float integral = IntegralCosNCos(N, N, nn[i]);
          //cout << "inc_eye = " << inc_eye << ",  transfIncDir = " << transfIncDir << ", length = " << transfIncDirLength <<  "\n";
	  //cout << "ccx = " << ccx[i] << ", ccy = " << ccy[i] << ", ccz = " << ccz[i] << ", nn = " << nn[i] << "\n";
          //cout << "integral = " << integral << ", coef = " << coefficient << ", tmp = " << integral * coefficient << ", lum = " << (d*integral * coefficient).Luminance() << "\n";
	  cos_lobes += integral * coefficient;
	}
      }
      c += d * (cos_lobes/ccx.size());
      //c += d * cos_lobes;*/
      c += d;
      retval = true;
    }
    //88 cout << "max: albedo for specular term,\t";
    //cout << "incoming = " << initDirInv << "cos_lobes = " << cos_lobes << ", d * cos_lobes = " << d * cos_lobes << "\n";
    //88 cout << "c = " << c << ", alb_lum = " << c.Luminance() << "\n";
    //else - incoming direction is under the surface
  //}


//#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";
      }
    }
#endif // _DEBUG_PHYSICAL_PLAUSIBILITY

    float alb_sum = 0;
    for(int i = 0; i < c.GetCount(); i++)
      alb_sum += c.Get(i);

    alb_sum /= c.GetCount();
    if(alb_sum < 1e-06)
    {
      WARNING_INFO("CSurfaceLafortune::GetAlbedoMax():\n\
Maximum albedo is equal to zero. This surface is a black body.");
      return false;
    }

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

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

  return false;
}
//========================================================================

// For a now it is done in this way and it seems to be reasonable to
// do it in such a way, however, who knows is it correct or not.
//========================================================================
float
CSurfaceLafortune::GetAverageDiffuseness(void) const
{
  assert(((_cx.size() == _cy.size()) && (_cx.size() == _cz.size())) &&
          (_cz.size() == _specnums.size()));

  bool retval = false;
  // mapping info for 2D - for textures - use the value exactly
  // in the middle of the texture
  static CMapPointInfo mapInfo2D(CMappingBase::EE_2D, 0.5);

  vector<float> ccx, ccy, ccz, nn;
  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    if (_cx[i]->GetType()==CMappingBase::EE_None)
      _cx[i]->GetValue(tmp1, CHitPointInfo::constMapInfo);
    else
      _cx[i]->GetValue(tmp1, mapInfo2D);

    if (_cy[i]->GetType()==CMappingBase::EE_None)
      _cy[i]->GetValue(tmp2, CHitPointInfo::constMapInfo);
    else
      _cy[i]->GetValue(tmp2, mapInfo2D);

    if (_cz[i]->GetType()==CMappingBase::EE_None)
      _cz[i]->GetValue(tmp3, CHitPointInfo::constMapInfo);
    else
      _cz[i]->GetValue(tmp3, mapInfo2D);

    // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
    // if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
    //    (tmp3.Get(0) > 1e-06))
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      _specnums[i]->GetValue(shin, mapInfo2D);
      nn.push_back(shin.Get(0));

      float check = pow(fabs(ccx[i]), nn[i]);
      float t1 = pow(fabs(ccy[i]), nn[i]);
      float t2 = pow(fabs(ccz[i]), nn[i]);

      if(check <= t1) check = t1;
      if(check <= t2) check = t2;
      if(check > (nn[i]+2.0)/(2.0*M_PI))
      {
        WARNING_INFO("CSurfaceLafortune::GetAverageDiffuseness():\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		return 0;
      }
    }
  }

  float diffuse=0,specular=0,n=0,refract=0;
  CValues c;


  if (_diffuse) {
    // SIG03 - it is not true
    // assert(_diffuse->GetType()==CMappingBase::EE_None);
    if (_diffuse->GetType()==CMappingBase::EE_None)
      _diffuse->GetValue(c, CHitPointInfo::constMapInfo);
    else
      _diffuse->GetValue(c, mapInfo2D);

    diffuse = c.Luminance();
  }

  if (_specular && ((int)ccx.size() > 0)) {
    // SIG03 - it is not true
    //assert(_specular->GetType()==CMappingBase::EE_None);
    //assert(_specnum->GetType()==CMappingBase::EE_None);

    _specular->GetValue(c, CHitPointInfo::constMapInfo);
    specular = c.Luminance();

    if(specular > 1e-06)
    {
      for(int i = 0; i < (int)ccx.size(); i++)
      {
        n += _nToDiffuseness(nn[i]);
      }
      n /= ccx.size();
    }
  }

  // result is the weighted average
  float res =
    ( diffuse + (specular * n) ) / (diffuse + specular);

  return res;
}
//========================================================================


// For a now it is done in this way and it seems to be reasonable to
// do it in such a way, however, who knows is it correct or not.
//========================================================================
float
CSurfaceLafortune::GetDiffuseness(const CHitPointInfo &hitInfo,
			      EDistributionType distType) const
{
  assert(((_cx.size() == _cy.size()) && (_cx.size() == _cz.size())) &&
          (_cz.size() == _specnums.size()));

  vector<float> ccx, ccy, ccz, nn;
  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
    hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
    hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

    // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
    //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
    //    (tmp3.Get(0) > 1e-06))
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
      nn.push_back(shin.Get(0));

      float check = pow(fabs(ccx[i]), nn[i]);
      float t1 = pow(fabs(ccy[i]), nn[i]);
      float t2 = pow(fabs(ccz[i]), nn[i]);

      if(check <= t1) check = t1;
      if(check <= t2) check = t2;
      if(check > (nn[i]+2.0)/(2.0*M_PI))
      {
        WARNING_INFO("CSurfaceLafortune::GetDiffuseness():\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		return 0;
      }
    }
  }
  float diffuse=0,specular=0,n=0,refract=0;
  CValues c(CValues::EE_Weigth3);

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

  if (((distType & EE_REFLECTED_GLOSSY) && _specular) && ((int)ccx.size() > 0)) {
    hitInfo.GetTextureValue(c, *_specular, CMappingBase::EE_KS);
    specular = c.Luminance();

    if(specular > 1e-06)
    {
      for(int i = 0; i < (int)ccx.size(); i++)
      {
        n += _nToDiffuseness(nn[i]);
      }
      n /= ccx.size();
    }
  }

  float res=
    ( diffuse + (specular * n) ) / (diffuse + specular);
  return res;
}
//========================================================================

//========================================================================
int
CSurfaceLafortune::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  //const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  // we would have to be spectral material to be able to handle
  // more specular indexes
  //assert(shininess.GetCount() == 1);
  if((kslocal[0] < 1e-5 && kslocal[1] < 1e-5) && kslocal[2] < 1e-5)
    return 0;

  //vector<CValues> cx(CValues::EE_Weigth1), cy(CValues::EE_Weigth1);
  //vector<CValues> cz(CValues::EE_Weigth1);
  //assert((_cx.size() > 0 && _cy.size() > 0) &&
  //       (_cz.size() > 0 && _specnums.size() > 0));
  assert(((_cx.size() == _cy.size()) && (_cz.size() == _cy.size())) &&
         (_cx.size() == _specnums.size()));

  vector<float> ccx, ccy, ccz, nn;

  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
    hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
    hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

    //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
    //    (tmp3.Get(0) > 1e-06))
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
      nn.push_back(shin.Get(0));
    }
  }
  // the number of samples generated at the end of this function
  int num = 0;
  //int countSamples = 0;

  //if(ccx <= 1e-05 && ccy <= 1e-05 && ccz <= 1e-05)
  //{
  //  cout << "At least one parameters for Lafortune model should be > 0\n";
  //  return 0;
  //}
  //cout << "specRays: incDir = " << (hitInfo.GetRay())->GetDir() << ", normal = " << hitInfo.GetNormal() << "\n";
  // Choose one of the cosine lobes to sample
  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  // We can reuse specular_rnd[0] for cosine lobe selection, because it is
  // already rescaled to fit into the range [0..1]. However, all these
  // schemes of selection should be revised for Halton random numbers.
  if(ccx.size() > 0)
  {
    double choose_lobe = specular_rnd[0];
    int index_lobe = 0, ccx_size = (int)ccx.size();

    //88 cout << "#cosines = " << ccx_size << ", choose_lobe = " << choose_lobe << "\t";
    //if((int)ccx.size() > 0)
    for(int i = 0; i < ccx_size; i++)
    {
      //88 cout << "i = " << i << ", comp = " << (double)(i+1)/(double)ccx_size << ", ";
      if(choose_lobe <= (double)(i+1)/(double)ccx_size)
      {
        index_lobe = i;
        break;
      }
    }
    //cout << "\nchoose lobe = " << choose_lobe << ", index = " << index_lobe << "\n";
    //88 cout << "index = " << index_lobe << "\n";
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    CVector3D N = hitInfo.GetNormal();
    //if((N[0] > 0.12828 &&  N[0] < 0.12829) && (N[1] > -0.975852 && N[1] < -0.97585) && (N[2] > -0.176798 && N[2] < -0.176796)) cout << "normal = " << N << "\n";
    //CVector3D ttt = N;
    CVector3D anisAxis = hitInfo.GetAnisAxis();
    //cout << ", normal2 = " << N << "\n";
    //cout << "anisDirLaf = " << anisAxis << "\n";
    bool neg_norm = false;
    if (hitInfo.GetDotProd() > 0.0)
    {
      N = -N;
      neg_norm = true;
      //anisAxis[0] = -anisAxis[0];
      //anisAxis[1] = -anisAxis[1];
      //anisAxis[2] = -anisAxis[2];
    }
    // now holds DotProd(N, incomingRay) < 0.0
    CVector3D dirInit = (hitInfo.GetRay())->GetDir();
    CVector3D initDirInv = -dirInit;  // Incident direction

    // diffuseness of the bounce
    //CValues shininess;
    //hitInfo.GetTextureValue(shininess, *nn[index_lobe], CMappingBase::EE_KSI);
    float n = nn[index_lobe];//.Get(0); //shininess.Get(0);

    // Compute diffuseness
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    float diffuseness = 0;
    for(int i = 0; i < (int)ccx.size(); i++)
    {
      diffuseness += _nToDiffuseness(nn[i]);
    }
    diffuseness /= ccx.size();
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    // Lafortune model perform all the computations for incident and
    // exitant directions defined in local coordinate system at the
    // surface.

    // 1. Transform incident direction to the local coordinate system
    //----------------------------------------------------------------------
    // 1.1. Make local coordinate system (locAxisX, locAxisY, locAxisZ)
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector3D locAxisX, locAxisY, locAxisZ;
    MakeLocalCS(N, anisAxis, locAxisX, locAxisY, locAxisZ);
    //MakeLocalCS(N, N, locAxisX, locAxisY, locAxisZ);
    //MakeLocalCS(N, locAxisX, locAxisY, locAxisZ);
    //if(neg_norm) locAxisY = -locAxisY;
    CVector3D inc_eye;

    // 1.2. Calculate vector's coordinates in local coordinate system
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    inc_eye[0] = DotProd(locAxisX, initDirInv);
    inc_eye[1] = DotProd(locAxisY, initDirInv);
    inc_eye[2] = DotProd(locAxisZ, initDirInv);

    // 1.3. Transform incident vector inside the local CS according
    // to the coefficients Cx, Cy, Cz. ('incident' vector now is eye
    // vector. Due to reciprocity we can use computed 'exitant' vector
    // as light-incident vector afterwards.)
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector3D transfIncDir(inc_eye[0]*ccx[index_lobe],
                           inc_eye[1]*ccy[index_lobe],
	  		   inc_eye[2]*ccz[index_lobe]);
    float transfIncDirLength = Magnitude(transfIncDir);
    //cout << "incGlob = " << initDirInv << ", inc_loc(untransf) = " << inc_eye << "\n";
    //cout << "transfIncDir_S = " << transfIncDir << ", n = " << n << ", transfIncDirLen = " << transfIncDirLength << "\t";
    float coefficient = 0;
    if(n > 1e-05) coefficient = pow(transfIncDirLength, n);
    else coefficient = M_1_PI; //transfIncDirLength;
    transfIncDir.Normalize();
    //cout << ", now transfIncDir_S = " << transfIncDir << "\n";

    // 1.4. Calculate transformed incident vector's coords in global CS
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CVector3D incTransfGlobal = transfIncDir.x * locAxisX +
                       transfIncDir.y * locAxisY + transfIncDir.z * locAxisZ;
    //cout << "incTransfGlobal = " << incTransfGlobal << ", length incGlob = " << Magnitude(incTransfGlobal) << "\n";

    // 1.5. Make new coordinate with Z-Axis equal to the reflected
    // vector of the transformed incoming direction
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    const CVector3D &R = 2.0*DotProd(N, initDirInv)*N - initDirInv;
    CVector3D U = ArbitraryNormal(incTransfGlobal); // U is perpendicular to Rr - X-axis
    CVector3D V = CrossProd(incTransfGlobal, U);    // V is perpendicular to Rr and V - Y-axis

    // this value is required to evaluate probability density
    // for outgoing direction
    int num_rejected = 0;
    int num_common = 0;

    for(int i = maxRays; i--; num++) {
      float outCos = -1;
      CVector3D dirRes(0,0,1);
      float cos_ang_Ia_R = 1.0;//, coefficient = 1.0;
      CVector3D local_norm(0,0,1);
      CVector3D ex_light(0,0,1);
      bool rejected = false;

      //CVector3D initDirInv = incident;  // Incident direction in local CS

      // If it is classical PathTracing, only one sample is generated
      // for a given hit point => we should artificially compute the
      // fraction of rejected samples. Samples will be generated until
      // the outCos will became positive
      CVector2D h;

      while(outCos < 0)
      {
        if( !rejected )
          h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
        else
          h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

        //rejected = false;

        // Generate random direction on spherical digon; density is proportional
        // to (cos(alpha))^n; alpha - angle from off-normal (Compendium, 38).
        float sqrt_1_mpow = sqrt(1.0 - pow(h.yy,(2.0f/(n + 1.0f))));
        CVector3D localDir(cos(2.0 * M_PI * h.xx) * sqrt_1_mpow,
                           sin(2.0 * M_PI * h.xx) * sqrt_1_mpow,
                           pow(h.yy, (1.0f/(n + 1.0f))));

        // Change the basis from the VUN coordinate system (local coordinate
        // system where normal Rr is Z-axis) to the World coordinate system
        //
        // Note, that mirror direction Rr is now transformed according
        // to the coefficients Cx, Cy, Cz
        dirRes = localDir.x * U + localDir.y * V + localDir.z * incTransfGlobal;
	//cout << "h = " << h << ", localDir = " << localDir << ", dirRes = " << dirRes << "\n";

        // This checking should be exactly for angle between N, dirRes
        // (see Compendium, p.19)
        outCos = DotProd(N, dirRes);

        if(outCos < 0)
        {
          rejected = true;
        }
      }

      ex_light[0] = DotProd(locAxisX, dirRes);
      ex_light[1] = DotProd(locAxisY, dirRes);
      ex_light[2] = DotProd(locAxisZ, dirRes);

      // Cs(incoming_transformed_local)
      //coefficient = pow(transfIncDirLength, n);

      // angle between the transformed incident direction and outgoing direction
      cos_ang_Ia_R = DotProd(ex_light, transfIncDir);
      //cout << "dirRes = " << dirRes << ", ex_light = " << ex_light << ", cosine = " << cos_ang_Ia_R << "\n";

      // the ray is sampled in this direction
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

      /* Probability density is proportional to Cos^n but not is now normalized:
         pdf = [(n+1)/2/Pi](cos(_alpha))^n/Integral((cos(_alpha))^n, _alpha)

         Importance sampling for Lafortune model is little bit tricky

         BRDF = ps * (Cx*ux*vx + Cy*uy*vy + Cz*uz*vz)^n
            (u - incident direction, v - exitant direction)

         BRDF = ps * ||u^T*M||^n * {u^T*M*v/||u^T*M||}^n
              = ps * Cs(u) * [_u*v]^n = ps * Cs(u) * (cos(_alpha))^n

         where _u = u^T*M*v/||u^T*M|| - transformed and normalized version of
         the incident direction, and _alpha is its angle with v.
           _u = (Cx*ux, Cy*uy, Cz*uz)^T/sqrt{(Cx*ux)^2 + (Cy*uy)^2 + (Cz*uz)^2)

         weight = ps * Cs(u) * (cos(_alpha))^n/ pdf
                = ps * Cs(u) * (cos(_alpha))^n/
	          {((n+1)/2/Pi)*(cos(_alpha))^n/Int((cos(_alpha))^n, _alpha)}
                = ps * Cs(u) * Int((cos(_alpha))^n, _alpha) * (2*Pi)/(n+1);
      */

      // Normalization factor for PDF
      float integral = IntegralCosN(transfIncDir,  local_norm, n);

      // Because the result for cosine lobes should be scaled by their number in
      // order to save plausibility, this chosen cosine lobe also should be scaled
      // by this number
      contrib->weigth = kslocal*(1.0/ccx_size) * coefficient * integral * ((2.0*M_PI)/(n+1.0));
      contrib->prob = ((n+1.0)/(2.0*M_PI))*pow(cos_ang_Ia_R, n)/integral;
      contrib->interaction = EE_REFLECTED_SPECULAR;
      contrib->diffuseness = diffuseness;
      //cout << "kslocal = " << kslocal << ", coef = " << coefficient << ", pow = " << pow(cos_ang_Ia_R, n) << ", brdf = " << kslocal * coefficient * pow(cos_ang_Ia_R, n) << "\n";
      //cout << "weight = " << contrib->weigth << ", prob = " << contrib->prob << "\n";

      // Add remaining cosine lobes (to supplement BRDF function)
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      float cos_lobes = 0;

      for(int i = 0; i < ccx_size; i++)
      {
        //cout << "i = " << i << ", index_lobe now = " << index_lobe << ", ";
        if(i != index_lobe)
        {
	  float tmp = ccx[i]*inc_eye[0]*ex_light[0] +
	              ccy[i]*inc_eye[1]*ex_light[1] +
	              ccz[i]*inc_eye[2]*ex_light[2];

          // If n~0 we should provide coefficient 1/M_PI to keep physical
          // plausibility (because for n = 0 Lafortune -> diffuse).
          if(nn[i] < 1e-06)
	  {
            cos_lobes += M_1_PI;
  	  }
	  else
	  {
	    if(tmp > 1e-05) cos_lobes += pow(tmp, nn[i]);
	  }
          //cout << "tmp = " << tmp << ", nn = " << nn[i] << ", cos_lobes = " << cos_lobes << "\n";
        }
      }
      contrib->weigth += kslocal * (cos_lobes/(contrib->prob*ccx_size));
      if ( _maxLafortune){
        //max Lafortune model [13]
        float incEyeCos  = fabs(hitInfo.GetDotProd());
        float lightCos = fabs(DotProd(hitInfo.GetNormal(), dirRes));
        float max_cos = Max(incEyeCos, lightCos);
        contrib->weigth *= (1.0/max_cos);
      }
      //contrib->weigth += kslocal * (cos_lobes/contrib->prob);
      //88 cout << "now cos_lobes (spec rays)= " << cos_lobes << ", weight = " << contrib->weigth << "\t";

      // Add diffuse term (to supplement BRDF function)
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CValues kd_local(CValues::EE_Weigth3);
      if (_diffuse) {
        hitInfo.GetTextureValue(kd_local, *_diffuse, CMappingBase::EE_KD);

        if(kd_local.Luminance() > 1e-05) {
          kd_local *= 1.0/(M_PI*contrib->prob);
          contrib->weigth += kd_local;
        }
      }
      //88 cout << "diffuse = " << kd_local << ", after diffuse weight = " << contrib->weigth << "\n";

      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


      // set the reflected ray, direction is normalized here
      contrib->weigth *= outCos;
      //cout << "outCos = " << outCos << ", after cosine = " << contrib->weigth << "\n";
      //cout << "probability = " << contrib->prob << ", weight = " << contrib->weigth << "\n";
      assert(Magnitude(dirRes) < 1.0+1e-6 && Magnitude(dirRes) > 0);
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      //cout << "dirRes = " << dirRes << ", weight = " << contrib->weigth << "\n";

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

      // the number of valid samples taken
      //countSamples++;
  } // for i
}
 // return countSamples;
  return num;
}
//========================================================================


/*
  kd_local is albedo => it should be divided by PI to be used
  as the value of BRDF.
*/
//========================================================================
int
CSurfaceLafortune::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  //const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  CVector3D initDirInv = -dirInit;  // Incident direction in local CS
  CVector3D N = hitInfo.GetNormal();
  CVector3D anisAxis = hitInfo.GetAnisAxis();

  vector<float> ccx, ccy, ccz, nn;
  //assert((_cx.size() > 0 && _cy.size() > 0) &&
  //       (_cz.size() > 0 && _specnums.size() > 0));
  assert(((_cx.size() == _cy.size()) && (_cz.size() == _cy.size())) &&
         (_cx.size() == _specnums.size()));

  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
    hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
    hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

    //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
    //   (tmp3.Get(0) > 1e-06))
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
      nn.push_back(shin.Get(0));
    }
  }

  if (hitInfo.GetDotProd() > 0.0)
    N = -N;

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

  //float n = shininess.Get(0);
  //88 cout << "========================================\n";
  // number of rays created
  int num = 0;

  for(int i = maxRays; i--; num++) {
    //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 * 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, lit.dir) / pdf =
    //         = pd * (n+2)/(2*Pi) * cos(normal, lit.dir) * (incZ*outZ)^n /
    //           (cos(normal, ray.dir)/Pi)
    //         = pd*((n+2)/2)*(incZ*outZ)^n
    //
    // 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);
    // This gives strange result -> to much light

    contrib->weigth = kd_local; // * M_1_PI;

    // Add specular term
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // Transform incident/exitant directions to the local coordinate system

    //CVector3D transfLightDir(ex_light[0]*_cx, ex_light[1]*_cy, ex_light[2]*_cz);
    //CVector3D transfLightDir(inc_eye[0]*ccx, inc_eye[1]*ccy, inc_eye[2]*ccz);
    //float transfLightDirLength = Magnitude(transfLightDir);
    //transfLightDir.Normalize();

    //CValues specularIndex(CValues::EE_Weigth3);
    //if (_specnum)
    //  hitInfo.GetTextureValue(specularIndex, *_specnum, CMappingBase::EE_KSI);
    //else
    // supposing that if not specified, surface is diffuse
    //specularIndex.Set(1.0);

    //----------------------------------------------------------------------
    CValues specularTerm(CValues::EE_Weigth3);
    specularTerm = CValues(0,0,0);

    // if surface has specular component
    if (_specular && ((int)ccx.size() > 0)) {
      CVector3D locAxisX, locAxisY, locAxisZ;
      MakeLocalCS(N, anisAxis, locAxisX, locAxisY, locAxisZ);

      CVector3D inc_eye(0,0,1), ex_light(0,0,1);

      inc_eye[0] = DotProd(locAxisX, initDirInv);
      inc_eye[1] = DotProd(locAxisY, initDirInv);
      inc_eye[2] = DotProd(locAxisZ, initDirInv);

      ex_light[0] = DotProd(locAxisX, dirRes);
      ex_light[1] = DotProd(locAxisY, dirRes);
      ex_light[2] = DotProd(locAxisZ, dirRes);

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

      if(specularTerm.Luminance() > 1e-05)
      {
        // the cosine angle between the ideally reflected direction and
        // the input direction
        float outCos = DotProd(N, dirRes);  // Light source is under surface

        // It can happen that angle(transfLightDir, ex_light) become more
	// than Pi/2 => sample should be rejected
	//float dott = DotProd(transfLightDir, ex_light);
        //(dott >= 0) ? dott = dott : dott = 0;

        if (outCos > 0.0) {
          float cos_lobes = 0;
          for(int i = 0; i < (int)ccx.size(); i++)
          {
            float tmp = ccx[i]*inc_eye[0]*ex_light[0] +
	                ccy[i]*inc_eye[1]*ex_light[1] +
	                ccz[i]*inc_eye[2]*ex_light[2];

            // If n~0 we should provide coefficient 1/M_PI to keep physical
	    // plausibility (because for n = 0 Lafortune -> diffuse).
            if(nn[i] < 1e-06)
	    {
              cos_lobes += M_1_PI;
	    }
	    else
	    {
	      if(tmp > 1e-06) cos_lobes += pow(tmp, nn[i]);
	    }
          }
          //specularTerm *= cos_lobes/contrib->prob;
          if ( _maxLafortune){
            //max Lafortune model [13]
            float incEyeCos  = fabs(hitInfo.GetDotProd());
            float lightCos = fabs(DotProd(hitInfo.GetNormal(), dirRes));
            float max_cos = Max(incEyeCos, lightCos);
             contrib->weigth += (1.0/max_cos) * specularTerm * outCos * (cos_lobes/(contrib->prob*ccx.size()));
          } else
            contrib->weigth += specularTerm * outCos * (cos_lobes/(contrib->prob*ccx.size()));
          //contrib->weigth += specularTerm * outCos * (cos_lobes/contrib->prob);
	}
      }
    }
    //88 cout << "contrib diffuse rays = " << contrib->weigth << "\n";

    //contrib->weigth = kd_local;

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

  } // for i

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


//    This function set weight = BRDF/probability
// where probability includes PDF and probability of choosing the
// BRDF model.
//    This multiplication done to reduce the time of computations.
// PDF is usually proportional to the BRDF, so it is possible
// to use this coeficient of proportionality as weight instead
// of compute both BRDF and PDF and then make division. 24.11.03 P.K.
//========================================================================
int
CSurfaceLafortune::SpawnStochasticRays(
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo,
  CSamplingInfo        &samplingInfo,
  EDistributionType     distType) const
{
  //cout << "1p\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.
  //refraction_rnd.clear();
  diffuse_rnd.clear();
  specular_rnd.clear();
  //!!!!!!!!!!!!!!!!!!!!!!!! K.P.!!!!!!!!!!!!!!!!!!!!!!!!!!

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

  assert(((_cx.size() == _cy.size()) && (_cz.size() == _cy.size())) &&
         (_cx.size() == _specnums.size()));

  //_ccx.clear(); _ccy.clear(); _ccz.clear(); _nn.clear();

  //88 cout << "SStR: _cx.size() = " << _cx.size() << "\n";
  vector<float> ccx, ccy, ccz, nn;
  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
    hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
    hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

    // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
    //if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
    //    (tmp3.Get(0) > 1e-06))
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
      nn.push_back(shin.Get(0));

      float check = pow(fabs(ccx[i]), nn[i]);
      float t1 =    pow(fabs(ccy[i]), nn[i]);
      float t2 =    pow(fabs(ccz[i]), nn[i]);

      if(check <= t1) check = t1;
      if(check <= t2) check = t2;
      if(check > (nn[i]+2.0)/(2.0*M_PI))
      {
        WARNING_INFO("CSurfaceLafortune::SpawnStochasticRays():\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		return 0;
      }
    }
  }

  // diffuse component
  CValues kd_local(CValues::EE_Weigth3);
  // specular component
  CValues ks_local(CValues::EE_Weigth3), ks_inc_local(CValues::EE_Weigth3);
  // refractive component (ideal refraction)
  //CValues kr_local(CValues::EE_Weigth3);

  // compute diffuse, specular, and refractive albedos
  float akd = 0.0f, aks = 0.0f, aks_inc = 0.0f; //, akr = 0.0f;
  EDistributionType tmp_distType;

  hitInfo.GetNormal(); // << "\n";

  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)ccx.size() > 0)) {
    hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
    tmp_distType = EE_REFLECTED_GLOSSY;
    GetAlbedo(ks_inc_local, hitInfo, tmp_distType);
    aks = ks_local.Luminance();
    aks_inc = ks_inc_local.Luminance();
  }

  //if (_refract && (distType & EE_REFRACTED_PERFECT)) {
  //   hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
  //   akr = kr_local.Luminance();
  // }

  //float albedo_whole = aks + akd; // + akr;
  float albedo_whole = aks_inc + akd; // + akr;
  //cout << "albedo_whole = " << albedo_whole << "\n";

  if(albedo_whole < 1.e-05) {
	WARNING_INFO("Albedo of this surface is almost zero => it is a black body!");
    return false;
  }
  
  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CSurfaceLafortune::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 refract / reflect diffusely / specularly
  int numRefract = 0;
  int numDiffuse = 0;
  int numSpecular = 0;

  /* This is the rule for choosing the interaction:
   *
   *  +----------------------------------------------------+
   *  |   refract   |  reflect diffusely |    specularly   |
   *  +----------------------------------------------------+
   *              (s1)                  (s2)              (s3)
   */

  int count = maxRays;
  while(count--)
  {
     float rnd_num = samplingInfo.Sample(0);
     //cout << "before scaling rnd_num: " << rnd_num << "\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;
     }

    /* This is the rule for choosing the interaction:
     *
     *  +----------------------------------------------------+
     *  |   refract   |  reflect diffusely |    specularly   |
     *  +----------------------------------------------------+
     *  |<-- aplha1 ->|<----- alpha2 ----->|<---- alpha3 --->|
     */
     float alpha1 = akd/albedo_whole;
     //float alpha2 = aks/albedo_whole;
     float alpha2 = aks_inc/albedo_whole;
     //float alpha3 = aks/albedo_whole;
     //cout << "alpha1 = " << alpha1 << ", alpha2 = " << alpha2 << ", cx_size = " << ccx.size() << ", rnd = " << rnd_num << "\n";

     /*if(rnd_num < alpha1 && alpha1 > 1e-5) {
       // Reuse random number for generation of one component of the random
       // direction. Scale it to fall into the range [0..1)
       refraction_rnd.push_back(rnd_num/alpha1);
       refraction_rnd.push_back(samplingInfo.Sample(1));
       numRefract++;
     }
     else {*/
       if(rnd_num < alpha1 && alpha1 > 1e-5) {
         // 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 {
         //if((alpha2 > 1e-5) && ((int)ccx.size() > 0))
	 //{
	    // 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));
            //cout << "specular_rnd[0] = " << specular_rnd[0] << ", specular_rnd[1] = " << specular_rnd[1] << "\n";
	    numSpecular++;

	 //}/
       }
     //}
     // For a now it does not work for Halton sequence.
     samplingInfo.IncDimension(2);
     samplingInfo.RecoverSeed();
     //cout << "7p\n";
  }
  //cout << "spec = " << numSpecular << ", dif = " << numDiffuse << "\n";

  int numRays = 0;

  //if (numRefract)
  //   numRays += SpawnIdealRefractedRay(kr_local, carray, hitInfo);

  // If there is refraction indexRefracted = 1
  //int indexRefracted = carray.GetUsedItems();

  //CValues shininess;
  //hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
  //cout << "getCount = " << shininess.GetCount() << "\n";

  if (numDiffuse) {
    numRays += SpawnStochasticDiffuseRays(kd_local, carray,
					  hitInfo, samplingInfo,
					  numDiffuse);
  }
  // The number of rays for diffuse interaction
  int indexDiffuse = carray.GetUsedItems(); // - indexRefracted;

  if (numSpecular) {
    numRays += SpawnStochasticSpecularRays(ks_local, carray,
					   hitInfo, samplingInfo,
					   numSpecular);
            //88 cout << "sssssssssssssssssssssssssssppppppppppppppppppppp\n";
  }
  //cout << "numRefract = " << numRefract << ", numDiffuse = " << numDiffuse << ", numSpecular = " << numSpecular << "\n";
  //cout << "really = " << numRays << "\n";

  int indexMax = carray.GetUsedItems();

  for(int i = 0; i < indexMax; i++) {
    CContributionItem* ci = carray.GetItem(i);
    /*float mult;
    if(i < indexDiffuse) {
      mult = akd;
    }
    else {
      mult = aks;
    }
    mult /= albedo_whole;
    //ci->prob *= mult;
    //ci->weigth *= (1.0/(mult * (float)maxRays));*/
    ci->weigth *= (1.0/(float)maxRays);
  } // for

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

  // recover the seed for the next dimension to be sampled
  samplingInfo.RecoverSeed();
  //cout << "2\n";

  // Be careful with this number of rays!!!!! Initial Phong model
  // have no possibility to add the transmitivity to the object.
  // We suppose to use deterministic refractive rays. But this causes
  // some difficultiues: only one ray will be produced for given hit
  // point does not matter which part of the light is refracted.
  // There should be shot maxRays number of rays, but it will be shot
  // less number of rays because of it is not necessary to shoot the
  // same ray several times - it is just wasting of time.
  // Therefore, if the result sum of all rays' contributions for a
  // given pixel will be divided by maxRays the pixel color will be
  // dimmer than it should be => the sum of all rays' contributions for
  // a given pixel should be divided by numRays -> actual number of rays.
  // 06/12/03/ K.P.

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


//========================================================================
bool
CSurfaceLafortune::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &dirOut,
				    EDistributionType distType) const
{
  CValues res(bsdf.GetType());
  bool retval = false;
  CVector3D normal = hitInfo.GetNormal();
  CVector3D anisAxis = hitInfo.GetAnisAxis();
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  CVector3D initDirInv = -dirInit;  // Incident eye direction in local CS

  vector<float> ccx, ccy, ccz, nn;
  //assert((_cx.size() > 0 && _cy.size() > 0) && _cz.size() > 0);
  assert(((_cx.size() == _cy.size()) && (_cz.size() == _cy.size())) &&
         (_cx.size() == _specnums.size()));

  for(int i = 0; i < (int)_cx.size(); i++)
  {
    CValues tmp1(CValues::EE_Weigth1), tmp2(CValues::EE_Weigth1);
    CValues tmp3(CValues::EE_Weigth1);

    hitInfo.GetTextureValue(tmp1, *_cx[i], CMappingBase::EE_CX);
    hitInfo.GetTextureValue(tmp2, *_cy[i], CMappingBase::EE_CY);
    hitInfo.GetTextureValue(tmp3, *_cz[i], CMappingBase::EE_CZ);

    // check for 'empty' cosine lobes (for them Cx, Cy, Cz are zeros)
    /*if(((tmp1.Get(0) > 1e-06) || (tmp2.Get(0) > 1e-06)) ||
        (tmp3.Get(0) > 1e-06))*/
    {
      ccx.push_back(tmp1.Get(0));
      ccy.push_back(tmp2.Get(0));
      ccz.push_back(tmp3.Get(0));

      CValues shin;
      hitInfo.GetTextureValue(shin, *_specnums[i], CMappingBase::EE_KSI);
      nn.push_back(shin.Get(0));

      float check = pow(fabs(ccx[i]), nn[i]);
      float t1 = pow(fabs(ccy[i]), nn[i]);
      float t2 = pow(fabs(ccz[i]), nn[i]);

      if(check <= t1) check = t1;
      if(check <= t2) check = t2;
      if(check > (nn[i]+2.0)/(2.0*M_PI))
      {
        WARNING_INFO("CSurfaceLafortune::EvaluateBSDF():\n\
For each cosine lobe absolute value of coefficients Cx, Cy, Cz should satisfy the condition :\n coefficient^shininess <= (shininess + 2) / (2*Pi)");
		return false;
      }
    }
  }

  // Transform incident/exitant directions to the local coordinate system
  CVector3D locAxisX, locAxisY, locAxisZ;
  MakeLocalCS(normal, anisAxis, locAxisX, locAxisY, locAxisZ);

  CVector3D inc_eye, ex_light, nl;

  inc_eye[0] = DotProd(locAxisX, initDirInv);
  inc_eye[1] = DotProd(locAxisY, initDirInv);
  inc_eye[2] = DotProd(locAxisZ, initDirInv);

  ex_light[0] = DotProd(locAxisX, dirOut);
  ex_light[1] = DotProd(locAxisY, dirOut);
  ex_light[2] = DotProd(locAxisZ, dirOut);

  nl[0] = DotProd(locAxisX, normal);
  nl[1] = DotProd(locAxisY, normal);
  nl[2] = DotProd(locAxisZ, normal);
  //cout << "====================\n";
  //cout << "anisAxis = " << anisAxis << ", locAxisZ = " << locAxisZ << ", locAxisX = " << locAxisX <<", locAxisY = " << locAxisY <<"\n";
  //cout << "normal = " << normal << ", loc_norm = " << nl << ", inc_eye = " << initDirInv << ", loc_eye = " << inc_eye << "\n";
  //cout << "ex_lit = " << dirOut << ", loc_lit = " << ex_light << "\n";

  // evaluating diffuse term
  //----------------------------------------------------------------------
  if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
    // We don't check the angle between normal and outgoing direction.
    // Maybe we should (not to illuminate backfaces), but in that case
    // the dot product must be evaluated here and in the renderer as well.
    CValues diffuseTerm(res.GetType());
    hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);

    if(diffuseTerm.Luminance() > 1e-06)
    {
      res += diffuseTerm * M_1_PI;
      retval = true;
    }
  }

  // evaluating specular term
  //----------------------------------------------------------------------
  CValues specularTerm(res.GetType());
  float ratio = DotProd(normal, dirOut);  // cos(light,normal)

  // For a now shininess is single-component color (non-RGB)
  // if surface has specular component
  if (((distType & EE_REFLECTED_GLOSSY) && _specular) && ((int)ccx.size() > 0))
  {
    if (ratio > 0.0) {
      // get the value of specular coefficient
      hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KS);

      if(specularTerm.Luminance() > 1e-06)
      {
        //cout << "-------------------\n";
	//cout << "inc_eye = " << inc_eye << ", ex_light = " << ex_light << "\n";
        float cos_lobes = 0;
        for(int i = 0; i < (int)ccx.size(); i++)
        {

	  float tmp = ccx[i]*inc_eye[0]*ex_light[0] +
	              ccy[i]*inc_eye[1]*ex_light[1] +
	              ccz[i]*inc_eye[2]*ex_light[2];

	  //cout << "ccx = " << ccx[i] << ", ccy = " << ccy[i] << ", ccz = " << ccz[i] << ", n = " << nn[i] << "\n";
	  // If n~0 we should provide coefficient 1/M_PI to keep physical
	  // plausibility (because for n = 0 Lafortune -> diffuse).
	  if(nn[i] < 1e-6)
	  {
            cos_lobes += M_1_PI;
	  }
	  else
	  {
	    if(tmp > 1e-05) cos_lobes += pow(tmp, nn[i]);
	  }
          //cout << "tmp = " << tmp << ", cos_lobes = " << cos_lobes << "\n";
        }
        specularTerm *= cos_lobes/ccx.size(); // * coefficient;
        //specularTerm *= cos_lobes; // * coefficient;

        res += specularTerm;
        retval = true;
      }
    }
  } // specular term

  // 18-10-2011, NEKOLA : ???? 
  // at 1.) : this will not influence result in any way, because specular term is already included in res 
  // at 2.) : if maxlafortune is true, result should not be multiplied once again with non-max lafortune (by ratio) (else statement missing)
  /*if (_maxLafortune){
    float incEyeCos   = -(hitInfo.GetDotProd());
    float lightCos    = DotProd(normal, dirOut);
    float max_cos = Max(fabs(incEyeCos),fabs(lightCos));
    specularTerm *= (1.0f/max_cos);
  }*/

  // corrected version
  if (_maxLafortune){
    float incEyeCos   = -(hitInfo.GetDotProd());
	// light cos is already in ratio variable
    float max_cos = Max(fabs(incEyeCos),fabs(ratio));
    res *= 1.0f/max_cos;
  }

  if (retval) {
    bsdf = res * ratio;
    return true;
  }
  return false;
}
//========================================================================

void CSurfaceLafortune::AddLobe(CTexture *cx, CTexture *cy, CTexture *cz, CTexture *specnum)
{
	_cx.push_back(cx);
	_cy.push_back(cy);
	_cz.push_back(cz);
	_specnums.push_back(specnum);
	_currentLobesCount++;
}
//========================================================================

const char* 
CSurfaceLafortune::GetHelp() const
{
	return BrdfStrings::KLafortuneDesc;
}
//========================================================================

int 
CSurfaceLafortune::GetNumParameters() const
{
	return (3 + _currentLobesCount); // diffuse, specular (from Phong parent)
									 // + maxLafortune
									 // + lobes count (x,y,z,specnum = 1 weighted xyz texture) 	
}
//========================================================================

CValues
CSurfaceLafortune::GetParameter(int parameterIndex) const
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
	{
		return CSurfacePhongBase::GetParameter(parameterIndex);
	}
	else if (parameterIndex == 2)
	{
		CValues val(0,0,0);
		val.Set(_maxLafortune);
		return val;
	}
	else if (parameterIndex > 0 && parameterIndex < (_currentLobesCount+3))// lobes
	{
		CValues auxC(0,0,0);
		CValues c(CValues::EE_WeigthN);
		if ((parameterIndex-3) >= _cx.size())
			return CValues(0);
		_cx[parameterIndex-3]->GetValue(auxC,CHitPointInfo::constMapInfo); 
		c.Set(0,auxC.Get(0)); //cx of the lobe
		_cy[parameterIndex-3]->GetValue(auxC,CHitPointInfo::constMapInfo); 
		c.Set(1,auxC.Get(0)); //cy of the lobe
		_cz[parameterIndex-3]->GetValue(auxC,CHitPointInfo::constMapInfo); 
		c.Set(2,auxC.Get(0)); //cz of the lobe
		_specnums[parameterIndex-3]->GetValue(auxC,CHitPointInfo::constMapInfo); 
		c.Set(3,auxC.Get(0)); //specnum of the lobe
		return c;
	}
	// Problem with an index
	CValues c(3);
	c.Set(0.f);
	return c;	
}
//========================================================================

CValues
CSurfaceLafortune::GetParameterConstraints(int parameterIndex) const
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
		return CSurfacePhongBase::GetParameterConstraints(parameterIndex);
	else if (parameterIndex == 2)
	{
		CValues c(3);
		_setBasicConstraints(EE_Boolean,c);
		return c;
	}
	else // lobes
	{
		CValues c(12);
		c.Set(0,-1); //cx of the lobe - min
		c.Set(1,1);  //cx of the lobe - max
		c.Set(2,0);  //cx of the lobe - default
		c.Set(3,-1); //cy of the lobe - min
		c.Set(4,1);  //cy of the lobe - max
		c.Set(5,0);  //cy of the lobe - default
		c.Set(6,-1); //cz of the lobe - min
		c.Set(7,1);  //cz of the lobe - max
		c.Set(8,0.001);  //cz of the lobe - default - lobe shouldnt direct under the surface
		c.Set(9,0.001); //specnum of the lobe - min
		c.Set(10,MAX_SHININESS); //specnum of the lobe - max
		c.Set(11,10); //specnum of the lobe - default
		return c;
	}
	// Problem with an index
	CValues c(3);
	c.Set(0.f);
	return c;
	
}
//========================================================================

const char* 
CSurfaceLafortune::GetParameterDescription(int parameterIndex) const
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
		return CSurfacePhongBase::GetParameterDescription(parameterIndex);
	if (parameterIndex == 2)
		return BrdfStrings::KMaxLafDesc;
	if (parameterIndex > 2)
		return BrdfStrings::KLobeDesc;
	else return BrdfStrings::KInvalidParameter;
}
//========================================================================

CSurface::EParameterUIType
CSurfaceLafortune::GetParameterSyntax(int parameterIndex) const
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
		return CSurfacePhongBase::GetParameterSyntax(parameterIndex);
	if (parameterIndex == 2)
		return EE_Boolean;
	if (parameterIndex > 2)
		return EE_XYZW;
	else return EE_Unknown;
}
//========================================================================

const char*
CSurfaceLafortune::GetParameterTitle(int parameterIndex) const
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
		return CSurfacePhongBase::GetParameterTitle(parameterIndex);
	if (parameterIndex == 2)
		return BrdfStrings::KMaxLafTitle;
	if (parameterIndex > 2)
		return BrdfStrings::KLobeTitle;
	else return BrdfStrings::KInvalidParameter;
}
//========================================================================

void
CSurfaceLafortune::SetParameter(int parameterIndex, const CValues& value)
{
	if ((parameterIndex>=0) && (parameterIndex<2)) // parent parameter
	{
		CSurfacePhongBase::SetParameter(parameterIndex,value);
		return;
	}
	if (parameterIndex == 2) 
	{
		bool boolValue = static_cast<bool>(value[0]);
		SetMaxLafortune(boolValue);
	}
	if (parameterIndex > 2)
	{
		if (value.GetType() != CValues::EE_Weight4)
		{
			return;
		}
		CValues auxC(CValues::EE_Weigth1);
		auxC.Set(value.Get(0));
		CConstantTexture* newTexX = new CConstantTexture(auxC);
		auxC.Set(value.Get(1));
		CConstantTexture* newTexY = new CConstantTexture(auxC);
		auxC.Set(value.Get(2));
		CConstantTexture* newTexZ = new CConstantTexture(auxC);
		auxC.Set(value.Get(3));
		CConstantTexture* newTexSpecnum = new CConstantTexture(auxC);
		int lobeIndex = parameterIndex - 3;
		if ((unsigned int)(lobeIndex) < _cx.size() )
		{
			delete _cx[lobeIndex];
			_cx[lobeIndex] = newTexX;
			delete _cy[lobeIndex];
			_cy[lobeIndex] = newTexY;
			delete _cz[lobeIndex];
			_cz[lobeIndex] = newTexZ;
			delete _specnums[lobeIndex];
			_specnums[lobeIndex] = newTexSpecnum;
		}
		else // new lobe
		{
			_cx.push_back(newTexX);
			_cy.push_back(newTexY);
			_cz.push_back(newTexZ);
			_specnums.push_back(newTexSpecnum);
		}
	}
}
//========================================================================





__END_BRDFLIB_SOURCE

