// ===================================================================
// $Id$
//
// kurt_kalos.cpp
//    Implementation of Kurt-Kalos BRDF.
//
// Class: CSurfaceKurtKalosAnisotropic, CSurfaceKurtKalosIsotropic
//
// Initial coding by Jiri Nekola within diploma thesis project, 
// November 2011
// ===================================================================

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

__BEGIN_BRDFLIB_SOURCE


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

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

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

		if (_slope) 
		{
			indent(app, ind+1);
			app << "  slope    = \t";
			_slope->Describe(app, ind); app << "\n";
		}
	}

	app << flush;

	return;
}

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

	if (_diffuse)
	{
		ret->SetDiffuse(_diffuse);
	}
	if (_specular)
	{
		ret->_specular = _specular;

		if(_slope)
		{
			ret->_slope = _slope;
		}
	}

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

/*
  Albedo is the integral of BRDF multiplied by cos(Eye,N) over the hemisphere.
  Integral of Kurt-Kalos BRDF cannot be evaluated analitically. However, it
  can be evaluated using Pade approximation. Unfortunately, it requres a lot
  of computations. So, for a now we leave this albedo function in this way.

  Possibly later on, we'll provide this integral evaluation in the way
  similar to the Arvo's procedures (Arvo'1995 "Applications of Irradiance
  Tensors to the Simulation of Non-Lambertian Phenomena"). */
//========================================================================
bool
CSurfaceKurtKalosAnisotropic::GetAlbedo(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
  // Aldebo is a function of the angle bewtween the incoming direction
  // and normal, but in most cases only the maximum of this function
  // computed below is intresting
  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 & _microfacetReflectionType()) && _specular) && _slope) {
    hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
    c += d;
    retval = true;
  }

//#define _DEBUG_PHYSICAL_PLAUSIBILITY
#ifdef _DEBUG_PHYSICAL_PLAUSIBILITY
  // test if the albedo is smaller than 1.0 -- physically plausible
  if (c.Max() > 1.0) {
    if (_diffuse) {
      hitInfo.GetTextureValue(d, *_diffuse, CMappingBase::EE_KD);
      DEBUG << "diffuse = " << d << "\n";
    }
    if (_specular) {
      hitInfo.GetTextureValue(d, *_specular, CMappingBase::EE_KS);
      DEBUG << "specular = " << d << "\n";
    }
    cout << "Current parameter settings for Cook-Torrance model make it\n";
    cout << "physically unplausible\n";
    exit(0);
  }
#endif // _DEBUG_PHYSICAL_PLAUSIBILITY

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

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

  return false;
}

/*=============================================================================
			 kd * (1 - F(w0.h))             ks * F(w0.h) * D(h)
   BRDF  =  -------------------- + ------------------------------------
					Pi			   4 * (w0.h) * ((wi.n) * (w0.n))^alpha

   where 
	- (x.y) is dot product of x and y
	- kd is diffuse albedo
	- ks is specular albedo
	- w0 is outgoing(viewer) direction
	- wi is incoming(light) direction
	- n is normal
	- h is halfvector between wo and wi (normalized (w0+wi))
	- alpha is shadowing/masking normalization factor
	- D(h) is Beckmann distribution :
		- for anisotropic case :
			                    1
			D(h) = ----------------------------- * q(h)
			       Pi * mx * my * cos^4 (thetaH)

		                                 cos^2(phiH)    sin^2(phiH)
			q(h) = e^(-tan^2(thetaH) * (------------ + ------------))
										     mx^2		   my^2
		
		where 
		 - phiH, thetaH are polar coordinates of halfvector h
		 - mx, my are standard deviations (RMS) of the surface slope in 
		   principal directions of anisotropy

  (see reference [1], equation (13))
=============================================================================*/
// Direct lighting computation
//=============================================================================
bool
CSurfaceKurtKalosAnisotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirLit,
			 EDistributionType   distType) const
{
	bool retval 	    = false;
	float d           = Magnitude(dirLit);
	CVector3D N       = hitInfo.GetNormal();
	CVector3D dirEye  = (hitInfo.GetRay())->GetDir();
	float lightCos    = DotProd(dirLit, N);
	float incEyeCos   = -(hitInfo.GetDotProd());

	result.Set(0);

	// Provide 'zero' result if Light source or viewer is under the surface
	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		CVector3D anisAxis = hitInfo.GetAnisAxis();

		if(lightCos < WTINY) lightCos       = WTINY;
		if(incEyeCos < WTINY) incEyeCos  = WTINY;

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
			Normalize(dirLit);

		CValues F(CValues::EE_Weigth3);

		// Specular part
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) 
		{
			CTexture *slopeX = _slope->GetSlopeX();
			CTexture *slopeY = _slope->GetSlopeY();
			if(slopeX && slopeY)
			{
				CValues rho_s, slopeXc, slopeYc;
				hitInfo.GetTextureValue(rho_s,  *_specular, CMappingBase::EE_KS);

				hitInfo.GetTextureValue(slopeXc, *slopeX,   CMappingBase::EE_KSI);
				hitInfo.GetTextureValue(slopeYc, *slopeY,   CMappingBase::EE_KSI);
				float ax  = slopeXc.Get(0), ay = slopeYc.Get(0);
				if(ax < WTINY) ax = WTINY;
				if(ay < WTINY) ay = WTINY;
				float ax2 = ax * ax;
				float ay2 = ay * ay;
				float ax_ay_pi = M_PI * ax * ay;

				if(rho_s.Luminance()/ax_ay_pi > 1e-06)
				{
					CVector3D  H  =  Normalize(dirLit - (hitInfo.GetRay()->GetDir()));
					float cos_HN  =  DotProd(H, N);
					assert(cos_HN >= 0);

					if(cos_HN < WTINY)  cos_HN = WTINY;

					float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
					if(tan2_HN < 0) tan2_HN = 0;

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

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

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

					// q(h)
					float ex  = exp( -tan2_HN * (cos2_phi/ax2 + sin2_phi/ay2));

					// D(h)
					float D = ex / (ax_ay_pi * pow(cos_HN,4.f));

					float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),H);

					// compute fresnel equation for (w0.h)
					hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);
					
					CValues aux;
					hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
					float a = aux.Get(0);

					// 1 / (4 * (w0.h) * ((wi.n) * (w0.n))^alpha)
					float mult = 1.f / (4 * cosHV * pow(lightCos * incEyeCos,a));

			
					// Add all together
					result += rho_s * F * D * mult;
					retval = true;

	
				} // not black body
			} // both slopes exist
		} // specular exist

		// Diffuse part
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
			CValues rho_d;
			hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
			result += rho_d * M_1_PI * (1.0 - F.Luminance()); // coupled by eq. (13) in ref [1]
			retval = true;
		}
	}
	else // under the surface
	{
		retval = false;
	}
	result *= lightCos;
	return retval;
}

// Diffuse part for stochastic sampling
//========================================================================
int
CSurfaceKurtKalosAnisotropic::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  const CValues       &ks_local,
  CVector2D           &slopeXY,
  CVector2D           &slopeXY_sq,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
	CVector3D N      = hitInfo.GetNormal();
	CVector3D anisAxis = hitInfo.GetAnisAxis();
	CVector3D dirEye = (hitInfo.GetRay())->GetDir();
	float incEyeCos  = hitInfo.GetDotProd();

	// If light source/viewer is under the surface
	if (incEyeCos > 0.0)
		N = -N;

	incEyeCos = fabs(incEyeCos);

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

	int num = 0;			       // number of rays created
	float ax_ay_pi = M_PI * slopeXY.xx * slopeXY.yy;

	// For translation of half-vector into the local coordinate system
	CVector3D lx(1,0,0), ly(0,1,0), lz(0,0,1), hl(1,0,0);
	if(ks_local.Luminance()/ax_ay_pi > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	// Compute directions for maxRays outgoing rays
	//----------------------------------------------------------------------
	for(int i = maxRays; i--; num++) 
	{
		// For Halton sequence random numbers generation
		//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 computed sample dirN 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 weight 
		//         = 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
		// visualisation
		//------------------------------------------------------------------------
		float lightCos = DotProd(dirRes, N);
		assert(lightCos >= -WTINY);

		if(incEyeCos < WTINY) incEyeCos = WTINY;
		if(lightCos  < WTINY) lightCos  = WTINY;

		CVector3D  H  =  Normalize(dirRes - (hitInfo.GetRay()->GetDir()));
		float cos_HN  =  DotProd(H, N);
		if(cos_HN < WTINY)  cos_HN = WTINY;

		float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
		if(tan2_HN < 0) tan2_HN = 0;

		// Now let us get phi polar angle for half-vector
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		hl[0] = DotProd(H, lx);
		hl[1] = DotProd(H, ly);
		hl[2] = DotProd(H, lz);
		assert(Magnitude(hl) > 1.0-1e-06 && Magnitude(hl) < 1.0+1e-06);

		// make projection onto XZ plane
		float hx2 = hl[0]*hl[0], hy2 = hl[1]*hl[1];
		float len_hprj2 = hx2 + hy2;
		//float len_hprj  = sqrt(len_hprj2);
		float cos2_phi   = hx2 / len_hprj2;
		float sin2_phi   = hy2 / len_hprj2;
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		// q(h)
		float ex  = exp( -tan2_HN * (cos2_phi/slopeXY_sq.xx + 
			                         sin2_phi/slopeXY_sq.yy));

		// D(h)
		float D = ex / (ax_ay_pi * pow(cos_HN,4.f));

		float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),H);

		// compute fresnel equation for (w0.h)
		CValues F(CValues::EE_Weigth3);
		hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);
		
		CValues aux;
		hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
		float a = aux.Get(0);

		// 1 / (4 * (w0.h) * ((wi.n) * (w0.n))^alpha)
		float mult = 1.f / (4 * cosHV * pow(lightCos * incEyeCos,a));


		// Add all together
		float k = F.Luminance() * D * mult;
		k /= contrib->prob;

		// add specular part to computed weight
		contrib->weigth += ks_local * k;

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

	} // for i

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

/* Specular part for stochastic sampling which is implemented in a similar
   to importance sampling scheme for Ward model way. (see ref [1])
   
							1
   pdf = ------------------------------------------ * q(h)
          4 * Pi * mx * my * cos^3(thetaH) * (w0.h)

			   ks * F(w0.h) * (wi.n)^(1-alpha)   
   weight w = ---------------------------------
                 cos(thetaH) * (w0.n)^alpha
   
   for parameters description, see EvaluateBSDF method commentary

*/
//========================================================================
int
CSurfaceKurtKalosAnisotropic::SpawnStochasticSpecularRays(
  const CValues       &ks_local,
  CVector2D           &slopeXY,
  CVector2D           &slopeXY_sq,
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
    // the number of samples generated at the end of this function
	int num = 0;

	float incEyeCos    = hitInfo.GetDotProd();
	CVector3D        N = hitInfo.GetNormal();
	CVector3D anisAxis = hitInfo.GetAnisAxis();

	// now holds DotProd(N, incomingRay) < 0.0
	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	CVector3D dirInit  = (hitInfo.GetRay())->GetDir();

	if (incEyeCos > 0.0)
	{
		N = -N;
	}

	// Make new basis with R as Z axis
	CVector3D V = CrossProd(N, anisAxis);    // V is perpendicular to R and V - Y-axis

	// diffuseness of the bounce
	//!!!!!!!!!!!!!!!!!!!!!!    it is incorrect     !!!!!!!!!!!!!!!!!!!!!!
	//!!!!!!!!!!!!!!!!!!!!!! it is just temporarily !!!!!!!!!!!!!!!!!!!!!!
	float diffuseness = (slopeXY.xx + slopeXY.yy) / 2;//_nToDiffuseness(shininess.Sum() / shininess.GetCount());
	
	float anisPart = 1;

	// for outgoing direction
	//float multi = M_PI * a2;
	//cout << "3ps\n";

	float ax_ay_pi_4 = 4.f * M_PI * slopeXY.xx * slopeXY.yy;

	for(int i = maxRays; i--; num++) 
	{
		// the ray is sampled in this direction
		CContributionItem *contrib = carray.AllocateNext();
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);

		float lightCos = -1, cos_HI = -1;
		CVector3D dirRes(0,0,1);
		if((slopeXY.xx >= 1e-03) || (slopeXY.yy >= 1e-03))
		{
			CVector3D hvector(0,0,1);
			float cos_HN = 0, cos2_HN = 0, sin_HN = 1, tan2_HN = 0;
			float sin_phi, cos_phi;
			float phi_HN;
			bool rejected = false;

			while(lightCos < 0) 
			{
				CVector2D rn; // random numbers
				if( !rejected )
					rn.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
				else
					rn.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

				// according to (16),(17) in ref [1] :
				//									-log(1-s)
				// delta = atan(sqrt(-------------------------------------))
				//					  cos^2(phiH)/mx^2 + sin^2(phiH)/my^2
				//
				// where alfa and beta are the standard deviations of the surface slope and
				//
				// phiH = atan2(my*sin(2*pi*t),mx*cos(2*pi*t))
				//
				// where s and t are independent random numbers uniformly distributed in [0, 1),
				//		
				float pi_2_t	=  2.0 * M_PI * rn.yy;
				phi_HN			=  atan2(slopeXY.yy * sin(pi_2_t),slopeXY.xx * cos(pi_2_t));
				cos_phi			=  cos(phi_HN);
				sin_phi			=  sqrtf(1-cos_phi*cos_phi);
				anisPart		=  (cos_phi*cos_phi)/slopeXY_sq.xx + (sin_phi*sin_phi)/slopeXY_sq.yy;

				// we dont need compute exactly delta ( = HN), it is sufficient to compute tan^2(delta), from 
				// which is easy to compute sin(delta) and cos(delta) (less expensive than atan(sqrt()) and then
				// back to tan^2
				// 					  -log(1-s)							-log(1-s)
				// delta = atan(sqrt(-----------)) ---> tan^2(delta) = ----------- 
				//					  anisPart							anisPart
				tan2_HN = (-log(1.0-rn.xx))/anisPart;
				assert(tan2_HN >= 0);
				cos2_HN  = (1.0/(1.0 + tan2_HN));
				cos_HN   = sqrt(cos2_HN);
				sin_HN   = sqrt(1.0 - cos2_HN);

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

				CVector3D localDir(sin_HN * cos_phi, sin_HN * sin_phi, cos_HN);

				// Change the basis from the VUN coordinate system (local coordinate
				// system where normal R is Z-axis) to the World coordinate system
				hvector = localDir.x * anisAxis + localDir.y * V + localDir.z * N;

				hvector.Normalize();


				cos_HI = -DotProd(hvector,dirInit);

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

				if(cos_HI >= 0)
				{
					// dirRes = 2*cos(HI)*h - dirInit, but dirInit is negative in our case
					// (points toward the surface), so its - (-dirInit) = + dirInit
					dirRes = (2.0 * cos_HI) * hvector + dirInit; 

					assert(Magnitude(dirRes) >= 1.0-1e-04 && Magnitude(dirRes) <= 1.0+1e-04);
					lightCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

					if(lightCos < 0)
					{
						rejected = true;
					}
				}
				else //(lightCos < 0)
				{
					rejected = true;
				}
			}

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(lightCos >= 0);

			float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),hvector);
			float k = 1.f / (ax_ay_pi_4 * pow(cos_HN,3.f) * cosHV);

			// q(h)
			float q = exp( -tan2_HN * anisPart );

			contrib->prob = k * q;

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

			CValues aux;
			hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
			float a = aux.Get(0);

			CValues F(CValues::EE_Weigth3);
			// compute fresnel equation for (w0.h)
			hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);

			float w = (F.Luminance() * pow(lightCos,1-a)) / (cos_HN * pow(incEyeCos,a));

			contrib->weigth = ks_local * w;

			// set the reflected ray, direction is normalized here
			contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
			// type of interaction
			contrib->interaction = EE_REFLECTED_SPECULAR;
			// diffuseness of the bounce
			contrib->diffuseness = diffuseness;
		}
		else // both slopes lower than 1e-3 -> mirror reflection
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			lightCos = DotProd(dirRes, N);
			contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
			// type of interaction
			contrib->interaction = EE_REFLECTED_MIRROR;
			// diffuseness of the bounce
			contrib->diffuseness = 0.0;

			// exp(-tan(0)/slopeXY) = 1, cos(H,N) = 1 =>
			contrib->prob = 1.0 / (slopeXY.xx * slopeXY.yy * M_PI);

			// For mirror cosine of incident equal to the cosine of
			// reflected ray => sqrt(lightCos/incEyeCos) =>
			contrib->weigth = ks_local * 0.25;
		}

		if(_diffuse && (kd_local.Luminance() > 1e-05)) 
		{
		   //kd_local *= 1.0/(M_PI*contrib->prob);
		   contrib->weigth += kd_local * (lightCos/(M_PI*contrib->prob));
		}

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

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


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

	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);
	CVector2D slopeXY(0,0), slopeXY_sq(0,0);

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

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

	if (_specular && (distType & EE_REFLECTED_SPECULAR)) 
	{
		hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
		aks = ks_local.Luminance();

		if(!_slope)
		{
			WARNING_INFO("CSurfaceKurtKalosAnisotropic::SpawnStochasticRays() :\n\
Slope structure should be set to compute specular part of Ward model.");
		  return 0;
		}

		CTexture *slopeX = _slope->GetSlopeX();
		CTexture *slopeY = _slope->GetSlopeY();
		if((!slopeX) || (!slopeY))
		{
			WARNING_INFO("CSurfaceKurtKalosAnisotropic::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope', which consists\n\
of 'slopeX' and 'slopeY' terms.");
			return 0;
		}
		CValues slopeXc(0,0,0), slopeYc(0,0,0);
		hitInfo.GetTextureValue(slopeXc, *slopeX, CMappingBase::EE_KSI);
		hitInfo.GetTextureValue(slopeYc, *slopeY, CMappingBase::EE_KSI);
		float slopeXf = slopeXc.Get(0);
		float slopeYf = slopeYc.Get(0);
		(slopeXf < WTINY) ? slopeXf = WTINY : slopeXf = slopeXf;
		(slopeYf < WTINY) ? slopeYf = WTINY : slopeYf = slopeYf;

		slopeXY.xx    = slopeXf;
		slopeXY_sq.xx = slopeXf * slopeXf;
		slopeXY.yy    = slopeYf;
		slopeXY_sq.yy = slopeYf * slopeYf;
		tmp_ratio     = aks/(slopeXY.xx * slopeXY.yy * 4.0 * M_PI);
	} // specular exists

	float albedo_whole = akd + aks;
  

	if(albedo_whole < 1e-05) 
	{
		WARNING_INFO("CSurfaceKurtKalosAnisotropic::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
		return 0;
	}
	if( albedo_whole >= 1.0 ) 
	{
		WARNING_INFO("CSurfaceKurtKalosAnisotropic::SpawnStochasticRays():\nAlbedo > 1.0\n\
To keep the energy conservation requirement the total albedo\n\
(albedo_specularity + albedo_diffuse) 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) || (tmp_ratio <= 1e-06)) 
		{
			// 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;

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

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

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

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

	samplingInfo.IncDimension(2);

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

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

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

const char* 
CSurfaceKurtKalosAnisotropic::GetHelp() const
{
	return BrdfStrings::KKurtKalosAnisoDesc;
}
//========================================================================

int 
CSurfaceKurtKalosAnisotropic::GetNumParameters() const
{
	// diffuse + fresnel term + alpha normalizator + slopeX + slopeY
	return 5; 
}
//========================================================================

CValues
CSurfaceKurtKalosAnisotropic::GetParameter(int parameterIndex) const
{
	CValues c(4);

	switch (parameterIndex)
	{
	case 0 :
		_diffuse->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 1 :
		{
			CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_specular);
			c.SetRepresentation(CValues::EE_WeigthN);
			CValues aux(CValues::EE_Weigth3);
			tex->GetColor(aux);
			c.Set(0,aux.Get(0));
			c.Set(1,aux.Get(1));
			c.Set(2,aux.Get(2));
			tex->GetIndexR(aux);
			c.Set(3,aux.Get(0));
			tex->GetExtinction(aux);
			c.Set(4,aux.Get(0));
			bool schlick;
			tex->GetApproxSchlick(schlick);
			c.Set(5,schlick);
			break;
		}
	case 2:
		_alpha->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 3:
		_slope->GetSlopeX()->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	case 4:
		_slope->GetSlopeY()->GetValue(c,CHitPointInfo::constMapInfo);
		break;
	}

	return c;
}

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

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

	switch (parameterIndex)
	{
	case 0 : // diffuse
		_setBasicConstraints(EE_ColorRGB,c);
		break;
	case 1 : // specular
		_setBasicConstraints(EE_Fresnel,c);
		break;
	case 2 : // alpha
		c = CValues(3);
		c.Set(0,0); // min 
		c.Set(1,10); // max 
		c.Set(2,0.2f); // default 
		break;
	case 3: // slope X
	case 4: // slope Y
		c = CValues(3);
		c.Set(0,0.001); // min 
		c.Set(1,M_PI/2); // max 
		c.Set(2,0.5f); // default 
		break;
	}

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

const char* 
CSurfaceKurtKalosAnisotropic::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseDesc;
	case 1 :
		return BrdfStrings::KFresnelDesc;
	case 2 :
		return BrdfStrings::KAlphaDesc;
	case 3 :
		return BrdfStrings::KSlopeXDesc;
	case 4 :
		return BrdfStrings::KSlopeYDesc;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CSurfaceKurtKalosAnisotropic::GetParameterSyntax(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 : // diffuse
		return EE_ColorRGB;
	case 1 : // specular
		return EE_Fresnel;
	case 2 : // alpha
	case 3 : // slope X
	case 4 : // slope Y
		return EE_Float;
	default :
		return EE_Unknown;
	}
}
//========================================================================

const char*
CSurfaceKurtKalosAnisotropic::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
	case 0 :
		return BrdfStrings::KDiffuseTitle;
	case 1 :
		return BrdfStrings::KFresnelTitle;
	case 2 :
		return BrdfStrings::KAlphaTitle;
	case 3 :
		return BrdfStrings::KSlopeXTitle;
	case 4 :
		return BrdfStrings::KSlopeYTitle;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CSurfaceKurtKalosAnisotropic::SetParameter(int parameterIndex, const CValues& value)
{
	switch (parameterIndex)
	{
	case 0 :{
				CConstantTexture* tex = new CConstantTexture(value);
				SetDiffuse(tex);
				break;
			}
	case 1 :{
				if (!_specular)
				{
					_specular = new CFresnelTexture(CValues(CValues::EE_Weigth3));
				}
				CFresnelTexture* tex = dynamic_cast<CFresnelTexture*>(_specular);
				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;
			}
	case 2 :{
				CConstantTexture* tex = new CConstantTexture(value);
				SetAlpha(tex);
				break;
			}
	case 3 :
	case 4 :{
				float w = 0.5f;
				CValues slope_value(value[0],value[0],value[0]);
				CConstantTexture* newTex = new CConstantTexture(slope_value);
				if (!_slope)
				{
					_slope = new CSlopeMicrofacets();
				}
				_slope->Set(newTex,parameterIndex-3);
				_slope->SetWeight(w,parameterIndex-3);
				break;
			}
				
	}
}










/**********************************************************************
 *
 *             CSurfaceKurtKalosIsotropic methods
 *
 **********************************************************************/

/*=============================================================================
			 kd * (1 - F(w0.h))             ks * F(w0.h) * D(h)
   BRDF  =  -------------------- + ------------------------------------
					Pi			   4 * (w0.h) * ((wi.n) * (w0.n))^alpha

   where 
	- (x.y) is dot product of x and y
	- kd is diffuse albedo
	- ks is specular albedo
	- w0 is outgoing(viewer) direction
	- wi is incoming(light) direction
	- n is normal
	- h is halfvector between wo and wi (normalized (w0+wi))
	- alpha is shadowing/masking normalization factor
	- D(h) is Beckmann distribution :
		- for anisotropic case :
			                    1
			D(h) = ---------------------------- * q(h)
			        Pi * mx^2 * cos^4 (thetaH)

			q(h) = e^(-tan^2(thetaH) / mx^2 )						 
		
		where 
		 - phiH, thetaH are polar coordinates of halfvector h
		 - mx is standard deviation (RMS) of the surface slope

  (see reference [1], equation (13))
=============================================================================*/
// Direct lighting computation
//=============================================================================
bool
CSurfaceKurtKalosIsotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirLit,
			 EDistributionType   distType) const
{
	bool retval 	    = false;
	float d           = Magnitude(dirLit);
	CVector3D N       = hitInfo.GetNormal();
	CVector3D dirEye  = (hitInfo.GetRay())->GetDir();
	float lightCos    = DotProd(dirLit, N);
	float incEyeCos   = -(hitInfo.GetDotProd());

	result.Set(0);

	// Provide 'zero' result if Light source or viewer is under the surface
	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		CVector3D anisAxis = hitInfo.GetAnisAxis();

		if(lightCos < WTINY) lightCos       = WTINY;
		if(incEyeCos < WTINY) incEyeCos  = WTINY;

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
			Normalize(dirLit);

		CValues F(CValues::EE_Weigth3);

		// Specular part
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) 
		{
			CTexture *slopeT = _slope->GetSlopeX();
			if(slopeT)
			{
				CValues rho_s, slopeVal;
				hitInfo.GetTextureValue(rho_s,  *_specular, CMappingBase::EE_KS);

				hitInfo.GetTextureValue(slopeVal, *slopeT,   CMappingBase::EE_KSI);
				float ax  = slopeVal.Get(0);
				if(ax < WTINY) ax = WTINY;
				float ax2 = ax * ax;
				float ax_2_pi = M_PI * ax2;

				if(rho_s.Luminance()/ax_2_pi > 1e-06)
				{
					CVector3D  H  =  Normalize(dirLit - (hitInfo.GetRay()->GetDir()));
					float cos_HN  =  DotProd(H, N);
					assert(cos_HN >= 0);

					if(cos_HN < WTINY)  cos_HN = WTINY;

					float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
					if(tan2_HN < 0) tan2_HN = 0;

					// q(h)
					float ex  = exp( -tan2_HN / ax2);

					// D(h)
					float D = ex / (ax_2_pi * pow(cos_HN,4.f));

					float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),H);

					// compute fresnel equation for (w0.h)
					hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);
					
					CValues aux;
					hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
					float a = aux.Get(0);

					// 1 / (4 * (w0.h) * ((wi.n) * (w0.n))^alpha)
					float mult = 1.f / (4 * cosHV * pow(lightCos * incEyeCos,a));

			
					// Add all together
					result += rho_s * F * D * mult;
					retval = true;

	
				} // not black body
			} // both slopes exist
		} // specular exist

		// Diffuse part
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) {
			CValues rho_d;
			hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
			result += rho_d * M_1_PI * (1.0 - F.Luminance()); // coupled by eq. (13) in ref [1]
			retval = true;
		}
	}
	else // under the surface
	{
		retval = false;
	}
	result *= lightCos;
	return retval;
}

// Diffuse part for stochastic sampling
//========================================================================
int
CSurfaceKurtKalosIsotropic::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  const CValues       &ks_local,
  float				  &slope,
  float				  &slope_sq,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
	CVector3D N      = hitInfo.GetNormal();
	CVector3D anisAxis = hitInfo.GetAnisAxis();
	CVector3D dirEye = (hitInfo.GetRay())->GetDir();
	float incEyeCos  = hitInfo.GetDotProd();

	// If light source/viewer is under the surface
	if (incEyeCos > 0.0)
		N = -N;

	incEyeCos = fabs(incEyeCos);

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

	int num = 0;			       // number of rays created
	float ax_2_pi = M_PI * slope_sq;

	// For translation of half-vector into the local coordinate system
	CVector3D lx(1,0,0), ly(0,1,0), lz(0,0,1), hl(1,0,0);
	if(ks_local.Luminance()/ax_2_pi > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	// Compute directions for maxRays outgoing rays
	//----------------------------------------------------------------------
	for(int i = maxRays; i--; num++) 
	{
		// For Halton sequence random numbers generation
		//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 computed sample dirN 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 weight 
		//         = 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
		// visualisation
		//------------------------------------------------------------------------
		float lightCos = DotProd(dirRes, N);
		assert(lightCos >= -WTINY);

		if(incEyeCos < WTINY) incEyeCos = WTINY;
		if(lightCos  < WTINY) lightCos  = WTINY;

		CVector3D  H  =  Normalize(dirRes - (hitInfo.GetRay()->GetDir()));
		float cos_HN  =  DotProd(H, N);
		if(cos_HN < WTINY)  cos_HN = WTINY;

		float tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
		if(tan2_HN < 0) tan2_HN = 0;

		// q(h)
		float ex  = exp( -tan2_HN / slope_sq);

		// D(h)
		float D = ex / (ax_2_pi * pow(cos_HN,4.f));

		float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),H);

		// compute fresnel equation for (w0.h)
		CValues F(CValues::EE_Weigth3);
		hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);
		
		CValues aux;
		hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
		float a = aux.Get(0);

		// 1 / (4 * (w0.h) * ((wi.n) * (w0.n))^alpha)
		float mult = 1.f / (4 * cosHV * pow(lightCos * incEyeCos,a));


		// Add all together
		float k = F.Luminance() * D * mult;
		k /= contrib->prob;

		// add specular part to computed weight
		contrib->weigth += ks_local * k;

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

	} // for i

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

/* Specular part for stochastic sampling which is implemented in a similar
   to importance sampling scheme for Ward model way. (see ref [1])
   
							1
   pdf = ------------------------------------------ * q(h)
          4 * Pi * mx^2 * cos^3(thetaH) * (w0.h)

			   ks * F(w0.h) * (wi.n)^(1-alpha)   
   weight w = ---------------------------------
                 cos(thetaH) * (w0.n)^alpha
   
   for parameters description, see EvaluateBSDF method commentary

*/
//========================================================================
int
CSurfaceKurtKalosIsotropic::SpawnStochasticSpecularRays(
  const CValues       &ks_local,
  float				  &slope,
  float               &slope_sq,
  const CValues       &kd_local,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
    // the number of samples generated at the end of this function
	int num = 0;

	float incEyeCos    = hitInfo.GetDotProd();
	CVector3D        N = hitInfo.GetNormal();
	CVector3D anisAxis = hitInfo.GetAnisAxis();

	// now holds DotProd(N, incomingRay) < 0.0
	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	CVector3D dirInit  = (hitInfo.GetRay())->GetDir();

	if (incEyeCos > 0.0)
	{
		N = -N;
	}

	// Make new basis with R as Z axis
	CVector3D V = CrossProd(N, anisAxis);    // V is perpendicular to R and V - Y-axis

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

	// for outgoing direction
	//float multi = M_PI * a2;
	//cout << "3ps\n";

	float ax_2_pi_4 = 4.f * M_PI * slope_sq;

	for(int i = maxRays; i--; num++) 
	{
		// the ray is sampled in this direction
		CContributionItem *contrib = carray.AllocateNext();
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);

		float lightCos = -1, cos_HI = -1;
		CVector3D dirRes(0,0,1);
		if(slope >= 1e-03)
		{
			CVector3D hvector(0,0,1);
			float cos_HN = 0, cos2_HN = 0, sin_HN = 1, tan2_HN = 0;
			float sin_phi, cos_phi;
			float phi_HN;
			bool rejected = false;

			while(lightCos < 0) 
			{
				CVector2D rn; // random numbers
				if( !rejected )
					rn.SetValue(specular_rnd[2*num], specular_rnd[2*num+1]);
				else
					rn.SetValue(samplingInfo.Sample(0), samplingInfo.Sample(1));

				// according to (16),(17) in ref [1] :
				//									
				// delta = atan(sqrt(-log(1-s) * mx^2))
				//
				// where mx is the standard deviation of the surface slope and
				//
				// phiH = 2*pi*t;
				//
				// where s and t are independent random numbers uniformly distributed in [0, 1),
				//		 
				phi_HN			=  2.0 * M_PI * rn.yy;
				cos_phi			=  cos(phi_HN);
				sin_phi			=  sin(phi_HN);

				// we dont need compute exactly delta ( = HN), it is sufficient to compute tan^2(delta), from 
				// which is easy to compute sin(delta) and cos(delta) (less expensive than atan(sqrt()) and then
				// back to tan^2
				// 					  							
				// delta = atan(sqrt(-log(1-s)*alfa^2)) ---> tan^2(delta) = -log(1-s)*alfa^2
				//
				tan2_HN = (-log(1.0-rn.xx))*slope_sq;
				assert(tan2_HN >= 0);
				cos2_HN  = (1.0/(1.0 + tan2_HN));
				cos_HN   = sqrt(cos2_HN);
				sin_HN   = sqrt(1.0 - cos2_HN);

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

				CVector3D localDir(sin_HN * cos_phi, sin_HN * sin_phi, cos_HN);

				// Change the basis from the VUN coordinate system (local coordinate
				// system where normal R is Z-axis) to the World coordinate system
				hvector = localDir.x * anisAxis + localDir.y * V + localDir.z * N;

				hvector.Normalize();


				cos_HI = -DotProd(hvector,dirInit);

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

				if(cos_HI >= 0)
				{
					// dirRes = 2*cos(HI)*h - dirInit, but dirInit is negative in our case
					// (points toward the surface), so its - (-dirInit) = + dirInit
					dirRes = (2.0 * cos_HI) * hvector + dirInit; 

					assert(Magnitude(dirRes) >= 1.0-1e-04 && Magnitude(dirRes) <= 1.0+1e-04);
					lightCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

					if(lightCos < 0)
					{
						rejected = true;
					}
				}
				else //(lightCos < 0)
				{
					rejected = true;
				}
			}

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(lightCos >= 0);

			float cosHV = DotProd(-hitInfo.GetRay()->GetDir(),hvector);
			float k = 1.f / (ax_2_pi_4 * pow(cos_HN,3.f) * cosHV);

			// q(h)
			float q = exp( -tan2_HN / slope_sq );

			contrib->prob = k * q;

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

			CValues aux;
			hitInfo.GetTextureValue(aux,  *_alpha, CMappingBase::EE_KAN);
			float a = aux.Get(0);

			CValues F(CValues::EE_Weigth3);
			// compute fresnel equation for (w0.h)
			hitInfo.GetTextureValue(F,  *_specular, CMappingBase::EE_KFR, cosHV);

			float w = (F.Luminance() * pow(lightCos,1-a)) / (cos_HN * pow(incEyeCos,a));

			contrib->weigth = ks_local * w;

			// set the reflected ray, direction is normalized here
			contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
			// type of interaction
			contrib->interaction = EE_REFLECTED_SPECULAR;
			// diffuseness of the bounce
			contrib->diffuseness = diffuseness;
		}
		else // both slopes lower than 1e-3 -> mirror reflection
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			lightCos = DotProd(dirRes, N);
			contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
			// type of interaction
			contrib->interaction = EE_REFLECTED_MIRROR;
			// diffuseness of the bounce
			contrib->diffuseness = 0.0;

			// exp(-tan(0)/slope) = 1, cos(H,N) = 1 =>
			contrib->prob = 1.0 / (slope_sq * M_PI);

			// For mirror cosine of incident equal to the cosine of
			// reflected ray => sqrt(lightCos/incEyeCos) =>
			contrib->weigth = ks_local * 0.25;
		}

		if(_diffuse && (kd_local.Luminance() > 1e-05)) 
		{
		   //kd_local *= 1.0/(M_PI*contrib->prob);
		   contrib->weigth += kd_local * (lightCos/(M_PI*contrib->prob));
		}

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

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


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

	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);
	float slopef, slope_2f;

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

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

	if (_specular && (distType & EE_REFLECTED_SPECULAR)) 
	{
		hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);
		aks = ks_local.Luminance();

		if(!_slope)
		{
			WARNING_INFO("CSurfaceKurtKalosIsotropic::SpawnStochasticRays() :\n\
Slope structure should be set to compute specular part of Ward model.");
		  return 0;
		}

		CTexture *slope = _slope->GetSlopeX();
		if(!slope)
		{
			WARNING_INFO("CSurfaceKurtKalosIsotropic::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope' parameter.");
			return 0;
		}
		CValues slopeC(0,0,0);
		hitInfo.GetTextureValue(slopeC, *slope, CMappingBase::EE_KSI);
		slopef = slopeC.Get(0);
		(slopef < WTINY) ? slopef = WTINY : slopef = slopef;

		slope_2f = slopef * slopef;
		tmp_ratio     = aks/(slope_2f * 4.0 * M_PI);
	} // specular exists

	float albedo_whole = akd + aks;
  

	if(albedo_whole < 1e-05) 
	{
		WARNING_INFO("CSurfaceKurtKalosIsotropic::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
		return 0;
	}
	if( albedo_whole >= 1.0 ) 
	{
		WARNING_INFO("CSurfaceKurtKalosIsotropic::SpawnStochasticRays():\nAlbedo > 1.0\n\
To keep the energy conservation requirement the total albedo\n\
(albedo_specularity + albedo_diffuse) 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) || (tmp_ratio <= 1e-06)) 
		{
			// 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;

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

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

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

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

	samplingInfo.IncDimension(2);

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

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

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

const char* 
CSurfaceKurtKalosIsotropic::GetHelp() const
{
	return BrdfStrings::KKurtKalosIsoDesc;
}
//========================================================================

int 
CSurfaceKurtKalosIsotropic::GetNumParameters() const
{
	// diffuse + fresnel term + alpha normalizator + slope
	return 4; 
}
//========================================================================

CValues
CSurfaceKurtKalosIsotropic::GetParameter(int parameterIndex) const
{
	if (parameterIndex<4) // parent parameter
		return CSurfaceKurtKalosAnisotropic::GetParameter(parameterIndex);

	// for sure
	return CValues(3);
}

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

CValues
CSurfaceKurtKalosIsotropic::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex < 4) // parent parameter
		return CSurfaceKurtKalosAnisotropic::GetParameterConstraints(parameterIndex);

	return CValues(4);
}
//========================================================================

const char* 
CSurfaceKurtKalosIsotropic::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
		case 1 : 
		case 2 : return CSurfaceKurtKalosAnisotropic::GetParameterDescription(parameterIndex);
		case 3 : return BrdfStrings::KSlopeDesc;
		default : return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CSurfaceKurtKalosIsotropic::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex < 4) // parent parameter
	{
		return CSurfaceKurtKalosAnisotropic::GetParameterSyntax(parameterIndex);
	}
	return EE_Unknown;
}
//========================================================================

const char*
CSurfaceKurtKalosIsotropic::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
		case 1 : 
		case 2 : return CSurfaceKurtKalosAnisotropic::GetParameterTitle(parameterIndex);
		case 3 : return BrdfStrings::KSlopeTitle;
		default : return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CSurfaceKurtKalosIsotropic::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex < 4) // parent parameter
	{
		CSurfaceKurtKalosAnisotropic::SetParameter(parameterIndex, value);
	}
}

__END_BRDFLIB_SOURCE
