#include <stdio.h>
#include <time.h>
#include <math.h>

#define OUTFILE "mandelbrot.bmp"

////////////////////////////////////////////////////////////////////////////////////

#define WIDTH 1024
#define HEIGHT 768

#define CENTRE_X 0.282
#define CENTRE_Y 0.01
#define ZOOM 150000

#define ITERATIONS 7000

// Plotting functions and parameters...

#define bailoutr(n) bailoutsqrt(n, 12)   // Try instead: bailoutlinear(n, 1.5)
#define bailoutg(n) bailoutsqrt(n, 3)  // Try instead: bailoutlinear(n, 0.35)
#define bailoutb(n) 0

// Subsamples, where we look "inside" each pixel for values and average them.
// Useful for despeckling the image, may "seem" to blur it slightly though...
// These values must be at least one.

#define SUBSAMPLES_X 2
#define SUBSAMPLES_Y 2

// Colours for the set itself...

#define IN_SET_R 0
#define IN_SET_G 0
#define IN_SET_B 0

////////////////////////////////////////////////////////////////////////////////////

struct imagemap {
   int width;
   int height;
   unsigned char * pixels;     // this will point to a malloc of unsigned chars.
   };

void setrgb (struct imagemap * bitmapstruct, int x, int y, int red, int green, int blue);
void cls (struct imagemap * bitmapstruct, int red, int green, int blue);
void drawbmp (struct imagemap * bitmapstruct, char * filename);
int bailoutlinear(int i, float multiplier);
int bailoutsqrt(int i, float multiplier);
int bailoutlog(int i, float multiplier);
int bailoutsquared(int i, float multiplier);
void init_imagemap(struct imagemap * bitmapstruct, int width, int height);
int getr(struct imagemap * bitmapstruct, int x, int y);
int getg(struct imagemap * bitmapstruct, int x, int y);
int getb(struct imagemap * bitmapstruct, int x, int y);

////////////////////////////////////////////////////////////////////////////////////

int main (void) 
{

double x; double r; double nextr;
double y; double s; double nexts;
int n; int xtoplot; int ytoplot;
double startx; double endx;
double starty; double endy;
double dx; double dy;
double dx_over_width;
int red, green, blue;
int tmpred, tmpgreen, tmpblue;
int subxnumber, subynumber;
double subx, suby;

struct imagemap bitmap;

init_imagemap(&bitmap, WIDTH, HEIGHT);
cls(&bitmap, IN_SET_R, IN_SET_G, IN_SET_B);

startx = CENTRE_X - ((double) WIDTH / (ZOOM * 2));
endx = CENTRE_X + ((double) WIDTH / (ZOOM * 2));

starty = CENTRE_Y - ((double) HEIGHT / (ZOOM * 2));
endy = CENTRE_Y + ((double) HEIGHT / (ZOOM * 2));

printf("\n Mandelbrot...\n From: %.10f, %.10f\n   To: %.10f, %.10f\n\n", startx, starty, endx, endy);
printf(" 0%%       20%%       40%%       60%%       80%%       100%%\n ");

dx = endx - startx;
dy = endy - starty;
dx_over_width = dx / WIDTH;

/*

- The actual screen pixel being plotted is in variables xtoplot and ytoplot.

- Coordinates we're Mandelbrotting are in variables x and y.
- We are interested in coordinates from (startx, starty) to (endx, endy).

- What distance a pixel covers has been calculated above and is dx_over_width.

- We do subsampling - for each pixel, "subpixels" are examined, and the colour is averaged.
- subxnumber and subynumber keep track of which subsample we are on,
    ie the first (which is 0), second (1) etc.
- subx and suby store the actual coordinates of the subpixel, for the Mandelbrot iteration.
    These will be x or y plus (the distance a subpixel covers, multiplied by subxnumber or subynumber).

*/

x = startx;
for (xtoplot = 0; xtoplot < WIDTH; xtoplot++)
{
   if (xtoplot % (WIDTH / 50) == 0)
   {
      printf(".");
      fflush(stdout);
   }
   
   y = starty;
   for (ytoplot = 0; ytoplot < HEIGHT; ytoplot++)
   {
   
#ifndef NO_OPTIM
   
      // To optimise a bit, we skip past bits in the main cardioid + it's major bulb
      // which are definitely in the set...

      while 
      (
         (x >  -1.2 && x <=  -1.1 && y >  -0.1 && y < 0.1)
      || (x >  -1.1 && x <=  -0.9 && y >  -0.2 && y < 0.2)
      || (x >  -0.9 && x <=  -0.8 && y >  -0.1 && y < 0.1)
      || (x > -0.69 && x <= -0.61 && y >  -0.2 && y < 0.2)
      || (x > -0.61 && x <=  -0.5 && y > -0.37 && y < 0.37)
      || (x >  -0.5 && x <= -0.39 && y > -0.48 && y < 0.48)
      || (x > -0.39 && x <=  0.14 && y > -0.55 && y < 0.55)
      || (x >  0.14 && x <   0.29 && y > -0.42 && y < -0.07)
      || (x >  0.14 && x <   0.29 && y >  0.07 && y < 0.42)
      )
      {
         // Uncomment to see optimised areas...
         // setrgb(&bitmap, xtoplot, ytoplot, 255, 255, 255, GFX_REPLACE);
         
         ytoplot++;
         y += dx_over_width;
         
         if (ytoplot >= HEIGHT - 1)
         {
            break;
         }
      }
      
#endif
      
      red = 0;
      green = 0;
      blue = 0;
      
      for (subxnumber = 0; subxnumber < SUBSAMPLES_X; subxnumber++)
      {
         subx = x + (subxnumber * dx_over_width / SUBSAMPLES_X);
         
         for (subynumber = 0; subynumber < SUBSAMPLES_Y; subynumber++)
         {
            suby = y + (subynumber * dx_over_width / SUBSAMPLES_Y);

            r = subx; s = suby;       // r = 0; s = 0;  also works (just adds an iteration)
            for (n = 0; n <= ITERATIONS; n++)
            {
               nextr = ((r * r) - (s * s)) + subx;
               nexts = (2 * r * s) + suby;
               r = nextr;
               s = nexts;
               
               if (n == ITERATIONS)
               {
                  red += IN_SET_R; green += IN_SET_G; blue += IN_SET_B;
                  break;
               } else if ((r * r) + (s * s) > 4) {
                  tmpred = bailoutr(n); tmpgreen = bailoutg(n); tmpblue = bailoutb(n);
                  
                  // These lines kinda alter the effect of the subsampling... arguably nicer without...
                  // if (tmpred > 255) tmpred = 255;
                  // if (tmpgreen > 255) tmpgreen = 255;
                  // if (tmpblue > 255) tmpblue = 255;
                  
                  red += tmpred; green += tmpgreen; blue += tmpblue;
                  break;
               }
            }
         }
      }
      red = red / (SUBSAMPLES_X * SUBSAMPLES_Y);
      green = green / (SUBSAMPLES_X * SUBSAMPLES_Y);
      blue = blue / (SUBSAMPLES_X * SUBSAMPLES_Y);
      setrgb(&bitmap, xtoplot, ytoplot, red, green, blue);
      
      y += dx_over_width;
   }
   x += dx_over_width;
}

drawbmp(&bitmap, OUTFILE);
printf("\n\n Saved to %s. Time elapsed: %d seconds.\n", OUTFILE, clock() / CLOCKS_PER_SEC);

return 0;
}

void cls (struct imagemap * bitmapstruct, int red, int green, int blue)
{
   int x;
   int y;
   
   for (x = 0; x < bitmapstruct->width; x++)
   {
      for (y = 0; y < bitmapstruct->height; y++)
      {
         setrgb(bitmapstruct, x, y, red, green, blue);
      }
   }
   return;
}

void setrgb (struct imagemap * bitmapstruct, int x, int y, int red, int green, int blue)
{
   if (x < 0 || x >= bitmapstruct->width || y < 0 || y >= bitmapstruct->height) return;
   if (red > 255) red = 255;
   if (red < 0) red = 0;
   if (green > 255) green = 255;
   if (green < 0) green = 0;
   if (blue > 255) blue = 255;
   if (blue < 0) blue = 0;
   
   bitmapstruct->pixels[(x * 3) + 0 + (y * bitmapstruct->width * 3)] = red;
   bitmapstruct->pixels[(x * 3) + 1 + (y * bitmapstruct->width * 3)] = green;
   bitmapstruct->pixels[(x * 3) + 2 + (y * bitmapstruct->width * 3)] = blue;
   
   return;
}

void drawbmp (struct imagemap * bitmapstruct, char * filename) {

   unsigned int headers[13];
   FILE * outfile;
   int extrabytes;
   int paddedsize;
   int x; int y; int n;

   extrabytes = 4 - ((bitmapstruct->width * 3) % 4);   // How many bytes of padding to add to each
                                                       // horizontal line - the size of which must
                                                       // be a multiple of 4 bytes.
   if (extrabytes == 4) extrabytes = 0;

   paddedsize = ((bitmapstruct->width * 3) + extrabytes) * bitmapstruct->height;

   // Headers...
   // Note that the "BM" identifier in bytes 0 and 1 is NOT included in these "headers".
                     
   headers[0]  = paddedsize + 54;      // bfSize (whole file size)
   headers[1]  = 0;                    // bfReserved (both)
   headers[2]  = 54;                   // bfOffbits
   headers[3]  = 40;                   // biSize
   headers[4]  = bitmapstruct->width;  // biWidth
   headers[5]  = bitmapstruct->height; // biHeight

   // Would have biPlanes and biBitCount in position 6, but they're shorts.
   // It's easier to write them out separately (see below) than pretend
   // they're a single int, especially with endian issues...

   headers[7]  = 0;                    // biCompression
   headers[8]  = paddedsize;           // biSizeImage
   headers[9]  = 0;                    // biXPelsPerMeter
   headers[10] = 0;                    // biYPelsPerMeter
   headers[11] = 0;                    // biClrUsed
   headers[12] = 0;                    // biClrImportant

   outfile = fopen(filename, "wb");

   //
   // Headers begin...
   // When printing ints and shorts, we write out 1 character at a time to avoid endian issues.
   //

   fprintf(outfile, "BM");

   for (n = 0; n <= 5; n++)
   {
      fprintf(outfile, "%c", headers[n] & 0x000000FF);
      fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
      fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
      fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
   }

   // These next 4 characters are for the biPlanes and biBitCount fields.

   fprintf(outfile, "%c", 1);
   fprintf(outfile, "%c", 0);
   fprintf(outfile, "%c", 24);
   fprintf(outfile, "%c", 0);

   for (n = 7; n <= 12; n++)
   {
      fprintf(outfile, "%c", headers[n] & 0x000000FF);
      fprintf(outfile, "%c", (headers[n] & 0x0000FF00) >> 8);
      fprintf(outfile, "%c", (headers[n] & 0x00FF0000) >> 16);
      fprintf(outfile, "%c", (headers[n] & (unsigned int) 0xFF000000) >> 24);
   }

   //
   // Headers done, now write the data...
   //

   for (y = bitmapstruct->height - 1; y >= 0; y--)     // BMP image format is written from bottom to top...
   {
      for (x = 0; x <= bitmapstruct->width - 1; x++)
      {
         // Also, it's written in (b,g,r) format...
   
         fprintf(outfile, "%c", getb(bitmapstruct, x, y));
         fprintf(outfile, "%c", getg(bitmapstruct, x, y));
         fprintf(outfile, "%c", getr(bitmapstruct, x, y));
      }
      if (extrabytes)      // See above - BMP lines must be of lengths divisible by 4.
      {
         for (n = 1; n <= extrabytes; n++)
         {
            fprintf(outfile, "%c", 0);
         }
      }
   }

   fclose(outfile);
   return;
}

////////////////////////////////////////////////////////////////////////////////////

// Some functions here that can be used as the bailout functions (see top)...

int bailoutlinear(int i, float multiplier)
{
   return i * multiplier;
}


int bailoutsqrt(int i, float multiplier)
{
   if (i < 0) return 0;
   return sqrt(i) * multiplier;
}


int bailoutlog(int i, float multiplier)
{
   if (i < 0) return 0;
   return log(i) * multiplier;
}


int bailoutsquared(int i, float multiplier)
{
   if (i < 0) return 0;
   return i * i * multiplier;
}

////////////////////////////////////////////////////////////////////////////////////

void init_imagemap(struct imagemap * bitmapstruct, int width, int height)
{
   bitmapstruct->width = width;
   bitmapstruct->height = height;
   bitmapstruct->pixels = (unsigned char *) malloc(bitmapstruct->width * bitmapstruct->height * 3 * sizeof(char));
   if (bitmapstruct->pixels == NULL)
   {
      printf("   Memory allocation failed. Quitting.\n");
      exit(1);
   }
   return;
}

int getr(struct imagemap * bitmapstruct, int x, int y)
{
   if (x < 0 || x >= bitmapstruct->width || y < 0 || y >= bitmapstruct->height)
   {
      return 0;
   }
   return (int) bitmapstruct->pixels[(x * 3) + (y * bitmapstruct->width * 3)];
}

int getg(struct imagemap * bitmapstruct, int x, int y)
{
   if (x < 0 || x >= bitmapstruct->width || y < 0 || y >= bitmapstruct->height)
   {
      return 0;
   }
   return (int) bitmapstruct->pixels[(x * 3) + 1 + (y * bitmapstruct->width * 3)];
}

int getb(struct imagemap * bitmapstruct, int x, int y)
{
   if (x < 0 || x >= bitmapstruct->width || y < 0 || y >= bitmapstruct->height)
   {
      return 0;
   }
   return (int) bitmapstruct->pixels[(x * 3) + 2 + (y * bitmapstruct->width * 3)];
}