/* Halftone pattern generator.

   ./halftone <width> <height> > output.pgm
                                            */
#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

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;
   double angle, scale, min_x, max_x;
   Transform halftone, gradient;
   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));

   /* 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). */
   scale = (double)rand() / RAND_MAX * (MAX_DOT_COUNT - MIN_DOT_COUNT)
         + MIN_DOT_COUNT;
   scale /= width > height ? width : height;
   angle = ((double)rand() / RAND_MAX) * 0.5 * PI;
   halftone.a = cos(angle) * scale;
   halftone.d = sin(angle) * scale;
   halftone.c = (double)rand() / RAND_MAX;
   halftone.b = -halftone.d;
   halftone.e = halftone.a;
   halftone.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.c = gradient.d = gradient.e = gradient.f = 0;

   /* Start by picking a random angle without scaling or translation. */
   angle = (double)rand() / RAND_MAX * 2.0 * PI;
   gradient.a = cos(angle);
   gradient.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);
   min_x = max_x = ApplyTransformX(&p2, &gradient);

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

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

   p1.y = 0;
   p2 = ApplyTransform(&p1, &halftone);
   UpdateRange(ApplyTransformX(&p2, &gradient), &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.a *= scale;
   gradient.b *= scale;
   gradient.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;
         p2 = ApplyTransform(&p1, &halftone);
         if( WithinRadiusOfNearbyDots(&p2, &gradient) )
         {
            fputc(255, stdout);
         }
         else
         {
            fputc(0, stdout);
         }
      }
   }
   return 0;
}