/* Distortion illusion generator.

   ./distortion <width> <height> > output.pgm
                                              */
#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_SQUARE_COUNT  20.0
#define MAX_SQUARE_COUNT  50.0
#define SQUARE_CORNER     0.25

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

/* Check if coordinate within sub-square are near the corners,
   returns 1 if so.                                            */
static int IsNearCorner(double x, double y)
{
   return x + y < SQUARE_CORNER ||
          x + y > 2.0 - SQUARE_CORNER ||
          y - x > 1.0 - SQUARE_CORNER ||
          x - y > 1.0 - SQUARE_CORNER;
}

int main(int argc, char **argv)
{
   int width, height, x, y, i, parity, phase[2][2], period[2][2], d, w;
   double angle, scale;
   Transform square;
   XY p1, p2, p3;

   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 square coordinates. */
   scale = (double)rand() / RAND_MAX * (MAX_SQUARE_COUNT - MIN_SQUARE_COUNT)
         + MIN_SQUARE_COUNT;
   scale /= width > height ? width : height;
   angle = ((double)rand() / RAND_MAX) * 0.5 * PI;
   square.a = cos(angle) * scale;
   square.d = sin(angle) * scale;
   square.c = (double)rand() / RAND_MAX * 4.0;
   square.b = -square.d;
   square.e = square.a;
   square.f = (double)rand() / RAND_MAX * 2.0;

   /* Generate wave offsets for corners. */
   for(x = 0; x < 2; x++)
   {
      for(y = 0; y < 2; y++)
      {
         phase[x][y] = (int)((double)rand() / RAND_MAX * 8.0);
         period[x][y] = (int)((double)rand() / RAND_MAX * 6.0) + 2;
      }
   }

   /* 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, &square);
         p3.x = floor(p2.x);
         p3.y = floor(p2.y);

         /* If pixel is not near the corner, set parity based on
            square coordinates.                                  */
         if( !IsNearCorner(p2.x - p3.x, p2.y - p3.y) )
         {
            parity = (int)(p3.x + p3.y) % 2;
            fputc(parity ? 255 : 0, stdout);
            continue;
         }

         /* Determine which diagonal line the pixels are on. */
         p3.x = round(p2.x);
         p3.y = round(p2.y);
         d = (int)(p3.x + p3.y) % 2 != 0;

         /* Apply wave function specific to that diagonal line,
            such that parity alternates every few corners along that line. */
         parity = 0;
         for(i = 0; i < 2; i++)
         {
            w = (int)fmod(phase[i][d] + p3.x - i * p3.y, period[i][d]);
            if( w < 0 )
               w += period[i][d];
            parity ^= w > period[i][d] / 2 ? 1 : 0;
         }

         /* Render corner as a 2x2 checkerboard.

            This step can be skipped to debug the wave function bits.
            Doing so causes the corners to be colored with solid
            colors in a pattern that kind of resembles certain circuit
            layouts, which is also kind of neat, but not quite the
            effect we were going for.                                  */
         parity ^= (p2.x - p3.x > fabs(p2.y - p3.y) ||
                    p2.x - p3.x < -fabs(p2.y - p3.y)) ? 1 : 0;

         fputc(parity ? 255 : 0, stdout);
      }
   }
   return 0;
}