优化曼德博集合分形

Question

优化曼德博集合分形

5

这是一个输出曼德博集合分形图案的代码，保存为 .ppm 文件。如何对其进行优化？

#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;
}

我只翻译图像高度的一半，因为在Photoshop中我可以复制另一半。我曾考虑过对数幂，但尝试后发现只适用于整数...

- user7312333

2

如果您已经有可用的代码并且想要改进它，那么这似乎是一个适合 https://codereview.stackexchange.com 的问题。 - Jesper Juhl

是的，这是输入内容：“512 512 512 -1.5 0.7 -1.0 1.0”。我想要将其与数据一起使用，例如：“40000 40000 40000 -1.5 0.7 -1.0 1.0”，甚至更多，这时时间部分就变得棘手了。 - user7312333

输入：512 512 512 -1.5 0.7 -1.0 1.0处理时间为0.315秒。 - user7312333

你使用了哪个CPU和哪些标志位？ - Anty

@Anty 我的处理器是 i7 4710 HQ 2.50GHz。 - user7312333

显示剩余3条评论

3个回答

1

所以这就是热循环：

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;

我知道如何在快速CPU指令中实现非整数幂，但这并不能解决复数问题。使用std::complex也没有帮助。你不会因为非内联而获得任何优势，也无法随时应用优化策略。所以我能做的最好的方法是：

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

是的，我知道将一个整数测试从循环中取出并没有带来太多的好处。我只找到了另一个优化器，您需要检查它是否真的更快：

int i = max_iterations;
double zr = 0.0, zi = 0.0;
do {
    double tempr = zr * zr - zi * zi + cr;
    double tempi = zr * zi;
    zi = tempi + tempi + ci;
    zr = tempr;
} while (--i && zr * zr + zi * zi < 4.0);
return max_iterations - i;

那就是全部内容。

- Joshua

你使用了工具来生成这些优化吗？似乎编译器应该能够自动执行第二个优化。 - Matt Messersmith

@mwm314，我使用gnuc++作为我的编译器，但它不具备自动优化功能。 - user7312333

@mwm314：我还没有见过一个编译器聪明到能够找到第一个问题。我见过几个编译器，如果在循环后丢弃 i，它们可以工作，但如果使用 i，则无法工作。至于第二个问题，不行。我使用领域知识知道它是安全的（由于双精度浮点数中累积舍入的方式，通常情况下是不安全的）。 - Joshua

在第一个循环中，zr和zi都是0.0。因此，在循环后，zr = cr和zi = ci。最好避免整个循环，而不仅仅是测试。 - Goswin von Brederlow

1

在findMandelbrot中，您在循环测试中使用表达式zr * zr和zi * zi，但如果测试成功，则重新计算相同的两个表达式。因此，一个明显的改进可能是使用类似以下内容的缓存...

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

- G.M.

我不知道为什么，但是使用你的实现方式得到了0.46，而使用我的实现方式得到了0.38。 - user7312333

编译器可能足够聪明，可以“缓存”这些临时变量。 - Matt Messersmith

@mwm314 确实。就我所知，对于最初指定的输入参数，我确实看到了5-10%的持续改进。但我同意这可能高度依赖于编译器（我在Linux上使用g++ 7.1.1）。如果编译器正在缓存这些值，那么编写明确缓存它们的代码可能会使事情变得更加混乱。 - G.M.

我没有看到加速，我的机器上情况更糟。不过我使用的是g++ 4.8.4。 - Matt Messersmith

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Goswin von Brederlow · Accepted Answer

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

其中一种方法是针对你的CPU甚至GPU进行代码优化。在SSE、AVX和OpenCL下的Mandelbrot中，展示了令人印象深刻的加速效果，将内部循环优化了近1000倍，快了3个数量级。

但是还有其他优化方法。你已经提到的第一个：Mandelbrot集在y=0处镜像。因此，只需要一半即可。还有一些更简单的方法来避免运行内部循环。如果你在维基百科上浏览Mandelbrot页面并向下滚动到Optimizations，你会看到“心形/球茎检查”。这是一个简单的检查，用于检查主要部分的苹果形状或直接位于其左侧的圆中的点。对于涵盖许多点的海报而言，这是个好办法。

我见过另一种加速方法，这种方法在生成预览或仅以黑白形式轮廓呈现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次迭代，并将网格分成4个半大小的网格。您需要进行插值来填充丢失的点。重复此过程直到达到最终分辨率。

即使它只能工作几次迭代，潜在收益也非常大。假设您想要40000×40000像素的图像。如果SOI在第一次细分之前工作10次循环（80次迭代），那么通过计算1040个点并进行一些插值和检查，您就可以节省80*40000*40000 = 128_000_000_000 次迭代。或者说加速了123_076_923倍，是8个数量级的加速。但仅限于前80次迭代。随着网格被分割越来越多，加速效果会逐渐减少。但每个节省下来的计算量都是有用的。这种方法的好处在于可以平滑地进行着色或将其映射到高度上。其他方法只能将迭代映射到颜色带上。