// ===================================================================
// $Id: rnd.h,v 1.3 2003/10/30 16:36:35 havran Exp $
//
// rnd.h:
//     A set of functions to generate random numbers according to
//     various distributions.
//
// Class: RND
//
// 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 Jaroslav Krivanek (oct 2000) under
// supervision of Vlastimil Havran.
//
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


#ifndef __RND_H__
#define __RND_H__

// BRDFLIB headers
#include "configh.h"
#include "errorg.h"
#include "vector2.h"
#include "sphcoord.h"


__BEGIN_BRDFLIB_HEADER

const float RAND_MAX_FLT = (float)RAND_MAX;

const float INV_RAND_MAX_FLT = (float) (1.0 / ((double)RAND_MAX+1.0));

//const double INV_RAND_MAX_DOUBLE = 1.0 / ((double)RAND_MAX+1.0);

const float TWO_TIMES_INV_RAND_MAX_FLT = (2.0 / ((double)RAND_MAX_FLT+1.0));

class RND
{
public:
  /**
     Returns a pseudo-random integer between 0 and RAND_MAX, int.
     Range is [0; RAND_MAX]
  */
  static inline int randi(void)
  { return rand(); }

  /**
     Returns a pseudo-random float between 0 and 1. Range is [0;1), float.
  */
  static inline float rand01f(void)
  {
    return RandomValue();
  }

  /**
     Returns a pseudo-random float between 0 and 1. Range is [0;1), double.
  */
  static inline double rand01d(void)
  {
    return (double) RandomValue();
  }
  
  /**
     Returns a pseudo-random float between -1 and 1. Range is [-1;1), float.
  */
  static inline float rand11f(void)
  {
    // return -1.0f + (float)rand() * TWO_TIMES_INV_RAND_MAX_FLT ;
    return -1.0f + (float)RandomValue() * 2.0;
  }
  
  /**
     Returns two pseudo-random floats between 0 and 1.
     Range is [0;1)x[0;1), float.
  */
  static inline CVector2D rand01v(void)
  {
    return CVector2D(rand01f(), rand01f());
  }   

  /**
     Returns two pseudo-random floats between -1 and 1.
     Range is [-1;1)x[-1;1), float.
  */
  static inline CVector2D rand11v(void)
  {
    return CVector2D(rand11f(), rand11f());
  }   
  
  /**
     Returns a pseudo-random float between 0 and b.
     Range is [0;b), float
  */
  static inline float rand0bf(float b)
  {
    return RandomValue() * b ;
  }
     
  /**
     Returns a vector of two pseudo-random floats between 0 and b.
     Range is [0;b) x [0;b), float.
  */
  static inline CVector2D rand0bv(float b)
  {
    return CVector2D(rand0bf(b), rand0bf(b));
  }
     
  /**
     Returns a pseudo-random float between 'a' and 'b'.
     Range is [a;b) x [a;b), float.
  */
  static inline float randabf(float a, float b)
  {
    float range = fabs(a - b);
    return RandomValue() * range + ((a < b) ? a : b);
  }

  /**
     Generate random point on unit circle with density p(x) = 1/PI.

     @return   the generated point as CVector2D
  */
  static inline CVector2D disc_const(void)
  {
    CVector2D v;
    do {v = rand11v();}
    while(v.xx * v.xx + v.yy * v.yy > 1.0f);
    return v;
  }

  /**
     Generate random point in unit disk in 2D using concentric map
     with density p(x) = 1/A.

     @param v  two random numbers between 0 and 1.
     @return   the generated point as CVector2D in 2D spherical coordinates
               (R,FI), radius R in <0,1>, angle FI in <0,2*PI>
  */
  static inline CVector2D
  disk_const(const CVector2D & w = rand01v())
  {
    assert(w.xx>=0.0);
    assert(w.xx<=1.0);
    assert(w.yy>=0.0);
    assert(w.yy<=1.0);

    float vx = -1.0 + 2.0*w.xx;
    float vy = -1.0 + 2.0*w.yy;
    float r, phi;

    if (vx > -vy) {
      if (vx > vy) {
	r = vx; phi = M_PI_4 * (vy / vx);
	// This modification is required to get range <0;2*M_PI>
	if (phi<0.0)
	  phi += 2.0f * M_PI;
      }
      else {
	r = vy; phi = M_PI_4 * (2.0 - (vx / vy) );
      }
    }
    else {
      // It holds: vx <= -vy
      if (vx < vy) {
	r = -vx; phi = M_PI_4 * (4.0 + (vy/vx) );
      }
      else {
	r = -vy;
	if (vy != 0.0)
	  phi = M_PI_4 * (6.0 - (vx/vy) );
	else
	  phi = 0.0;
      }
    }
    
    assert(r<=1.00001);
    assert(phi>=-0.000001);
    assert(phi<= 2.0*M_PI+0.000001);    
    return CVector2D(r, phi);
  }
  
  /**
     Generate random point on unit sphere.

     @param  v(in)        two random numbers between 0 and 1
     @return the generated point as CVector3D
  */
  static inline CVector3D sph_const(const CVector2D &v = rand01v())
  {
    float tmp = 2.0 * sqrt(v.xx * (1.0 - v.xx));
    float tmp2 = 2.0 * M_PI * v.yy;
    return CVector3D(cos(tmp2) * tmp, sin(tmp2) * tmp, 1.0 - 2.0 * v.xx);
  }

  /**
     Generate random direction on unit hemisphere with density
     p(theta) = 1/2*PI .

     @param  result(out)  resulting direction
     @param  v(in)        two random numbers between 0 and 1

     @return the probability density at generated direction
             (which is always 1/2*PI for this method)
  */
  static inline float
  hemisph_const(CVector3D &result, const CVector2D &v = rand01v())
  {
    float tmp = sqrt(1.0 - v.xx * v.xx);
    float tmp2 =  2.0 * M_PI * v.yy;
    result.SetValue( cos( tmp2 ) * tmp, sin( tmp2 ) * tmp, v.xx );
    return M_PI_2;
  }

  /**
     The same as above in normalized spherical coordinates.
  */
  static inline float
  hemisph_const(CSphCoordNorm &result,const CVector2D &v = rand01v())
  {
    result.SetValue(acos(v.xx), 2.0 * M_PI * v.yy);
    return M_PI_2;
  }

  /**
     Generate random direction on unit hemisphere with uniform density.

     @param result(out) resulting direction
     @param v(in)       two random numbers discribing the radius and angle
                        in two dimensions, assuming uniform density on a disk
			v.xx .. radius in <0,1>, v.yy .. angle in <0,2*PI>

     @return the probability density at generated direction
  */
  static inline float
  hemisph_from_disk(CVector3D &result, const CVector2D &v)
  {
    float tmp2 = 1.0 - v.xx * v.xx;
    float tmp3 = v.xx * sqrt(1.0 + tmp2);
    result.SetValue( cos(v.yy) * tmp3, sin(v.yy) * tmp3, tmp2);
    return M_PI_2;
  }
  
  /**
     Generate random direction on hemisphere with theta<theta_max with density
     p(theta) = 1/(2*PI*(1-cos(theta_max)) .

     @param  result(out)  resulting direction
     @param  cos_t_max    cosinus of maximaml theta
     @param  v(in)        two random numbers between 0 and 1

     @return the probability density at generated direction
             (which is always 1/(2*PI*(1-cos(theta_max)) for this method)
  */
  static inline float
  hemisph_const_clamp(CVector3D &result, float cos_t_max,
		      const CVector2D &v = rand01v())
  {
    float tmp = sqrt(1 - sqr (1 - v.xx * (1.0 - cos_t_max) )  );
    float tmp2 =  2.0 * M_PI * v.yy;
    result.SetValue( cos(tmp2) * tmp, sin(tmp2) * tmp,
		     1.0 - v.xx * (1.0 - cos_t_max) );
    return M_PI_2/(1.0 - cos_t_max);
  }

  /**
     Generate random direction on unit hemisphere with density
     p(theta) = cos(theta) / PI .

     @param result(out) resulting direction
     @param v(in)       two random numbers between 0 and 1

     @return the probability density at generated direction
  */
  static inline float
  hemisph_cos(CVector3D &result, const CVector2D &v = rand01v())
  {
    float tmp = sqrt(1.0 - v.xx);
    float tmp2 = 2.0*M_PI*v.yy;
    result.SetValue( cos(tmp2) * tmp, sin(tmp2) * tmp, sqrt(v.xx) );
    return sqrt(v.xx) / M_PI;
  }

  /**
     The same as above in normalized spherical coordinates.
  */
  static inline float
  hemisph_cos(CSphCoordNorm &result, const CVector2D &v = rand01v())
  {
    float sq = sqrt(v.xx);
    result.SetValue(acos(sq), 2.0 * M_PI * v.yy);
    return sq/M_PI;
  }
  
  /**
     Generate random direction on unit hemisphere with density
     p(theta) = (n+1)/(2*PI) * cos^n(theta) .

     @param result(out) resulting direction
     @param n(in)       exponent in density
     @param v(in)       two random numbers between 0 and 1

     @return the probability density at generated direction
  */
  static inline float
  hemisph_cosn(CVector3D &result, float n, const CVector2D &v = rand01v())
  {
    float tmp2 = pow((float)(v.xx), (float)(1.0/(n + 1.0)));
    float tmp3 = sqrt(1.0 - tmp2 * tmp2);
    float tmp4 =  2.0*M_PI*v.yy;
    result.SetValue( cos(tmp4) * tmp3, sin(tmp4) * tmp3, tmp2);
    return pow(tmp2, n) * (n + 1.0)/(2.0*M_PI);
  }
};


__END_BRDFLIB_HEADER

#endif // __RND_H__

