Mandelbrot set 以parallel_for_實現

凤凰_1發表於2024-06-19
Parallel

我們將以繪製曼德布羅集合為例，展示如何從常規的順序程式碼輕鬆地修改程式碼以實現並行化計算。

曼德布羅特集理論：

曼德布羅特集的定義是以數學家本諾·曼德布羅特的名字命名的，由阿德里安·杜瓦迪命名。它因其在數學領域之外的形象表示而聞名，因為它是一個類分形的例子，一個在每個尺度上重複顯示模式的數學集合（更進一步說，曼德布羅特集是自相似的，因為整個形狀可以在不同的尺度上反覆看到）。對於更深入的介紹，您可以檢視相應的維基百科文章。在這裡，我們將只介紹繪製曼德布羅特集的公式（來自上述維基百科文章）。

Mandelbrot集合是複平面上 的值集，對於這些值，在二次對映的迭代下，從0開始的軌道保持收斂。 二次迭代對映表示式如下：

這就是說，複數c是Mandelbrot集合的一部分，如果從z0=0開始並反覆應用迭代，無論n有多大，zn的絕對值仍然保持有界。這也可以表示為

虛擬碼

🔍發現奇蹟，每一畫素藏著無限秘密！在五彩光闌的影像宇宙裡，我們一起用最簡單的逃逸時間演算法，探索複數世界的無盡邊界。🌈Mandelbrot集的每一點色彩，都是遞迴關係的精彩演繹。我們一起來看看你的世界裡，哪些畫素會逃逸，哪些會在迭代中沉澱。🌟 加入我們的旅程，讓數字藝術綻放在你的螢幕上。🎨 不斷迭代，細節漸顯，時間彷彿也在這一花一世界中流轉。🎯 專注每個細節，渲染出屬於你的Mandelbrot世界地圖。🚀 跟著我，一起繪製未知，記錄每一次迭代的成長軌跡。🎓

For each pixel (Px, Py) on the screen, do:
{
    x0 = scaled x coordinate of pixel (scaled to lie in the Mandelbrot X scale (-2, 1));
    y0 = scaled y coordinate of pixel (scaled to lie in the Mandelbrot Y scale (-1, 1));
    x = 0.0;
    y = 0.0
    iteration = 0
    max_iteration = 1000
    while (x*x + y*y < 2*2 AND iteration < max_iteration) {
        xtemp = x*x - y*y + x0
        y = 2*x*y + y0
        x = xtemp
        iteration = iteration + 1
    }
    color = palette[iteration]
    plot(Px, Py, color)
}

為了把理論與虛擬碼聯絡起來，我們將曼德布羅特表示式寫入如下表示式：

在這個圖中，我們回憶一下實數部分位於x軸上，虛數部分位於y軸上。您可以看到，如果我們放大特定位置，整個形狀可以重複出現。

實現程式碼如下：

//1.Escape time algorithm implementation
int mandelbrot(const complex<float> &z0, const int max)
{
    complex<float> z = z0;
    for (int t = 0; t < max; t++)
    {
        if (z.real()*z.real() + z.imag()*z.imag() > 4.0f) return t;
        z = z*z + z0;
    }
    return max;
}
//2.Sequential Mandelbrot implementation
int mandelbrotFormula(const complex<float> &z0, const int maxIter=500)
{
    int value = mandelbrot(z0, maxIter);
    if(maxIter - value == 0)
    {
        return 0;
    }

    return cvRound(sqrt(value / (float) maxIter) * 255);
}
//3.Parallel Mandelbrot implementation
class ParallelMandelbrot : public ParallelLoopBody
{
public:
    ParallelMandelbrot (Mat &img, const float x1, const float y1, const float scaleX, const float scaleY)
        : m_img(img), m_x1(x1), m_y1(y1), m_scaleX(scaleX), m_scaleY(scaleY)
    {  }
    virtual void operator ()(const Range& range) const CV_OVERRIDE
    {
        for (int r = range.start; r < range.end; r++)
        {
            int i = r / m_img.cols;
            int j = r % m_img.cols;

            float x0 = j / m_scaleX + m_x1;
            float y0 = i / m_scaleY + m_y1;

            complex<float> z0(x0, y0);
            uchar value = (uchar) mandelbrotFormula(z0);
            m_img.ptr<uchar>(i)[j] = value;
        }
    }

    ParallelMandelbrot& operator=(const ParallelMandelbrot &) {
        return *this;
    };

private:
 Mat &m_img;
 float m_x1;
 float m_y1;
 float m_scaleX;
 float m_scaleY;
};
//4.示例
int main()
{
    Mat mandelbrotImg(4800, 5400, CV_8U);
    float x1 = -2.1f, x2 = 0.6f;
    float y1 = -1.2f, y2 = 1.2f;
    float scaleX = mandelbrotImg.cols / (x2 - x1);
    float scaleY = mandelbrotImg.rows / (y2 - y1);

    cv::parallel_for_(Range(0, mandelbrotImg.rows*mandelbrotImg.cols), [&](const Range& range)
    {
        for (int r = range.start; r < range.end; r++)
        {
            int i = r / mandelbrotImg.cols;
            int j = r % mandelbrotImg.cols;

            float x0 = j / scaleX + x1;
            float y0 = i / scaleY + y1;

            complex<float> z0(x0, y0);
            uchar value = (uchar) mandelbrotFormula(z0);
            mandelbrotImg.ptr<uchar>(i)[j] = value;
        }
    });
    string winname="mandelbrotImg";
    namedWindow(winname,WINDOW_FREERATIO);
    imshow(winname,mandelbrotImg);
    waitKey();
    return 0;
}

相關文章