Mandelbrot集的内部距离估计算法 [英] Interior Distance Estimate algorithm for the Mandelbrot set
问题描述
我正在寻找一种方法来估计从Mandelbrot集内部的点到Mandelbrot集的边界的距离，以用于GLSL着色器.
至少可以保证下限的东西很多.
我敢打赌， 1./length(z)达到了 float 的精度，如果有任何区别，请使用 double 和 dvec2 而不是 float，vec2 .如果可以，那么我会忽略 length(z)的太小的值.
或者，您可以
只需查找//***这就是我在代码中添加的注释，它是添加到标准mandelbrot渲染中以渲染距离的内容.ps我的(x，y)是您的 z 和(cx，cy)是您的 c
无论如何，距离仍然是高度非线性的，并取决于位置
[Edit1]非各向同性尺度
黑点是您可以将阈值大小设置为 1e-20 的阈值...现在，我添加了水平线以显示距离标度的分布(因为我不知道非等向性和非线性...)这是输出:
并对片段的一部分进行着色(在 for 循环之后):
ar = 1.0-(2.0 * nn/ar);aa = 10.0 * ar;//每单位10条水平线aa- = floor(aa);如果(abs(aa)< 0.05)col = vec4(0.0,1.0,0.0,1.0);//水平线的宽度和颜色else col = vec4(ar，ar，ar，1.0); 如您所见，它与边界不是很平行，但是在局部仍是恒定"的.(在分形的局部特征中，水平线与每个等距线等距)，因此，如果使用梯度(微分)，则结果将只是非常粗略的估计(但应该起作用).如果足够，您应该做的是:
-
计算查询位置的非线性距离，并在全部"路径中计算出距离它很少的点
d.方向. -
选择距离原点距离变化最大的那个邻居
-
重新缩放估计的距离，这样它们的减法将为您提供
d.然后使用重新缩放的第一个距离作为输出.
放入片段代码(使用8个邻居)时:
<代码>//片段#version 450核心均匀dvec2 p0 = vec2(0.0,0.0);//鼠标位置< -1，+ 1>均匀双倍缩放= 1.000；//飞涨 [-]均匀整数n = 100;//迭代[-]在平滑vec2 p32中；vec4 col;double mandelbrot_distance(double cx，double cy){//到边界的距离int i，j;双x，y，q，xx，yy，ar = 0.0，aa，nn = 0.0;对于(x = 0.0，y = 0.0，xx = 0.0，yy = 0.0，i = 0;(i 0.0){d0 = mandelbrot_distance(cx-d，cy);如果(d0> 0.0)d0 = abs(d0-dd);d1 = mandelbrot_distance(cx + d，cy);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx，cy-d);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx，cy + d);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx-e，cy-e);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx + e，cy-e);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx-e，cy + e);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}d1 = mandelbrot_distance(cx + e，cy + e);如果(d1＞ 0.0){d1 = abs(d1-dd);如果(d0< d1)d0 = d1;}dd * = d/d0;}dd * = zoom;//仅用于小细节的可视化，实际距离不应由此缩放col = vec4(dd，dd，dd，1.0);} 结果在这里:
如您所见，它现在更正确了(但由于上述非等向性，非常接近边界是不准确的).8个邻居在圆形斑点中产生8条对角线样的线型.如果要摆脱它们，则应扫描该位置周围的整个圆圈，而不是仅扫描8个点.
还存在一些白点(它们与精度无关)，我认为是某些情况，即选定的 d 远方邻居跨越了mandelbrot边缘的斑点与原始斑点不同.可以将其过滤掉...(您知道在同一方向上的 d/2 距离应该是一半，如果不是，则位于不同的斑点)
但是，即使是8个邻居也很慢.因此，为了获得更高的准确性，我建议您进行2次射线投射"测试.方法.
I'm looking for a way to estimate the distance to the boundary of the Mandelbrot set from a point inside of it for use in a GLSL shader.
This page links to various resources online touching on the subject of interior distance estimation such as the underlying mathematical formula, a Haskell implementation, some other blogs, forum posts and a C99 implementation, but I got the impression that they are all either very complex to implement or very computationally heavy to run.
After many hours of trying, I managed to make this code that runs in 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);
}
It does produce a gradient in every bulb, but it is clear that it's not usable as a distance estimator by itself because values in smaller bulbs have inconsistent magnitude (brightness) compared to bigger bulbs. So it's clear that a parameter is missing but I don't know what it is.
I don't require a perfect solution nor one that converges into a perfect solution like in this image.
Something that at least guarantees a lower bound is plenty.
My bet is that 1./length(z) is hitting the precision of float try to use double and dvec2 instead of float,vec2 if it makes any difference. If it does then I would ignore too small values of length(z).
Alternatively you can render just the boundary into texture in one pass and then just scan neighbors in all directions until boundary found returning the ray length. (may require some morphology operators before safe use)
This can be speed up with another pass where you "flood" fill incrementing distance into texture until its filled (better done on CPU side as you need R/W access to the same texture) its similar to A* filling however your precision will be limited by texture resolution.
If I ported my mandlebrot from link above to your computation ported to doubles and added the threshold:
// 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
}
I got these outputs:
Just look for // *** this is what I added comment in the code that is what is added to the standard mandelbrot rendering to render distance instead. ps my (x,y) is your z and (cx,cy) is your c
anyway the distance is still highly nonlinear and depends on the position
[Edit1] non-isotropic scale
The black dot is the thresholds size you can lover it to 1e-20 ... Now I added level lines to show the distribution of distance scale (as I did not know how non isotropic and non linear it is...) here the output:
And coloring part of fragment (after the for loop):
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);
As you can see its not very parallel to border but still locally "constant" (the level lines are equidistant to each in local feature of fractal) so if gradient (derivate) used the result will be just very rough estimate (but should work). If that is enough what you should do is:
compute non linear distance for queried position and few point
ddistant to it in "all" directions.pick that neighbor that has bigest change in distance to original point
rescale the estimated distances so their substraction will give you
d. then use the first distance rescaled as output.
When put to fragment code (using 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);
}
here the result:
As you can see its now much more correct (but very close to border is inaccurate due to non isotropy mentioned above). The 8 neighbors produces the 8 diagonal like lines pattern in the circular blobs. If you want to get rid of them you should scan whole circle around the position instead of just 8 points.
Also there are still some white dots (they are not accuracy related) I think they are cases when the selected d distant neighbor is across the mandelbrot edge in different blob than original. That could be filtered out ... (you know d/2 distance in the same direction should be half if not you are in different blob)
However even 8 neighbors are pretty slow. So I for more accuracy I would recommend to go for the 2 pass "ray casting" method instead.
这篇关于Mandelbrot集的内部距离估计算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
