// ===================================================================
// $Id: texture.h,v 1.8 2003/06/17 08:44:15 havran Exp $
//
// texture.h
//     Header file for CTexture class hierarchy.
//
// Class: CTextureList, CTexture, CConstantTexture
//
// Licence and usage conditions in November 2003:
// This file can be used and modified only with written permission of
// Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
// This library is granted to be used by Polina Kondratieva within
// her MPI HIWI job working on BRDFLIB library. The usage condition
// may change in the next version of the library. At the moment,
// the library development is supported by Max-Planck-Institut for
// Informatik, Stuhlsatzenhausweg 85, 66123 Saarbruecken Germany
// See WWW: http://www.mpi-sb.mpg.de
// REPLACEMENT_STRING
//
// Initial coding by Vlastimil Havran + Petr Mladek + Jiri Liska, 1999
//
// Given for the use to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


#ifndef __TEXTURE_H__
#define __TEXTURE_H__

// BRDFLIB headers
#include "configh.h"
#include "color.h"
#include "mappibas.h"

// standard headers
#include <iostream>
#include <vector>


__BEGIN_BRDFLIB_HEADER

// forward declarations
class CMapPointInfo;
class CTexture;

// Reference:
//  1. Complex refraction:     http://www.luxpop.com
//  2. Shirley, "Realistic Ray Tracing", 2003
//  3. McCool, "Reflectance lectures", 2002
//  4. Neumann, A. Ph.D. "Construction of Bidirectional Reflection Distribution
//     Functions", 2001
//  5. Glassner, "Principles of Digital Image Synthesis", Vol. 2
//  6. Dutre, P. "Global Illumination Compendium", 2001
//     http://www.graphics.cornell.edu/~phil/GI
//  7. Schlick, C. "An Inexpensive BRDF Model for Physically-based
//     Rendering", 1994
 	  
// ----------------------------------------------------------------
// CTextureList
//  List of textures; derives from CSimpleList and adds a function
//  to search for a given color and return nonzero if there.
// ----------------------------------------------------------------
class CTextureList:
  public vector<CTexture *>
{
public:
  // returns true, if the texture is already present in the list
  bool TextureThere(CTexture *probe);
};

// --------------------------------------------------------------------------
// This is general class for textures
// most important is metod GetValue which comes in for the mapInfo and 
// returns appropriate value
// --------------------------------------------------------------------------
class CTexture 
{
public:
  // this identifikator allows to identify real class of the object CImage *
  // 0x01xxxx classes which calculate mapInfo on the texture value
  // 0x02xxxx classes which operate with the value from the one texture
  // 0x04xxxx classes which operate with the values from the two textures
  // 0x08xxxx classes which operate with the values from many textures
  // 0x10xxxx classes which change mapInfo
  enum EClassId {
    // unknown, not recognized texture
    EE_UnknownTexture = 0,
    
    // classes which calculate mapInfo on the texture value
    EE_ConstantTexture = 0x010000,
    EE_FresnelTexture = 0x020000,
    EE_GeneralImageTexture,
    EE_InterpolatedTexture,
    EE_MipMapTexture,
    EE_NearestNeighbourTexture,
    EE_RipMapTexture,
    EE_SummedAreaTexture,
  };

private:
  // these must not be used
  CTexture & operator=(const CTexture &texture);
  // copy constructor
  CTexture(const CTexture &);
protected:
  // if the type (dimensionality) of the texture is static then is true
  // by other words if the demand on the dimensionality of the mapping
  // is static then this variable is true
  bool _isTypeStatic;
public:
  // constructor
  CTexture(bool isTypeStatic = true) : _isTypeStatic(isTypeStatic) { }

  // explicitely virtual destructor
  virtual ~CTexture() { }

  // returns the type (dimensionality) of the texture (2D,3D, ...).
  // It corresponds to the mapping that has to be performed to access the data
  // from this texture. So the dimensionality of mapping performed has to be
  // larger or greater than the dimensionality of this texture.
  virtual CMappingBase::EDimType GetType() const = 0;

  // returns true if the type of the texture is static
  // by other words it returns true if the demand on the dimensionality
  // of the mapping is static
  bool IsTypeStatic() const { return _isTypeStatic; }

  // returns identification of the class
  virtual EClassId GetClassId() const = 0;

  // returns one value for specifed mapInfo
  virtual void GetValue(CValues &value, const CMapPointInfo &mapInfo) const = 0;

  // Returns one value for specifed array of mapInfos.
  // The class mapInfos defines area instead of only one point.
  // This function can better solve antialiasing.
  virtual void GetValue(CValues &value, const CMapPointInfo mapInfo[],
			int number) const = 0;

  // returns refraction parameters (index of refraction, coefficient of
  // refraction) if their were set (for computing of Fresnel reflection)
  virtual bool GetRefractionParams(CValues &ior, CValues &extinction,
                                   bool &approx, const CMapPointInfo &) const = 0;

  // Returns true is this texture is of FresnelTexture type
  virtual bool IsFresnel() = 0;

  // Describes the texture to a given stream
  virtual void Describe(ostream &app, int ind, float multiplier = 1.0) const {
    indent(app, ind); multiplier = multiplier; app << "Not yet implemented";
  }
};

// --------------------------------------------------------------------------
// this is class for constant texture, this texture needn't mapping
// --------------------------------------------------------------------------
class CConstantTexture:
  public CTexture
{
public:
  // constructor
  CConstantTexture(const CValues &value):
    CTexture(), _value(value) { }

  CConstantTexture(const CConstantTexture &masterTexture):
    CTexture(),
    _value(masterTexture._value) { }

  // returns the type (dimensionality) of the texture (2D,3D, ...).
  virtual CMappingBase::EDimType GetType() const { return CMappingBase::EE_None; }

  // returns identification of the class
  virtual EClassId GetClassId() const { return EE_ConstantTexture; }

  // set the new color to constant texture
  void SetColor(const CValues &newColor) {
    _value = newColor;
  }

  // set the new color to constant texture
  void GetColor(CValues &color) {
    color = _value;
  }

  // returns value for specifed mapInfo
  virtual void GetValue(CValues &value, const CMapPointInfo &) const {
    value = _value;
  }

  virtual void GetValue(CValues &value, const CMapPointInfo [], int) const {
    value = _value;
 }

  // return value of complex index of refraction (index of refraction and
  // coeffitient of extinction) for specified mapInfo. Added by KP
  virtual bool GetRefractionParams(CValues &ior, CValues &extinction,
                                   bool &approx, const CMapPointInfo &) const {
    ior        = CValues(0, 0, 0);
    extinction = CValues(0, 0, 0);
    approx     = false;
    return false;
  }

  // return true if texture is of Fresnel type, otherwise false
  virtual bool IsFresnel() {return false;}

  CConstantTexture & operator=(const CConstantTexture &texture) {
    if (this != &texture)
      this->_value = texture._value;
    return *this;
  }

  // Describes the texture to a given stream
  virtual void Describe(ostream &app, int ind, float multiplier = 1.0) const {
    indent(app, ind); app << _value * multiplier;
  }
protected:
  // value of this texture for each texture coordinates
  CValues _value;
};

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


// ==========================================================================
// This is class for Fresnel texture, and for a now it needn't mapping
// The value of the texture depends only on the angle between the incoming
// light direction and normal. By default complex refraction coefficient [1]
// corresponds to that of ideal mirror (n = 0, k > 0)
// There are two types of materials: conductors and dielectrics (isotropic
// case).
//  * For conductors index of refraction is complex and all the transmitted
//    light is adsorbed.
//  * For dielectrics index of refraction is real and the most part of the
//    light is transmitted. Brewster angle is the critical angle between the
//    incomning direction and normal. If light transmitts through the
//    dielectric media towards the outside (outside media has refractive
//    index less than dielectric) under the Brewster angle, all the light is
//    reflected along the media boundary. If angle is greater than Brewster
//    angle, the total internal reflection is occurs - light is reflected
//    in mirror direction into the dielectric.
//
// abridgements: ior     - index of refraction
//               rel.ior - relative index of refraction = ior2/ior1
//
// added by KP, March 2004
// ==========================================================================
class CFresnelTexture:
  public CTexture
{
public:
  // constructor
  CFresnelTexture(const CValues &value, const CValues ior,
		  const CValues extinction, bool approx):
    CTexture(),
    _value(value),
    _index(ior),
    _extinction(extinction),
    _approxSchlick(approx)
    {
      assert(ior.Get(0) >= 0 && extinction.Get(0) >= 0);
    }

  CFresnelTexture(const CValues &value, const CValues ior,
		  const CValues extinction):
    CTexture(),
    _value(value),
    _index(ior),
    _extinction(extinction),
    _approxSchlick(false)
    {
      assert(ior.Get(0) >= 0 && extinction.Get(0) >= 0);
    }

  CFresnelTexture(const CValues &value, const CValues ior):
    CTexture(),
    _value(value),
    _index(ior),
    _extinction(CValues(0,0,0)),
    _approxSchlick(false)
    {
      assert(ior.Get(0) >= 0);
    }

  CFresnelTexture(const CValues &value):
    CTexture(),
    _value(value),
    _index(CValues(0,0,0)),
    _extinction(CValues(0,0,0)),
    _approxSchlick(false)
    { }

  CFresnelTexture(const CFresnelTexture &masterTexture):
    CTexture(),
    _value(masterTexture._value) { }

  // returns the type (dimensionality) of the texture (2D,3D, ...).
  virtual CMappingBase::EDimType GetType() const { return CMappingBase::EE_None; }

  // returns identification of the class
  virtual EClassId GetClassId() const { return EE_FresnelTexture; }

  // set the new color to constant texture
  void SetColor(const CValues &newColor) {
    _value = newColor;
  }

  // set the new rel.ior (real part) for constant texture
  void SetIndexR(const CValues &newIor) {
    _index = newIor;
    assert(newIor.Get(0) >= 0);
  }

  // set the new extinction coefficient for constant texture
  void SetExtinction(const CValues &newExtinction) {
    _extinction = newExtinction;
    assert(newExtinction.Get(0) >= 0);
  }

  // set the new new value for bool variable _approxSchlick
  void SetApproxSchlick(bool newApprox) {
    _approxSchlick = newApprox;
  }

  // get the color from constant texture
  void GetColor(CValues &color) const {
    color = _value;
  }

  // get the rel.ior (real part)
  void GetIndexR(CValues &ior) const {
    ior = _index;
  }

  // get the extinction coefficient
  void GetExtinction(CValues &extinction) const {
    extinction = _extinction;
  }

  // use Schlick approximation of Fresnel formula or not
  void GetApproxSchlick(bool &approx) const {
    approx = _approxSchlick;
 }

 // returns complex index of refraction (index of refraction and
 // coeffitient of extinction)
 virtual bool GetRefractionParams(CValues &ior, CValues &extinction,
                                  bool &approx, const CMapPointInfo &) const {
    if(_index.Get(0) > -1e-06) 
    {
      ior        = _index;
      extinction = _extinction;
      approx     = _approxSchlick;
      return true;
    }
    ior        = CValues(0, 0, 0);
    extinction = CValues(0, 0, 0);
    approx     = false;
    return false;
 }

 // return true if texture is of Fresnel type, otherwise false
 virtual bool IsFresnel()
 {
    if(_index.Get(0) > -1e-06)
      return true;
    return false;
 }

  // Reflectance of real metals and dielectrics (isotropic case!) depends on
  // angle between the incoming light direction and normal and is defined by
  // the Fresnel equation [2]. The maximum of albedo is achieved:
  //    for dielectics at normal incidence,
  //    for metal at grazing angles [3].
  //----------------------------------------------------------------------
  // Returns the value of albedo at normal incidence
  //----------------------------------------------------------------------
  virtual void GetValue(CValues &value, const CMapPointInfo &mapInfo) const {
    // For metals
    CValues indexa(0,0,0), extinctiona(0,0,0);
    float index, extinction;
    bool approx = false;
    bool use_fresnel = GetRefractionParams(indexa, extinctiona, approx, mapInfo);
    index       = indexa.Get(0);
    extinction  = extinctiona.Get(0);

    if(use_fresnel && (index > 1e-06 && extinction > 1e-06))
    {
      // See [4], p. 17
      float sqr_ex = extinction*extinction;
      float reflectance = ((index-1.0)*(index-1.0) + sqr_ex) /
                          ((index+1.0)*(index+1.0) + sqr_ex);

      value = reflectance * _value;
    }
    else
    {
      // For dielectrics - it is value of transmitted light
      if(index > 1e-06)
      {
        float reflectance = ((index-1.0)*(index-1.0)) /
                            ((index+1.0)*(index+1.0));
        value = (1.0-reflectance) * _value;
      }
      else // For ideal mirror
      {
        value = _value;
      }
    }
  }

  // Returns the value of albedo for specified cosinus and sinus of the angle
  // between the incoming light direction and normal and sinus of the angle
  // between the refracted ray direction and the inversed normal.
  //    (acos(angle) is time consuming operation, moreover, explicit value
  //    of the angle is not necessary. In this way it is better to use the
  //    values of cosinus and sinus directly).
  // First CMapPointInfo specifies incoming angle, second - texture mapping
  // for some particular surface
  //----------------------------------------------------------------------
  virtual void GetValue(CValues &value, const CMapPointInfo mapInfo[],
                        int size) const
  {
    assert(size > 1);

    // Actualy, incoming angle should fit into the range [0..90]. The user's
    // application is responsible for checking for "under-the-surface" ray 
    // hits. Therefore, here absolute value of cosine could be taken.
    float cosi = fabs((mapInfo[0].coord)[0]);
    float sini = fabs((mapInfo[0].coord)[1]);
    float sinr = fabs((mapInfo[0].coord)[2]);

    CValues indexa(0,0,0), extinctiona(0,0,0);
    float index      = 0, extinction = 0;
    bool approx      = false;
    bool use_fresnel = GetRefractionParams(indexa, extinctiona, approx,
                                           mapInfo[1]);
    index            = indexa.Get(0);
    extinction       = extinctiona.Get(0);

    if(use_fresnel)
    {
	  assert(cosi <= 1.0+WTINY); // cosinus of incoming angle
      assert(sini <= 1.0+WTINY); // sinus of incoming angle
     
      // Compute reflectance(theta_inc) according to Fresnel equation
      //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      if(!_approxSchlick)
      {
        // without taking into account Fresnel effect (for normal incidence)
        //------------------------------------------------------------------
        if(cosi >= (1.0 - 1e-06))
        {
	  GetValue(value, mapInfo[0]);
          assert((value[0] <= _value[0] && value[1] <= _value[1])
	          && value[2] <= _value[2]);
	  return;
        }
        else
        {
          float ind2 = index * index;
          float ext2 = extinction * extinction;

          // For metals (for unpolarized light) the reflected part of the light
	  // is computed as following (see [5], p. 734., [6], p. 26 - there is a
	  // typo in (59))
	  /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	          a^2 + b^2 - 2*a*cos(Ti) + cos(Ti)  
	     Fs = ---------------------------------
	          a^2 + b^2 + 2*a*cos(Ti) + cos(Ti)
	          
		        a^2 + b^2 - 2*a*sin(Ti)*tan(Ti) + sin^2(Ti)*tan^2(Ti)  
	     Fp = Fs * ------------------------------------------------------
	                a^2 + b^2 + 2*a*sin(Ti)*tan(Ti) + sin^2(Ti)*tan^2(Ti)  
		  
	     Fr = (Fs + Fp)/2       		     // For unpolarized light
	     
	     Ti - angle of incidence, Fs - perpendicular reflectivity
	     component, Fp - parallel reflectivity component 
	     
	     2a^2 = sqrt[(n^2 - k^2 - sin^2(Ti))^2 + 4*n^2*k^2] +
	                 (n^2 - k^2 - sin^2(Ti))
	     
	     2b^2 = sqrt[(n^2 - k^2 - sin^2(Ti))^2 + 4*n^2*k^2] -
	                 (n^2 - k^2 - sin^2(Ti))
	  
	     n, k - relative index of refraction and coeffitient of extinction
	     (n2/n1 = n + j*k       - complex index of refraction)
	  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
          if(index > 1e-06 && extinction > 1e-06)
          { 
	    if(cosi >= 1e-06)
	    {
              float sini2 = sini*sini;
              float tani = sini/cosi;
              float tani2 = tani * tani;

              float n2k2sin2 =(ind2 - ext2 - sini2);
	      float sqr_part = sqrt(pow(n2k2sin2, 2.0f) + 4.0*ind2*ext2);
	      float a2 = (sqr_part + n2k2sin2)/2.0f;
	      float b2 = (sqr_part - n2k2sin2)/2.0f;
	      float a_2 = 2.0*sqrt(a2);

	      float a2b2cos2 = a2 + b2 + cosi*cosi;
	      float a_cos_2 = a_2 * cosi;
	      float a_sin_tan_2 = a_2 * sini * tani;
	      float a2b2sin2tan2 = a2 + b2 + sini2*tani2;

	      float fs = (a2b2cos2 - a_cos_2)/(a2b2cos2 + a_cos_2);
	      float fp = fs * (a2b2sin2tan2 - a_sin_tan_2)/
	                      (a2b2sin2tan2 + a_sin_tan_2);
              float fresnel = (fs + fp)/2.0f;

	      value = fresnel * _value;
            }
	    else
	    {
              // At grazing angles limit of fresnel function goes to 1 =>
	      value = _value;
	    }
          }
          else
          {
            // For dielectrics (Compendium, p. 26)
	    /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	           [ cos(Tr)/cos(Ti) - 1/n ]^2     [ sin(Ti-Tr) ]^2
	      Fs = ---------------------------- =  ----------------
	           [ cos(Tr)/cos(Ti) + 1/n ]^2     [ sin(Ti+Tr) ]^2
	  
	           [ cos(Ti)/cos(Tr) - 1/n ]^2     [ tan(Ti-Tr) ]^2
	      Fp = ---------------------------- =  ----------------
	           [ cos(Ti)/cos(Tr) + 1/n ]^2     [ tan(Ti+Tr) ]^2
	       
	      Tr - angle of refraction
	      
	      Fr = ( Fs + Fp ) / 2
	    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
            if(index > 1e-06)
            {
	      // For the case if n_inc > n_trans and if angle is greater or
	      // equal than Brewster's angle, the total reflection happens. In
	      // this case all the light is reflected and reflectance is 
	      // maximum and constant.
	      if((sini >= 1.0/index) && (index < 1.0))
	      {
                value = _value;
	      }
	      else // Part of the light is refracted, another part is reflected
              {
                float cosr = sqrt(1.0 - sinr*sinr);
	        float ratio =  cosr / fabs(cosi);
	        float inv_index = 1.0 / index;
	        float inv_ratio = 1.0 / ratio;

	        float fs = pow(((ratio - inv_index)/(ratio + inv_index)),2.0f);
	        float fp = pow(((inv_ratio - inv_index)/
		                (inv_ratio + inv_index)), 2.0f);
                float fresnel = (fs + fp)/2.0f;
   	        value = (1.0-fresnel) * _value;
              }
	    } // if(_index > 1e-06)
            else
            {
              // For ideal mirror
	      value = _value;
            } // else(_index > 1e-06)
          }// else(_index > 1e-06 && _extinction > 1e-06)
        } // if(fabs(cosi) >= 1.0-1e-06)
        //------------------------------------------------------------------
      }
      else // Compute reflectance(theta_inc) using Schlick's approximation, [7]
      {    // See [7]; [8], p. 175

        // R(theta) = R0 + (1 - R0)*pow((1 - cos(theta)), 5)
        // Get reflectance at normal incidence -> R0
        GetValue(value, mapInfo[0]);

        // reflectance(theta) = R(theta) * coefficient
        // coefficient = _value
        // value = _value * R0 => we can rewrite reflectance formulae as:
        // reflectance(theta) = _value*(R0+(1-R0)*((1-cos(theta))^5)) =
        //                    =  value+(_value-value)*((1-cos(theta))^5)
        CValues tmp(CValues::EE_Weigth3);

        float cos_term = pow((1.0f - cosi), 5.0f);

        tmp = value + (_value - value)*cos_term;
        value = tmp;
      }
      assert((value[0] <= _value[0] && value[1] <= _value[1])
              && value[2] <= _value[2]);
    }
    else // else for if(use_fresnel)
    {
      std::cout << "Parameters are set incorrectly\n";
      exit(0);
    }
    //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  }

  CFresnelTexture & operator=(const CFresnelTexture &texture) {
    if (this != &texture)
    {
      this->_value         = texture._value;
      this->_index         = texture._index;
      this->_extinction    = texture._extinction;
      this->_approxSchlick = texture._approxSchlick;
    }
    return *this;
  }

  // Describes the texture to a given stream
  virtual void Describe(ostream &app, int ind, float multiplier = 1.0) const {
    indent(app, ind); app << "\tmax_value\t = " << _value * multiplier << "\n";
    indent(app, ind); app << "\textinction\t = " << _extinction << "\n";
    indent(app, ind); app << "\tindex of refraction\t = " << _index<< "\n";
    indent(app, ind); app << "\tapproximate by Schlick?\t -> ";
    if(_approxSchlick) app << "yes" << "\n";
    else               app << "no" << "\n";
  }

protected:
  CValues _value;          //value of this texture for each texture coordinates
  CValues _index;          // Real part of complex rel.ior
  CValues _extinction;     // Imaginary part of complex rel.ior
  bool    _approxSchlick;  // Use Schlick's aproximation or compute
                           //     reflectance directly by Fresnel's equation
};

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

__END_BRDFLIB_HEADER

#endif // __TEXTURE_H__
