使用OpenGL进行Mandelbrot集的可视化

Mandelbrot集是哪一集？？

Mandelbrot集不是哪一集！！ 啊不对……
Mandelbrot集是哪一集！！ 好像也不对……
Mandelbrot集是数集！！ 因此……他不是一集而是数集？？……

因此这个M...dem...集究竟是什么啊？？

Mandelbrot集是一个数集

Mandelbrot集\(\mathbb{M}\)（简称曼集）是一个由二元复数构成的集合，也就是一个复数集：
\[ \mathbb{M}\subset\mathbb{C} \]

也就是说，曼集的元素都是复数，也就是以下的形式：
\[ \forall m\in \mathbb{M},m=a+b\mathrm{i},\left(a,b\in \mathbb{R}\right) \]

为何叫这么复杂的名字？？

本华·曼德勃罗特（Benoit B.Mandelbrot）是分形学（Fractals） 的鼻祖，是的，Fractals这个词就是他提出的。
不过，这个集合自己并非他本人提出的，而是一个叫作阿德里安·杜阿迪的法国数学家提出的，至于为何以曼德勃罗特命名，根据wiki的说法是提出者为了向曼德勃罗特致敬而放弃用本身的名字命名。

这是个什么样子的集合？？

一个复数\(z\)在或不在曼集中，取决于以下数列是否收敛：
\[ \zeta:z,f(z;c),f(f(z;c);c),\cdots,f\circ^n(z;c) \tag{1} \]
\[ z,c\in\mathbb{C} \]
其中，\(f(z;c)=z^2+c\)，\(c\)是一个复参数。
若是把\(f\circ^n(z;c)\)展开写的话就是：
\[ \begin{split} \zeta_n&=f\circ^n(z;c)=z_n^2+nc \\ z_n&=z_{n-1}^2+c \end{split} \tag{2} \]
因为在其中迭代的前一个\(z\)都是以平方的形式出现，这也就意味着
\[ \left|\Re\left(z_n\right)\right|\ge\left|\Re\left(z_{n-1}\right)\right|\ge 0 \\ \left|\Im\left(z_n\right)\right|\ge\left|\Im\left(z_{n-1}\right)\right|\ge 0 \]
这样一来，若\(z_n\)收敛，则\(z_m(m<n)\)必然收敛。

综上描述，有下结论：
\[ z\in\mathbb{M} \Leftrightarrow \lim_{n\rightarrow\infty}\zeta_n\in\mathbb{C} \tag{3} \]

式3即曼集断定的充要条件或定义。

为了让问题变得简单，这里能够令初始值\(z=0\)，这样式2变成了仅与参数\(c\)相关的形式：
\[ \begin{split} z_0 &=0 \\ z_1 &=c \\ z_2 &= z_1^2+c = c^2+c \\ z_3 &= z_2^2+c = c^4+2c^3+c^2+c \\ &\cdots \\ \end{split} \]

复数？？迭代？？还收敛？？听起来好像很复杂……

这个……是的，这些东西听起来很不直观，以致于我也花了很久看了大量资料才搞清楚是什么意思。

可是，这件事自己的中心思想并不复杂，其无非作了这样的一件事情：

1. 给定一个复数\(c\)
2. 计算\(z_+=z^2+c\)（\(z_+\)会成为下一次的\(z\)）
3. 就这样，反复执行2，直到山穷水尽之时
   1. \(z_+\)仍然没有飞升为\(\infty\)，那\(c\in\mathbb{M}\)
   2. 不然存在某一次\(z_+=\infty\)则\(c\notin\mathbb{M}\)

山穷水尽之时？？那我恐怕等不到那一天了……

没必要沮丧，咱们甚至没打算去推导出\(z_n\)的完整展开形式，并且，计算机也没法处理须要运行无数次才能解决的问题（停机问题）。

可是，在有限的，可预见的次数（迭代数）内，进行这样的断定仍是可行的而且能得出结果的。

可是这样的话确定会有纰漏吧

说的对，由于只要迭代数\(n\)有限，\(z_n\)是总能求出肯定值。但正由于如此，迭代数的增加才会产生实际的意义。这就好像，没有人能求出\(\pi\)的全部小数位数，可是每当有人多求出正确的一位的时候，其结果与实际值就总会更加接近。也就是说，只要范围有限，咱们求得的曼集就总比实际的曼集大，但随着次数的增多，结果集合会向曼集逐渐地靠近。

所幸，在有限的次数内，咱们尽管不能当即断言哪些数必定属于曼集，但咱们能断言哪些数必定不属于曼集，那正是：

还没上高速公路就要起飞的，那必然是要玩完的

在有限次计算就体现出明显的发散趋势的时候，那基本也就抢救无效了

好吧好吧，那有哪些值是能够放弃治疗的？

根据迭代函数的定义，\(f(z;c)=z^2+c\)，咱们对函数的范数（也就是模，或者称为复数的绝对值）进行考察，至少若是\(f\)收敛，其模与带入z的模至少是相等的（固然，原则上幅角也应当相等），：
\[ |f(z;c)|=\left|z^2+c\right|=|z| \\ |z|\ge0 \]

根据绝对值三角不等式（这个对复数域一样成立，实际上就是向量绝对值三角不等式的移植），若\(|z|=0\)则\(|c|\)必然为\(0\)，不然
\[ |z|^2-|c|=|z^2|-|c|\le|f(z;c)|=|z|\le|z^2|+|c|=|z|^2+|c| \\ 1-|c|\le|z|\le1+|c| \]
因为迭代初始条件\(z_0=0\)，所以通过n次迭代后的\(z\)必定是
\[ z_n=\sum_{m=0}^{(n-2)^2}{k_mc^{m}} \\ c\le|z_n|\le\sum_{m=0}^{(n-2)^2}{|k_m||c^{m}|} \]
上面的不等数右侧产生了一个很熟悉的形式——幂级数，若是要幂级数\(\sum_{m=0}^{(n-2)^2}{|k_m||c^{m}|}\)收敛，则\(|c|<1\)是必要条件（因为\(k_m\)的大小这里没有多加考虑，即便\(|c|<1\)知足，发散的可能一样存在，这个在后续的迭代过程当中会进一步的被发掘出来），不然根据幂级数的性质级数必然发散，这样获得了一个必要阈值：\(|c_\Theta|=1\)，将这个阈值带入，咱们就能获得\(|z|\)的基础发散断定值：
\[ |z|\ge2\Rightarrow \lim_{n\rightarrow\infty}{\zeta_n}=\infty \]

在后面的程序中，咱们将使用这个值做为发散断定，意思就是说，若是算出\(|z|\ge 2\)了，那就不用考虑了，你不可能收敛的，至于剩下的值，我只假设但不断言它必定在集合内，所以我暂且将它认为是集合内的，若是通过更多的计算发现也达到这个红线条件了，那么一样不论这个值挣扎了多久，OUT！

可是，复数的收敛并不限于模收敛，当复数被整理为辐角形式\(r\angle\theta\)时，若\(\theta\)没有稳定值即使\(r\)已经收敛仍然不能断言复数收敛（毕竟啊，那个复平面上的向量若是一直在原地打转，份量就总在以一个有界值变化，宛如论证的丧钟在转动），可是对于这个函数而言根本不会不存在这种情形！

哇，断言的这么轻松？？何出此言？？

其实很容易，找到这样的向量，利用迭代函数平方相加后保持模不变便可，所幸的是，这样的值很是有限，由于首先一个条件就将范围急剧缩减：
\[ |c|=1 \]

是的，只有\(|c|=1\)的时候，平方计算后才可能保证模不变（依据辐角计算法若\(z=r_z\angle\theta\)，则\(z^2=r_z^2\angle{2\theta}\)），固然我是指仅一次，由于还另外须要一件事保证\(z^2+c=c^2+c\)后的模长也是\(|c|\)或者说是\(1\)：
\[ \cos{<c^2,c>}=-\frac{1}{2} \]
也就是说\(\arg{z^2}-\arg{c}=\frac{2k\pi}{3}\)才能够（这个式子写的不是很严谨，但你当成是两个夹角为\(120^\circ\)就好），固然你会反驳：

\(c\)是随意取得，也就是说幅角就能够是任意值，固然，\(z^2\)的取值决定于\(c\)，也就是\(z^2\)能够和\(c\)构造连续的函数映射，那么你根本保证不了\(z^2\)必定可以避开这些值啊你个代数白痴！！

啊是的，可是你很快就会意识到一个事实：

这个丧钟，转不久的！！

什么？？

确实，当知足上述两个条件以后，获得的\(z_+\)的结果模仍然是\(1\)，只是转了个角度，可是架不住这个函数是迭代的啊，此次获得的\(z_+\)很明显位于这\(z^2\)和\(c\)的角平分线处……

而后，夹角条件被打破了！！

就像雪崩同样，这个微妙的平衡很快就被打破了，而且立刻开始急剧的发散（或收敛）。

如丸走坂
- Ramp Rollerball

这样一来，经过屡次的迭代运算，不合格的点一点一点被筛出，若是次数足够多，这个集合最后造成的图形被称为Mandelbrot分形。
大概长出来是这样的：

这个M..d什么b什么的分形有什么特色么？？

做为一个分形，它最大的特色之一就是：

德罗斯特效应（自类似性，Droste Effect）

若是对这个图形局部放大，你会看到这个与最开始看到一致的子图性，并且：

放大一次有，一直放大一直有！放大放大再放大！每个子部都看的清清楚楚！

德罗斯特效应本质是递归式的图形体现，由于迭代函数经过\(z_+=f(z;c)\)进行了递推，这样的话，也就意味着，分形具备无限细微的结构，并且跟总体是类似的。

而具备无限精细的结构也就意味着：它具备有限的面积，而周长却极可能无穷大

固然因为曼集的边缘性质很是复杂，这里就再也不详细讨论他的周长是否收敛了（面积收敛是一看就能看出来的），毕竟这篇文章是为了说明如何经过OpenGL实现这个集合的可视化，上面其实已经花了大量的篇幅去推导一个阈值了。

Mandelbrot集的兄弟——Julia集简介

以前是以\(c\)做为参数进行迭代，获得了曼集，那么，若是把\(c\)和\(z\)交换一下，让\(z\)作参数，就获得了Mandelbrot集的兄弟——Julia集（如下简称朱集），以法国数学家加斯顿·朱利亚（Gaston Julia）命名。

— Oh,Julia~
— Oh,Mandelbrot~
— W...What??

一样，由于递归仍然存在，Julia集获得的也会是Julia分形：

迭代函数没变，所以它的发散断定也是和曼集是一致的，这里仅仅介绍一下这个亲戚，并不打算实现这个Julia集的可视化。

OpenGL可视化实现

其实上面的函数都已经推导完毕的状况下，而且也知道集合是如何生成的以后，代码实现就容易多了，绘制曼集的代码以下（只考虑整数点）：

基本绘制（二值化）

void DrawMandelbrot(int maxIter)
{
    complex<double> z(0,0); // 位于标准库<complex>头文件中
    complex<double> c;
    complex<double> f(0,0);
    bool diverge = false;
    glBegin(GL_POINTS);
    for (GLint x = -400; x < 400; x++)
    {
        diverge = false;
        for (GLint y = -300; y < 300; y++)
        {
            c = { static_cast<double>(x)/150,static_cast<double>(y)/150 };
            int iter;
            for (iter = 0; iter < maxIter; iter++)
            {
                z = z*z + c;
                if ( abs(z) >= 2 ) //发散断定
                {
                    diverge = true;
                    break;
                }
            }
            if (!diverge)
            {
                glColor3f(0.0f, 0.0f, 0.0f); //集合内黑色
            }
            else
            {
                glColor3f(1.0f, 1.0f, 1.0f); //集合外白色
                diverge = false;
            }
            z = { 0,0 };
            glVertex2i(x, y);
        }
    }
    glEnd();
}

效果图：

有问题，若是发散断定值不取2会有什么影响？？

问得好，实际上我正想说，其实断定值最好取2，若是实在没法取得的话，也至少请保证它大于2
若是小于2，则有些原本应该在曼集内的点被剔除，这样造成的曼集图形是不完整的。

而若是大于2，当迭代次数很是大的时候，实际上获得的结果与2的时候相差不大（细微的内部可能会有所差异），由于次数很是大的时候，该发散的值通常早都散得没影了，用再大的有限值根本拦不住。可是若是迭代次数很是小的时候，获得的曼集会包含大量冗余的点，由于阈值没能称职地起到约束做用。

迭代着色

固然，不管如何，咱们可让上面的图更加漂亮一些，好比咱们以迭代多少次就发散了做为着色的依据，这里对上面的图形进行着色处理。
实际上只须要把集合外的着色代码里面的颜色份量改为与迭代数iter相关的形式：

//...
else
{
    glColor3f(static_cast<float>(iter) / maxIter, (0.5f*iter) / maxIter, 0.15f*iter / maxIter);
    diverge = false;
}

绘制结果：

值得注意的是，因为迭代次数是一个离散的值，所以整个图片的着色显得并非那般的丝滑和连续，在其余文章中还有提到过对数化将使着色更加柔和漂亮的方法，这里再也不阐述。

当指定参数maxIter的次数不一样，造成的图形也不一样，而且，随着迭代次数的增大， 形状逐渐趋于稳定，更接近于标准的曼集图样：