/* Circular moire pattern generator.

   ./circular_moire_dither_fs <width> <height> > output.pgm

   Forked from circular_moire.c to use Atkinson error diffusion.

   https://en.wikipedia.org/wiki/Atkinson_dithering
                                                                */
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#ifdef _WIN32
   #include<fcntl.h>
   #include<io.h>
#endif

#define CIRCLE_COUNT  3

int main(int argc, char **argv)
{
   int width, height, x, y, i, y0, y1, y2;
   double cx[CIRCLE_COUNT], cy[CIRCLE_COUNT], scale[CIRCLE_COUNT], a;
   int *row_error[3], o, e;

   if( argc != 3 ||
       (width = atoi(argv[1])) <= 0 ||
       (height = atoi(argv[2])) <= 0 )
   {
      return printf("%s <width> <height>\n", *argv);
   }
   if( width >= 0x8000 || height >= 0x8000 )
      return !puts("Output size too large.");

   row_error[0] = (int*)calloc(width + 3, sizeof(int));
   row_error[1] = (int*)calloc(width + 3, sizeof(int));
   row_error[2] = (int*)malloc((width + 3) * sizeof(int));
   if( row_error[0] == NULL || row_error[1] == NULL )
      return !puts("Not enough memory");

   #ifdef _WIN32
      setmode(STDOUT_FILENO, O_BINARY);
   #endif

   srand(time(NULL));

   for(i = 0; i < CIRCLE_COUNT; i++)
   {
      cx[i] = (double)rand() / RAND_MAX * width;
      cy[i] = (double)rand() / RAND_MAX * height;
      scale[i] = ((double)rand() / RAND_MAX * 30.0 + 6.0)
                  / (width < height ? width : height);
   }

   printf("P5\n%d %d\n255\n", width, height);
   for(y = 0; y < height; y++)
   {
      /* Clear the last error row at the start of each scanline. */
      y0 = y % 3;
      y1 = (y + 1) % 3;
      y2 = (y + 2) % 3;
      memset(row_error[y2], 0, (width + 3) * sizeof(int));

      /* Render scanline. */
      for(x = 0; x < width; x++)
      {
         a = 0;
         for(i = 0; i < CIRCLE_COUNT; i++)
            a += sin(hypot(x - cx[i], y - cy[i]) * scale[i]) + 1.0;
         a /= 2.0 * CIRCLE_COUNT;

         /* i = intended grayscale level. */
         i = (int)(a * 255) + row_error[y0][x + 1] / 8;

         /* o = output grayscale level. */
         o = i > 127 ? 255 : 0;
         fputc(o, stdout);

         /* Propagate error. */
         e = i - o;
         row_error[y0][x + 2] += e;
         row_error[y0][x + 3] += e;
         row_error[y1][x    ] += e;
         row_error[y1][x + 1] += e;
         row_error[y1][x + 2] += e;
         row_error[y2][x    ] += e;
      }
   }

   free(row_error[0]);
   free(row_error[1]);
   free(row_error[2]);
   return 0;
}