/*
 * Coherent noise function over 1 or 2 dimensions.
 */
#include "dmlib.h"
#include <math.h>

#define B	(0x100)
#define BM	(0x0ff)
#define NP	(12)
#define N	(0x1000)
#define NM 	(0x0fff)


#define DM_PERLIN_SETUP(i, b0, b1, r0, r1) {	\
        t = (vec[i] + N);	\
        b0 = (((int) t) & BM);	\
        b1 = ((b0 + 1) & BM);	\
        r0 = (t - ((int) t));	\
        r1 = (r0 - 1.0f);	\
        }
        

#define DM_PERLIN_AT2(rx,ry) ((rx) * q[0] + (ry) * q[1])


static int p[B + B + 2];
static DMFloat g2[B + B + 2][2];
static DMFloat g1[B + B + 2];


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

    DM_PERLIN_SETUP(0, bx0, bx1, rx0, rx1);
    DM_PERLIN_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 = DMM_S_CURVE(rx0);
    sy = DMM_S_CURVE(ry0);

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

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

    return DMM_LERP(sy, a, b);
}


static void dmPerlinNormalize2(DMFloat v[2])
{
    DMFloat s;

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

    v[0] /= s;
    v[1] /= s;
}


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

    srand(32);
    
    for (i = 0; i < B; i++)
    {
        p[i] = i;
        g1[i] = (DMFloat) ((rand() % (B + B)) - B) / B;

        for (j = 0; j < 2; j++)
            g2[i][j] = (DMFloat) ((rand() % (B + B)) - B) / B;

        dmPerlinNormalize2(g2[i]);
    }

    while (--i)
    {
        k = p[i];
        p[i] = p[j = rand() % 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];
    }
}


/* 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.
 */
DMFloat dmPerlinNoise2D(DMFloat x, DMFloat y, DMFloat alpha, DMFloat beta, int n)
{
    int i;
    DMFloat val, sum = 0;
    DMFloat p[2], scale = 1;

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

    return (sum);
}