Lode 的计算机图形学教程

Julia 集与 Mandelbrot 集

目录

返回目录

简介

Julia 集和 Mandelbrot 集是复平面上相当著名的集合，能生成那些美丽的无限细节图像。它们如此漂亮，甚至有人用它们创作艺术品。从定义上说，Mandelbrot 集并不是真正的分形，但它具有半自相似性，仍然展示出无限细节，因此通常也被称为分形。

Julia 集的研究由 Gaston Julia 本人于 1917 年完成，但他当时没有计算机来实际绘制它。这个话题没有引起太多关注，直到 1970 年代，Benoît Mandelbrot 在计算机上绘制了 Julia 集，并发现了 Mandelbrot 集。

本教程只涵盖理论基础（关于 Julia 集和 Mandelbrot 集还有很多内容可以讲），以及绘制它们的代码。此外还包含一些非常实用的分形查看器代码，可以实时缩放和移动。

Julia 集

本章首先尝试解释 Julia 集背后的一些数学公式，具体来说是二次 Julia 集。需要具备复数知识，如有需要可以阅读关于复数的附录。

那么如何生成如此美丽的分形？简而言之：对每个像素，在复平面上迭代 znew = zold² + c，直到它离开以原点为圆心、半径为 2 的圆。迭代次数即为像素的颜色。

屏幕将代表复平面的一部分，在以原点为圆心、半径为 2 的圆内。对于一个像素，x 坐标代表其复坐标的实部，y 坐标代表虚部。

对于 Julia 集，对每个像素应用一个迭代复函数。该函数为 newz = oldz² + c，z 和 c 都是复数。z 最初是像素的坐标，然后在每次迭代中不断更新：每次迭代，上一次迭代的"newz"作为"oldz"使用。

持续迭代该函数，根据初始条件（像素），z 要么趋向无穷，要么永远留在以复平面原点为圆心、半径为 2 的圆内。永远留在圆内的点就属于 Julia 集。因此持续迭代函数，直到 z 到原点 (0,0) 的距离大于 2。还需设置最大迭代次数（例如 256），否则计算机会陷入无限循环。

像素的颜色值将等于在 z 到原点的距离大于 2 之前需要迭代的次数。公式中的常数 c 可以是任意值，只要它也在半径为 2 的圆内。不同的 c 值会产生不同的 Julia 集。有些 Julia 集是连通的，有些不是。Mandelbrot 集是所有能生成连通 Julia 集的点 c 的集合。

以下是计算示例：

首先可以为函数选择常数 c，选择哪个值将决定分形的形状。本例取 c = (-0.5, 0.5)，即实部为 -0.5，虚部为 0.5。

假设我们正在计算 256×256 屏幕上像素 (256, 192) 的颜色。首先将坐标变换到 -1 到 1 之间（如果在分形中缩放或移动则需要不同的变换）：坐标变为 (1, 0.5)，即 p = 1 + 0.5i。

现在第一次应用函数：

z = p² + c
  = (1 + 0.5i)² - 0.5 + 0.5i
  = 1 + 2*0.5i + 0.25*i² - 0.5 + 0.5i
  = 0.5 + 1.5i - 0.25
  = 0.25 + 1.5i

所以 z = (0.25, 1.5)，z 到原点的距离 = sqrt(0.25*0.25 + 1.5*1.5) = 1.52069...，仍小于 2。

现在将计算出的 z 再次代入函数计算下一个 z：

z = (0.25 + 1.5i)² - 0.5 + 0.5i
  = 0.0625 + 2*0.375i - 2.25 - 0.5 + 0.5i
  = -2.6875 + 1.25i

距离现在为 8.78515625，超出了半径为 2 的圆，因此该点在 Julia 集之外。迭代次数仅为 2，所以像素的颜色值为 2，计算完成。某些初始值会给出超过 256 次迭代，根据设置的最大迭代次数，可以选择停止或继续。

迭代次数越多，深度缩放时 Julia 集看起来越详细，但需要的计算量也越多。数值精度越高，可以缩放的深度越大而不会出现像素化。

以下是绘制 Julia 集的程序源代码。将其放入 main.cpp 文件的 "int main(int argc, char *argv[])" 函数中。代码中的蓝色注释解释了大部分代码。代码不使用复数，而是使用普通浮点数，实部和虚部简单地分别计算，就像手工计算一样。

int main(int argc, char *argv[])
{
  screen(400, 300, 0, "Julia Set"); //make larger to see more detail!

  //each iteration, it calculates: new = old*old + c, where c is a constant and old starts at current pixel
  double cRe, cIm;           //real and imaginary part of the constant c, determinate shape of the Julia Set
  double newRe, newIm, oldRe, oldIm;   //real and imaginary parts of new and old
  double zoom = 1, moveX = 0, moveY = 0; //you can change these to zoom and change position
  ColorRGB color; //the RGB color value for the pixel
  int maxIterations = 300; //after how much iterations the function should stop

  //pick some values for the constant c, this determines the shape of the Julia Set
  cRe = -0.7;
  cIm = 0.27015;

  //loop through every pixel
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    //calculate the initial real and imaginary part of z, based on the pixel location and zoom and position values
    newRe = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
    newIm = (y - h / 2) / (0.5 * zoom * h) + moveY;
    //i will represent the number of iterations
    int i;
    //start the iteration process
    for(i = 0; i < maxIterations; i++)
    {
      //remember value of previous iteration
      oldRe = newRe;
      oldIm = newIm;
      //the actual iteration, the real and imaginary part are calculated
      newRe = oldRe * oldRe - oldIm * oldIm + cRe;
      newIm = 2 * oldRe * oldIm + cIm;
      //if the point is outside the circle with radius 2: stop
      if((newRe * newRe + newIm * newIm) > 4) break;
    }
    //use color model conversion to get rainbow palette, make brightness black if maxIterations reached
    color = HSVtoRGB(ColorHSV(i % 256, 255, 255 * (i < maxIterations)));
    //draw the pixel
    pset(x, y, color);
  }
  //make the Julia Set visible and wait to exit
  redraw();
  sleep();
  return 0;
}

迭代次数参数用于 HSV 颜色模型的"色相"。色相的优点是它是循环的，因此无论最大迭代次数是多少，基于色相的调色板都能生成漂亮的连续值。

结果如下：


Julia 集探索器

你可以修改上面代码中的"zoom"、"moveX"和"moveY"值来缩放到特定位置，但更好的做法是在程序运行时实时操作，例如用方向键移动，用小键盘 + 和 - 键缩放。更好的是还能用小键盘方向键实时改变 cRe 和 cIm 的值来改变 Julia 集的形状。

编写能做到这些的程序非常简单，只需使用 SDL 按键来改变这些变量的值，并在循环中重新绘制 Julia 集即可！以下代码实现了所有这些功能，还在屏幕上显示所有变量的值，这样你就能确切知道 Julia 集中漂亮部分的坐标。还可以修改"maxIterations"等值……代码中的注释会再次解释一切。没有新的计算机图形学代码，只是改变参数的输入按键：

int main(int argc, char *argv[])
{
  screen(320, 240, 0, "Julia  Explorer");

  //each iteration, it calculates: new = old*old + c, where c is a constant and old starts at current pixel
  double cRe, cIm;           //real and imaginary part of the constant c, determines shape of the Julia Set
  double newRe, newIm, oldRe, oldIm;   //real and imaginary parts of new and old
  double zoom=1, moveX=0, moveY=0; //you can change these to zoom and change position
  ColorRGB color; //the RGB color value for the pixel
  int maxIterations=128; //after how much iterations the function should stop

  double time, oldTime, frameTime; //current and old time, and their difference (for input)
  int showText=0;

  //pick some values for the constant c, this determines the shape of the Julia Set
  cRe = -0.7;
  cIm = 0.27015;

  //begin the program loop
  while(!done())
  {
    //draw the fractal
    for(int y = 0; y < h; y++)
    for(int x = 0; x < w; x++)
    {
      //calculate the initial real and imaginary part of z, based on the pixel location and zoom and position values
      newRe = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
      newIm = (y - h / 2) / (0.5 * zoom * h) + moveY;
      //i will represent the number of iterations
      int i;
      //start the iteration process
      for(i = 0; i < maxIterations; i++)
      {
        //remember value of previous iteration
        oldRe = newRe;
        oldIm = newIm;
        //the actual iteration, the real and imaginary part are calculated
        newRe = oldRe * oldRe - oldIm * oldIm + cRe;
        newIm = 2 * oldRe * oldIm + cIm;
        //if the point is outside the circle with radius 2: stop
        if((newRe * newRe + newIm * newIm) > 4) break;
      }
      //use color model conversion to get rainbow palette, make brightness black if maxIterations reached
      color = HSVtoRGB(ColorHSV(i % 256, 255, 255 * (i < maxIterations)));
      //draw the pixel
      pset(x, y, color);
    }

    //print the values of all variables on screen if that option is enabled
    if(showText <= 1)
    {
      print("X:", 1, 1, RGB_White, 1); print(moveX, 17, 1, RGB_White, 1);
      print("Y:", 1, 9, RGB_White, 1); print(moveY, 17, 9, RGB_White, 1);
      print("Z:", 1, 17, RGB_White, 1); print(zoom, 17, 17, RGB_White, 1);
      print("R:", 1, 25, RGB_White, 1); print(cRe, 17, 25, RGB_White, 1);
      print("I:", 1, 33, RGB_White, 1); print(cIm, 17, 33, RGB_White, 1);
      print("N:", 1, 41, RGB_White, 1); print(maxIterations, 17, 41, RGB_White, 1);
    }
    //print the help text on screen if that option is enabled
    if(showText == 0)
    {
      print("Arrows move (X,Y), Keypad +,- zooms (Z)", 1, h - 33, RGB_White, 1);
      print("Keypad arrows change shape (R,I)     ", 1, h - 25, RGB_White, 1);
      print("Keypad *,/ changes iterations (N)    ", 1, h - 17, RGB_White, 1);
      print("a to z=presets (qwerty), F1=cycle texts", 1, h - 9, RGB_White, 1);
    }
    redraw();

    //get the time and old time for time dependent input
    oldTime = time;
    time = getTicks();
    frameTime = time - oldTime;
    readKeys();
    //ZOOM keys
    if(keyDown(SDLK_KP_PLUS))  {zoom *= pow(1.001, frameTime);}
    if(keyDown(SDLK_KP_MINUS)) {zoom /= pow(1.001, frameTime);}
    //MOVE keys
    if(keyDown(SDLK_DOWN))  {moveY += 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_UP))  {moveY -= 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_RIGHT)) {moveX += 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_LEFT))  {moveX -= 0.0003 * frameTime / zoom;}
    //CHANGE SHAPE keys
    if(keyDown(SDLK_KP2)) {cIm += 0.0002 * frameTime / zoom;}
    if(keyDown(SDLK_KP8)) {cIm -= 0.0002 * frameTime / zoom;}
    if(keyDown(SDLK_KP6)) {cRe += 0.0002 * frameTime / zoom;}
    if(keyDown(SDLK_KP4)) {cRe -= 0.0002 * frameTime / zoom;}
    //keys to change number of iterations
    if(keyPressed(SDLK_KP_MULTIPLY)) {maxIterations *= 2;}
    if(keyPressed(SDLK_KP_DIVIDE))   {if(maxIterations > 2) maxIterations /= 2;}
    //key to change the text options
    if(keyPressed(SDLK_F1)) {showText++; showText %= 3;}
  }
}

现在你可以探索每种可能的 Julia 集的所有细节！用小键盘数字找到漂亮的形状，然后用方向键移动到 Julia 集的边界或有趣的位置，开始缩放以查看更多细节。缩放后可以按小键盘"*"键查看更多细节。

注意：在当前版本的 gcc 中，如果混用 float 和 double，pow 函数无法正常工作，因此请确保全部使用 double 或全部使用 float。

以下是一些截图：

几种不同的形状：


缩放：


增加迭代次数：


Mandelbrot 集

缩放时，在 Julia 集中会不断看到相同的细节——毕竟它是分形。Mandelbrot 集不是完全自相似的，只是半自相似，因此在 Mandelbrot 集中缩放时会出现更多惊喜。

Mandelbrot 集代表所有能使 Julia 集连通的复数点 c，即所有包含原点的 Julia 集。生成 Mandelbrot 集时使用与 Julia 集相同的迭代函数，只是这次 c 代表像素的位置，z 从 (0,0) 开始。

以下代码与 Julia 集绘制器非常相似，只是输出带彩虹调色板的 Mandelbrot 集。修改的部分用粗体标出。

int main(int argc, char *argv[])
{
  screen(400, 300, 0, "Mandelbrot Set"); //make larger to see more detail!

  //each iteration, it calculates: newz = oldz*oldz + p, where p is the current pixel, and oldz stars at the origin
  double pr, pi;           //real and imaginary part of the pixel p
  double newRe, newIm, oldRe, oldIm;   //real and imaginary parts of new and old z
  double zoom = 1, moveX = -0.5, moveY = 0; //you can change these to zoom and change position
  ColorRGB color; //the RGB color value for the pixel
  int maxIterations = 300;//after how much iterations the function should stop

  //loop through every pixel
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    //calculate the initial real and imaginary part of z, based on the pixel location and zoom and position values
    pr = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
    pi = (y - h / 2) / (0.5 * zoom * h) + moveY;
    newRe = newIm = oldRe = oldIm = 0; //these should start at 0,0
    //"i" will represent the number of iterations
    int i;
    //start the iteration process
    for(i = 0; i < maxIterations; i++)
    {
      //remember value of previous iteration
      oldRe = newRe;
      oldIm = newIm;
      //the actual iteration, the real and imaginary part are calculated
      newRe = oldRe * oldRe - oldIm * oldIm + pr;
      newIm = 2 * oldRe * oldIm + pi;
      //if the point is outside the circle with radius 2: stop
      if((newRe * newRe + newIm * newIm) > 4) break;
    }
    //use color model conversion to get rainbow palette, make brightness black if maxIterations reached
    color = HSVtoRGB(ColorHSV(i % 256, 255, 255 * (i < maxIterations)));
     //draw the pixel
     pset(x, y, color);
  }
  //make the Mandelbrot Set visible and wait to exit
  redraw();
  sleep();
  return 0;
}

是不是很漂亮：


Mandelbrot 集探索器

以下是允许你在 Mandelbrot 集中移动和缩放的完整程序代码。

int main(int argc, char *argv[])
{
  screen(320, 240, 0, "Mandelbrot Explorer");

  //each iteration, it calculates: new = old*old + c, where c is a constant and old starts at current pixel
  double pr, pi;           //real and imaginary part of the pixel p
  double newRe, newIm, oldRe, oldIm;   //real and imaginary parts of new and old
  double zoom = 1, moveX = -0.5, moveY = 0; //you can change these to zoom and change position
  ColorRGB color; //the RGB color value for the pixel
  int maxIterations = 128; //after how much iterations the function should stop

  double time, oldTime, frameTime; //current and old time, and their difference (for input)
  int showText = 0;

  //begin main program loop
  while(!done())
  {
    //draw the fractal
    for(int y = 0; y < h; y++)
    for(int x = 0; x < w; x++)
    {
      //calculate the initial real and imaginary part of z, based on the pixel location and zoom and position values
      pr = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
      pi = (y - h / 2) / (0.5 * zoom * h) + moveY;
      newRe = newIm = oldRe = oldIm = 0; //these should start at 0,0
      //i will represent the number of iterations
      int i;
      //start the iteration process
      for(i = 0; i < maxIterations; i++)
      {
        //remember value of previous iteration
        oldRe = newRe;
        oldIm = newIm;
        //the actual iteration, the real and imaginary part are calculated
        newRe = oldRe * oldRe - oldIm * oldIm + pr;
        newIm = 2 * oldRe * oldIm + pi;
        //if the point is outside the circle with radius 2: stop
        if((newRe * newRe + newIm * newIm) > 4) break;
      }
      //use color model conversion to get rainbow palette, make brightness black if maxIterations reached
      color = HSVtoRGB(ColorHSV(i % 256, 255, 255 * (i < maxIterations)));
      //draw the pixel
      pset(x, y, color);
    }

    //print the values of all variables on screen if that option is enabled
    if(showText <= 1)
    {
      print("X:", 1, 1, RGB_White, 1); print(moveX, 17, 1, RGB_White, 1);
      print("Y:", 1, 9, RGB_White, 1); print(moveY, 17, 9, RGB_White, 1);
      print("Z:", 1, 17, RGB_White, 1); print(zoom, 17, 17, RGB_White, 1);
      print("N:", 1, 25, RGB_White, 1); print(maxIterations, 17, 25, RGB_White, 1);
    }
    //print the help text on screen if that option is enabled
    if(showText == 0)
    {
      print("Arrows move (X,Y), Keypad +,- zooms (Z)", 1, h - 25, RGB_White, 1);
      print("Keypad *,/ changes iterations (N)    ", 1, h - 17, RGB_White, 1);
      print("a to z=presets (qwerty), F1=cycle texts", 1, h - 9, RGB_White, 1);
    }
    redraw();

    //get the time and old time for time dependent input
    oldTime = time;
    time = getTicks();
    frameTime = time - oldTime;
    readKeys();
    //ZOOM keys
    if(keyDown(SDLK_KP_PLUS))  {zoom *= pow(1.001, frameTime);}
    if(keyDown(SDLK_KP_MINUS)) {zoom /= pow(1.001, frameTime);}
    //MOVE keys
    if(keyDown(SDLK_DOWN))  {moveY += 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_UP))  {moveY -= 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_RIGHT)) {moveX += 0.0003 * frameTime / zoom;}
    if(keyDown(SDLK_LEFT))  {moveX -= 0.0003 * frameTime / zoom;}
    //keys to change number of iterations
    if(keyPressed(SDLK_KP_MULTIPLY)) {maxIterations *= 2;}
    if(keyPressed(SDLK_KP_DIVIDE))   {if(maxIterations > 2) maxIterations /= 2;}
    //key to change the text options
    if(keyPressed(SDLK_F1)) {showText++; showText %= 3;}
  }
  return 0;
}

通过移动和缩放可以得到非常漂亮的图像。图片上保留了参数值，以便你看到生成这些图片所需的坐标、缩放级别和最大迭代次数：

Mandelbrot 集的这部分称为海马谷（Seahorse Valley）：


深度缩放至海马谷的大小（放大 1779 倍）：


集合另一侧的细节：


几个"小 Mandelbrot"（minibrot），第一个只放大了 8463 倍，第二个放大了 419622325484 倍！这样的小 Mandelbrot 只有在允许足够多次迭代时才会可见。


这是使用更多迭代次数的效果：右图与左图相同，但最大迭代次数是左图的两倍：


这是在象谷（elephant valley）侧面进行的深度缩放，使用了非常高的迭代次数。象谷是 Mandelbrot 集右侧 X 轴上的尖锐形状。


这张图片是后来添加的，因为它非常漂亮：



这张图片放大倍数极大，已达到 64 位浮点数的数值极限，开始出现像素化。它放大了超过 10^18 倍，要查看更深的细节，需要更高精度的处理器或模拟无限精度数值：


如果你喜欢，这里有一张 1280×1024 像素的 Mandelbrot 集图片。

最后编辑：2004 年

版权所有 (c) 2004-2007 Lode Vandevenne，保留所有权利。