/* Trail path generator.

   ./trail <width> <height> > output.pgm
                                         */
#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 PI           3.14159265358979323846264338327950288419716939937510
#define EPSILON      1e-8
#define POINT_COUNT  16


typedef struct { double x, y; } XY;

/* Fisher-Yates shuffle. */
static void Shuffle(int size, int *elements)
{
   int i, j, tmp;

   for(i = size - 1; i > 1; i--)
   {
      j = (int)((double)rand() / RAND_MAX * i);
      tmp = elements[i];
      elements[i] = elements[j];
      elements[j] = tmp;
   }
}

/* Move a point toward another point. */
static void MoveToward(double max_velocity,
                       double acceleration,
                       XY target,
                       XY input_position,
                       XY *output_position,
                       XY *velocity)
{
   double dx = target.x - input_position.x;
   double dy = target.y - input_position.y;
   double d = hypot(dx, dy);

   if( d < EPSILON )
      d = EPSILON;
   dx *= acceleration / d;
   dy *= acceleration / d;

   velocity->x += dx;
   velocity->y += dy;
   d = hypot(velocity->x, velocity->y);
   if( d > max_velocity )
   {
      velocity->x *= max_velocity / d;
      velocity->y *= max_velocity / d;
   }
   output_position->x = input_position.x + velocity->x;
   output_position->y = input_position.y + velocity->y;
}

/* Compute location a point that is some distance along the segment from
   source to target.                                                     */
static void Project(XY source, XY target, double distance, XY *output)
{
   double dx = target.x - source.x;
   double dy = target.y - source.y;
   double d = hypot(dx, dy);

   if( d < EPSILON )
      d = EPSILON;
   output->x = source.x + dx * distance / d;
   output->y = source.y + dy * distance / d;
}

/* Minimum/maximum. */
static int Min(int a, int b) { return a < b ? a : b; }
static int Max(int a, int b) { return a > b ? a : b; }
static int Min3(int a, int b, int c) { return Min(Min(a, b), c); }
static int Max3(int a, int b, int c) { return Max(Max(a, b), c); }

/* Cross product of two vectors. */
static double CrossProduct(double ax, double ay, double bx, double by)
{
   return ax * by - bx * ay;
}

/* Barycentric algorithm for drawing triangles.

   This is slow in theory but fast enough for what we do. */
static void DrawTriangle(int width, int height,
                         XY a, XY b, XY c, unsigned char *pixels)
{
   int x, y;
   double c1, c2, c3;

   for(y = Max(0, Min3((int)a.y, (int)b.y, (int)c.y));
       y <= Min(height - 1, Max3((int)a.y, (int)b.y, (int)c.y));
       y++)
   {
      for(x = Max(0, Min3((int)a.x, (int)b.x, (int)c.x));
          x <= Min(width - 1, Max3((int)a.x, (int)b.x, (int)c.x));
          x++)
      {
         /* One possible optimization here is to skip all these cross
            product computation if the pixel is already covered.  That
            turned out to be a wash because the pixels don't overlap
            all that much in practice.                                 */
         c1 = CrossProduct(x - a.x, y - a.y, x - b.x, y - b.y);
         c2 = CrossProduct(x - b.x, y - b.y, x - c.x, y - c.y);
         c3 = CrossProduct(x - c.x, y - c.y, x - a.x, y - a.y);
         if( (c1 < 0 && c2 < 0 && c3 < 0) || (c1 > 0 && c2 > 0 && c3 > 0) )
            pixels[y * width + x] = 255;
      }
   }
}

int main(int argc, char **argv)
{
   int width, height, short_edge, frame, i, j;
   int visit_order[POINT_COUNT];
   int visited_quadrant[4], current_quadrant, previous_quadrant;
   unsigned char *pixels;
   XY target, position[2], extension[2], velocity, origin;
   double max_velocity, acceleration, a;

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

   if( (pixels = (unsigned char*)calloc(width * height, 1)) == NULL )
      return !puts("Not enough memory");

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

   srand(time(NULL));

   /* Initialize origin to be horizontally on either the left or right 1/3
      of the screen, and vertically somewhere in the middle 1/3.           */
   origin.x = rand() > RAND_MAX / 2 ? 2 * width / 3.0 : width / 3.0;
   origin.y = ((double)rand() / RAND_MAX + 1.0) / 3.0 * height;

   /* Initialize random position and velocity. */
   position[0].x = (double)rand() / RAND_MAX * width;
   position[0].y = (double)rand() / RAND_MAX * height;
   short_edge = width < height ? width : height;
   max_velocity = short_edge * 0.03;
   acceleration = max_velocity * 0.1;
   a = (double)rand() / RAND_MAX * PI * 2.0;
   velocity.x = max_velocity * cos(a);
   velocity.y = max_velocity * sin(a);

   /* Initialize first extension point. */
   Project(position[0], origin, max_velocity, &extension[0]);

   /* Randomize order to visit each area. */
   for(i = 0; i < POINT_COUNT; i++)
      visit_order[i] = i;
   Shuffle(POINT_COUNT, visit_order);

   /* Run simulation. */
   frame = 0;
   for(i = 0; i < POINT_COUNT; i++)
   {
      /* Set target to be the center of the area we want to visit. */
      j = visit_order[i];
      target.x = (j % 4) * width / 4 + width / 8;
      target.y = (j / 4) * height / 4 + height / 8;

      /* Reset quadrant visit counts. */
      memset(visited_quadrant, 0, sizeof(visited_quadrant));
      previous_quadrant = -1;

      /* Generate a trail of triangle strips while moving toward a
         target, but only do this for a finite number of steps before
         moving on to the next target.                                */
      for(j = 0; j < 1000; j++)
      {
         /* Stop early if we are already circling around the target,
            but couldn't get close enough.                           */
         current_quadrant = target.x > position[frame].x
            ? target.y > position[frame].y ? 0 : 1
            : target.y > position[frame].y ? 2 : 3;
         if( previous_quadrant != current_quadrant )
         {
            previous_quadrant = current_quadrant;
            ++visited_quadrant[current_quadrant];
            if( visited_quadrant[current_quadrant] > 2 )
               break;
         }

         /* Stop early if we got close enough to target. */
         if( hypot(target.x - position[frame].x,
                   target.y - position[frame].y) <= short_edge / 6 )
         {
            break;
         }

         /* Move position toward target. */
         MoveToward(max_velocity,
                    acceleration,
                    target,
                    position[frame],
                    &position[frame ^ 1],
                    &velocity);

         /* Project extension toward origin. */
         Project(position[frame],
                 origin,
                 max_velocity,
                 &extension[frame ^ 1]);

         /* Draw triangles.
            E[f^1]         P[f+1]
                 +---------+
                 |       / |
                 |     /   |
                 |   /     |
                 | /       |
                 +---------+
              E[f]         P[f]     */
         DrawTriangle(width, height,
                      position[frame],
                      position[frame ^ 1],
                      extension[frame],
                      pixels);
         DrawTriangle(width, height,
                      extension[frame],
                      extension[frame ^ 1],
                      position[frame ^ 1],
                      pixels);
         frame ^= 1;
      }
   }

   /* Write pixels. */
   printf("P5\n%d %d\n255\n", width, height);
   for(i = 0; i < width * height; i++)
      fputc((int)pixels[i], stdout);
   free(pixels);
   return 0;
}