File:From decomposition to checkerboard.gif

From Wikimedia Commons, the free media repository
Jump to navigation Jump to search
From_decomposition_to_checkerboard.gif ‎(600 × 600 pixels, file size: 1.84 MB, MIME type: image/gif, looped, 26 frames, 13 s)

Summary[edit]

Description
English: From decomposition to checkerboard. It is insipired by the program Mandel by Wolf Jung. See algorithm 9 : zeros of qn(c)[1]
Date
Source own work with help of Wolf Jung
Author Adam majewski

Licensing[edit]

I, the copyright holder of this work, hereby publish it under the following license:
w:en:Creative Commons

attribution share alike

This file is licensed under the Creative Commons Attribution-Share Alike 4.0 International license.
You are free:
  • to share – to copy, distribute and transmit the work
  • to remix – to adapt the work
Under the following conditions:
  • attribution – You must attribute the work in the manner specified by the author or licensor (but not in any way that suggests that they endorse you or your use of the work).
  • share alike – If you alter, transform, or build upon this work, you may distribute the resulting work only under the same or similar license to this one.

See also[edit]

  • Cauliflower = parabolic Julia set

  • binary escape time showing filled -in Julia set

  • Julia set

  • Decomposition

  • Binary decomposition of whole dynamical plane with circle Julia set c = 0

  • Binary decomposition along internal ray 0

  • Modified binary decomposition of whole dynamical plane with circle Julia set. Image with src code

  • Another modification of binary decomposition

  • Modified decomposition of dynamical plane with dendrite Julia set. Image with src code

  • Modified decomposition of dynamical plane with basilica Julia set. Image with src code

Summary[edit]

This animated gif file shows parabolic Julia set called Cauliflower.

Gif consist of 26 frames numbered from 0 to 25 with the number of iterations ( n).

making animated gif[edit]

  • first find exterior ( escaping set) and filled-in julia set ( non-escaping or bounded set) using simple escape time ( see ComputeFatouComponents function in c src code). This image is saved in data array. It is common for all frames and this image is not saved, but one can do it by uncomment SaveArray2PGMFile( data,999); line
  • then make all pgm images ( see ComputeColorZero funtion in C src code and zero array)
  • convert pgm images to frames of animated gif using bash src code

Frames and plane[edit]

Each frame of animated gif shows rectangle part of the dynamic plane : 

/* world ( double) coordinate = dynamic plane */
static   const double ZxMin=-1.5;
static  const double ZxMax=1.5;
static  const double ZyMin=-1.5;
static  const double ZyMax=1.5;

for different number of iterations ( n) which is displayed in left upper corner

Arrays[edit]

There are 3 arrays :

  • data for main components of dynamic plane ( escaping and non-escaping set). See escape time algorithm and ComputeColorZero funcion
  • edge for boundaries of components : Julia set, but also for boundaries between other components. See Sobel filter and ComputeBoundariesIn function
  • zero for components of Zero algorithm. See ComputeColorZero function
// memmory 1D array 
unsigned char *data;
unsigned char *edge;
unsigned char *zero;

Each array can be :

  • saved as a pgm image ( see SaveArray2PGMFile function)
  • copied to other array ( see CopyBoundariesTo function )
  • used to save information about colour of pixels ( shades of gray = char = integer from 0 to 255 )

In each array there 3 types of coordinate are :

  • i = index of 1D array
  • ix and iy = indices of virtual 2D array ( image ). These are integer = pixel = screeen coordinate
  • zx and zy = world or double coordinate ( of dynamic plane )

There are also functions for conversion :

  • i = Give_i(ix,iy)
  • Zx = GiveZx(ix) and Zy = GiveZy(iy)

Colouring algorithms[edit]

There are 3 colouring algorithms:

Program flow :

  • ComputeZerosOfQnc function calls PlotPointZero function for each element of the array zero
  • function PlotPointZero calls ComputeColorZero function and computes colour of each pixel
  • Zero algorithm is implemented in the ComputeColorZero function

ComputeColorZero function :

  • takes as the input :
    • integer coordinate of the point ( ix, iy)
    • maximal number of iterations (n=iMax)
  • converts integer to double coordinate of pixel ( from screen to world, here dynamic z-plane ) : z = Zx+Zy*I
  • makes (iMax) of iterations iterations : zn = fn(z0)
  • checks zn and returns colour ( char number iColor )
// check position near fixed point after n iterations 
  if ( Zx>0 && Zy>0) iColor = 150; //re(z_n) > 0 and im(z_n) > 0 (first quadrant)
  if ( Zx<0 && Zy>0) iColor = 170; //re(z_n) < 0 and im(z_n) > 0 (second)
  if ( Zx<0 && Zy<0) iColor = 190; //re(z_n) < 0 and im(z_n) < 0 (third)
  if ( Zx>0 && Zy<0) iColor = 200; //re(z_n) > 0 and im(z_n) < 0 (fourth).

Note that in Zero algorithm :

the number of iterations[edit]

The maximal number of iterations is displayed in left upper corner of each frame.

In the source code it is denoted by n variable : 

 int n;

Note that for number of iterations ( n) :

  • n<10 both exterior ( escaping set) and interior ( non-escaping = bounded set) is coloured with the same algorithm ( see ComputeColorZero funtion)
  • n>=10 only interior is coloured with ComputeColorZero funcion. Exterior is coloured with simple solid color ( iColorOfExterior )
  • number of details is proportional to n. See what happens near the cusps of Julia set, especially when N>=10 and all exterior has the same colour
  • for a big n only 2 colours stay inside !!

Look at the code : 

iColor = ComputeColorZero(ix, iy, n);
  if (data[i]==iColorOfInterior) A[i] =  iColor+20; // interior 
     else {if (n<10) A[i]= iColor; else A[i]= iColorOfExterior;} // exterior , only for low n

C source code[edit]

The src code was formatted with Emacs 

/*

  Adam Majewski
  adammaj1 aaattt o2 dot pl  // o like oxygen not 0 like zero 
  fraktal.republika.pl

 https://gitlab.com/adammajewski/zero_algorithm_cauliflower_in_c

  c console progam 

  gcc c.c -lm -Wall -march=native 
  time ./a.out

  gcc c.c -lm -Wall -march=native -fopenmp

  time ./a.out

  From program Mandel by Wolf Jung http://www.mndynamics.com/indexp.html
  On parameter plane for complex quadratic map algorithm " 9 is showing the location of centers / p.p. The period is set with q. (Use 9 also to check the displayed period, which might be inaccurate.)"

*/

#include <stdio.h>
#include <stdlib.h> // malloc
#include <string.h> // strcat
#include <math.h> // M_PI; needs -lm also 
#include <complex.h>
#include <omp.h>

/* --------------------------------- global variables and consts ------------------------------------------------------------ */

#define iPeriodChild 0 // Period of secondary component joined by root point with the parent component 
int iPeriodParent = 1; // main cardioid of Mandelbrot set

// virtual 2D array and integer ( screen) coordinate
// Indexes of array starts from 0 not 1 
//unsigned int ix, iy; // var
static unsigned int ixMin = 0; // Indexes of array starts from 0 not 1
static unsigned int ixMax ; //
static unsigned int iWidth ; // horizontal dimension of array

static unsigned int iyMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iyMax ; //

static unsigned int iHeight = 5000; //  
// The size of array has to be a positive constant integer 
static unsigned int iSize ; // = iWidth*iHeight; 

// memmory 1D array 

unsigned char *data;
unsigned char *edge;
unsigned char *zero;

// unsigned int i; // var = index of 1D array
//static unsigned int iMin = 0; // Indexes of array starts from 0 not 1
static unsigned int iMax ; // = i2Dsize-1  = 
// The size of array has to be a positive constant integer 
// unsigned int i1Dsize ; // = i2Dsize  = (iMax -iMin + 1) =  ;  1D array with the same size as 2D array

/* world ( double) coordinate = dynamic plane */
static   const double ZxMin=-1.5;
static  const double ZxMax=1.5;
static  const double ZyMin=-1.5;
static  const double ZyMax=1.5;
static  double PixelWidth; // =(ZxMax-ZxMin)/ixMax;
static  double PixelHeight; // =(ZyMax-ZyMin)/iyMax;
static  double ratio ;
 

// complex numbers of parametr plane 
double Cx; // c =Cx +Cy * i
double Cy;
double complex c; // parameter of function fc(z)=z^2 + c

static unsigned long int iterMax  = 10000; //iHeight*100;

static double ER = 2.0; // Escape Radius for bailout test 
static double ER2;

/* colors = shades of gray from 0 to 255 */
// 8 bit color = int number from 0 to 255

// Arrays are 0-indexed, so the first array element is at index = 0, and the highest is =(size_of_array – 1) 
//unsigned char iColorInterior[2][iPeriodChild]={{255,231}, {123,99}}; /* shades of gray used in image */
unsigned char iColors[4]={255,231,180, 160}; //
static unsigned char iColorOfExterior = 245;
static unsigned char iColorOfInterior = 200;
unsigned char iColorOfUnknown = 133;
//int iNumberOfUknknown = 0;

//static double TwoPi=2*M_PI;

/* ------------------------------------------ functions -------------------------------------------------------------*/

//------------------complex numbers -----------------------------------------------------

// from screen to world coordinate ; linear mapping
// uses global cons
double GiveZx(unsigned int ix)
{ return (ZxMin + ix*PixelWidth );}

// uses globaal cons
double GiveZy(unsigned int iy)
{ return (ZyMax - iy*PixelHeight);} // reverse y axis

/* -----------  array functions = drawing -------------- */

/* gives position of 2D point (ix,iy) in 1D array  ; uses also global variable iWidth */
unsigned int Give_i(unsigned int ix, unsigned int iy)
{ return ix + iy*iWidth; }

// plots raster point (ix,iy) 
int iDrawPoint(unsigned char A[], unsigned int ix, unsigned int iy, unsigned char iColor)
{ 

  /* i =  Give_i(ix,iy) compute index of 1D array from indices of 2D array */
  A[Give_i(ix,iy)] = iColor;

  return 0;
}

// draws point to memmory array data
// uses complex type so #include <complex.h> and -lm 
int dDrawPoint(unsigned char A[], complex double point,unsigned char iColor )
{

  unsigned int ix, iy; // screen coordinate = indices of virtual 2D array
  //unsigned int i; // index of 1D array
  
  ix = (creal(point)- ZxMin)/PixelWidth; 
  iy = (ZyMax - cimag(point))/PixelHeight; // inverse Y axis 
  iDrawPoint(A, ix, iy, iColor);
  return 0;
}

//;;;;;;;;;;;;;;;;;;;;;;  setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

int setup(int ParentPeriod, int ChildPeriod)
{

  
  
  
  
  
  printf("setup\n");

  

  
  Cx=0.25;
  Cy=0.0;
  c=Cx+Cy*I;
  
  

  /* 2D array ranges */
  if (!(iHeight % 2)) iHeight+=1; // it sholud be even number (variable % 2) or (variable & 1)
  iWidth = iHeight;
  iSize = iWidth*iHeight; // size = number of points in array 
  // iy
  iyMax = iHeight - 1 ; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].
  //ix
  
  ixMax = iWidth - 1;

  /* 1D array ranges */
  // i1Dsize = i2Dsize; // 1D array with the same size as 2D array
  iMax = iSize-1; // Indexes of array starts from 0 not 1 so the highest elements of an array is = array_name[size-1].

  /* Pixel sizes */
  PixelWidth = (ZxMax-ZxMin)/ixMax; //  ixMax = (iWidth-1)  step between pixels in world coordinate 
  PixelHeight = (ZyMax-ZyMin)/iyMax;
  ratio = ((ZxMax-ZxMin)/(ZyMax-ZyMin))/((float)iWidth/(float)iHeight); // it should be 1.000 ...
  
  

  // for numerical optimisation in iteration
  ER2 = ER * ER;
  
  
  /* create dynamic 1D arrays for colors ( shades of gray ) */
  data = malloc( iSize * sizeof(unsigned char) );
  edge = malloc( iSize * sizeof(unsigned char) );
  zero = malloc( iSize * sizeof(unsigned char) );

  if (edge== NULL || data == NULL || zero==NULL)
    {
      fprintf(stderr," Could not allocate memory");
      getchar(); 
      return 1;
    }

  
  

 
  

  
   
  
  printf(" end of setup \n");
  
  return 0;

} // ;;;;;;;;;;;;;;;;;;;;;;;;; end of the setup ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

unsigned char ComputeColor(unsigned int ix, unsigned int iy, int IterationMax)
{ 
  // check behavour of z under fc(z)=z^2+c
  // using 1 target set:
  // 1. exterior or circle (center at origin and radius ER ) 
  // as a target set containing infinity = for escaping points ( bailout test)
  // for points of exterior of julia set
  

  double Zx2, Zy2;
  int i=0;  // number of the iteration = fc(z)
  
 
  double Zx, Zy;
  
 
 
  
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  

  // if not inside target set around attractor ( alfa fixed point )
  while (1 )
    { // then iterate 
      
      Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
      // bailout test 
      if (Zx2 + Zy2 > ER2) return iColorOfExterior; // if escaping stop iteration
       
      // if not escaping or not attracting then iterate = check behaviour
      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx; 
      //
      i+=1;
      if (i > IterationMax) break; 
      
     
       
                  
	  
    }
    
      
      
      
    
   
  return   iColorOfInterior;   //
}

// plots raster point (ix,iy) 
int PlotPoint(unsigned char A[] , unsigned int ix, unsigned int iy, int IterationMax)
{
  unsigned i; /* index of 1D array */
  unsigned char iColor;
  

  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  iColor = ComputeColor(ix, iy, IterationMax);
  A[i] = iColor;

  return 0;
}

// fill array 
// uses global var :  ...
// scanning complex plane 
int ComputeFatouComponents(unsigned char A[], int IterationMax )
{
  unsigned int ix, iy; // pixel coordinate 

  //printf("compute image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(ixMax , iyMax, IterationMax)
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d \r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPoint(A, ix, iy, IterationMax ) ; //  
    } 
   
  return 0;
}

int ComputeBoundariesIn(unsigned char A[]) // compute in A, but save to edge
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
  /* sobel filter */
  unsigned char G, Gh, Gv; 
  // boundaries are in edge array ( global var )
 
  // clear all pixels
  for(i=1;i<iSize;++i) edge[i]=iColorOfExterior; 
 
 
  printf(" find boundaries in A array using  Sobel filter\n");   
  // #pragma omp parallel for schedule(dynamic) private(i,iY,iX,Gv,Gh,G) shared(iyMax,ixMax, ER2)
  for(iY=1;iY<iyMax-1;++iY){ 
    for(iX=1;iX<ixMax-1;++iX){ 
      Gv= A[Give_i(iX-1,iY+1)] + 2*A[Give_i(iX,iY+1)] + A[Give_i(iX-1,iY+1)] - A[Give_i(iX-1,iY-1)] - 2*A[Give_i(iX-1,iY)] - A[Give_i(iX+1,iY-1)];
      Gh= A[Give_i(iX+1,iY+1)] + 2*A[Give_i(iX+1,iY)] + A[Give_i(iX-1,iY-1)] - A[Give_i(iX+1,iY-1)] - 2*A[Give_i(iX-1,iY)] - A[Give_i(iX-1,iY-1)];
      G = sqrt(Gh*Gh + Gv*Gv);
      i= Give_i(iX,iY); /* compute index of 1D array from indices of 2D array */
      if (G==0) {edge[i]=255;} /* background */
      else {edge[i]=0;}  /* boundary */
    }
  }
 
   
 
  return 0;
}

int CopyBoundariesTo(unsigned char A[]) // copy boundaries from edge to A
{
 
  unsigned int iX,iY; /* indices of 2D virtual array (image) = integer coordinate */
  unsigned int i; /* index of 1D array  */
 
 
  printf("copy boundaries from edge array to data array \n");
  for(iY=1;iY<iyMax-1;++iY)
    for(iX=1;iX<ixMax-1;++iX)
      {i= Give_i(iX,iY); if (edge[i]==0) A[i]=0;}
 
 
 
  return 0;
}

// Check Orientation of image : mark first quadrant 
// it should be in the upper right position
// uses global var :  ...
int CheckOrientation(unsigned char A[] )
{
  unsigned int ix, iy; // pixel coordinate 
  double Zx, Zy; //  Z= Zx+ZY*i;
  unsigned i; /* index of 1D array */
  for(iy=iyMin;iy<=iyMax;++iy) 
    {
      Zy = GiveZy(iy);
      for(ix=ixMin;ix<=ixMax;++ix) 
	{

	  // from screen to world coordinate 
	  Zx = GiveZx(ix);
	  i = Give_i(ix, iy); /* compute index of 1D array from indices of 2D array */
	  if (Zx>0 && Zy>0) A[i]=255-A[i];   // check the orientation of Z-plane by marking first quadrant */

	}
    }
   
  return 0;
}

unsigned char ComputeColorZero(unsigned int ix, unsigned int iy, int iMax)
{ 
  

  double Zx2, Zy2;
  int i=0;  // number of the iteration = fc(z)
  
  unsigned char iColor;
  double Zx, Zy;
  
 
 
  
  
  // from screen to world coordinate 
  Zx = GiveZx(ix);
  Zy = GiveZy(iy);
  

  
  while (i <iMax)   
    { // then iterate 
      
      Zx2 = Zx*Zx; 
      Zy2 = Zy*Zy;
       
             
      
      // new z : Z(n+1) = Zn * Zn  + C
      Zy = 2*Zx*Zy + Cy; 
      Zx = Zx2 - Zy2 + Cx; 
      //
      i+=1;
    }

  // check position near fixed point after n iterations 
  if ( Zx>0 && Zy>0) iColor = 150; //re(z_n) > 0 and im(z_n) > 0 (first quadrant)
  if ( Zx<0 && Zy>0) iColor = 170; //re(z_n) < 0 and im(z_n) > 0 (second)
  if ( Zx<0 && Zy<0) iColor = 190; //re(z_n) < 0 and im(z_n) < 0 (third)
  if ( Zx>0 && Zy<0) iColor = 200; //re(z_n) > 0 and im(z_n) < 0 (fourth).
  //
  return iColor;  
}

// plots raster point (ix,iy) 
int PlotPointZero(unsigned char A[] , unsigned int ix, unsigned int iy, int n)
{
  unsigned i; /* index of 1D array */
  unsigned char iColor;
  

  i = Give_i(ix,iy); /* compute index of 1D array from indices of 2D array */
  iColor = ComputeColorZero(ix, iy, n);
  if (data[i]==iColorOfInterior) A[i] =  iColor+20; // interior 
  else {if (n<10) A[i]= iColor; else A[i]= iColorOfExterior;} // exterior , only for low n 

  return 0;
}

// fill array 
// uses global var :  ...
// scanning complex plane 
int ComputeZerosOfQnc(unsigned char A[], int n)
{
  unsigned int ix, iy; // pixel coordinate 

  //printf("compute image \n");
  // for all pixels of image 
#pragma omp parallel for schedule(dynamic) private(ix,iy) shared(ixMax , iyMax, n)
  for(iy = iyMin; iy<=iyMax; ++iy) 
    { printf(" %d z %d \r", iy, iyMax); //info 
      for(ix= ixMin; ix<=ixMax; ++ix) PlotPointZero(A, ix, iy, n ) ; //  
    } 
   
  return 0;
}
 

// save "A" array to pgm file 
int SaveArray2PGMFile( unsigned char A[], double k)
{
  
  FILE * fp;
  const unsigned int MaxColorComponentValue=255; /* color component is coded from 0 to 255 ;  it is 8 bit color file */
  char name [30]; /* name of file */
  sprintf(name,"%.0f", k); /*  */
  char *filename =strcat(name,".pgm");
  char *comment="# ";/* comment should start with # */

  /* save image to the pgm file  */      
  fp= fopen(filename,"wb"); /*create new file,give it a name and open it in binary mode  */
  fprintf(fp,"P5\n %s\n %u %u\n %u\n",comment,iWidth,iHeight,MaxColorComponentValue);  /*write header to the file*/
  fwrite(A,iSize,1,fp);  /*write A array to the file in one step */
  printf("File %s saved. \n", filename);
  fclose(fp);

  return 0;
}

int info()
{
  // diplay info messages
  printf("Numerical approximation of parabolic Julia set for fc(z)= z^2 + c \n");
  printf("iPeriodParent = %d \n", iPeriodParent);
  printf("iPeriodOfChild  = %d \n", iPeriodChild);
  printf("parameter c = ( %f ; %f ) \n", Cx, Cy); 
  
  printf("Image Width = %f \n", ZxMax-ZxMin);
  printf("PixelWidth = %f \n", PixelWidth);
  printf("Maximal number of iterations = iterMax = %ld \n", iterMax);
  printf("ratio of image  = %f ; it should be 1.000 ...\n", ratio);
  return 0;
}

/* -----------------------------------------  main   -------------------------------------------------------------*/
int main()
{
  setup(iPeriodParent, iPeriodChild);
  int n; 

 
   
    
 
  printf("components of Fatou set  : \n");
  ComputeFatouComponents(data, iterMax);
  //SaveArray2PGMFile( data, ); 
  printf("done \n");

  
  
  for(n = 0; n<=25; ++n)
    { 
      ComputeZerosOfQnc(zero, n);
      ComputeBoundariesIn(zero); // from zero to edge
      CopyBoundariesTo(zero); // copy from edge to zero
      SaveArray2PGMFile( zero, n); 
    }
  
 
 

  printf(" allways free memory  to avoid buffer overflow \n");
 
  free(data);
  free(edge);
  free(zero);

  
  info();

  return 0;
}

Bash ang Image Magic src code[edit]

 #!/bin/bash
 
# script file for BASH 
# which bash
# save this file as g.sh
# chmod +x g.sh
# ./g.sh
 

 
# for all ppm files in this directory
for file in *.pgm ; do
  # b is name of file without extension
  b=$(basename $file .pgm)
  # convert from pgm to gif and add text ( level ) using ImageMagic
  convert $file -resize 600x600 -pointsize 80 -annotate +500+100 $b ${b}.gif
  echo $file
done
 
# convert gif files to animated gif
convert  -delay 50  -loop 0 %d.gif[0-25] a.gif
 
echo OK
  1. ↑ Mandel: software for real and complex dynamics by Wolf Jung

File history

Click on a date/time to view the file as it appeared at that time.

Date/TimeThumbnailDimensionsUserComment
current10:47, 28 November 2015Thumbnail for version as of 10:47, 28 November 2015600 × 600 (1.84 MB)Adam majewski (talk | contribs)User created page with UploadWizard
  • You cannot overwrite this file.

File usage on other wikis

Metadata

This file contains additional information such as Exif metadata which may have been added by the digital camera, scanner, or software program used to create or digitize it. If the file has been modified from its original state, some details such as the timestamp may not fully reflect those of the original file. The timestamp is only as accurate as the clock in the camera, and it may be completely wrong.
Retrieved from "https://commons.wikimedia.org/w/index.php?title=File:From_decomposition_to_checkerboard.gif&oldid=295940048"