File:Mandelbrot Islands of Consciousness.png

/* 
   c program: console, 1-file
   samm.c
   
   algorithms:
   - escape time
   - DEM/M = distance estimation, Milnor version 
   
   
   "we have an estimate dn for the distance from z0 to M, 
   where n is the first index for which |zn|>ER. 
   ...
   we compared this estimate to a fixed size s and set the pixel 
   * if dn<s  set gray = dn/s
   * white otherwise" A Cheritat
   
   
  
   
   --------------------------------
   1. draws Mandelbrot set for complex quadratic polynomial 
   Fc(z)=z*z +c
   using samm = Stripe Average Method/Coloring  
   -------------------------------         
   2. technique of creating ppm file is  based on the code of Claudio Rocchini
   http://en.wikipedia.org/wiki/Image:Color_complex_plot.jpg
   create 24 bit color graphic file ,  portable pixmap file = PPM 
   see http://en.wikipedia.org/wiki/Portable_pixmap
   to see the file use external application ( graphic viewer)
   -----
   it is example  for : 
   https://www.math.univ-toulouse.fr/~cheritat/wiki-draw/index.php/Mandelbrot_set
   -------------
   http://www.fractalforums.com/programming/dense-image/
   https://mathoverflow.net/questions/187842/is-there-an-almost-dense-set-of-quadratic-polynomials-which-is-not-in-the-inte
   -------------
   compile : 

 
 
   gcc d.c -lm -Wall 
 
 
   ./a.out
   
   convert samm.ppm -resize 800x800 samm.png

   -------- git --------
   
   
   cd existing_folder
   git init
   git remote add origin git@gitlab.com:adammajewski/mandelbrot_wiki_ACh.git
   git add samm.c
   git commit -m "samm + linear interpolation "
   git push -u origin master


   https://www.mathplanet.com/education/geometry/points,-lines,-planes-and-angles/finding-distances-and-midpoints
   middle point between 2 points: = (z1 + z2) / 2
   distance between 2 points = |z2-z1| = abs(z2-z1)
 
 
*/
#include <stdio.h>
#include <math.h>
#include <complex.h> // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c

 
#define LENGTH 5


int example = 1;

//  A viewport is a polygon viewing region in computer graphics


/* screen ( integer) coordinate ,  remember about aspect ratio !!!*/
int iX,iY;
const int iXmax = 8000; 
const int iYmax = 8000;

double width2; // = iWidth/2.0
double height2; // = iHeight/2.0


/* world ( double) coordinate = parameter plane, remember about aspect ratio !!! */
double Cx,Cy;

// c = -0.355298979690605  +1.178016112830515 i    period = 0
// c = 0.213487557597331  +0.604701020309395 i    period = 0
double CxMin;
double CxMax;
double CyMin;
double CyMax;
// 
complex double Center;
double Magnification;
double Radius;

/* */
double PixelWidth; //=(CxMax-CxMin)/iXmax;
double PixelHeight; // =(CyMax-CyMin)/iYmax;

/* color component ( R or G or B) is coded from 0 to 255 */
/* it is 24 bit color RGB file */
const int MaxColorComponentValue=255; 
FILE * fp;
char *filename="g8000_01.ppm"; // 
char *comment="# ";/* comment should start with # */
        
static unsigned char color[3]; // 24-bit rgb color


const int IterationMax=320000; //  N in wiki 


/* bail-out value for the bailout test for exaping points
   radius of circle centered ad the origin, exterior of such circle is a target set  */
const double EscapeRadius=2000; // = ER big !!!!
const double eps=0.01;

// DEM parameters
double thickness = 0.1;// multiplier 
double MinBoundaryWidth; //       
        
// number of pixels        
int iBoundaryPixels=0;
int iInteriorPixels=0;
int iExteriorPixels=0;
int iUnknownPixels=0;
int iAllPixels=0;
double density=0.0; // = iBoundary/all


// ------------------- examples -------------------
// all tables shoud have the same examples


/*
  {xmin,xmax, ymin,ymax}
  e[0 ] = center (-1.747370489, 0.004733564), magnification 4.42685e+06, 
  https://commons.wikimedia.org/wiki/File:Mandelbrot_Islands_of_Consciousness.jpg
  CxMin =  -1.747370789;    
  CxMax =  -1.747370189  ;
  CyMin = 0.004733264;
  CyMax = 0.004733864; 

  it depends of image size in pixels

  usually 640/480 = 4/3 =1,3333333 =
  the aspect ratio, or more precisely the Display Aspect Ratio (DAR) – the aspect ratio of the image as displayed; for TV, 
  DAR was traditionally 4:3 (a.k.a. Fullscreen), 
  with 16:9 (a.k.a. Widescreen) now the standard for HDTV

*/
double corners[LENGTH][4]={
  {-2.5,1.5,-1.5,1.5}, // standard viewport
  {-1.7473706,  -1.7473700,   0.0047332,   0.0047338}, //
  {-0.118608341441624798, 
   -0.118608339455450601, 
   -0.649678595428078265, 
   -0.649678593938447624} //https://commons.wikimedia.org/wiki/File:Fractal_08.jpg
};

/* -1.747370489, 0.004733564), magnification 4.42685e+06 
   https://dinkydauset.deviantart.com/art/Density-near-the-cardoid-662181713


*/
double CenterMag[LENGTH][3]={
  {-0.5, 0.0, 0.666666667},
  {-1.747370489, 0.004733564,4.42685e+06}
};

/* -1.747370489, 0.004733564), radsius 4.42685e+06 */
double CenterRadius[LENGTH][3]={
  {-0.5,0.0, 1.5},
  {-1.747370489, 0.004733564,0.0000002259}
};


       
        
double complex C_fromCorners(int iX, int iY){
  double Cx,Cy;
  
  Cy=CyMax - iY*PixelHeight; // reverse Y axis
  //if (fabs(Cy)<PixelHeight) Cy = 0.0; 
  Cx=CxMin + iX*PixelWidth;
   
  return Cx+Cy*I;
}


double complex C_fromCenterMag(int iX, int iY){
  double x,y;
  
  y = (iY - height2) / height2;// reverse Y axis
  //if (fabs(Cy)<PixelHeight) Cy = 0.0; 
  x= (iX - width2) / height2; ;
   
  return Center + (x - I * y)/Magnification;
}

 



int SetPlaneFromCorners(int i){
  if ((i >= LENGTH) || (i<0)) return 1; // error : bad input
  //
  CxMin = corners[i][0];
  CxMax = corners[i][1];
  CyMin = corners[i][2];
  CyMax = corners[i][3];
  //
  Radius = fabs(CyMax-CyMin)/2.0;
  Magnification = 1.0/Radius;
  
  // http://www.purplemath.com/modules/midpoint.htm
  Center = (CxMax+CxMin)/2.0 + I*(CyMax+CyMin)/2.0;
  
  
  //
  
  
  return 0; 
 


}

int SetPlaneFromCenterMag(int i){
  
  if ((i >= LENGTH) || (i<0)) return 1; // error : bad input
  
  Center = CenterMag[i][0] + CenterMag[i][1]*I;
  Magnification = CenterMag[i][2];
  //
  Radius = 1.0/Magnification;
  
  //
  width2  = iXmax/2.0;
  height2 = iYmax/2.0;
  
  complex double C1 = C_fromCenterMag(0,0);
  complex double C2 = C_fromCenterMag(iXmax,iYmax);
  CxMin = creal(C1);
  CxMax = creal(C2);
  CyMin = cimag(C2); // !!!!
  CyMax = cimag(C1);
  
  
  return 0; 
 


}


int SetPlaneFromCenterRadius(int i){
  
  if ((i >= LENGTH) || (i<0)) return 1; // error : bad input
  
  Center = CenterRadius[i][0] + CenterRadius[i][1]*I;
  Radius = CenterRadius[i][2];
  //
  Magnification = 1.0/Radius;
  
  //
  width2  = iXmax/2.0;
  height2 = iYmax/2.0;
  
  complex double C1 = C_fromCenterMag(0,0);
  complex double C2 = C_fromCenterMag(iXmax,iYmax);
  CxMin = creal(C1);
  CxMax = creal(C2);
  CyMin = cimag(C2); // !!!!
  CyMax = cimag(C1);
  
  
  return 0; 
}

// -------------------- end examples ------------------------




/*
  input :
  - complex number
  - intege
  output =  estimator dn
   
*/
double Give_D(double complex C , int iMax)
{
  int i=0; // iteration 
   
   
  double complex Z= C; // initial value for iteration Z0 is not 0 because of dZ  
  double R; // =radius = cabs(Z)
  
  complex double dZ = 1.0; // derivative with respect to z 
  double complex dC = 0; // derivative with respect to c 
  // dn in A Cheritat notation
  double de; // = 2 * z * log(cabs(z)) / dc =  estimated distance to the Mandelbrot set
   
    
  // iteration = computing the orbit
  for(i=0;i<iMax;i++)
    { 
      dZ = 2.0*dZ*Z;
      dC = 2.0*dC*Z + 1.0; 
      Z=Z*Z+C; // https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/qpolynomials
      
      //      
      R = cabs(Z);
      if(R > EscapeRadius)  // 
	{ de = 2.0 * R * log(R) / cabs(dC) ; // 
	  if (de < MinBoundaryWidth) { iBoundaryPixels +=1; return de/MinBoundaryWidth; }//  boundary = exterior near boundary
	  else   { iExteriorPixels +=1; return FP_ZERO;} //  exterior far from boundary
              
	} // 
      //    
      if(cabs(dZ) < eps) {iInteriorPixels +=1; return -1.0; } 
      
    } // for(i=0
   
   
  // if (i == iMax)  
  iUnknownPixels +=1; 
  return FP_INFINITE;
    
 
}
 
 
 
int SetColor(unsigned char color[3], unsigned char r, unsigned char g, unsigned char b) {

  //  gradient
  color[0]= r;  /* Red*/
  color[1]= g;  /* Green */ 
  color[2]= b;  /* Blue */
 
  return 0;


}

/* 
 
   input = complex number
   output 
   - color array as Formal parameters
   - int = return code
*/
int compute_color(complex double c, unsigned char color[3]){
 
  double d ;
  unsigned char b;
     
  // compute 
  d = Give_D( c, IterationMax);
  
  if (d >0)     { b =  (unsigned char )(255.0*d) % 255; SetColor(color,   b,   b,   b); return 0;} // boundary = gray gadient   
  
  if (d == FP_ZERO)          {SetColor(color, 255, 255, 255); return 0;} // exterior = white
  if (d < 0.0)          {SetColor(color, 255, 178,   0); return 0;} // interior = yelllow 
  if (d == FP_INFINITE) {SetColor(color, 255,   0,   0); return 0;} // unknown  = red
                      
 
    
 
 
 
  return 0;
}
 
 
 
void setup(int example ){
  printf(" render = compute and write image data bytes to the file \n");
  printf("example %d \n", example);
  
  SetPlaneFromCorners(example);
  //SetPlaneFromCenterMag(example);
  
  //
  PixelWidth=(CxMax-CxMin)/iXmax;
  PixelHeight=(CyMax-CyMin)/iYmax;
  
  
  
  // DEM
  MinBoundaryWidth = thickness*PixelWidth;
  
        
         
  /*create new file,give it a name and open it in binary mode  */
  fp= fopen(filename,"wb"); /* b -  binary mode */
  /*write ASCII header to the file*/
  fprintf(fp,"P6\n %s\n %d\n %d\n %d\n",comment,iXmax,iYmax,MaxColorComponentValue);
 
}
 



void info(){

  double distortion;
  // widt/height
  double PixelsAspectRatio = (double)iXmax/iYmax;  // https://en.wikipedia.org/wiki/Aspect_ratio_(image) 
  
 
  printf(" CxMin = %.16f  Re(Center) = %.16f  CxMax = %.16f \n", CxMin, creal(Center), CxMax);
  printf(" CyMin = %.16f  Im(Center) = %.16f  CyMax = %.16f \n", CyMin, cimag(Center), CyMax);   
 
  printf("Magnification = %.16f = %e \n", Magnification, Magnification);
  printf("Radius = %.16f = %e \n", Radius, Radius );
  double WorldAspectRatio = (CxMax-CxMin)/(CyMax-CyMin);
  printf("PixelsAspectRatio = %.16f \n", PixelsAspectRatio );
  printf("WorldAspectRatio = %.16f \n", WorldAspectRatio );
  distortion = PixelsAspectRatio - WorldAspectRatio;
  printf("distortion = %.16f ( it should be zero !)\n", distortion );
  printf("\n");
  printf("bailut value = Escape Radius = %.0f \n", EscapeRadius);
  printf("IterationMax = %d \n", IterationMax);
 
  iAllPixels = iXmax*iYmax;
  if (iAllPixels -(iExteriorPixels+iBoundaryPixels+iInteriorPixels+iUnknownPixels)==0) {
    printf("iAllPixels = %d \n", iAllPixels);
    printf("iBoundaryPixels = %d pixels = %f fraction of all pixels \n", iBoundaryPixels, (double)iBoundaryPixels/iAllPixels); 
    printf("iInteriorPixels = %d pixels = %f fraction of all pixels \n", iInteriorPixels, (double)iInteriorPixels/iAllPixels);   
    printf("iExteriorPixels = %d pixels = %f fraction of all pixels \n", iExteriorPixels, (double)iExteriorPixels/iAllPixels);
    printf("iUnknownPixels = %d pixels = %f fraction of all pixels \n", iUnknownPixels, (double)iUnknownPixels/iAllPixels);
    printf("density = %f = (iBoundaryPixels + iInteriorPixels)/iAllPixels \n", (double)(iBoundaryPixels+iInteriorPixels)/iAllPixels);    
  }
  else printf("error in counting pixels");
  
  
  

  // file  
  printf("file %s saved. \n", filename);

}
 



 
void close(){
 
  fclose(fp);
  info(); 

}
 
 
 
 
 
// ************************************* main ************************* 
int main()
{
        
  complex double c;
        
        
 
  setup(example);      
        
        
  
 
  for(iY=0;iY<iYmax;iY++)
    for(iX=0;iX<iXmax;iX++)
      { // compute pixel coordinate        
	//c = C_fromCenterMag(iX, iY);  
	c = C_fromCorners(iX,iY);
	/* compute  pixel color (24 bit = 3 bytes) */
	compute_color(c,color);         
	/*write color to the file*/
	fwrite(color,1,3,fp);
      }
        
  
  
  close();
  
         
  return 0;
}