// generate Barnsley fern
// refs:
// https://en.wikipedia.org/wiki/Barnsley_fern
// http://www.home.aone.net.au/~byzantium/ferns/fractal.html
//
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "ppm.h"

typedef struct { double x, y; } Point;

#define WIDTH 400 // image width

#define ITER 1000000 // number of iterations

Point P[ITER]; // all of the generated points

double params[4][7] = {
//    0,	1,	2,	3,	4,	5,	6
//    a,	b,	c,	d,	e,	f,	p
#if 0 //  Barnsley
    { 0,	0,	0,	0.16,	0,	0,	0.01 },
    { 0.85,	0.04,	-0.04,	0.85,	0,	1.60,	0.85 },
    { 0.20,	-0.26,	0.23,	0.22,	0,	1.60,	0.07 },
    { -0.15,	0.28,	0.26,	0.24,	0,	0.44,	0.07 }
#else
#if 0 // alternate 1
    { 0,	0,	0,	0.25,	0,	-0.4,	0.02 },
    { 0.95,	0.005,	-0.005,	0.93,	-0.002,	0.5,	0.84 },
    { 0.035,	-0.2,	0.16,	0.04,	-0.09,	0.02,	0.07 },
    { -0.04,	0.2,	0.16,	0.04,	0.083,	0.12,	0.07  },
#else // alternate 2
    { 0,	0,	0,	0.25,	0,	-0.14,	0.02 },
    { 0.85,	0.02,	-0.02,	0.83,	0,	1,	0.84 },
    { 0.09,	-0.28,	0.3,	0.11,	0,	0.6,	0.07 },
    { -0.09,	0.28,	0.3,	0.09,	0,	0.7,	0.07 },
#endif
#endif
};

// create next point
//
Point f( Point p, double g [])
{
  Point r = { g[0]*p.x + g[1]*p.y + g[4], g[2]*p.x + g[3]*p.y + g[5] };

  return r;
}

int main( int argc, char *argv[])
{
  srand(time(0));

 // make params probability cumulative
 //
  for( int i = 1; i < 4; ++i) params[i][6] += params[i-1][6];

#if 0 // debug
  for( int i = 0; i < 4; ++i) {
    for( int j = 0; j < 7; ++j) printf(" %g", params[i][j]);
    putchar( '\n');
  }
#endif


  int k; double xmin, xmax, ymin, ymax, r;

  xmin = xmax = ymin = ymax = P[0].x = P[0].y = 0;

  for( int i = 1; i < ITER; ++i)
  {
    r = rand()/(RAND_MAX+1.0);

    k = 3; for( int i = 0; i < 3; ++i) if( r <= params[i][6]) { k = i; break; }

    P[i] = f( P[i-1], params[k]);

    xmin = fmin( xmin, P[i].x);
    xmax = fmax( xmax, P[i].x);
    ymin = fmin( ymin, P[i].y);
    ymax = fmax( ymax, P[i].y);
  }

  // printf( "%g %g %g %g\n", xmin, xmax, ymin, ymax);

  double dx = xmax - xmin, dy = ymax - ymin;

  ppm a = ppm_new( WIDTH*(dy/dx), WIDTH);

  memset( a.r[0], 255, a.rows*a.cols);
  memset( a.g[0], 255, a.rows*a.cols);
  memset( a.b[0], 255, a.rows*a.cols);

  double xscale = dx/(a.cols-1), yscale = dy/(a.rows-1);

  for( int i = 1; i < ITER; ++i)
  {
    int row = (a.rows-1) - (P[i].y - ymin)/yscale; // flip vertically
    int col = (P[i].x - xmin)/xscale;

    a.r[row][col] = a.b[row][col] = 0;
  }

  ppm_write(a);
}