/* Coherent noise function over 1, 2 or 3 dimensions */
/* (copyright Ken Perlin) */

#include "config.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "perlin.h"

/* random is not portable to all platforms */
#define random g_random_int

static int p[B + B + 2];
static double g3[B + B + 2][3];
static double g2[B + B + 2][2];
static double g1[B + B + 2];
static int start = 1;

double
noise1 (double arg)
{
  int       bx0, bx1;
  double    rx0, rx1, sx, t, u, v, vec[1];

  vec[0] = arg;
  if (start)
    {
      start = 0;
      perlin_init ();
    }

  setup (0, bx0, bx1, rx0, rx1);

  sx = s_curve (rx0);
  u = rx0 * g1[p[bx0]];
  v = rx1 * g1[p[bx1]];

  return (lerp (sx, u, v));
}

double
noise2 (double vec[2])
{
  int       bx0, bx1, by0, by1, b00, b10, b01, b11;
  double    rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
  int       i, j;

  if (start)
    {
      start = 0;
      perlin_init ();
    }

  setup (0, bx0, bx1, rx0, rx1);
  setup (1, by0, by1, ry0, ry1);

  i = p[bx0];
  j = p[bx1];

  b00 = p[i + by0];
  b10 = p[j + by0];
  b01 = p[i + by1];
  b11 = p[j + by1];

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

  q = g2[b00];
  u = at2 (rx0, ry0);
  q = g2[b10];
  v = at2 (rx1, ry0);
  a = lerp (sx, u, v);

  q = g2[b01];
  u = at2 (rx0, ry1);
  q = g2[b11];
  v = at2 (rx1, ry1);
  b = lerp (sx, u, v);

  return lerp (sy, a, b);
}

double
noise3 (double vec[3])
{
  int       bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
  double    rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
  int       i, j;

  if (start)
    {
      start = 0;
      perlin_init ();
    }

  setup (0, bx0, bx1, rx0, rx1);
  setup (1, by0, by1, ry0, ry1);
  setup (2, bz0, bz1, rz0, rz1);

  i = p[bx0];
  j = p[bx1];

  b00 = p[i + by0];
  b10 = p[j + by0];
  b01 = p[i + by1];
  b11 = p[j + by1];

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

  q = g3[b00 + bz0];
  u = at3 (rx0, ry0, rz0);
  q = g3[b10 + bz0];
  v = at3 (rx1, ry0, rz0);
  a = lerp (t, u, v);

  q = g3[b01 + bz0];
  u = at3 (rx0, ry1, rz0);
  q = g3[b11 + bz0];
  v = at3 (rx1, ry1, rz0);
  b = lerp (t, u, v);

  c = lerp (sy, a, b);

  q = g3[b00 + bz1];
  u = at3 (rx0, ry0, rz1);
  q = g3[b10 + bz1];
  v = at3 (rx1, ry0, rz1);
  a = lerp (t, u, v);

  q = g3[b01 + bz1];
  u = at3 (rx0, ry1, rz1);
  q = g3[b11 + bz1];
  v = at3 (rx1, ry1, rz1);
  b = lerp (t, u, v);

  d = lerp (sy, a, b);

  return lerp (sz, c, d);
}

void
normalize2 (double v[2])
{
  double    s;

  s = sqrt (v[0] * v[0] + v[1] * v[1]);
  v[0] = v[0] / s;
  v[1] = v[1] / s;
}

void
normalize3 (double v[3])
{
  double    s;

  s = sqrt (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
  v[0] = v[0] / s;
  v[1] = v[1] / s;
  v[2] = v[2] / s;
}

void
perlin_init (void)
{
  int       i, j, k;

  for (i = 0; i < B; i++)
    {
      p[i] = i;
      g1[i] = (double) ((random () % (B + B)) - B) / B;

      for (j = 0; j < 2; j++)
        g2[i][j] = (double) ((random () % (B + B)) - B) / B;
      normalize2 (g2[i]);

      for (j = 0; j < 3; j++)
        g3[i][j] = (double) ((random () % (B + B)) - B) / B;
      normalize3 (g3[i]);
    }

  while (--i)
    {
      k = p[i];
      p[i] = p[j = random () % B];
      p[j] = k;
    }

  for (i = 0; i < B + 2; i++)
    {
      p[B + i] = p[i];
      g1[B + i] = g1[i];
      for (j = 0; j < 2; j++)
        g2[B + i][j] = g2[i][j];
      for (j = 0; j < 3; j++)
        g3[B + i][j] = g3[i][j];
    }
}

/* --- My harmonic summing functions - PDB --------------------------*/

/*
   In what follows "alpha" is the weight when the sum is formed.
   Typically it is 2, As this approaches 1 the function is noisier.
   "beta" is the harmonic scaling/spacing, typically 2.
*/

double
PerlinNoise1D (double x, double alpha, double beta, int n)
{
  int       i;
  double    val, sum = 0;
  double    p, scale = 1;

  p = x;
  for (i = 0; i < n; i++)
    {
      val = noise1 (p);
      sum += val / scale;
      scale *= alpha;
      p *= beta;
    }
  return (sum);
}

double
PerlinNoise2D (double x, double y, double alpha, double beta, int n)
{
  int       i;
  double    val, sum = 0;
  double    p[2], scale = 1;

  p[0] = x;
  p[1] = y;
  for (i = 0; i < n; i++)
    {
      val = noise2 (p);
      sum += val / scale;
      scale *= alpha;
      p[0] *= beta;
      p[1] *= beta;
    }
  return (sum);
}

double
PerlinNoise3D (double x, double y, double z, double alpha, double beta, int n)
{
  int       i;
  double    val, sum = 0;
  double    p[3], scale = 1.0;

  if (z < 0.0000)
    return PerlinNoise2D (x, y, alpha, beta, n);
  p[0] = x;
  p[1] = y;
  p[2] = z;

  for (i = 0; i < n; i++)
    {
      val = noise3 (p);
      sum += val / scale;
      scale *= alpha;
      p[0] *= beta;
      p[1] *= beta;
      p[2] *= beta;
    }
  return (sum);
}