/* Fractal pattern generator.

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

   https://en.wikipedia.org/wiki/Iterated_function_system
                                                          */
#include<math.h>
#include<stdio.h>
#include<stdlib.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.         */
#define ITERATIONS 16

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

/* 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 coefficients, 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[0] * x + output[2] * y + output[4]
   Y(x,y) = output[1] * x + output[3] * y + output[5]
                                                                  */
static void GetCoefficients(double ax, double ay,
                            double bx, double by,
                            double cx, double cy,
                            double *output)
{
   /* (0,0) -> A */
   output[4] = ax;
   output[5] = ay;

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

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

int main(int argc, char **argv)
{
   int width, height, points, out_of_bounds, i, j, k;
   unsigned char *pixels;
   double vertex[3][2], rule[RULE_COUNT][6], fx, fy, nx, ny;

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

   /* 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][j] = vertex[i][j % 2] + (double)rand() / RAND_MAX - 0.5;
         fx = GetParallelogramArea(rule[i][0], rule[i][1],
                                   rule[i][2], rule[i][3],
                                   rule[i][4], rule[i][5]);
      } while( fx < 0.3 || fx > 0.7 );
      GetCoefficients(rule[i][0], rule[i][1],
                      rule[i][2], rule[i][3],
                      rule[i][4], rule[i][5],
                      rule[i]);
   }

   /* 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; i++)
      {
         k = (int)((double)rand() / (RAND_MAX + 1u) * RULE_COUNT);
         nx = fx * rule[k][0] + fy * rule[k][2] + rule[k][4];
         ny = fx * rule[k][1] + fy * rule[k][3] + rule[k][5];
         fx = nx;
         fy = ny;
      }

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