/* Fractal pattern generator.

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

   This is a fork of fractal.c that generates rules in a more
   constrained way.
                                                              */

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

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

   /* For each of those vertices, generate 3 vertices that are axis-aligned.
      This means only translation and scaling are applied by each rule, no
      rotation or skew.

      Usually what we get from this are variants of Sierpinski triangle.
                                                                            */
   for(i = 0; i < 3; i++)
   {
      fx = (double)rand() / RAND_MAX * 0.6 + 0.1 * (i + 2);
      fy = (double)rand() / RAND_MAX * 0.6 + 0.1 * (i + 2);
      rule[i].x[0] = fx;
      rule[i].x[4] = fy;
      rule[i].x[1] = rule[i].x[3] = 0;
      rule[i].x[2] = vertex[i][0] - fx * 0.5;
      rule[i].x[5] = vertex[i][1] - fy * 0.5;
   }

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