此算法有 2 个版本

外部
内部

与动态平面和 Julia 集版本：DEM/J 进行比较

外部 DEM/M

曼德勃罗集的距离估计方法 (DEM/M) 估计点 $c\,$ （在曼德勃罗集的外部）到曼德勃罗集最近点的距离。^[1]^[2]

名称

（方向）距离估计公式

它可用于从 BOF 创建黑白图像：^[3]

第 84 页的 41 图
第 85 页的 43 图
第 188 页的无编号图

类型

解析 DE = DEM = 真实 DE
Claude DE = 假 DE (FDE)：当 g>0 时 FDE=1/(g*log(2))，当 g=0 时未定义（即在内部）^[4]

数学公式

数学公式：^[5]^[6] ^[7]

distance(c_{0},\sigma M)=\lim _{n\to \infty }2{\frac {|z_{n}|}{|z'_{n}|}}\ln |z_{n}|

其中

z_{n}=f_{c}^{n}(z=0.0)

z_{n}'={\frac {d}{dc}}f_{c}^{n}(z_{0}).

是 $f_{c}^{n}(z_{0})$ 关于 c 的一阶导数，它可以通过从以下开始的迭代找到

z_{0}'={\frac {d}{dc}}f_{c}^{0}(z_{0})=1

然后在每个连续步骤中替换

z_{n+1}'={\frac {d}{dc}}f_{c}^{n+1}(z_{0})=2\cdot {}f_{c}^{n}(z)\cdot {\frac {d}{dc}}f_{c}^{n}(z_{0})+1=2\cdot z_{n}\cdot z_{n}'+1.

逃逸半径

"你需要将“If mag > 2”更改为增加逃逸半径。附件是一个比较，显示当半径增加时，“带子”会减少，没有问题可以将其增加到像 R = 1e6 这样的大小，因为一旦它逃逸了 2，它就会快速增长。半径越大，距离估计越准确。” Claude ^[8]
基本上，在测试幅度平方与逃逸值时
- 如果你的逃逸值是 10^4，那么你的 DE 估计值精确到两位小数，例如，如果你的 DE 是 0.01，那么它只精确到 +/-0.0001，
- 如果你的逃逸值是 10^6，那么你的 DE 估计值精确到三位小数（David Makin 在 fractalforums.com 上说）

"对于外部距离估计，你需要一个大的逃逸半径，例如 100100。迭代次数限制是任意的，在有限的限制下，某些像素将始终被归类为“未知”。" Claude Heiland-Allen^[9]

比较

分形的美丽
分形图像科学，第 198 页，
由 Robert P. Munafo 编写的距离估计器^[10]

由 Claude Heiland-Allen 编写的伪代码^[11]

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

代码 ^[12]

      double _Complex m_exterior_distance(int N, double R, double _Complex c)
{
  double _Complex z = 0;
  double _Complex dc = 0;
  for (int n = 0; n < N; ++n)
  {
    if (cabs(z) > R)
      return 2 * z * log(cabs(z)) / dc;
    dc = 2 * z * dc + 1;
    z = z * z + c;
  }
  return -1;
}

来自 BOF 的代码分析

"分形的美丽" 为右侧图像中显示的距离估计提供了几乎正确的计算机程序。该方法没有获得支持的一个可能原因是该程序中的过程存在严重缺陷： $z_{k}$ 的计算是在 $z'_{k}$ 计算之前（并完成），而不是之后，它应该是在之后（ $z'_{k+1}$ 使用 $z_{k}$ ，而不是 $z_{k+1}$ ）。为了避免重新计算 $z_{k}$ （k = 0, 1, 2, ...），此序列被保存在一个数组中。使用此数组， $z'_{k+1}=2z_{k}z'_{k}+1$ 被计算到最后一个迭代次数，并指出可能发生溢出。如果发生溢出，则该点被认为属于边界（逃逸条件）。如果未发生溢出，则可以执行距离计算。除了发生溢出并非事实之外，该方法使用了不必要的存储和迭代重复，使其变得不必要地更慢且吸引力更低。本书中的以下评论也不是令人鼓舞的：“事实证明，图像对各种选择非常敏感”（逃逸半径、最大迭代次数、溢出、边界厚度和放大系数）。是这种无稽之谈让人们对使用和推广该方法失去了所有兴趣吗？" Gertbuschmann

所以

在循环中，导数应该在新的 z 之前计算

"我终于生成了让我满意的 DEM（距离估计方法）图像。事实证明，我的代码中存在一些错误。新的代码运行速度相当快，即使计算距离估计值需要额外的数学运算。

S of F I（“分形图像科学”）中的算法使用从白色到黑色的锐利截止，当像素足够靠近集合时，但我发现这会导致图像看起来参差不齐。我已经实现了一个非线性颜色渐变，效果很好。

对于 DE 值为 threshold>=DE>=0 的像素，我将 DE 值缩放到 1>=DE>=0，然后进行计算：gradient_value = 1 - (1-DE)^n，（“^n”表示 n 次方。我希望我能使用上标！）

"n" 是一个调整因子，它允许我更改渐变曲线的形状。值“gradient_value”决定像素的颜色。如果它为 0，则像素以起始颜色（对于黑白图像，为白色）着色。如果“gradient_value”为 1，则像素获得结束颜色（例如，

黑色) 对于 n>1 的值，这将导致远离集合的像素颜色发生快速变化，并且随着 DE 接近 0，值的变化速度会变慢。对于 n<1 的值，它会导致远离集合的像素颜色发生缓慢变化，并且靠近集合的像素颜色发生快速变化。对于 n=1，我得到一个线性

颜色梯度。非线性梯度允许我使用 DE 值来反锯齿我的绘图。通过调整我的阈值和调整值，我可以为各种图像获得良好的效果。”Duncan C ^[13]

“我们上面列出的所有 DE 公式都只是近似值——在 n→∞ 的极限情况下有效，而且一些公式也只适用于靠近分形边界点的点。当您开始渲染这些结构时，这一点会变得非常明显——您通常会遇到噪点和伪影。

将 DE 估计值乘以小于 1 的数字可以用来减少噪点（这就是 Fragmentarium 中的“ Fudge Factor”）。

另一种常见的方法是 **过采样**——或以较大的尺寸渲染图像并缩小尺寸。”Mikael Hvidtfeldt Christensen ^[14]

示例代码

ultrafractal

两种算法，在两个循环中

Inigo Quilez 在 shadertoy : Distance estimation to the Mandelbrot set by Inigo Quilez 中提供的曼德布罗特集合距离估计方法

这是一个 Java 函数。它在一个循环中同时计算 : 迭代 $z_{n}\,$ 和导数 $z'_{n}\,$ 。如果（在动态平面中）临界点

没有逃逸到无穷远（跳出条件），那么它属于曼德布罗特集合的内部，并具有最大颜色值，
逃逸则它可能在外部或边界附近。它的颜色与点 c 与曼德布罗特集合中最近点的距离“成比例”。它也使用颜色循环（（int）R % maxColor）。

// Java function by Evgeny Demidov from http://www.ibiblio.org/e-notes/MSet/DEstim.htm
// based on code by Robert P. Munafo from http://www.mrob.com/pub/muency/distanceestimator.html
  public int iterate(double cr, double ci, double K, double f) {
    double Cr,Ci, I=0, R=0,  I2=I*I, R2=R*R, Dr=0,Di=0,D;   int n=0;
    if (f == 0){ Cr = cr; Ci = ci;}
    else{ Cr = ci; Ci = cr;}
    do  {
      D = 2*(R*Dr - I*Di) +1;  Di = 2*(R*Di + I*Dr);  Dr = D;
      I=(R+R)*I+Ci;  R=R2-I2+Cr;  R2=R*R;  I2=I*I;  n++;
    } while ((R2+I2 < 100.) && (n < MaxIt) );
    if (n == MaxIt) return maxColor; // interior
    else{ // boundary and exterior 
     R = -K*Math.log(Math.log(R2+I2)*Math.sqrt((R2+I2)/(Dr*Dr+Di*Di))); // compute distance
     if (R < 0) R=0;
     return (int)R % maxColor; }; 
 }

这是一个 cpp 函数。它将颜色的整数索引作为输出。

/*
 this function is  from program mandel ver 5.10 by Wolf Jung
 see file mndlbrot.cpp   
 http://www.mndynamics.com/indexp.html
 
It is checked first (in pixcolor)
that the point escapes before calling this function.  
So we do not compute the derivative and the logarithm 
when it is not escaping.  
This would not be a good idea if most points escape anyway.

*/
int mndlbrot::dist(double a, double b, double x, double y)
{  uint j; 
   double xp = 1, yp = 0; // zp = xp+yp*i =  1 ; derivative 
   double nz, nzp; 
   
  // iteration 
  for (j = 1; j <= maxiter; j++)
   {  // zp
      nz = 2*(x*xp - y*yp); 
      yp = 2*(x*yp + y*xp); 
      xp = nz; //zp = 2*z*zp; on the dynamic plane 
      // if sign is positive it is parameter plane, if negative it is dynamic plane.
      if (sign > 0) xp++; //zp = 2*z*zp + 1 ; on the parameter plane
      //z = z*z + c; 
      nz = x*x - y*y + a; 
      y = 2*x*y + b; 
      x = nz; 
      //
      nz = x*x + y*y; 
      nzp = xp*xp + yp*yp;
      // 2 conditions for stopping the iterations
      if (nzp > 1e40 || nz > bailout) break;
   }

   // 5 types of points but 3 colors  
   /*  not escaping, rare 
       Is it not simply interior ???
       This should not be necessary.  I do not know why I included it,
       because in this case  pixcolor should not call dist.  If you do not
       have pixcolor before,  you should return 0 here. */
   if (nz < bailout) return 1; // not escaping, rare
   /*  If The Julia set is disconnected and the orbit of z comes close to
    z = 0 before escaping,  nzp will be small  */
   if (nzp < nz) return 10; // includes escaping through 0
    
   // compute estimated distance  = x 
   x = log(nz); 
   x = x*x*nz / (temp[1]*nzp); //4*square of dist/pixelwidth
   if (x < 0.04) return 1; // exterior but very close to the boundary
   if (x < 0.24) return 9; // exterior but close to the boundary
   return 10; //exterior : escaping and not close to the boundary
} //dist

两种算法，在一个循环中

ińigo quilez 在 c++ 中提供的分形距离渲染方法^[15]

Function GiveDistance(xy2,eDx,eDy:extended):extended;

begin
    result:=2*log2(sqrt(xy2))*sqrt(xy2)/sqrt(sqr(eDx)+sqr(eDy));
  end;

//------------------------------------

Function PointIsInCardioid (Cx,Cy:extended):boolean;
 //Hugh Allen
 // http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip
  var DeltaX,DeltaY:extended;
      //
      PDeltyX,PDeltyY:extended;
      //
      ZFixedX,ZFixedY:extended;

  begin
     result:=false;
     // cardioid checkig - thx to Hugh Allen
     //sprawdzamy Czy punkt  C jest w głównej kardioidzie
     //Cardioid in squere :[-0.75,0.4] x [ -0.65,0.65]
     if InRange(Cx,-0.75,0.4)and InRange(Cy,-0.65,065) then
      begin
        // M1= all C for which Fc(z) has  attractive( = stable) fixed point
        // znajdyjemy punkt staly z: z=z*z+c
        // czyli rozwiazujemy rownanie kwadratowe
        // zmiennej zespolonej o wspolczynnikach zespolonych
        // Z*Z - Z + C = 0
        //Delta:=1-4*a*c; Delta i C sa liczbami zespolonymi
        DeltaX:=1-4*Cx;
        DeltaY:=-4*Cy;
        // Pierwiastek zespolony z delty
        CmplxSqrt(DeltaX,DeltaY,PDeltyX,PDeltyY);
        // obliczmy punkt staly jeden z dwóch, ten jest prawdopodobnie przycišgajšcy
        ZFixedX:=1.0-PDeltyX; //0.5-0.5*PDeltyX;
        ZFixedY:=PDeltyY; //-0.5*PDeltyY;
        // jesli punkt stały jest przycišgajšcy
        // to należy do M1
        If (ZfixedX*ZFixedX + ZFixedY*ZFixedY)<1.0
          then result:=true;

         
          // ominięcie iteracji M1 przyspiesza 3500 do 1500 msek
        end; // if InRange(Cx ...

  end;
//------------------------------------
Function PointIsInComponent (Cx,Cy:extended):boolean;
//Hugh Allen
// http://homepages.borland.com/ccalvert/Contest/MarchContest/Fractal/Source/HughAllen.zip

  var Dx:extended;
  begin
    result:=false;
    // czy punkt C nalezy do koła na lewo od kardioidy
    // circle: center = -1.0  and radius 1/4
    dx:=Cx+1.0;
    if (Dx*Dx+Cy*Cy) < 0.0625
      then result:=true;

  end;

//------------------------------
Procedure DrawDEM_DazibaoTrueColor; 
// draws Mandelbrot set in black and its complement in true color
// see   http://ibiblio.org/e-notes/MSet/DEstim.htm
// by  Evgeny Demidov
//
// see also
//http://www.mandelbrot-dazibao.com/Mset/Mset.htm
// translation ( with modification) of Q-Basic program:
//  http://www.mandelbrot-dazibao.com/Mset/Mdb3.bas
//
// see also my page http://republika.pl/fraktal/mset_dem.html

var iter:integer;
     iY,iX:integer;
     eCy ,eCx:extended; // C:=eCx + eCy*i
     eX,eY:extended;    // Zn:=eX+eY*i
     eTempX,eTempY:extended;
     eX2,eY2:extended;  //x2:=eX*eX;  y2:=eY*eY;
     eXY2:extended;  //  xy2:=x2+y2;
     eXY4:extended;
     eTempDx,eTempDy:extended;
     eDx,eDy:extended; // derivative
     distance:extended;
     color:TColor;

begin
  //compute bitmap
  for iY:= iYmin to iYMax do
    begin
      eCy:=Convert_iY_to_eY(iY);
      for iX:= iXmin to iXmax do
        begin
          eCx:=Convert_iX_to_eX(iX);
          If not PointIsInCardioid (eCx,eCy) and not PointIsInComponent(eCx,eCy)
            then
              begin
                //  Z0:=0+0*i
                eX:=0;
                eY:=0;
                eTempX:=0;
                eTempY:=0;
                //
                eX2:=0;
                eY2:=0;
                eXY2:=0;
                //
                eDx:=0;
                eDy:=0;
                eTempDx:=0;
                eTempDy:=0;
                //
                iter:=0;
                // iteration of Z ; Z= Z*z +c
                while ((iter<IterationMax) and (eXY2<=BailOut2)) do
                  begin
                    inc(iter);
                    //
                    eTempY:=2*eX*eY + eCy;
                    eTempX:=eX2-eY2 + eCx;
                    //
                    eX2:=eTempX*eTempX;
                    eY2:=eTempY*eTempY;
                    //
                    eTempDx:=1+2*(eX*eDx-eY*eDy);
                    eTempDy:=2*(eX*eDY+eY*eDx);
                    //
                    eXY2:=eX2+eY2;
                    //
                    eX:=ETempX;
                    eY:=eTempY;
                    //
                    eDx:=eTempDx;
                    eDy:=eTempDy;
                  end; // while

                // drawing procedure
                if (iter<IterationMax)
                then
                  begin
                    distance:= GiveDistance(eXY2,eDx,eDy);
                    color:=Rainbow(1,500,Abs(Round(100*Log10(distance)) mod 500));
                    with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := GetBValue(color);
                          G := GetGValue(color);
                          R := GetRValue(color);
                          //A := 0;
                        end; // with FirstLine[Y*LineLength+X]

                end // if (iter<IterationMax) then
              else  with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := 0;
                          G := 0;
                          R := 0;
                          //A := 0;
                        end;
             //--- end of drawing procedure ---
            end //  If not PointIsInCardioid ... then
          else with Bitmap1.FirstLine[iY*Bitmap1.LineLength+iX] do
                        begin
                          B := 0;
                          G := 0;
                          R := 0;
                          //A := 0;
                        end;
         //--- If not PointIsInCardioid ...
        end; // for iX
      end; // for iY

end;
//------------------------------------------------------

修改/优化

可以通过使用来改进 DEM 图像的创建

柯比 1/4 定理^[16]
基于距离的灰度图^[17]^[18]
分形论坛 : 跳出条件与解析 DE 精度之间的关系
DEM/M 中的权重
如果你看到主天线有双重轮廓
- 提高分辨率并不能解决问题
- 在高度上增加一个像素，或更改 CyMax 或 CyMin（增加一个像素大小）。

均衡

距离估计均衡^[19]

程序 exrdeeq^[20]

读取 DEX 和 DEY 通道
输出 Y 通道
DE 小于 0 则变为 0
DE 大于 1 则变为 1
否则 Y 为 DE 在排序数组中的索引，缩放到 0 到 1 的范围内

exrdeeq in.exr out.exr

示例

Duncan Champney 绘制的 B of F 地图 43 DEM^[21]
该图以 -0.7436636773584340, 0.1318632143960234i 为中心，宽度约为 6.25E-11
Duncan Champney 绘制的微笑万花筒 ^[22]

height=800 
width=800 
max_iterations=10000
center_r=-1.74920463346
center_i=-0.000286846601561
r_width=3.19E-10

内部距离估计

估计极限周期点 $c\,$ （在曼德布罗特集合的内部）到曼德布罗特集合边界的距离。

描述 :^[23]

书籍 : The Science of Fractal Images, Springer Verlag, Tokyo, Springer, New York, 1988 by Heinz-Otto Peitgen , D. Saupe, page 294
Albert Lobo Cusidó 关于曼德布罗特集合的内部和外部距离边界
R Munafo 的内部距离估计方法
mandelbrot-numerics 库 : m_d_interior by Claude

‎

// https://mathr.co.uk/mandelbrot/book-draft/#interior-distance
// Claude Heiland-Allen
#include <complex.h>
#include <math.h>

double cnorm(double _Complex z)
{
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

double m_interior_distance
    (double _Complex z0, double _Complex c, int p)
{
    double _Complex z = z0;
    double _Complex dz = 1;
    double _Complex dzdz = 0;
    double _Complex dc = 0;
    double _Complex dcdz = 0;
    for (int m = 0; m < p; ++m)
    {
        dcdz = 2 * (z * dcdz + dz * dc);
        dc = 2 * z * dc + 1;
        dzdz = 2 * (dz * dz + z * dzdz);
        dz = 2 * z * dz;
        z = z * z + c;
    }
    return (1 - cnorm(dz))
        / cabs(dcdz + dzdz * dc / (1 - dz));
}

double m_distance(int N, double R, double _Complex c)
{
    double _Complex dc = 0;
    double _Complex z = 0;
    double m = 1.0 / 0.0;
    int p = 0;
    for (int n = 1; n <= N; ++n)
    {
        dc = 2 * z * dc + 1;
        z = z * z + c;
        if (cabs(z) > R)
            return 2 * cabs(z) * log(cabs(z)) / cabs(dc);
        if (cabs(z) < m)
        {
            m = cabs(z);
            p = n;
            double _Complex z0 = m_attractor(z, c, p);
            double _Complex w = z0;
            double _Complex dw = 1;
            for (int k = 0; k < p; ++k)
            {
                dw = 2 * w * dw;
                w = w * w + c;
            }
            if (cabs(dw) <= 1)
                return m_interior_distance(z0, c, p);
        }
    }
    return 0;
}

数学描述

该估计值由给出

distance(c,\sigma M)=d(c,z_{0},p)={\frac {1-\mid {{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}\mid ^{2}}{\mid {{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})+{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0}){\frac {{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})}{1-{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})}}}\mid }}

其中

$p$ 是周期 = 周期轨道的长度

$c$ 是要估计的参数平面中的点

$f_{c}(z)$ 是复二次多项式 $f_{c}(z)=z^{2}+c$

$f_{c}^{p}(z_{0})$ 是 $p$ -次迭代 $f_{c}(z)$

$z_{0}$ 是构成周期轨道（极限环） $\{z_{0},\dots ,z_{p-1}\}$ 的 $p$ 个点中的任意一个

在 $z_{0},c,p$ 处计算的 $f_{c}^{p}$ 的 4 个导数

${\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})$

${\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$

关于 z 的一阶偏导数

它必须通过应用链式法则递归地计算

起点： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

第一次迭代： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

第二次迭代： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

第p次迭代的一般规则： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

关于c的一阶偏导数

它必须通过应用链式法则递归地计算

起始点： ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

第一次迭代： ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{0}(z_{0})+1=1\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

第二次迭代： ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{1}(z_{0})+1=2z_{0}^{2}+2c+1\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

p 次迭代的一般规则： ${\begin{cases}{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}f_{c}^{p}(z_{0})+1\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

关于 z 的二阶偏导数

它必须计算

与： $f_{c}^{n}$ 和 ${\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})$ 一起
递归地应用链式法则

起点： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=1\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=0\\f_{c}^{0}(z_{0})=z_{0}\end{cases}}$

第一次迭代： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=2*f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})=2z_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0}))^{2}+f_{c}^{0}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{0}(z_{0})\}=2\\f_{c}^{1}(z_{0})=(f_{c}^{0}(z_{0}))^{2}+c=z_{0}^{2}+c\end{cases}}$

第二次迭代： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{2}(z_{0})=2*f_{c}^{1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{1}(z_{0})=4z_{0}^{3}+4cz_{0}\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=\\f_{c}^{2}(z_{0})=(f_{c}^{1}(z_{0}))^{2}+c=(z_{0}^{2}+c)^{2}+c=z_{0}^{4}+2*c*z_{0}^{2}+c^{2}+c\end{cases}}$

p 迭代的一般规则： ${\begin{cases}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2*f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\\{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2\{({\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))^{2}+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})\}\\f_{c}^{p}(z_{0})=f_{c}^{p-1}(z_{0})^{2}+c\end{cases}}$

二阶混合偏导数

${\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p}(z_{0})=2({\frac {\partial }{\partial {c}}}f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0})+f_{c}^{p-1}(z_{0})*{\frac {\partial }{\partial {c}}}{\frac {\partial }{\partial {z}}}f_{c}^{p-1}(z_{0}))$

算法

步骤

选择 c
检查给定 c 的临界点是否为周期性的（内部）或非周期性的（外部或接近边界或周期过大）
- 如果非周期性，则不要使用此算法（使用外部版本）
- 如果周期性，则计算周期和周期点 $z_{0}$
使用 $c,z_{0},p$ 计算 p 迭代的距离

计算周期和周期点

计算距离

代码

Albert Lobo Cusidó 的 Java 代码^[24]
Claude 的 Haskell 代码^[25] 和图像^[26]
tit_toinou 的 processing 代码^[27] 和描述^[28]

问题

内部距离估计有两个实际问题

首先，我们需要精确找到 $z_{0}$ ，
其次，我们需要精确找到 $p$ 。

$z_{0}$ 的问题是，理论上，通过迭代 $P_{c}(z)$ 来收敛到 $z_{0}$ 需要无限次的操作。

周期的问题是，有时由于舍入误差，会错误地将周期识别为真实周期的整数倍（例如，检测到 86 的周期，而真实周期只有 43=86/2）。在这种情况下，距离被高估，即报告的半径可能包含曼德尔布罗特集合之外的点。

对于内部距离估计，您需要周期，然后需要一些（可能 10 个左右，通常就足够了）牛顿步来找到极限环。

内部检查代码绝对要求参考点位于岛的核心中（不是它的任何子圆盘，当然也不是随机点）

优化

与外部情况类似，一旦找到b，我们就知道距离c 的b/4 内的所有点都在曼德尔布罗特集合内。
使用距离估计的自适应超采样 Claude Heiland-Allen 使用距离估计的自适应超采样

参考资料

Claude Heiland-Allen 2020-04-23 的混合逃逸时间分形距离估计
Milnor "曼德尔布罗特集合中的自相似性和毛发"，1989 年的几何与拓扑学中的计算机。

↑ A Cheritat wiki-draw: Mandelbrot_set
↑ fractalforums : m-set-distance-estimation
↑ fractalforums discussion : How are B&W images from "Beauty of Fractals" created?
↑ fractalforums.org : m-brot-distance-estimation-versus-claudes-fake-de
↑ Milnor (在单个复变量中的动力学，附录 G)
↑ Heinz-Otto Peitgen（编辑、贡献者），Dietmar Saupe（编辑、贡献者），Yuval Fisher（贡献者），Michael McGuire（贡献者），Richard F. Voss（贡献者），Michael F. Barnsley（贡献者），Robert L. Devaney（贡献者），Benoit B. Mandelbrot（贡献者） : 分形图像的科学。施普林格；第一版（1988 年 7 月 19 日），第 199 页
↑ ińigo quilez 的分形距离渲染
↑ fractalforums : mandelbrot-distance-estimation-problem
↑ fractalforums.org : list-of-numbers-with-fixed-digit-length-in-the-mandelbrot-set
↑ Robert P. Munafo 的距离估计器
↑ math.stackexchange 问题：如何绘制具有连接细丝的曼德尔布罗特集合
↑ 曼德尔布罗特书籍
↑ Fractal forum : 如何从 "Beauty of Fractals" 中创建黑白图像？
↑ Mikael Hvidtfeldt Christensen 的距离估计 3D 分形（V）：曼德尔球体和不同的 DE 近似
↑ ińigo quilez 的分形距离渲染
↑ Claude Heiland-Allen 的距离估计
↑ ińigo quilez 的分形距离渲染
↑ fractalforums 画廊，作者 Pauldelbrot
↑ fractalforums.org : histogram-de-coloring
↑ Claude Heiland-Allen 的距离估计均衡
↑ Duncan Champney 的 B of F 地图 43 DEM
↑ Duncan Champney 的 Smily Kaleidoscope
↑ 曼德博集合的内部和外部距离界限，作者 Albert Lobo Cusidó
↑ 曼德博集合的内部和外部距离界限，作者 Albert Lobo Cusidó
↑ 曼德博集合中的内部坐标，作者 Claude
↑ 分形论坛 : 曼德博集合内部点的着色
↑ ttoinou 编写的 Processing 中的 MandelbrotDE
↑ 分形论坛 : 使用距离和梯度着色的经典曼德博集合

[1] A Cheritat wiki-draw: Mandelbrot_set

[2] ractalforums : m-set-distance-estimation

[3] ractalforums discussion : How are B&W images from "Beauty of Fractals" created?

[4] ractalforums.org : m-brot-distance-estimation-versus-claudes-fake-de

[5] Milnor (在单个复变量中的动力学，附录 G)

[6] Heinz-Otto Peitgen（编辑、贡献者），Dietmar Saupe（编辑、贡献者），Yuval Fisher（贡献者），Michael McGuire（贡献者），Richard F. Voss（贡献者），Michael F. Barnsley（贡献者），Robert L. Devaney（贡献者），Benoit B. Mandelbrot（贡献者） : 分形图像的科学。施普林格；第一版（1988 年 7 月 19 日），第 199 页

[7] ńigo quilez 的分形距离渲染

[8] ractalforums : mandelbrot-distance-estimation-problem

[9] ractalforums.org : list-of-numbers-with-fixed-digit-length-in-the-mandelbrot-set

[10] Robert P. Munafo 的距离估计器

[11] th.stackexchange 问题：如何绘制具有连接细丝的曼德尔布罗特集合

[12] 曼德尔布罗特书籍

[13] Fractal forum : 如何从 "Beauty of Fractals" 中创建黑白图像？

[14] Mikael Hvidtfeldt Christensen 的距离估计 3D 分形（V）：曼德尔球体和不同的 DE 近似

[15] ńigo quilez 的分形距离渲染

[16] Claude Heiland-Allen 的距离估计

[17] ńigo quilez 的分形距离渲染

[18] ractalforums 画廊，作者 Pauldelbrot

[19] ractalforums.org : histogram-de-coloring

[20] Claude Heiland-Allen 的距离估计均衡

[21] Duncan Champney 的 B of F 地图 43 DEM

[22] Duncan Champney 的 Smily Kaleidoscope

[23] 曼德博集合的内部和外部距离界限，作者 Albert Lobo Cusidó

[24] 曼德博集合的内部和外部距离界限，作者 Albert Lobo Cusidó

[25] 曼德博集合中的内部坐标，作者 Claude

[26] 分形论坛 : 曼德博集合内部点的着色

[27] ttoinou 编写的 Processing 中的 MandelbrotDE

[28] 分形论坛 : 使用距离和梯度着色的经典曼德博集合

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]