/* Fractal pattern generator.

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

   This is an optimized version of fractal_slow.c.  The optimization
   works by pre-computing transformation matrices to reduce the number of
   multiplications at each point.

   Before:                             After:
   --------------------------------    -----------------------------------
                                       Pre-compute 3**8 rules
   Repeat a few million times:         Repeat a few million times:
      Pick a random starting point        Pick a random starting point
      Repeat for 16 iterations:           Repeat for 2 iterations:
         Select random rule among 3          Select random rule among 3**8
         Apply rule to point                 Apply rule to point
      Draw transformed point              Draw transformed point

   Thus we get a 8x speedup, at the cost of some extra memory, and the
   resulting output will be slightly different because we made fewer
   random selections, and maybe also some difference due to precision
   loss in chaining transformation matrices.  But 8x speedup, it's
   totally worth it.
                                                                           */

#include<assert.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#ifdef _WIN32
   #include<fcntl.h>
   #include<io.h>
#endif

/* Number of rules is fixed, because the way we generate randomized
   parameters is hardcoded to a specific number of rules.

   Honestly we already have enough trouble generating good looking
   fractals with the rule count being fixed.  Letting the rule count
   be variable would take a lot more testing.                        */
#define RULE_COUNT 3

/* Number of iterations is fixed.  It's tricky to dynamically figure
   out when the coordinates would converge, so we just run a fixed
   number of steps.

   Higher number of iterations tend to lead to fractals that are
   more detailed because rules that causes expansion gets to expand
   deeper, but after some threshold it's mostly just noise.

   This must be an integer multiple of SUB_ITERATIONS, below.       */
#define ITERATIONS 16

/* Number of iterations to expand in a single step.

   That is, instead of applying the transformations like this:
      repeat 16 times:
         select 1 random rule out of 3
         apply rule

   We do this:
      repeat 2 times:
         select 1 random transformation out of 3**8
         apply transformation

   This means we do way fewer computations per point, at the cost of
   having to pre-compute and store 3**8 transformations.

   Note that the "8" iteration count is important, not just because it
   saves memory, but also because on systems where rand() only returns
   15 bits (e.g. MingW), setting an iteration count at 10 or above
   would cause some of the transformations to be unused.               */
#define SUB_ITERATIONS 8

/* Maximum number of points to land on the same pixel before concluding
   that we are getting diminishing returns from trying to draw more points.

   The optimal value was empirically determined to be around 80.  Lower
   values noticeably increases the amount of noise in output, and higher
   values only offers miniscule improvement in noise.                      */
#define MAX_OVERLAP 80

/* Top 6 elements of a transformation matrix.
   [ x[0] x[1] x[2]
     x[3] x[4] x[5]
     0.0  0.0  1.0  ]                         */
typedef struct { double x[6]; } Transform;

/* Compute a*b, storing result in b. */
static void ApplyTransformation(const Transform *a, Transform *b)
{
   const Transform r =
   {
      {
         a->x[0] * b->x[0] + a->x[1] * b->x[3],
         a->x[0] * b->x[1] + a->x[1] * b->x[4],
         a->x[0] * b->x[2] + a->x[1] * b->x[5] + a->x[2],

         a->x[3] * b->x[0] + a->x[4] * b->x[3],
         a->x[3] * b->x[1] + a->x[4] * b->x[4],
         a->x[3] * b->x[2] + a->x[4] * b->x[5] + a->x[5]
      }
   };
   memcpy(b, &r, sizeof(r));
}

/* Compute a*point. */
static void ApplyTransformationToPoint(const Transform *a, double *x, double *y)
{
   const double nx = a->x[0] * *x + a->x[1] * *y + a->x[2];
   const double ny = a->x[3] * *x + a->x[4] * *y + a->x[5];
   *x = nx;
   *y = ny;
}

/* Recursively compute all RULE_COUNT**SUB_ITERATIONS transformations. */
static void GenerateTransformations(int depth,
                                    const Transform *rule,
                                    const Transform *parent,
                                    Transform **output)
{
   Transform child;
   int i;

   if( depth == SUB_ITERATIONS )
   {
      memcpy(*output, parent, sizeof(*parent));
      ++*output;
      return;
   }

   for(i = 0; i < RULE_COUNT; i++)
   {
      memcpy(&child, parent, sizeof(child));
      ApplyTransformation(&rule[i], &child);
      GenerateTransformations(depth + 1, rule, &child, output);
   }
}

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

/* Compute area of parallelogram from 3 corner vertices. */
static double GetParallelogramArea(double ax, double ay,
                                   double bx, double by,
                                   double cx, double cy)
{
   return fabs(CrossProduct(bx - ax, by - ay, cx - ax, cy - ay));
}

/* Compute transformation, such that when applied to a right triangle
   with corners at (0,0), (1,0), and (0,1), the transformed triangle
   would touch the points A+B+C.

   F(0,0) -> A
   F(1,0) -> B
   F(0,1) -> C

   X(x,y) = output->x[0] * x + output->x[1] * y + output->x[2]
   Y(x,y) = output->x[3] * x + output->x[4] * y + output->x[5]        */
static void GetTransform(double ax, double ay,
                         double bx, double by,
                         double cx, double cy,
                         Transform *output)
{
   /* (0,0) -> A */
   output->x[2] = ax;
   output->x[5] = ay;

   /* (1,0) -> B */
   output->x[0] = bx - ax;
   output->x[3] = by - ay;

   /* (0,1) -> C */
   output->x[1] = cx - ax;
   output->x[4] = cy - ay;
}

int main(int argc, char **argv)
{
   int width, height, points, out_of_bounds, transform_count, i, j, k;
   unsigned char *pixels;
   double vertex[3][2], fx, fy;
   Transform rule[RULE_COUNT], *transforms, *t;
   const Transform identity = {{1, 0, 0,  0, 1, 0}};

   assert(pow(RULE_COUNT, SUB_ITERATIONS) <= RAND_MAX);

   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");
   transform_count = (int)pow(RULE_COUNT, SUB_ITERATIONS);
   transforms = (Transform*)malloc(transform_count * sizeof(Transform));
   if( transforms == NULL )
      return printf("Not enough memory for %d transforms\n", transform_count);

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

   srand(time(NULL));

   /* Generate triangle vertices until we have one that covers a
      sufficient area.  The area requirement here is to ensure that the
      vertices are sufficiently spread apart.

      This is potentially an infinite loop because we will keep generating
      random numbers until the area requirement is satisfied.  But assuming
      that the numbers we get from rand() is truly random, the passing rate
      should be ~0.006 according to fractal_area_stats.c, so we will take
      our chances.                                                          */
   do
   {
      for(i = 0; i < 3; i++)
      {
         for(j = 0; j < 2; j++)
            vertex[i][j] = 0.8 * (double)rand() / RAND_MAX + 0.1;
      }
      fx = GetParallelogramArea(vertex[0][0], vertex[0][1],
                                vertex[1][0], vertex[1][1],
                                vertex[2][0], vertex[2][1]);
   } while( fx < 0.4 );

   /* Pick three random points near each of those vertices.

      Also enforce a range check on areas covered by these parallelograms.
      We don't want the area to be too small because that means
      fractal points will converge toward thin lines, and we don't
      want the area to be too large because that will end up covering
      most of the canvas with noise.

      An area range between [0.3 .. 0.7] appears to produce the best
      results, with a good mix of pixels and blank space.  Lowering the
      thresholds resulted in too much blank space, and raising the
      thresholds causes most of the canvas to be covered by noise.

      Passing rate for the area test is ~0.14 according to
      fractal_area_stats.c.                                               */
   for(i = 0; i < 3; i++)
   {
      do
      {
         for(j = 0; j < 6; j++)
            rule[i].x[j] = vertex[i][j % 2] + (double)rand() / RAND_MAX - 0.5;
         fx = GetParallelogramArea(rule[i].x[0], rule[i].x[1],
                                   rule[i].x[2], rule[i].x[3],
                                   rule[i].x[4], rule[i].x[5]);
      } while( fx < 0.3 || fx > 0.7 );
      GetTransform(rule[i].x[0], rule[i].x[1],
                   rule[i].x[2], rule[i].x[3],
                   rule[i].x[4], rule[i].x[5],
                   &rule[i]);
   }

   /* Compute all transforms. */
   t = transforms;
   GenerateTransformations(0, rule, &identity, &t);

   /* Draw points.

      We make a decision on whether adding more points is seeing
      diminishing returns based on where it lands:

      - If points are landing on top of pixels that have already been
        covered previously, stop after a pixel has gotten enough
        overlap.

      - Otherwise, stop after we saw enough points land outside of
        image area.                                                   */
   for(points = out_of_bounds = 0;
       out_of_bounds < width * height * MAX_OVERLAP;)
   {
      fx = (double)rand() / RAND_MAX;
      fy = (double)rand() / RAND_MAX;
      for(i = 0; i < ITERATIONS / SUB_ITERATIONS; i++)
      {
         k = (int)((double)rand() / (RAND_MAX + 1u) * transform_count);
         ApplyTransformationToPoint(&transforms[k], &fx, &fy);
      }

      /* Convert fractal coordinates back to integer image
         coordinates.  Note the use of round() here, which avoids
         truncation artifacts near the top and left edges.        */
      i = (int)round(fx * width);
      j = (int)round(fy * height);
      if( i < 0 || j < 0 || i >= width || j >= height )
      {
         out_of_bounds++;
         continue;
      }
      k = j * width + i;
      if( pixels[k] >= MAX_OVERLAP )
         break;
      pixels[k]++;
      points++;
   }

   fprintf(stderr, "width = %d, height = %d, points = %d, out_of_bounds = %d\n",
           width, height, points, out_of_bounds);

   /* Normalize and write pixels. */
   j = 1;
   for(i = 0; i < width * height; i++)
      j = j < (int)pixels[i] ? (int)pixels[i] : j;
   printf("P5\n%d %d\n255\n", width, height);
   for(i = 0; i < width * height; i++)
      fputc((int)pixels[i] * 255 / j, stdout);
   free(pixels);
   free(transforms);
   return 0;
}