// ===================================================================
// $Id:$
//
// rndgauss.h:
//     A set of functions to generate random numbers according to
//     bivariate normal (=Gaussian) distribution + special distribution
//     on disk combining it with bivariate normal distribution
//
// Initial coding by Vlastimil Havran, January 2007

#ifndef __RNDGAUSS_H__
#define __RNDGAUSS_H__

//#include "configh.h"
//#include "basrnd.h"
//#include "errorg.h"
#include <assert.h>
#include "Vector2.h"
//#include "sphcoord.h"

namespace GtpVisibilityPreprocessor {

// Function: GaussianOn2D
// Here we generate the Gaussian distribution on the 2D domain
// (=bivariate normal distribution) given the uniformly distributed random
// numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
// We specify also standard deviation (sigma) of the normal distribution.
// The result has a mean value <0.0f, 0.0f>. We use Box-Muller transform
// in polar form.
// Please note that in the distance from mean value this is true:
// the distance = 1 x sigma ... we have 68.25% of generated points
//              = 2 x sigma ... we have 95.45% of generated points
//              = 3 x sigma ... we have 99.73% of generated points
//              = 4 x sigma ... we have 99.99% of generated points

inline void
GaussianOn2D(const Vector2 &input, // input (RND in interval [0-1)x[0-1))
			 float sigma,            // standard deviation of gaussian distribution
			 Vector2 &outputVec)   // resulting vector according to gaussian distr.
{
  // Handling zero
  float xf = input.xx;
  if (xf < 1e-20f)
    xf = 1e-20f;
  assert( xf <= 1.0f);
  assert( (input.yy >= 0.f) && (input.yy <= 1.0f));

  // Here we have to use natural logarithm function !
  float t = -2.0f * log(xf);
  assert(t > 0.f);
  t = sigma * sqrt(t);
  float phi = 2.0f * M_PI * input.yy;
  outputVec.xx = t * cos(phi);
  outputVec.yy = t * sin(phi);

  return;
}

// Function: PolarGaussianOnDisk
// Here we generate the Gaussian distribution on the 2D domain on the part
// of 2D disk. The input of the function is tuple of random
// numbers in the range [0-1)x[0-1) (bivariate uniform distribution).
// The result of this function 'outputVec' (in general in the range
// [-infinity,+infinity]x[-infinity,+infinity], so the point generated
// according to proposed distribution. The algorithm does not use rejection
// sampling. We have to specify the radius of 2D circle
// where the distribution mean value is zero for one variable !!!!
// Sigma is specified independently, but if you want to have small probability
// of generating the original point, then consider:
// the radius = 1 x sigma ... 68.25% points in the distance from sigma from the circle
//            = 2 x sigma ... 95.45% points in the distance from sigma from the circle
//            = 3 x sigma ... 99.73% points in the distance from sigma from the circle
//            = 4 x sigma ... 99.99% points in the distance from sigma from the circle
// Recommendation: use sigma = radius/3.3f or so (according to visualization)

// This is only helper function
inline float
ComputeSigmaFromRadius(float radius)
{
  // For standard deviation = 3 we have achieve probability 0.4% that the proximity
  // of the center of the circle will be covered again by a generated point.
  return radius / 3.3f;
}

inline void
PolarGaussianOnDisk(const Vector3 &input, // input (RND in interval [0-1)x[0-1))
					float sigma,    // standard deviation of gaussian distribution
					float radius,   // standard deviation of gaussian distribution
					Vector2 &outputVec) // result
{
  // Handling zero
  float xf = input.x;
  if (xf < 1e-20f)
    xf = 1e-20f;
  assert( xf <= 1.0f);
  assert( (input.y >= 0.f) && (input.y <= 1.0f));

  // Here we have to use natural logarithm function !
  float t = -2.0f * log(xf);
  assert(t > 0.f);
  t = sigma * sqrt(t);
  float phi = 2.0f * M_PI * input.y;

#if 1
  // This is what Jiri Bittner wanted originally,
  // another idea how to implement this distribution
  // First, we generate point as for standard polar
  // inmplementation
  outputVec.xx = t * cos(phi);
  outputVec.yy = t * sin(phi);

  // Then we generate a point on 2D circle at the distance
  // of radius
  float r = input.z; //RND::rand01f();
  float phi2 = 2.0f * M_PI * r;
  Vector2 c(cos(phi2),sin(phi2));
  outputVec.xx += radius * c.xx;
  outputVec.yy += radius * c.yy;
#endif
#if 0
  // This is what Jiri Bittner wanted originally
  // however, I cannot make it running correctly !!!
  // Do not use it at the moment !!! VH 11/1/2007
  float r = RND::rand01f();  
  if (r < 0.5f) {
    outputVec.xx = (radius + t) * cos(phi);
    outputVec.yy = (radius + t) * sin(phi);
  }
  else {
    // going inside the cirle
    outputVec.xx = (radius - t) * cos(phi);
    outputVec.yy = (radius - t) * sin(phi);    
  }
#endif
  
  return;
}

}

#endif // __RND_H__