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

   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 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;
         p = (double)rand() / RAND_MAX > a;
         if( p )
            memset(buffer1 + (s * 4), 0, 4);
         else
            memset(buffer2 + (s * 4), 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;
}