Boolean Escape time

Algorithm: for every point z of dynamical plane (z-plane) compute iteration number ( last iteration) for which magnitude of z is greater than escape radius. If last_iteration=max_iteration then point is in filled-in Julia set, else it is in its complement (attractive basin of infinity ). Here one has 2 options, so it is named boolean algorithm.

if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else color=WHITE;  /* unbounded orbits = exterior of Filled-in Julia set  */

In theory this method is for drawing Filled-in Julia set and its complement ( exterior), but when c is Misiurewicz point ( Filled-in Julia set has no interior) this method draws nothing. For example for c=i . It means that it is good for drawing interior of Filled-in Julia set.

; common lisp
(loop for y from -2 to 2 by 0.05 do
      (loop for x from -2 to 2 by 0.025 do
		(let* ((z (complex x y))
                   	(c (complex -1 0))
                   	(iMax 20)
			(i 0))

		(loop  	while (< i iMax ) do 
			(setq z (+ (* z z) c))
			(incf i)
			(when (> (abs z) 2) (return i)))
			

           (if (= i iMax) (princ (code-char 42)) (princ (code-char 32)))))
      (format t "~%"))

PPM file with raster graphic

Filled-in Julia set for c= =-1+0.1*i. Image and C source code

Integer escape time = Level Sets of the Basin of Attraction of Infinity = Level Sets Method= LSM/J

• 
  level sets of escape time for C=0 and 0.9<Z<1.5
• 

Escape time measures time of escaping to infinity ( infinity is superattracting point for polynomials). Time is measured in steps ( iterations = i) needed to escape from circle of given radius ( ER= Escape Radius).

One can see few things:

• this is discontinuous function
• i is iMax for z in Filled-in Julia set
• i=0 for x0>ER
• this is nonlinear function

Level sets here are sets of points with the same escape time. Here is algorithm of choosing color in black & white version.

  if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else   /* unbounded orbits = exterior of Filled-in Julia set  */
        if ((LastIteration%2)==0) /* odd number */
           then color=BLACK; 
           else color=WHITE;

Here is the c function which:

• uses complex double type numbers
• computes 8 bit color ( shades of gray)
• checks both escape and attraction test

unsigned char ComputeColorOfLSM(complex double z){

 int nMax = 255;
  double cabsz;
  unsigned char iColor;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	cabsz = cabs(z);
    	if (cabsz > ER) break; // esacping
    	if (cabsz< PixelWidth) break; // fails into finite attractor = interior
  			
   
     	z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  
  iColor = 255 - 255.0 * ((double) n)/20; // nMax or lower walues in denominator
  
  
  return iColor;
}

Normalized iteration count (real escape time or fractional iteration or Smooth Iteration Count Algorithm (SICA))

• 
  integer escape time = LSM
• 
  real escape time
• 
  Escape time for C=0 and 0.5<Z<2.5

Math formula :

${\displaystyle \nu (z)=\lim \limits _{i\to \infty }(i-\log _{2}\log _{2}|z_{i}|)\ .}$

Maxima version :

GiveNormalizedIteration(z,c,E_R,i_Max):=
/* */
block(
 [i:0,r],
 while abs(z)<E_R and i<i_Max
   do (z:z*z + c,i:i+1),
 r:i-log2(log2(cabs(z))),
 return(float(r))
)$

In Maxima log(x) is a natural (base e) logarithm of x. To compute log2 use :

log2(x) := log(x) / log(2);

description:

• FF : julia-smooth-colouring-how-to-do
• stefan bion : fraktal-generator/colormapping/
• FF smooth-histogram-rendering/
• FF creating-a-good-palette-using-bezier-interpolation/

Level Curves of escape time Method = eLCM/J

• 
  edge detection Image and C code
• 
  preimages of circle

These curves are boundaries of Level Sets of escape time ( eLSM/J ). They can be drawn using these methods:

• edge detection of Level Curves ( =boundaries of Level sets).
  • Algorithm based on paper by M. Romera et al.^[7]
  • Sobel filter
• drawing lemniscates = curves ${\displaystyle L_{n}=\{z:abs(z_{n})=ER\}\,}$, see explanation and source code
• drawing circle ${\displaystyle L_{0}=\{z:abs(z)=ER\}\,}$ and its preimages. See this image, explanation and source code
• method described by Harold V. McIntosh^[8]

/* Maxima code : draws lemniscates of Julia set */
c: 1*%i;
ER:2;
z:x+y*%i;
f[n](z) := if n=0 then z else (f[n-1](z)^2 + c);
load(implicit_plot); /* package by Andrej Vodopivec */ 
ip_grid:[100,100];
ip_grid_in:[15,15];
implicit_plot(makelist(abs(ev(f[n](z)))=ER,n,1,4),[x,-2.5,2.5],[y,-2.5,2.5]);

Components of Interior of Filled Julia set ( Fatou set)

• use limited color ( palette = list of numbered colors)
• find period of attracting cycle
• find one point of attracting cycle
• compute number of iteration after when point reaches the attractor
• color of component=iteration % period^[9]
• use edge detection for drawing Julia set

• 
  Components of Fatou set for period 4
• 
  Components of Fatou set for period 3
• 
  In case of Siegel disc critical orbit is a boundary Siegel disc compponent. All other componnats are preimages of this component

Internal Level Sets

See :

• algorithm 0 of program Mandel by Wolf Jung

• 
  c = 0.5*i. Image and source code
• 
  c = -0.584895. Image and source code

Binary decomposition

Here color of pixel ( exterior of Julia set) is proportional to sign of imaginary part of last iteration .

Main loop is the same as in escape time.

In other words target set is decompositioned in 2 parts ( binary decomposition) :

${\displaystyle T_{b}+=\{z:abs(z_{n})>ER~~{\mbox{and}}~~im(z_{n})>0\}\,}$

${\displaystyle T_{b}-=\{z:abs(z_{n})>ER~~{\mbox{and}}~~im(z_{n})<=0\}\,}$

Algorithm in pseudocode ( Im(Zn) = Zy ) :

if (LastIteration==IterationMax)
   then color=BLACK;   /* bounded orbits = Filled-in Julia set */
   else   /* unbounded orbits = exterior of Filled-in Julia set  */
        if (Zy>0) /* Zy=Im(Z) */
           then color=BLACK; 
           else color=WHITE;

Modified decomposition

Here exterior of Julia set is decompositioned into radial level sets.

It is because main loop is without bailout test and number of iterations ( iteration max) is constant.

It creates radial level sets. See also video by bryceguy72^[10] and by FreymanArt^[11]

  for (Iteration=0;Iteration<8;Iteration++)
 /* modified loop without checking of abs(zn) and low iteration max */
    {
    Zy=2*Zx*Zy + Cy;
    Zx=Zx2-Zy2 +Cx;
    Zx2=Zx*Zx;
    Zy2=Zy*Zy;
   };
  iTemp=((iYmax-iY-1)*iXmax+iX)*3;        
  /* --------------- compute  pixel color (24 bit = 3 bajts) */
 /* exterior of Filled-in Julia set  */
 /* binary decomposition  */
  if (Zy>0 ) 
   { 
    array[iTemp]=255; /* Red*/
    array[iTemp+1]=255;  /* Green */ 
    array[iTemp+2]=255;/* Blue */
  }
  if (Zy<0 )
   {
    array[iTemp]=0; /* Red*/
    array[iTemp+1]=0;  /* Green */ 
    array[iTemp+2]=0;/* Blue */    
  };

It is also related with automorphic function for the group of Mobius transformations ^[12]

How to use distance

One can use distance for colouring  :

• only Julia set ( boundary of filled Julia set)
• boundary and exterior of filled Julia set.

Here is first example :

 if (LastIteration==IterationMax) 
     then { /*  interior of Julia set, distance = 0 , color black */ }
     else /* exterior or boundary of Filled-in Julia set  */
          {  double distance=give_distance(Z0,C,IterationMax);
             if (distance<distanceMax)
                 then { /*  Julia set : color = white */ }
                 else  { /*  exterior of Julia set : color = black */}
          }

Here is second example ^[16]

 if (LastIteration==IterationMax or distance < distanceMax) then ... // interior by ETM/J and boundary by DEM/J
 else .... // exterior by real escape time

Zoom

DistanceMax is smaller than pixel size. During zooming pixel size is decreasing and DistanceMax should also be decreased to obtain good picture. It can be made by using formula :

${\displaystyle DistanceMax={\frac {PixelSize}{n}}\,}$

where ${\displaystyle n\leq 1\,}$

One can start with n=1 and increase n if picture is not good. Check also iMax !!

DistanceMax may also be proportional to zoom factor ${\displaystyle mag\,}$:^[17]

${\displaystyle DistanceMax={\sqrt {\frac {thick}{1000*mag}}}}$

where thick is image width ( in world units) and mag is a zoom factor.

Examples of code

For cpp example see mndlbrot::dist from mndlbrot.cpp in src code of program mandel ver 5.3 by Wolf Jung.

C function using complex type :

unsigned char ComputeColorOfDEMJ(complex double z){
// https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Julia_set#DEM.2FJ


  
  int nMax = iterMax;
  complex double dz = 1.0; //  is first derivative with respect to z.
  double distance;
  double cabsz;
	
  int n;

  for (n=0; n < nMax; n++){ //forward iteration
	cabsz = cabs(z);
    	if (cabsz > 1e60 || cabs(dz)> 1e60) break; // big values
    	if (cabsz< PixelWidth) return iColorOfInterior; // falls into finite attractor = interior
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  if (distance <distanceMax) {//printf(" distance = %f \n", distance); 
  		return iColorOfBoundary;}
  
  
  return iColorOfExterior;

 
}

C function using double type:

/*based on function  mndlbrot::dist  from  mndlbrot.cpp
 from program mandel by Wolf Jung (GNU GPL )
 http://www.mndynamics.com/indexp.html  */
double jdist(double Zx, double Zy, double Cx, double Cy ,  int iter_max)
 { 
 int i;
 double x = Zx, /* Z = x+y*i */
         y = Zy, 
         /* Zp = xp+yp*1 = 1  */
         xp = 1, 
         yp = 0, 
         /* temporary */
         nz,  
         nzp,
         /* a = abs(z) */
         a; 
 for (i = 1; i <= iter_max; i++)
  { /* first derivative   zp = 2*z*zp  = xp + yp*i; */
    nz = 2*(x*xp - y*yp) ; 
    yp = 2*(x*yp + y*xp); 
    xp = nz;
    /* z = z*z + c = x+y*i */
    nz = x*x - y*y + Cx; 
    y = 2*x*y + Cy; 
    x = nz; 
    /* */
    nz = x*x + y*y; 
    nzp = xp*xp + yp*yp;
    if (nzp > 1e60 || nz > 1e60) break;
  }
 a=sqrt(nz);
 /* distance = 2 * |Zn| * log|Zn| / |dZn| */
 return 2* a*log(a)/sqrt(nzp); 
 }

Delphi function :

function Give_eDistance(zx0,zy0,cx,cy,ER2:extended;iMax:integer):extended;

var zx,zy ,  // z=zx+zy*i
    dx,dy,  //d=dx+dy*i  derivative :  d(n+1)=  2 * zn * dn
    zx_temp,
    dx_temp,
    z2,  //
    d2, //
    a // abs(d2)
     :extended;
    i:integer;
begin
  //initial values
  // d0=1
  dx:=1;
  dy:=0;
  //
  zx:=zx0;
  zy:=zy0;
  // to remove warning : variables may be not initialised ?
  z2:=0;
  d2:=0;

  for i := 0 to iMax - 1 do
    begin
      // first derivative   d(n+1) = 2*zn*dn  = dx + dy*i;
      dx_temp := 2*(zx*dx - zy*dy) ;
      dy := 2*(zx*dy + zy*dx);
      dx := dx_temp;
      // z = z*z + c = zx+zy*i
      zx_temp := zx*zx - zy*zy + Cx;
      zy := 2*zx*zy + Cy;
      zx := zx_temp;
      //
      z2:=zx*zx+zy*zy;
      d2:=dx*dx+dy*dy;
      if ((z2>1e60) or (d2 > 1e60)) then  break;
      
    end; // for i
   if (d2 < 0.01) or (z2 < 0.1)  // when do not use escape time
    then  result := 10.0
    else
      begin
        a:=sqrt(z2);
        // distance = 2 * |Zn| * log|Zn| / |dZn|
        result := 2* a*log10(a)/sqrt(d2);
      end;

end;  //  function Give_eDistance

Matlab code by Jonas Lundgren^[18]

function D = jdist(x0,y0,c,iter,D)
%JDIST Estimate distances to Julia set by potential function
% by Jonas Lundgren http://www.mathworks.ch/matlabcentral/fileexchange/27749-julia-sets
% Code covered by the BSD License http://www.mathworks.ch/matlabcentral/fileexchange/view_license?file_info_id=27749

% Escape radius^2
R2 = 100^2;

% Parameters
N = numel(x0);
M = numel(y0);
cx = real(c);
cy = imag(c);
iter = round(1000*iter);

% Create waitbar
h = waitbar(0,'Please wait...','name','Julia Distance Estimation');
t1 = 1;

% Loop over pixels
for k = 1:N/2
    x0k = x0(k);
    for j = 1:M
        % Update distance?
        if D(j,k) == 0
            % Start values
            n = 0;
            x = x0k;
            y = y0(j);
            b2 = 1;                 % |dz0/dz0|^2
            a2 = x*x + y*y;         % |z0|^2
            % Iterate zn = zm^2 + c, m = n-1
            while n < iter && a2 <= R2
                n = n + 1;
                yn = 2*x*y + cy;
                x = x*x - y*y + cx;
                y = yn;
                b2 = 4*a2*b2;       % |dzn/dz0|^2
                a2 = x*x + y*y;     % |zn|^2
            end
            % Distance estimate
            if n < iter
                % log(|zn|)*|zn|/|dzn/dz0|
                D(j,k) = 0.5*log(a2)*sqrt(a2/b2);
            end
        end
    end
    % Lap time
    t = toc;
    % Update waitbar
    if t >= t1
        str = sprintf('%0.0f%% done in %0.0f sec',200*k/N,t);
        waitbar(2*k/N,h,str)
        t1 = t1 + 1;
    end
end

% Close waitbar
close(h)

Maxima function :

 GiveExtDistance(z0,c,e_r,i_max):=
 /* needs z in exterior of Kc */
 block(
 [z:z0,
 dz:1,
 cabsz:cabs(z),
 cabsdz:1, /* overflow limit */
 i:0],
 while  cabsdz < 10000000 and  i<i_max
  do 
   (
    dz:2*z*dz,
    z:z*z + c,
    cabsdz:cabs(dz),
    i:i+1
   ),
  cabsz:cabs(z), 
  return(2*cabsz*log(cabsz)/cabsdz)
 )$

Exterior of circle

This is typical target set. It is exterior of circle with center at origin ${\displaystyle z=0\,}$ and radius =ER :

${\displaystyle T_{ER}=\{z:abs(z)>ER\}\,}$

Radius is named escape radius ( ER ) or bailout value. Radius should be greater than 2.

Exterior of square

Here target set is exterior of square of side length ${\displaystyle s\,}$ centered at origin

${\displaystyle T_{s}=\{z:abs(re(z))>s~~{\mbox{or}}~~abs(im(z))>s\}\,}$