novus log

(2010.5.03補充: 本篇程式碼有一些不理想的地方，已修正於第二版，這裡僅保留方便網友理解流程。更進階的 render 法見 Nebulabrot)

對上次 Julia Set 有興趣的人應該也會喜歡這個。

Buddhabrot 是 Mandelbrot Set 的一種特殊render法，叫這個名字的原因應該不用解釋，真佩服最初想出這些點子的人。

其實我不確定這樣寫對不對，好像參數改一點點就會跑出差蠻多的圖形。編譯很簡單，這已經是非常懶人包了，只要隨便弄個視窗的 HDC 餵 DrawBuddha 即可，懶得弄的話用 GetDC(0) 取得desktop DC也可以，再不然把 SetPixel 改成寫入圖片應該也不難。

逐點畫的方法很笨，改天再來用好一點的演算法來畫 Nebulabrot 好了

進一步的資訊可參考

http://www.superliminal.com/fractals/bbrot/bbrot.htm

http://www.linas.org/art-gallery/mandel/mandel.html

http://www.mrob.com/pub/muency/buddhabrot.html (這裡的pseudo code不完整，他沒有過濾掉會發散的點)

#include <windows.h>
#include <algorithm>

using namespace std;

template <typename T> struct Complex
{
    Complex(T rVal = 0, T iVal = 0) : real(rVal), imag(iVal) {}
    T real, imag;
};

const int ITER_MAX = 4096;
const int PLOT_WIDTH = 500, PLOT_HEIGHT = 500;

const double REAL_MAX = 1, REAL_MIN = -2;
const double IMAG_MAX = 1.5, IMAG_MIN = -1.5;

const double REAL_STEP = (REAL_MAX - REAL_MIN) / PLOT_WIDTH;
const double IMAG_STEP = (IMAG_MIN - IMAG_MAX) / PLOT_HEIGHT;

int trace[PLOT_WIDTH * PLOT_HEIGHT];

void PlotData(HDC hdc)
{
    int maxTrace = *max_element(trace, trace + PLOT_WIDTH * PLOT_HEIGHT);
    for (int y = 0; y < PLOT_HEIGHT; ++y) {
        int offset = y * PLOT_WIDTH;
        for (int x = 0; x < PLOT_WIDTH; ++x) {
            int d = (double)trace[x + offset] / maxTrace * 1024;
            if (d > 255) d = 255;
            SetPixel(hdc, x, y, RGB(d, d, d));
        }
    }
}

bool Calc(const Complex<double> &c, bool updateTrace)
{
    Complex<double> z = c;

    for (int i = 0; i < ITER_MAX; i++) {
        double rr = z.real * z.real;
        double ii = z.imag * z.imag;

        if (rr + ii >= 4) return false;

        if (updateTrace) {
            int px = (z.real - REAL_MIN) / REAL_STEP;
            int py = (IMAG_MAX - z.imag) / -IMAG_STEP;

            if (0 <= px && px < PLOT_WIDTH &&
                0 <= py && py < PLOT_HEIGHT) {
                ++trace[px + py * PLOT_WIDTH];
            }
        }

        z.imag = 2 * z.real * z.imag + c.imag;
        z.real = rr - ii + c.real;
    }
    return true;
}

void DrawBuddha(HDC hdc)
{
    Complex<double> c(REAL_MIN, IMAG_MAX);

    memset(trace, 0, sizeof(trace));

    for (int y = 0; y < PLOT_HEIGHT; ++y) {
        for (int x = 0; x < PLOT_WIDTH; ++x) {
            if (!Calc(c, false)) {
                Calc(c, true);
            } else {   // 下面這行可砍掉，不影響結果
                SetPixel(hdc, x, y, RGB(0, 0, 0));  
            }

            c.real += REAL_STEP;
        }
        c.imag += IMAG_STEP;
        c.real = REAL_MIN;
    }
    PlotData(hdc);
}