/* Halftone pattern generator.

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

   Forked from halftone.c to generate multiple layers of halftone dots.
                                                                        */
#include<assert.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
#ifdef _WIN32
   #include<fcntl.h>
   #include<io.h>
#endif

#define PI             3.14159265358979323846264338327950288419716939937510
#define MIN_DOT_COUNT  30.0
#define MAX_DOT_COUNT  120.0
#define LAYER_COUNT    3

typedef struct { double x, y; } XY;
typedef struct { double a, b, c, d, e, f; } Transform;

/* Apply transformations to point. */
static double ApplyTransformX(XY *input, Transform *t)
{
   return input->x * t->a + input->y * t->b + t->c;
}
static double ApplyTransformY(XY *input, Transform *t)
{
   return input->x * t->d + input->y * t->e + t->f;
}
static XY ApplyTransform(XY *input, Transform *t)
{
   XY output;
   output.x = ApplyTransformX(input, t);
   output.y = ApplyTransformY(input, t);
   return output;
}

/* Update value range. */
static void UpdateRange(double input, double *range_min, double *range_max)
{
   if( *range_min > input ) *range_min = input;
   if( *range_max < input ) *range_max = input;
}

/* Determine if a point is within the radius of a particular target,
   returns 1 if so.                                                  */
static int WithinDotRadius(XY *point, XY *target, Transform *gradient)
{
   double r = ApplyTransformX(target, gradient);
   double dx = point->x - target->x;
   double dy = point->y - target->y;
   return dx * dx + dy * dy < r * r;
}

/* Determine if a point is within the radius of nearby halftone dots,
   returns 1 if so.                                                   */
static int WithinRadiusOfNearbyDots(XY *point, Transform *gradient)
{
   XY target;

   target.x = floor(point->x);
   target.y = floor(point->y);
   if( WithinDotRadius(point, &target, gradient) )
      return 1;

   target.y = ceil(point->y);
   if( WithinDotRadius(point, &target, gradient) )
      return 1;

   target.x = ceil(point->x);
   if( WithinDotRadius(point, &target, gradient) )
      return 1;

   target.y = floor(point->y);
   return WithinDotRadius(point, &target, gradient);
}

int main(int argc, char **argv)
{
   int width, height, x, y, i;
   double dot_count, angle, scale, min_x, max_x;
   Transform halftone[LAYER_COUNT], gradient[LAYER_COUNT];
   XY p1, p2;

   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.");

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

   srand(time(NULL));

   dot_count = (double)rand() / RAND_MAX * (MAX_DOT_COUNT - MIN_DOT_COUNT)
             + MIN_DOT_COUNT;
   for(i = 0; i < LAYER_COUNT; i++)
   {
      /* Generate transformation from screen coordinates to halftone
         coordinates.

         Halftone dots are centered at integer intervals in halftone
         coordinates, so the center between any two dots is at most sqrt(2).

         Halftone dot counts are adjusted such that we get more dots
         at lower layers.  Different dot frequencies at each layer looks
         nicer because it makes the output less periodic.                   */
      scale = pow(1.2, i) * dot_count / (width > height ? width : height);
      angle = ((double)rand() / RAND_MAX) * 0.5 * PI;
      halftone[i].a = cos(angle) * scale;
      halftone[i].d = sin(angle) * scale;
      halftone[i].c = (double)rand() / RAND_MAX;
      halftone[i].b = -halftone[i].d;
      halftone[i].e = halftone[i].a;
      halftone[i].f = (double)rand() / RAND_MAX;

      /* Generate transformation from halftone coordinates to gradient
         coordinates.  The X value in gradient space is used to set
         halftone dot sizes.                                           */
      gradient[i].c = gradient[i].d = gradient[i].e = gradient[i].f = 0;

      /* Start by picking a random angle without scaling or translation. */
      angle = (double)rand() / RAND_MAX * 2.0 * PI;
      gradient[i].a = cos(angle);
      gradient[i].b = -sin(angle);

      /* Apply transforms to screen coordinates to get the range of X
         values in gradient coordinates.

         Note that the range is (0,0) .. (width,height), as opposed
         to the screen area of (0,0) .. (width-1,height-1).  This means
         the 4 points tested always enclose a rectangular area that is at
         least one pixel wide on each side, which means the min and max X
         values will never be equal and we won't have a divide by zero.   */
      p1.x = p1.y = 0;
      p2 = ApplyTransform(&p1, &halftone[i]);
      min_x = max_x = ApplyTransformX(&p2, &gradient[i]);

      p1.y = height;
      p2 = ApplyTransform(&p1, &halftone[i]);
      UpdateRange(ApplyTransformX(&p2, &gradient[i]), &min_x, &max_x);

      p1.x = width;
      p2 = ApplyTransform(&p1, &halftone[i]);
      UpdateRange(ApplyTransformX(&p2, &gradient[i]), &min_x, &max_x);

      p1.y = 0;
      p2 = ApplyTransform(&p1, &halftone[i]);
      UpdateRange(ApplyTransformX(&p2, &gradient[i]), &min_x, &max_x);

      assert(min_x < max_x);

      /* Update the gradient transforms so that output falls within
         [0, 0.5 * sqrt(2)] range.                                  */
      scale = 0.5 * sqrt(2) / (max_x - min_x);
      gradient[i].a *= scale;
      gradient[i].b *= scale;
      gradient[i].c = -min_x * scale;
   }

   /* Render pixels. */
   printf("P5\n%d %d\n255\n", width, height);
   for(y = 0; y < height; y++)
   {
      p1.y = y;
      for(x = 0; x < width; x++)
      {
         p1.x = x;
         for(i = 0; i < LAYER_COUNT; i++)
         {
            p2 = ApplyTransform(&p1, &halftone[i]);
            if( WithinRadiusOfNearbyDots(&p2, &gradient[i]) )
            {
               /* Assign the colors such that the topmost layer gets
                  the brightest color.                               */
               fputc((LAYER_COUNT - i) * 255 / LAYER_COUNT, stdout);
               break;
            }
         }
         if( i == LAYER_COUNT )
            fputc(0, stdout);
      }
   }
   return 0;
}