此演算法有 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]