// Noise.cpp

#include "stdafx.h"
#include "Noise.h"
#include "..\..\Global.h"
#include "..\..\Twister.h"

//_____________________________________________________________________________
// One Dimensional:
double Noise::NoiseFunction(DWORD x) {x=(x<<13)^x; return Signed(((x*(x*x*15731+789221)+1376312589)&LONG_MAX)/double(LONG_MAX));} // returns in the interval [0,1]

// Mosaic: Each coordinate gives a specific number.
// Mosaic(x/17,y/7,0.5); will create random rectangles 17 by 7 pixels in size
double Noise::Mosaic(DWORD x, DWORD y, DWORD z) {return NoiseFunction(x + 1301*y + 314159*z);}

// One Dimensional Perlin Noise:
double Noise::SmoothedNoise(DWORD x) {return NoiseFunction(x)/2 + NoiseFunction(x-1)/4 + NoiseFunction(x+1)/4;}
double Noise::InterpolatedNoise(double x) {DWORD X=Floor(x); return Interpolate(x-X, SmoothedNoise(X), SmoothedNoise(X+1));}
double Noise::Get(double x, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result +=Amplitude * InterpolatedNoise(Frequency * x);
  }
  return Result;
}
double Noise::Turbulence(double x, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result+=fabs(Amplitude * InterpolatedNoise(Frequency * x));
  }
  return Result;
}
//_____________________________________________________________________________
// Two Dimensional Perlin Noise:
double Noise::        NoiseFunction(DWORD x, DWORD y) {return NoiseFunction(x+57*y);} // returns in the interval [-1,1]
double Noise::SmoothedNoise(DWORD x, DWORD y) {
  return (NoiseFunction(x-1, y-1)+NoiseFunction(x+1, y-1)+NoiseFunction(x+1, y+1)+NoiseFunction(x-1, y+1))/16
       + (NoiseFunction(x-1, y  )+NoiseFunction(x+1, y  )+NoiseFunction(x  , y-1)+NoiseFunction(x  , y+1))/8
       +  NoiseFunction(x,y)/4;
}
double Noise::InterpolatedNoise(double x, double y) {
  DWORD X=Floor(x), Y=Floor(y);
  return Interpolate(y-Y, Interpolate(x-X, SmoothedNoise(X,Y  ), SmoothedNoise(X+1,Y  )),
                          Interpolate(x-X, SmoothedNoise(X,Y+1), SmoothedNoise(X+1,Y+1)));
}
double Noise::Get(double x, double y, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result +=Amplitude * InterpolatedNoise(Frequency * x, Frequency * y);
  }
  return Result;
}
double Noise::Turbulence(double x, double y, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result+=fabs(Amplitude * InterpolatedNoise(Frequency * x, Frequency * y));
  }
  return Result;
}
//_____________________________________________________________________________
// Three Dimensional Perlin Noise:
double Noise::InterpolatedNoise(double x, double y, double z) { // Returns interval [-1,1]
  int X=Floor(x); // Locate the unit Voxel:
  int Y=Floor(y); // Use as a coordinate in the p array
  int Z=Floor(z);
  x-=X; y-=Y; z-=Z; // x,y,z are now in the Interval [0,1] within the unit Voxel.
  double u=x*x*x*(x*(x*6-15)+10); // Triangle to normal curve: 6t^5 - 15t^4 + 10t^3
  double v=y*y*y*(y*(y*6-15)+10);
  double w=z*z*z*(z*(z*6-15)+10);
  const BYTE p[]={151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
/*Twice to wrap*/ 151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180};
  X&=0xFF; // Use as a coordinate in the p array
  Y&=0xFF;
  Z&=0xFF;
  int A=p[X  ]+Y, AA=p[A]+Z, AB=p[A+1]+Z; // Hash coordinates
  int B=p[X+1]+Y, BA=p[B]+Z, BB=p[B+1]+Z; // of the 8 cube corners
  return Interpolate(w, Interpolate(v, Interpolate(u, Grad(p[AA  ], x  , y  , z  ), // Interpolate coordinates from each corner of the cube:
                                                      Grad(p[BA  ], x-1, y  , z  )),
                                       Interpolate(u, Grad(p[AB  ], x  , y-1, z  ),
                                                      Grad(p[BB  ], x-1, y-1, z  ))),
                        Interpolate(v, Interpolate(u, Grad(p[AA+1], x  , y  , z-1),
                                                      Grad(p[BA+1], x-1, y  , z-1)),
                                       Interpolate(u, Grad(p[AB+1], x  , y-1, z-1),
                                                      Grad(p[BB+1], x-1, y-1, z-1))));
}
double Noise::Grad(int Hash, double x, double y, double z) {
  int h=Hash & 0xF; // Generate 12 gradient directions from least significan nybble
  double u=(h<8 ? x : y);
  double v=(h<4 ? y : ((h==12)||(h==14) ? x : z));
  return ((h&1)==0 ? u : -u) + ((h&2)==0 ? v : -v);
}

double Noise::Get(double x, double y, double z, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result+=Amplitude * InterpolatedNoise(Frequency * x, Frequency * y, Frequency * z);
  }
  return Result;
}
double Noise::Turbulence(double x, double y, double z, int Octaves, double Persistence) { // Persistence [0-1]. 1=All octaves added equally.
  double Result=0;
  DWORD Frequency=1; // Can't see amplitudes less than 1/256.
  for(double Amplitude=1; (Octaves--) && (Amplitude>1/256.); Amplitude*=Persistence, Frequency<<=1) {
    Result+=fabs(Amplitude * InterpolatedNoise(Frequency * x, Frequency * y, Frequency * z));
  }
  return Result;
}

double Noise::Stripes(double x, double y, double z, int Octaves, double Persistence) {return sin(y-3*Turbulence(x,y,z,Octaves,Persistence));} // Persistence [0-1]. 1=All octaves added equally.
double Noise::Marble (double x, double y, int Octaves, double Persistence) {return sin(50*(x+Get(x*5, y*5,Octaves,Persistence)/2));} // Persistence [0-1]. 1=All octaves added equally.