#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;
      w = image.width;
      h = image.height;
      for(x = 0; x < 2;)
         buffer[x++] = (png_bytep)malloc(s = w * h * 4);
      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 = (double)(w > h ? w : 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.

                      rnd() returns a random number in [-0.5, 0.5] range.
                      Setting threshold to zero here will result in
                      balanced distribution of dots for the two output
                      images, which is not as interesting as an unbalanced
                      threshold here.                                      */
                   0.25
                 : /* 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))
                      /* Perlin noise results in a number in the range
                        of [-1, 1], while rnd() returns a range of [-0.5, 0.5].

                        To get more dispersed dots, we should multiply
                        rnd() here by 2, but by not doing that, we get
                        more of a clustering behavior in the [-1, -0.5]
                        and [0.5, 1] ranges, which looks nicer for
                        our use case.

                        We can get larger clusters by reducing the
                        rnd() range further, but we will need to adjust
                        the fixed threshold for C mode accordingly.
                        It also doesn't look as good without some
                        scattered dots.                                 */
                      > rnd()
            ] + 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;
}