Fractals/Computer graphic techniques/2D/optimisation

Symmetry

Mirror symmetry

Parameter plane and the Mandelbrot set is symmetric with respect to the x-axis in the plane :^[1]

colour(x,y) = colour(x,-y)

Here is a C code that shows how to divide symmetric image (with axis in the middle of its height) into 3 parts :

axis of symmetry
2 symmetric parts : above and below axis

// fill array using mirror symmetry of image 
// uses global var :  ...
int FillArray(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
// draw axis of symmetry
iy = iyAxisOfSymmetry; 
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
 
// draw symmetric parts : above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ix, iyAxisOfSymmetry - iyAbove , Color ); 
   }   
 return 0;
}

See also Pascal code ( Delphi) for general case ^[2]

Rotational symmetry

Dynamical plane and the Julia set have rotational symmetry. It can be used to speed up computing :

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry; 
 
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}

↑Jump back a section

Complex numbers

Computations without explicit definition of complex numbers or without complex type are usually faster.

↑Jump back a section

Parallel computing

In Mandelbrot and Julia sets each point/pixel is calculated independently, so it can be easly divided into an smaller tasks and carried out simultaneously.^[3]

Types of parallel computing : Wikipedia description^[4]

Multi-core (computing)
- Intel MIC^[5]
- MPI ^[6]
- OpenMP^[7]^[8]^[9]
General-purpose computing on graphics processing units (GPGPU):
- OpenCl^[10]
- CUDA^[11]
- WebGl ^[12]^[13]

Vectorisation

C

gcc vector extension ^[14]^[15]^[16]^[17]
Intel® Advanced Vector Extensions^[18]^[19]

R

This image is made with use of vectorised code. It's creation time is very short.

Vectorisation ^[20]^[21]

Here is R code of above image with comments: ^[22]

## based on the R code by Jarek Tuszynski
## http://commons.wikimedia.org/wiki/File:Mandelbrot_Creation_Animation_%28800x600%29.gif

iMax=20; # number of frames and iterations

## world coordinate ( float )
cxmin = -2.2;
cxmax = 1.0;
cymax = 1.2;
cymin = -1.2;

## screen coordinate ( integer)

dx=800; dy=600 ;           # define grid size

## sequences of  values = rows and columns
cxseq = seq(cxmin, cxmax, length.out=dx);
cyseq = seq(cymin, cymax, length.out=dy);
## sequences of  values = rows * columns
cxrseq = rep(cxseq, each = dy);
cyrseq = rep(cyseq,  dx);
C = complex( real=cxrseq, imag = cyrseq); # complex vector
C = matrix(C,dy,dx)       # convert from vector to matrix
# Now C is a matrix of c points coordinates (cx,cy)

# allocate memory for all the frames
F = array(0, dim = c(dy,dx,iMax)) # dy*dx*iMax array filled with zeros

Z = 0                     # initialize Z to zero
for (i in 1:iMax)         # perform iMax iterations
{
  Z = Z^2+C               # iteration; uses n point to compute n+1 point
  
  F[,,i]  = exp(-abs(Z))              # an omitted index is used to represent an entire column or row
     # store magnitude of the complex number, scale it using an exponential function to emphasize fine detail, 
} 


# color and save F array to the file 
library(caTools)          # load library with write.gif function
# mapped to color palette jetColors . Dark red is a very low number, dark blue is a very high number
jetColors = colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan", "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
write.gif(F, "Mandelbrot.gif", col=jetColors, delay=100, transparent=0)

Multicore

Image and src code using OpenMP and symmetry

"I always create as many worker threads as I have cores. More than that and your system spends too much time task switching" - Duncan C^[23]

Using OpenMP :

// fill array using symmetry of image 
// uses global var :  ...
int FillArraySymmetric(unsigned char data[] )
{
 
 unsigned char Color; // gray from 0 to 255 
 
printf("axis of symmetry \n"); 
iy = iyAxisOfSymmetry; 
#pragma omp parallel for schedule(dynamic) private(ix,Color) shared(ixMin,ixMax, iyAxisOfSymmetry)
for(ix=ixMin;ix<=ixMax;++ix) PlotPoint(ix, iy, GiveColor(ix, iy)); 
 
 
/*
The use of ‘shared(variable, variable2) specifies that these variables should be shared among all the threads.
The use of ‘private(variable, variable2)’ specifies that these variables should have a seperate instance in each thread.
*/
 
#pragma omp parallel for schedule(dynamic) private(iyAbove,ix,iy,Color) shared(iyAboveMin, iyAboveMax,ixMin,ixMax, iyAxisOfSymmetry)
 
// above and below axis 
for(iyAbove = iyAboveMin; iyAbove<=iyAboveMax; ++iyAbove) 
  {printf(" %d from %d\n", iyAbove, iyAboveMax); //info 
  for(ix=ixMin; ix<=ixMax; ++ix) 
 
  { // above axis compute color and save it to the array
    iy = iyAxisOfSymmetry + iyAbove;
    Color = GiveColor(ix, iy);
    PlotPoint(ix, iy, Color ); 
    // below the axis only copy Color the same as above without computing it 
    PlotPoint(ixMax-ix, iyAxisOfSymmetry - iyAbove , Color ); 
   } 
}  
 return 0;
}