/* Trail path generator.

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

   This is forked from trail.c to produce more than one trail.
                                                               */
#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
#define TRAIL_COUNT  3


typedef struct { double x, y; } XY;

typedef struct
{
   /* Order for visiting each of the target points. */
   int visit_order[POINT_COUNT];

   /* Count of number of quadrants around each target point
      that was visited, used for counting loops.            */
   int visited_quadrant[4];

   /* Quadrant indices. */
   int current_quadrant;
   int previous_quadrant;

   /* Trail origin. */
   XY origin;

   /* Trail state. */
   XY position[2], extension[2], velocity, target;

   /* Simulation step counters. */
   int outer_step;
   int inner_step;
} Trail;

/* 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, unsigned char color,
                         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] = color;
      }
   }
}

/* Check if we need to advance the outer loop counter, return 1 if so. */
static int DoneWithInnerLoop(int short_edge, int frame, Trail *trail)
{
   /* Stop if we ran enough steps.

      We will always trip over this condition on first iteration. */
   if( trail->inner_step > 1000 )
      return 1;

   /* Stop early if we are already circling around the same target,
      but couldn't get close enough.                                */
   trail->current_quadrant = trail->target.x > trail->position[frame].x
      ?  trail->target.y > trail->position[frame].y ? 0 : 1
      :  trail->target.y > trail->position[frame].y ? 2 : 3;
   if( trail->previous_quadrant != trail->current_quadrant )
   {
      trail->previous_quadrant = trail->current_quadrant;
      trail->visited_quadrant[trail->current_quadrant]++;
      if( trail->visited_quadrant[trail->current_quadrant] > 2 )
         return 1;
   }

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

int main(int argc, char **argv)
{
   int width, height, short_edge, frame, completed_trails, i, j, t;
   unsigned char *pixels;
   double max_velocity, acceleration, a;
   Trail trail[TRAIL_COUNT];

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

   short_edge = width < height ? width : height;
   max_velocity = short_edge * 0.03;
   acceleration = max_velocity * 0.1;

   /* Initial trail states. */
   for(i = 0; i < TRAIL_COUNT; i++)
   {
      /* Randomize visit order. */
      for(j = 0; j < POINT_COUNT; j++)
         trail[i].visit_order[j] = j;
      Shuffle(POINT_COUNT, trail[i].visit_order);

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

      /* Initialize random position and velocity. */
      trail[i].position[0].x = (double)rand() / RAND_MAX * width;
      trail[i].position[0].y = (double)rand() / RAND_MAX * height;
      a = (double)rand() / RAND_MAX * PI * 2.0;
      trail[i].velocity.x = max_velocity * cos(a);
      trail[i].velocity.y = max_velocity * sin(a);

      /* Initialize first extension point. */
      Project(trail[i].position[0],
              trail[i].origin,
              max_velocity,
              &trail[i].extension[0]);
   
      /* Set inner step counter to overflow such that outer loop
         states will be initialized on first iteration.          */
      trail[i].outer_step = -1;
      trail[i].inner_step = 0x7fffffff;
   }

   /* Run simulation for all trails in lockstep.
      This is so that the trail pixels would interleave. */
   frame = 0;
   do
   {
      completed_trails = 0;
      for(t = 0; t < TRAIL_COUNT; t++)
      {
         if( trail[t].outer_step >= POINT_COUNT )
         {
            completed_trails++;
            continue;
         }
         if( DoneWithInnerLoop(short_edge, frame, &trail[t]) )
         {
            trail[t].outer_step++;
            if( trail[t].outer_step >= POINT_COUNT )
            {
               completed_trails++;
               continue;
            }
            trail[t].inner_step = 0;

            /* Set target to be the center of the area we want to visit. */
            i = trail[t].visit_order[trail[t].outer_step];
            trail[t].target.x = (i % 4) * width / 4 + width / 8;
            trail[t].target.y = (i / 4) * height / 4 + height / 8;

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

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

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

         /* Draw triangles.
            E[f^1]         P[f+1]
                 +---------+
                 |       / |
                 |     /   |
                 |   /     |
                 | /       |
                 +---------+
              E[f]         P[f]     */
         DrawTriangle(width, height,
                      (unsigned char)((t + 1) * 255 / TRAIL_COUNT),
                      trail[t].position[frame],
                      trail[t].position[frame ^ 1],
                      trail[t].extension[frame],
                      pixels);
         DrawTriangle(width, height,
                      (unsigned char)((t + 1) * 255 / TRAIL_COUNT),
                      trail[t].extension[frame],
                      trail[t].extension[frame ^ 1],
                      trail[t].position[frame ^ 1],
                      pixels);
      }
      frame ^= 1;
   } while( completed_trails < TRAIL_COUNT );

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