#include <math.h>
#include <stdlib.h>
#include <string.h>

#include "log.h"
#include "rng.h"

#include "perlin.h"

#define NOISE_MAX_OCTAVES         128
#define NOISE_MAX_DIMENSIONS      4
#define NOISE_DEFAULT_HURST       0.5
#define NOISE_DEFAULT_LACUNARITY  2.

#define LERP(a, b, x)   (a + x * (b - a))
#define ABS(a)          ((a)<0?-(a):(a))
#define CLAMP(a, b, x)  ((x) < (a) ? (a) : ((x) > (b) ? (b) : (x)))

typedef void* noise_t;

/* Used internally. */
typedef struct {
  int ndim;
  unsigned char map[256]; /* Randomized map of indexes into buffer. */
  float buffer[256][NOISE_MAX_DIMENSIONS];  // Random 256 x ndim buffer. */
  /* Fractal stuff. */
  float H;
  float lacunarity;
  float exponent[NOISE_MAX_OCTAVES];
} perlin_data_t;

static float* noise_genNebulae(const int w, const int h, const int n, float rug);
static noise_t noise_new(int dimensions, float hurst, float lacunarity);
/* Basic perlin noise. */
static float noise_get(noise_t noise, float* f);
/* Fractional brownian motion. */
/*static float noise_fbm(noise_t noise, float* f, float octaves);*/
/* Turbulence. */
static float noise_turbulence(noise_t noise, float* f, float octaves);
static void noise_delete(noise_t noise);

static float lattice(perlin_data_t* data, int ix, float fx, int iy,
    float fy, int iz, float fz, int iw, float fw) {

  int   n[4] = { ix, iy, iz, iw };
  float f[4] = { fx, fy, fz, fw };
  int nindex = 0;
  int i;
  float value = 0;

  for(i = 0; i < data->ndim; i++)
    nindex = data->map[(nindex + n[i]) & 0xFF];
  for(i = 0; i < data->ndim; i++)
    value += data->buffer[nindex][i] * f[i];

  return value;
}

#define DEFAULT_SEED    0x15687436
#define DELTA           1e-6f
#define SWAP(a, b, t)   t = a; a = b; b = t

#define FLOOR(a)        ((int) a - (a < 0 && a != (int)a))
#define CUBIC(a)        (a * a * (3 - 2 * a))

static void normalize(perlin_data_t* data, float* f) {
  float magnitude = 0;
  int i;
  for(i = 0; i < data->ndim; i++)
    magnitude += f[i] * f[i];
  magnitude = 1 / sqrtf(magnitude);
  for(i = 0; i < data->ndim; i++)
    f[i] *= magnitude;
}

static noise_t noise_new(int ndim, float hurst, float lacunarity) {
  perlin_data_t* data=(perlin_data_t*)calloc(sizeof(perlin_data_t), 1);
  int i, j;
  unsigned char tmp;
  float f = 1;
  data->ndim = ndim;
  for(i = 0; i < 256; i++) {
    data->map[i] = (unsigned char) i;
    for(j = 0; j < data->ndim; j++)
      data->buffer[i][j] = RNGF()-0.5;
    normalize(data, data->buffer[i]);
  }

  while(--i) {
    j = RNG(0, 255);
    SWAP(data->map[i], data->map[j], tmp);
  }

  data->H = hurst;
  data->lacunarity = lacunarity;
  for(i = 0; i < NOISE_MAX_OCTAVES; i++) {
    /*exponent[i] = powf(f, -H); */
    data->exponent[i] = 1.0f / f;
    f *= lacunarity;
  }
  return (noise_t)data;
}

static float noise_get(noise_t noise, float *f )
{
  perlin_data_t* data = (perlin_data_t*) noise;
  int n[NOISE_MAX_DIMENSIONS];     /* Indexes to pass to lattice function */
  int i;
  float r[NOISE_MAX_DIMENSIONS];   /* Remainders to pass to lattice function */
  float w[NOISE_MAX_DIMENSIONS];   /* Cubic values to pass to interpolation function */
  float value;

  for(i=0; i<data->ndim; i++) {
    n[i] = FLOOR(f[i]);
    r[i] = f[i] - n[i];
    w[i] = CUBIC(r[i]);
  }

  switch(data->ndim) {
    case 1:
      value = LERP(lattice(data,n[0], r[0],0,0,0,0,0,0),
          lattice(data,n[0]+1, r[0]-1,0,0,0,0,0,0),
          w[0]);
      break;
    case 2:
      value = LERP(LERP(lattice(data,n[0], r[0], n[1], r[1],0,0,0,0),
                 lattice(data,n[0]+1, r[0]-1, n[1], r[1],0,0,0,0),
                 w[0]),
              LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1,0,0,0,0),
                 lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1,0,0,0,0),
                 w[0]),
              w[1]);
      break;
    case 3:
      value = LERP(LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2], r[2],0,0),
                  lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2], r[2],0,0),
                  w[0]),
                 LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2], r[2],0,0),
                  lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2], r[2],0,0),
                  w[0]),
                 w[1]),
              LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2]+1, r[2]-1,0,0),
                  lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2]+1, r[2]-1,0,0),
                  w[0]),
                 LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2]+1, r[2]-1,0,0),
                  lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1,0,0),
                  w[0]),
                 w[1]),
              w[2]);
      break;
    case 4:
    default:
      value = LERP(LERP(LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2], r[2], n[3], r[3]),
                     lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2], r[2], n[3], r[3]),
                     w[0]),
                  LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2], r[2], n[3], r[3]),
                     lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2], r[2], n[3], r[3]),
                     w[0]),
                  w[1]),
                  LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2]+1, r[2]-1, n[3], r[3]),
                     lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2]+1, r[2]-1, n[3], r[3]),
                     w[0]),
                  LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2]+1, r[2]-1,0,0),
                     lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1, n[3], r[3]),
                     w[0]),
                  w[1]),
                 w[2]),
              LERP(LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2], r[2], n[3]+1, r[3]-1),
                     lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2], r[2], n[3]+1, r[3]-1),
                     w[0]),
                  LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2], r[2], n[3]+1, r[3]-1),
                     lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2], r[2], n[3]+1, r[3]-1),
                     w[0]),
                  w[1]),
                  LERP(LERP(lattice(data,n[0], r[0], n[1], r[1], n[2]+1, r[2]-1, n[3]+1, r[3]-1),
                     lattice(data,n[0]+1, r[0]-1, n[1], r[1], n[2]+1, r[2]-1, n[3]+1, r[3]-1),
                     w[0]),
                  LERP(lattice(data,n[0], r[0], n[1]+1, r[1]-1, n[2]+1, r[2]-1,0,0),
                     lattice(data,n[0]+1, r[0]-1, n[1]+1, r[1]-1, n[2]+1, r[2]-1, n[3]+1, r[3]-1),
                     w[0]),
                  w[1]),
                 w[2]),
              w[3]);
      break;
  }
  return CLAMP(-0.99999f, 0.99999f, value);
}

#if 0
float noise_fbm(noise_t noise, float* f, float octaves) {
  float tf[NOISE_MAX_DIMENSIONS];
  perlin_data_t* data = (perlin_data_t*) noise;
  /* Init locals. */
  float value = 0;
  int i, j;
  memcpy(tf, f, sizeof(float) * data->ndim);

  /* Inner loop for spectral construction, where the fractal is build. */
  for(i = 0; i < (int)octaves; i++) {
    value += noise_get(noise, tf) * data->exponent[i];
    for(j = 0; j < data->ndim; j++) tf[j] *= data->lacunarity;
  }

  /* Take care of remainder in octaves. */
  octaves -= (int)octaves;
  if(octaves > DELTA)
    value += octaves * noise_get(noise, tf) * data->exponent[i];
  return CLAMP(-0.99999f, 0.99999f, value);
}
#endif

static float noise_turbulence(noise_t noise, float* f, float octaves) {
  float tf[NOISE_MAX_DIMENSIONS];
  perlin_data_t* data = (perlin_data_t*) noise;
  /* Init locals. */
  float value = 0;
  int i, j;
  memcpy(tf, f, sizeof(float) * data->ndim);

  /* Inner loop of spectral construction, where the fractal is built. */
  for(i = 0; i < (int)octaves; i++) {
    value += ABS(noise_get(noise, tf)) * data->exponent[i];
    for(j = 0; j < data->ndim; j++) tf[j] *= data->lacunarity;
  }

  /* Take care of remainders in octaves. */
  octaves -= (int)octaves;
  if(octaves > DELTA)
    value += octaves * ABS(noise_get(noise, tf)) * data->exponent[i];
  return CLAMP(-0.99999f, 0.99999f, value);
}

void noise_delete(noise_t noise) {
  free((perlin_data_t*)noise);
}

static float* noise_genNebulae(const int w, const int h, const int n, float rug) {
  int x, y, z;
  float f[3];
  float octaves;
  float hurst;
  float lacunarity;
  noise_t noise;
  float* nebulae;;
  float value;

  octaves = 3.;
  hurst = NOISE_DEFAULT_HURST;
  lacunarity = NOISE_DEFAULT_LACUNARITY;

  noise = noise_new(2, hurst, lacunarity);

  nebulae = malloc(sizeof(float)*w*h*n);
  if(nebulae == NULL) {
    WARN("Out of memory!");
    return NULL;
  }

  for(z = 0; z < n; z++) {
    for(y = 0; y < h; y++) {
      for(x = 0; x < w; x++) {
        
        f[0] = rug * (float)x / (float)w;
        f[1] = rug * (float)y / (float)h;
        f[2] = rug * (float)z / (float)n;

        value = noise_turbulence(noise, f, octaves);

        value = value + 0.3;
        nebulae[z*w*h + y*w+x] = (value < 1.) ? value : 1.;
      }
    }
  }

  noise_delete(noise);

  return nebulae;
}

glTexture* noise_genCloud(const int w, const int h, double rug) {
  int i;
  float* map;
  SDL_Surface* sur;
  uint32_t* pix;
  glTexture* tex;
  double c;

  map = noise_genNebulae(w, h, 1, rug);

  sur = SDL_CreateRGBSurface(SDL_SWSURFACE, w, h, 32, RGBMASK);
  pix = sur->pixels;

  /* Convert from mapping to actual colours. */
  SDL_LockSurface(sur);
  for(i = 0; i < h*w; i++) {
    c = map[i];
    pix[i] = RMASK + BMASK + GMASK + (AMASK & (uint32_t)(AMASK*c));
  }
  SDL_UnlockSurface(sur);

  free(map);

  tex = gl_loadImage(sur);

  return tex;
}