/* Perlin noise generator.

   ./perlin_noise <width> <height> > output.pgm

   https://en.wikipedia.org/wiki/Perlin_noise
                                                */
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#ifdef _WIN32
   #include<fcntl.h>
   #include<io.h>
#endif

/* Number of grid divisions along the long edge of the image. */
#define GRID_SIZE 16

/* Random gradients. */
static double gradient[GRID_SIZE + 1][GRID_SIZE + 1][2];

/* Dot product of two vectors. */
static double DotProduct(double ax, double ay, double bx, double by)
{
   return ax * bx + ay * by;
}

/* Interpolate between two gradient values.

   https://en.wikipedia.org/wiki/Smoothstep
                                            */
static double Interpolate(double a, double b, double x)
{
   if( x <= 0.0 ) return a;
   if( x >= 1.0 ) return b;

#if 0
   return (b - a) * (3.0 - x * 2.0) * x * x + a;
#else
   return (b - a) * ((x * (x * 6.0 - 15.0) + 10.0) * x * x * x) + a;
#endif
}

int main(int argc, char **argv)
{
   int width, height, x, y, grid_x, grid_y, p;
   double cell_size, d, dx, dy, dot1, dot2, p1, p2;

   if( argc != 3 ||
       (width = atoi(argv[1])) <= 0 ||
       (height = atoi(argv[2])) <= 0 )
   {
      return printf("%s <width> <height>\n", *argv);
   }
   if( width >= 0x8000 || height >= 0x8000 )
      return !puts("Output size too large.");

   #ifdef _WIN32
      setmode(STDOUT_FILENO, O_BINARY);
   #endif

   srand(time(NULL));

   /* Determine size of each grid cell such that the long edge of the
      output image is covered by exactly the number of grid cells we have. */
   cell_size = width > height ? (double)width / GRID_SIZE
                              : (double)height / GRID_SIZE;
   if( cell_size < 1 )
      cell_size = 1;

   /* Generate random gradients. */
   for(y = 0; y <= GRID_SIZE; y++)
   {
      for(x = 0; x <= GRID_SIZE; x++)
      {
         do
         {
            dx = (double)rand() / RAND_MAX - 0.5;
            dy = (double)rand() / RAND_MAX - 0.5;
            d = hypot(dx, dy);
         } while( d < 1e-6 );
         gradient[y][x][0] = dx / d;
         gradient[y][x][1] = dy / d;
      }
   }

   /* Render points. */
   printf("P5\n%d %d\n255\n", width, height);
   for(y = 0; y < height; y++)
   {
      dy = (double)y / cell_size;
      grid_y = (int)floor(dy);
      for(x = 0; x < width; x++)
      {
         dx = (double)x / cell_size;
         grid_x = (int)floor(dx);

         dot1 = DotProduct(gradient[grid_y][grid_x][0],
                           gradient[grid_y][grid_x][1],
                           dx - grid_x,
                           dy - grid_y);
         dot2 = DotProduct(gradient[grid_y][grid_x + 1][0],
                           gradient[grid_y][grid_x + 1][1],
                           dx - (grid_x + 1),
                           dy - grid_y);
         p1 = Interpolate(dot1, dot2, dx - grid_x);

         dot1 = DotProduct(gradient[grid_y + 1][grid_x][0],
                           gradient[grid_y + 1][grid_x][1],
                           dx - grid_x,
                           dy - (grid_y + 1));
         dot2 = DotProduct(gradient[grid_y + 1][grid_x + 1][0],
                           gradient[grid_y + 1][grid_x + 1][1],
                           dx - (grid_x + 1),
                           dy - (grid_y + 1));
         p2 = Interpolate(dot1, dot2, dx - grid_x);

         p = (int)((Interpolate(p1, p2, dy - grid_y) * 0.5 + 0.5) * 255.0);
         fputc(p < 0 ? 0 : p > 255 ? 255 : p, stdout);
      }
   }
   return 0;
}