// ===================================================================
//
// coupledModel.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 "coupledModel.h"
#include "rnd.h"
#include "texture.h"
#include <iostream>

#define constant 21/(20*M_PI)

__BEGIN_BRDFLIB_SOURCE

// kd_local is albedo => it should be divided by PI to be used
//   as the value of BRDF.
//========================================================================
int
CSurfaceCoupledModel::SpawnStochasticDiffuseRays(
  const CValues       &kd_local,
  const CValues       &shininess,
  CContributionArray  &carray,
  const CHitPointInfo &hitInfo,
  CSamplingInfo       &samplingInfo,
  int                 maxRays) const
{
  CVector3D N = hitInfo.GetNormal();
  if (hitInfo.GetDotProd() > 0.0)
    N = -N;

  int num = 0;			        // number of rays created
  float theta, fi, lightCos;

  // Compute data for maxRays outgoing rays
  //----------------------------------------------------------------------
  for(int i = maxRays; i--; num++) {
     // allocate next contribution
    CContributionItem *contrib = carray.AllocateNext();
    fi = 2*M_PI*diffuse_rnd[2*num];
    float eps = diffuse_rnd[2*num + 1];
    float upper = (38295.9*pow(eps,3) + 13002.4*pow(eps,2) + 154.4 * eps + 0.0114813);
    float lower = (16520.2*pow(eps,3) + 33596.4*pow(eps,2) + 1483.57*eps + 1);
    theta = acos(upper / lower);
    //theta = M_PI*diffuse_rnd[2*num+1];

    /*float cosTheta = cos(theta);
    float temp = 7/2*pow(cosTheta,3) - 21/4*pow(cosTheta,4) + 21/5*pow(cosTheta,5)
                - 7/4*pow(cosTheta,6) + 3/10*pow(cosTheta,7);
    if (diffuse_rnd[2*num + 1] + 0.01 < temp || diffuse_rnd[2*num+1] - 0.01 > temp){
     cout << "need Newton iteration method" << endl;
    }*/


    CVector3D dirRes;
    float one = 1;
    PolarToCartesian(N, fi, theta, one, dirRes);

    lightCos = DotProd(N, dirRes);
    // Reject only samples which are outside the hemisphere
    if(lightCos < 0.0){
       //rejected = true;
       dirRes = -dirRes;
       lightCos = DotProd(N, dirRes);
    }

    // probability is set to value cos(theta)/M_PI
    contrib->prob = 100;//(1-pow(1-cos(theta),5))*cos(theta)*sin(theta);
    // 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
    contrib->weigth.SetRepresentation(CValues::EE_Weigth3);
    contrib->weigth = 21/(20*M_PI)*(1-shininess.Get(0))*
          (1-pow(1-lightCos,5))/sin(theta) *kd_local;

    // Additional (specular) component in order to provide nice and simple
    // visualization
    //------------------------------------------------------------------------
    CValues ks_local(CValues::EE_Weigth3);
    if (_specular) {
      hitInfo.GetTextureValue(ks_local, *_specular, CMappingBase::EE_KS);

      if(ks_local.Luminance() > 1e-5) {
        float alphaCos = DotProd(dirRes,hitInfo.GetReflectedDir());

	// This checking is necessary, otherwise, samples, which are outside
	// the cosine lobe whould be incorrectly involved into the model
	if(alphaCos > 0.0)
	{
          CValues shininess;
          hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
          float n = shininess.Get(0);
          float lightCos = DotProd(N, dirRes);

  	  ks_local *= (n+2.0)*lightCos*pow(alphaCos,n)/(2.0*M_PI*contrib->prob);
          contrib->weigth += ks_local;
	}
      }
    }

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

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


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

  diffuse_rnd.clear();
  specular_rnd.clear();

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

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

  // compute diffuse, specular, and refractive albedos
  float akd = 0.0f, aks = 0.0f, 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("CSurfaceCoupledModel::SpawnStochasticRays() :\nAlbedo of this surface is < 1e-05. It is a black body!");
    return 0;
  }
  if( albedo_whole >= 1.0 ) {
	WARNING_INFO("CSurfaceCoupledModel::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   |
   *  +----------------------------------------------------+
   *                        (s2)              (s3)
   */

  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 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   |
     *  +----------------------------------------------------+
     *  |<----- alpha2 ----->|<---- alpha3 --->|
     */
     float alpha1 = akd/albedo_whole;
     float alpha2 = aks/albedo_whole;

     if(rnd_num < alpha1 && alpha1 > 1e-6) {
       // 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);
       diffuse_rnd.push_back(samplingInfo.Sample(1));
       numDiffuse++;
     }
     else {
       if(rnd_num < (alpha1+alpha2) && alpha2 > 1e-6) {
         // 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();
  }

  int numRays = 0;

  CValues shininess;
  hitInfo.GetTextureValue(shininess, *_specnum, CMappingBase::EE_KSI);
  // The number of rays for diffuse interaction
  //----------------------------------------------------------------------
  if (numDiffuse) {
    numRays += SpawnStochasticDiffuseRays(kd_local, shininess, carray,
					  hitInfo, samplingInfo,
					  numDiffuse);
  }

  int indexDiffuse = carray.GetUsedItems() ;

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

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

    }
    mult /= albedo_whole;
  } // for*/
  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
CSurfaceCoupledModel::EvaluateBSDF(CValues &bsdf,
				    const CHitPointInfo &hitInfo,
				    const CVector3D &dirLit,
				    EDistributionType distType) const
{
  bool retval       = false;
  CVector3D normal  = hitInfo.GetNormal();
  CVector3D dirInit = (hitInfo.GetRay())->GetDir();
  float incEyeCos   = -(hitInfo.GetDotProd());
  float lightCos    = DotProd(normal, dirLit);

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

    // now get the R0 constant
    CValues specularIndex(CValues::EE_Weigth3);
    if (_specnum)
      hitInfo.GetTextureValue(specularIndex, *_specnum,
                         CMappingBase::EE_KSI);
    else
      // supposing that if not specified, surface is diffuse
      specularIndex.Set(0.2);

    // Diffuse part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if ((distType & EE_REFLECTED_DIFFUSE) && _diffuse) {
      CValues diffuseTerm(bsdf.GetType());
      hitInfo.GetTextureValue(diffuseTerm, *_diffuse, CMappingBase::EE_KD);
	  // k = 21 / 20*pi*(1-specnum) = constant / (1-specnum)
	  // Rm = diffuse component
	  // diffuse part : k * Rm * (1 - (1-cos(in))^5) * (1 - (1-cos(out))^5)
      bsdf = (constant / (1-specularIndex.Get(0))) *(1-pow(1-lightCos,5))
            *(1-pow(1-incEyeCos,5)) * diffuseTerm;
	  // edited 16-10-2011 Nekola : constant 21/(20*M_1_PI) is wrong and is not corrresponding to 
	  // given equation k = 21 / 20*pi*(1-specnum) = constant / (1-specnum)
	  // in that way 21/(20*(1/pi)) = 21/(20/pi) = 21pi/20 and not 21/20pi 
	  // correct constant k should be 21/(20*M_PI)

      retval = true;
    }

    // Specular part
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    CValues specularTerm(bsdf.GetType());

    // if surface has specular component
    if ((distType & EE_REFLECTED_GLOSSY) && _specular) {

      //if (hitInfo.GetReflectedDir() == -dirInit) ???? 16-10-2011, Nekola
	  {
        // get the value of specular component
        hitInfo.GetTextureValue(specularTerm, *_specular, CMappingBase::EE_KS);

        specularTerm *= specularIndex.Get(0) + pow(1-incEyeCos,5)*(1-specularIndex.Get(0))
                        *Max(incEyeCos, lightCos);

        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* 
CSurfaceCoupledModel::GetHelp() const
{
	return BrdfStrings::KCoupledModelDesc;
}
//========================================================================


int 
CSurfaceCoupledModel::GetNumParameters() const
{
	return 3; // diffuse,specular,specnum
}
//========================================================================

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

CValues
CSurfaceCoupledModel::GetParameterConstraints(int parameterIndex) const
{
	if (parameterIndex == 2) // specnum is different from phong brdf (possible values between 0 and 1)
	{
		CValues c(9);
		c.Set(0,0); // min R
		c.Set(1,0.99f); // max R
		c.Set(2,0.1f); // default R
		c.Set(3,0); // min G
		c.Set(4,0.99f); // max G
		c.Set(5,0.1f); // default G
		c.Set(6,0); // min B
		c.Set(7,0.99f); // max B
		c.Set(8,0.1f); // default B
		return c;
	}
	// this just sets aside phong parameters that are not used in this model
	if (parameterIndex >= 0 && parameterIndex < 3)
		return CSurfacePhongBase::GetParameterConstraints(parameterIndex);
	else 
		return CValues(0,0,0);
}
//========================================================================

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

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

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

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