// ===================================================================
// $Id$
//
// ward.h
//    (Hemi)spherical harmonics representation of BRDFs.
//
// Class:
//
// Initial coding by Jaroslav Krivanek, January 2004.
// Edited by Jiri Nekola, October 2011

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

__BEGIN_BRDFLIB_SOURCE

typedef vector<double> rnd_array;
rnd_array CSurfaceMicrofacet::diffuse_rnd;
rnd_array CSurfaceMicrofacet::specular_rnd;




/**********************************************************************
 *
 *             CSurfaceWard methods
 *
 **********************************************************************/
// describes the surface to a given stream
//========================================================================
void
CSurfaceWard::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;
}

/*
  Integral of this function cannot be evaluated analitically. However, it can
  be evaluated using Taylor's series. 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
CSurfaceWard::GetAlbedo(CValues &albedo,
                          const CHitPointInfo &hitInfo,
                          EDistributionType   distType) const
{
  //CVector3D N       = hitInfo->GetNormal();
  //CVector3D dirInit = hitInfo->GetRay()->GetDir();

   // As default for all the Microfacet-based models. Actually, for each
  // model aldebo is a function of the angle bewtween the incoming
  // direction and normal.
  CValues c(CValues::EE_Weigth3); // result - zero by default
  CValues d(CValues::EE_Weigth3);
  bool retval = false;

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

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

  if (retval) {
    albedo = c;
    //cout << "albedo from GetAlbedo = " << albedo << "\n";
    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;
//  GetAlbedoMax(albedo, hitInfo, distType);
}

int
CSurfaceWard::SpawnStochasticRays(CContributionArray  &carray,
							const CHitPointInfo &hitInfo,
                            CSamplingInfo       &samplingInfo,
							EDistributionType	distType) const
{
	int maxRays = samplingInfo.GetMaxSamples();
	assert(maxRays > 0);

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

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

	// diffuse component
	CValues kd_local(CValues::EE_Weigth3);
	// specular component
	CValues ks_local(CValues::EE_Weigth3);
	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();//*(n+2.0)/(2.0*M_PI))*

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


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

	samplingInfo.IncDimension(2);

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

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

int 
CSurfaceWard::GetNumParameters() const
{
	return 4; // diffuse,specular, slope x(from CSurfaceMicrofacet), slope y
}
//========================================================================

CValues
CSurfaceWard::GetParameter(int parameterIndex) const
{
	CValues c(0,0,0);
	if (parameterIndex<3) // parent parameter
		return CSurfaceMicrofacet::GetParameter(parameterIndex);
	else if (parameterIndex == 3)
	{
		_slope->GetSlopeX()->GetValue(c,CHitPointInfo::constMapInfo);
	}
	return c;
}
//========================================================================

CValues
CSurfaceWard::GetParameterConstraints(int parameterIndex) const
{
	
    if (parameterIndex<2) // parent parameter
	  return CSurfaceMicrofacet::GetParameterConstraints(parameterIndex);

    if (parameterIndex == 2 || parameterIndex == 3) // slope x,y
      {
	CValues c(3);
	// usually, slope is weighted xyz parameter, but weight is not used in wards brdf 
	// at all, so its better to not confuse the user with slope's weight if it is not
	// used for computation anyway
	c.Set(0,0.001f); // min 
	c.Set(1,1); // max 
	c.Set(2,0.5f); // default 
	return c;
      }
    return CValues(0.f);
}
//========================================================================

const char* 
CSurfaceWard::GetParameterDescription(int parameterIndex) const
{
	if (parameterIndex<2) // parent parameter
		return CSurfaceMicrofacet::GetParameterDescription(parameterIndex);

	switch (parameterIndex)
	{
	case 2 :
		return BrdfStrings::KSlopeXDesc;
	case 3 :
		return BrdfStrings::KSlopeYDesc;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CSurfaceWard::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex<2) // parent parameter
		return CSurfaceMicrofacet::GetParameterSyntax(parameterIndex);
	switch (parameterIndex)
	{
	// usually, slope is weighted xyz parameter, but weight as well as y,z are not used 
	// in wards brdf at all, so its better to not confuse the user with these redundant
	// values, if they are not used for computation anyway
	case 2 : //slope[0] - x 
		return EE_Float;
	case 3 : //slope[1] - y
		return EE_Float;
	default :
		return EE_Unknown;
	}
}
//========================================================================

const char*
CSurfaceWard::GetParameterTitle(int parameterIndex) const
{
	if (parameterIndex<2) // parent parameter
		return CSurfaceMicrofacet::GetParameterTitle(parameterIndex);
	switch (parameterIndex)
	{
	case 2 :
		return BrdfStrings::KSlopeXTitle;
	case 3 :
		return BrdfStrings::KSlopeYTitle;
	default :
		return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CSurfaceWard::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex<2) // parent parameter
		return CSurfaceMicrofacet::SetParameter(parameterIndex, value);

	switch (parameterIndex)
	{
	case 2 :
		{
			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,0);
			_slope->SetWeight(w,0);
			break;
		}
	case 3 :
		{
			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,1);
			_slope->SetWeight(w,1);
			break;
		}
	}
}
//========================================================================



/**********************************************************************
 *
 *             CSurfaceWardAnisotropic methods
 *
 **********************************************************************/
//========================================================================

#define WTINY 1e-06
//========================================================================
bool
CSurfaceWardAnisotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		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_4 = 4.* M_PI * ax * ay;

				if(rho_s.Luminance()/ax_ay_pi_4 > 1e-06)
				{
					CVector3D  H  =  Normalize(dirIn - (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);
					hl.Normalize();
					//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;
					//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

					float ex  = exp( -tan2_HN * (cos2_phi/ax2 + sin2_phi/ay2));
					float k;
					if ( _max)
					{
						k = ex / (Max(fabs(lightCos),fabs(incEyeCos)) * ax_ay_pi_4);
					}
					else 
					{
						k =  ex / ( sqrt(lightCos*incEyeCos) * ax_ay_pi_4);
					}
					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Original Ward anisotropic BRDF :
//
//			 1	    exp(-tan^2(delta)*(cos^2(phi)/a^2 + sin^2(phi)/b^2)
//  fW =  ------ * ---------------------------------------------------
//		   Pi*a*b			4*sqrt(cos(thetaLight)*cos(thetaViewer))
//
//  where a,b - are the standard deviations of the surface slope, and
//  delta - is the angle between the half-vector and normal
//
//	PDF of such a BRDF by Arne Dür ("An Improved Normalization for the Ward Reflectance
//  Model", 2005)
//
//		  fW
//	PDF = --- 
//		   w
//
//										(cos(thetaI)+cos(thetaV))^3
// where w = ------------------------------------------------------------------------------------
//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
//
// where V = outgoing direction (viewer) and I = incident direction(light)
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceWardAnisotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_ay_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_ay_pi_4 > 1e-06) 
		{

			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);
			hl.Normalize();
			//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;
			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

			float ex  = exp( -tan2_HN * (cos2_phi/slopeXY_sq.xx +
                                   sin2_phi/slopeXY_sq.yy));
			float multi;
			if (_max)
			{
				//multi = Max(lightCos,incEyeCos);	// wrong, Max should be in denominator!
				multi = 1.f / Max(lightCos,incEyeCos);
			}
			else
			{
				// same fault as in SpawnStochasticSpecularRays !
				// 1/sqrt(cos(Out,N)*cos(Inc,N)) is not equal to sqrt(cos(Inc,N)/cos(Out,N))
				//multi = sqrt(lightCos/incEyeCos);
				multi = 1.f / sqrt(lightCos*incEyeCos);
			}
			float k   = ex * multi/( ax_ay_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceWardAnisotropic::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
	//CVector3D U = anisAxis; // U is perpendicular to R - X-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";

	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 outCos = -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(outCos < 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 corrected Ward sampling method :
				//									-log(1-s)
				// delta = atan(sqrt(-------------------------------------))
				//					 cos^2(phi)/alfa^2 + sin^2(phi)/beta^2
				//
				// where alfa and beta are the standard deviations of the surface slope and
				//
				// phi = atan2(beta*sin(2*pi*t),alfa*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);
					outCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

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

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= 0);
			//	PDF by Arne Dür ("An Improved Normalization for the Ward Reflectance Model", 2005)
			//
			//		  fW
			//	PDF = --
			//		  w
			//
			//										(cos(thetaI)+cos(thetaV))^3
			// where w = ------------------------------------------------------------------------------------
			//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
			//
			// where V = outgoing direction (viewer) and I = incident direction(light) and fW is Ward's BRDF

			float k = 1.f / (4 * M_PI * slopeXY.xx * slopeXY.yy);
			float brdf = (k * exp(-tan2_HN*anisPart));
			if (_max)
			{
				brdf /= Max(fabs(outCos),fabs(incEyeCos));
			}
			else
			{
				brdf /= sqrt(outCos*incEyeCos);
			}

			// l,v sins
			float incEyeSin = sqrt(1-incEyeCos*incEyeCos);
			float outSin = sqrt(1-outCos*outCos);

			// now get phi angles from output direction (dirRes) and input direction (eye)
			CVector3D incEyeDir = -hitInfo.GetRay()->GetDir();
			// we will compute phi angle as angle between the projection of the vector and some vector
			// laying in the plane with normal N -> e.g. anisAxis is such a vector
			// projection of the vector A onto the plane alfa with normal B is A||B
			// A||B = B × (A×B / |B|) / |B|
			// in the case both vectors are normalized
			// A||B = B × (A×B)
			CVector3D proj = CrossProd(N,CrossProd(incEyeDir,N));
			// now get cos of the mentioned angle

			float cosPhiInc = DotProd(proj,anisAxis);
			// now we have to compute sin of this angle..but be careful, we cannot compute sin(phi) as sqrt(1-cos^2(phi)),
			// because phi from cosinus is always in the range (0,pi> (dotproduct returns cos of the smallest angle between
			// given vectors so e.g. for 220° it returns cos for 360-220 = 140°). we have to compute sin(phi) carefully :
			// angle is between 0-pi, if cross product between the vectors is normal of the plane (right hand rule)
			// angle s between pi-2*pi(0), if cross product between the vectors is negative normal of the plane
			// so if the cross product between the given vectors on the plane is plane's negative normal, we have to invert
			// the sign of the sin(phi)
			float sinPhiInc = sqrt(1-cosPhiInc*cosPhiInc);
			// if angle between the vector perpendicular to projected vector and anisAxis is larger than some epsilon 
			// it means this perpendicular vector is negative normal and phi angle is in the range pi,2pi and we have to
			// invert the sin(phi)
			CVector3D perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiInc = -sinPhiInc;
			}

			// same for light direction (dirRes)
			proj = CrossProd(N,CrossProd(dirRes,N));
			float cosPhiLight = DotProd(proj,anisAxis);
			float sinPhiLight = sqrt(1-cosPhiLight*cosPhiLight);
			perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiLight = -sinPhiLight;
			}

			// cos(x-y) = cos(x)*cos(y) + sin(x)*sin(y)
			float cos_phii_phil = cosPhiInc*cosPhiLight + sinPhiInc*sinPhiLight;

			float w = pow(incEyeCos+outCos,3);
			w /= 4.f * incEyeCos * (1.f + outCos*incEyeCos + outSin*incEyeSin*cos_phii_phil);
			contrib->prob = brdf / w;

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

			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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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(outCos/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 * (outCos/(M_PI*contrib->prob));
		}

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

	} // for i

	
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceWardAnisotropic::GetHelp() const
{
	return BrdfStrings::KAnisoWardDesc;
}
//========================================================================

int 
CSurfaceWardAnisotropic::GetNumParameters() const
{
	return CSurfaceWard::GetNumParameters() + 1; // max ward modification as an addition to the ward base
}
//========================================================================

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

	CValues c(4);

	if (parameterIndex == 4)
	{
		c.Set(0,GetMaxWard());
	}
	return c;
}
//========================================================================

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

	CValues c(4);

	if (parameterIndex == 4)
	{
		_setBasicConstraints(EE_Boolean,c);
	}
	return c;
}
//========================================================================

const char* 
CSurfaceWardAnisotropic::GetParameterDescription(int parameterIndex) const
{
	if (parameterIndex < 4) // parent parameter
	{
		return CSurfaceWard::GetParameterDescription(parameterIndex);
	}

	if (parameterIndex == 4)
	{
		return BrdfStrings::KMaxWardDesc;
	}
		
	return BrdfStrings::KInvalidParameter;
}
//========================================================================

CSurface::EParameterUIType
CSurfaceWardAnisotropic::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex < 4) // parent parameter
	{
		return CSurfaceWard::GetParameterSyntax(parameterIndex);
	}
	
	if (parameterIndex == 4)
	{
		return EE_Boolean;
	}

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

const char*
CSurfaceWardAnisotropic::GetParameterTitle(int parameterIndex) const
{
	if (parameterIndex < 4) // parent parameter
	{
		return CSurfaceWard::GetParameterTitle(parameterIndex);
	}
	if (parameterIndex == 4)
	{
		return BrdfStrings::KMaxWardTitle;
	}
	
	return BrdfStrings::KInvalidParameter;
}
//========================================================================

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

	if (parameterIndex == 4)
	{
		bool boolValue = static_cast<bool>(value[0]);
		SetMaxWard(boolValue);
	}
}
//========================================================================




/**********************************************************************
 *
 *             CSurfaceWardIsotropic methods
 *
 **********************************************************************/

//========================================================================
bool
CSurfaceWardIsotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) 
		{
			CTexture *slope = _slope->GetSlopeX();
			if(slope)
			{
				CValues rho_s, slopeC;
				hitInfo.GetTextureValue(rho_s,  *_specular, CMappingBase::EE_KS);

				hitInfo.GetTextureValue(slopeC, *slope,   CMappingBase::EE_KSI);
				float ax  = slopeC.Get(0);
				if(ax < WTINY) ax = WTINY;
				float ax2 = ax * ax;
				float ax_2_pi_4 = 4.* M_PI * ax2;

				if(rho_s.Luminance()/ax_2_pi_4 > 1e-06)
				{
					CVector3D  H  =  Normalize(dirIn - (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;

					float ex  = exp( -tan2_HN / ax2 );
					float k;
					if ( _max)
					{
						k = ex / (Max(fabs(lightCos),fabs(incEyeCos)) * ax_2_pi_4);
					}
					else 
					{
						k =  ex / ( sqrt(lightCos*incEyeCos) * ax_2_pi_4);
					}
					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

int
CSurfaceWardIsotropic::SpawnStochasticRays(CContributionArray  &carray,
							const CHitPointInfo &hitInfo,
                            CSamplingInfo       &samplingInfo,
							EDistributionType	distType) const
{
	int maxRays = samplingInfo.GetMaxSamples();
	assert(maxRays > 0);

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

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

	// diffuse component
	CValues kd_local(CValues::EE_Weigth3);
	// specular component
	CValues ks_local(CValues::EE_Weigth3);
	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();//*(n+2.0)/(2.0*M_PI))*

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

		CTexture *slope = _slope->GetSlopeX();
		if(!slope)
		{
			WARNING_INFO("CSurfaceWardIsotropic::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope'");
			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("CSurfaceWardIsotropic::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
		return 0;
	}
	if( albedo_whole >= 1.0 ) 
	{
		WARNING_INFO("CSurfaceWardIsotropic::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;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Original Ward isotropic BRDF :
//
//			 1				exp(-tan^2(delta)*a^2)
//  fW =  ------ * ------------------------------------------
//		   Pi*a^2	4*sqrt(cos(thetaLight)*cos(thetaViewer))
//
//  where a,b - are the standard deviations of the surface slope, and
//  delta - is the angle between the half-vector and normal
//
//	PDF of such a BRDF by Arne Dür ("An Improved Normalization for the Ward Reflectance
//  Model", 2005)
//
//		  fW
//	PDF = --- 
//		   w
//
//										(cos(thetaI)+cos(thetaV))^3
// where w = ------------------------------------------------------------------------------------
//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
//
// where V = outgoing direction (viewer) and I = incident direction(light)
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceWardIsotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_2_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_2_pi_4 > 1e-06) 
		{

			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;

			float ex  = exp( -tan2_HN / slope_sq );
			float multi;
			if (_max)
			{
				//multi = Max(lightCos,incEyeCos);	// wrong, Max should be in denominator!
				multi = 1.f / Max(lightCos,incEyeCos);
			}
			else
			{
				// same fault as in SpawnStochasticSpecularRays !
				// 1/sqrt(cos(Out,N)*cos(Inc,N)) is not equal to sqrt(cos(Inc,N)/cos(Out,N))
				//multi = sqrt(lightCos/incEyeCos);
				multi = 1.f / sqrt(lightCos*incEyeCos);
			}
			float k   = ex * multi/( ax_2_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceWardIsotropic::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
	//CVector3D U = anisAxis; // U is perpendicular to R - X-axis

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

	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 outCos = -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(outCos < 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 corrected Ward sampling method, for isotropic version stands:
				//									
				// delta = atan(sqrt(-log(1-s) * alfa^2))
				//
				// where alfa is the standard deviation of the surface slope and
				//
				// phi = (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);
					outCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

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

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= 0);
			//	PDF by Arne Dür ("An Improved Normalization for the Ward Reflectance Model", 2005)
			//
			//		  fW
			//	PDF = --
			//		  w
			//
			//										(cos(thetaI)+cos(thetaV))^3
			// where w = ------------------------------------------------------------------------------------
			//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
			//
			// where V = outgoing direction (viewer) and I = incident direction(light) and fW is Ward's BRDF

			float brdf;
			float ax_2_pi_4 = 4 * M_PI * slope_sq;
			float ex = exp(-tan2_HN / slope_sq);
			if ( _max)
			{
				brdf = ex / (Max(fabs(outCos),fabs(incEyeCos)) * ax_2_pi_4);
			}
			else 
			{
				brdf =  ex / ( sqrt(outCos*incEyeCos) * ax_2_pi_4);
			}

			// l,v sins
			float incEyeSin = sqrt(1-incEyeCos*incEyeCos);
			float outSin = sqrt(1-outCos*outCos);

			// now get phi angles from output direction (dirRes) and input direction (eye)
			CVector3D incEyeDir = -hitInfo.GetRay()->GetDir();
			// we will compute phi angle as angle between the projection of the vector and some vector
			// laying in the plane with normal N -> e.g. anisAxis is such a vector
			// projection of the vector A onto the plane alfa with normal B is A||B
			// A||B = B × (A×B / |B|) / |B|
			// in the case both vectors are normalized
			// A||B = B × (A×B)
			CVector3D proj = CrossProd(N,CrossProd(incEyeDir,N));
			// now get cos of the mentioned angle

			float cosPhiInc = DotProd(proj,anisAxis);
			// now we have to compute sin of this angle..but be careful, we cannot compute sin(phi) as sqrt(1-cos^2(phi)),
			// because phi from cosinus is always in the range (0,pi> (dotproduct returns cos of the smallest angle between
			// given vectors so e.g. for 220° it returns cos for 360-220 = 140°). we have to compute sin(phi) carefully :
			// angle is between 0-pi, if cross product between the vectors is normal of the plane (right hand rule)
			// angle s between pi-2*pi(0), if cross product between the vectors is negative normal of the plane
			// so if the cross product between the given vectors on the plane is plane's negative normal, we have to invert
			// the sign of the sin(phi)
			float sinPhiInc = sqrt(1-cosPhiInc*cosPhiInc);
			// if angle between the vector perpendicular to projected vector and anisAxis is larger than some epsilon 
			// it means this perpendicular vector is negative normal and phi angle is in the range pi,2pi and we have to
			// invert the sin(phi)
			CVector3D perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiInc = -sinPhiInc;
			}

			// same for light direction (dirRes)
			proj = CrossProd(N,CrossProd(dirRes,N));
			float cosPhiLight = DotProd(proj,anisAxis);
			float sinPhiLight = sqrt(1-cosPhiLight*cosPhiLight);
			perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiLight = -sinPhiLight;
			}

			// cos(x-y) = cos(x)*cos(y) + sin(x)*sin(y)
			float cos_phii_phil = cosPhiInc*cosPhiLight + sinPhiInc*sinPhiLight;

			float w = pow(incEyeCos+outCos,3);
			w /= 4.f * incEyeCos * (1.f + outCos*incEyeCos + outSin*incEyeSin*cos_phii_phil);
			contrib->prob = brdf / w;

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

			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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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 / (slope_sq * M_PI);

			// For mirror cosine of incident equal to the cosine of
			// reflected ray => sqrt(outCos/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 * (outCos/(M_PI*contrib->prob));
		}

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

	} // for i

	
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceWardIsotropic::GetHelp() const
{
	return BrdfStrings::KIsoWardDesc;
}
//========================================================================

int 
CSurfaceWardIsotropic::GetNumParameters() const
{
	// isotropic ward has only one slope parameter, not two
	return CSurfaceWardAnisotropic::GetNumParameters() - 1; 
}
//========================================================================

CValues
CSurfaceWardIsotropic::GetParameter(int parameterIndex) const
{
	if (parameterIndex<3) // parent parameter
		return CSurfaceWardAnisotropic::GetParameter(parameterIndex);

	if (parameterIndex == 3)
	{
		// skip slope Y
		return CSurfaceWardAnisotropic::GetParameter(parameterIndex+1);
	}

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

CValues
CSurfaceWardIsotropic::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
		return CSurfaceWardAnisotropic::GetParameterConstraints(parameterIndex);

	if (parameterIndex == 3)
	{
		// skip slope Y
		return CSurfaceWardAnisotropic::GetParameterConstraints(parameterIndex+1);
	}
	return CValues(4);
}
//========================================================================

const char* 
CSurfaceWardIsotropic::GetParameterDescription(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
		case 1 : return CSurfaceWardAnisotropic::GetParameterDescription(parameterIndex);
		case 2 : return BrdfStrings::KSlopeDesc;
		// skip slope Y
		case 3 : return CSurfaceWardAnisotropic::GetParameterDescription(parameterIndex+1);

		default : return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

CSurface::EParameterUIType
CSurfaceWardIsotropic::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
	{
		return CSurfaceWardAnisotropic::GetParameterSyntax(parameterIndex);
	}
	
	if (parameterIndex == 3)
	{
		// skip slope Y
		return CSurfaceWardAnisotropic::GetParameterSyntax(parameterIndex+1);
	}

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

const char*
CSurfaceWardIsotropic::GetParameterTitle(int parameterIndex) const
{
	switch (parameterIndex)
	{
		case 0 :
		case 1 : return CSurfaceWardAnisotropic::GetParameterTitle(parameterIndex);
		case 2 : return BrdfStrings::KSlopeTitle;
		// skip slope Y
		case 3 : return CSurfaceWardAnisotropic::GetParameterTitle(parameterIndex+1);

		default : return BrdfStrings::KInvalidParameter;
	}
}
//========================================================================

void
CSurfaceWardIsotropic::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex < 3) // parent parameter
	{
		CSurfaceWardAnisotropic::SetParameter(parameterIndex, value);
	}

	if (parameterIndex == 3)
	{
		// skip slope Y
		CSurfaceWardAnisotropic::SetParameter(parameterIndex+1, value);
	}
}
//========================================================================




/**********************************************************************
 *
 *             CSurfaceWardDuerAnisotropic methods
 *
 **********************************************************************/

//========================================================================
bool
CSurfaceWardDuerAnisotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		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_4 = 4.* M_PI * ax * ay;

				if(rho_s.Luminance()/ax_ay_pi_4 > 1e-06)
				{
					//  exp(-tan^2(theta_half)/slopeX^2) / (4 * Pi * slopeX^2)
					//                                  / sqrt(cos(theta_in)*cos(theta_out))
					//if(lightCos      < WTINY) lightCos       = WTINY;
					//if(incEyeCos < WTINY) incEyeCos  = WTINY;

					CVector3D  H  =  Normalize(dirIn - (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);
					hl.Normalize();
					//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;
					//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

					float ex = exp( -tan2_HN * (cos2_phi/ax2 + sin2_phi/ay2));
					float k = ex / ( (lightCos*incEyeCos) * ax_ay_pi_4);
					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Ward-Dür anisotropic BRDF :
//
//			 1	    exp(-tan^2(delta)*(cos^2(phi)/a^2 + sin^2(phi)/b^2)
//  fWD =  ------ * ---------------------------------------------------
//		   Pi*a*b			4*cos(thetaLight)*cos(thetaViewer)
//
//  where a,b - are the standard deviations of the surface slope, and
//  delta - is the angle between the half-vector and normal
//
//	PDF of such a BRDF by Arne Dür ("An Improved Normalization for the Ward Reflectance
//  Model", 2005)
//
//		  fWD
//	PDF = --- 
//		   w
//
//										(cos(thetaI)+cos(thetaV))^3
// where w = ------------------------------------------------------------------------------------
//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
//
// where V = outgoing direction (viewer) and I = incident direction(light)
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceWardDuerAnisotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_ay_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_ay_pi_4 > 1e-06) 
		{

			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);
			hl.Normalize();
			//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;
			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

			float ex  = exp( -tan2_HN * (cos2_phi/slopeXY_sq.xx +
                                   sin2_phi/slopeXY_sq.yy));
			float multi = 1.f / (lightCos*incEyeCos);
			float k   = ex * multi/( ax_ay_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
		//if (_refract && (distType & EE_REFRACTED_PERFECT)) {
		//  hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
		//  contrib->weigth += ks_local;
		//}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceWardDuerAnisotropic::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;
	//int countSamples = 0;

	float incEyeCos    = hitInfo.GetDotProd();
	//cout << "incEysCos = " << incEyeCos;
	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 N as Z axis
	CVector3D V = CrossProd(N, anisAxis);   

	float anisPart = 1;


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

	// sampling according to paper : 
	// David Geisler-Moroder, Arne Dür. A New Ward BRDF Model with Bounded Albedo. Comput. Graph. Forum, 2010: 1391~1398   
 

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

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

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

				// according to corrected Ward sampling method :
				//									-log(1-s)
				// delta = atan(sqrt(-------------------------------------))
				//					 cos^2(phi)/alfa^2 + sin^2(phi)/beta^2
				//
				// where alfa and beta are the standard deviations of the surface slope and
				//
				// phi = atan2(beta*sin(2*pi*t),alfa*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;
				float phi_HN	=  atan2f(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
				h= localDir.x * anisAxis + localDir.y * V + localDir.z * N;

				h.Normalize();


				cos_HI = -DotProd(h,dirInit);

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

				if(cos_HI >= 0)
				{
					// make final direction vector
					// 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) * h + dirInit;

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

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

		
			// Set the value of BSDF type
			contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
			
			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= -WTINY);

			// constant part in BSDF denominator
			float ax_ay_pi_4 = 4.* M_PI * slopeXY.xx * slopeXY.yy;

			// WEIGHT .........................................................................
			//
			//							                  cos^3(theta_l + theta_v) 
			// weight = -------------------------------------------------------------------------------------------------
			//			(4 * cos(theta_v) * (1.0 + cos(theta_l)*cos(theta_v) + sin(theta_l)*sin(theta_v)*cos(phi_v-phi_l)
			//
			// where v = viewer (input direction), l = light (generated output direction)

			// numerator of the weight
			float w_numer = pow(incEyeCos + outCos, 3.0f);

			// l,v sins
			float incEyeSin = sqrt(1-incEyeCos*incEyeCos);
			float outSin = sqrt(1-outCos*outCos);
			
			// For translation vectors 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_4 > 1e-06) 
			{
				MakeLocalCS(N, anisAxis, lx, ly, lz);
			}
			// now get phi angles from output direction (dirRes) and input direction (eye)
			CVector3D incEyeDir = -hitInfo.GetRay()->GetDir();
			// we will compute phi angle as angle between the projection of the vector and some vector
			// laying in the plane with normal N -> e.g. anisAxis is such a vector
			// projection of the vector A onto the plane alfa with normal B is A||B
			// A||B = B × (A×B / |B|) / |B|
			// in the case both vectors are normalized
			// A||B = B × (A×B)
			CVector3D proj = CrossProd(N,CrossProd(incEyeDir,N));
			// now get cos of the mentioned angle

			float cosPhiInc = DotProd(proj,anisAxis);
			// now we have to compute sin of this angle..but be careful, we cannot compute sin(phi) as sqrt(1-cos^2(phi)),
			// because phi from cosinus is always in the range (0,pi> (dotproduct returns cos of the smallest angle between
			// given vectors so e.g. for 220° it returns cos for 360-220 = 140°). we have to compute sin(phi) carefully :
			// angle is between 0-pi, if cross product between the vectors is normal of the plane (right hand rule)
			// angle s between pi-2*pi(0), if cross product between the vectors is negative normal of the plane
			// so if the cross product between the given vectors on the plane is plane's negative normal, we have to invert
			// the sign of the sin(phi)
			float sinPhiInc = sqrt(1-cosPhiInc*cosPhiInc);
			// if angle between the vector perpendicular to projected vector and anisAxis is larger than some epsilon 
			// it means this perpendicular vector is negative normal and phi angle is in the range pi,2pi and we have to
			// invert the sin(phi)
			CVector3D perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiInc = -sinPhiInc;
			}

			// same for light direction (dirRes)
			proj = CrossProd(N,CrossProd(dirRes,N));
			float cosPhiLight = DotProd(proj,anisAxis);
			float sinPhiLight = sqrt(1-cosPhiLight*cosPhiLight);
			perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiLight = -sinPhiLight;
			}

			// cos(x-y) = cos(x)*cos(y) + sin(x)*sin(y)
			float cos_phii_phil = cosPhiInc*cosPhiLight + sinPhiInc*sinPhiLight;
			
			// denominator of the weight
			float w_denom = 4.0 * incEyeCos * (1.0 + incEyeCos*outCos + incEyeSin*outSin*cos_phii_phil);
			
			// -----------------------------------------------------------------------------------

			// PDF of the Ward-Dür :
			// PDF = fWD() / weight = Ward's specular BRDF / weight
			// see EvaluateBSDF() method for more info about BRDF
			
			cos_HN = DotProd(h,N);
			tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;

			float ex  = exp( -tan2_HN * anisPart);
			float multi = 1.f / (outCos*incEyeCos);
			float k   = (ex * multi) / ax_ay_pi_4; // == Ward Dür BRDF


			// weight = ks * w
			contrib->weigth = ks_local * (w_numer/w_denom);

			// pdf = brdf / w
			contrib->prob = (k * w_denom) / w_numer;

			// 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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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(outCos/incEyeCos) =>
			contrib->weigth = ks_local * 0.25;
		}

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

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

	} // for i
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceWardDuerAnisotropic::GetHelp() const
{
	return BrdfStrings::KAnisoWardDuerDesc;
}

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




/**********************************************************************
 *
 *             CSurfaceWardDuerIsotropic methods
 *
 **********************************************************************/
//========================================================================
/*
  Integral of this function cannot be evaluated analitically. However, it can
  be evaluated using Taylor's series. 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
CSurfaceWardDuerIsotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) 
		{
			CTexture *slope = _slope->GetSlopeX();
			if(slope)
			{
				CValues rho_s, slopeC;
				hitInfo.GetTextureValue(rho_s,  *_specular, CMappingBase::EE_KS);

				hitInfo.GetTextureValue(slopeC, *slope,   CMappingBase::EE_KSI);
				float ax  = slopeC.Get(0);
				if(ax < WTINY) ax = WTINY;
				float ax2 = ax * ax;
				float ax_2_pi_4 = 4.* M_PI * ax2;

				if(rho_s.Luminance()/ax_2_pi_4 > 1e-06)
				{
					CVector3D  H  =  Normalize(dirIn - (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;

					float ex  = exp( -tan2_HN / ax2 );
					float k = ex / ((lightCos*incEyeCos) * ax_2_pi_4);

					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

int
CSurfaceWardDuerIsotropic::SpawnStochasticRays(CContributionArray  &carray,
							const CHitPointInfo &hitInfo,
                            CSamplingInfo       &samplingInfo,
							EDistributionType	distType) const
{
	int maxRays = samplingInfo.GetMaxSamples();
	assert(maxRays > 0);

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

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

	// diffuse component
	CValues kd_local(CValues::EE_Weigth3);
	// specular component
	CValues ks_local(CValues::EE_Weigth3);
	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();//*(n+2.0)/(2.0*M_PI))*

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

		CTexture *slope = _slope->GetSlopeX();
		if(!slope)
		{
			WARNING_INFO("CSurfaceWardDuerIsotropic::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope'");
			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("CSurfaceWardDuerIsotropic::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
		return 0;
	}
	if( albedo_whole >= 1.0 ) 
	{
		WARNING_INFO("CSurfaceWardDuerIsotropic::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;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Ward-Dür isotropic BRDF :
//
//			 1			  exp(-tan^2(delta)*a^2)
//  fW =  ------ * --------------------------------------
//		   Pi*a^2	4*(cos(thetaLight)*cos(thetaViewer))
//
//  where a,b - are the standard deviations of the surface slope, and
//  delta - is the angle between the half-vector and normal
//
//	PDF of such a BRDF by Arne Dür ("An Improved Normalization for the Ward Reflectance
//  Model", 2005)
//
//		  fW
//	PDF = --- 
//		   w
//
//										(cos(thetaI)+cos(thetaV))^3
// where w = ------------------------------------------------------------------------------------
//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
//
// where V = outgoing direction (viewer) and I = incident direction(light)
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceWardDuerIsotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_2_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_2_pi_4 > 1e-06) 
		{

			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;

			float ex  = exp( -tan2_HN / slope_sq );
			float multi = 1.f / (lightCos*incEyeCos);

			float k   = ex * multi/( ax_2_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceWardDuerIsotropic::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
	//CVector3D U = anisAxis; // U is perpendicular to R - X-axis

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

	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 outCos = -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(outCos < 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 corrected Ward sampling method, for isotropic version stands:
				//									
				// delta = atan(sqrt(-log(1-s) * alfa^2))
				//
				// where alfa is the standard deviation of the surface slope and
				//
				// phi = (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);
					outCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

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

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= 0);
			//	PDF by Arne Dür ("An Improved Normalization for the Ward Reflectance Model", 2005)
			//
			//		  fW
			//	PDF = --
			//		  w
			//
			//										(cos(thetaI)+cos(thetaV))^3
			// where w = ------------------------------------------------------------------------------------
			//			 4*cos(thetaV)*(1 + cos(thetaI)*cos(thetaV) + sin(thetaI)*sin(thetaV)*cos(phiV-phiI))
			//
			// where V = outgoing direction (viewer) and I = incident direction(light) and fW is Ward's BRDF

			
			float ax_2_pi_4 = 4 * M_PI * slope_sq;
			float ex = exp(-tan2_HN / slope_sq);
			float brdf = ex / ((outCos*incEyeCos) * ax_2_pi_4);

			// l,v sins
			float incEyeSin = sqrt(1-incEyeCos*incEyeCos);
			float outSin = sqrt(1-outCos*outCos);

			// now get phi angles from output direction (dirRes) and input direction (eye)
			CVector3D incEyeDir = -hitInfo.GetRay()->GetDir();
			// we will compute phi angle as angle between the projection of the vector and some vector
			// laying in the plane with normal N -> e.g. anisAxis is such a vector
			// projection of the vector A onto the plane alfa with normal B is A||B
			// A||B = B × (A×B / |B|) / |B|
			// in the case both vectors are normalized
			// A||B = B × (A×B)
			CVector3D proj = CrossProd(N,CrossProd(incEyeDir,N));
			// now get cos of the mentioned angle

			float cosPhiInc = DotProd(proj,anisAxis);
			// now we have to compute sin of this angle..but be careful, we cannot compute sin(phi) as sqrt(1-cos^2(phi)),
			// because phi from cosinus is always in the range (0,pi> (dotproduct returns cos of the smallest angle between
			// given vectors so e.g. for 220° it returns cos for 360-220 = 140°). we have to compute sin(phi) carefully :
			// angle is between 0-pi, if cross product between the vectors is normal of the plane (right hand rule)
			// angle s between pi-2*pi(0), if cross product between the vectors is negative normal of the plane
			// so if the cross product between the given vectors on the plane is plane's negative normal, we have to invert
			// the sign of the sin(phi)
			float sinPhiInc = sqrt(1-cosPhiInc*cosPhiInc);
			// if angle between the vector perpendicular to projected vector and anisAxis is larger than some epsilon 
			// it means this perpendicular vector is negative normal and phi angle is in the range pi,2pi and we have to
			// invert the sin(phi)
			CVector3D perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiInc = -sinPhiInc;
			}

			// same for light direction (dirRes)
			proj = CrossProd(N,CrossProd(dirRes,N));
			float cosPhiLight = DotProd(proj,anisAxis);
			float sinPhiLight = sqrt(1-cosPhiLight*cosPhiLight);
			perp = CrossProd(anisAxis,proj);
			perp.Normalize();
			if (DotProd(perp,N) < 0) // > 90°
			{
				sinPhiLight = -sinPhiLight;
			}

			// cos(x-y) = cos(x)*cos(y) + sin(x)*sin(y)
			float cos_phii_phil = cosPhiInc*cosPhiLight + sinPhiInc*sinPhiLight;

			float w = pow(incEyeCos+outCos,3);
			w /= 4.f * incEyeCos * (1.f + outCos*incEyeCos + outSin*incEyeSin*cos_phii_phil);
			contrib->prob = brdf / w;

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

			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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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 / (slope_sq * M_PI);

			// For mirror cosine of incident equal to the cosine of
			// reflected ray => sqrt(outCos/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 * (outCos/(M_PI*contrib->prob));
		}

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

	} // for i

	
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceWardDuerIsotropic::GetHelp() const
{
	return BrdfStrings::KIsoWardDuerDesc;
}
//========================================================================

int 
CSurfaceWardDuerIsotropic::GetNumParameters() const
{
	// isotropic ward-dür has only one slope parameter, not two
	return CSurfaceWardDuerAnisotropic::GetNumParameters() - 1; 
}
//========================================================================

CValues
CSurfaceWardDuerIsotropic::GetParameter(int parameterIndex) const
{
	if (parameterIndex<3) // parent parameter
		return CSurfaceWardDuerAnisotropic::GetParameter(parameterIndex);

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

CValues
CSurfaceWardDuerIsotropic::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
		return CSurfaceWardDuerAnisotropic::GetParameterConstraints(parameterIndex);

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

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

CSurface::EParameterUIType
CSurfaceWardDuerIsotropic::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
	{
		return CSurfaceWardDuerAnisotropic::GetParameterSyntax(parameterIndex);
	}

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

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

void
CSurfaceWardDuerIsotropic::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex < 3) // parent parameter
	{
		CSurfaceWardDuerAnisotropic::SetParameter(parameterIndex, value);
	}
}
//========================================================================




/**********************************************************************
 *
 *             CSurfaceMoroderDuerAnisotropic methods
 *
 **********************************************************************/

//========================================================================
bool
CSurfaceMoroderDuerAnisotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		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_4 = 4.* M_PI * ax * ay;

				if(rho_s.Luminance()/ax_ay_pi_4 > 1e-06)
				{
					CVector3D  H  =  Normalize(dirIn - (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;

					float cos_LH  =  DotProd(dirIn, H);
					assert(cos_LH >= 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);
					hl.Normalize();
					//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;
					//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

					float ex = exp( -tan2_HN * (cos2_phi/ax2 + sin2_phi/ay2));
					float k = ex / ( pow(cos_LH,2) * pow(cos_HN,4) * ax_ay_pi_4);
					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Moroder-Dür anisotropic BRDF :
//
//			 1	    exp(-tan^2(delta)*(cos^2(phi)/a^2 + sin^2(phi)/b^2)
//  fMD =  ------ * ---------------------------------------------------
//		   Pi*a*b			   4*cos^2(lh)*cos^4(delta)
//
//  where a,b - are the standard deviations of the surface slope, and
//  delta is the angle between the half-vector and normal and lh is angle between
//  half-vector and light direction (== angle between light direction and viewer
//  direction / 2)
//
//	PDF of such a BRDF according ot the paper (using original Ward sampling method) is
//
//		  fMD
//	PDF = --- 
//		   w
//
//					     2
// where w = --------------------------
//			 1 + thetaViewer/thetaLight
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceMoroderDuerAnisotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_ay_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_ay_pi_4 > 1e-06) 
		{

			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;

			float cos_LH = DotProd(dirRes,H);

			// 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);
			hl.Normalize();
			//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;
			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

			float ex  = exp( -tan2_HN * (cos2_phi/slopeXY_sq.xx +
                                   sin2_phi/slopeXY_sq.yy));
			float multi = 1.f / (pow(cos_LH,2) * pow(cos_HN,4));
			float k   = ex * multi/( ax_ay_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
		//if (_refract && (distType & EE_REFRACTED_PERFECT)) {
		//  hitInfo.GetTextureValue(kr_local, *_refract, CMappingBase::EE_KR);
		//  contrib->weigth += ks_local;
		//}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceMoroderDuerAnisotropic::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;
	//int countSamples = 0;

	float incEyeCos    = hitInfo.GetDotProd();
	//cout << "incEysCos = " << incEyeCos;
	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 N as Z axis
	CVector3D V = CrossProd(N, anisAxis);   

	float anisPart = 1;


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

	// sampling according to paper : 
	// David Geisler-Moroder, Arne Dür. A New Ward BRDF Model with Bounded Albedo. Comput. Graph. Forum, 2010: 1391~1398   
 

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

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

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

				// according to corrected Ward sampling method :
				//									-log(1-s)
				// delta = atan(sqrt(-------------------------------------))
				//					 cos^2(phi)/alfa^2 + sin^2(phi)/beta^2
				//
				// where alfa and beta are the standard deviations of the surface slope and
				//
				// phi = atan2(beta*sin(2*pi*t),alfa*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;
				float phi_HN	=  atan2f(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
				h = localDir.x * anisAxis + localDir.y * V + localDir.z * N;

				h.Normalize();


				cos_HI = -DotProd(h,dirInit);

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

				if(cos_HI >= 0)
				{
					// make final direction vector
					// 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) * h + dirInit;

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

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

		
			// Set the value of BSDF type
			contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
			
			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= -WTINY);

			// constant part in BSDF denominator
			float ax_ay_pi_4 = 4.* M_PI * slopeXY.xx * slopeXY.yy;

			// WEIGHT --------------------------------------------------------------------------
			//				   2
			// w = --------------------------
			//	   1 + thetaViewer/thetaLight

			float w = 2.f / (1.f + incEyeCos/outCos);
			
			// ----------------------------------------------------------------------------------
			// PDF of the Moroder-Dür :
			// PDF = fMD() / weight = Moroder-Dür specular BRDF / weight
			// see EvaluateBSDF() method for more info about BRDF
			
			cos_HN = DotProd(h,N);
			tan2_HN =  1.0f/(cos_HN * cos_HN) - 1.0;
			float cos_LH = DotProd(dirRes,h);

			float ex  = exp( -tan2_HN * anisPart);
			float multi = 1.f / (pow(cos_LH,2) * pow(cos_HN,4));
			float brdf   = (ex * multi) / ax_ay_pi_4; // == Moroder-Dür BRDF


			// weight = ks * w
			contrib->weigth = ks_local * w;

			// pdf = brdf / w
			contrib->prob = brdf / 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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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(outCos/incEyeCos) =>
			contrib->weigth = ks_local * 0.25;
		}

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

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

	} // for i
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceMoroderDuerAnisotropic::GetHelp() const
{
	return BrdfStrings::KAnisoMoroderDuerDesc;
}





/**********************************************************************
 *
 *             CSurfaceMoroderDuerIsotropic methods
 *
 **********************************************************************/
//========================================================================
/*
  Integral of this function cannot be evaluated analitically. However, it can
  be evaluated using Taylor's series. 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
CSurfaceMoroderDuerIsotropic::EvaluateBSDF(CValues &result,
			 const CHitPointInfo &hitInfo,
			 const CVector3D     &dirIn,
			 EDistributionType   distType) const
{
	bool retval = false;
	float d = Magnitude(dirIn);
	CVector3D N          =  hitInfo.GetNormal();
	CVector3D anisAxis   =  hitInfo.GetAnisAxis();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06);

	float lightCos       =  DotProd(dirIn, N);
	float incEyeCos  =  -(hitInfo.GetDotProd());

	if(lightCos >= -1e-06 && incEyeCos >= -1e-06)
	{
		if(lightCos < 0) lightCos = 0;
		else
		{
			if(lightCos < WTINY) lightCos       = WTINY;
		}

		if(incEyeCos < 0) incEyeCos = 0;
		else
		{
			if(incEyeCos < WTINY) incEyeCos  = WTINY;
		}

		if ( (d < 1.0-1e-06) || (d > 1.0+1e-06) )
		{
			Normalize(dirIn);
		}

		if(_diffuse && (distType & EE_REFLECTED_DIFFUSE)) 
		{
		  CValues rho_d;
		  hitInfo.GetTextureValue(rho_d, *_diffuse, CMappingBase::EE_KS);
		  result = rho_d * M_1_PI;
		  retval = true;
		}

		if((_specular && (distType & EE_REFLECTED_GLOSSY)) && _slope) 
		{
			CTexture *slope = _slope->GetSlopeX();
			if(slope)
			{
				CValues rho_s, slopeC;
				hitInfo.GetTextureValue(rho_s,  *_specular, CMappingBase::EE_KS);

				hitInfo.GetTextureValue(slopeC, *slope,   CMappingBase::EE_KSI);
				float ax  = slopeC.Get(0);
				if(ax < WTINY) ax = WTINY;
				float ax2 = ax * ax;
				float ax_2_pi_4 = 4.* M_PI * ax2;

				if(rho_s.Luminance()/ax_2_pi_4 > 1e-06)
				{
					CVector3D  H  =  Normalize(dirIn - (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;

					float cos_LH  =  DotProd(dirIn, H);
					assert(cos_LH >= 0);

					float ex  = exp( -tan2_HN / ax2 );
					float k = ex / ( pow(cos_LH,2) * pow(cos_HN,4) * ax_2_pi_4);

					result += rho_s * k;
					retval = true;
				} // lum > 1e-6
			} // slopes exist
		} // specular exist

		if(retval)
		{
			result *= lightCos;
		}
	} // not under the surface
	else // under the surface
	{
		retval = false;
	}
	return retval;
}

int
CSurfaceMoroderDuerIsotropic::SpawnStochasticRays(CContributionArray  &carray,
							const CHitPointInfo &hitInfo,
                            CSamplingInfo       &samplingInfo,
							EDistributionType	distType) const
{
	int maxRays = samplingInfo.GetMaxSamples();
	assert(maxRays > 0);

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

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

	// diffuse component
	CValues kd_local(CValues::EE_Weigth3);
	// specular component
	CValues ks_local(CValues::EE_Weigth3);
	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("CSurfaceMoroderDuerIsotropic::SpawnStochasticRays() :\n\
Slope structure should be set to compute specular part of Ward model.");
		  return 0;
		}

		CTexture *slope = _slope->GetSlopeX();
		if(!slope)
		{
			WARNING_INFO("CSurfaceMoroderDuerIsotropic::SpawnStochasticRays() :\n\
Glossy term requires structure 'slope'");
			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("CSurfaceMoroderDuerIsotropic::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
		return 0;
	}
	if( albedo_whole >= 1.0 ) 
	{
		WARNING_INFO("CSurfaceMoroderDuerIsotropic::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;
}

//-------------------------------------------------------------------------------------------------
// Proportional sampling is implemented according to sampling scheme in
// "A NewWard BRDF Model with Bounded Albedo" paper, written by David 
// Geisler-Moroder and Arne Dür in 2010. Original sampling scheme from Ward92
// was incorrect (arctan for phi angle was omitted) and then lately corrected by
// Walter05 and Arne Dur. Final sampling method was nicely recapitulated in the
// mentioned paper describing bounded-albedo Ward model :
//
// Moroder-Dür isotropic BRDF :
//
//			 1	      exp(-tan^2(delta)*a^2)
//  fMD =  ------ * --------------------------
//		   Pi*a*b	 4*cos^2(lh)*cos^4(delta)
//
//  where a is the standard deviations of the surface slope, and
//  delta is the angle between the half-vector and normal
//  lh is the angle between light direction and half vector
//
//	PDF of such a BRDF by Arne Dür ("An Improved Normalization for the Ward Reflectance
//  Model", 2005)
//
//		  fW
//	PDF = --- 
//		   w
//
//					     2
// where w = --------------------------
//			 1 + thetaViewer/thetaLight
//
//--------------------------------------------------------------------------------------------------

//========================================================================
int
CSurfaceMoroderDuerIsotropic::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();

	assert(fabs(DotProd(N, anisAxis)) < 1e-06); // angle ~ pi/2

	float incEyeCos    = hitInfo.GetDotProd();

	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

	// number of rays created
	int num = 0;
	float ax_2_pi_4 = 4.* 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_4 > 1e-06) 
	{
		MakeLocalCS(N, anisAxis, lx, ly, lz);
	}

	for(int i = maxRays; i--; num++) 
	{
		// allocate next contribution
		CContributionItem *contrib = carray.AllocateNext();
		// generate the sample in the cannonical space
		CVector3D dirN;
		// Generate a sample in 2D, range [0;1)x[0;1)
		CVector2D vv(diffuse_rnd[2*num], diffuse_rnd[2*num+1]);

		// probability is set to value cos(theta)/M_PI
		contrib->prob = RND::hemisph_cos(dirN, vv);
		// transform to real coordinates
		CVector3D dirRes = dirN.x * U + dirN.y * V + dirN.z * N ;
		// set the reflected ray
		contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
		// set the flags concerning ray type
		contrib->SetFlags(CContributionItem::EE_ReflectedRay);
		// type of interaction
		contrib->interaction = EE_REFLECTED_DIFFUSE;
		// diffuseness of the bounce
		contrib->diffuseness = 1.0;
		// Set the value of BSDF type
		// compute = BRDF * cos(normal, ray.dir) / pdf =
		//         = pd/Pi * cos(normal, ray.dir) /(cos(normal, ray.dir)/Pi)
		//         = pd
		// This means that we should not multiply the weight by the cosine
		// of the angle between the incoming light direction and normal
		contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
		// Cosine between the incoming angle and normal is already taken into
		// account!
		contrib->weigth = kd_local; // * M_1_PI;


		if(ks_local.Luminance()/ax_2_pi_4 > 1e-06) 
		{

			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;

			float cos_LH = DotProd(dirRes,H);

			float ex  = exp( -tan2_HN / slope_sq );
			float multi = 1.f / (pow(cos_LH,2) * pow(cos_HN,4));

			float k   = ex * multi/( ax_2_pi_4 );
			k *= 1.0/contrib->prob;

			contrib->weigth += ks_local * k;
   
		}
#ifdef __TRAVERSAL_STATISTICS
    world->_statistics->IncReflectedRays();
#endif

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

int
CSurfaceMoroderDuerIsotropic::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
	//CVector3D U = anisAxis; // U is perpendicular to R - X-axis

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

	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 outCos = -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(outCos < 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 corrected Ward sampling method, for isotropic version stands:
				//									
				// delta = atan(sqrt(-log(1-s) * alfa^2))
				//
				// where alfa is the standard deviation of the surface slope and
				//
				// phi = (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);
					outCos = DotProd(N, dirRes);  // (angle betwrrn the light ray and normal)

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

			incEyeCos = fabs(incEyeCos);
			if(incEyeCos < WTINY) incEyeCos = WTINY;
			assert(outCos >= 0);
			
			
			// ----------------------------------------------------------------------------------
			// PDF of the Moroder-Dür :
			// PDF = fMD() / weight = Moroder-Dür specular BRDF / weight
			// see EvaluateBSDF() method for more info about BRDF

			// WEIGHT --------------------------
			//				   2
			// w = --------------------------
			//	   1 + thetaViewer/thetaLight
			float w = 2.f / (1.f + incEyeCos/outCos);

			float cos_LH = DotProd(dirRes,hvector);

			float ax_2_pi_4 = 4 * M_PI * slope_sq;
			float ex = exp(-tan2_HN / slope_sq);
			float multi = 1.f / (pow(cos_LH,2) * pow(cos_HN,4));
			float brdf = (ex * multi) / ax_2_pi_4;

			contrib->prob = brdf / w;

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

			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 
		{
			// set the reflected ray, direction is normalized here
			dirRes = hitInfo.GetReflectedDir();
			outCos = 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 / (slope_sq * M_PI);

			// For mirror cosine of incident equal to the cosine of
			// reflected ray => sqrt(outCos/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 * (outCos/(M_PI*contrib->prob));
		}

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

	} // for i

	
  assert(num == maxRays);
  return num;
}

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

const char* 
CSurfaceMoroderDuerIsotropic::GetHelp() const
{
	return BrdfStrings::KIsoMoroderDuerDesc;
}
//========================================================================

int 
CSurfaceMoroderDuerIsotropic::GetNumParameters() const
{
	// isotropic ward-dür has only one slope parameter, not two
	return CSurfaceMoroderDuerAnisotropic::GetNumParameters() - 1; 
}
//========================================================================

CValues
CSurfaceMoroderDuerIsotropic::GetParameter(int parameterIndex) const
{
	if (parameterIndex<3) // parent parameter
		return CSurfaceMoroderDuerAnisotropic::GetParameter(parameterIndex);

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

CValues
CSurfaceMoroderDuerIsotropic::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
		return CSurfaceMoroderDuerAnisotropic::GetParameterConstraints(parameterIndex);

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

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

CSurface::EParameterUIType
CSurfaceMoroderDuerIsotropic::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex < 3) // parent parameter
	{
		return CSurfaceMoroderDuerAnisotropic::GetParameterSyntax(parameterIndex);
	}

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

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

void
CSurfaceMoroderDuerIsotropic::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex < 3) // parent parameter
	{
		CSurfaceMoroderDuerAnisotropic::SetParameter(parameterIndex, value);
	}
}
//========================================================================


__END_BRDFLIB_SOURCE
