/* 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);

}