// ===================================================================
// $Id: mgfsurf.cpp,v 1.4 2003/10/30 16:36:36 havran Exp $
//
// mgfsurf.cpp
//     Implementation of MGF surface.
//
// 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 Vlastimil Havran, 2002.
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.


/// !!!!!!!!!!!!!!!!!!!
// THIS FILE IS NOT TESTED AND SHOULD BE CHECKED AND VERIFIED
// Note by Vlastimil Havran, November 21, 2003 !!!!
// !!!!!!!!!!!!!!!!!!!!!

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

__BEGIN_BRDFLIB_SOURCE

// This defines if max-Phong or normal Phong shading model is used
//#define MAX_PHONG

// When this macro is used, classical normals as described originally
// by Phong are used to evaluate the specular component. If not define,
// Blinn variant with H vector is used
#define PHONG_CLASSIC_NORMALS

/**********************************************************************
 *
 *             CSurfaceMGF methods
 *
 **********************************************************************/


CSurface *
CSurfaceMGF::Dup() const
{
  CSurfaceMGF *ret = new CSurfaceMGF();

  if (_diffuseReflection)
    ret->_diffuseReflection = _diffuseReflection;
  if (_specularReflection)
    ret->_specularReflection = _specularReflection;
  if (_specnumReflection)
    ret->_specnumReflection = _specnumReflection;
  if (_diffuseTransmission)
    ret->_diffuseTransmission = _diffuseTransmission;
  if (_specularTransmission)
    ret->_specularTransmission = _specularTransmission;
  if (_specnumTransmission)
    ret->_specnumTransmission = _specnumTransmission;
  if (_iorReal)
    ret->_iorReal = _iorReal;
  if (_iorImag)
    ret->_iorImag = _iorImag;

  ret->_supportedType = _supportedType;
  return ret;
}


// For given input direction and distribution type
// computes value of albedo to variable 'albedo'
bool
CSurfaceMGF::GetAlbedoMax(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

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

    if ((distType & EE_REFLECTED_SPECULAR) && _specularReflection) {
      hitInfo.GetTextureValue(d, *_specularReflection, CMappingBase::EE_KS);
      c += d * fabs(hitInfo.GetDotProd());
      retval = true;
    }
  }
  else {
    if ((distType & EE_REFRACTED_DIFFUSE) && _diffuseTransmission) {
      hitInfo.GetTextureValue(d, *_diffuseTransmission, CMappingBase::EE_KDT);
      c += d;
      retval = true;
    }

    if ((distType & EE_REFRACTED_SPECULAR) && _specularTransmission) {
      hitInfo.GetTextureValue(d, *_specularTransmission, CMappingBase::EE_KST);
      c += d * fabs(hitInfo.GetDotProd());
      retval = true;
    }
  }

  // test if the albedo is smaller than 1.0 -- physically plausible ?
  assert(c.Max() <= 1.0);

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

  return false; // black body - no albedo
}

// !!!!!!!!!!!!!!!!!!!!!!!!!!!!! TODO
// For given input direction and distribution type
// computes value of albedo to variable 'albedo'
bool
CSurfaceMGF::GetAlbedo(
  CValues &albedo,
  const CHitPointInfo &hitInfo,
  EDistributionType   distType) const
{
  CValues c(CValues::EE_Weigth3); // zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

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

    if ((distType & EE_REFLECTED_SPECULAR) && _specularReflection) {
      hitInfo.GetTextureValue(d, *_specularReflection, CMappingBase::EE_KS);
      c += d * fabs(hitInfo.GetDotProd());
      retval = true;
    }
  }
  else {
    if ((distType & EE_REFRACTED_DIFFUSE) && _diffuseTransmission) {
      hitInfo.GetTextureValue(d, *_diffuseTransmission, CMappingBase::EE_KDT);
      c += d;
      retval = true;
    }

    if ((distType & EE_REFRACTED_SPECULAR) && _specularTransmission) {
      hitInfo.GetTextureValue(d, *_specularTransmission, CMappingBase::EE_KST);
      c += d * fabs(hitInfo.GetDotProd());
      retval = true;
    }
  }

  // test if the albedo is smaller than 1.0 -- physically plausible ?
  assert(c.Max() <= 1.0);

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

  return false; // black body - no albedo
}


float
CSurfaceMGF::GetAverageDiffuseness(void) const
{
  float diffuseRefl = 0, specularRefl = 0, nRefl = 0,
        diffuseTrans = 0, specularTrans = 0, nTrans = 0;

  // the auxiliary variable
  CValues c;
  
  if (_diffuseReflection) {
    // ?????????? why this assert is required - works only for
    // constant texture ?
    assert(_diffuseReflection->GetType() == CMappingBase::EE_None);
    _diffuseReflection->GetValue(c, CHitPointInfo::constMapInfo);
    diffuseRefl = c.Luminance();
  }
  
  if ((_specularReflection) && (_specnumReflection)) {
    // ????????? can be EE_KS a EE_KSI
    assert(_specularReflection->GetType()==CMappingBase::EE_None);
    assert(_specnumReflection->GetType()==CMappingBase::EE_None);
    
    _specularReflection->GetValue(c, CHitPointInfo::constMapInfo);
    specularRefl = c.Luminance();
    
    _specnumReflection->GetValue(c, CHitPointInfo::constMapInfo);
    nRefl = c.Luminance();
  }

  if (_diffuseTransmission) {
    // ????????
    assert(_diffuseTransmission->GetType() == CMappingBase::EE_None);
    _diffuseTransmission->GetValue(c, CHitPointInfo::constMapInfo);
    diffuseTrans = c.Luminance();
  }
  
  if (_specularTransmission) {
    // ????????? can be EE_KST a EE_KSIT
    assert(_specularTransmission->GetType()==CMappingBase::EE_None);
    assert(_specnumTransmission->GetType()==CMappingBase::EE_None);
    
    _specularTransmission->GetValue(c, CHitPointInfo::constMapInfo);
    specularTrans = c.Luminance();
    
    _specnumTransmission->GetValue(c, CHitPointInfo::constMapInfo);
    nTrans = c.Luminance();
  }
  
  // result is the weighted average
  float res =
    ( diffuseRefl + specularRefl * _nToDiffuseness(nRefl) +
      diffuseTrans + specularTrans * _nToDiffuseness(nTrans))
    / (diffuseRefl + specularRefl + diffuseTrans + specularTrans);

  return res;
}


// for given input direction compute the average diffuseness of the bounce
float
CSurfaceMGF::GetDiffuseness(const CHitPointInfo &hitInfo,
			    EDistributionType distType) const
{
  float diffuseRefl = 0, specularRefl = 0, nRefl = 0,
        diffuseTrans = 0, specularTrans = 0, nTrans = 0;

  // the auxiliary variable
  CValues c;
  
  if ((distType & EE_REFLECTED_DIFFUSE) && _diffuseReflection) {
    hitInfo.GetTextureValue(c, *_diffuseReflection, CMappingBase::EE_KD);
    diffuseRefl = c.Luminance();
  }
  
  if ((distType & EE_REFLECTED_SPECULAR) && _specularReflection) {
    hitInfo.GetTextureValue(c, *_specularReflection, CMappingBase::EE_KS);
    specularRefl = c.Luminance();
    
    hitInfo.GetTextureValue(c, *_specnumReflection, CMappingBase::EE_KSI);
    nRefl = c.Luminance();
  }

  if ((distType & EE_REFRACTED_DIFFUSE) && _diffuseTransmission) {
    hitInfo.GetTextureValue(c, *_diffuseTransmission, CMappingBase::EE_KDT);
    diffuseTrans = c.Luminance();
  }
  
  if ((distType & EE_REFRACTED_SPECULAR) &&_specularTransmission) {
    hitInfo.GetTextureValue(c, *_specularTransmission, CMappingBase::EE_KST);
    specularTrans = c.Luminance();

    hitInfo.GetTextureValue(c, *_specnumTransmission, CMappingBase::EE_KSIT);
    nTrans = c.Luminance();
  }

  float denom = (diffuseRefl + specularRefl + diffuseTrans + specularTrans);
  
  // result is the weighted average
  float res = 0.0;

  if (denom > 0.0)
    res = ( diffuseRefl + specularRefl * _nToDiffuseness(nRefl) +
	    diffuseTrans + specularTrans * _nToDiffuseness(nTrans))
      / denom ;

  return res;
}

// The value of brdf is ks_local * cosn^n / max_cos.
int
CSurfaceMGF::SpawnStochasticSpecularRays(
  const CValues       &kslocal,
  const CValues       &shininess,
  const CVector3D     &direction,
  const bool          &transmissionFlag,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  // we would have to be spectral material to be able to handle
  // more specular indexes
  assert(shininess.GetCount() == 1);

  // the number of samples generated at the end of this function
  int num = 0;

  CVector3D N = hitInfo.GetNormal();
  if (DotProd(N, direction) < 0.0)
    N = -N;
  
  // Direction of the ray below corresponds to 
  // hitInfo.GetReflectedDir() or to ideally ReflectedRay.
  const CVector3D &R = direction;
  CVector3D U = ArbitraryNormal(R);  // U is perpendicular to R
  CVector3D V = CrossProd(R, U);     // ortonormal base with R and U
  CVector3D dirN;

  // 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);
  // this value is required for BRDF evaluation
  float multi2 = 1.0/multi * (n + 2.0)/(M_PI * 2.0);

  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 / phiInt;
    r2w = 1.0 / thetaInt;
    // compute the number of rays now
    maxRays = thetaInt * phiInt;
  }
  
  float ratio = 1.0 / (phiInt * thetaInt);
  float r1 = 0, r2 = 0;

  ERayType rayType = REFLECTED_RAY;

  if (transmissionFlag)
    rayType = REFRACTED_RAY;
  
  // loop in the phi angle
  for(int ir1 = phiInt; ir1-- ; r1 += r1w, r2 = 0 ) {
    // loop in the theta angle
    for(int ir2 = thetaInt ; ir2-- ; r2 += r2w ) {
      CContributionItem *contrib = carray.AllocateNext();
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);

      // generate random tuple
      CVector2D h(aux_rnd, samplingInfo.Sample(1));
      CVector2D v(r2 + h.xx * r2w , r1 + h.yy * r1w);
      
      // generate canonical direction
      float 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
      CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * R ;

      // check if it is valid direction
      if (DotProd( N , dirRes ) < 0.0) {
	// This can happen because the specular lobe can reach below
	// the plane of the surface. This is clear disadvantage of
	// Phong model.
#ifdef _DEBUG
	/* 
	   DEBUG << "CSurfaceMGF::SpawnStochasticSpecularRays: "
	   " reflection turned to refraction and rejected." << endl;
	*/
#endif
	carray.DeleteLastItem();
	continue;
      }

      // compute the probability density
      contrib->prob = multi *                     // (n+1)/2/M_PI 
	pow( cos_ang_V_R, n);                     // = cos^n(theta)

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

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

      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, rayType);

      // set the flags concerning ray type and type of interaction
      if (transmissionFlag) {
	contrib->SetFlags(CContributionItem::EE_ReflectedRay);
	contrib->interaction = EE_REFLECTED_GLOSSY;
      }
      else {
	contrib->SetFlags(CContributionItem::EE_RefractedRay);
	contrib->interaction = EE_REFRACTED_GLOSSY;
      }

      // diffuseness of the bounce
      contrib->diffuseness = diffuseness;  
      
#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif
      num++;

      // generate a new sample for the next time
      if (num < maxRays)
	aux_rnd = samplingInfo.Sample(0);
    } // for ir2
  } // for ir1

  return num;  
}


/*
  kd_local is the value got from the texture divided by PI - thus it is
  direcly used as the value of BRDF.
*/
int
CSurfaceMGF::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  const bool          &transmissionFlag,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  // for transmission we have to change the normal orientation
  if (transmissionFlag)
    N = -N;
  
  CVector3D U = ArbitraryNormal(N);     // U is perpendicular to N
  CVector3D V = CrossProd(N, U);        // orthonormal base with N and U

  CValues kd_weighted(kd_local * ( M_1_PI / maxRays));

  ERayType rayType = REFLECTED_RAY;
  if (transmissionFlag)
    rayType = REFRACTED_RAY;

  // for all rays
  for(int i = maxRays; i--; ) {
    // allocate next contribution
    CContributionItem *contrib = carray.AllocateNext();
    // generate the sample in the cannonical space
    CVector3D dirN;
    // the random number in 2D
    CVector2D vv(aux_rnd, samplingInfo.Sample(1));
    
    // probability density 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, rayType);

    // set the flags concerning ray type and type of interaction
    if (transmissionFlag) {
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      contrib->interaction = EE_REFLECTED_DIFFUSE;
    }
    else {
      contrib->SetFlags(CContributionItem::EE_RefractedRay);
      contrib->interaction = EE_REFRACTED_DIFFUSE;      
    }
    // diffuseness of the bounce
    contrib->diffuseness = 1.0;
    // Set the value of BSDF type
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    // the BRDF value for generated ray direction
    contrib->weigth = kd_weighted;

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

    // generate a new sample for the next time
    if (i > 1)
      aux_rnd = samplingInfo.Sample(0);
  } // for
  return maxRays;
}  

int
CSurfaceMGF::SpawnDeterministicRays(
  CContributionArray  &contribution,
  const CHitPointInfo &hitInfo,
  EDistributionType distributionType) const
{
// whether we allow deterministic reflection for ray tracing stuff
#define ALLOW_DETERMINISTIC_REFLECTION
  
#ifndef ALLOW_DETERMINISTIC_REFLECTION
  WARNING << "CSurfaceMGF should be used only in physically based renderers. "
	  << "For classical raytracing, use CSurfacePhongClassical. " << endl;
  return 0;
#else
  // This is from CSurfacePhongClassical - for testing ray tracing over MGF scenes
  int numRays = 0;

  CVector3D normal = hitInfo.GetNormal();
  
  // a reflected ray should be spawned ?
  if (_specularReflection && (distributionType & EE_REFLECTED) ) {
    // according to SPD and other source, one reflected ray should be spawned
    // when specular component is not zero
    CContributionItem *ctr = contribution.AllocateNext();
    ctr->weigth.SetRepresentation(CValues::EE_Weigth3);
    
    // compute the normal and ideally reflected ray 
    const CVector3D &dir = hitInfo.GetReflectedDir();
    
    // set the reflected ray 
    ctr->ray.Init(hitInfo.GetPoint(), dir, REFLECTED_RAY);

    // the probability density of shooting ray to this direction | coming from
    // given direction
    ctr->prob = 1.0;      
    
    // evaluating diffuse term
    if (_diffuseReflection) {
      // get value of diffuse coefficient for intersection point
      hitInfo.GetTextureValue(ctr->weigth, *_diffuseReflection,
			      CMappingBase::EE_KD);
    }
    else
      // the value is set to zero at beginning
      ctr->weigth.Set(0.0);
    
    // evaluating specular term
    CValues specularTerm(CValues::EE_Weigth3);
	
    // get value of specular coefficient for intersection point
    hitInfo.GetTextureValue(specularTerm, *_specularReflection,
			    CMappingBase::EE_KS);
	
    // specular term is not required to multiply by cos^N angle, since
    // the angle is always zero - deterministic ideal reflection.
	
    // now perform summing of diffuse and specular component
    ctr->weigth += specularTerm;

    // set the flags concerning ray type
    ctr->SetFlags(CContributionItem::EE_ReflectedRay);
    ctr->interaction = EE_REFLECTED_MIRROR;

    // does the reflection make sense
    const float epsilon = 1e-5;
    if (ctr->weigth.AbsSum() < epsilon ) {
      // reflection with zero contribution
      contribution.DeleteLastItem();
    }
    else {
#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncReflectedRays();
#endif // __TRAVERSAL_STATISTICS
      numRays++; // non-zero reflection
    }
  } // if(_specular)

  // a refracted ray should be spawned ?
  if ( (distributionType & EE_REFRACTED) &&
       (_specularTransmission) ) {
    // the item for refracted ray
    CContributionItem *ctr = contribution.AllocateNext();
    ctr->weigth.SetRepresentation(CValues::EE_Weigth3);
    
    // get value of refraction coefficient for intersection point
    hitInfo.GetTextureValue(ctr->weigth, *_specularTransmission,
			    CMappingBase::EE_KST);

    // direction of refracted ray to be computed
    CVector3D dir;
    
    // evaluating the index of refraction
    CValues indexOfRefractionReal(CValues::EE_Weigth1);
    CValues indexOfRefractionImag(CValues::EE_Weigth1);
    if ((_iorReal) && (_iorImag)) {
      // get value of index of refraction  coefficient for intersection point
      hitInfo.GetTextureValue(indexOfRefractionReal, *_iorReal, CMappingBase::EE_KRI);
      hitInfo.GetTextureValue(indexOfRefractionImag, *_iorImag, CMappingBase::EE_KRI);

      // currently, rays are not bended according to spectral wavelegth
      assert(indexOfRefractionReal.GetType() == CValues::EE_Weigth1);
      assert(indexOfRefractionImag.GetType() == CValues::EE_Weigth1);
      
      // here is the point to be modified for spectral-dependent effects
      float gamma, sqrterm;
      
      float dot = hitInfo.GetDotProd();
      float ior = sqrt( sqr(indexOfRefractionReal.Get(0)) +
		        sqr(indexOfRefractionImag.Get(0)) );
	
      if (dot < 0.0) {
	// from outside into the inside of object
	gamma = 1.0 / ior;
      }
      else {
	// from the inside to outside of object
	gamma = ior;
	dot = -dot;
	normal = -normal;
      }
      sqrterm = 1.0 - gamma * gamma * (1.0 - dot * dot);
      
      // Check for total internal reflection, do nothing if it applies.
      if (sqrterm > 0.0) {
	sqrterm = dot * gamma + sqrt(sqrterm);
	dir = -sqrterm * normal + hitInfo.GetRay()->GetDir() * gamma;
      }
      else  {
	// total internal reflection occured
	// DEBUG << "Total reflection in Phong" << endl;
	contribution.DeleteLastItem();
	return numRays;
      }
    }
    else { // _index
      // assuming the light is not bent at all.
      dir = hitInfo.GetRay()->GetDir();
    }

#ifdef _DEBUG
    float d = SqrMagnitude(dir);
    if ( (d < 0.99) ||
	 (d > 1.01) ) {
      WARNING << "direction in CSurfacePhong::SpawnDeterministicRay "
	      << " is not normalized\n";
    }
#endif // _DEBUG
    // compute the refracted ray
    ctr->ray.Init(hitInfo.GetPoint(), dir,
		  REFRACTED_RAY);
    
    // the probability that the ray will be shot in this direction
    ctr->prob = 1.0;
    
    // set the flags concerning ray type
    ctr->SetFlags(CContributionItem::EE_RefractedRay);
    // !!!!!!!! where is set ctr->weigth for refracted ray ????????
    if (ctr->weigth.AbsSum() < 1e-5 )
      contribution.DeleteLastItem(); // zero reflection component
    else {
#ifdef __TRAVERSAL_STATISTICS
      world->_statistics->IncRefractedRays();
#endif // __TRAVERSAL_STATISTICS
      numRays++; // non-zero reflection
    }
  }
  return numRays;
#endif // ALLOW_DETERMINISTIC_REFLECTION
}


// Main function shooting rays
int
CSurfaceMGF::SpawnStochasticRays(
  CContributionArray   &carray,
  const CHitPointInfo  &hitInfo,
  CSamplingInfo        &samplingInfo,
  EDistributionType     distType) const
{
  int maxRays = samplingInfo.GetMaxSamples();
  assert(maxRays > 0);

  // diffuse component for reflection
  CValues kd_local_refl(CValues::EE_Weigth3);
  // specular component for reflection
  CValues ks_local_refl(CValues::EE_Weigth3);

  // diffuse component for reflection
  CValues kd_local_trans(CValues::EE_Weigth3);
  // specular component for reflection
  CValues ks_local_trans(CValues::EE_Weigth3);
  
  // compute diffuse and specular albedos for reflection
  float akd = 0.0f, aks = 0.0f;
  // compute diffuse and specular albedos for reflection
  float akdt = 0.0f, akst = 0.0f;

  // Validate hitInfo.Dotprod() function for future use
  CVector3D N = hitInfo.GetNormal();

  // REFLECTION
  // get the values for diffuse and specular reflection to 'akd' and 'aks'
  if (_diffuseReflection && (distType & EE_REFLECTED_DIFFUSE)) {
    hitInfo.GetTextureValue(kd_local_refl, *_diffuseReflection,
			    CMappingBase::EE_KD);
    akd = kd_local_refl.Luminance();
  }
  if (_specularReflection && (distType & EE_REFLECTED_SPECULAR)) {
    hitInfo.GetTextureValue(ks_local_refl, *_specularReflection,
			    CMappingBase::EE_KS);
    aks = ks_local_refl.Luminance();
  }

  // TRANSMISSION
  // get the values for diffuse and specular transmission to 'akdt' and 'akst'
  if (_diffuseTransmission && (distType & EE_REFRACTED_DIFFUSE)) {
    hitInfo.GetTextureValue(kd_local_trans, *_diffuseTransmission,
			    CMappingBase::EE_KDT);
    akdt = kd_local_trans.Luminance();
  }
  if (_specularTransmission && (distType & EE_REFRACTED_SPECULAR)) {
    hitInfo.GetTextureValue(ks_local_trans, *_specularTransmission,
			    CMappingBase::EE_KST);
    aks = ks_local_trans.Luminance();
  }
  
  if ( (akd <= 1e-7) &&
       (aks <= 1e-7) &&
       (akdt <= 1e-7) &&
       (akst <= 1e-7) )
    return 0; // no rays -- this surface is black body

  //  DEBUG << (int)distType << endl;

  // now the akd, aks and akr is known

  // Compute the number of rays to reflect diffusely / specularly
  // or trasmitt diffusely / specular
  int numDiffuseRefl = 0;
  int numSpecularRefl = 0;
  int numDiffuseTrans = 0;
  int numSpecularTrans = 0;

  //  DEBUG << "akd: " << akd << " " << kd_local_refl << endl;
  //  DEBUG << "aks: " << aks << " " << ks_local_refl << endl;
  //  DEBUG << "akdt: " << akdt << " " << kd_local_trans << endl;
  //  DEBUG << "akst: " << akst << " " << ks_local_trans << endl;

  // perform kind of normalization to get the exact number of rays
  float scale = 1.0 / (akd + aks + akdt + akst);
  akd *= scale;
  aks *= scale;
  akdt *= scale;
  akst *= scale;
  
  int count = maxRays; 
  while (count--) {
    float rnd_num = samplingInfo.Sample(0);
    if (rnd_num < akd) {
      numDiffuseRefl++;
      SetAuxRnd(rnd_num, 0.0, akd);
      continue;
    }
    rnd_num -= akd;
    if (rnd_num < aks) {
      numSpecularRefl++;
      SetAuxRnd(rnd_num, 0.0, aks);
      continue;
    }
    rnd_num -= aks;
    if (rnd_num < akdt) {
      numDiffuseTrans++;
      SetAuxRnd(rnd_num, 0.0, akdt);
      continue;
    }
    SetAuxRnd(rnd_num, akdt, akdt+akst);
    numSpecularTrans++;
  } // while

  /*
  DEBUG << endl;
  DEBUG << "numDiffuseRefl:  " << numDiffuseRefl  << endl;
  DEBUG << "numSpecularRefl: " << numSpecularRefl << endl;
  DEBUG << "numDiffuseTrans:  " << numDiffuseTrans  << endl;
  DEBUG << "numSpecularTrans: " << numSpecularTrans << endl;
  */

  int numRays = 0;

  // For the rays that goes to the same halfspace where the input ray
  // come from -> reflection.
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  // ideally reflected ray
  const CVector3D &R = hitInfo.GetRay()->ReflectRay(N);
  assert(DotProd(R, N) >= 0.0);
  
  if (numDiffuseRefl) {
    bool transmissionFlag = false;
    numRays += SpawnStochasticDiffuseRays(kd_local_refl, transmissionFlag,
					  carray, hitInfo,
					  samplingInfo, numDiffuseRefl);
  }
  if (numSpecularRefl) {
    CValues shininess;
    hitInfo.GetTextureValue(shininess, *_specnumReflection, CMappingBase::EE_KSI);
    bool transmissionFlag = false;
    numRays += SpawnStochasticSpecularRays(ks_local_refl, shininess, R,
					   transmissionFlag, carray,
					   hitInfo,
					   samplingInfo,numSpecularRefl);
  }

  // now shoot rays for transmission.
  if (numDiffuseTrans) {
    bool transmissionFlag = true;
    numRays += SpawnStochasticDiffuseRays(kd_local_trans, transmissionFlag,
					  carray, hitInfo,
					  samplingInfo, numDiffuseRefl);
  }

  // Specular refraction - the most difficult case
  // recompute the R - the direction of ideally refracted ray
  if (numSpecularTrans) {
    float dot   = hitInfo.GetDotProd();
    float gamma, sqrterm;

    CValues iorReal(CValues::EE_Weigth1);
    if (_iorReal)
      // get value of index of refraction  coefficient for intersection point
      hitInfo.GetTextureValue(iorReal, *_iorReal, CMappingBase::EE_KRI);
    else {
      iorReal.Set(1.0); // undefined - no other choice
    }
    assert(iorReal.GetType() == CValues::EE_Weigth1);

    CValues iorImag(CValues::EE_Weigth1);
    if (_iorImag)
      // get value of index of refraction  coefficient for intersection point
      hitInfo.GetTextureValue(iorImag, *_iorImag, CMappingBase::EE_KRI);
    else {
      iorImag.Set(1.0); // undefined - no other choice
    }
    assert(iorImag.GetType() == CValues::EE_Weigth1);
  
    float ior = sqrt(sqr(iorReal.MeanValue()) + sqr(iorImag.MeanValue()));
    if (dot < 0.0) {
      // from outside into the inside of object
      gamma = 1.0 / ior;
    }
    else {
      // from the inside to outside of object
      gamma = ior;
      dot = -dot;
      N = -N;
    }
    sqrterm = 1.0 - gamma * gamma * (1.0 - dot * dot);
  
    sqrterm = dot * gamma + sqrt(sqrterm);

    // test if we did not get the internal reflection
    if (sqrterm >= 0.0) {
      // so ideally reflected ray exist
      CVector3D dirRefracted = -sqrterm * N + hitInfo.GetRay()->GetDir() *
	gamma;  
  
      CValues shininess;
      hitInfo.GetTextureValue(shininess, *_specnumTransmission,
			      CMappingBase::EE_KSIT);
      bool transmissionFlag = true;
      numRays += SpawnStochasticSpecularRays(ks_local_trans, shininess,
					     dirRefracted, transmissionFlag,
					     carray, hitInfo,
					     samplingInfo,
					     numSpecularTrans);
    }
  } // if (numSpecularTrans)
  
  //  DEBUG << "returning " << numRays << endl;
  return numRays;
}
  

// Given input and output direction, evaluate BSDF
bool
CSurfaceMGF::EvaluateBSDF(CValues &bsdf,
			  const CHitPointInfo &hitInfo,
			  const CVector3D &dirOut,
			  EDistributionType distType) const
{

  CValues res(bsdf.GetType()); // has zero value as default
  bool retval = false;

  // REFLECTION - does it occur now ?

  // the cosine between incoming ray and normal
  float dotProd = hitInfo.GetDotProd();
  // the cosine between outgoing ray and normal
  float cosN_L = DotProd(hitInfo.GetNormal(), dirOut);

  // Taking into account two-sideness of our surface
  // CVector3D normal = hitInfo.GetNormal();
  
  if ( ((dotProd < 0.0) && (cosN_L > 0.0)) ||
       ((dotProd > 0.0) && (cosN_L < 0.0)) ) {

    // reflection occurs for sure
    if (cosN_L > 0.0) {
      // evaluating diffuse term for reflection
      if ((distType & EE_REFLECTED_DIFFUSE) && _diffuseReflection) {
	// We don't check the angle between normal and outgoing direction.
	// Maybe we should (not to illuminate backfaces), but in that case
	// the dot product must be evaluated here and in the renderer as well.
	CValues diffuseTerm(res.GetType());
	hitInfo.GetTextureValue(diffuseTerm, *_diffuseReflection,
				CMappingBase::EE_KD);
	res += diffuseTerm * M_1_PI;
	retval = true;
      }

      // evaluating specular term for reflection
      CValues specularTerm(res.GetType());

      // if surface has specular component
      if ((distType & EE_REFLECTED_GLOSSY) && _specularReflection) {
	// now get the index of imperfection for specular reflection
	CValues specularIndex(CValues::EE_Weigth3);
	if (_specnumReflection)
	  hitInfo.GetTextureValue(specularIndex, *_specnumReflection,
				  CMappingBase::EE_KSI);
	else
	  // supposing that if not specified, surface is diffuse
	  specularIndex.Set(1.0);

	// now perform summing of diffuse and specular component
	float ratio;

#ifdef PHONG_CLASSIC_NORMALS
	// evaluation of the specular component according to Phong's
	// original paper
	ratio = DotProd(hitInfo.GetReflectedDir(), dirOut);
#else
	// evaluation of the specular component according to Blinn's article
	CVector3D hVector = 0.5 * (normal + dirOut);
	ratio = DotProd(hVector, normal);
#endif

	// only if the specular component makes sense
	if (ratio > 0.0) {
#ifdef MAX_PHONG
	  // max-Phong model
	  float cos_in  = fabs(dotProd);
	  float cos_out = fabs(cosN_L);
	  double max_cos = Max(cos_in, cos_out);
#endif
	  // get the value of specular component
	  hitInfo.GetTextureValue(specularTerm, *_specularReflection,
				  CMappingBase::EE_KS);

	  if (specularTerm.GetCount() == specularIndex.GetCount()) {
	    for (int i = specularTerm.GetCount(); i--;)
	      specularTerm.Mul(i, pow(ratio, specularIndex.Get(i))
#ifdef MAX_PHONG
			       / max_cos
#endif
			       );
	  }
	  else {
	    assert(specularIndex.GetCount() == 1);
	    // fr = ks * cos^n / max_cos
	    specularTerm *=
	      (
	       pow(ratio, specularIndex.Get(0))
#ifdef MAX_PHONG
	       / max_cos
#endif
	       );
	  }
	  // fr is composed of specular and diffuse term
	  res += specularTerm;
	  retval = true;
	} // ratio
      } // specular term for reflection
    } // cosN_L > 0
  }
  else { // TRANSMISSION - it must occur now
    if (cosN_L < 0.0) {
      // the transmission occurs
    
      // cosine of outgoing vector
      cosN_L = fabs(cosN_L);
    
      // evaluating diffuse term for transmission
      if ((distType & EE_REFLECTED_DIFFUSE) && _diffuseTransmission) {
	// We don't check the angle between normal and outgoing direction.
	// Maybe we should (not to illuminate backfaces), but in that case
	// the dot product must be evaluated here and in the renderer as well.
	CValues diffuseTerm(res.GetType());
	hitInfo.GetTextureValue(diffuseTerm, *_diffuseTransmission,
				CMappingBase::EE_KDT);
	res += diffuseTerm * M_1_PI;
	retval = true;
      }

      // evaluating specular term for transmission
      CValues specularTerm(res.GetType());

      // if surface has specular component
      if ((distType & EE_REFLECTED_GLOSSY) && _specularTransmission) {
	// now get the index of imperfection for specular reflection
	CValues specularIndex(CValues::EE_Weigth3);
	if (_specnumTransmission)
	  hitInfo.GetTextureValue(specularIndex, *_specnumTransmission,
				  CMappingBase::EE_KSIT);
	else
	  // supposing that if not specified, surface is diffuse
	  specularIndex.Set(1.0);

	// now perform summing of diffuse and specular component
	CVector3D idealRefractedDir;

	CVector3D N = hitInfo.GetNormal();
	float dot   = hitInfo.GetDotProd();
	float gamma, sqrterm;

	CValues iorReal(CValues::EE_Weigth1);
	if (_iorReal)
	  // get value of index of refraction  coefficient for intersection
	  // point
	  hitInfo.GetTextureValue(iorReal, *_iorReal, CMappingBase::EE_KRI);
	else
	  iorReal.Set(1.0);
	CValues iorImag(CValues::EE_Weigth1);
	if (_iorImag)
	  // get value of index of refraction  coefficient for intersection
	  // point
	  hitInfo.GetTextureValue(iorImag, *_iorImag, CMappingBase::EE_KRI);
	else
	  iorImag.Set(1.0);

	float ior = sqrt(sqr(iorReal.Get(0)) + sqr(iorImag.Get(0)));

	if (dot < 0.0) {
	  // from outside into the inside of object
	  gamma = 1.0 / ior;
	}
	else {
	  // from the inside to outside of object
	  gamma = ior;
	  dot = -dot;
	  N = -N;
	}
	sqrterm = 1.0 - gamma * gamma * (1.0 - dot * dot);
  
	// Check for total internal reflection, do nothing if it applies.
	if (sqrterm < 0.0)  {
	   // total internal reflection, but value computed correctly!
	  goto finish_bsdf;
	}
	else {
	  sqrterm = dot * gamma + sqrt(sqrterm);
	  idealRefractedDir = -sqrterm * N + hitInfo.GetRay()->GetDir()*gamma;
	}	

	// the dot product between outgoing direction and some vector
	float ratio;
#ifdef PHONG_CLASSIC_NORMALS
	// evaluation of the specular component according to Phong's
	// original paper
	ratio = DotProd(idealRefractedDir, dirOut);	
#else
	// evaluation of the specular component according to Blinn's article
	CVector3D hVector = 0.5 * (normal + dirOut);
	ratio = DotProd(hVector, normal);
#endif

	if (ratio > 0.0) {
#ifdef MAX_PHONG
	  // max-Phong model
	  float cos_in  = fabs(dotProd);
	  float cos_out = fabs(cosN_L);
	  double max_cos = Max(cos_in, cos_out);
#endif
	  // get the value of specular component
	  hitInfo.GetTextureValue(specularTerm, *_specularTransmission,
				  CMappingBase::EE_KST);

	  if (specularTerm.GetCount() == specularIndex.GetCount()) {
	    for (int i = specularTerm.GetCount(); i--;)
	      specularTerm.Mul(i, pow(ratio, specularIndex.Get(i))
#ifdef MAX_PHONG
			       / max_cos
#endif
			       );
	  }
	  else {
	    assert(specularIndex.GetCount() == 1);
	    // fr = ks * cos^n / max_cos
	    specularTerm *=
	      (
	       pow(ratio, specularIndex.Get(0))
#ifdef MAX_PHONG
	       / max_cos
#endif
	       );
	  }
	  // fr is composed of specular and diffuse term
	  res += specularTerm;
	  retval = true;
	} // ratio > 0.0
      } // specular transmission possible
    } // cosine has appropriate 
  } // transmission  


finish_bsdf:
  if (retval) {  
    bsdf = res * cosN_L;
    // DEBUG << "bsdf: " << bsdf << endl;
    return true; // OK - value computed
  }

  // for sure set the value to zero
  bsdf.Set(0.0);  
  return false; // problem/error
}

const char* 
CSurfaceMGF::GetHelp() const
{
	return BrdfStrings::K_MGF_Desc;
}

int 
CSurfaceMGF::GetNumParameters() const
{
	return 8; /* diffuseReflection, specularReflection,specnumReflection, diffuseTransmition,
				 specularTransmission, specnumTransmission, iorReal, iorImag*/
}

CValues
CSurfaceMGF::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	switch (parameterIndex)
	{
	case 0 : //diffuse refl
		_diffuseReflection->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 1 : //specular refl
		_specularReflection->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 2 : //specnum refl
		_specnumReflection->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 3 : //diffuse refr
		_diffuseTransmission->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 4 : //specular refr
		_specularTransmission->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 5 : //specnum refr
		_specnumTransmission->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 6 : //ior real
		_iorReal->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 7 : //ior img
		_iorImag->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	default :
		break;
	}
	return c;
}

CValues
CSurfaceMGF::GetParameterConstraints(int parameterIndex) const
{
	CValues c(3);
	switch (parameterIndex)
	{
	case 0 : //diffuse refl
	case 1 : //specular refl
	case 3 : //diffuse refr
	case 4 : //specular refr
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 2 : //specnum refl
	case 5 : //specnum refr
		c.Set(0,0); // min
		c.Set(1,MAX_SHININESS); // max
		c.Set(2,1); // default
		break;
	case 6 : //ior real
	case 7 : //ior img
		c.Set(0,1); // min
		c.Set(1,MAX_IOR); // max
		c.Set(2,1); // default
		break;
	default :
		break;
	}
	return c;
}

const char* 
CSurfaceMGF::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseReflDesc;
	case 1 :
		return BrdfStrings::KSpecularReflDesc;
	case 2 :
		return BrdfStrings::KShininessReflDesc;
	case 3 :
		return BrdfStrings::KDiffuseRefrDesc;
	case 4 :
		return BrdfStrings::KSpecularRefrDesc;
	case 5 :
		return BrdfStrings::KShininessRefrDesc;
	case 6 :
		return BrdfStrings::KIORDesc;
	case 7 :
		return BrdfStrings::KExtinctionCoefDesc;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}

CSurface::EParameterUIType
CSurfaceMGF::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 : //diffuse refl
		return EE_ColorRGB;
	case 1 : //specular refl
		return EE_ColorRGB;
	case 2 : //specnum refl
		return EE_Float;
	case 3 : //diffuse refr
		return EE_ColorRGB;
	case 4 : //specular refr
		return EE_ColorRGB;
	case 5 : //specnum refr
		return EE_Float;
	case 6 : //ior real
		return EE_Float;
	case 7 : //ior img
		return EE_Float;
	default :
		return EE_Unknown;
	}
}

const char*
CSurfaceMGF::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseReflTitle;
	case 1 :
		return BrdfStrings::KSpecularReflTitle;
	case 2 :
		return BrdfStrings::KShininessReflTitle;
	case 3 :
		return BrdfStrings::KDiffuseRefrTitle;
	case 4 :
		return BrdfStrings::KSpecularRefrTitle;
	case 5 :
		return BrdfStrings::KShininessRefrTitle;
	case 6 :
		return BrdfStrings::KIORTitle;
	case 7 :
		return BrdfStrings::KExtinctionCoefTitle;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}

void
CSurfaceMGF::SetParameter(int parameterIndex, const CValues& value)
{
	CConstantTexture* newTex = new CConstantTexture(value);
	switch (parameterIndex)
	{
	case 0 :
		SetDiffuseReflection(newTex);
		return;
	case 1 :
		SetSpecularReflection(newTex);
		return;
	case 2 :
		SetSpecnumReflection(newTex);
		return;
	case 3 :
		SetDiffuseTransmission(newTex);
		return;
	case 4 :
		SetSpecularTransmission(newTex);
		return;
	case 5 :
		SetSpecnumTransmission(newTex);
		return;
	case 6 :
		SetIORReal(newTex);
		return;
	case 7 :
		SetIORImag(newTex);
		return;
	default :
		break;
	}
}

__END_BRDFLIB_SOURCE
