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

#define WAVE_COUNT 3
#define PI 3.14159265358979323846264338327950288419716939937510

int main(int argc, char **argv)
{
   png_image image;
   png_bytep buffer1, buffer2;
   double a, rx[WAVE_COUNT], ry[WAVE_COUNT], phase[WAVE_COUNT];
   double scale[WAVE_COUNT], cx[WAVE_COUNT], cy[WAVE_COUNT];
   int s, i, p, x, y;
   int *error[3], e[3];

   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);
   for(i = 0; i < 3; i++)
   {
      e[i] = i;
      error[i] = (int*)calloc(image.width + 3, sizeof(int));
   }
   if( buffer1 == NULL || buffer2 == NULL ||
       error[0] == NULL || error[1] == NULL || error[2] == 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 wave parameters. */
   for(i = 0; i < WAVE_COUNT; i++)
   {
      a = (double)rand() / RAND_MAX * 2.0 * PI;
      scale[i] = ((double)rand() / RAND_MAX * 40.0 + 6.0) /
                 (image.width < image.height ? image.width : image.height);
      rx[i] = cos(a) * scale[i];
      ry[i] = sin(a) * scale[i];
      phase[i] = (double)rand() / RAND_MAX * 2.0 * PI;

      cx[i] = (double)rand() / RAND_MAX * image.width;
      cy[i] = (double)rand() / RAND_MAX * image.height;
   }

   for(y = s = 0; y < (int)image.height; y++)
   {
      for(x = 0; x < (int)image.width; x++, s++)
      {
         a = 0;
         for(i = 0; i < WAVE_COUNT; i++)
         {
            if( sizeof('c') > 1 )
            {
               /* C mode: linear moire. */
               a += sin(x * rx[i] + y * ry[i] + phase[i]) + 1.0;
            }
            else
            {
               /* C++ mode: circular moire. */
               a += sin(hypot(x - cx[i], y - cy[i]) * scale[i]) + 1.0;
            }
         }
         a /= 2.0 * WAVE_COUNT;

         /* Add pixel gray level plus accumulated error for
            current pixel, then render pixel.               */
         p = (int)(a * 255) + error[e[0]][x + 1] / 8;

         if( p > 127 )
            memset(buffer1 + (s * 4), 0, 4);
         else
            memset(buffer2 + (s * 4), 0, 4);

         /* Reset previously accumulated error here for next scanline. */
         error[e[0]][x + 1] = 0;

         /* Propagate current error to other pixels.

            This uses Atkinson dithering, as opposed to the smoother
            Floyd-Steinberg.  Atkinson generates some splotchy
            artifacts which are sometimes undesirable, but for our
            purposes we actually want those, since the highly diffused
            dot patterns from Floyd-Steinberg don't produce good masks.

            https://en.wikipedia.org/wiki/Atkinson_dithering
                                                                       */

         p = p - (p > 127 ? 255 : 0);
         error[e[0]][x + 2] += p;
         error[e[0]][x + 3] += p;
         error[e[1]][x] += p;
         error[e[1]][x + 1] += p;
         error[e[1]][x + 2] += p;
         error[e[2]][x + 1] += p;
      }
      for(i = 0; i < 3; i++)
         e[i] = (e[i] + 1) % 3;
   }

   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;
}