// ===================================================================
//
// anisotropicPhong.cpp
//     Implementation of Anisotropic Phong shading model (Asikhmin&Shirley).
//
// Licence and usage conditions in November 2003:
// This file can be used and modified only with written permission of
// Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
// This library is granted to be used by Polina Kondratieva within
// her MPI HIWI job working on BRDFLIB library. The usage condition
// may change in the next version of the library. At the moment,
// the library development is supported by Max-Planck-Institut for
// Informatik, Stuhlsatzenhausweg 85, 66123 Saarbruecken Germany
// See WWW: http://www.mpi-sb.mpg.de
// REPLACEMENT_STRING
//
// Initial coding by Hana Truskova & Michael Smutny
// as a term work for X36RSO 2008/2009 at Prague CTU FEL
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


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

__BEGIN_BRDFLIB_SOURCE

// This variable is used to store intermediate result of random number
// generation if Selection Roulette is used.

typedef vector<double> rnd_array;

#define FAKTOR 1.217391304 //28/23
#define M_1_8PI M_1_PI/8


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

  int num = 0;			    // the number of samples generated
  CVector3D dirInit = -hitInfo.GetRay()->GetDir();  //normalized ray direction
  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;
  const CVector3D &Rm = hitInfo.GetReflectedDir();
  CVector3D R = Rm;

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

  float nu = shininess.Get(0);
  float nv = shininess.Get(1);


  // Compute directions for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
    float lightCos = -1;
    CVector3D dirRes(0.f);
    //bool rejected = false;
    float theta, fi;
    CVector3D h;
    // Rejection sampling for Phong model. It can happens, that random
    // direction will be generated under the surface. In this case we should
    // reject the sample and make the new one (in other case sample will be
    // considered as absorbed, but it contradicts with the fact that
    // exactly maxRays must be spawned by this kind of surface interaction).
    while(lightCos < 0) {
      // first we generate half vector h
      //TODO: ensure sampling over all of the hemisphere - flip psi along the PI/2 axis
      //to ensure stratification
       if ( specular_rnd[2*num] > 0.25 && specular_rnd[2*num] < 0.5 )
          specular_rnd[2*num] = 1 - 4*(0.5 - specular_rnd[2*num]);
        else {
       if ( specular_rnd[2*num] > 0.5 && specular_rnd[2*num] < 0.75)
          specular_rnd[2*num] = 1 - 4*(0.75 - specular_rnd[2*num]);
       else {
        if ( specular_rnd[2*num] > 0.75)
          specular_rnd[2*num] = 1 - 4*(1 - specular_rnd[2*num]);
        else
          specular_rnd[2*num] = 1 - 4*(0.25 - specular_rnd[2*num]);
       }
       }
       // uniform distribution
       //fi = 2.0*M_PI*specular_rnd[2*num];
       //theta = M_PI*specular_rnd[2*num+1];
       fi = atan(sqrt((nu+1)/(nv+1))*tan(M_PI_2*specular_rnd[2*num]));
       theta = acos(pow((float)(1-specular_rnd[2*num+1]),(1/(1+nu*pow(cos(fi),2)+nv*pow(sin(fi),2)))));
       fi *= 4;
       theta *= 2;
       float one = 1.0;
       PolarToCartesian(N, fi, theta, one, h);
       h.Normalize();

       dirRes = CVector3D(-dirInit + 2*DotProd(dirInit,h)*h);
       dirRes.Normalize();


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

      //float exp = nu*pow(cos(fi),2)+nv*pow(sin(fi),2);
      //float upper = sqrt((nu+1)*(nv+1))/M_PI;
      //float next = pow(DotProd(N,h),exp);
      //float lower = (2*DotProd(dirInit,h));
      // compute the probability density ot the outgoing direction
      contrib->prob = 1;//upper*next / lower; //seems to be in percentage

      // Set the value of BSDF type
      contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
      /* compute = value_of_BRDF / probability_density_of_generated direction
         or more precisely:
         compute = BRDF / pdf
         Note, ps is the albedo which is given as input parameter. */
      //contrib->weigth = lightCos / contrib->prob * kslocal;//
      contrib->weigth = 1;//1/Max(DotProd(N,dirInit),DotProd(N,dirRes))
                      //*(kslocal + (1-kslocal)*pow(1-DotProd(dirInit,h),5));

      // set the reflected ray, direction is normalized here
      contrib->ray.Init(hitInfo.GetPoint(), dirRes, REFLECTED_RAY);
      // set the flags concerning ray type
      contrib->SetFlags(CContributionItem::EE_ReflectedRay);
      // type of interaction
      contrib->interaction = EE_REFLECTED_SPECULAR;

      // Additional (diffuse) components in order to provide nice and simple
      // vizualisation
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      CValues kd_local(CValues::EE_Weigth3);
      if (_diffuse) {
        hitInfo.GetTextureValue(kd_local, *_diffuse, CMappingBase::EE_KD);
      }

      if(kd_local.Luminance() > 1e-05) {
        kd_local *= 1.0/(M_PI*contrib->prob);
        contrib->weigth += kd_local;
      }
      // Cosine between the incoming ray and the normal does not taken
      // into account.
      contrib->weigth *= lightCos;

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

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


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

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

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

  // diffuse component
  CValues kd_local(CValues::EE_Weigth3);
  // specular component
  CValues ks_local(CValues::EE_Weigth3), ks_inc_local(CValues::EE_Weigth3);
  // refractive component (ideal refraction)

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

  hitInfo.GetNormal();
  EDistributionType tmp_distType;

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

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


  // Complete albedo of physically correct Phong model
  //----------------------------------------------------------------------
  float albedo_whole =  akd + aks;

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

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

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

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

     // This RND number is from the User's program. Actually, it should be
     // fit into the range 0..albedo_whole if some event besides absorbtion
     // was choosen. Therefore, we shoul make RND fall into the [0..1)
     if(count >= maxRays-1)
     {
       if(rnd_num <= albedo_whole)
         rnd_num /= albedo_whole;  // make first RND fall into the [0..1)
     }

    /* This is the rule for choosing the interaction:
     *
     *  +----------------------------------------------------+
     *  |  reflect diffusely |    specularly   |
     *  +----------------------------------------------------+
     *  |<----- alpha1 ----->|<---- alpha2 --->|
     */
     float alpha1 = akd/albedo_whole;
     float alpha2 = aks/albedo_whole;

     if(rnd_num < alpha1 && alpha1 > 1e-5) {
       // Reuse random number for generation of one component of the random
       // direction. Scale it to fall into the range [0..1)
       diffuse_rnd.push_back(rnd_num/alpha1);
       // Choose another component of direction randomly
       diffuse_rnd.push_back(samplingInfo.Sample(1));
       numDiffuse++;
     }
     else {
          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();
  }

  int numRays = 0;

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

  //int indexDiffuse = carray.GetUsedItems() ;

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

  samplingInfo.IncDimension(2);

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

  // WARNING:
  // Be careful with this number of rays!!!!! Initial Phong model have no
  // possibility to add the transmitivity to the object. We suppose to use
  // deterministic refractive rays. But this causes some difficultiues:
  // only one ray will be produced for given hit point does not matter
  // which part of the light is refracted.
  //
  // There should be shot maxRays number of rays, but it will be shot less
  // number of rays because of it is not necessary to shoot the same ray
  // several times - it is just wasting of time. Therefore, if the result
  // sum of all rays' contributions for a given pixel will be divided by
  // maxRays the pixel color will be dimmer than it should be => the sum
  // of all rays' contributions for a given pixel should be divided by
  // numRays -> actual number of rays.      06/12/03/ KP.

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

// Direct lighting computation
//========================================================================
bool
CSurfaceAnisotropicPhong::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &k1,
				    EDistributionType distType) const
{
  bool retval       = false;
  CVector3D n  = hitInfo.GetNormal(); //normal
  CVector3D k2  = -((hitInfo.GetRay())->GetDir());
  CVector3D h = (k2 + k1)/2;
  h.Normalize();
  float incEyeCos   = -(hitInfo.GetDotProd());
  float lightCos    = DotProd(n, k1);

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

    float dotNK1 = DotProd(n,k1);
    float dotNK2 = DotProd(n,k2);

    // Diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      CValues diffuseTerm(bsdf.GetType());
      hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);
      
      //28/23*PI*Rd*(1-Rs)*(1-(1-(n*k1)/2)^5)*(1-(1-(n*k2)/2)^5)
      //k1...normalized vector to light
      //k2...normalized vector to viewer
      bsdf = FAKTOR * M_1_PI
            *(1-pow((1-dotNK1/2),5))*(1-pow((1-dotNK2/2),5)) *  diffuseTerm;
      //multiplication 1-Rs is in the following block
      retval = true;
    }

    // Specular part
    // sqrt((nu + 1)*(nv + 1))/M_8_PI * (n*h)^()...
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CValues specularTerm(bsdf.GetType());

    // if surface has specular component
    if ((distType & EE_REFLECTED_GLOSSY) && _specular) {
      CValues white(CValues::EE_Weigth3);
      white.Set(0, 1); white.Set(1, 1);  white.Set(2, 1);
      hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KS);
      if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse){
          bsdf*= (white-specularTerm); // from the diffuse block
      }

      float nu, nv;
      // now get the index of imperfection for specular reflection
      CValues specularIndex(CValues::EE_Weigth3);
      if (_specnum){
        hitInfo.GetTextureValue(specularIndex, *_specnum,
	                         CMappingBase::EE_KSI);
        nu = specularIndex.Get(0);
        nv = specularIndex.Get(1);
        retval = true;
      }
      else {
        // supposing that if not specified, surface is diffuse
          if (retval) bsdf *= lightCos;
          return retval;
      }

      CVector3D u = ArbitraryNormal(n); // ortonormal basis - n,u,v
      CVector3D v = CrossProd(n, u);
      
      float dotNH = DotProd(n,h);
      //float dotNU = DotProd(n,u);
      //float dotNV = DotProd(n,v);
      float dotHV = DotProd(h,v);
      float dotHU = DotProd(h,u);
      float dotHK1 = DotProd(h,k1);
      specularTerm = sqrt((nu+1)*(nv+1))*M_1_8PI*pow(dotNH,(nu*pow(dotHU,2)+nv*pow(dotHV,2)/(1-pow(dotNH,2))))
        /(dotHK1*max(DotProd(n,k1),DotProd(n,k2)))*(specularTerm+(white-specularTerm)*pow(1-(dotHK1),5));


      bsdf += specularTerm;
      retval = true;
    } // specular term

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

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

const char* 
CSurfaceAnisotropicPhong::GetHelp() const
{
	return BrdfStrings::KAnisoPhongDesc;
}
//========================================================================


int 
CSurfaceAnisotropicPhong::GetNumParameters() const
{
	return 3; // diffuse,specular,specnum // only used parameters from phong base
}
//========================================================================

CValues
CSurfaceAnisotropicPhong::GetParameter(int parameterIndex) const
{
	// this just sets aside phong parameters that are not used in this anisotropic kind
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::GetParameter(parameterIndex);
	else 
		return CValues(0,0,0);
}

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

CValues
CSurfaceAnisotropicPhong::GetParameterConstraints(int parameterIndex) const
{
	// this just sets aside phong parameters that are not used in this anisotropic kind
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::GetParameterConstraints(parameterIndex);
	else 
		return CValues(0,0,0);
}

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

const char* 
CSurfaceAnisotropicPhong::GetParameterDescription(int parameterIndex) const
{
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::GetParameterDescription(parameterIndex);
	else
		return BrdfStrings::KInvalidParameter;
}
//========================================================================

CSurface::EParameterUIType
CSurfaceAnisotropicPhong::GetParameterSyntax(int parameterIndex) const
{
	if (parameterIndex >= 0 && parameterIndex < 2)
		return CSurfacePhongBase::GetParameterSyntax(parameterIndex);
	if (parameterIndex == 2)
		return EE_XY; // specnum anisotropic
	else
		return EE_Unknown;
}
//========================================================================

const char*
CSurfaceAnisotropicPhong::GetParameterTitle(int parameterIndex) const
{
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::GetParameterTitle(parameterIndex);
	else
		return BrdfStrings::KInvalidParameter;
}
//========================================================================

void
CSurfaceAnisotropicPhong::SetParameter(int parameterIndex, const CValues& value)
{
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::SetParameter(parameterIndex,value);
}
//========================================================================
