/*
 * DMLib
 * -- Perlin noise functionality
 * Programmed and designed by Matti 'ccr' Hamalainen
 * (C) Copyright 2012 Tecnic Software productions (TNSP)
 */
#include "dmperlin.h"
#include "dmcurves.h"
#include <math.h>


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

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


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

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

    i = ctx->p[bx0];
    j = ctx->p[bx1];

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

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

    q = ctx->g2[b00];
    u = DM_PERLIN_AT2(rx0, ry0);

    q = ctx->g2[b10];
    v = DM_PERLIN_AT2(rx1, ry0);
    a = DMM_LERP(sx, u, v);

    q = ctx->g2[b01];
    u = DM_PERLIN_AT2(rx0, ry1);

    q = ctx->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 = sqrt(v[0] * v[0] + v[1] * v[1]);
    v[0] /= s;
    v[1] /= s;
}


int dmPerlinInit(DMPerlinContext *ctx, const int seed)
{
    int i, j;

    if (ctx == NULL)
        return DMERR_NULLPTR;

    memset(ctx, 0, sizeof(*ctx));

    if (seed < 0)
        return DMERR_INVALID_ARGS;

    srand(seed);

    for (i = 0; i < DM_PLNS_B; i++)
    {
        ctx->p[i] = i;
        ctx->g1[i] = (DMFloat) ((rand() % (DM_PLNS_B + DM_PLNS_B)) - DM_PLNS_B) / DM_PLNS_B;

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

        dmPerlinNormalize2(ctx->g2[i]);
    }

    while (--i)
    {
        int k = ctx->p[i];
        ctx->p[i] = ctx->p[j = rand() % DM_PLNS_B];
        ctx->p[j] = k;
    }

    for (i = 0; i < DM_PLNS_B + 2; i++)
    {
        ctx->p[DM_PLNS_B + i] = ctx->p[i];
        ctx->g1[DM_PLNS_B + i] = ctx->g1[i];

        for (j = 0; j < 2; j++)
            ctx->g2[DM_PLNS_B + i][j] = ctx->g2[i][j];
    }

    return DMERR_OK;
}


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

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

    return sum;
}