Mandelbrot 集的内部距离估计算法

Question

我正在寻找一种方法来估计从 Mandelbrot 集内部的一点到其边界的距离，以便在 GLSL 着色器中使用。

This page 链接到各种涉及内部距离估计主题的在线资源，例如底层数学公式、Haskell 实现、其他一些博客、论坛帖子和 C99 实现，但我的印象是它们都是要么实现起来非常复杂，要么运行起来计算量很大。

经过数小时的尝试，我设法使这段代码在 Shadertoy 中运行：

void mainImage( out vec4 fragColor, in vec2 fragCoord ) {

    float zoom = 1.;
    vec2 c = vec2(-0.75, 0.0) + zoom * (2.*fragCoord-iResolution.xy)/iResolution.y;

    vec2 z = c;

    float ar = 0.; // average of reciprocals
    float i;
    for (i = 0.; i < 1000.; i++) {
        ar += 1./length(z);
        z = vec2(z.x * z.x - z.y * z.y, 2.0 * z.x * z.y) + c;
    }
    ar = ar / i;

    fragColor = vec4(vec3(2. / ar), 1.0);
}

它确实会在每个灯泡中产生渐变，但很明显，它本身不能用作距离估计器，因为与较大灯泡相比，较小灯泡中的值具有不一致的幅度（亮度）。所以很明显缺少一个参数，但我不知道它是什么。

我不需要完美的解决方案，也不需要像this image 那样收敛成完美解决方案的解决方案。

至少可以保证下限的东西很多。

Answer 1

我敢打赌，1./length(z) 达到了float 的精度，如果有任何不同，请尝试使用double 和dvec2 而不是float,vec2。如果是这样，那么我会忽略太小的 length(z) 值。

或者您可以render 只在一次传递中将边界转换为纹理，然后just scan neighbors in all directions 直到找到返回光线长度的边界。 （安全使用前可能需要一些形态学算子）

这可以通过另一个通道来加快速度，在该通道中，您将增量距离“填充”到纹理中直到填充（最好在 CPU 端完成，因为您需要对同一纹理进行 R/W 访问）它类似于 A* filling 但您的精度会受到纹理分辨率的限制。

如果我将我的 mandlebrot 从上面的链接移植到你的计算移植到双精度并添加阈值：

// Fragment
#version 450 core
uniform dvec2 p0=vec2(0.0,0.0);     // mouse position <-1,+1>
uniform double zoom=1.000;          // zoom [-]
uniform int  n=100;                 // iterations [-]
in smooth vec2 p32;
out vec4 col;

vec3 spectral_color(float l)        // RGB <0,1> <- lambda l <400,700> [nm]
    {
    float t;  vec3 c=vec3(0.0,0.0,0.0);
         if ((l>=400.0)&&(l<410.0)) { t=(l-400.0)/(410.0-400.0); c.r=    +(0.33*t)-(0.20*t*t); }
    else if ((l>=410.0)&&(l<475.0)) { t=(l-410.0)/(475.0-410.0); c.r=0.14         -(0.13*t*t); }
    else if ((l>=545.0)&&(l<595.0)) { t=(l-545.0)/(595.0-545.0); c.r=    +(1.98*t)-(     t*t); }
    else if ((l>=595.0)&&(l<650.0)) { t=(l-595.0)/(650.0-595.0); c.r=0.98+(0.06*t)-(0.40*t*t); }
    else if ((l>=650.0)&&(l<700.0)) { t=(l-650.0)/(700.0-650.0); c.r=0.65-(0.84*t)+(0.20*t*t); }
         if ((l>=415.0)&&(l<475.0)) { t=(l-415.0)/(475.0-415.0); c.g=             +(0.80*t*t); }
    else if ((l>=475.0)&&(l<590.0)) { t=(l-475.0)/(590.0-475.0); c.g=0.8 +(0.76*t)-(0.80*t*t); }
    else if ((l>=585.0)&&(l<639.0)) { t=(l-585.0)/(639.0-585.0); c.g=0.84-(0.84*t)           ; }
         if ((l>=400.0)&&(l<475.0)) { t=(l-400.0)/(475.0-400.0); c.b=    +(2.20*t)-(1.50*t*t); }
    else if ((l>=475.0)&&(l<560.0)) { t=(l-475.0)/(560.0-475.0); c.b=0.7 -(     t)+(0.30*t*t); }
    return c;
    }

void main()
    {
    int i,j;
    dvec2 pp,p;
    double x,y,q,xx,yy,mu,cx,cy;
    p=dvec2(p32);

    pp=(p/zoom)-p0;         // y (-1.0, 1.0)
    pp.x-=0.5;              // x (-1.5, 0.5)
    cx=pp.x;                // normal
    cy=pp.y;
/*
    // single pass mandelbrot integer escape
    for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
        {
        q=xx-yy+cx;
        y=(2.0*x*y)+cy;
        x=q;
        xx=x*x;
        yy=y*y;
        }
    float f=float(i)/float(n);
    f=pow(f,0.2);
    col=vec4(spectral_color(400.0+(300.0*f)),1.0);
*/
    // distance to boundary
    double ar=0.0,aa,nn=0.0;              // *** this is what I added
    for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
        {
        aa=length(dvec2(x,y));            // *** this is what I added
        if (aa>1e-3){ ar+=1.0/aa; nn++; } // *** this is what I added
        q=xx-yy+cx;
        y=(2.0*x*y)+cy;
        x=q;
        xx=x*x;
        yy=y*y;
        }
    ar=ar/nn;                             // *** this is what I added
    col=vec4(vec3(1.0-(2.0/ar)),1.0);     // *** this is what I added
    }

我得到了这些输出：

只需在代码中查找 // *** this is what I added 注释，该注释是添加到标准 mandelbrot 渲染以渲染距离的内容。 ps 我的(x,y) 是你的z 和(cx,cy) 是你的c

无论如何，距离仍然是高度非线性的，并且取决于位置

[Edit1]非各向同性比例

黑点是阈值大小，您可以将其设置为1e-20 ...现在我添加了水平线以显示距离比例的分布（因为我不知道它是多么非各向同性和非线性... ) 这里的输出：

并对片段的一部分进行着色（在for 循环之后）：

ar=1.0-(2.0*nn/ar);
aa=10.0*ar; // 10 level lines per unit
aa-=floor(aa);
if (abs(aa)<0.05) col=vec4(0.0,1.0,0.0,1.0); // width and color of level line
 else             col=vec4(ar,ar,ar,1.0);

正如您所看到的，它与边界不是很平行，但仍然局部“恒定”（水平线与分形的局部特征中的每条线等距）所以如果使用梯度（导数），结果将只是非常粗略的估计（但应该管用）。如果这已经足够了，你应该做的是：

1. 计算查询位置的非线性距离，以及在“所有”方向上距离它的几个点d。

2. 选择与原点距离变化最大的邻居

3. 重新调整估计的距离，这样它们的减法就会给你d。然后使用第一个重新缩放的距离作为输出。

当放入片段代码时（使用 8-neighbors）：

// Fragment
#version 450 core
uniform dvec2 p0=vec2(0.0,0.0);     // mouse position <-1,+1>
uniform double zoom=1.000;          // zoom [-]
uniform int  n=100;                 // iterations [-]
in smooth vec2 p32;
out vec4 col;

double mandelbrot_distance(double cx,double cy)
    {
    // distance to boundary
    int i,j;
    double x,y,q,xx,yy,ar=0.0,aa,nn=0.0;
    for (x=0.0,y=0.0,xx=0.0,yy=0.0,i=0;(i<n)&&(xx+yy<4.0);i++)
        {
        aa=length(dvec2(x,y));
        if (aa>1e-20){ ar+=1.0/aa; nn++; }
        q=xx-yy+cx;
        y=(2.0*x*y)+cy;
        x=q;
        xx=x*x;
        yy=y*y;
        }
    return 1.0-(2.0*nn/ar);
    }
void main()
    {
    dvec2 pp,p;
    double cx,cy,d,dd,d0,d1,e;
    p=dvec2(p32);

    pp=(p/zoom)-p0;         // y (-1.0, 1.0)
    pp.x-=0.5;              // x (-1.5, 0.5)
    cx=pp.x;                // normal
    cy=pp.y;
    d =0.01/zoom;           // normalization distance
    e =sqrt(0.5)*d;
    dd=mandelbrot_distance(cx,cy);
    if (dd>0.0)
        {
        d0=mandelbrot_distance(cx-d,cy  ); if (d0>0.0)  d0=abs(d0-dd);
        d1=mandelbrot_distance(cx+d,cy  ); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx  ,cy-d); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx  ,cy+d); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx-e,cy-e); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx+e,cy-e); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx-e,cy+e); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        d1=mandelbrot_distance(cx+e,cy+e); if (d1>0.0){ d1=abs(d1-dd); if (d0<d1) d0=d1; }
        dd*=d/d0;
        }
    dd*=zoom; // just for visualization of small details real distance should not be scaled by this
    col=vec4(dd,dd,dd,1.0);
    }

结果如下：

正如您所看到的，它现在更加正确（但由于上述非各向同性，非常接近边界是不准确的）。 8 个邻居在圆形斑点中产生 8 个类似对角线的线条图案。如果您想摆脱它们，您应该围绕该位置扫描整个圆圈，而不是仅扫描 8 个点。

还有一些白点（它们不准确）我认为它们是在Selected d遥远邻居跨越Mandelbrot边缘的情况而不是原版。那可以被过滤掉......（你知道d/2在同一方向的距离应该是一半，如果不是你在不同的blob）

但是即使是 8 个邻居也很慢。因此，为了获得更高的准确性，我建议改用 2 遍“光线投射”方法。