#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<png.h>

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

/* Gradient unit vectors. */
static double gradient[GRID_SIZE + 1][GRID_SIZE + 1][2];

static double DotProduct(double ax, double ay, double bx, double by)
{
   return ax * bx + ay * by;
}

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

int main(int argc, char **argv)
{
   png_image image;
   png_bytep buffer1, buffer2;
   int s, x, y, p, grid_x, grid_y;
   double d, dx, dy, cell_size, dot1, dot2, p1, p2;

   if( argc != 4 )
      return printf("%s <input.png> <output1.png> <output2.png>\n", *argv);

   memset(&image, 0, sizeof(image));
   image.version = PNG_IMAGE_VERSION;
   if( !png_image_begin_read_from_file(&image, argv[1]) )
      return printf("Error reading %s\n", argv[1]);

   image.format = PNG_FORMAT_RGBA;
   s = image.width * image.height * 4;
   buffer1 = (png_bytep)malloc(s);
   buffer2 = (png_bytep)malloc(s);
   if( buffer1 == NULL || buffer2 == NULL )
      return puts("Out of memory");
   if( !png_image_finish_read(&image, NULL, buffer1, 0, NULL) )
      return printf("Error loading %s\n", argv[1]);

   memcpy(buffer2, buffer1, s);
   srand(time(NULL));

   /* Generate 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;
      }
   }
   cell_size = image.width > image.height ? (double)image.width / GRID_SIZE
                                          : (double)image.height / GRID_SIZE;

   for(y = 0; y < (int)image.height; y++)
   {
      dy = (double)y / cell_size;
      grid_y = (int)floor(dy);
      for(x = 0; x < (int)image.width; x++)
      {
         if( sizeof('c') > 1 )
         {
            /* Use white noise pattern in C mode. */
            p = (rand() > RAND_MAX / 2);
         }
         else
         {
            /* Use Perlin noise in C++ mode.
               https://en.wikipedia.org/wiki/Perlin_noise */
            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 = Interpolate(p1, p2, dy - grid_y) * 0.5 + 0.5 >=
                (double)rand() / RAND_MAX;
         }

         s = (y * image.width + x) * 4;
         if( p )
            memset(buffer1 + s, 0, 4);
         else
            memset(buffer2 + s, 0, 4);
      }
   }

   if( !png_image_write_to_file(&image, argv[2], 0, buffer1, 0, NULL) )
      return printf("Error writing %s\n", argv[2]);
   if( !png_image_write_to_file(&image, argv[3], 0, buffer2, 0, NULL) )
      return printf("Error writing %s\n", argv[3]);
   return 0;
}