/*
 * PerlinNoise.hpp
 *
 *  Created on: Jul 11, 2008
 *      Author: pete
 */

#ifndef PERLINNOISE_HPP_
#define PERLINNOISE_HPP_

#include <math.h>


class PerlinNoise
{
public:
  float noise_lq[256][256];
  float noise_lq_lite[32][32];
  float noise_mq[256][256];
  float noise_hq[256][256];
  float noise_perlin[512][512];
  float noise_lq_vol[32][32][32];
  float noise_hq_vol[32][32][32];


  PerlinNoise();
  virtual ~PerlinNoise();

private:

  static inline float noise( int x)
  {
    x = (x<<13)^x;
    return (((x * (x * x * 15731 + 789221) + 1376312589) & 0x7fffffff) / 2147483648.0);
  }

  static inline float noise(int x, int y)
  {
    int n = x + y * 57;
    return noise(n);
  }

  static inline float noise(int x, int y, int z)
  {
    int n = x + y * 57 + z * 141;
    return noise(n);
  }

  static inline float cos_interp(float a, float b, float x)
  {
    float ft = x * 3.1415927;
    float f = (1 - cos(ft)) * .5;

    return  a*(1-f) + b*f;
  }

  static inline float cubic_interp(float v0, float v1, float v2, float v3, float x)
  {
    float P = (v3 - v2) - (v0 - v1);
    float Q = (v0 - v1) - P;
    float R = v2 - v0;

    return P*pow(x,3) + Q * pow(x,2) + R*x + v1;
  }

  static inline float InterpolatedNoise(float x, float y)
  {
    int integer_X = int(x);
    float fractional_X = x - integer_X;

    int integer_Y = int(y);
    float fractional_Y = y - integer_Y;

    float a0 = noise(integer_X - 1,      integer_Y - 1);
    float a1 = noise(integer_X,          integer_Y - 1);
    float a2 = noise(integer_X + 1,      integer_Y - 1);
    float a3 = noise(integer_X + 2,      integer_Y - 1);

    float x0 = noise(integer_X - 1,      integer_Y);
    float x1 = noise(integer_X,          integer_Y);
    float x2 = noise(integer_X + 1,      integer_Y);
    float x3 = noise(integer_X + 2,      integer_Y);

    float y0 = noise(integer_X + 0,      integer_Y + 1);
    float y1 = noise(integer_X,        integer_Y + 1);
    float y2 = noise(integer_X + 1,    integer_Y + 1);
    float y3 = noise(integer_X + 2,    integer_Y + 1);

    float b0 = noise(integer_X - 1,      integer_Y + 2);
    float b1 = noise(integer_X,          integer_Y + 2);
    float b2 = noise(integer_X + 1,      integer_Y + 2);
    float b3 = noise(integer_X + 2,      integer_Y + 2);

    float i0 = cubic_interp(a0 , a1, a2, a3, fractional_X);
    float i1 = cubic_interp(x0 , x1, x2, x3, fractional_X);
    float i2 = cubic_interp(y0 , y1, y2, y3, fractional_X);
    float i3 = cubic_interp(b0 , b1, b2, b3, fractional_X);

    return cubic_interp(i0, i1 , i2 , i3, fractional_Y);

  }

  static inline float perlin_octave_2d(float x,float y, int width, int seed, float period)
  {
    float freq=1/(float)(period);

    int num=(int)(width*freq);
    int step_x=(int)(x*freq);
    int step_y=(int)(y*freq);
    float zone_x=x*freq-step_x;
    float zone_y=y*freq-step_y;
    int box=step_x+step_y*num;
    int noisedata=(box+seed);

    float u=cubic_interp(noise(noisedata-num-1),noise(noisedata-num),noise(noisedata-num+1),noise(noisedata-num+2),zone_x);
    float a=cubic_interp(noise(noisedata-1),noise(noisedata),noise(noisedata+1),noise(noisedata+2),zone_x);
    float b=cubic_interp(noise(noisedata+num -1),noise(noisedata+num),noise(noisedata+1+num),noise(noisedata+2+num),zone_x);
    float v=cubic_interp(noise(noisedata+2*num -1),noise(noisedata+2*num),noise(noisedata+1+2*num),noise(noisedata+2+2*num),zone_x);

    float value=cubic_interp(u,a,b,v,zone_y);

    return value;
  }


  static inline float perlin_octave_2d_cos(float x,float y, int width, int seed, float period)
  {
    float freq=1/(float)(period);

    int num=(int)(width*freq);
    int step_x=(int)(x*freq);
    int step_y=(int)(y*freq);
    float zone_x=x*freq-step_x;
    float zone_y=y*freq-step_y;
    int box=step_x+step_y*num;
    int noisedata=(box+seed);

    float a=cos_interp(noise(noisedata),noise(noisedata+1),zone_x);
    float b=cos_interp(noise(noisedata+num),noise(noisedata+1+num),zone_x);

    float value=cos_interp(a,b,zone_y);

    return value;
  }


  static inline float perlin_octave_3d(float x,float y, float z,int width, int seed, float period)
  {
    float freq=1/(float)(period);

    int num=(int)(width*freq);
    int step_x=(int)(x*freq);
    int step_y=(int)(y*freq);
    int step_z=(int)(z*freq);
    float zone_x=x*freq-step_x;
    float zone_y=y*freq-step_y;
    float zone_z=z*freq-step_z;

    int boxB=step_x+step_y+step_z*num;
    int boxC=step_x+step_y+step_z*(num+1);
    int boxD=step_x+step_y+step_z*(num+2);
    int boxA=step_x+step_y+step_z*(num-1);

    float u,a,b,v,noisedata,box;

    box = boxA;
    noisedata=(box+seed);
    u=cubic_interp(noise(noisedata-num-1),noise(noisedata-num),noise(noisedata-num+1),noise(noisedata-num+2),zone_x);
    a=cubic_interp(noise(noisedata-1),noise(noisedata),noise(noisedata+1),noise(noisedata+2),zone_x);
    b=cubic_interp(noise(noisedata+num -1),noise(noisedata+num),noise(noisedata+1+num),noise(noisedata+2+num),zone_x);
    v=cubic_interp(noise(noisedata+2*num -1),noise(noisedata+2*num),noise(noisedata+1+2*num),noise(noisedata+2+2*num),zone_x);
    float A=cubic_interp(u,a,b,v,zone_y);

    box = boxB;
    noisedata=(box+seed);
    u=cubic_interp(noise(noisedata-num-1),noise(noisedata-num),noise(noisedata-num+1),noise(noisedata-num+2),zone_x);
    a=cubic_interp(noise(noisedata-1),noise(noisedata),noise(noisedata+1),noise(noisedata+2),zone_x);
    b=cubic_interp(noise(noisedata+num -1),noise(noisedata+num),noise(noisedata+1+num),noise(noisedata+2+num),zone_x);
    v=cubic_interp(noise(noisedata+2*num -1),noise(noisedata+2*num),noise(noisedata+1+2*num),noise(noisedata+2+2*num),zone_x);
    float B=cubic_interp(u,a,b,v,zone_y);

    box = boxC;
    noisedata=(box+seed);
    u=cubic_interp(noise(noisedata-num-1),noise(noisedata-num),noise(noisedata-num+1),noise(noisedata-num+2),zone_x);
    a=cubic_interp(noise(noisedata-1),noise(noisedata),noise(noisedata+1),noise(noisedata+2),zone_x);
    b=cubic_interp(noise(noisedata+num -1),noise(noisedata+num),noise(noisedata+1+num),noise(noisedata+2+num),zone_x);
    v=cubic_interp(noise(noisedata+2*num -1),noise(noisedata+2*num),noise(noisedata+1+2*num),noise(noisedata+2+2*num),zone_x);
    float C=cubic_interp(u,a,b,v,zone_y);

    box = boxD;
    noisedata=(box+seed);
    u=cubic_interp(noise(noisedata-num-1),noise(noisedata-num),noise(noisedata-num+1),noise(noisedata-num+2),zone_x);
    a=cubic_interp(noise(noisedata-1),noise(noisedata),noise(noisedata+1),noise(noisedata+2),zone_x);
    b=cubic_interp(noise(noisedata+num -1),noise(noisedata+num),noise(noisedata+1+num),noise(noisedata+2+num),zone_x);
    v=cubic_interp(noise(noisedata+2*num -1),noise(noisedata+2*num),noise(noisedata+1+2*num),noise(noisedata+2+2*num),zone_x);
    float D=cubic_interp(u,a,b,v,zone_y);

    float value =cubic_interp(A,B,C,D,zone_z);

    return value;
  }


  static inline float perlin_noise_2d(int x, int y,  int width, int octaves, int seed, float persistance, float basePeriod)
  {
    float p = persistance;
    float val = 0.0;

    for (int i = 0; i<octaves;i++)
    {
      val += perlin_octave_2d_cos(x,y,width,seed,basePeriod) * p;

      basePeriod *= 0.5;
      p *= persistance;
    }
    return val;
  }

  static inline float perlin_noise_3d(int x, int y, int z, int width, int octaves, int seed, float persistance, float basePeriod)
  {
    float p = persistance;
    float val = 0.0;

    for (int i = 0; i<octaves;i++)
    {
      val += perlin_octave_3d(x,y,z,width,seed,basePeriod) * p;

      basePeriod *= 0.5;
      p *= persistance;
    }
    return val;
  }
};

#endif /* PERLINNOISE_HPP_ */