// ===================================================================
// $Id: phong.cpp,v 1.22 2003/11/04 15:27:18 havran Exp $
//
// phong.cpp
//     Implementation of Phong shading model.
//
// Licence and usage conditions in November 2003:
// This file can be used and modified only with written permission of
// Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
// This library is granted to be used by Polina Kondratieva within
// her MPI HIWI job working on BRDFLIB library. The usage condition
// may change in the next version of the library. At the moment,
// the library development is supported by Max-Planck-Institut for
// Informatik, Stuhlsatzenhausweg 85, 66123 Saarbruecken Germany
// See WWW: http://www.mpi-sb.mpg.de
// REPLACEMENT_STRING
//
// Initial coding by Nicholas Wilt,
// rewritten by Vlastimil Havran & Jaroslav Krivanek, 2000.
// Changed by KP, December 2003
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


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

__BEGIN_BRDFLIB_SOURCE

// This variable is used to store intermediate result of random number
// generation if Selection Roulette is used.
double
CSurfacePhongBase::aux_rnd = 0.0;
typedef vector<double> rnd_array;
rnd_array CSurfacePhongBase::diffuse_rnd;
rnd_array CSurfacePhongBase::specular_rnd;
rnd_array CSurfacePhongBase::refraction_rnd;

/**********************************************************************
 *
 *             CSurfacePhongBase methods
 *
 **********************************************************************/

// describes the surface to a given stream
//========================================================================
void
CSurfacePhongBase::Describe(ostream &app, int ind) const
{
  indent(app, ind);
  app << "CSurfacePhong:\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";
  }
  if (_specnum) {
    indent(app, ind+2);
    app << "  specnum = \t";
    _specnum->Describe(app, ind, 1.0/128.0); app << "\n";
  }
  if (_refract) {
    indent(app, ind+2);
    app << "  refract = \t";
    _refract->Describe(app, ind+4); app << "\n";
  }
  if (_index) {
    indent(app, ind+2);
    app << "  index = \t";
    _index->Describe(app, ind+4); app << "\n";
  }
  app << flush;

  return;
}


const char* 
CSurfacePhongBase::GetHelp() const
{
	return BrdfStrings::KPhongPhysicalDesc;
}

int 
CSurfacePhongBase::GetNumParameters() const
{
	return 9; // diffuse,specular,specnum,refract,index,extinction,maxPhong,normalsBlinn,retroReflective
}

CValues
CSurfacePhongBase::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	switch (parameterIndex)
	{
	case 0 : //diffuse
		_diffuse->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 1 : //specular
		_specular->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 2 : //specnum
		_specnum->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 3 : //refract
		_refract->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 4 : //index
		_index->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 5 : //extinction
		_extinction->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 6 : //maxPhong
		c.Set(_maxPhong);
		break;
	case 7 : //normalsBlinn
		c.Set(_normalsBlinn);
		break;
	case 8 : //retroReflective
		c.Set(_retroReflective);
		break;
	default :
		break;
	}
	return c;
}

CValues
CSurfacePhongBase::GetParameterConstraints(int parameterIndex) const
{
	CValues c(3);
	switch (parameterIndex)
	{
	case 0 : //diffuse
	case 1 : //specular
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 2 : //specnum
		_setBasicConstraints(EE_RGB_Components,c);
		break;
	case 3 : //refract
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 4 : //index
	case 5 : //extinction
		c.Set(0,1); // min
		c.Set(1,MAX_IOR); // max
		c.Set(2,1); // default
		break;
	case 6 : //maxPhong
	case 7 : //normalsBlinn
	case 8 : //retroReflective
		_setBasicConstraints(EE_Boolean,c);
		break;
	default :
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	}
	return c;
}

const char* 
CSurfacePhongBase::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseDesc;
	case 1 :
		return BrdfStrings::KSpecularDesc;
	case 2 :
		return BrdfStrings::KShininessDesc;
	case 3 :
		return BrdfStrings::KRefractionDesc;
	case 4 :
		return BrdfStrings::KIORDesc;
	case 5 :
		return BrdfStrings::KExtinctionCoefDesc;
	case 6 :
		return BrdfStrings::KMaxPhongDesc;
	case 7 :
		return BrdfStrings::KBlinnNormalsDesc;
	case 8 :
		return BrdfStrings::KRetroReflectiveDesc;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}

CSurface::EParameterUIType
CSurfacePhongBase::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 : //diffuse
		return EE_ColorRGB;
	case 1 : //specular
		return EE_ColorRGB;
	case 2 : //specnum
		return EE_RGB_Components;
	case 3 : //refract
		return EE_ColorRGB;
	case 4 : //index
		return EE_XYZ;
	case 5 : //extinction
		return EE_XYZ;
	case 6 : //maxPhong
		return EE_Boolean;
	case 7 : //normalsBlinn
		return EE_Boolean;
	case 8 : //retroReflective
		return EE_Boolean;
	default :
		return EE_Unknown;
	}
}

const char*
CSurfacePhongBase::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseTitle;
	case 1 :
		return BrdfStrings::KSpecularTitle;
	case 2 :
		return BrdfStrings::KShininessTitle;
	case 3 :
		return BrdfStrings::KRefractionTitle;
	case 4 :
		return BrdfStrings::KIORTitle;
	case 5 :
		return BrdfStrings::KExtinctionCoefTitle;
	case 6 :
		return BrdfStrings::KMaxPhongTitle;
	case 7 :
		return BrdfStrings::KBlinnNormalsTitle;
	case 8 :
		return BrdfStrings::KRetroReflectiveTitle;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}

void
CSurfacePhongBase::SetParameter(int parameterIndex, const CValues& value)
{
	CConstantTexture* newTex = 0;
	bool boolValue;
	if ((parameterIndex >= 0) && (parameterIndex < 6))
	{
		newTex = new CConstantTexture(value);
	}
	if ((parameterIndex >= 6) && (parameterIndex < 9))
	{
		boolValue = static_cast<bool>(value[0]);
	}
	switch (parameterIndex)
	{
	case 0 :
		SetDiffuse(newTex);
		break;
	case 1 :
		SetSpecular(newTex);
		break;
	case 2 :
		SetShininess(newTex);
		break;
	case 3 :
		SetRefract(newTex);
		break;
	case 4 :
		SetIOR(newTex);
		break;
	case 5 :
		SetExtinction(newTex);
		break;
	case 6 :
		SetMaxPhong(boolValue);
		break;
	case 7 :
		SetBlinnNormals(boolValue);
		break;
	case 8 :
		SetRetroReflective(boolValue);
		break;
	default :
		break;
	}
}
//========================================================================

/**
   Returns the transparency at the specified point, possibly taking the
   incoming direction of the ray into account
*/
//========================================================================
bool
CSurfacePhongBase::GetTransparency(CValues &transparency,
			       const CHitPointInfo &hitInfo) const
{
  transparency.Set(0.0);

  if (_refract) {
    hitInfo.GetTextureValue(transparency, *_refract, CMappingBase::EE_KR);
    if (transparency.AbsSum() > 0.0)
      return false; // transparent surface - maybe
    return true;    // opaque surface
  }

  transparency.Set(0.0);
  return true; // opaque surface
}
//========================================================================

// Get default color of the surface at the specified hit point
//========================================================================
bool
CSurfacePhongBase::GetDefaultColor(CValues &defaultColor,
			       const CHitPointInfo &hitInfo) const
{
  // let us consider _diffuse color to be the default one
  if (_diffuse) {
    hitInfo.GetTextureValue(defaultColor, *_diffuse, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}
//========================================================================


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

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

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

// Get maximum of albedo
//========================================================================
bool
CSurfacePhongBase::GetAlbedoMax(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

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

  // compute albedo for specular part
  //----------------------------------------------------------------------
  if ((distType & _phongReflectionType()) && _specular) {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
    c += d;
    retval = true;
  }

  if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
    c += d;
    retval = true;
  }

//#define _DEBUG_PHYSICAL_PLAUSIBILITY

#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#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

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

/* Albedo for Phong reflectance model is different for different incoming
   light directions.*/
//========================================================================
bool
CSurfacePhongBase::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);

  CVector3D N = hitInfo.GetNormal();
  float cosine_inc = hitInfo.GetDotProd();

  bool retval = false;

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

  // Specular part of albedo. Albedo is an integral over spherical digon
  // of BRDF function multiplied by cosine of the angle between the
  // outgoing direction and normal. The approximation of integrals of this
  // type was adopted from [2], (31), (original paper [10])
  if ((distType & _phongReflectionType()) && _specular) {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);

    CValues shininess;
    hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
    float nnn = shininess.Get(0);

	if (cosine_inc < 0.0 && d.Max() > WTINY) // avoid heavy numerical computation when specular component is (0,0,0)
    {
      if(fabs(cosine_inc) <= (1.0 - 1e-06))
      {
        CVector3D Rm = hitInfo.GetReflectedDir();

        float integral = IntegralCosNCos(Rm, N, nnn);
        float tmpt = (nnn+2.0)*integral/(2.0*M_PI);
        assert(tmpt >= -1e-06 && tmpt <= 1.0+1e-6);
        d *= tmpt;
      }
      c += d;
      retval = true;
    }
    //else - incoming direction is under the surface
  }

  // Refractive part of albedo
  if ((distType & EE_REFRACTED_PERFECT) && _refract) {
    hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
    // For physically valid BSDF (which represents either reflection from
    // the boundary of an opaque object or the interaction between two
    // non-adsording media), the total_albedo <= 1 ([7], p.184)
    // Total albedo is definded as the sum of all albedos ([8], Ch.2, Ph.D)
    // 24.11.03 K.P.
    CValues ior(CValues::EE_Weigth1);
    // get value of index of refraction  coefficient for intersection point
    if (_index)
      hitInfo.GetTextureValue(ior, *_index, CMappingBase::EE_KRI);
    else
      ior.Set(1.0);

    // For simplicity we assume here that ks(theta_inc) = const and
    // kr(theta_inc) = const. Actually, these parameters depend on refract-
    // tion coefficient and defined by the Fresnel equation
    assert(ior.GetType() == CValues::EE_Weigth1);
    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 || c.Min() < 0.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#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

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


// Get average diffuseness
//========================================================================
float
CSurfacePhongBase::GetAverageDiffuseness(void) const
{
  float diffuse=0,specular=0,n=0,refract=0;
  CValues c;

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

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

    _specnum->GetValue(c, CHitPointInfo::constMapInfo);
    n = c.Luminance();
  }

  if (_refract) {
    // SIG03 - to be removed
    //assert(_refract->GetType()==CMappingBase::EE_None);

    _refract->GetValue(c,CHitPointInfo::constMapInfo);
    refract=c.Luminance();
  }

  // result is the weighted average
  float res =
    ( diffuse + (specular * _nToDiffuseness(n)) ) / (diffuse + specular);
  res *= (refract >= 1.0) ? 1.0 : 1.0 - refract;
  return res;
}
//========================================================================


// Get value of diffuseness for particular incoming direction
//========================================================================
float
CSurfacePhongBase::GetDiffuseness(const CHitPointInfo &hitInfo,
			      EDistributionType distType) const
{
  float diffuse=0,specular=0,n=0,refract=0;
  CValues c;

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

  if ( (distType & EE_REFLECTED_GLOSSY) && _specular) {
    hitInfo.GetTextureValue(c, *_specular, CMappingBase::EE_KS);
    specular = c.Luminance();
    hitInfo.GetTextureValue(c, *_specnum, CMappingBase::EE_KSI);
    n = c.Luminance();
  }

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

  if ( (refract >= 1.0)  ||
       (diffuse + specular < 1e-5) )
    return 0;

  float res=
    ( diffuse + (specular * _nToDiffuseness(n)) ) / (diffuse + specular);
  return res * (1.0 - refract);
}
//========================================================================

// Scale the color 'res' with inverted albedo of refractive part
//========================================================================
inline void
CSurfacePhongBase::ScaleWithTransmission(CValues &res,
				     const CHitPointInfo &hitInfo) const
{
  // If transmission is non-zero, we must multiply the computed
  // value of brdf by (1-kr). Otherwise the transparent object looks bad.
  CValues kr(CValues::EE_Weigth3);
  // get the value of transmission
  hitInfo.GetTextureValue(kr, *_refract, CMappingBase::EE_KR);

  // multiply the resulting BSDF by max(0, 1-kr) in every wavelength
  res.ScaleOneMinusMultiplier(kr);

  return;
}







/**********************************************************************
 *
 *             CSurfacePhong methods
 *
 **********************************************************************/

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



//========================================================================
// Specular part for stochastic sampling
//
// there is no description of importance sampling of original Phong model.
// Physically plausible Phong BRDF sampling is used instead.
//========================================================================
int
CSurfacePhong::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  // TODO: difraction => RGB values for shininess
  //assert(shininess.GetCount() == 1);
  if((kslocal[0] < 1e-5 && kslocal[1] < 1e-5) && kslocal[2] < 1e-5)
    return 0;

  int num = 0;			    // the number of samples generated

  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  // now holds DotProd(N, incomingRay) < 0.0
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  const CVector3D &R = hitInfo.GetReflectedDir();
  // Make new basis with R as Z axis
  CVector3D U = ArbitraryNormal(R); // U is perpendicular to R - X-axis
  CVector3D V = CrossProd(R, U);    // V is perpendicular to R and V - Y-axis

  // diffuseness of the bounce
  float diffuseness = _nToDiffuseness(shininess.Sum() / shininess.GetCount());

  float n = shininess.Get(0);
  // value, required to evaluate probability density for outgoing direction
  float multi = (n + 1.0f)/(M_PI * 2.0f);

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    float lightCos = -1;
    CVector3D dirRes;
    float cos_ang_V_R;
    bool rejected = false;
    // Rejection sampling for Phong model. It can happens, that random
    // direction will be generated under the surface. In this case we should
    // reject the sample and make the new one (in other case sample will be
    // considered as absorbed, but it contradicts with the fact that
    // exactly maxRays must be spawned by this kind of surface interaction).
    while(lightCos < 0) {
       // TODO: Halton sequence
       //if (num) {
         // increase the seed not the generate the same sample again
         //samplingInfo.IncSeed();
         // generate a new sample for this time
         //aux_rnd = samplingInfo.Sample(0);
       //}
       CVector2D h;
       if( !rejected )
         h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
       else
         h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

       // If n->0 we get the distribution close to diffuse distribution and
       // if n->inf we get a ray close to the input (reflected) ray [2], p.19
       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 R is Z-axis) to the World coordinate system
       dirRes = localDir.x * U + localDir.y * V + localDir.z * R;

	   cos_ang_V_R = DotProd(R, dirRes);
       lightCos = DotProd(N, dirRes);

       // Reject only samples which are outside the hemisphere
       if(lightCos < 0.0)
          rejected = true;
      }
      // the ray is sampled in this direction
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

      // compute the probability density
      
      contrib->prob = multi * pow(cos_ang_V_R, n);

      // this value is required for BRDF evaluation
      float multi2 = ((n + 2.0f)/(n + 1.0f));

      // Set the value of BSDF type
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      /* compute = value_of_BRDF / probability_density_of_generated direction
         or more precisely:
         compute = BRDF / pdf =
                 = [ps*(n+2)/(2*Pi)*cos^n]/[(n+1)*cos^n/(2*Pi)] =
                 = ps * (n+2)/(n+1)
         Note, ps is the albedo which is given as input parameter. */
      contrib->weigth = kslocal * multi2;

      // set the reflected ray, direction is normalized here
      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_SPECULAR;
      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;

      // Additional (diffuse) components in order to provide nice and simple
      // vizualisation
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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;
      }
      // Cosine between the incoming ray and the normal does not taken
      // into account.
      contrib->weigth *= lightCos;

#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
  } // for i

  // return number of generated samples
  assert(num == maxRays);
  return num;
}
//========================================================================
// kd_local is albedo => it should be divided by PI to be used
//   as the value of BRDF.
//========================================================================
int
CSurfacePhong::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  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

  int num = 0;			        // number of rays created

  // Compute data for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    // TODO: Halton sequence
    //if (num) {
      // increase the seed not the generate the same sample again
      //samplingInfo.IncSeed();
      // generate a new sample for this time
      //aux_rnd = samplingInfo.Sample(0);
    //}
     // 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, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // This means that we should not multiply the weight by the cosine
    // of the angle between the incoming light direction and normal
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = kd_local;

    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
    CValues ks_local(CValues::EE_Weigth3);
    if (_specular) {
      if(_specular->IsFresnel())
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
      else
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);

      if(ks_local.Luminance() > 1e-5) {
        float alphaCos = DotProd(dirRes,hitInfo.GetReflectedDir());

	// This checking is necessary, otherwise, samples, which are outside
	// the cosine lobe whould be incorrectly involved into the model
	if(alphaCos > 0.0)
	{
          CValues shininess;
		  float incEyeCos   = -(hitInfo.GetDotProd());
          hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
          float n = shininess.Get(0);
          float lightCos = DotProd(N, dirRes);

  	      ks_local *= lightCos*pow(alphaCos,n)/(M_PI*incEyeCos*contrib->prob);
          contrib->weigth += ks_local;
	}
      }
    }

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

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

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



// Gather all parts for stochastic sampling together. Scaling approach
// of random values for stochastic methods was adopted from [9]
//========================================================================
int
CSurfacePhong::SpawnStochasticRays(
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo,
  CSamplingInfo        &samplingInfo,
  EDistributionType     distType) const
{
	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);

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

	// validate hitInfo.dotprod
	hitInfo.GetNormal();

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

	// compute albedo for specular part
	//----------------------------------------------------------------------
	if (_specular && (distType & EE_REFLECTED_SPECULAR)) {
		hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
		aks = ks_local.Luminance();
	}

	// Complete albedo of Phong with mirror reflection at grazing angles model
	//----------------------------------------------------------------------
	float albedo_whole = akd + aks;

	if(albedo_whole < 1e-05) {
	WARNING_INFO("CSurfacePhong::SpawnStochasticRays() :\n\
Albedo of this surface < 1e-05 => it is the black body!");
	return 0;
	}

	if( albedo_whole >= 1.0 ) {
		WARNING_INFO("CSurfacePhong::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);

		// 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 chosen. Therefore, we should 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));
			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++;
		}
		
		// For a now it does not work for Halton sequence.
		samplingInfo.IncDimension(2);
		samplingInfo.RecoverSeed();
	} // while count

	int numRays = 0;


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

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

	// The number of rays for specular interaction
	//----------------------------------------------------------------------
	if (numSpecular) {
	CValues shininess;
	hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
	numRays += SpawnStochasticSpecularRays(ks_local, shininess, carray,
					   hitInfo, samplingInfo,
					   numSpecular);
	}
	int indexMax = carray.GetUsedItems();
	assert(indexMax > 0);

	// Scaling for distributed ray tracing
	/*
	for(int i = 0; i < indexMax; i++) 
	{
	CContributionItem* ci = carray.GetItem(i);
	float mult;

	if(i < indexDiffuse) 
	{
		mult = akd;
		assert(numDiffuse > 0);
		ci->weigth *= (1.0 / ((float)numDiffuse));
	}
	else 
	{
		mult = aks;
		assert(numSpecular > 0);
		ci->weigth *= (1.0 / ((float)numSpecular));
	}
	mult /= albedo_whole;
	} // for*/
	samplingInfo.IncDimension(2);

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

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

// Direct lighting computation
//
//						  kd	  ks * cos(alfa)^n
// original Phong BRDF = ---- + --------------------
//						  Pi	 Pi * cos(thetaInc)
//
// where alfa is angle between the viewer direction and direction of
// perfectly reflected ray, n is shininess of the material and thetaInc
// is theta angle between light direction nad surface normal
//
//========================================================================
bool
CSurfacePhong::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &dirLit,
				    EDistributionType distType) const
{
  bool retval       = false;
  CVector3D  normal  = hitInfo.GetNormal();
  const CVector3D& dirInit = (hitInfo.GetRay())->GetDir();
  float incEyeCos   = -(hitInfo.GetDotProd());
  float lightCos    = DotProd(normal, dirLit);

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

    // Diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      CValues diffuseTerm(bsdf.GetType());
      hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);
      bsdf = diffuseTerm * M_1_PI;
      retval = true;
    }

    // Specular part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CValues specularTerm(bsdf.GetType());

    // if surface has specular component
    if ((distType & EE_REFLECTED_GLOSSY) && _specular) {
      // now get the index of imperfection for specular reflection
      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);

      // the cosine angle between the ideally reflected direction and
      // the viewing (eye) direction, or between the halfvector and normal
      float ratio = DotProd(hitInfo.GetReflectedDir(), dirLit);

      // It is possible that dirLit goes out of cosine lobe, in this case
      // direct illumination should be 0.
      if (ratio > 0.0) {
        // get the value of specular component
        bool use_fresnel = _specular->IsFresnel();
        // Specular texture for phong can be fresnel texture to modelate metalic materials
        hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KS);

        // shininess (coeffitient of imperfection) can have 3 components (RGB)
	    if (specularTerm.GetCount() == specularIndex.GetCount()) {
	    for (int i = specularTerm.GetCount(); i--;)
	       specularTerm.Mul(i, pow(ratio, specularIndex.Get(i)) * 1.f/(M_PI * incEyeCos));
        }
        else {
		  assert(specularIndex.GetCount() == 1);
		  specularTerm *= (pow(ratio, specularIndex.Get(0)) * 1.f/(M_PI * incEyeCos));
        }
        bsdf += specularTerm;
        retval = true;
      }
    } // specular term


  // multiply the result by cosine between the light direction and normal
  if (retval) {
      bsdf *= lightCos;
    }
  }
  else
  {
    retval = false;
  }

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

const char* 
CSurfacePhong::GetHelp() const
{
	return BrdfStrings::KPhongDesc;
}
//========================================================================

int 
CSurfacePhong::GetNumParameters() const
{
	return 3; // diffuse,specular,specnum
}
//========================================================================

CValues
CSurfacePhong::GetParameter(int parameterIndex) const
{
	if (parameterIndex>2)
		return CValues(0,0,0);
	return CSurfacePhongBase::GetParameter(parameterIndex);
}

//========================================================================
CValues
CSurfacePhong::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex>2)
		return CValues(0,0,0);
	return CSurfacePhongBase::GetParameterConstraints(parameterIndex);
}
//========================================================================

const char* 
CSurfacePhong::GetParameterDescription(int parameterIndex) const
{
	if (parameterIndex>2)
		return BrdfStrings::KInvalidParameter;
	return CSurfacePhongBase::GetParameterDescription(parameterIndex);
}
//========================================================================

CSurface::EParameterUIType
CSurfacePhong::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 : //diffuse
	case 1 : //specular
		return EE_ColorRGB;
	case 2 : //specnum
		return EE_Float;
	default :
		return EE_Unknown;
	}
}
//========================================================================

const char*
CSurfacePhong::GetParameterTitle(int parameterIndex) const
{
	if (parameterIndex>2)
		return BrdfStrings::KInvalidParameter;
	return CSurfacePhongBase::GetParameterTitle(parameterIndex);
}
//========================================================================

void
CSurfacePhong::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex>2)
		return;
	CSurfacePhongBase::SetParameter(parameterIndex,value);
}
//========================================================================






/**********************************************************************
 *
 *             CSurfacePhongPhysical methods
 *
 **********************************************************************/
int
CSurfacePhongPhysical::SpawnIdealRefractedRay(const CValues& kr_local,
					      CContributionArray  &carray,
					      const CHitPointInfo &hitInfo)
  const
{
  CContributionItem *contrib = carray.AllocateNext();
  contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

  CVector3D N = hitInfo.GetNormal();
  float dot   = hitInfo.GetDotProd();
  CVector3D dirInit = hitInfo.GetRay()->GetDir();
  CVector3D invInit = -dirInit;

  float gamma, sqrterm;

  CValues ior(CValues::EE_Weigth1);
  if (_index)
    // get value of index of refraction  coefficient for intersection point
    hitInfo.GetTextureValue(ior, *_index, CMappingBase::EE_KRI);
  else {
    ior.Set(1.0);
#if 0
    DEBUG << "CSurfacePhongPhysical::SpawnIdealRefractedRay - ior undefined"
	  << endl;
#endif
  }

  assert(ior.GetType() == CValues::EE_Weigth1);

  // Incident ray - eyeRay (from the eye to the surface) and dot(N,eye) < 0
  // (if normal is oriented in opposite to the eyeRay direction). In this
  // case eyeRay falls from the air onto the object and refracts in the
  // object. And gamma should be equal to n2 (object's index of refraction)
  // Otherwise, if dot(N,eye) < 0 => gamma = 1/n2
  CVector3D dir(0,0,1);
  float sinrefr = (1.0 - dot*dot);

  // Compute direction of outcoming (to the light source) ray
  //---------------------------------------------------------------------------
  if(fabs(dot) > 1e-10 && fabs(dot) < 1.0f - 1e-10)
  {
    if(dot <= 0.0)
      gamma = 1.0/ior.Get(0);
    else
    {
      gamma = ior.Get(0);
      N = -N;
    }

    // compute the cosine between the refracted ray and the normal
    //--------------------------------------------------------------------
    // eta        = n2/n1 - relative index of refraction
    // gamma      = 1/eta
    // sin_refr   = sin_inc/eta = sin_inc * gamma
    // cos_refr^2 = 1 - sin_refr^2 = 1 - sin_inc^2 * gamma^2 =
    //            = 1 - gamma^2 * (1 - dot_inc * dot_inc)
    sqrterm = 1.0 - gamma * gamma * sinrefr;

    // the direction of refracted ray

    // Check for total internal reflection (cos(N,incoming) > 0 in this case:
    // ray falls onto the internal side of the object's surface and object's
    // index of refraction is more than outside media index of refraction.
    // Total interreflection happens when sinus of incoming angle is more
    // than 1/eta - Brewster's angle(). Send ray in mirror direction if
    // it happens.
    if (sqrterm <= 0.0)  {
      // Incoming angle is close to Brewster angle.  Refracted ray goes
      // along the boundary of two media.
      if(sqrterm > -1e-10)
      {
        // Ref[le/ra]cted ray is on the boundary of two media [11], p.176
        dir = (dirInit + fabs(dot)*N)*(1.0/sqrt(sinrefr));
        // set the reflected ray, direction is normalized here
        contrib->interaction = EE_REFRACTED_PERFECT;
        // set the reflected ray, direction is normalized here
        contrib->ray.Init(hitInfo.GetPoint(), dir , REFRACTED_RAY);
        contrib->SetFlags(CContributionItem::EE_RefractedRay);
        float cos_inc_out = DotProd(dir, dirInit);
        assert(cos_inc_out > 0);
      }
      else
      { // Total internal reflection (incoming angle > Brewster's angle)
        dir = N*2.0*fabs(dot) + dirInit;
        contrib->interaction = EE_REFLECTED_PERFECT;
        // set the reflected ray, direction is normalized here
        contrib->ray.Init(hitInfo.GetPoint(), dir , REFLECTED_RAY);
        contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      }
      assert(dot >= 0.0);
    }
    else {
      // Ordinary refraction. There are the following cases:
      // * ray falls from outside the object onto the object's surface,
      //   object coefficient of refraction is more than outside surface's
      //   coeffictient of refraction.
      // * ray falls from the inside the object onto the internal side if the
      //   object's surface. Sinus of incoming angle is less than Brewster
      //   angle.
      sqrterm = -fabs(dot) * gamma + sqrt(sqrterm);

      // Transmitted ray direction: [1], p.578
      dir = -sqrterm * N + gamma*dirInit;
      contrib->interaction = EE_REFRACTED_PERFECT;
      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dir, REFRACTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_RefractedRay);

      float cos_out_inc = DotProd(dir, dirInit) - WTINY;
      assert(cos_out_inc > 0);
      float out_inc_a   = acos(cos_out_inc);
      float lightCos     = fabs(DotProd(N, dir)) - WTINY;
      float out_a       = acos(lightCos);
      float inc_a       = acos(fabs(dot));
      assert(out_inc_a < (out_a + inc_a));
    }
  }
  else // Boundary cases: the angle between the incoming direction and
  {    // the normal is 0 or 90 (normal incidence or grazing angle)
    dir = dirInit;

    if(fabs(dot) <= 1e-6) // Incoming angle is 90 -> ray is grazing the object
    {
      contrib->interaction = EE_REFLECTED_PERFECT;
      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dir, REFLECTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
    }
    else
    { // Incomig angle is 0 -> transmission without refraction
      contrib->interaction = EE_REFRACTED_PERFECT;
      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dir, REFRACTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_RefractedRay);
    }
  }
  float len = Magnitude(dir);
  assert(len <= 1.0 + 1e-06 && len >= 1.0 - 1e-06);

  // diffuseness of the bounce
  contrib->diffuseness = 0.0;
  // the reflection probability is 1 for this type of reflection
  contrib->prob = 1.0;

  // In [12] there is a note, that
  //       BRDF = ps*Dirac()/cos(normal, light) =>
  // In this case to avoid multiplication and division by the same term
  // (weight = BRDF * cos(normal, light)) we should not divide here by
  // cosine and should not multiply by cosine in the PathTracing
  // correspondingly.  05/12/03 KP.
  contrib->weigth = kr_local;

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

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


// Specular part for stochastic sampling
//
// Coefitient ks_local is not multiplied by (n+2.0)/(2.0*PI) => it should be
// taken into account here. The value of specular part of BSDF is
//     ks_local * (cos(phi))^n * (n+2.0)/(2.0*PI),
// where phi is the angle between ideally reflected ray and outgoing ray.
//========================================================================
int
CSurfacePhongPhysical::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  // TODO: difraction => RGB values for shininess
  //assert(shininess.GetCount() == 1);
  if((kslocal[0] < 1e-5 && kslocal[1] < 1e-5) && kslocal[2] < 1e-5)
    return 0;

  int num = 0;			    // the number of samples generated

  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  // now holds DotProd(N, incomingRay) < 0.0
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  const CVector3D &R = hitInfo.GetReflectedDir();
  // Make new basis with R as Z axis
  CVector3D U = ArbitraryNormal(R); // U is perpendicular to R - X-axis
  CVector3D V = CrossProd(R, U);    // V is perpendicular to R and V - Y-axis

  // diffuseness of the bounce
  float diffuseness = _nToDiffuseness(shininess.Sum() / shininess.GetCount());

  float n = shininess.Get(0);
  // value, required to evaluate probability density for outgoing direction
  float multi = (n + 1.0f)/(M_PI * 2.0f);

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    float lightCos = -1;
    CVector3D dirRes;
    float cos_ang_V_R;
    bool rejected = false;
    // Rejection sampling for Phong model. It can happens, that random
    // direction will be generated under the surface. In this case we should
    // reject the sample and make the new one (in other case sample will be
    // considered as absorbed, but it contradicts with the fact that
    // exactly maxRays must be spawned by this kind of surface interaction).
    while(lightCos < 0) {
       // TODO: Halton sequence
       //if (num) {
         // increase the seed not the generate the same sample again
         //samplingInfo.IncSeed();
         // generate a new sample for this time
         //aux_rnd = samplingInfo.Sample(0);
       //}
       CVector2D h;
       if( !rejected )
         h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
       else
         h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

       // If n->0 we get the distribution close to diffuse distribution and
       // if n->inf we get a ray close to the input (reflected) ray [2], p.19
       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 R is Z-axis) to the World coordinate system
       dirRes = localDir.x * U + localDir.y * V + localDir.z * R;

       // If cosine lobe is computed according to the Blinn article [4], in
       // SpawnStochasticSpecularRay function this fact should be also
       // taken into account
       if (!_normalsBlinn) {
         // evaluation of the specular component according to Phong
	 cos_ang_V_R = DotProd(R, dirRes);
       }
       else {
         CVector3D hVector = 0.5 * (-dirInit + dirRes);
         hVector.Normalize();
         cos_ang_V_R = DotProd(hVector, N);
       }
       lightCos = DotProd(N, dirRes);

       // Reject only samples which are outside the hemisphere
       if(lightCos < 0.0)
          rejected = true;
      }
      // the ray is sampled in this direction
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

      // compute the probability density
      if ( _retroReflective)
        contrib->prob = multi * pow(fabs(DotProd(-dirInit,dirRes)),n);
      else
        contrib->prob = multi * pow(cos_ang_V_R, n);

      // this value is required for BRDF evaluation
      float multi2 = ((n + 2.0f)/(n + 1.0f));

      // Set the value of BSDF type
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      /* compute = value_of_BRDF / probability_density_of_generated direction
         or more precisely:
         compute = BRDF / pdf =
                 = [ps*(n+2)/(2*Pi)*cos^n]/[(n+1)*cos^n/(2*Pi)] =
                 = ps * (n+2)/(n+1)
         Note, ps is the albedo which is given as input parameter. */
      contrib->weigth = kslocal * multi2;

      if (_maxPhong || _retroReflective) {
         // if max-Phong model is used divide by cosine
         float incEyeCos  = fabs(hitInfo.GetDotProd());
         float lightCos = fabs(DotProd(N, dirRes));
         float max_cos = Max(incEyeCos, lightCos);
         contrib->weigth *= (1.0f/max_cos);
      }
      // set the reflected ray, direction is normalized here
      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_SPECULAR;
      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;

      // Additional (diffuse) components in order to provide nice and simple
      // vizualisation
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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;
      }
      // Cosine between the incoming ray and the normal does not taken
      // into account.
      contrib->weigth *= lightCos;

#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
  } // for i

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


// kd_local is albedo => it should be divided by PI to be used
//   as the value of BRDF.
//========================================================================
int
CSurfacePhongPhysical::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  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

  int num = 0;			        // number of rays created

  // Compute data for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    // TODO: Halton sequence
    //if (num) {
      // increase the seed not the generate the same sample again
      //samplingInfo.IncSeed();
      // generate a new sample for this time
      //aux_rnd = samplingInfo.Sample(0);
    //}
     // 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, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // This means that we should not multiply the weight by the cosine
    // of the angle between the incoming light direction and normal
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = kd_local;

    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
    CValues ks_local(CValues::EE_Weigth3);
    if (_specular) {
      if(_specular->IsFresnel())
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
      else
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);

      if(ks_local.Luminance() > 1e-5) {
        float alphaCos = DotProd(dirRes,hitInfo.GetReflectedDir());

	// This checking is necessary, otherwise, samples, which are outside
	// the cosine lobe whould be incorrectly involved into the model
	if(alphaCos > 0.0)
	{
          CValues shininess;
          hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
          float n = shininess.Get(0);
          float lightCos = DotProd(N, dirRes);

  	  ks_local *= (n+2.0)*lightCos*pow(alphaCos,n)/(2.0*M_PI*contrib->prob);
          contrib->weigth += ks_local;
	}
      }
    }

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

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

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


// Specular part for stochastic sampling - multithreaded optimized variant
// (it uses local random arrays)
//========================================================================
int
CSurfacePhongPhysical::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  const rnd_array     &rnd_nums,
  int                 maxRays) const
{
  if((kslocal[0] < 1e-5 && kslocal[1] < 1e-5) && kslocal[2] < 1e-5)
    return 0;

  int num = 0;			    // the number of samples generated

  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  // now holds DotProd(N, incomingRay) < 0.0
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  const CVector3D &R = hitInfo.GetReflectedDir();
  // Make new basis with R as Z axis
  CVector3D U = ArbitraryNormal(R); // U is perpendicular to R - X-axis
  CVector3D V = CrossProd(R, U);    // V is perpendicular to R and V - Y-axis

  // diffuseness of the bounce
  float diffuseness = _nToDiffuseness(shininess.Sum() / shininess.GetCount());

  float n = shininess.Get(0);
  // value, required to evaluate probability density for outgoing direction
  float multi = (n + 1.0f)/(M_PI * 2.0f);

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    float lightCos = -1;
    CVector3D dirRes;
    float cos_ang_V_R;
    bool rejected = false;
    // Rejection sampling for Phong model. It can happens, that random
    // direction will be generated under the surface. In this case we should
    // reject the sample and make the new one (in other case sample will be
    // considered as absorbed, but it contradicts with the fact that
    // exactly maxRays must be spawned by this kind of surface interaction).
    while(lightCos < 0) {
       // TODO: Halton sequence
       //if (num) {
         // increase the seed not the generate the same sample again
         //samplingInfo.IncSeed();
         // generate a new sample for this time
         //aux_rnd = samplingInfo.Sample(0);
       //}
       CVector2D h;
       if( !rejected )
         h.SetValue(rnd_nums[2*num], rnd_nums[2*num+1]);
       else
         h.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

       // If n->0 we get the distribution close to diffuse distribution and
       // if n->inf we get a ray close to the input (reflected) ray [2], p.19
       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 R is Z-axis) to the World coordinate system
       dirRes = localDir.x * U + localDir.y * V + localDir.z * R;

       // If cosine lobe is computed according to the Blinn article [4], in
       // SpawnStochasticSpecularRay function this fact should be also
       // taken into account
       if (!_normalsBlinn) {
         // evaluation of the specular component according to Phong
	 cos_ang_V_R = DotProd(R, dirRes);
       }
       else {
         CVector3D hVector = 0.5 * (-dirInit + dirRes);
         hVector.Normalize();
         cos_ang_V_R = DotProd(hVector, N);
       }
       lightCos = DotProd(N, dirRes);

       // Reject only samples which are outside the hemisphere
       if(lightCos < 0.0)
          rejected = true;
      }
      // the ray is sampled in this direction
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

      // compute the probability density
      if ( _retroReflective)
        contrib->prob = multi * pow(fabs(DotProd(-dirInit,dirRes)),n);
      else
        contrib->prob = multi * pow(cos_ang_V_R, n);

      // this value is required for BRDF evaluation
      float multi2 = ((n + 2.0f)/(n + 1.0f));

      // Set the value of BSDF type
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      /* compute = value_of_BRDF / probability_density_of_generated direction
         or more precisely:
         compute = BRDF / pdf =
                 = [ps*(n+2)/(2*Pi)*cos^n]/[(n+1)*cos^n/(2*Pi)] =
                 = ps * (n+2)/(n+1)
         Note, ps is the albedo which is given as input parameter. */
      contrib->weigth = kslocal * multi2;

      if (_maxPhong || _retroReflective) {
         // if max-Phong model is used divide by cosine
         float incEyeCos  = fabs(hitInfo.GetDotProd());
         float lightCos = fabs(DotProd(N, dirRes));
         float max_cos = Max(incEyeCos, lightCos);
         contrib->weigth *= (1.0f/max_cos);
      }
      // set the reflected ray, direction is normalized here
      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_SPECULAR;
      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;

      // Additional (diffuse) components in order to provide nice and simple
      // vizualisation
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      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;
      }
      // Cosine between the incoming ray and the normal does not taken
      // into account.
      contrib->weigth *= lightCos;

#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
  } // for i

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

// Diffuse rays spawning - multithreading optimized variant
//========================================================================
int
CSurfacePhongPhysical::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  const rnd_array     &rnd_nums,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  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

  int num = 0;			        // number of rays created

  // Compute data for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    // TODO: Halton sequence
    //if (num) {
      // increase the seed not the generate the same sample again
      //samplingInfo.IncSeed();
      // generate a new sample for this time
      //aux_rnd = samplingInfo.Sample(0);
    //}
     // 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(rnd_nums[2*num], rnd_nums[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, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // This means that we should not multiply the weight by the cosine
    // of the angle between the incoming light direction and normal
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = kd_local;

    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
    CValues ks_local(CValues::EE_Weigth3);
    if (_specular) {
      if(_specular->IsFresnel())
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
      else
        hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);

      if(ks_local.Luminance() > 1e-5) {
        float alphaCos = DotProd(dirRes,hitInfo.GetReflectedDir());

	// This checking is necessary, otherwise, samples, which are outside
	// the cosine lobe whould be incorrectly involved into the model
	if(alphaCos > 0.0)
	{
          CValues shininess;
          hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
          float n = shininess.Get(0);
          float lightCos = DotProd(N, dirRes);

  	  ks_local *= (n+2.0)*lightCos*pow(alphaCos,n)/(2.0*M_PI*contrib->prob);
          contrib->weigth += ks_local;
	}
      }
    }

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

  } // for i
  return maxRays;
}


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

  // this is multithreading optimized spawning - more than one
  // thread can access random numbers array = this method,
  // because each run of this method and all the called methods 
  // has its own local arrays with random numbers
  rnd_array local_rnd_diffuse;
  rnd_array local_rnd_specular;
  rnd_array local_rnd_refraction;

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

  //hitInfo.GetNormal();
  EDistributionType tmp_distType;

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

  // compute albedo for specular part
  //----------------------------------------------------------------------
  if (_specular && (distType & EE_REFLECTED_SPECULAR)) {
    if(_specular->IsFresnel())
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
    else
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
    tmp_distType = EE_REFLECTED_SPECULAR;
    aks = ks_local.Luminance();
  }

  // compute albedo for refractive part
  //----------------------------------------------------------------------
  if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
    tmp_distType = EE_REFRACTED_PERFECT;
    akr = kr_local.Luminance();
  }
  // Complete albedo of physically correct Phong model
  //----------------------------------------------------------------------
  float albedo_whole = akr + akd + aks;

  if(albedo_whole < 1e-05) {
	WARNING_INFO("CSurfacePhongPhysical::SpawnStochasticRays() :\n\
Albedo of this surface < 1e-05. It is the black body!");
    return 0;
  }
  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CSurfacePhongPhysical::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)");
  }

  // 'Distribute' (precalculate random numbers for each type of interaction:
  // diffuse, specular, refraction) maxRays rays of different type according
  // to the corresponding albedos
  //----------------------------------------------------------------------
  int numRefract = 0, numDiffuse = 0, numSpecular = 0;

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

  int count = maxRays;
  while(count--)
  {  //    We assume, that the Russian roulette is used to choose between
     // the continuation or terninamtion of the random walk. We need to get
     // the 3 random numbers: first one to choose the type of interaction,
     // two the rest to generate the random direction for the choosen type
     // of interaction.
     //    1) Choose the type of interaction
     //    In the first iteration random number should be taken from the user's
     // program (pushedValue). Note, that this random number is LESS than the
     // total albedo because it was already decided that the walk will be
     // continued (if random number is more than the total albedo the
     // absorbtion should be occured in the user's program and this function
     // SpawnStochasticRays should NOT be called). This random number should
     // be scaled by the total albedo to fall into the range [0..1).
     //    For the next iterations it is impossible to reuse that number from
     // the users's program and we should generate this number locally. If we
     // use standard random generator, random number falls into the range
     // [0..1) and should not be scaled.
     //    Note that the albedos to choose the type of interaction should
     // be scaled in this case to sum to one (see [9], p.35 for details)
     //    2) Choose random direction
     //    It is possible to reuse the random number for choosing the
     // type of interaction for generation of the first component of the
     // random direction. It should be rescaled according to [9], p.35. Second
     // component of the random direction should be generated from the scratch.
     float rnd_num = samplingInfo.Sample(0);

     // 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)
     }

    /* This is the rule for choosing the interaction:
     *
     *  +----------------------------------------------------+
     *  |   refract   |  reflect diffusely |    specularly   |
     *  +----------------------------------------------------+
     *  |<-- aplha1 ->|<----- alpha2 ----->|<---- alpha3 --->|
     */
     float alpha1 = akr/albedo_whole;
     float alpha2 = akd/albedo_whole;
     float alpha3 = 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)
       local_rnd_refraction.push_back(rnd_num/alpha1);
       local_rnd_refraction.push_back(samplingInfo.Sample(1));
       numRefract++;
     }
     else {
       if(rnd_num < (alpha1+alpha2) && alpha2 > 1e-6) {
         // Reuse random number for generation of one component of the random
         // direction. Scale it to fall into the range [0..1)
         local_rnd_diffuse.push_back((rnd_num-alpha1)/alpha2);
         // Choose another component of direction randomly
	     local_rnd_diffuse.push_back(samplingInfo.Sample(1));
         numDiffuse++;
       }
       else {
	 // Reuse random number for generation of one component of the random
         // direction. Scale it to fall into the range [0..1)
         local_rnd_specular.push_back((rnd_num-alpha1-alpha2)/alpha3);
         // Choose another component of direction randomly
	     local_rnd_specular.push_back(samplingInfo.Sample(1));
         numSpecular++;
       }
     }
     // For a now it does not work for Halton sequence.
     samplingInfo.IncDimension(2);
     samplingInfo.RecoverSeed();
  }

  int numRays = 0;

  // The number of rays for refraction
  //----------------------------------------------------------------------
  if (numRefract)
    numRays += SpawnIdealRefractedRay(kr_local, carray, hitInfo);

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

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

  int indexDiffuse = carray.GetUsedItems() - indexRefracted;

  // The number of rays for specular interaction
  //----------------------------------------------------------------------
  if (numSpecular) {
    CValues shininess;
    hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
    numRays += SpawnStochasticSpecularRays(ks_local, shininess, carray,
					hitInfo, samplingInfo, local_rnd_specular,
					   numSpecular);
  }
  int indexMax = carray.GetUsedItems();

  // Scaling for distributed ray tracing
  /*
  for(int i = 0; i < indexMax; i++) {
    CContributionItem* ci = carray.GetItem(i);
    float mult;
    if(i < indexRefracted) {
      mult = akr;
    }
    else {
      if(i < (indexRefracted + indexDiffuse)) {
        mult = akd;
        assert(numDiffuse > 0);
        ci->weigth *= (1.0 / ((float)numDiffuse));
      }
      else {
        mult = aks;
        assert(numSpecular > 0);
	ci->weigth *= (1.0 / ((float)numSpecular));
      }
    }
    mult /= albedo_whole;
  } // for*/
  samplingInfo.IncDimension(2);

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

  // WARNING:
  // 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/ KP.

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

// Direct lighting computation
//========================================================================
bool
CSurfacePhongPhysical::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &dirLit,
				    EDistributionType distType) const
{
  bool retval       = false;
  CVector3D  normal  = hitInfo.GetNormal();
  const CVector3D& dirInit = (hitInfo.GetRay())->GetDir();
  float incEyeCos   = -(hitInfo.GetDotProd());
  float lightCos    = DotProd(normal, dirLit);

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

    // Diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      CValues diffuseTerm(bsdf.GetType());
      hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);
      bsdf = diffuseTerm * M_1_PI;
      retval = true;
    }

    // Refraction should not be taken into account for direct illumination
    // computation

    // Specular part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CValues specularTerm(bsdf.GetType());

    // if surface has specular component
    if ((distType & EE_REFLECTED_GLOSSY) && _specular) {
      // now get the index of imperfection for specular reflection
      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);

      // the cosine angle between the ideally reflected direction and
      // the viewing (eye) direction, or between the halfvector and normal
      float ratio;

      if (!_normalsBlinn) {
        // evaluation of the specular component according to Phong
        if ( _retroReflective)
          ratio = fabs(DotProd(dirLit, -dirInit));
        else
          ratio = DotProd(hitInfo.GetReflectedDir(), dirLit);
      }
      else {
        // evaluation of the specular component according to Blinn [4]
        // Halfvector - is the bisector vector between the incoming and
        // outgoing directions K.P. 07/12/03
        CVector3D hVector = -dirInit + dirLit;
        hVector.Normalize();
        if ( _retroReflective)
          ratio = fabs(DotProd(dirLit, -dirInit));
        else
          ratio = DotProd(hVector, normal);
      }

      // It is possible that dirLit goes out of cosine lobe, in this case
      // direct illumination should be 0.
      if (ratio > 0.0) {
        // get the value of specular component
        bool use_fresnel = _specular->IsFresnel();
        // Specular texture for phong can be fresnel texture to modelate metalic materials
        if(use_fresnel)
          hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KFM);
        else
          hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KS);

        // shininess (coeffitient of imperfection) can have 3 components (RGB)
	if (specularTerm.GetCount() == specularIndex.GetCount()) {
	  for (int i = specularTerm.GetCount(); i--;)
	    specularTerm.Mul(i, pow(ratio, specularIndex.Get(i))
	                        * (specularIndex.Get(i)+2.0f)/(M_PI*2.0f));
        }
        else {
	  assert(specularIndex.GetCount() == 1);
	  specularTerm *= (pow(ratio, specularIndex.Get(0))*
	                       (specularIndex.Get(0)+2.0f)/(M_PI*2.0f));
        }
        // compute BSDF by so called max-Phong formula proposed in [3]
	if (_maxPhong || _retroReflective) {
          float max_cos = Max(incEyeCos, lightCos);
          specularTerm *= (1.0f/max_cos);
        }
        bsdf += specularTerm;
        retval = true;
      }
    } // specular term


  // multiply the result by cosine between the light direction and normal
  if (retval) {
      bsdf *= lightCos;
    }
  }
  else
  {
    retval = false;
  }

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

const char* 
CSurfacePhongPhysical::GetHelp() const
{
	return BrdfStrings::KPhongPhysicalDesc;
}
//========================================================================

int 
CSurfacePhongPhysical::GetNumParameters() const
{
	return 8; // diffuse,specular,specnum,refraction,ior,maxPhong,normalsBlinn,retroReflective
}
//========================================================================

CValues
CSurfacePhongPhysical::GetParameter(int parameterIndex) const
{
	if (parameterIndex>4) // we skip 5th parameter of phong model == extinction coefficient
		parameterIndex++;
	return CSurfacePhongBase::GetParameter(parameterIndex);
}

//========================================================================
CValues
CSurfacePhongPhysical::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex>4) // we skip 5th parameter of phong model == extinction coefficient
		parameterIndex++;
	return CSurfacePhongBase::GetParameterConstraints(parameterIndex);
}
//========================================================================

const char* 
CSurfacePhongPhysical::GetParameterDescription(int parameterIndex) const
{
	if (parameterIndex>4) // we skip 5th parameter of phong model == extinction coefficient
		parameterIndex++;
	return CSurfacePhongBase::GetParameterDescription(parameterIndex);
}
//========================================================================

CSurface::EParameterUIType
CSurfacePhongPhysical::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 : //diffuse
		return EE_ColorRGB;
	case 1 : //specular
		return EE_ColorRGB;
	case 2 : //specnum
		return EE_Float; // isotropic
	case 3 : //refract
		return EE_ColorRGB;
	case 4 : //index
		return EE_Float; // isotropic
	case 5 : //maxPhong
		return EE_Boolean;
	case 6 : //normalsBlinn
		return EE_Boolean;
	case 7 : //retroReflective
		return EE_Boolean;
	default :
		return EE_Unknown;
	}
}
//========================================================================

const char*
CSurfacePhongPhysical::GetParameterTitle(int parameterIndex) const
{
	if (parameterIndex>4) // we skip 5th parameter of phong model == extinction coefficient
		parameterIndex++;
	return CSurfacePhongBase::GetParameterTitle(parameterIndex);
}
//========================================================================

void
CSurfacePhongPhysical::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex>4) // we skip 5th parameter of phong model == extinction coefficient
		parameterIndex++;
	CSurfacePhongBase::SetParameter(parameterIndex,value);
}
//========================================================================

/**********************************************************************
 *
 *             CSurfacePhongGrazingMirror methods
 *
 *
 **********************************************************************/
/* This model [6] is physically plausible:
     1. Albedo of Phong model is changed, but is still < 1.
     2. Reciprocity of Phong model is not violated (for ordinary function
        division by cosine of incoming angle makes function assymetric,
	but for Mirror, incoming angle is always the same as outgoing
	mirror-reflected angle).

   Actually, albedo for specular part of Phong model differs with incoming
   angle (it decreases for greater incoming angles). For
   SurfacePhongGrazingMirror model rejected samples are sent into the mirror
   direction and their weight is equal to albedo of specular part. In this
   way, albedo is approximately the same and equal to the sum of albedo of
   refractive, diffuse, and specular parts.  */
//========================================================================
bool
CSurfacePhongGrazingMirror::GetAlbedo(CValues &albedo,
                                      const CHitPointInfo &hitInfo,
                                      EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

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

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

  // Refractive part of albedo
  if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
    c += d;
    retval = true;
  }

//#define _DEBUG_PHYSICAL_PLAUSIBILITY
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#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

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

// Specular part for stochastic sampling
//========================================================================
int
CSurfacePhongGrazingMirror::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  assert(shininess.GetCount() == 1);
  if((kslocal[0] < 1e-5 && kslocal[1] < 1e-5) && kslocal[2] < 1e-5)
    return 0;

  int num = 0;                      // the number of samples to be generated

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

  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  const CVector3D &R = hitInfo.GetReflectedDir();
  CVector3D Rm = R;
  // Make new basis with R as Z axis
  CVector3D U = ArbitraryNormal(R); // U is perpendicular to R - X-axis
  CVector3D V = CrossProd(R, U);    // V is perpendicular to R and V - Y-axis

  // diffuseness of the bounce
  float diffuseness = _nToDiffuseness(shininess.Sum() / shininess.GetCount());

  float n = shininess.Get(0);
  // this value is required to evaluate probability density
  // for outgoing direction
  float multi = (n + 1.0f)/(M_PI * 2.0f);

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    float lightCos = -1;
    CVector3D dirRes;
    float cos_ang_V_R;
    bool rejected = false;

    // TODO: Halton sequence
    //if (num) {
    // increase the seed not the generate the same sample again
    //samplingInfo.IncSeed();
    // generate a new sample for this time
    //aux_rnd = samplingInfo.Sample(0);
    //}

    CVector2D h;
    h.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);

    // See [2], p.19
    // If n->0 we get the distribution close to diffuse distribution and
    // if n->infinity we get a ray close to the input (reflected) ray,
    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 R is Z-axis) to the World coordinate system
    dirRes = localDir.x * U + localDir.y * V + localDir.z * R;

    // Compute cosine of corresponding angle (between the viewing and
    // ideally reflected directions or between the half-vector and the normal)
    // according to Phong or Blinn-Phong [4] formulae
    if (!_normalsBlinn) {
      cos_ang_V_R = DotProd(R, dirRes);
    }
    else {
      CVector3D hVector = -dirInit + dirRes;
      hVector.Normalize();
      cos_ang_V_R = DotProd(hVector, N);
    }
    // Reject only samples which are outside the upper hemisphere
    lightCos = DotProd(N, dirRes);

    if(lightCos < 0.0)
    {
      rejected = true;
    }
    // the BRDF is sampled in this direction
    CContributionItem *contrib = carray.AllocateNext();
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

    // If sample was not rejected, send the rays according to Phong model
    if(!rejected)
    {
      if ( _retroReflective)
        contrib->prob = multi * pow(fabs(DotProd(-dirInit, dirRes)),n);
      else
        contrib->prob = multi * pow(cos_ang_V_R, n);

      // this value is required for BRDF evaluation
      float multi2 = ((n + 2.0f)/(n + 1.0f));

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

      if (_maxPhong || _retroReflective) {
         // if max-Phong model [3]
         float incEyeCos  = fabs(hitInfo.GetDotProd());
         float lightCos = fabs(DotProd(N, dirRes));
         float max_cos = Max(incEyeCos, lightCos);
         contrib->weigth *= (1.0f/max_cos);
      }
      // set the reflected ray, direction is normalized here
      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_SPECULAR;
      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;

      // weight = BSDF * lightCos = ks * pow(cos(Mirror,Eye),n) * lightCos
      contrib->weigth *= lightCos;
    }
    else
    {
      // send the rays in mirror-reflection direction
      contrib->prob = 1;
      // weight = BSDF * lightCos = ks/lightCos * lightCos = ks,  [12]
      contrib->weigth = kslocal;
      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), Rm, REFLECTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      // type of interaction
      contrib->interaction = EE_REFLECTED_MIRROR;
      // diffuseness of the bounce
      contrib->diffuseness = 0;

      lightCos = DotProd(Rm, N);
    }
    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
     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 * lightCos;
       }
     }

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

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

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

  // 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), 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, akr = 0.0f;

  // validate hitInfo.dotprod
  hitInfo.GetNormal();

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

  // compute albedo for specular part
  //----------------------------------------------------------------------
  if (_specular && (distType & EE_REFLECTED_SPECULAR)) {

    if(_specular->IsFresnel())
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KFM);
    else
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
    aks = ks_local.Luminance();
  }

  // compute albedo for refractive part
  //----------------------------------------------------------------------
  if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
    akr = kr_local.Luminance();
  }

  // Complete albedo of Phong with mirror reflection at grazing angles model
  //----------------------------------------------------------------------
  float albedo_whole = akr + akd + aks;

  if(albedo_whole < 1e-05) {
	WARNING_INFO("CSurfacePhongGrazingMirror::SpawnStochasticRays() :\n\
Albedo of this surface < 1e-05 => it is the black body!");
    return 0;
  }
  
  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CSurfacePhongGrazingMirror::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)");
  }
  // 'Distribute' (precalculate random numbers for each type of interaction:
  // diffuse, specular, refraction) maxRays rays of different type according
  // to the corresponding albedos
  //----------------------------------------------------------------------
  int numRefract = 0, numDiffuse = 0, numSpecular = 0;
  int count = maxRays;
  while(count--)
  {
     float rnd_num = samplingInfo.Sample(0);
     if(count >= maxRays-1)
     {
       if(rnd_num <= albedo_whole)
         rnd_num /= albedo_whole;  // make first RND fall into the [0..1)
     }

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

     if(rnd_num < alpha1 && alpha1 > 1e-6) {
       refraction_rnd.push_back(rnd_num/alpha1);
       refraction_rnd.push_back(samplingInfo.Sample(1));
       numRefract++;
     }
     else {
       if(rnd_num < (alpha1+alpha2) && alpha2 > 1e-6) {
         diffuse_rnd.push_back((rnd_num-alpha1)/alpha2);
	 diffuse_rnd.push_back(samplingInfo.Sample(1));
         numDiffuse++;
       }
       else {
         specular_rnd.push_back((rnd_num-alpha1-alpha2)/alpha3);
         specular_rnd.push_back(samplingInfo.Sample(1));
         numSpecular++;
       }
     }
     // For a now it does not work for Halton sequence.
     samplingInfo.IncDimension(2);
     samplingInfo.RecoverSeed();
  }
  int numRays = 0;

  // The number of rays for refraction
  //----------------------------------------------------------------------
  if (numRefract)
    numRays += SpawnIdealRefractedRay(kr_local, carray, hitInfo);

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

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

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

  // The number of rays for specular interaction
  //----------------------------------------------------------------------
  if (numSpecular) {
    CValues shininess;
    hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
    numRays += SpawnStochasticSpecularRays(ks_local, shininess, carray,
					   hitInfo, samplingInfo,
					   numSpecular);
  }
  int indexMax = carray.GetUsedItems();
  assert(indexMax > 0);

  // Scaling for distributed ray tracing
  /*
  for(int i = 0; i < indexMax; i++) {
    CContributionItem* ci = carray.GetItem(i);
    float mult;
    if(i < indexRefracted) {
      mult = akr;
    }
    else {
      if(i < (indexRefracted + indexDiffuse)) {
        mult = akd;
        assert(numDiffuse > 0);
        ci->weigth *= (1.0 / ((float)numDiffuse));
      }
      else {
        mult = aks;
        assert(numSpecular > 0);
	ci->weigth *= (1.0 / ((float)numSpecular));
      }
    }
    mult /= albedo_whole;
  } // for*/
  samplingInfo.IncDimension(2);

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

   // WARNING:
  // 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/ KP.
  assert(indexMax == maxRays);
  return numRays;
}
//========================================================================

const char* 
CSurfacePhongGrazingMirror::GetHelp() const
{
	return BrdfStrings::KPhongGrazingMirrorDesc;
}
//========================================================================

/**********************************************************************
 *
 *             CSurfacePhongPhysicalDR methods
 *
 *    (deterministic reflection)
 *
 **********************************************************************/
/*
  See CSurfacePhongPhysical::SpawnStochasticRays for comments
*/
//========================================================================
/* If specular reflected rays are spawn deterministically => EvaluateBSDF
   function of CSurfacePhongPhysicalDR class should differ from the
   correspondent function of CSurfacePhongPhysical class: specular
   term should not be taken into account at all because for Ideal
   specular interaction emissive power of the Light is taken into
   account when calculating indirect illumination. K.P. 07/12/03
*/
//========================================================================
bool
CSurfacePhongPhysicalDR::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);

  CVector3D N = hitInfo.GetNormal();

  bool retval = false;

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

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

  // Refractive part of albedo
  if ((distType & EE_REFRACTED_PERFECT) && _refract) {
    hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
    c += d;
    retval = true;
  }

//#define _DEBUG_PHYSICAL_PLAUSIBILITY

#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    if (_refract) {
      hitInfo.GetTextureValue(d, *_refract, CMappingBase::EE_KR);
      DEBUG << "refr = " << d << endl;
    }
  }
#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

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

// Direct lighting computation
//=============================================================================
bool
CSurfacePhongPhysicalDR::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &dirLit,
				    EDistributionType distType) const
{
  bool retval = false;
  float lightCos    = DotProd(dirLit, hitInfo.GetNormal());
  float incEyeCos   = -(hitInfo.GetDotProd());

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

    // evaluating diffuse term
    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      CValues diffuseTerm(bsdf.GetType());
      hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);
      bsdf = diffuseTerm * M_1_PI;
      retval = true;
    }

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


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

  // TODO
  // 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();

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

  CValues kd_local(CValues::EE_Weigth3);	// diffuse component
  CValues ks_local(CValues::EE_Weigth3);	// specular component
  CValues kr_local(CValues::EE_Weigth3);	// refractive component

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

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

  // compute albedo for specular part
  //----------------------------------------------------------------------
  if (_specular && (distType & EE_REFLECTED_PERFECT)) {
    hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KM);
    aks = ks_local.Luminance();
  }

  // compute albedo for refractive part
  //----------------------------------------------------------------------
  if (_refract && (distType & EE_REFRACTED_PERFECT)) {
    hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
    akr = kr_local.Luminance();
  }

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

  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CSurfacePhongPhysicalDR::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)");
  }

  if ( albedo_whole < 1e-5) {
    WARNING_INFO("CSurfacePhongPhysicalDR::SpawnStochasticRays()\n\
Albedo of this surface < 1e-05 => it is the black body!");
    samplingInfo.IncDimension(2);
    samplingInfo.RecoverSeed();
    return 0; // no rays -- this surface is black body
  }
  // 'Distribute' (precalculate random numbers for each type of interaction:
  // diffuse, specular, refraction) maxRays rays of different type according to the
  // corresponding albedos
  //----------------------------------------------------------------------
  int numRefract = 0, numDiffuse = 0, numSpecular = 0;
  int count = maxRays;
  while(count--)
  {
     float rnd_num = samplingInfo.Sample(0);
     if(count >= maxRays-1)
     {
       if(rnd_num <= albedo_whole)
         rnd_num /= albedo_whole;  // make first RND fall into the [0..1)
     }

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

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

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

  int indexRefracted = carray.GetUsedItems();

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

  int indexDiffuse = carray.GetUsedItems() - indexRefracted;

  // The number of rays for specular interaction
  //----------------------------------------------------------------------
  if (numSpecular)
    numRays += SpawnIdealReflectedRay(ks_local, carray, hitInfo);

  int indexMax = carray.GetUsedItems();

  // Scaling for distributed ray tracing
  /*
  for(int i = 0; i < indexMax; i++) {
    CContributionItem* ci = carray.GetItem(i);
    float mult;
    if(i < indexRefracted) {
      mult = akr;
    }
    else {
      if(i < (indexRefracted + indexDiffuse)) {
        mult = akd;
        ci->weigth *= (1.0/((float)numDiffuse));
      }
      else {
        mult = aks;
      }
    }
    mult /= albedo_whole;
    ci->prob *= mult;
    ci->weigth *= (1.0/mult);
  } // for*/
  samplingInfo.IncDimension(2);

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


// Spawn the ray deterministically in mirror-reflection direction //========================================================================
int
CSurfacePhongPhysicalDR::SpawnIdealReflectedRay(
  const CValues        &ks_local,
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo) const
{
  CContributionItem *ctr = carray.AllocateNext();
  ctr->weigth.SetRepresentation(CValues::EE_Weigth3);

  const CVector3D &reflectedDir = hitInfo.GetReflectedDir();

  // set the reflected ray
  ctr->ray.Init(hitInfo.GetPoint(), reflectedDir, REFLECTED_RAY, true);

  // the probability density of shooting ray to this direction
  ctr->prob = 1.0;
  ctr->diffuseness = 0.0;

  // this is BRDF for given input and output direction multiplied
  // by cosine between the incoming direction and normal
  ctr->weigth = ks_local;

  if(_maxPhong || _retroReflective)
  {
    // max-Phong model [3]
    float incEyeCos  = fabs(hitInfo.GetDotProd());
    float lightCos = fabs(DotProd(hitInfo.GetNormal(), reflectedDir));
    float max_cos = Max(incEyeCos, lightCos);
    ctr->weigth *= (1.0/max_cos);
  }
  // set the flags concerning ray type
  ctr->SetFlags(CContributionItem::EE_ReflectedRay);
  ctr->interaction = EE_REFLECTED_PERFECT;

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

  return 1; // number of rays in this case is one
}

// Diffuse part for stochastic sampling
//========================================================================
int
CSurfacePhongPhysicalDR::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();

  // If light source/viewer is under the surface
  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

  // number of rays created
  int num = 0;

  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    // TODO
    //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, ray.dir) / pdf =
    //         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
    //         = pd
    // This means that we should not multiply the weight by the cosine
    // of the angle between the incoming light direction and normal
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = kd_local;

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

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

const char* 
CSurfacePhongPhysicalDR::GetHelp() const
{
	return BrdfStrings::KPhongPhysicalDRDesc;
}

__END_BRDFLIB_SOURCE

