//
// Noise.cc
//
// original implementation due to Ken Perlin, via Ken Musgrave
//
// $Id: noise.cc,v 1.1.1.1 2001/02/28 00:28:37 cstolte Exp $
//

#include <sgl/noise.h>
#include <sgl/mathfunc.h>
#include <sgl/vecmath3.h>
#include <assert.h>

extern int p_precomputed[];
extern double g_precomputed[][3];
const u_int hash_size = 0x100;

//#define FILTER_ARGS

// utility functions

static void
setup(double pos, int &hash0, int &hash1, double &r0, double &r1)
{
  const int big_random_constant = 0x100000;

  double tmp = pos + big_random_constant;
  int int_tmp = int(tmp);

  hash0 = int_tmp & (hash_size-1);
  hash1 = (hash0+1) & (hash_size-1);
  r0 = tmp-int_tmp;
  r1 = r0-1.;
}

static double 
at2(double rx, double ry, double *q)
{
  return rx*q[0]+ry*q[1];
}

static double 
s_curve(double val)
{
  return val*val*(3-2*val);
}

static double 
at3(double rx, double ry, double rz, double *q)
{
  return rx*q[0]+ry*q[1]+rz*q[2];
}

#ifdef FILTER_ARGS
static void 
filter_arg(double *arg)
{
  assert(fabs(*arg) < 1.e13);
  const u_int maxval = 2147483647;
  const u_int two_maxval = 4294967294;

  if (*arg < 0) *arg = -*arg;

  while (*arg > maxval || *arg < 0)
    {
      if (*arg > maxval)
	  *arg = two_maxval - *arg;
      else
	  *arg = -*arg;
    }
}
#endif /* FILTER_ARGS */

// implementations of noise functions

double
Noise(double x)
{
    return Noise2(x, 0.5525221);
}

double 
Noise2(double x, double y)
{
  int bx0,bx1,by0,by1;
  double rx0,rx1,ry0,ry1;

#ifdef FILTER_ARGS
  filter_arg(&x);
  filter_arg(&y);
#endif /* FILTER_ARGS */

  setup(x, bx0, bx1, rx0, rx1);
  setup(y, by0, by1, ry0, ry1);

  int hashLL = p_precomputed[p_precomputed[bx0]+by0];
  int hashLR = p_precomputed[p_precomputed[bx1]+by0];
  int hashUL = p_precomputed[p_precomputed[bx0]+by1];
  int hashUR = p_precomputed[p_precomputed[bx1]+by1];

  double weightLL = at2(rx0, ry0, g_precomputed[hashLL]);
  double weightLR = at2(rx1, ry0, g_precomputed[hashLR]);
  double weightUL = at2(rx0, ry1, g_precomputed[hashUL]);
  double weightUR = at2(rx1, ry1, g_precomputed[hashUR]);

  double sx = s_curve(rx0);
  double lower_weight = lerp(sx, weightLL, weightLR);
  double upper_weight = lerp(sx, weightUL, weightUR);
  double sy = s_curve(ry0);
  double alt = lerp(sy, lower_weight, upper_weight);

  return alt;
}

double 
DNoise2(double x, double y, double *dx,double *dy)
{
  int bx0, bx1, by0, by1;
  double rx0, rx1, ry0, ry1;

#ifdef FILTER_ARGS
  filter_arg(&x);
  filter_arg(&y);
#endif /* FILTER_ARGS */

  setup(x, bx0, bx1, rx0, rx1);
  setup(y, by0, by1, ry0, ry1);

  int hashLL = p_precomputed[p_precomputed[bx0]+by0];
  int hashLR = p_precomputed[p_precomputed[bx1]+by0];
  int hashUL = p_precomputed[p_precomputed[bx0]+by1];
  int hashUR = p_precomputed[p_precomputed[bx1]+by1];

  double weightLL = at2(rx0, ry0, g_precomputed[hashLL]);
  double weightLR = at2(rx1, ry0, g_precomputed[hashLR]);
  double weightUL = at2(rx0, ry1, g_precomputed[hashUL]);
  double weightUR = at2(rx1, ry1, g_precomputed[hashUR]);

  double sx = s_curve(rx0);
  double lower_weight = lerp(sx, weightLL, weightLR);
  double upper_weight = lerp(sx, weightUL, weightUR);
  double sy = s_curve(ry0);
  double alt = lerp(sy, lower_weight, upper_weight);

  const double delta = 0.001;
  const double delta_inv = 1000.;

  double sxd = s_curve(rx0+delta);
  double syd = s_curve(ry0+delta);
  double dx0 = lerp(sxd, weightLL, weightLR) - lower_weight;

  double dx1= lerp(sxd, weightUL, weightUR) - upper_weight;
  *dx = lerp(sy, dx0, dx1) * delta_inv;
  *dy = (lerp(syd, lower_weight, upper_weight) - alt) * delta_inv;

  return alt;
}

double 
Noise3(double x, double y, double z)
{
  int bx0, bx1, by0, by1, bz0, bz1;
  double rx0, rx1, ry0, ry1, rz0, rz1;
  double a,b,c,d;

#ifdef FILTER_ARGS
  filter_arg(&x);
  filter_arg(&y);
  filter_arg(&z);
#endif /* FILTER_ARGS */

  setup(x,bx0,bx1,rx0,rx1);
  setup(y,by0,by1,ry0,ry1);
  setup(z,bz0,bz1,rz0,rz1);

  int hashLL = p_precomputed[p_precomputed[bx0]+by0];
  int hashLR = p_precomputed[p_precomputed[bx1]+by0];
  int hashUL = p_precomputed[p_precomputed[bx0]+by1];
  int hashUR = p_precomputed[p_precomputed[bx1]+by1];

  double weightBLL = at3(rx0, ry0, rz0, g_precomputed[hashLL+bz0]);
  double weightBLR = at3(rx1, ry0, rz0, g_precomputed[hashLR+bz0]);
  double weightBUL = at3(rx0, ry1, rz0, g_precomputed[hashUL+bz0]);
  double weightBUR = at3(rx1, ry1, rz0, g_precomputed[hashUR+bz0]);
  double weightTLL = at3(rx0, ry0, rz1, g_precomputed[hashLL+bz1]);
  double weightTLR = at3(rx1, ry0, rz1, g_precomputed[hashLR+bz1]);
  double weightTUL = at3(rx0, ry1, rz1, g_precomputed[hashUL+bz1]);
  double weightTUR = at3(rx1, ry1, rz1, g_precomputed[hashUR+bz1]);

  double sx = s_curve(rx0);
  double sy = s_curve(ry0);
  double sz = s_curve(rz0);

  a = lerp(sx, weightBLL, weightBLR);
  b = lerp(sx, weightBUL, weightBUR);
  c = lerp(sy, a, b);
  a = lerp(sx, weightTLL, weightTLR);
  b = lerp(sx, weightTUL, weightTUR);
  d = lerp(sy, a, b);

  return lerp(sz,c,d);
}

double 
Noise3(const Tuple3d &p)
{
    return Noise3(p.x, p.y, p.z);
}

Vector3
DNoise3(double x, double y, double z)
{
  int bx0, bx1, by0, by1, bz0, bz1;
  double rx0, rx1, ry0, ry1, rz0, rz1;
  double sx, sy, sz;
  double a0, a1, a2, b0, b1, b2, c0, c1, c2, d0, d1, d2;
  double u0, u1, u2, v0, v1, v2;

#ifdef FILTER_ARGS
  filter_arg(&x);
  filter_arg(&y);
  filter_arg(&z);
#endif // FILTER_ARGS

  setup(x, bx0, bx1, rx0, rx1);
  setup(y, by0, by1, ry0, ry1);
  setup(z, bz0, bz1, rz0, rz1);

  int hashLL = p_precomputed[p_precomputed[bx0]+by0];
  int hashLR = p_precomputed[p_precomputed[bx1]+by0];
  int hashUL = p_precomputed[p_precomputed[bx0]+by1];
  int hashUR = p_precomputed[p_precomputed[bx1]+by1];

  sx = s_curve(rx0);
  sy = s_curve(ry0);
  sz = s_curve(rz0);

  u0 = at3(rx0, ry0, rz0, g_precomputed[hashLL+bz0]);
  u1 = at3(rx0, ry0, rz0, g_precomputed[hashLL+bz0+1]);
  u2 = at3(rx0, ry0, rz0, g_precomputed[hashLL+bz0+2]);
  v0 = at3(rx1, ry0, rz0, g_precomputed[hashLR+bz0]);
  v1 = at3(rx1, ry0, rz0, g_precomputed[hashLR+bz0+1]);
  v2 = at3(rx1, ry0, rz0, g_precomputed[hashLR+bz0+2]);
  a0 = lerp(sx, u0, v0);
  a1 = lerp(sx, u1, v1);
  a2 = lerp(sx, u2, v2);

  u0 = at3(rx0, ry1, rz0, g_precomputed[hashUL+bz0]);
  u1 = at3(rx0, ry1, rz0, g_precomputed[hashUL+bz0+1]);
  u2 = at3(rx0, ry1, rz0, g_precomputed[hashUL+bz0+2]);
  v0 = at3(rx1, ry1, rz0, g_precomputed[hashUR+bz0]);
  v1 = at3(rx1, ry1, rz0, g_precomputed[hashUR+bz0+1]);
  v2 = at3(rx1, ry1, rz0, g_precomputed[hashUR+bz0+2]);
  b0 = lerp(sx, u0, v0);
  b1 = lerp(sx, u1, v1);
  b2 = lerp(sx, u2, v2);

  c0= lerp(sy, a0, b0);
  c1= lerp(sy, a1, b1);
  c2= lerp(sy, a2, b2);

  u0 = at3(rx0, ry0, rz1, g_precomputed[hashLL+bz1]);
  u1 = at3(rx0, ry0, rz1, g_precomputed[hashLL+bz1+1]);
  u2 = at3(rx0, ry0, rz1, g_precomputed[hashLL+bz1+2]);
  v0 = at3(rx1, ry0, rz1, g_precomputed[hashLR+bz1]);
  v1 = at3(rx1, ry0, rz1, g_precomputed[hashLR+bz1+1]);
  v2 = at3(rx1, ry0, rz1, g_precomputed[hashLR+bz1+2]);
  a0 = lerp(sx, u0, v0);
  a1 = lerp(sx, u1, v1);
  a2 = lerp(sx, u2, v2);

  u0 = at3(rx0,ry1,rz1,g_precomputed[hashUL+bz1]);
  u1 = at3(rx0,ry1,rz1,g_precomputed[hashUL+bz1+1]);
  u2 = at3(rx0,ry1,rz1,g_precomputed[hashUL+bz1+2]);
  v0 = at3(rx1,ry1,rz1,g_precomputed[hashUR+bz1]);
  v1 = at3(rx1,ry1,rz1,g_precomputed[hashUR+bz1+1]);
  v2 = at3(rx1,ry1,rz1,g_precomputed[hashUR+bz1+2]);
  b0 = lerp(sx,u0,v0);
  b1 = lerp(sx,u1,v1);
  b2 = lerp(sx,u2,v2);

  d0 = lerp(sy,a0,b0);
  d1 = lerp(sy,a1,b1);
  d2 = lerp(sy,a2,b2);

  return Vector3(c0 + sz * (d0 - c0),
		 c1 + sz * (d1 - c1),
		 c2 + sz * (d2 - c2));
#if 0
  // why oh why does g++ give errors with this? cpp's output
  // checks out just fine--it's same as the above fix, just with
  // more (all balanced, etc) parenthesis.. --mmp

  return Tuple3d((lerp(sz,c0,d0)),
		 (lerp(sz,c1,d1)),
		 (lerp(sz,c2,d2)));
#endif
}

Vector3
DNoise3(const Tuple3d &p)
{
    return DNoise3(p.x, p.y, p.z);
}

// functions built on top of noise

double 
FBm(double x, double y, double z, double omega, double lambda, double octaves)
{
  double omega_prime = 1;
  double ret = 0;

  for (int i = 0; i < int(octaves); ++i)
    {
      ret += omega_prime * Noise3(x, y, z);
      omega_prime *= omega;
      x *= lambda;
      y *= lambda;
      z *= lambda;
    }
  ret += (octaves-int(octaves)) * omega_prime * Noise3(x, y, z);
  return ret;
}

double
FBm(const Tuple3d &p, double omega, double lambda, double octaves)
{
    return FBm(p.x, p.y, p.z, omega, lambda, octaves);
}

double 
Chaos(double x, double y, double z, double octaves)
{
  return FBm(x, y, z, .5, 2., octaves);
}

double 
Chaos(const Tuple3d &p, double octaves)
{
  return FBm(p.x, p.y, p.z, .5, 2., octaves);
}

Vector3
VfBm(double x, double y, double z, double omega, double lambda, double octaves)
{
  double omega_prime = 1;
  Vector3 ret(0, 0, 0);