数学科普

自绘Mandelbrot集合（七）

作者：安迁

十六、连续的迭代次数

还剩下的问题有两个。一是靠近右边斜线那里的孔雀羽眼看起来和左边那些有点不同，它们的中心是黑的。这其实是我们在第十节中遇到的老问题。我们把迭代次数上限定得太低，以至于许多不属于Mandelbrot集合的点也被算在集合内部，于是被画成了黑色。

第二个问题是，如果观察左边大色块部分，会发现颜色呈不连续的条带状。这是因为迭代次数总是一个整数，所以逃逸次数（即第一次使得zn超出逃逸半径的迭代次数n）函数是离散的。当然这还是一个口味问题，也许有人就是觉得平滑过渡的颜色没有不连续的条带状颜色带漂亮。但是如果我们希望消除这些条带，希望颜色能够平滑过渡，那该怎么处理？

如果有两个复数c1和c2，都花了n次迭代逃出规定的逃逸半径（我们前面都取逃逸半径为2，于是其平方为4），按照前面离散的迭代次数的定义，它们逃出逃逸半径的速度是一样的，于是在调色板中取到的颜色也完全相同。我们希望有一种方法更细致地来区分出它们逃逸速度的快慢。

常见的处理方法是使用所谓的“连续的逃逸次数”函数来取代我们前面使用的离散的逃逸次数函数。先设定足够大的逃逸半径，对一个复数c，令n是第一次使得zn超出逃逸半径的迭代次数，定义它的“连续的逃逸次数”为 　　d = ｎ+１ - log(log(|zn|)) /log(2) 虽然这个定义取决于预先设定的逃逸半径，但可以证明，当这个预先定义的逃逸半径趋向于无穷大时，对于固定的复数c，上面这个式子趋向于一个固定值。所以如果我们只要确定一个较大的逃逸半径，我们就能近似地计算出它的“连续的逃逸次数”，其值不取决于所选取的具体的逃逸半径值。

假设在逃逸半径为2时，对于某复数c，离散的逃逸次数为n，也就是说|zn|>2但并不超过很多，那么如果我们取比较大的逃逸半径比如说为2k，其中k是一个整数，那么从zn开始还需要多少次才能超出这个大的逃逸半径？|zn+1|=|zn2+c|，如果c的模较小的话，|zn+1|≈|zn2|=|zn|2；类似地|zn+k|约等于|zn|2k，也就是大约还要k次迭代才能超出逃逸半径2k。按照上面的“连续的逃逸次数”的定义， 　　d = ｎ+ k +１ - log(log(|zn+k|)) /log(2) 　　　≈ ｎ+ k +１ - log(log(|zn|2k)) /log(2) 　　　= n + k + 1 - log(2klog(|zn|))/log(2) 　　　= n + k + 1 - k - log(|zn|)/log(2) 　　　= n + 1 - log(|zn|)/log(2) 　　　≈ n + 1 - log(2)/log(2) 　　　= n 所以这个“连续的逃逸次数”又是可以和原来的离散的逃逸次数可以相比的，可以看作是原来的定义的精细化。

按照上述公式，将逃逸半径的平方ESCAPERADIUS增加到128.0，将迭代次数上限增加到10000，并将在源程序的第(7)部分中使用上述“连续的逃逸次数”函数： Java版本

    // (1)可调参数
    double ESCAPERADIUS = 128.0;
    int MAXITERNUMBER = 10000;

    ……

    // (7)迭代计算
    double iter(double cx, double cy) {
        double x = 0;
        double y = 0;
        double newx;
        double newy;

        double smodz = 0;
        int i = 0;
        while (i < MAXITERNUMBER) {
            newx = x * x - y * y + cx;
            newy = 2 * x * y + cy;
            x = newx;
            y = newy;
            i++;

            smodz = x * x + y * y;
            if (smodz >= ESCAPERADIUS) {
                double d = i + 1 - Math.log(Math.log(smodz) * 0.5) / Math.log(2);
                return d;
            }
        }

        return -1.0;
    }

JavaScript版本

    // (1)可调参数
    var ESCAPERADIUS = 128.0;
    var MAXITERNUMBER = 10000;

    ……

    // (7)迭代计算
    function iter(cx, cy) {
        var x = 0;
        var y = 0;
        var newx;
        var newy;

        var smodz = 0;
        var i = 0;
        while (i < MAXITERNUMBER) {
            newx = x * x - y * y + cx;
            newy = 2 * x * y + cy;
            x = newx;
            y = newy;
            i++;

            smodz = x * x + y * y;
            if (smodz >= ESCAPERADIUS) {
                var d = i + 1 - Math.log(Math.log(smodz) * 0.5) / Math.log(2);
                return d;
            }
        }

        return -1.0;
    }

注意smodz计算的其实是前面所说的|zn|的平方，所以在log中需要乘以0.5。

运行得到（因为增大了逃逸半径和迭代次数上限，所需要的计算时间更长了，如有需要请耐心等待几分钟）：

-0.743030 + 0.126433i @ 0.016110 /0.75
-0.743030 + 0.126433i @ 0.016110 /0.75
我们终于得到了符合要求的图像。

十七、完整的源码

为了防止在前面修改过程中可能有操作失误，本节给出修改完成后最终完整的源码。懒得看上述一步步修改过程的读者也可以拿最后这个完整的版本来测试。容易验证最终两个版本的源程序长度均为一百五十行左右。 Java版本

package mandelbrot;

import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;

import javax.imageio.ImageIO;

public class Mandelbrot {
    // (1)可调参数
    double ESCAPERADIUS = 128.0;
    int MAXITERNUMBER = 10000;

    // -0.743030 + 0.126433i @ 0.016110 /0.75
    double OX = -0.743030;
    double OY = 0.126433;
    double WIDTH = 0.016110;
    double RATIO = 0.75;

    int IMAGEWIDTH = 800;

    // (2)计算所得参数
    int IMAGEHEIGHT = (int) (IMAGEWIDTH * RATIO);
    double PIXELSIZE = WIDTH / IMAGEWIDTH;
    double COFFSET = IMAGEWIDTH % 2 == 0 ? (IMAGEWIDTH / 2) - 0.5 : (IMAGEWIDTH / 2);
    double ROFFSET = IMAGEHEIGHT % 2 == 0 ? (IMAGEHEIGHT / 2) - 0.5 : (IMAGEHEIGHT / 2);

    // (3)图象缓存
    BufferedImage image = new BufferedImage(IMAGEWIDTH, IMAGEHEIGHT, BufferedImage.TYPE_INT_RGB);

    // (4)程序入口
    public static void main(String[] args) {
        new Mandelbrot().draw();
    }

    // (5)主程序
    void draw() {
        for (int row = 0; row < IMAGEHEIGHT; row++) {
            for (int col = 0; col < IMAGEWIDTH; col++) {
                int color = calcColor(col, row);
                drawColor(col, row, color);
            }
        }

        saveImage();
    }

    // (6)计算颜色
    int calcColor(int col, int row) {
        double cx = (col - COFFSET) * PIXELSIZE + OX;
        double cy = (row - ROFFSET) * PIXELSIZE + OY;
        int r = 0;
        int g = 0;
        int b = 0;
        for (int i = -1; i <= 1; i++) {
            for (int j = -1; j <= 1; j++) {
                double d = iter(cx + i * PIXELSIZE / 3, cy + j * PIXELSIZE / 3);
                int c = getColor(d);
                r += (c >> 16) & 0xFF;
                g += (c >> 8) & 0xFF;
                b += c & 0xFF;
            }
        }

        r /= 9;
        g /= 9;
        b /= 9;
        return (r << 16) | (g << 8) | b;
    }

    // (7)迭代计算
    double iter(double cx, double cy) {
        double x = 0;
        double y = 0;
        double newx;
        double newy;

        double smodz = 0;
        int i = 0;
        while (i < MAXITERNUMBER) {
            newx = x * x - y * y + cx;
            newy = 2 * x * y + cy;
            x = newx;
            y = newy;
            i++;

            smodz = x * x + y * y;
            if (smodz >= ESCAPERADIUS) {
                double d = i + 1 - Math.log(Math.log(smodz) * 0.5) / Math.log(2);
                return d;
            }
        }

        return -1.0;
    }

    // (8)调色盘
    int getColor(double d) {
        if (d >= 0) {
            double k = 0.021 * (d - 1 + Math.log(Math.log(128)) / Math.log(2));
            k = Math.log(1 + k) - 29.0 / 400;

            k = k - Math.floor(k);
            k *= 400;
            if (k < 63.0) {
                return interpolation(k / 63.0, 0x000764, 0x206BCB);
            } else if (k < 167.0) {
                return interpolation((k - 63.0) / (167.0 - 63.0), 0x206BCB, 0xEDFFFF);
            } else if (k < 256.0) {
                return interpolation((k - 167.0) / (256.0 - 167.0), 0xEDFFFF, 0xFFAA00);
            } else if (k < 342.0) {
                return interpolation((k - 256.0) / (342.0 - 256.0), 0xFFAA00, 0x310230);
            } else {
                return interpolation((k - 342.0) / (400.0 - 342.0), 0x310230, 0x000764);
            }
        } else {
            return 0x000000;
        }
    }

    int interpolation(double f, int c0, int c1) {
        int r0 = (c0 >> 16) & 0xFF;
        int g0 = (c0 >> 8) & 0xFF;
        int b0 = c0 & 0xFF;
        int r1 = (c1 >> 16) & 0xFF;
        int g1 = (c1 >> 8) & 0xFF;
        int b1 = c1 & 0xFF;
        int r = (int) ((1 - f) * r0 + f * r1 + 0.5);
        int g = (int) ((1 - f) * g0 + f * g1 + 0.5);
        int b = (int) ((1 - f) * b0 + f * b1 + 0.5);
        return (r << 16) | (g << 8) | b;
    }

    // (9)在图像像素(col, row)处画上颜色rgb
    void drawColor(int col, int row, int rgb) {
        image.setRGB(col, IMAGEHEIGHT - row - 1, rgb);
    }

    // (10)保存图像
    void saveImage() {
        try {
            ImageIO.write(image, "png", new File("c:\\temp\\image.png"));
        } catch (IOException e) {
            e.printStackTrace();
        }
    }
}

JavaScript版本

<!DOCTYPE html>
<html>
<head>
<meta charset="UTF-8">
<title>Mandelbrot</title>
<script>
// (4)程序入口
window.onload = function() {
    // (1)可调参数
    var ESCAPERADIUS = 128.0;
    var MAXITERNUMBER = 10000;

    // -0.743030 + 0.126433i @ 0.016110 /0.75
    var OX = -0.743030;
    var OY = 0.126433;
    var WIDTH = 0.016110;
    var RATIO = 0.75;

    var IMAGEWIDTH = 800;

    // (2)计算所得参数
    var IMAGEHEIGHT = IMAGEWIDTH * RATIO;
    var PIXELSIZE = WIDTH / IMAGEWIDTH;
    var COFFSET = IMAGEWIDTH % 2 == 0 ? (IMAGEWIDTH / 2) - 0.5 : (IMAGEWIDTH / 2);
    var ROFFSET = IMAGEHEIGHT % 2 == 0 ? (IMAGEHEIGHT / 2) - 0.5 : (IMAGEHEIGHT / 2);

    // (3)图象缓存
    var canvas = document.getElementById("image");
    var context = canvas.getContext("2d");
    canvas.width = IMAGEWIDTH;
    canvas.height = IMAGEHEIGHT;
    var imagedata = context.createImageData(IMAGEWIDTH, IMAGEHEIGHT);

    // (5)主程序
    for (row = 0; row < IMAGEHEIGHT; row++) {
        for (col = 0; col < IMAGEWIDTH; col++) {
            var color = calcColor(col, row);
            drawColor(col, row, color);
        }
    }

    saveImage();

    // (6)计算颜色
    function calcColor(col, row) {
        var cx = (col - COFFSET) * PIXELSIZE + OX;
        var cy = (row - ROFFSET) * PIXELSIZE + OY;
        var r = 0;
        var g = 0;
        var b = 0;
        for (i = -1; i <= 1; i++) {
            for (j = -1; j <= 1; j++) {
                var d = iter(cx + i * PIXELSIZE / 3, cy + j * PIXELSIZE / 3);
                var c = getColor(d);
                r += (c >> 16) & 0xFF;
                g += (c >> 8) & 0xFF;
                b += c & 0xFF;
            }
        }

        r /= 9;
        g /= 9;
        b /= 9;
        return (r << 16) | (g << 8) | b;
    }

    // (7)迭代计算
    function iter(cx, cy) {
        var x = 0;
        var y = 0;
        var newx;
        var newy;

        var smodz = 0;
        var i = 0;
        while (i < MAXITERNUMBER) {
            newx = x * x - y * y + cx;
            newy = 2 * x * y + cy;
            x = newx;
            y = newy;
            i++;

            smodz = x * x + y * y;
            if (smodz >= ESCAPERADIUS) {
                var d = i + 1 - Math.log(Math.log(smodz) * 0.5) / Math.log(2);
                return d;
            }
        }

        return -1.0;
    }

    // (8)调色盘
    function getColor(d) {
        if (d >= 0) {
            var k = 0.021 * (d - 1 + Math.log(Math.log(128)) / Math.log(2));
            k = Math.log(1 + k) - 29.0 / 400;

            k = k - Math.floor(k);
            k *= 400;
            if (k < 63.0) {
                return interpolation(k / 63.0, 0x000764, 0x206BCB);
            } else if (k < 167.0) {
                return interpolation((k - 63.0) / (167.0 - 63.0), 0x206BCB, 0xEDFFFF);
            } else if (k < 256.0) {
                return interpolation((k - 167.0) / (256.0 - 167.0), 0xEDFFFF, 0xFFAA00);
            } else if (k < 342.0) {
                return interpolation((k - 256.0) / (342.0 - 256.0), 0xFFAA00, 0x310230);
            } else {
                return interpolation((k - 342.0) / (400.0 - 342.0), 0x310230, 0x000764);
            }
        } else {
            return 0x000000;
        }
    }

    function interpolation(f, c0, c1) {
        var r0 = (c0 >> 16) & 0xFF;
        var g0 = (c0 >> 8) & 0xFF;
        var b0 = c0 & 0xFF;
        var r1 = (c1 >> 16) & 0xFF;
        var g1 = (c1 >> 8) & 0xFF;
        var b1 = c1 & 0xFF;
        var r = Math.floor((1 - f) * r0 + f * r1 + 0.5);
        var g = Math.floor ((1 - f) * g0 + f * g1 + 0.5);
        var b = Math.floor ((1 - f) * b0 + f * b1 + 0.5);
        return (r << 16) | (g << 8) | b;
    }

    // (9)在图像像素(col, row)处画上颜色rgb
    function drawColor(col, row, rgb) {
        var pindex = ((IMAGEHEIGHT - row - 1) * IMAGEWIDTH + col) * 4;
        imagedata.data[pindex] = (rgb>>16) & 0xFF;
        imagedata.data[pindex+1] = (rgb>>8) & 0xFF;
        imagedata.data[pindex+2] = rgb & 0xFF;
        imagedata.data[pindex+3] = 255;
    }

       // (10)保存图像
    function saveImage() {
        context.putImageData(imagedata, 0, 0);
    }
};

</script>
</head>
<body>
<canvas id="image" width="0" height="0"></canvas>
</body>
</html>

十八、进一步的话题

关于Mandelbrot集合图像绘制的具体编程的讨论，本文到此已经全部完成了。下面给希望进一步了解这方面内容的读者提供一些建议。

本文介绍的按照逃逸次数来绘制Mandelbrot集合图像的方法也许是比较简单的，但并不是最精细的。维基百科中还介绍了其他的算法，最为精细的是“距离估计法”（Distance Estimation Method）。它通过Koebe定理来直接估计一个点到Mandelbrot集合的距离，而非象逃逸次数法那样间接地通过逃逸次数来估计。这对于希望直接看到Mandelbrot集合图像，而非通过看到其临近结构从而间接看到它的人来说更有用。这一点在观察所谓的“细丝”结构，以及验证Mandelbrot集合图像确为连通集合，而需要非黑即白地画出Mandelbrot集合图像时尤为明显。

-1.770119403588788 + 0.00411094010822248 @ 1e-9 /1 逃逸次数法绘图
-1.770119403588788 + 0.00411094010822248 @ 1e-9 /1 逃逸次数法绘图
-1.770119403588788 + 0.00411094010822248 @ 1e-9 /1 距离估计法绘图
-1.770119403588788 + 0.00411094010822248 @ 1e-9 /1 距离估计法绘图
当然，要了解距离估计法绘图则需要更进一步的数学工具，在此就不具体介绍了。

本文的程序没有用户交互功能，只能绘制指定区域坐标的图像。如果事先不知道区域坐标而想自己探索Mandelbrot集合有趣的地点的话，本程序就很不实用了。好在网上有许多此类程序，比如英文维基百科中介绍的Web Mandelbrot网页就是一个很好的工具，在本文的写作过程中帮了我不少忙。具体链接请见参考文献。

本文的分形图像基本都是宽度为800个像素，最后的反锯齿也只用3x3超级采样。但分形图像以大和精细为美，如果真要作装饰（比如电脑桌面），需要增加宽度、逃逸半径、迭代上限以及超级采样次数的设置。此时的计算量将会大大增加。一般应该使用专业软件绘制，但是自己动手也未尝不可。我的建议是使用低设置来测试调色板和图像的区域，确定以后再选择高设置来绘制最终结果。另外我们注意到，在第(5)部分主程序的循环，以及第(6)部分计算颜色中，我们计算出一个像素中心点的逃逸次数后立刻使用它来计算出像素颜色并绘制。但其实一个像素中心点的逃逸次数和它的颜色是独立的，后者随着不同调色板而变，但前者总是一样的。所以在象Java这样有保存中间计算结果功能（比如使用Java序列化技术）的语言中，对固定的区域坐标，我们完全可以先计算好每个像素中心点逃逸次数并保存。在以后使用不同调色板时就无需再重复计算，只需直接读取已完成的计算结果即可，这可以省下许多计算时间。而JavaScript这样运行速度为弱项的语言来绘制大而精细的图像则不太实用。

试用本文的程序来绘制这个区域坐标的图像：-0.744539860355908380 + 0.121723773894424824i @ 6.0e-15 /1。我们得到下面这个模糊的图像：

-0.744539860355908380 + 0.121723773894424824i @ 6.0e-15 /1
-0.744539860355908380 + 0.121723773894424824i @ 6.0e-15 /1
这是因为64比特双精度浮点运算到了它的能力极限，这个极限大约在画面宽度为1e-14处。所有只使用64比特双精度浮点运算的分形绘图程序都有这个问题，包括上面所说的Web Mandelbrot网页。要突破这个极限，就必须使用更高精度的浮点数运算。理论上我们有任意精度的浮点运算工具，但是如果对每一点都进行这样的运算，计算时间会极其漫长。有一种微扰算法，可以只需对一点进行高精度浮点运算，而对其周围点只进行普通双精度浮点运算而得到令人满意的结果，第三节中提到的Kalles Fraktaler软件就是利用这个算法来达到令人乍舌的Mandelbrot集合超深放大的图像。

更多的对程序的优化方法可参考维基百科Mandelbrot Set条“优化”一节。

参考文献：

[1] Javier Barrallo, and Damien M. Jones, Coloring Algorithms for Dynamical Systems in the Complex Plane, Visual Mathematics 1.4 (1999)

[2] Karol Guciek, Web Mandelbrot guciek.github.io/web_mandelbrot.html

[3] 维基百科Mandelbrot set条目 https://en.wikipedia.org/wiki/Mandelbrot_set

<（六）
（待续）