c++ - 优化 Mandelbrot 分形

Question

这是一个在 .ppm 文件中输出 mandelbrot 分形的代码。我该如何优化它？

#include<bits/stdc++.h>
using namespace std;

int findMandelbrot(double cr, double ci, int max_iterations)
{
    int i = 0;
    double zr = 0.0, zi = 0.0;
    while (i < max_iterations && zr * zr + zi * zi < 4.0)
    {
        double temp = zr * zr - zi * zi + cr;
        zi = 2.0 * zr * zi + ci;
        zr = temp;
        ++i;
    }
    return i;
}

double mapToReal(int x, int imageWidth, double minR, double maxR)
{
    double range = maxR - minR;
    return x * (range / imageWidth) + minR;
}

double mapToImaginary(int y, int imageHeight, double minI, double maxI)
{
    double range = maxI - minI;
    return y * (range / imageHeight) + minI;
}

int main()
{
    ifstream f("input.txt");
    int imageWidth, imageHeight, maxN;
    double minR, maxR, minI, maxI;

    if (!f)
    {
        cout << "Could not open file!" << endl;
        return 1;
    }

    f >> imageWidth >> imageHeight >> maxN;
    f >> minR >> maxR >> minI >> maxI;

    ofstream g("output_image.ppm");
    g << "P3" << endl;
    g << imageWidth << " " << imageHeight << endl;
    g << "255" << endl;


    double start = clock();

    for (int i = 0; i < imageHeight; i++)
    {
        for (int j = 0; j < imageWidth; j++)
        {
            double cr = mapToReal(j, imageWidth, minR, maxR);
            double ci = mapToImaginary(i, imageHeight, minI, maxI);

            int n = findMandelbrot(cr, ci, maxN);

            int r = ((int)sqrt(n) % 256);
            int gr = (2*n % 256);
            int b = (n % 256);

            g << r << " " << gr << " " << b << " ";
        }
        g << endl;

        if(i == imageHeight / 2) break;
    }

    cout << "Finished!" << endl;

    double stop = clock();

    cout << (stop-start)/CLOCKS_PER_SEC;
    return 0;
}

我一直走到 imageHeight/2 因为在 photoshop 中我可以只复制另一半。我在考虑对数幂，但尝试了一些东西并且只适用于整数......

Answer 1

有很多方法可以优化 Mandelbrot 分形。

一种方法是针对您的 CPU 甚至 GPU 优化代码。令人印象深刻的加速显示在:Mandelbrot with SSE, AVX and OpenCL .这将内部循环优化了将近 1000 倍。速度提高了 3 个数量级。

但是还有其他的优化方式。您已经提到的第一个:mandelbrot 集在 y=0 处镜像。所以你只需要一半。有一些更简单的方法可以避免运行内部循环。如果您向下滚动 Wikipedia on Mandelbrot到优化，您会看到“心形/灯泡检查”。这是对苹果形主要部分或直接在其左侧的圆中的点的简单检查。对于可以涵盖很多点的海报。

我见过的另一种加速方法是在生成预览时使用距离估计(也在维基百科上)方法或仅使用黑白设置的 mandelbrot 轮廓。计算图像中的随机点，如果在 mandelbrot 集之外，则绘制一个圆，其半径由距离估计方法给出。该圆圈中的任何东西都不是 mandelbrot 集的一部分。这可以快速覆盖 mandelbrot 集之外的许多点。

还有一些近似结果的方法可能不会产生完美的图像，但通常已经足够好了。例如:

• 边界追踪

计算边界后的点，其中像素需要 N 和 N+1 次迭代才能逃逸。 N+1 和 N+2 相同。这两个边界之间的所有内容都需要 N 次迭代。

• 分而治之的盒子

计算矩形的边界(从整个图像开始)。如果边界全部进行 N 次迭代，则矩形内部进行 N 次迭代并被填充。否则将矩形分成 4 份并重复每份。通过计算 2 或 3 像素宽的边框可能会改善结果，但保存的会更少。

• 猜测

以低分辨率计算图像，然后将分辨率加倍以保持计算点。检查图像，如果 5x5 个原始点具有相同的 N，则在内部 3x3 点周围填充一个矩形(也适用于 3x3、7x7、9x9 等点)。未填写的分数由您计算。然后重复整个过程，直到你得到最终的解决方案。

• 单轨道迭代

这是最难做到的，我只见过它的一种实现方式。这个想法是靠近在一起的点在迭代中表现相同。如果您迭代一个 3x3 点的网格(首先覆盖整个图像)，您可以使用牛顿插值法对中间点的值进行插值。这工作得很好——直到它不工作。

因此，除了 3x3 点网格之外，您还迭代了 4 个测试点，即每个网格单元的中间。对这 13 个点进行 8 次迭代，然后从网格点中插入 4 个测试点。如果计算结果和插值结果相差太大(这是困难的部分)，您将丢弃最后 8 次迭代并将网格 segmentation 为 4 个半尺寸网格。您插入的缺失点。重复直到达到最终分辨率。

即使它只适用于几次迭代，潜在的好处也是巨大的。假设您想要 40000 x 40000 像素的图像。如果 SOI 在第一次 segmentation 之前工作 10 个循环(80 次迭代)，那么您可以通过计算 1040 和一些插值和检查来节省 80 * 40000 * 40000 = 128_000_000_000 次迭代。或者加速 123_076_923，8 个数量级。但仅限于前 80 次迭代。随着网格被越来越多地划分，加速变得越来越小。但是，节省的每一点都会累加起来。

此方法的优点在于适用于平滑/连续着色或映射到高度。其他方法仅适用于将迭代映射到色带。