Fractals/Iterations in the complex plane/MandelbrotSetExterior
Colouring of exterior of Mandelbrot set can be :
- non-smooth :
- Boolean/binary Escape Time Method ( bETM/M )
- discrete = Level Set Method = LSM/M = iETM/M
- Smooth :
- radial measures
- Real Escape Time Method( rETM/M )
- Distance Estimation Method( DEM/M )
- radius of complex potential = Complex Potential Method ( CPM/M )
- radial measures
- angular measures
One can also draw curves :
- external rays
- equipotential lines ( closed curves - quasi circles)
Similar projects:
- Mandelbrot Notebook by Claude Heiland-Allen
- different drawing techniques and algorithms by Arnaud Cheritat
- Linas Vepstas' Art Gallery:
Escape time or dwell[edit | edit source]
Here for given point c on parameter plane one checks how critical point behaves on dynamical plane under forward iteration. If you change initial point you will get different result [5]
To draw given plane one needs to check/scan (all) its points. See here for more details ( optimisation) Read definitions first.
How to find the number of iterations required to escape the mandelbrot set ?[edit | edit source]
- math.stackexchange question: is-there-an-equation-for-the-number-of-iterations-required-to-escape-the-mandelb
- math.stackexchange question: a-way-to-determine-the-ideal-number-of-maximum-iterations-for-an-arbitrary-zoom?
Boolean escape time[edit | edit source]
This algorithm answers the question: “For which values of c will the Julia fractal, J(c), be line-like and for which dust-like?”[6]
Here complex plane consists of 2 sets : Mandelbrot set and its complement :
ASCI graphic ( on screen)[edit | edit source]
// http://mrl.nyu.edu/~perlin/
main(k){float i,j,r,x,y=-16;while(puts(""),y++<15)for(x
=0;x++<84;putchar(" .:-;!/>)|&IH%*#"[k&15]))for(i=k=r=0;
j=r*r-i*i-2+x/25,i=2*r*i+y/10,j*j+i*i<11&&k++<111;r=j);}
-- Haskell code by Ochronus
-- http://blog.mostof.it/mandelbrot-set-in-ruby-and-haskell/
import Data.Complex
mandelbrot a = iterate (\z -> z^2 + a) a !! 500
main = mapM_ putStrLn [[if magnitude (mandelbrot (x :+ y)) < 2 then '*' else ' '
| x <- [-2, -1.9685 .. 0.5]]
| y <- [1, 0.95 .. -1]]
; common lisp
(loop for y from -1.5 to 1.5 by 0.05 do
(loop for x from -2.5 to 0.5 by 0.025 do
(let* ((c (complex x y)) ; parameter
(z (complex 0 0))
(iMax 20) ; maximal number of iterations
(i 0)) ; iteration number
(loop while (< i iMax ) do
(setq z (+ (* z z) c)) ; iteration
(incf i)
(when (> (abs z) 2) (return i)))
; color of pixel
(if (= i iMax) (princ (code-char 42)) ; inside M
(princ (code-char 32))))) ; outside M
(format t "~%")) ; new line
Comparison programs in various languages [7][8]
Graphic file ( PPM )[edit | edit source]
Here are various programs for creating pbm file [9]
C[edit | edit source]
This is complete code of C one file program.
- It makes a ppm file which consists an image. To see the file (image) use external application ( graphic viewer).
- Program consists of 3 loops:
- iY and iX, which are used to scan rectangle area of parameter plane
- iterations.
For each point of screen (iX,iY) it's complex value is computed c=cx+cy*i.
For each point c is computed iterations of critical point
It uses some speed_improvement. Instead of checking :
sqrt(Zx2+Zy2)<ER
it checks :
(Zx2+Zy2)<ER2 // ER2 = ER*ER
It gives the same result but is faster.
/*
c program:
--------------------------------
1. draws Mandelbrot set for Fc(z)=z*z +c
using Mandelbrot algorithm ( boolean escape time )
-------------------------------
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)
*/
#include <stdio.h>
#include <math.h>
int main()
{
/* screen ( integer) coordinate */
int iX,iY;
const int iXmax = 800;
const int iYmax = 800;
/* world ( double) coordinate = parameter plane*/
double Cx,Cy;
const double CxMin=-2.5;
const double CxMax=1.5;
const double CyMin=-2.0;
const double CyMax=2.0;
/* */
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="new1.ppm";
char *comment="# ";/* comment should start with # */
static unsigned char color[3];
/* Z=Zx+Zy*i ; Z0 = 0 */
double Zx, Zy;
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
/* */
int Iteration;
const int IterationMax=200;
/* bail-out value , radius of circle ; */
const double EscapeRadius=2;
double ER2=EscapeRadius*EscapeRadius;
/*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);
/* compute and write image data bytes to the file*/
for(iY=0;iY<iYmax;iY++)
{
Cy=CyMin + iY*PixelHeight;
if (fabs(Cy)< PixelHeight/2) Cy=0.0; /* Main antenna */
for(iX=0;iX<iXmax;iX++)
{
Cx=CxMin + iX*PixelWidth;
/* initial value of orbit = critical point Z= 0 */
Zx=0.0;
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
/* */
for (Iteration=0;Iteration<IterationMax && ((Zx2+Zy2)<ER2);Iteration++)
{
Zy=2*Zx*Zy + Cy;
Zx=Zx2-Zy2 +Cx;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
/* compute pixel color (24 bit = 3 bytes) */
if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
else
{ /* exterior of Mandelbrot set = white */
color[0]=255; /* Red*/
color[1]=255; /* Green */
color[2]=255;/* Blue */
};
/*write color to the file*/
fwrite(color,1,3,fp);
}
}
fclose(fp);
return 0;
}
Integer escape time = LSM/M = dwell bands[edit | edit source]
Here color is proportional to last iteration ( of final_n, final iteration).[11]
This is also called Level Set Method ( LSM )
C[edit | edit source]
Difference between Mandelbrot algorithm and LSM/M is in only in part instruction, which computes pixel color of exterior of Mandelbrot set. In LSM/M is :
if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set = LSM */
else if ((Iteration%2)==0)
{ /* even number = black */
color[0]=0; /* Red */
color[1]=0; /* Green */
color[2]=0; /* Blue */
}
else
{/* odd number = white */
color[0]=255; /* Red */
color[1]=255; /* Green */
color[2]=255; /* Blue */
};
Here is C function whithout explicit complex numbers, only doubles:
int GiveEscapeTime(double C_x, double C_y, int iMax, double _ER2)
{
int i;
double Zx, Zy;
double Zx2, Zy2; /* Zx2=Zx*Zx; Zy2=Zy*Zy */
Zx=0.0; /* initial value of orbit = critical point Z= 0 */
Zy=0.0;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
for (i=0;i<iMax && ((Zx2+Zy2)<_ER2);i++)
{
Zy=2*Zx*Zy + C_y;
Zx=Zx2-Zy2 +C_x;
Zx2=Zx*Zx;
Zy2=Zy*Zy;
};
return i;
}
here a short code with complex numbers:
// https://gitlab.com/adammajewski/mandelbrot_wiki_ACh/blob/master/betm.c
int iterate(double complex C , int iMax)
{
int i;
double complex Z= 0.0; // initial value for iteration Z0
for(i=0;i<iMax;i++)
{
Z=Z*Z+C; // https://stackoverflow.com/questions/6418807/how-to-work-with-complex-numbers-in-c
if(cabs(Z)>EscapeRadius) break;
}
return i;
}
C++[edit | edit source]
Here is C++ function which can be used to draw LSM/M :
int iterate_mandel(complex C , int imax, int bailout)
{
int i;
std::complex Z(0,0); // initial value for iteration Z0
for(i=0;i<=imax-1;i++)
{
Z=Z*Z+C; // overloading of operators
if(abs(Z)>bailout)break;
}
return i;
}
I think that it can't be coded simpler (it looks better than pseudocode), but it can be coded in other way which can be executed faster .
Here is faster code :
// based on cpp code by Geek3
inline int fractal(double cx, double cy, int max_iters)
// gives last iteration
{
double zx = 0, zy = 0;
if (zx * zx + zy * zy > 4) return(0); // it=0
for (int it = 1; it < max_iters; it++)
{ double zx_old = zx;
zx = zx * zx - zy * zy;
zy = 2 * zx_old * zy;
zx += cx;
zy += cy;
if (zx * zx + zy * zy > 4.0) return(it);
}
return(max_iters);
}
A touch more optimised :
// optimised from cpp code by Geek3
inline int fractal(double cReal, double cImg, int max_iters)
// gives last iteration
{
double zReal = 0, zImg = 0, zReal2 = 0, zImg2 = 0;
//iteration zero is always 0^2+0^2, it will never escape
for (int it = 1; it < max_iters; it++)
{ //because we have zReal^2 and zImg^2 pre-calculated
//we can calculate zImg first
//then we don't need to calculate/store the "old" zReal
zImg = (2 * zReal * zImg ) + cImg;
zReal = zReal2 - zImg2 + cReal;
// calculate next iteration: zReal^2 and zImg^2
// they are used twice so calculate them once
zReal2 = zReal * zReal;
zImg2 = zImg * zImg;
if (zReal2 + zImg2 > 4.0) return(it);
}
return(max_iters);
}
See also :
GLSL[edit | edit source]
Java[edit | edit source]
//Java code by Josef Jelinek
// http://java.rubikscube.info/
int mandel(double cx, double cy) {
double zx = 0.0, zy = 0.0;
double zx2 = 0.0, zy2 = 0.0;
int iter = 0;
while (iter < iterMax && zx2 + zy2 < 4.0) {
zy = 2.0 * zx * zy + cy;
zx = zx2 - zy2 + cx;
zx2 = zx * zx;
zy2 = zy * zy;
iter++;
}
return iter;
}
Java Script[edit | edit source]
Here is JavaScript function which does not give last iteration but LastIteration modulo maxCol. It makes colour cycling ( if maxCol < maxIt ).
function iterate(Cr,Ci) {
// JavaScript function by Evgeny Demidov
// http://www.ibiblio.org/e-notes/html5/fractals/mandelbrot.htm
var I=0, R=0, I2=0, R2=0, n=0;
if (R2+I2 > max) return 0;
do { I=(R+R)*I+Ci; R=R2-I2+Cr; R2=R*R; I2=I*I; n++;
} while ((R2+I2 < max) && (n < maxIt) );
if (n == maxIt) return maxCol; else return n % maxCol;
}
Above functions do not use explicit definition of complex number.
Khan Academy[edit | edit source]
Lisp program[edit | edit source]
Whole Lisp program making ASCII graphic based on code by Frank Buss [12] [13]
; common lisp
(loop for y from -1.5 to 1.5 by 0.1 do
(loop for x from -2.5 to 0.5 by 0.04 do
(let* ((i 0)
(z (complex x y))
(c z))
(loop while (< (abs
(setq z (+ (* z z) c)))
2)
while (< (incf i) 32))
(princ (code-char (+ i 32))))) ; ASCII chars <= 32 contains non-printing characters
(format t "~%"))
MathMap plugin for Gimp[edit | edit source]
filter mandelbrot (gradient coloration) c=ri:(xy/xy:[X,X]*1.5-xy:[0.5,0]); z=ri:[0,0]; # initial value z0 = 0 # iteration of z iter=0; while abs(z)<2 && iter<31 do z=z*z+c; # z(n+1) = fc(zn) iter=iter+1 end; coloration(iter/32) # color of pixel end
Pov-Ray[edit | edit source]
Pov-Ray has a built-in function mandel[14]
Wolfram Mathematica[edit | edit source]
Here is code by Paul Nylander
Level Curves of escape time Method = LCM/M[edit | edit source]
Lemniscates are boundaries of Level Sets of escape time ( LSM/M ). They can be drawn using :
- edge detection of Level sets.
- Algorithm described in paper by M. Romera et al.[15] This method is fast and allows looking for high iterations.
- boundary trace[16]
- drawing curves , see explanation and source code. This method is very complicated for iterations > 5.
Decomposition of target set for Mandelbrot set drawing[edit | edit source]
Decomposition is modification of escape time algorithm.
The target set is divided into parts (2 or more). Very large escape radius is used, for example ER = 12.
Binary decomposition of LSM/M[edit | edit source]
Here target set on dynamic plane is divided into 2 parts (binary decomposition = 2-decomposition ):
- upper half ( white)
- lower half (black)
Division of target set induces decomposition of level sets into parts ( cells, subsets):
- which is colored white,
- which is colored black.
External rays of angles (measured in turns):
can be seen as borders of subsets.
Difference between binary decomposition algorithm and Mandel or LSM/M is in only in part of instruction , which computes pixel color of exterior of Mandelbrot set. In binary decomposition is :
if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set = LSM */
else if (Zy>0)
{
color[0]=0; /* Red */
color[1]=0; /* Green */
color[2]=0; /* Blue */
}
else
{
color[0]=255; /* Red */
color[1]=255; /* Green */
color[2]=255; /* Blue */
};
also GLSL code from Fragmentarium :
#include "2D.frag"
#group Simple Mandelbrot
// maximal number of iterations
uniform int iMax; slider[1,100,1000] // increase iMax
// er2= er^2 wher er= escape radius = bailout
uniform float er2; slider[4.0,1000,10000] // increase er2
// compute color of pixel
vec3 color(vec2 c) {
vec2 z = vec2(0.0); // initial value
// iteration
for (int i = 0; i < iMax; i++) {
z = vec2(z.x*z.x-z.y*z.y,2*z.x*z.y) + c; // z= z^2+c
if (dot(z,z)> er2) // escape test
// exterior
if (z.x>0){ return vec3( 1.0);} // upper part of the target set
else return vec3(0.0); //lower part of the target set
}
return vec3(0.0); //interior
}
// zoomasm -- zoom video assembler
// (c) 2019,2020,2021,2022 Claude Heiland-Allen
// SPDX-License-Identifier: AGPL-3.0-only
// recommended KF bailout settings: linear smoothing, custom radius 25
vec3 colour(void)
{
if (getInterior())
{
return vec3(1.0, 0.0, 0.0);
}
bool decomp = getT() < 0.5;
return vec3(decomp ? 0.0 : 1.0);
}
Point c is plotting white or black if imaginary value of last iteration ( Zy) is positive or negative.[17]
nth-decomposition[edit | edit source]
This method is extension of binary decomposition.
The target set T = { z : |zn| > R } with a very large escape radius ( for example R = 12 ) is divided into more then 2 parts ( for example 8).[18]
Real Escape Time[edit | edit source]
Other names of this method/algorithm are :
- the fully-renormalized fractional iteration count ( by Linas Vepstas in 1997)[19]
- smooth iteration count for generalized Mandelbrot sets ( by Inigo Quilez in 2016)[20]
- continuous iteration count for the Mandelbrot set
- Normalized Iteration Count Algorithm
- Continuous coloring
- smooth colour gradient
- fractional iterations
- fractional escape time
Here color of exterior of Mandelbrot set is proportional not to Last Iteration ( which is integer number) but to real number :
Other methods and speedups
Colouring formula in Ultrafractal :[21]
smooth iter = iter + 1 + ( log(log(bailout)-log(log(cabs(z))) )/log(2)
where :
- log(log(bailout) can be precalculated
C[edit | edit source]
To use log2 function add :
#include <math.h>
at the beginning of program.
if (Iteration==IterationMax)
{ /* interior of Mandelbrot set = black */
color[0]=0;
color[1]=0;
color[2]=0;
}
/* exterior of Mandelbrot set */
else GiveRainbowColor((double)(Iteration- log2(log2(sqrt(Zx2+Zy2))))/IterationMax,color);
where :
- Zx2 = Zx*Zx
- Zy2 = Zy*Zy
Here is another version by Tony Finch[22]
while (n++ < max &&
x2+y2 < inf) {
y = 2*x*y + b;
x = x2-y2 + a;
y2 = y*y;
x2 = x*x;
}
nu = n - log(log(x2+y2)/2)/ log(2);
based on equation [23]
C++[edit | edit source]
// based on cpp code by Geek3 from http://en.wikibooks.org/wiki/File:Mandelbrot_set_rainbow_colors.png
sqrxy = x * x + y * y;
double m = LastIteration + 1.5 - log2(log2(sqrxy));
java[edit | edit source]
/**
Smooth coloring algorithm
https://gitlab.com/shreyas.siravara/mandelbrot-with-smooth-coloring/blob/master/Mandelbrot.java
Mandelbrot with Smooth Coloring by Shreyas Siravara
*/
double nsmooth = (iterations + 1 - Math.log(Math.log(Zn.getMagnitude())) / Math.log(ESCAPE_RADIUS));
double smoothcolor = nsmooth / MAX_ITERATIONS;
if (iterations < MAX_ITERATIONS) {
int rgb = Color.HSBtoRGB((float) (0.99f + 1.9 * smoothcolor), 0.9f, 0.9f);
g2d.setColor(new Color(rgb));
} else {
g2d.setColor(Color.black.darker());
}
Matemathica[edit | edit source]
Here is code by Paul Nylander. It uses different formula :
Python[edit | edit source]
Python code using mpmath library[24]
def mandelbrot(z):
c = z
for i in xrange(ITERATIONS):
zprev = z
z = z*z + c
if abs(z) > ESCAPE_RADIUS:
return ctx.exp(1j*(i + 1 - ctx.log(ctx.log(abs(z)))/ctx.log(2)))
return 0
Distance estimation DEM/M[edit | edit source]
Variants :
- exterior DEM/M
- interior DEM/M
Complex potential[edit | edit source]
See also[edit | edit source]
- https://web.archive.org/web/20071008112609/http://rgba.scenesp.org/iq/trastero/fieldlines/
- http://fraktal.republika.pl/mset_bottcher.html
References[edit | edit source]
- ↑ Mathematics of Divergent Fractals by
- ↑ jussi harkonen : coloring-techniques
- ↑ wikipedia : Orbit trap
- ↑ Mandelbrot Orbit Trap Rendering! Programming How-To Video by DKM101
- ↑ Java program by Dieter Röß showing result of changing initial point of Mandelbrot iterations
- ↑ Julia fractals in PostScript by Kees van der Laan, EUROTEX 2012 & 6CM PROCEEDINGS 47
- ↑ Fractal Benchmark by Erik Wrenholt
- ↑ 12-minute Mandelbrot: fractals on a 50 year old IBM 1401 mainframe
- ↑ Computer Language Benchmarks Game
- ↑ example-code-from-presentation-ways-of-seeing-julia-sets by ed Burke
- ↑ Computing the Mandelbrot set by Andrew Williams
- ↑ LIsp Program by Frank Buss
- ↑ Mandelbrot Set ASCII art at Bill Clementson's blog
- ↑ mandel function from 2.5.11.14 Fractal Patterns at Pov-Ray docs
- ↑ Drawing the Mandelbrot set by the method of escape lines. M. Romera et al.
- ↑ http://www.metabit.org/~rfigura/figura-fractal/math.html boundary trace by Robert Figura
- ↑ http://web.archive.org/20010415125044/www.geocities.com/CapeCanaveral/Launchpad/5113/fr27.htm%7C An open letter to Dr. Meech from Joyce Haslam in FRACTAL REPORT 27
- ↑ mandelbrot set n-th-decomposition
- ↑ linas.org : Renormalizing the Mandelbrot Escape
- ↑ I Quilez : mset_smooth
- ↑ fractalforums : What range/precision for fractional escape counts for Mandelbrot/Julia sets?
- ↑ Making Mandelbrot Set Movies by Tony Finch
- ↑ Linas Vepstas. Renormalizing the mandelbrot escape.
- ↑ mpmath Python library