#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

png_image image;
png_bytep buffer[2];
double a, rx[WAVE_COUNT], ry[WAVE_COUNT], phase[WAVE_COUNT],
       scale[WAVE_COUNT], cx[WAVE_COUNT], cy[WAVE_COUNT];
int s, i, p, x, y, *error[3], e;

double Rnd(double range)
{
   return (double)rand() / RAND_MAX * range;
}

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;
      for(s = image.width * image.height * 4; i < 3;
          error[i++] = (int*)calloc(image.width + 3, sizeof(int)))
         if( i < 2 )
             buffer[i] = (png_bytep)malloc(s);
      p = *buffer && buffer[1] &&
          *error && error[1] && error[2] &&
          png_image_finish_read(&image, NULL, *buffer, 0, NULL);
   }
   if( !p )
      return printf("Error reading %s\n", argv[1]);

   srand(time(NULL));

   /* Generate wave parameters. */
   for(memcpy(buffer[1], *buffer, s); x < WAVE_COUNT;
       phase[x++] = Rnd(2.0 * PI))
   {
      rx[x] = cos(a = Rnd(2.0 * PI)) *
              (scale[x] =
                 (Rnd(40.0) + 6.0) /
                 (image.width < image.height ? image.width : image.height));
      ry[x] = sin(a) * scale[x];

      cx[x] = Rnd(image.width);
      cy[x] = Rnd(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++)
            a +=
                sizeof('c') > 1
                ?  /* C mode: linear moire. */
                   sin(x * rx[i] + y * ry[i] + phase[i]) + 1.0
                :  /* C++ mode: circular moire. */
                   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 = y % 3][x + 1] / 8;

         memset(buffer[p > 127 ? 0 : 1] + s * 4, 0, 4);

         /* Reset previously accumulated error here for next scanline. */
         error[e][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
                                                                       */

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

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