This algorithm has 2 versions:

Compare it with version for parameter plane and Mandelbrot set : DEM/M

External distance estimation

Distance Estimation Method for Julia set ( DEM/J ) estimates distance from point z ( in the exterior of Julia set ) to nearest point in Julia set.

c = 0.255
c= -0.75+0.11
c=-0.1+0.651
c=-0.74543+0.11301*i
c=- 0.181502832839439 -0.582826014844503*I

For distance estimate it has been proved that the computed value differs from the true distance at most by a factor of 4:

   Koebe Quarter Theorem. The image of an injective analytic function f : D → C from the unit disk D onto a subset of the complex plane contains the disk whose center is f(0) and whose radius is |f′(0)|/4.^[1]

Math formula :

$distance(z_{0},K_{c},n)=2*|z_{n}|*{\frac {\log |z_{n}|}{|z'_{n}|}}\,$

where :

$z'_{n}\,$ is first derivative with respect to z.

This derivative can be found by iteration starting with

z'_{0}=1\,

and then :

$z'_{n}=2*z_{n-1}*z'_{n-1}\,$

Pseudocode and the code

The Beauty of Fractals
The Science of Fractal Images, page 198,
Distance Estimator by Robert P. Munafo^[2]

Pseudocode by Claude Heiland-Allen^[3]

foreach pixel c
  while not escaped and iteration limit not reached
    dz := 2 * z * dz + 1
    z := z^2 + c
  de := 2 * |z| * log(|z|) / |dz|
  d := de / pixel_spacing
  plot_grey(tanh(d))

First derivative wrt z
degree	Function f(z)	derivative wrt z
2	$z^{2}+c$	$2*z$
3	$z^{3}+c$	$3*z^{2}$
4	$z^{4}+c$	$4*z^{3}$
d	$z^{d}+c$	$d*z^{d-1}$

// ********************************************************************************************************************
/* -----------------------------------------  basic function ( iteration)  -------------------------------------------------------------*/
// ********************************************************************************************************************

// update with f function 
const char *f_description = "Numerical approximation of Julia set for f(z)= z^3 + c "; // without /n !!!
/* ------------------------------------------ functions -------------------------------------------------------------*/

// complex function
// upadte f_description also
complex double f(const complex double z0) {

  double complex z = z0;
  z = z*z*z + c;
  return  z;
}
	
complex double derivative_wrt_z(const complex double z0) {

  double complex z = z0;
  z = 3.0*z*z ;
  return  z;
}



/* ************************** DEM/J*****************************************
 
it can be used for 
* whole image thru Compute8BitColor function
* only drawing boundary thru 

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


	int nMax = IterMax_DEM;
	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

		if (cabs2(z)> ER2 || cabs(dz)> 1e60) break; // big values
  			
		dz = derivative_wrt_z(z) * dz; 
		z  = f(z); /* forward iteration : complex quadratic polynomial */ 
	}
  
	if (n==nMax) return iColorOfInterior; // not escaping
	
	
	// escaping and boundary
	cabsz = cabs(z);
	distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
	double g = tanh(distance / PixelWidth);
  
	return 255*g; // iColorOfExterior;
}

How to use distance

tanh

double g = tanh(distance / PixelWidth);
return 255*g; // iColorOfExterior;

distance max

Distance max is the old ( deprecated) method.

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 ^[4]

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

Zoom

distance max

Distance max is the old ( deprecated) method.

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 :

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

where $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 $mag\,$ :^[5]

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

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

One can use also tanh which gives more precise look:

distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
 double g = tanh(distance /PixelWidth);
 return 255*g; // iColorOfExterior;

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_DEM;
  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

    if (cabs2(z)> ER2 || cabs(dz)> 1e60) break; // big values
    
  			
    dz = 2.0*z * dz; 
    z = z*z +c ; /* forward iteration : complex quadratic polynomial */ 
  }
  
  if (n==nMax) return iColorOfInterior; // not escaping
  // escaping and boundary
  cabsz = cabs(z);
  distance = 2.0 * cabsz* log(cabsz)/ cabs(dz);
  double g = tanh(distance /PixelWidth);
  return 255*g; // 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^[6]

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)
 )$

shadertoy

Julia - Distance 1 by iq. Analytical distance to a Julia set z^2 + c where vec2 c = vec2( -0.745, 0.186 )
Julia - Distance 2 by iq = SDF for the Julia set of f(z) = z^3+C where const vec2 kC = vec2(0.105,0.7905);
Julia - Distance 3 by iq. A generic Julia set renderer, using distance to the set (Douady_Hubbard potential). In this case I'm using a rational function of order 6: f(z) = (z-(1+i)/10)(z-i)(z-1)^4 / ((z+1)(z-(1+i)) + c
Inigo Quilez :: articles :: distance to fractals - 2004

 // Julia - Distance 2 by iq
        // compute Julia set
        const float threshold = 64.0;
        const vec2  kC = vec2(0.105,0.7905);
        const int   kNumIte = 200;

        float it = 0.0;
        float dz2 = 1.0;
        float m2 = 0.0;
        for( int i=0; i<kNumIte; i++ )
        {
            // df(z)/dz = 3*z^2
            dz2 *= 9.0*dot2(vec2(z.x*z.x-z.y*z.y,2.0*z.x*z.y));
            // f(z) = z^3 + c
            z = vec2( z.x*z.x*z.x - 3.0*z.x*z.y*z.y, 3.0*z.x*z.x*z.y - z.y*z.y*z.y ) + kC;
            // check divergence
            it++;
            m2 = dot2(z);
            if( m2>threshold ) break;
        }
        
        // distance
        float d = 0.5 * log(m2) * sqrt(m2/dz2);
        // interation count
        float h = it - log2(log2(dot(z,z))/(log2(threshold)))/log2(3.0); // https://iquilezles.org/articles/msetsmooth
        
        // coloring
        vec3 tmp = vec3(0.0);
        if( it<(float(kNumIte)-0.5) )
        {
            #if COLOR_SCHEMA==0
            tmp = 0.5 + 0.5*cos( 5.6 + sqrt(h)*0.5 + vec3(0.0,0.15,0.2));
            tmp *= smoothstep(0.0,0.0005,d);
            tmp *= 1.2/(0.3+tmp);
            tmp = pow(tmp,vec3(0.4,0.55,0.6));
            #else
            tmp = vec3(0.12,0.10,0.09);
            tmp *= smoothstep(0.005,0.020,d);
            float f = smoothstep(0.0005,0.0,d);
            tmp += 3.0*f*(0.5+0.5*cos(3.5 + sqrt(h)*0.4 + vec3(0.0,0.5,1.0)));
            tmp = clamp(tmp,0.0,1.0);
			#endif
        }
        
        col += vec4(tmp*w,w);
	#if AA>1
    }
    col /= col.w;
	#endif

    return col.xyz;

Internal distance estimation

Colouring the Julia set by Gert Buschmann

Julia set drawn from distance estimation

In order to get a nice picture, we must also colour the Julia set, since otherwise the Julia set is only visible through the colouring of the Fatou domains, and this colouring changes vigorously near the Julia set, giving a muddy look (it is possible to avoid this by choosing the colour scale and the density carefully). A point z belongs to the Julia set if the iteration does not stop, that is, if we have reached the chosen maximum number of iterations, M. But as the Julia set is infinitely thin, and as we only perform calculations for regularly situated points, in practice we cannot colour the Julia set in this way. But happily there exists a formula that (up to a constant factor) estimates the distance from the points z outside the Julia set to the Julia set. This formula is associated to a Fatou domain, and the number given by the formula is the more correct the closer we come to the Julia set, so that the deviation is without significance.

The distance function is the function $\delta (z)=\phi (z)/|\phi '(z)|$ (see the section Julia and Mandelbrot sets for non-complex functions), whose equipotential lines must lie approximately regularly. In the formula appears the derivative $z'_{k}$ of $z_{k}$ with respect to z. But as $z_{k}=f(f(...f(z)))$ (the k-fold composition), $z'_{k}$ is the product of the numbers $f'(z_{i})$ (i = 0, 1, ..., k-1), and this sequence can be calculated recursively by $z'_{k+1}=f'(z_{k})z'_{k}$ and $z'_{0}=1$ (before the calculation of the next iteration $z_{k+1}=f(z_{k})$ ). In the three cases we have:

\delta (z)=

lim_k→∞

|z_{kr}-z*|/|z'_{kr}|

(non-super-attraction)

\delta (z)=

lim_k→∞

log|z_{kr}-z*||z_{kr}-z*|/|z'_{kr}|

(super-attraction)

\delta (z)=

lim_k→∞

log|z_{k}||z_{k}|/|z'_{k}|

(d ≥ 2 and z* = ∞)

If this number (calculated for the last iteration number kr - to be divided by r) is smaller that a given small number, we colour the point z dark-blue, for instance.

For more see Pictures_of_Julia_and_Mandelbrot_Sets

code

/*
gcc -std=c99 -Wall -Wextra -pedantic -O3 -o julia-de julia-de.c -lm
https://math.stackexchange.com/questions/1153052/interior-distance-estimate-for-julia-sets-getting-rid-of-spots
code by Claude Heiland-Allen
*/

#include <complex.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>

void hsv2rgb(double h, double s, double v, int *rp, int *gp, int *bp) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *rp = fmin(fmax(round(r * 255), 0), 255);
  *gp = fmin(fmax(round(g * 255), 0), 255);
  *bp = fmin(fmax(round(b * 255), 0), 255);
}

complex double julia_attractor(complex double c, int maxiters, int *period) {
  double epsilon = nextafter(2, 4) - 2;
  complex double z = c;
  double mzp = 1.0/0.0;
  int p = 0;
  for (int n = 1; n < maxiters; ++n) {
    double mzn = cabs(z);
    if (mzn < mzp) {
      mzp = mzn;
      p = n;
      complex double z0 = z;
      for (int i = 0; i < 64; ++i) {
        complex double f = z0;
        complex double df = 1;
        for (int j = 0; j < p; ++j) {
          df = 2 * f * df;
          f = f * f + c;
        }
        complex double z1 = z0 - (f - z0) / (df - 1);
        if (cabs(z1 - z0) <= epsilon) {
          z0 = z1;
          break;
        }
        if (isinf(creal(z1)) || isinf(cimag(z1)) || isnan(creal(z1)) || isnan(cimag(z1))) {
          break;
        }
        z0 = z1;
      }
      complex double w = z0;
      complex double dw = 1;
      for (int i = 0; i < p; ++i) {
        dw = 2 * w * dw;
        w = w * w + c;
      }
      if (cabs(dw) <= 1) {
        *period = p;
        return z0;
      }
    }
    z = z * z + c;
  }
  *period = 0;
  return 0;
}

double julia_exterior_de(complex double c, complex double z, int maxiters, double escape_radius) {
  complex double dz = 1;
  for (int n = 0; n < maxiters; ++n) {
    if (cabs(z) >= escape_radius) {
      return cabs(z) * log(cabs(z)) / cabs(dz);
    }
    dz = 2 * z * dz;
    z = z * z + c;
  }
  return 0;
}

double julia_interior_de(complex double c, complex double z, int maxiters, double escape_radius, double pixel_size, complex double z0, int period, bool superattracting, int *fatou) {
  complex double dz = 1;
  for (int n = 0; n < maxiters; ++n) {
    if (cabs(z) >= escape_radius) {
      *fatou = -1;
      return cabs(z) * log(cabs(z)) / cabs(dz);
    }
    if (cabs(z - z0) <= pixel_size) {
      *fatou = n % period;
      if (superattracting) {
        return cabs(z - z0) * fabs(log(cabs(z - z0))) / cabs(dz);
      } else {
        return cabs(z - z0) / cabs(dz);
      }
    }
    dz = 2 * z * dz;
    z = z * z + c;
  }
  *fatou = -2;
  return 0;
}

int main(int argc, char **argv) {
  int size = 512;
  double radius = 2;
  double escape_radius = 1 << 10;
  int maxiters = 1 << 13;
  if (! (argc > 2)) { return 1; }
  complex double c = atof(argv[1]) + I * atof(argv[2]);

  int period = 0;
  bool superattracting = false;
  complex double z0 = julia_attractor(c, maxiters, &period);
  if (period > 0) {
    double epsilon = nextafter(1, 2) - 1;
    complex double z = z0;
    complex double dz = 1;
    for (int i = 0; i < period; ++i) {
      dz = 2 * z * dz;
      z = z * z + c;
    }
    superattracting = cabs(dz) <= epsilon;
  }

  double pixel_size = 2 * radius / size;
  printf("P6\n%d %d\n255\n", size, size);
  for (int j = 0; j < size; ++j) {
    for (int i = 0; i < size; ++i) {
      double x = 2 * radius * ((i + 0.5) / size - 0.5);
      double y = 2 * radius * (0.5 - (j + 0.5) / size);
      complex double z = x + I * y;
      double de = 0;
      int fatou = -1;
      if (period > 0) {
        de = julia_interior_de(c, z, maxiters, escape_radius, pixel_size, z0, period, superattracting, &fatou);
      } else {
        de = julia_exterior_de(c, z, maxiters, escape_radius);
      }
      int r, g, b;
      hsv2rgb(fatou / (double) period, 0.25 * (0 <= fatou), tanh(de / pixel_size), &r, &g, &b);
      putchar(r);
      putchar(g);
      putchar(b);
    }
  }
  return 0;
}

references

Fractals/Iterations in the complex plane
	demj	Iterations_in_the_complex_plane/def_cqp#Riemann_map →

[1] Koebe quarter theorem in wikipedia

[2] Distance Estimator by Robert P. Munafo

[3] th.stackexchange question: how-to-draw-a-mandelbrot-set-with-the-connecting-filaments-visible

[4] Pictures of Julia and Mandelbrot sets by Gert Buschmann

[5] Pictures of Julia and Mandelbrot sets by Gert Buschmann

[6] Julia sets by Jonas Lundgren in Matlab

[1]

[2]

[3]

[4]

[5]

[6]

Fractals/Iterations in the complex plane/demj

Contents

External distance estimation

Pseudocode and the code

How to use distance

tanh

distance max

Zoom

distance max

Examples of code

shadertoy

Internal distance estimation

Colouring the Julia set by Gert Buschmann

code

references