namespace Noise;

using System;

// Used for the hashing example that's commented out
//using System.Text;
//using System.Security.Cryptography;

// Based on https://adrianb.io/2014/08/09/perlinnoise.html

public class PerlinNoise {
  // Hash lookup table as defined by Ken Perlin.  This is a randomly
  // arranged array of all numbers from 0-255 inclusive.
  // Works as the seed.
  private static readonly int[] permutation = { 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
  };

  // Later on, in the hash function, access similar to p[p[_]+255+1] is performed
  // p contains two copies of permutation, len 512 in total.
  // The max number p[n] can return is 255, and 255+255+1 is 511,
  // which is the largest index in the 512 long array.
  private static readonly int[] p;

  // Controls intentional noise repeating
  // Do note the noise will repeat each 256 anyway, due to the hashing function 
  // (to be exact, the preprocessing of the input for the hashing function)
  private static readonly int repeat = 0;

  static PerlinNoise() {
    p = new int[512];
    for(int i=0;i<512;i++) {
        p[i] = permutation[i%256];
    }
  }

  // TODO figure out how the vectors were chosen
  // TODO figure out whether there is bias due to a few vectors being repeated
  // TODO experiment with other vectors
  //
  // Choose a random vector from
  //
  // The 12 vectors
  // (1,1,0),(-1,1,0),(1,-1,0),(-1,-1,0),
  // (1,0,1),(-1,0,1),(1,0,-1),(-1,0,-1),
  // (0,1,1),(0,-1,1),(0,1,-1),(0,-1,-1)
  //
  // Then calculate the dot product. 
  // Take as an example the first one, (1,1,0)
  // The dot product would thus be 1x + 1y + 0z. Thus, x + y
  //
  // The equations below are just simplified dot products.
  //
  // We have 12 vectors but 16 cases, so some of the vectors repeat.
  //
  // The vectors represent the middles of edges of a cube. 
  // A cube has 12 edges, so 12 vectors.
  public static double grad(int hash, double x, double y, double z) {
    switch(hash & 0xF)
    {
      // Unique
      case 0x0: return  x + y; //  1,  1,  0 A
      case 0x1: return -x + y; // -1,  1,  0 B
      case 0x2: return -y + z; //  0, -1,  1 C
      case 0x3: return -x - y; // -1, -1,  0
      case 0x4: return  x + z; //  1,  0,  1
      case 0x5: return -x + z; // -1,  0,  1
      case 0x6: return  x - z; //  1,  0, -1
      case 0x7: return -x - z; // -1,  0, -1
      case 0x8: return  y + z; //  0,  1,  1
      case 0x9: return  x - y; //  1, -1,  0
      case 0xA: return  y - z; //  0,  1, -1
      case 0xB: return -y - z; //  0, -1, -1
      // Repetitions
      case 0xC: return  y + x; //  1,  1,  0 A
      case 0xD: return -y + z; //  0, -1,  1 C
      case 0xE: return  y - x; // -1,  1,  0 B
      case 0xF: return -y - z; //  0, -1,  1 C
      default: return 0; // never happens
    }
  }

  // TODO Unsure of other properties it needs to satisfy
  // TODO Unsure why it's required, but without it, 
  //      edges between cells/cubes become visible
  // TODO There are indeed other properties.
  //      1 / (1 + Math.Exp(-12 * (x - 0.5)))
  //      satisfies the 3 defined ones, but does not give good results.
  //
  // A magical easing function
  //
  // Intended range is 0..1 -> 0..1
  // It's a simple sigmoid that satisfies 
  // f(0) = 0
  // f(1) = 1
  //
  // Fade function as defined by Ken Perlin.  This eases coordinate values
  // so that they will ease towards integral values.  This ends up smoothing
  // the final output.
  public static double fade(double t) {
    // 6t^5 - 15t^4 + 10t^3
    return t * t * t * (t * (t * 6 - 15) + 10);
  }

  // Incrementing function that wraps around based on `repeat`
  public static int inc(int num) {
    num++;
    if (repeat > 0) num %= repeat;
    
    return num;
  }

  // Linear Interpolate
  public static double lerp(double a, double b, double x) {
    return a + x * (b - a);
  }

  // Hash function
  // Any hash function can be utilized (for example sha256 works), though this one is fast and simple.
  // This hash function performs 3 random "jumps" in the permutations array, 
  // seeded & offset by the coordinates and offset by the number it landed on.
  private static int hash(int xi, int yi, int zi) {
    return p[p[p[xi]+yi]+zi];
  }

  // Here is a Sha256 replacement, to show that the permutations are literally just a prng with nothing special
//private static int hash(int xi, int yi, int zi) {
//    string inputString = $"{xi},{yi},{zi}";
//    byte[] hashBytes = SHA256.HashData(Encoding.UTF8.GetBytes(inputString));
//    return BitConverter.ToInt32(hashBytes, 0);
//}

  // TODO figure out why n.0 always returns 0.5, and n.5 always returns numbers like 
  // 0.75, 0.5 or 0.375. Overall, the noise seems to break down when it is not "zoomed in" enough.
  public static double perlin(double x, double y, double z) {
    // TODO this is probably useless
    // If we have any repeat on, change the coordinates to their "local" repetitions
    if(repeat > 0) {
      x = x%repeat;
      y = y%repeat;
      z = z%repeat;
    }
    
    // Calculate the "unit cube" that the point asked will be located in
    // The left bound is ( |_x_|,|_y_|,|_z_| ) and the right bound is that
    // plus 1.  Next we calculate the location (from 0.0 to 1.0) in that cube.

    // This is just a singular corner of the cube.
    // Modulos with 256 - this causes integers to repeat every 256.
    // Thus, 0.25 = 256.25, 3.6 = 259.6...
    // And yes, this causes a visual repeating pattern. 
    // The edge is not visible under normal conditions, though.
    int xi = (int)x & 255;
    int yi = (int)y & 255;
    int zi = (int)z & 255;

    // 0..1 coordinates within the cube
    double xf = x-(int)x;
    double yf = y-(int)y;
    double zf = z-(int)z;

    // TODO figure out why the lerping coordinates are eased rather than anything else.
    // Coordinates used for the lerping
    double u = fade(xf);
    double v = fade(yf);
    double w = fade(zf);

    // Magic hash
    // Compute 8 hashes for each vertex of the cube, from each coordinate
    int aaa, aba, aab, abb, baa, bba, bab, bbb;
    aaa = hash(    xi,      yi,      zi);
    aba = hash(    xi,  inc(yi),     zi);
    aab = hash(    xi,      yi,  inc(zi));
    abb = hash(    xi,  inc(yi), inc(zi));
    baa = hash(inc(xi),     yi,      zi);
    bba = hash(inc(xi), inc(yi),     zi);
    bab = hash(inc(xi),     yi,  inc(zi));
    bbb = hash(inc(xi), inc(yi), inc(zi));

    // TODO figure out how the -1s work
    // TODO the u, v, and w being used in this order matters. Why?
    // TODO figure out why lerping, and this particular way works
    //
    // This is then used as the seed for the grad() function that uses it to
    // pick a random edge vector and compute the dot product between it and the input vector,
    // offset into 8 corners, where one is the origin corner of the cube
    // (This works fine because the origin corner of the cube is also the origin of the vectors)
    // Results in vectors pointing from the corners to the input point specified by the input vector.
    //
    // The gradient function calculates the dot product between a pseudorandom                  
    // gradient vector and the vector from the input coordinate to the 8                        
    // surrounding points in its unit cube.                                                     
    // This is all then lerped together as a sort of weighted average based on the faded (u,v,w)
    // values we made earlier.                                                                  
    double x1, x2, x3, x4, y1, y2;
    x1 = lerp(grad(aaa, xf  , yf  , zf  ),           
              grad(baa, xf-1, yf  , zf  ), u);                                     
    x2 = lerp(grad(aba, xf  , yf-1, zf  ),           
              grad(bba, xf-1, yf-1, zf  ), u);
    y1 = lerp(x1, x2, v);

    x3 = lerp(grad(aab, xf  , yf  , zf-1),
              grad(bab, xf-1, yf  , zf-1), u);
    x4 = lerp(grad(abb, xf  , yf-1, zf-1),
              grad(bbb, xf-1, yf-1, zf-1), u);
    y2 = lerp(x3, x4, v);

    // For convenience we bind the result to 0 - 1 (theoretical min/max before is [-1, 1])
    double res = (lerp (y1, y2, w)+1)/2;

    return res;
  }
}