渲染曼德尔布罗集的最快算法是什么?

Viv*_*d V 4 c c++ graphics optimization mandelbrot

我尝试了许多用于渲染曼德尔布罗特集的算法,包括朴素的逃逸时间算法,以及优化的逃逸时间算法。但是,是否有更快的算法可以像我们在 YouTube 上看到的那样有效地产生真正深度的缩放。另外,我很想了解一些关于如何提高 C/C++ 之外的精度的想法 double

Spe*_*tre 5

与普通 GPU 相比,即使是高端 CPU 也会慢得多。即使使用 GPU 上的简单迭代算法,您也可以实现实时渲染。因此,在 GPU 上使用更好的算法可以获得高变焦,但是对于您需要的任何合适的算法:

  • 多通道渲染,因为我们无法在 GPU 上自行修改纹理
  • 高精度浮点作为floats/doubles还不够。

这里有一些相关的 QA:

这可能会让你开始......

加快速度的一种方法是您可以使用分数转义,就像我在第一个链接中所做的那样。它提高了图像质量,同时保持较低的最大迭代次数。

第二个链接将让您近似了解分形的哪些部分进入和退出以及距离有多远。它不是很准确,但可以用来避免对“肯定在外部”的部分进行计算迭代。

下一个链接将向您展示如何获得更好的精度。

最后一个链接是关于扰动的想法是,您仅对某些参考点使用高精度数学,并使用它以低精度数学计算其邻居点,而不会失去精度。然而从未使用过,但看起来很有希望。

最后,一旦您实现了快速渲染,您可能希望实现以下目标:

下面是 GLSL 中 3* 64 位双精度用于单值的小示例:

// high precision float (very slow)
dvec3 fnor(dvec3 a)
    {
    dvec3 c=a;
    if (abs(c.x)>1e-5){ c.y+=c.x; c.x=0.0; }
    if (abs(c.y)>1e+5){ c.z+=c.y; c.y=0.0; }
    return c;
    }
double fget(dvec3 a){ return a.x+a.y+a.z; }
dvec3 fset(double a){ return fnor(dvec3(a,0.0,0.0)); }
dvec3 fadd(dvec3 a,double b){ return fnor(a+fset(b)); }
dvec3 fsub(dvec3 a,double b){ return fnor(a-fset(b)); }
dvec3 fmul(dvec3 a,double b){ return fnor(a*b); }
dvec3 fadd(dvec3 a,dvec3 b){ return fnor(a+b); }
dvec3 fsub(dvec3 a,dvec3 b){ return fnor(a-b); }
dvec3 fmul(dvec3 a,dvec3 b)
    {
    dvec3 c;
    c =fnor(a*b.x);
    c+=fnor(a*b.y);
    c+=fnor(a*b.z);
    return fnor(c);
    }
Run Code Online (Sandbox Code Playgroud)

所以每个 hi 精度值是dvec3... fnor 中的阈值可以更改为任何范围。您可以将其转换vec3float...

[Edit1]“快速”C++ 示例

好吧,我想尝试我的新SSD1306 驱动程序和我的 AVR32 MCU 来计算 Mandelbrot,所以我可以将速度与这个Arduino + 3D + Pong + Mandelbrot进行比较。我使用 AT32UC3A3256,约 66MHz 无 FPU 无​​ GPU 和 128x64x1bpp 显示屏。无外部存储器,仅内部 16+32+32 KByte。Naive Mandlebrot 速度太慢(每帧大约 2.5 秒),所以我搞砸了这样的事情(利用这个位置和视图的缩放是连续的):

  1. 将分辨率降低 2

    为抖动腾出空间,因为我的输出只是黑白

  2. 使用变量最大迭代n基于zoom

    更改使n最后一帧无效以强制完全重新计算。我知道这很慢,但在变焦范围之间的转换时只发生 3 次。

    从最后一帧开始的缩放计数看起来不太好,因为它不是线性的。

    可以使用最后的计数,但为此还需要记住用于迭代的复杂变量,这会占用太多内存。

  3. 记住最后一帧以及哪个x,y屏幕坐标映射到哪个 Mandelbrot 坐标。

  4. 在每一帧上计算屏幕坐标和曼德尔布罗特坐标之间的映射。

  5. 重新映射最后一帧以调整到新位置和缩放

    因此,只需查看#3、#4中的数据,如果我们在最后一帧和实际帧中都有相同的位置(接近像素大小的一半),则复制像素。并重新计算其余部分。

    如果您的视图平滑(因此位置和缩放在每帧的基础上不会发生太大变化),这将极大地提高性能。

我知道它的描述有点模糊,所以这里有一个 C++ 代码,您可以在其中推断出所有疑问:

// high precision float (very slow)
dvec3 fnor(dvec3 a)
    {
    dvec3 c=a;
    if (abs(c.x)>1e-5){ c.y+=c.x; c.x=0.0; }
    if (abs(c.y)>1e+5){ c.z+=c.y; c.y=0.0; }
    return c;
    }
double fget(dvec3 a){ return a.x+a.y+a.z; }
dvec3 fset(double a){ return fnor(dvec3(a,0.0,0.0)); }
dvec3 fadd(dvec3 a,double b){ return fnor(a+fset(b)); }
dvec3 fsub(dvec3 a,double b){ return fnor(a-fset(b)); }
dvec3 fmul(dvec3 a,double b){ return fnor(a*b); }
dvec3 fadd(dvec3 a,dvec3 b){ return fnor(a+b); }
dvec3 fsub(dvec3 a,dvec3 b){ return fnor(a-b); }
dvec3 fmul(dvec3 a,dvec3 b)
    {
    dvec3 c;
    c =fnor(a*b.x);
    c+=fnor(a*b.y);
    c+=fnor(a*b.z);
    return fnor(c);
    }
Run Code Online (Sandbox Code Playgroud)

这些lcd内容shade8x8可以在链接的 SSD1306 QA 中找到。但是您可以忽略它,它只是抖动并输出像素,因此您可以i直接输出(即使没有缩放到<0..16>.

这里预览(在电脑上,因为我懒得连接相机......):

预览

因此 64x32 Mandelbrot 像素显示为 128x64 抖动图像。在我的 AVR32 上,这可能比简单方法快 8 倍(可能是 3-4fps)...代码可能会更优化,但请记住 Mandelbrot 并不是唯一运行的东西,因为我在后台有一些 ISR 处理程序来处理LCD 和我的TTS 引擎都基于此,从那时起我对其进行了很多升级,并用于调试它(是的,它可以与渲染并行发音)。另外,我的内存不足,因为我的 3D 引擎占用了大约 11 KB 的空间(主要是深度缓冲区)。

预览是用以下代码完成的(在计时器内):

static float zoom=1.0;
mandelbrot_draw<128,64,1>(+0.37,-0.1,zoom);
zoom*=1.02; if (zoom>100000) zoom=1.0;
Run Code Online (Sandbox Code Playgroud)

对于非 AVR32 C++ 环境也可以使用以下命令:

//------------------------------------------------------------------------------------------
#ifndef _AVR32_compiler_h
#define _AVR32_compiler_h
//------------------------------------------------------------------------------------------
typedef int32_t  S32;
typedef int16_t  S16;
typedef int8_t   S8;
typedef uint32_t U32;
typedef uint16_t U16;
typedef uint8_t  U8;
//------------------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------------------
Run Code Online (Sandbox Code Playgroud)

[Edit2] GLSL 中更高的浮点精度

Mandelbrot 的主要问题是它需要将指数差很大的数字相加。对于+,-运算,我们需要对齐两个操作数的尾数,将它们添加为整数并标准化回科学计数法。然而,如果指数差异很大,则结果尾数需要的位数超出了 32 位浮点数的容纳范围,因此仅保留 24 个最高有效位。这会产生舍入误差,导致像素化。如果你查看float二进制的 32 位,你会看到:

float a=1000.000,b=0.00001,c=a+b; 

//012345678901234567890123456789 ... just to easy count bits
a=1111101000b                                         // a=1000
b=         0.00000000000000001010011111000101101011b  // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b  // not rounded result
c=1111101000.00000000000000b                          // c=1000 rounded to 24 bits of mantissa
Run Code Online (Sandbox Code Playgroud)

现在的想法是扩大尾数位数。最简单的技巧是使用 2 个浮点数而不是 1 个:

//012345678901234567890123 ... just to easy count bits
a=1111101000b                                         // a=1000
                       //012345678901234567890123 ... just to easy count     b=         0.0000000000000000101001111100010110101100b    // b=0.00000999999974737875
c=1111101000.00000000000000001010011111000101101011b  // not rounded result
c=1111101000.00000000000000b                          // c=1000 rounded to 24 
 +          .0000000000000000101001111100010110101100b
                           //012345678901234567890123 ... just to easy count bits
Run Code Online (Sandbox Code Playgroud)

所以结果的一部分在一个浮点数中,其余部分在另一个浮点数中......每个值的浮点数越多,尾数就越大。然而,在 GLSL 中,将大尾数按位精确划分为 24 位块会很复杂且缓慢(即使由于 GLSL 限制而可能)。相反,我们可以为每个浮点数选择一定范围的指数(就像上面的示例一样)。

因此,在示例中,每个 (float) 值都有 3 个浮点数 (vec3)。每个坐标代表不同的范围:

       abs(x) <= 1e-5
1e-5 < abs(y) <= 1e+5
1e+5 < abs(z) 
Run Code Online (Sandbox Code Playgroud)

所以value = (x+y+z)我们有 3*24 位尾数,但范围并不完全匹配 24 位。为此,指数范围应除以:

log10(2^24)=7.2247198959355486851297334733878
Run Code Online (Sandbox Code Playgroud)

而不是 10...例如这样的:

       abs(x) <= 1e-7
1e-7 < abs(y) <= 1e+0
1e+0 < abs(z) 
Run Code Online (Sandbox Code Playgroud)

此外,还必须选择范围,以便它们处理您使用的值的范围,否则将毫无意义。因此,如果您的数字<4具有范围是没有意义的>10^+5,那么首先您需要查看您拥有的值的范围,然后将其分解为指数范围(每个值有多少个浮点数)。

请注意,仍然会发生一些舍入(但比本机浮点数少得多)!

现在对这些数字进行操作比对普通浮点数进行操作稍微复杂一些,因为您需要将每个值处理为所有组件的括号总和,因此:

(x0+y0+z0) + (x1+y1+z1) = (x0+x1 + y0+y1 + z0+z1)
(x0+y0+z0) - (x1+y1+z1) = (x0-x1 + y0-y1 + z0-z1)
(x0+y0+z0) * (x1+y1+z1) = x0*(x1+y1+z1) + y0*(x1+y1+z1) + z0*(x1+y1+z1)
Run Code Online (Sandbox Code Playgroud)

并且不要忘记将值标准化回定义的范围。避免添加小和大(abs)值,因此避免x0+z0等等......

[Edit3] 新的 win32 演示 CPU 与 GPU

两个可执行文件都预设在同一位置,并缩放以显示doubles 何时开始四舍五入。我必须稍微升级一下坐标的计算方式,px,py因为 y 轴在这个位置开始偏离 10^9 左右(对于其他位置来说阈值仍然可能太大)

此处预览 CPU 与 GPU 的高变焦 (n=1600):

CPU 与 GPU

CPU 的 RT GIF 捕获(n=62++,GIF 缩小 4 倍):

CPU放大