// ===================================================================
// $Id: reflsurf.cpp,v 1.8 2003/11/04 17:11:11 havran Exp $
//
// reflsurf.cpp
//     Implementation basic reflectors.
//
// 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 Jaroslav Krivanek (oct 2000) under supervision
// of Vlastimil Havran.
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


// BRDFLIB headers
#include "config.h"
#include "reflsurf.h"
#include "rnd.h"
#include "texture.h"

__BEGIN_BRDFLIB_SOURCE


/****************************************************************************
 *
 *
 *             Ideal diffuse reflector
 *
 *
 ***************************************************************************/

int
CIdealDiffuseReflector::
SpawnStochasticRaysStratify(CContributionArray  &carray,
			    const CHitPointInfo &hitInfo,
			    EDistributionType   distType,
			    CSamplingInfo       &samplingInfo,
			    const CVector3D &N) const
{
  ERayType rayType =
    (distType == EE_REFLECTED_DIFFUSE) ?
    REFLECTED_RAY :
    REFRACTED_RAY;

  // build the coordinate system
  CVector3D U = ArbitraryNormal(N);     // U is perpendicular to N
  CVector3D V = CrossProd(N,U);         // orthonormal base with N and U
  CVector3D dirN;

  int maxRays = samplingInfo.GetMaxSamples();
  
//#define _STRATIFY
#ifdef _STRATIFY
  // original code with explicit stratification (ordinary stochastic
  // importance sampling, using Monte Carlo method)
  int   thetaInt, phiInt;
  float r1w, r2w;
  if (maxRays==1) {
    thetaInt = phiInt  = 1;
    r1w      = r2w     = 1.0;    
  }
  else {
    thetaInt = (int)(sqrt((double)maxRays / 2.0));
    if (thetaInt < 1)
      thetaInt = 1;
    phiInt = (int)((float)maxRays/(float)thetaInt);
    r1w = 1.0 / (float)phiInt;
    r2w = 1.0 / (float)thetaInt;
  }
  float r1 = 0.0, r2 = 0.0;

  CValues kd_local( kd(hitInfo) * (1.0f/ (thetaInt*phiInt)) );

  for(int ir1 = phiInt; ir1-- ; r1 += r1w, r2 = 0 ) {
    for(int ir2 = thetaInt ; ir2-- ; r2 += r2w) {
      CContributionItem *contrib = carray.AllocateNext();
      // sample the hemisphere in canonical coordinates
      // with density cos(Phi)/M_PI - in L. Szirmay-Kalos script it
      // is explained why, pp. 59.
      // Sets the probability for generated direction
      contrib->prob = RND::hemisph_cos(dirN,
				       CVector2D(r1 + RND::rand0bf(r1w),
						 r2 + RND::rand0bf(r2w)));
      // transform from local coordinates
      CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, rayType);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      // type of interaction
      contrib->interaction = distType;
      // diffuseness of the bounce
      contrib->diffuseness = 1.0;
      // Evaluate the value of BSDF for given input and generated
      // output direction
      /* weight = BSDF * cos(Light, N) / PDF =
	        = pd / M_PI * cos(Light, N) / (cos(Light, N)/M_PI) = pd    */
      contrib->weigth = kd_local;
    } // for ir1
  } // for ir2

  // return number of samples
  return phiInt * thetaInt;

#else
  // without explicit stratification
  CValues kd_local( kd(hitInfo) );

  for(int i = 0; i < maxRays ; i++) {
    CContributionItem *contrib = carray.AllocateNext();

    if (i)
      // increase the seed not to generate the same random sample again
      samplingInfo.IncSeed();

    // sample the hemisphere in canonical coordinates with density 
    // cos(Phi)/M_PI (see [1], pp. 59).
    // 
    // Sets the probability for generated direction
    CVector2D vv(samplingInfo.Sample(0), samplingInfo.Sample(1));

    contrib->prob = RND::hemisph_cos(dirN, vv);

    // transform from local coordinates
    CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
    // set the reflected ray, direction is normalized here
    contrib->ray.Init(hitInfo.GetPoint(), dirRes, rayType);
    // set the flags concerning ray type
    contrib->SetFlags(CContributionItem::EE_ReflectedRay);
    // type of interaction
    contrib->interaction = distType;
    // diffuseness of the bounce
    contrib->diffuseness = 1.0;
    // Evaluate the value of BSDF for given input and generated
    // output direction
    /* weight = BSDF * cos(Light, N) / PDF =
              = pd / M_PI * cos(Light, N) / (cos(Light, N)/M_PI) = pd    */
    contrib->weigth = kd_local ;
  } // for i

  return maxRays;
#endif // _STRATIFY

}


// Gather all parts (only diffuse here) for stochastic sampling together
//========================================================================
int
CIdealDiffuseReflector::SpawnStochasticRays(CContributionArray  &carray,
					    const CHitPointInfo &hitInfo,
					    CSamplingInfo       &samplingInfo,
					    EDistributionType distType) const
{
  int maxRays = samplingInfo.GetMaxSamples();
  assert(_kd);

  if ( !(distType & EE_REFLECTED_DIFFUSE) ||
       (maxRays < 1) )
    return 0; // no rays

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

  // build coordinate system
  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;

  // spawn the rays using stratification
  int spawned =
    SpawnStochasticRaysStratify(carray, hitInfo, EE_REFLECTED_DIFFUSE,
				samplingInfo, N);

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

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

#ifdef __TRAVERSAL_STATISTICS
  world->_statistics->reflectedRays += spawned;
#endif

  return spawned;
}

// Return default color for the surface at a given point
//========================================================================
bool
CIdealDiffuseReflector::GetDefaultColor(CValues &defaultColor,
					const CHitPointInfo &hitInfo) const
{
  // let us consider _diffuse color to be the default one
  if (_kd) {
    hitInfo.GetTextureValue(defaultColor, *_kd, CMappingBase::EE_KD);
    return false;
  }

  defaultColor.Set(0.0);
  return true;
}

const char* 
CIdealDiffuseReflector::GetHelp() const
{
	return BrdfStrings::KIdealDiffuseReflDesc;
}
//========================================================================

int 
CIdealDiffuseReflector::GetNumParameters() const
{
	return 1; // diffuse
}
//========================================================================

CValues
CIdealDiffuseReflector::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	switch (parameterIndex)
	{
		case 0 :
			_kd->GetValue(c,CHitPointInfo::constMapInfo);
			break;
	}
	return c;
}

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

CValues
CIdealDiffuseReflector::GetParameterConstraints(int parameterIndex) const
{
	CValues c(CValues::EE_WeigthN);
	switch (parameterIndex)
	{
		case 0 :
			_setBasicConstraints(EE_ColorRGB,c);
			break;
	}
	return c;
}
//========================================================================

const char* 
CIdealDiffuseReflector::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KDiffuseDesc;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CIdealDiffuseReflector::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 : // diffuse texture
			return EE_ColorRGB;
		default :
			return EE_Unknown;
	}
}
//========================================================================

const char*
CIdealDiffuseReflector::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 : // diffuse texture
			return BrdfStrings::KDiffuseTitle;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CIdealDiffuseReflector::SetParameter(int parameterIndex, const CValues& value)
{
	CConstantTexture* newTex = new CConstantTexture(value);
	switch (parameterIndex)
	{
		case 0 : // diffuse texture
			SetKd(newTex);
			break;
	}
}
//========================================================================


/****************************************************************************
 *
 *
 *             Ideal mirror reflector
 *
 *
 ***************************************************************************/
// This type of surfaces always produce only one (deterministic) ray -
// ideally reflected ray.
//========================================================================
int
CIdealMirrorReflector::SpawnDeterministicRays(
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  EDistributionType   distributionType) const
{
  if ( !(distributionType & EE_REFLECTED_MIRROR) )
    return 0;

  CContributionItem *contrib = carray.AllocateNext();

  hitInfo.GetNormal();
  // set the reflected ray, direction is normalized here
  contrib->ray.Init(hitInfo.GetPoint(), hitInfo.GetReflectedDir() ,
		    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.0;
  // the reflection probabity is 1 for this type of reflection
  contrib->prob = 1.0;
  // Get albedo of ideal mirror
  CValues km_local(CValues::EE_Weigth3);
  hitInfo.GetTextureValue(km_local, *_km, CMappingBase::EE_KM);
  
  // The weight corresponds to BSDF multiplied by cosine.
  // where cosine = cos_out = cos(N,Light) = cos_relflected_ray =>
  //
  //    BRDF * cos(normal, ray.dir)/pdf =
  //                     = [ps*Dirac()/cos(normal, refl_ray.dir)]*
  //                       cos(normal, refl_ray.dir)/1
  //                     = ps
  //  23/11/03 K.P.
  contrib->weigth = km_local; 
#ifdef __TRAVERSAL_STATISTICS
  world->_statistics->IncReflectedRays();
#endif
  return 1;
}

// Return default color for the surface at a given point
//========================================================================
bool
CIdealMirrorReflector::GetDefaultColor(CValues &defaultColor,
				       const CHitPointInfo &hitInfo) const
{
  // let us consider _km color to be the default one
  if (_km) {
    hitInfo.GetTextureValue(defaultColor, *_km, CMappingBase::EE_KM);
    return false;
  }

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

const char* 
CIdealMirrorReflector::GetHelp() const
{
	return BrdfStrings::KIdealMirrorReflDesc;
}
//========================================================================

int 
CIdealMirrorReflector::GetNumParameters() const
{
	return 1; // albedo
}
//========================================================================

CValues
CIdealMirrorReflector::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	switch (parameterIndex)
	{
		case 0 :
			_km->GetValue(c,CHitPointInfo::constMapInfo);
			break;
	}
	return c;
}
//========================================================================

CValues
CIdealMirrorReflector::GetParameterConstraints(int parameterIndex) const
{
	CValues c(CValues::EE_WeigthN);
	switch (parameterIndex)
	{
		case 0 :
			_setBasicConstraints(EE_ColorRGB,c);
			break;
	}
	return c;
}
//========================================================================

const char* 
CIdealMirrorReflector::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KSpecularDesc;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CIdealMirrorReflector::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 : // diffuse texture
			return EE_RGB_Components;
		default :
			return EE_Unknown;
	}
}
//========================================================================

const char*
CIdealMirrorReflector::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 : // diffuse texture
			return BrdfStrings::KAlbedoTitle;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CIdealMirrorReflector::SetParameter(int parameterIndex, const CValues& value)
{
	CConstantTexture* newTex = new CConstantTexture(value);
	switch (parameterIndex)
	{
		case 0 : // volume texture
			SetKm(newTex);
			break;
	}
}
//========================================================================



/****************************************************************************
 *
 *
 *             Ideal metal reflector
 *
 *
 ***************************************************************************/
// This type of surfaces always produce only one (deterministic) ray -
// ideally reflected ray.
//========================================================================
int
CIdealMetal::SpawnDeterministicRays(
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  EDistributionType   distributionType) const
{
  if ( !(distributionType & EE_REFLECTED_MIRROR) )
    return 0;

  CContributionItem *contrib = carray.AllocateNext();

  hitInfo.GetNormal();
  // set the reflected ray, direction is normalized here
  contrib->ray.Init(hitInfo.GetPoint(), hitInfo.GetReflectedDir() ,
		    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.0;
  // the reflection probabity is 1 for this type of reflection
  contrib->prob = 1.0;
  
  CValues km_local(CValues::EE_Weigth3);

  // Check whether it is Fresnel texture
  bool use_fresnel = _km->IsFresnel();

  if(use_fresnel) 
  {
    hitInfo.GetTextureValue(km_local, *_km, CMappingBase::EE_KFM);
  }
  else
  {
	WARNING_INFO("CIdealMetal::SpawnDeterministicRays() :\n\
albedo for CIdealMetal must be of CFresnelTexture type");
    return 0;
  }
  // The weight corresponds to BSDF multiplied by cosine.
  // where cosine = cos_out = cos(N,Light) = cos_relflected_ray =>
  //
  //    BRDF * cos(normal, ray.dir)/pdf =
  //                     = [ps*Dirac()/cos(normal, refl_ray.dir)]*
  //                       cos(normal, refl_ray.dir)/1
  //                     = ps
  //  23/11/03 K.P.
  contrib->weigth = km_local; 
  
#ifdef __TRAVERSAL_STATISTICS
  world->_statistics->IncReflectedRays();
#endif
  return 1;
}



// Return default color for the surface at a given point
//========================================================================
bool
CIdealMetal::GetDefaultColor(CValues &defaultColor,
				       const CHitPointInfo &hitInfo) const
{
  // let us consider _kr color to be the default one
  if (_km) {
    bool use_fresnel = _km->IsFresnel();

    // return value according to Fresnel formulae (EE_KFM mapping) if Fresnel
    // texture is used and according to EE_KM mapping otherwise
    if(use_fresnel)
      hitInfo.GetTextureValue(defaultColor, *_km, CMappingBase::EE_KFM);
    else
      hitInfo.GetTextureValue(defaultColor, *_km, CMappingBase::EE_KM);

    return false;
  }
  defaultColor.Set(0.0);
  return true;
}

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

const char* 
CIdealMetal::GetHelp() const
{
	return BrdfStrings::KIdealMirrorReflDesc;
}
//========================================================================

int 
CIdealMetal::GetNumParameters() const
{
	return 1; // Fresnel texture
}
//========================================================================

CValues
CIdealMetal::GetParameter(int parameterIndex) const
{
	CValues c(CValues::EE_WeigthN);
	if (!_km)
	{
		return c;
	}
	switch (parameterIndex)
	{
		case 0 :
			CFresnelTexture* tex = static_cast<CFresnelTexture*>(_km);
			CValues aux(CValues::EE_Weigth3);
			tex->GetColor(aux);
			c.Set(0,aux.Get(0));
			c.Set(1,aux.Get(1));
			c.Set(2,aux.Get(2));
			tex->GetIndexR(aux);
			c.Set(3,aux.Get(0));
			tex->GetExtinction(aux);
			c.Set(4,aux.Get(0));
			bool schlick;
			tex->GetApproxSchlick(schlick);
			c.Set(5,schlick);
			break;
	}
	return c;
}
//========================================================================

CValues
CIdealMetal::GetParameterConstraints(int parameterIndex) const
{
	CValues c(CValues::EE_WeigthN);
	switch (parameterIndex)
	{
		case 0 :
			_setBasicConstraints(EE_Fresnel,c);
			break;
	}
	return c;
}
//========================================================================

const char* 
CIdealMetal::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KFresnelDesc;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CIdealMetal::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return EE_Fresnel;
		default :
			return EE_Unknown;
	}
}
//========================================================================

const char*
CIdealMetal::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KFresnelTitle;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

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


/****************************************************************************
 *
 *
 *             Specular reflector
 *
 *
 ***************************************************************************/
// Gather all parts (only diffuse here) for stochastic sampling together
//========================================================================
int
CSpecularReflector::SpawnStochasticRays(CContributionArray  &carray,
					const CHitPointInfo &hitInfo,
					CSamplingInfo       &samplingInfo,
					EDistributionType   distType) const
{
  int maxRays = samplingInfo.GetMaxSamples();
  if(!(distType & EE_REFLECTED_SPECULAR))
    return 0;

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

  CVector3D dirEye = (hitInfo.GetRay())->GetDir();
  const CVector3D &R = hitInfo.GetReflectedDir();
  CVector3D U = ArbitraryNormal(R);  // U is perpendicular to R
  CVector3D V = CrossProd(R,U);      // ortonormal base with R and U

  CValues kslocal = ks(hitInfo), nlocal = n(hitInfo);
  // diffuseness of the bounce
  float diffuseness = _nToDiffuseness(nlocal.Sum() / nlocal.GetCount());

  // TODO: make possible difraction for Phong model (=> different values
  // of each of RGB component for shininess)
  assert(nlocal.GetCount() == 1);

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

#ifdef _STRATIFY
  // original code with explicit stratification (ordinary stochastic
  // importance sampling, using Monte Carlo method)
  int   thetaInt, phiInt;
  float r1w, r2w;
  CVector3D dirN;
  if (maxRays==1) {
    thetaInt = phiInt  = 1;
    r1w      = r2w     = 1.0;
  }
  else {
    thetaInt = (int)(sqrt((double)maxRays/2.0));
    if (thetaInt < 1)
      thetaInt = 1;
    phiInt = (int)((float)maxRays/(float)thetaInt);
    r1w = 1.0 / phiInt;
    r2w = 1.0 / thetaInt;
  }

  // the inverse proportional number of samples
  float ratio = 1.0/(thetaInt * phiInt);
  float r1 = 0, r2 = 0;

  // number of rays created

  for(int ir1 = phiInt; ir1-- ; r1 += r1w, r2 = 0 ) {
    for(int ir2 = thetaInt ; ir2-- ; r2 += r2w, num++ ) {
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      // generate random tuplet
      float lightCos = -1;
      CVector3D dirRes;
      float cos_ang_V_R;

      while(lightCos < 0) {
        if (num) {
	  // increase the seed not the generate the same sample again
          samplingInfo.IncSeed();
        }

        // Generate a sample in 2D, range [0;1)x[0;1)
        CVector2D vv(samplingInfo.Sample(0), samplingInfo.Sample(1));
        // stratified sampling !!!!
        CVector2D v(r2 + vv.xx * r2w , r1 + vv.yy * r1w);

        // generate canonical direction
        cos_ang_V_R  = pow(float(v.xx), float(1.0/(n + 1.0)) );
        float cos2_ang_V_R = sqr(cos_ang_V_R);
        float sin_ang_V_R = sqrt(1.0 - cos2_ang_V_R);

        dirN.SetValue( cos( 2.0 * M_PI * v.yy ) * sin_ang_V_R,
		       sin( 2.0 * M_PI * v.yy ) * sin_ang_V_R,
		       cos_ang_V_R );

        // transform the canonical direction to global coords
        dirRes = dirN.x * U + dirN.y * V + dirN.z * R ;

	// rejected samples if dirRes appears to be under the surface
	lightCos = DotProd(N , dirRes );
      } // while
      // compute the probability density of reflection in this direction
      // 	pdf = (n+1)/(2*M_PI) * cos(Mirror, Eye)^n
      contrib->prob = multi * pow(cos_ang_V_R, n); 

#ifdef MAX_PHONG
      // if max-Phong model [4] is used divide by cosine
      float cos_in  = fabs(hitInfo.GetDotProd());
      float cos_out = fabs(DotProd(N, dirRes));
      float max_cos = Max(cos_in, cos_out);
#endif

      // compute the BRDF for this direction
      contrib->weigth = kslocal * (ratio * lightCos / multi
#ifdef MAX_PHONG
				     / max_cos
#endif
		        );

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

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

    } // for ir2
  } // for ir1
  samplingInfo.IncDimension(2);

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

  return num;

#else
  // without explicit stratification
  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 new one (in other case sample will be 
    // considered as adsorbed, 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
       }*/

       CVector2D h(samplingInfo.Sample(0), samplingInfo.Sample(1));
 
       // See [3] 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 * V + localDir.y * U + localDir.z * R;

       // If cosine lobe is computed according to the Blinn [5], 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 * (-dirEye + dirRes);     // Blinn-Phong
         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
    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 weight
    contrib->weigth = kslocal * multi2;
    /* weight = value_of_BRDF / probability_density_of_generated direction
       or more precisely:
       weight = 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. */

    if (_maxPhong) {
      // if max-Phong model is used divide by cosine
      float cos_in  = fabs(hitInfo.GetDotProd());
      float cos_out = fabs(DotProd(N, dirRes));
      float max_cos = Max(cos_in, cos_out);
      contrib->weigth *= (1.0f/max_cos);  // (*)
    }
    contrib->weigth *= lightCos;  // (*)

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

#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif
    samplingInfo.IncDimension(2);

    // recover the seed for the next dimension to be sampled
    samplingInfo.RecoverSeed();
    // the number of valid samples taken
  } // for i

  assert(num == maxRays);
  return num;
#endif // _STRATIFY
}

// Direct lighting computation
//========================================================================
bool
CSpecularReflector::EvaluateBSDF(CValues &res,
				 const CHitPointInfo &hitInfo,
				 const CVector3D &dirLit,
				 EDistributionType distType) const
{
  if (distType & EE_REFLECTED_GLOSSY) {
    CVector3D normal = hitInfo.GetNormal();
    CVector3D dirEye = (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 >= -1e-05 && incEyeCos >= -1e-05)
    {
      float cosR_L;
      // If cosine lobe is computed according to the Blinn article, in
      // SpawnStochasticSpecularRay function this fact should be also
      // taken into account
      if (!_normalsBlinn) {
        // evaluation of the specular component according to Phong
        cosR_L = DotProd(hitInfo.GetReflectedDir(), dirLit);
		if (cosR_L < WTINY) cosR_L = WTINY;
      }
      else {
        // evaluation of the specular component according to Blinn [5]
        // According to Blinn's article Halfvector - is the half vector
        // between the incoming and outgoing directions K.P. 07/12/03
        CVector3D hVector = -dirEye + dirLit;
        hVector.Normalize();
        cosR_L = DotProd(hVector, normal);
      }

      if (cosR_L >= 0.0) {
        CValues nlocal = n(hitInfo);
        assert(nlocal.GetCount() == 1);

        float nn = nlocal.Get(0);
        res = ks(hitInfo) * (pow(cosR_L, nn) * (nn + 2.0) / (2.0*M_PI) );
        if (_maxPhong) {
           // if max-Phong model is used divide by cosine
           float max_cos = Max(incEyeCos, lightCos);
           res *= (1.0f/max_cos);  // (*)
        }
        res *= lightCos;
        return true; // OK, value is set
      }
	} // lightCos
  } // type
  
  return false; // value is not set
}


//========================================================================
float
CSpecularReflector::GetAverageDiffuseness(void) const
{
  // texture needs mapping - not yet implemented
  assert(_n->GetType() == CMappingBase::EE_None);
  CValues c;
  _n->GetValue(c, CHitPointInfo::constMapInfo);
  float n = c.MeanValue();

  return _nToDiffuseness(n);
}


//========================================================================
float
CSpecularReflector::GetDiffuseness(const CHitPointInfo &hitInfo,
				   EDistributionType distType) const
{
  if (distType & EE_REFLECTED_GLOSSY) {
    CValues c;
    hitInfo.GetTextureValue(c, *_n, CMappingBase::EE_KSI);
    float n = c.MeanValue();
    
    return _nToDiffuseness(n);
  }

  return 0; // not possible to evaluate
}


//========================================================================
bool
CSpecularReflector::GetDefaultColor(CValues &defaultColor,
				    const CHitPointInfo &hitInfo) const
{
  // let us consider _ks color to be the default one
  if (_ks) {
    hitInfo.GetTextureValue(defaultColor, *_ks, CMappingBase::EE_KS);
    return false;
  }

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

const char* 
CSpecularReflector::GetHelp() const
{
	return BrdfStrings::KSpecularReflectorDesc;
}
//========================================================================

int 
CSpecularReflector::GetNumParameters() const
{
	return 4; // specular, shininess, maxPhong, normalsBlinn
}
//========================================================================

CValues
CSpecularReflector::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	switch (parameterIndex)
	{
		case 0 :
			_ks->GetValue(c,CHitPointInfo::constMapInfo);
			break;
		case 1 :
			_n->GetValue(c,CHitPointInfo::constMapInfo);
			break;
		case 2 :
			c.Set(_maxPhong);
			break;
		case 3 :
			c.Set(_normalsBlinn);
			break;
	}
	return c;
}
//========================================================================

CValues
CSpecularReflector::GetParameterConstraints(int parameterIndex) const
{
	CValues c(3);
	switch (parameterIndex)
	{
		case 0 :
			_setBasicConstraints(EE_ColorRGB,c);
			break;
		case 1 :
			c.Set(0,0);
			c.Set(1,MAX_SHININESS);
			c.Set(2,1);
			break;
		case 2 :
		case 3 :
			_setBasicConstraints(EE_Boolean,c);
			break;
	}
	return c;
}
//========================================================================

const char* 
CSpecularReflector::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KSpecularDesc;
		case 1 :
			return BrdfStrings::KShininessDesc;
		case 2 :
			return BrdfStrings::KMaxPhongDesc;
		case 3 :
			return BrdfStrings::KBlinnNormalsDesc;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CSpecularReflector::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 : //specular
			return EE_ColorRGB;
		case 1 : //shininess
			return EE_Float;
		case 2 : //max phong
			return EE_Boolean;
		case 3 : //normals blinn
			return EE_Boolean;
		default :
			return EE_Unknown;
	}
}
//========================================================================

const char*
CSpecularReflector::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
			return BrdfStrings::KSpecularTitle;
		case 1 :
			return BrdfStrings::KShininessTitle;
		case 2 :
			return BrdfStrings::KMaxPhongTitle;
		case 3 :
			return BrdfStrings::KBlinnNormalsTitle;
		default :
			return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CSpecularReflector::SetParameter(int parameterIndex, const CValues& value)
{
	CConstantTexture* newTex = 0;
	bool boolValue = false;
	if ((parameterIndex >=0) && (parameterIndex < 2))
	{
		newTex = new CConstantTexture(value);
	}
	if ((parameterIndex >= 2) && (parameterIndex < 4))
	{
		boolValue = static_cast<bool>(value[0]);
	}
	
	switch (parameterIndex)
	{
		case 0 :{
					SetSpecular(newTex);
					break;
				}
		case 1 :{
					SetShininess(newTex);
					break;
				}
		case 2 :{
					SetMaxPhong(boolValue);
					break;
				}
		case 3 :{
					SetBlinnNormals(boolValue);
					break;
				}
	}
}
//========================================================================

__END_BRDFLIB_SOURCE

