此演算法有兩個版本

外部距離估計

Julia 集的距離估計方法 ( DEM/J ) 估計點 z ( 在 Julia 集外部 ) 到 Julia 集中最近點的距離。

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

對於距離估計，已證明計算值與真實距離最多相差 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]

數學公式

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

其中

$z'_{n}\,$ 是關於 z 的一階導數.

這個導數可以透過迭代找到，從

z'_{0}=1\,

然後

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

虛擬碼和程式碼

分形的美麗
分形影像的科學，第 198 頁，
羅伯特·P·穆納福的距離估計器^[2]

克勞德·海蘭德-艾倫的虛擬碼^[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))

關於 z 的第一個導數
度數	函式 f(z)	關於 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://wikibook.tw/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;
}

如何使用距離

tanh

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

最大距離

Distance max 是舊的（已棄用）方法。

可以使用距離來進行著色

僅限於 Julia 集（填充 Julia 集的邊界）
填充 Julia 集的邊界和外部。

這裡第一個例子

 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 */}
          }

這裡第二個例子 ^[4]

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

縮放

最大距離

Distance max 是舊的（已棄用）方法。

DistanceMax 小於畫素大小。在縮放時，畫素大小會減小，DistanceMax 也應該減小以獲得好的圖片。可以使用以下公式進行操作：

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

其中 $n\leq 1\,$

可以從 n=1 開始，如果圖片不好，可以增加 n。還要檢查 iMax !!

DistanceMax 也可能與縮放因子 $mag\,$ 成正比：^[5]

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

其中 thick 是影像寬度（以世界單位計），mag 是縮放因子。

還可以使用 tanh，它可以提供更精確的外觀

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

程式碼示例

有關 cpp 示例，請參見 src 程式碼中 mndlbrot.cpp 中的 mndlbrot::dist，該程式碼來自程式 mandel ver 5.3 by Wolf Jung。

使用複數型別的 C 函式

unsigned char ComputeColorOfDEMJ(complex double z){
// https://wikibook.tw/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 函式

/*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 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

Jonas Lundgren 編寫的 Matlab 程式碼^[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 函式

 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。到 Julia 集 z^2 + c 的解析距離，其中 vec2 c = vec2( -0.745, 0.186 )
Julia - Distance 2 by iq = f(z) = z^3+C 的 Julia 集的 SDF，其中 const vec2 kC = vec2(0.105,0.7905);
Julia - Distance 3 by iq。一個通用的 Julia 集渲染器，使用到該集合的距離（Douady_Hubbard 勢）。在本例中，我使用的是 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;

內部距離估計

Gert Buschmann 對 Julia 集的著色

為了獲得不錯的圖片，我們還必須對 Julia 集進行著色，因為否則 Julia 集僅透過對 Fatou 域的著色才能看到，而這種著色在 Julia 集附近會發生劇烈變化，導致外觀渾濁（可以透過仔細選擇顏色比例和密度來避免這種情況）。如果迭代沒有停止，則點z 屬於 Julia 集，也就是說，如果我們已經達到選擇的最大迭代次數 M。但是由於 Julia 集無限薄，而且我們僅對規律排列的點進行計算，在實際中我們無法透過這種方式對 Julia 集進行著色。但幸運的是，存在一個公式可以（直到一個常數因子）估計點z（位於 Julia 集外部）到 Julia 集的距離。該公式與 Fatou 域相關聯，並且該公式給出的數字越靠近 Julia 集越正確，因此偏差無關緊要。

距離函式 是函式 $\delta (z)=\phi (z)/|\phi '(z)|$ （參見非複數函式的 Julia 集和 Mandelbrot 集部分），其等勢線必須近似規則排列。公式中出現了 $z'_{k}$ 相對於z 的導數。但由於 $z_{k}=f(f(...f(z)))$ （k 重複合）， $z'_{k}$ 是 $f'(z_{i})$ （i = 0, 1, ..., k-1）的乘積，這個序列可以透過 $z'_{k+1}=f'(z_{k})z'_{k}$ 和 $z'_{0}=1$ （在計算下一個迭代 $z_{k+1}=f(z_{k})$ 之前）進行遞迴計算。在三種情況下，我們有

\delta (z)=

lim_k→∞

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

(非超吸引)

\delta (z)=

lim_k→∞

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

(超吸引)

\delta (z)=

lim_k→∞

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

(d ≥ 2 且 z* = ∞)

如果這個數字（根據最後一次迭代次數kr計算，併除以r）小於給定的一個小數，那麼我們就可以將點z染成深藍色。

更多資訊請參見Pictures_of_Julia_and_Mandelbrot_Sets

程式碼

/*
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;
}

參考文獻

分形/複平面上的迭代
	demj	Iterations_in_the_complex_plane/def_cqp#Riemann_map →

[1] 維基百科中的Koebe 四分之一定理

[2] Robert P. Munafo 的距離估計器

[3] th.stackexchange 問題：如何繪製具有連線細絲的曼德勃羅集合

[4] Gert Buschmann 繪製的朱利亞和曼德勃羅集影像

[5] Gert Buschmann 繪製的朱利亞和曼德勃羅集影像

[6] Jonas Lundgren 在 Matlab 中編寫的朱利亞集合

[1]

[2]

[3]

[4]

[5]

[6]