#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. */
double gradient[GRID_SIZE + 1][GRID_SIZE + 1][2];

png_image image;
png_bytep buffer[2];
int s, x, y, grid_x, grid_y, w, h;
double d, dx, dy;

double rnd()
{
   return (double)rand() / RAND_MAX - 0.5;
}

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

double Interpolate(double a, double b, double x)
{
   return x < 0
      ? a
      : x > 1 ? b
              : /* https://en.wikipedia.org/wiki/Smoothstep */
                (b - a) * (x * (x * 6 - 15) + 10) * x * x * x + a;
}

int main(int argc, char **argv)
{
   if( argc - 4 )
      return printf("%s <input.png> <output1.png> <output2.png>\n", *argv);

   image.version = PNG_IMAGE_VERSION;
   if( png_image_begin_read_from_file(&image, argv[1]) )
   {
      image.format = PNG_FORMAT_RGBA;
      s = (w = image.width) * (h = image.height) * 4;
      for(x = 0; x < 2; x++)
         buffer[x] = (png_bytep)malloc(s);
      y = *buffer && buffer[1] &&
          png_image_finish_read(&image, NULL, *buffer, 0, NULL);
   }
   if( !y )
      return printf("Error reading %s\n", argv[1]);

   memcpy(buffer[1], *buffer, s);
   srand(time(NULL));

   /* Generate gradients. */
   for(y = 0; y <= GRID_SIZE; y++)
      for(x = 0; x <= GRID_SIZE; gradient[y][x++][1] = dy / d)
      {
         for(dx = dy = 0; (d = hypot(dx, dy)) < 1e-6; dy = rnd())
            dx = rnd();
         gradient[y][x][0] = dx / d;
      }
   d = w > h ? (double)w / GRID_SIZE : (double)h / GRID_SIZE;

   for(y = s = 0; y < h; y++)
      for(grid_y = floor(dy = y / d), x = 0; x < w; x++, s += 4)
      {
         grid_x = floor(dx = x / d);

         memset(
            buffer
            [
               (sizeof('c') > 1
                 ? /* Compare with white noise level in C mode. */
                   0
                 : /* Compare with Perlin noise in C++ mode.
                      https://en.wikipedia.org/wiki/Perlin_noise */
                   Interpolate(
                      Interpolate(DotProduct(gradient[grid_y][grid_x][0],
                                             gradient[grid_y][grid_x][1],
                                             grid_x - dx,
                                             grid_y - dy),
                                  DotProduct(gradient[grid_y][grid_x + 1][0],
                                             gradient[grid_y][grid_x + 1][1],
                                             grid_x + 1 - dx,
                                             grid_y - dy),
                                  dx - grid_x),
                      Interpolate(DotProduct(gradient[grid_y + 1][grid_x][0],
                                             gradient[grid_y + 1][grid_x][1],
                                             grid_x - dx,
                                             grid_y + 1 - dy),
                                  DotProduct(gradient[grid_y + 1][grid_x + 1][0],
                                             gradient[grid_y + 1][grid_x + 1][1],
                                             grid_x + 1 - dx,
                                             grid_y + 1 - dy),
                                  dx - grid_x),
                      dy - grid_y)) > rnd() * 2 ? 0 : 1
            ] + s,
            0,
            4);
      }

   for(x = 0; x < 2; x++)
      if( !png_image_write_to_file(&image, argv[2 + x], 0, buffer[x], 0, NULL) )
         return printf("Error writing %s\n", argv[2 + x]);
   return 0;
}