#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<time.h>
#include<png.h>

/* Number of fractal rules. */
#define RULE_COUNT 3

/* Number of fractal iterations. */
#define ITERATIONS 16

/* Number of sub-iterations to apply in one step.

   Normal fractal works like this:
   1. Select a rule at random, out of rule_count choices.
   2. Apply rule.
   3. Repeat steps 1-2 multiple times.

   We can speed this up like this:
   0. Pre-compute explosion of rules_count**substep rules,
      i.e. all possible rules up to substep count.
   1. Select one rule out of rule_count**substep choices.
   2. Apply rule.
   3. Repeat steps 1-2 multiple times.

   This is faster because the number of iterations at step 3 only need
   to be ITERATION_COUNT / SUBSTEP_COUNT, so we reduce the number of
   transformations per input pixel.                                    */
#define SUBSTEP_COUNT 8

/* Maximum amount of overlap on the same pixel. */
#define MAX_OVERLAP 80

/* Base vertices for rules.  The transforms will be centered around
   these vertices.                                                  */
static double vertex[RULE_COUNT][2];

/* Transformation rule.  Basically 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;

/* Transformation rules. */
Transform base_rule[RULE_COUNT];
Transform rule[6561];   /* 6561 = 3 ** 8 */
Transform identity = {{1, 0, 0,  0, 1, 0}};
Transform *t = rule;

/* Count of each pixel value. */
static int color_count[MAX_OVERLAP + 1];

/* Apply transformation a to b, storing result in b. */
static void ApplyTransformation(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));
}

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

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

   for(i = 0; i < RULE_COUNT; i++)
   {
      memcpy(&child, parent, sizeof(child));
      ApplyTransformation(&base_rule[i], &child);
      GenerateTransformations(depth + 1, &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));
}

/* Convert vertices to transformation matrix.

   Assume that input specifies 3 vertices of a triangle, we want a
   transformation matrix such that

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

   Transformations are applied as:
   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 ConvertToTransform(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)
{
   png_image image;
   png_bytep buffer1, buffer2;
   double a, fx, fy, nx, ny;
   int s, i, j, k, out_of_bounds, *pixel;

   if( argc != 4 )
      return printf("%s <input.png> <output1.png> <output2.png>\n", *argv);

   memset(&image, 0, sizeof(image));
   image.version = PNG_IMAGE_VERSION;
   if( !png_image_begin_read_from_file(&image, argv[1]) )
      return printf("Error reading %s\n", argv[1]);

   image.format = PNG_FORMAT_RGBA;
   s = image.width * image.height;
   buffer1 = (png_bytep)malloc(s * 4);
   buffer2 = (png_bytep)calloc(s, 4);
   if( buffer1 == NULL || buffer2 == NULL )
      return puts("Out of memory");
   if( !png_image_finish_read(&image, NULL, buffer1, 0, NULL) )
      return printf("Error loading %s\n", argv[1]);

   srand(time(NULL));

   /* Initialize vertices in [0.1,0.9] range with sufficient spread. */
   do
   {
      for(i = 0; i < RULE_COUNT; i++)
      {
         for(j = 0; j < 2; j++)
            vertex[i][j] = 0.8 * (double)rand() / RAND_MAX + 0.1;
      }
      a = GetParallelogramArea(vertex[0][0], vertex[0][1],
                               vertex[1][0], vertex[1][1],
                               vertex[2][0], vertex[2][1]);
   } while( a < 0.4 );

   /* Pick random points near the generated vertices.  These will be
      used to initialize transformation rules.

      To ensure that each of these transformations produce interesting
      results, we also require that the parallelograms formed by the
      vertices cover sufficient area.                                  */
   for(i = 0; i < RULE_COUNT; i++)
   {
      if( sizeof('c') > 1 )
      {
         /* C mode: generate constrained axis-aligned rules. */
         fx = (double)rand() / RAND_MAX * 0.6 + 0.1 * (i + 2);
         fy = (double)rand() / RAND_MAX * 0.6 + 0.1 * (i + 2);
         base_rule[i].x[0] = fx;
         base_rule[i].x[4] = fy;
         base_rule[i].x[1] = base_rule[i].x[3] = 0.0;
         base_rule[i].x[2] = vertex[i][0] - fx * 0.5;
         base_rule[i].x[5] = vertex[i][1] - fy * 0.5;
      }
      else
      {
         /* C++ mode: generate free range random transforms that
            falls within expected area range.                    */
         do
         {
            for(j = 0; j < 6; j++)
            {
               base_rule[i].x[j] =
                  vertex[i][j % 2] + (double)rand() / RAND_MAX - 0.5;
            }
            a = GetParallelogramArea(base_rule[i].x[0], base_rule[i].x[1],
                                     base_rule[i].x[2], base_rule[i].x[3],
                                     base_rule[i].x[4], base_rule[i].x[5]);
         } while( a < 0.3 || a > 0.7 );

         /* Convert vertices to transformation matrix. */
         ConvertToTransform(base_rule[i].x[0], base_rule[i].x[1],
                            base_rule[i].x[2], base_rule[i].x[3],
                            base_rule[i].x[4], base_rule[i].x[5],
                            &base_rule[i]);
      }
   }

   /* Pre-compute all transformations. */
   GenerateTransformations(0, &identity, &t);

   /* Draw points.

      This is done by applying random rules iteratively to each starting
      point, and repeat the process until we appear to be getting
      diminishing returns, which is one of these conditions:
      - Many points are landing on pixels that were previously drawn.
      - Many points are drawn out of bounds.                             */
   pixel = (int*)buffer2;
   for(out_of_bounds = 0; out_of_bounds < s * MAX_OVERLAP;)
   {
      fx = (double)rand() / RAND_MAX;
      fy = (double)rand() / RAND_MAX;
      for(i = 0; i < ITERATIONS / SUBSTEP_COUNT; i++)
      {
         k = (int)((double)rand() / (RAND_MAX + 1u) * 6561);
         nx = fx * rule[k].x[0] + fy * rule[k].x[1] + rule[k].x[2];
         ny = fx * rule[k].x[3] + fy * rule[k].x[4] + rule[k].x[5];
         fx = nx;
         fy = ny;
      }

      /* Convert fractal coordinates from [0,1] range back to screen range. */
      i = (int)round(fx * image.width);
      j = (int)round(fy * image.height);
      if( i < 0 || j < 0 || i >= (int)image.width || j >= (int)image.height )
      {
         out_of_bounds++;
         continue;
      }
      k = j * image.width + i;
      if( pixel[k] >= MAX_OVERLAP )
         break;  /* Too many pixels landed on the same point. */
      pixel[k]++;
   }

   /* Count number of colors, and pick a threshold such that number of
      pixels above and below the threshold are roughly equal.          */
   memset(color_count, 0, sizeof(color_count));
   for(i = 0; i < s; i++)
      color_count[pixel[i]]++;
   j = 0;
   for(k = MAX_OVERLAP; k > 1; k--)
   {
      j += color_count[k];
      if( j >= s / 2 )
         break;
   }

   /* Move pixels according to threshold. */
   for(i = 0; i < s; i++)
   {
      if( pixel[i] >= k )
      {
         memcpy(buffer2 + (i * 4), buffer1 + (i * 4), 4);
         memset(buffer1 + (i * 4), 0, 4);
      }
      else
      {
         pixel[i] = 0;
      }
   }

   if( !png_image_write_to_file(&image, argv[2], 0, buffer1, 0, NULL) )
      return printf("Error writing %s\n", argv[2]);
   if( !png_image_write_to_file(&image, argv[3], 0, buffer2, 0, NULL) )
      return printf("Error writing %s\n", argv[3]);

   return 0;
}