问题描述
我正在寻找一种方法来估计 Mandelbrot 集边界与 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 那样收敛成一个完美的解决方案。
至少能保证下限的东西就足够了。
解决方法
我敢打赌，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,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,1.0); // width and color of level line
else col=vec4(ar,ar,1.0);
正如你所看到的，它与边界不是很平行，但仍然是局部“恒定的”（水平线与分形局部特征中的每条线等距）所以如果使用梯度（导数），结果将只是非常粗略的估计（但应该管用）。如果这足够了，您应该做的是：
-
计算查询位置的非线性距离以及在“所有”方向上距离它
d的几个点。 -
选择到原点距离变化最大的那个邻居
-
重新调整估计的距离，这样他们的减法就会给你
d。然后使用重新缩放的第一个距离作为输出。
当放入片段代码时（使用 8 个邻居）：
// Fragment
#version 450 core
uniform dvec2 p0=vec2(0.0,+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,ar=0.0,nn=0.0;
for (x=0.0,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,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,1.0);
}
结果如下：
正如您所看到的，它现在更加正确（但由于上述非各向同性，非常接近边界是不准确的）。 8 个邻居在圆形斑点中产生 8 条对角线图案。如果你想摆脱它们，你应该围绕位置扫描整个圆圈，而不是仅仅 8 个点。
还有一些白点（它们不准确）我认为它们是当所选d遥远邻居跨越Mandelbrot边缘的情况而不是原件。可以过滤掉......（如果不是在不同的blob中，您知道同一方向的d/2距离应该是一半）
然而，即使是 8 个邻居也很慢。因此，为了更准确，我建议改用 2 pass“光线投射”方法。