Lode 的计算机图形学教程

谢尔宾斯基分形

目录

简介

以瓦茨瓦夫·谢尔宾斯基(Waclaw Sierpinski)命名的分形有很多,他是一位生活在 1882 年至 1969 年间的波兰数学家。

这些分形包括谢尔宾斯基三角形、谢尔宾斯基地毯、谢尔宾斯基金字塔(谢尔宾斯基三角形的三维版本)以及谢尔宾斯基立方体(谢尔宾斯基地毯的三维版本)。本文将介绍其中的二维图形。

谢尔宾斯基三角形

谢尔宾斯基三角形,又称谢尔宾斯基垫片和谢尔宾斯基筛,可以按如下方式手工绘制:

从一个三角形开始。这是唯一一个正向的三角形,其余所有三角形都将是倒置的:


在这个三角形内部,画一个较小的倒置三角形。它的顶点应该恰好位于大三角形各边的中点:



现在,在 3 个向上的三角形中各画 3 个更小的三角形,同样将顶点置于向上三角形各边的中点:



现在有 9 个向上的三角形。在每一个三角形中,再次画更小的倒置三角形:



在 27 个向上的三角形中,再次画 27 个向下的三角形:



如此循环往复。

经过无限步骤之后,如果所有向上的三角形都被填充,就得到了谢尔宾斯基筛。

每一步都需要绘制更多的三角形。这是一个递归过程,用计算机也可以用同样的方式来绘制。

递归法

现在我们将用计算机来实现手工绘制的过程:编写一个三角形绘制函数,该函数会递归调用自身 3 次,直到达到 n 步递归为止。这个程序适用于任意初始三角形,不要求对称,唯一的条件是三角形的顶点位于屏幕范围内。

主函数负责初始化屏幕并调用 drawSierpinski 函数。drawSierpinski 函数本身只绘制一个三角形:初始的向上三角形。然后它会调用 subTriangle 函数,这才是真正的递归函数,负责绘制所有倒置的三角形。

subTriangle 函数根据传入的参数绘制一个倒置三角形的 3 个顶点。然后它会再次调用自身 3 次,以绘制 3 个更小的三角形。当然,这 3 个三角形需要使用新的顶点坐标,这些坐标需要计算得出。在下图中,如果黑色三角形是 subTriangle 函数已绘制的大三角形,那么 3 个红色三角形就是需要计算的新三角形:



大三角形的顶点为 a1、a2 和 a3。如图所示,其中一个小三角形的顶点为 b1、b2 和 b3。如果将这些点视为向量,则在已知点 a 的情况下,点 b 的计算公式为:

b3 = (a1 + a2) / 2,因为 b3 位于 a1 和 a2 的中点,即 a1 和 a2 的平均值!

b1 = b3 + (a1 - a3) / 2:仔细观察可以发现,点 b1 是点 b3 与向量 (a1 - a3) / 2 的和,除以 2 是因为小三角形对应边的长度是大三角形的一半。

最后,

b2 = b3 + (a2 - a3) / 2:这与上面的公式非常相似,只是换成了另一条边。

另外 2 个小三角形的处理方式类似。

在代码中,我们没有使用向量类,x 和 y 是分开的变量。由于向量加法的运算方式,我们只需用相同的公式分别对 x 和 y 各做一次相同的运算即可。

三角形顶点的坐标使用浮点数以获得更高的精度。

有了这些知识,就可以编写程序了,下面代码中的注释将解释其工作原理:

//Declaration of the drawSierpinski function. The coordinates are the 3 outer corners of the Sierpinski Triangle.
void drawSierpinski(float x1, float y1, float x2, float y2, float x3, float y3);
//Declaration of the subTriangle function, the coordinates are the 3 corners, and n is the number of recursions.
void subTriangle(int n, float x1, float y1, float x2, float y2, float x3, float y3);

//depth is the number of recursive steps
int depth = 7;

//The main function sets up the screen and then calls the drawSierpinski function
int main(int argc, char *argv[])
{
  screen(640, 480, 0, "Sierpinski Triangle");
  cls(RGB_White); //Make the background white
  drawSierpinski(10, h - 10, w - 10, h - 10, w / 2, 10); //Call the sierpinski function (works with any corners inside the screen)
  //After drawing the whole thing, redraw the screen and wait until the any key is pressed
  redraw();
  sleep();
  return(0);
}

//This function will draw only one triangle, the outer triangle (the only not upside down one), and then start the recursive function
void drawSierpinski(float x1, float y1, float x2, float y2, float x3, float y3)
{
    //Draw the 3 sides of the triangle as black lines
    drawLine(int(x1), int(y1), int(x2), int(y2), RGB_Black);
    drawLine(int(x1), int(y1), int(x3), int(y3), RGB_Black);
    drawLine(int(x2), int(y2), int(x3), int(y3), RGB_Black);

    //Call the recursive function that'll draw all the rest. The 3 corners of it are always the centers of sides, so they're averages
    subTriangle
    (
      1, //This represents the first recursion
      (x1 + x2) / 2, //x coordinate of first corner
      (y1 + y2) / 2, //y coordinate of first corner
      (x1 + x3) / 2, //x coordinate of second corner
      (y1 + y3) / 2, //y coordinate of second corner
      (x2 + x3) / 2, //x coordinate of third corner
      (y2 + y3) / 2  //y coordinate of third corner
    );
}

//The recursive function that'll draw all the upside down triangles
void subTriangle(int n, float x1, float y1, float x2, float y2, float x3, float y3)
{
  //Draw the 3 sides as black lines
  drawLine(int(x1), int(y1), int(x2), int(y2), RGB_Black);
  drawLine(int(x1), int(y1), int(x3), int(y3), RGB_Black);
  drawLine(int(x2), int(y2), int(x3), int(y3), RGB_Black);

  //Calls itself 3 times with new corners, but only if the current number of recursions is smaller than the maximum depth
  if(n < depth)
  {
    //Smaller triangle 1
    subTriangle
    (
      n+1, //Number of recursions for the next call increased with 1
      (x1 + x2) / 2 + (x2 - x3) / 2, //x coordinate of first corner
      (y1 + y2) / 2 + (y2 - y3) / 2, //y coordinate of first corner
      (x1 + x2) / 2 + (x1 - x3) / 2, //x coordinate of second corner
      (y1 + y2) / 2 + (y1 - y3) / 2, //y coordinate of second corner
      (x1 + x2) / 2, //x coordinate of third corner
      (y1 + y2) / 2  //y coordinate of third corner
    );
    //Smaller triangle 2
    subTriangle
    (
      n+1, //Number of recursions for the next call increased with 1
      (x3 + x2) / 2 + (x2 - x1) / 2, //x coordinate of first corner
      (y3 + y2) / 2 + (y2 - y1) / 2, //y coordinate of first corner
      (x3 + x2) / 2 + (x3 - x1) / 2, //x coordinate of second corner
      (y3 + y2) / 2 + (y3 - y1) / 2, //y coordinate of second corner
      (x3 + x2) / 2, //x coordinate of third corner
      (y3 + y2) / 2  //y coordinate of third corner
    );
    //Smaller triangle 3
    subTriangle
    (
      n+1, //Number of recursions for the next call increased with 1
      (x1 + x3) / 2 + (x3 - x2) / 2, //x coordinate of first corner
      (y1 + y3) / 2 + (y3 - y2) / 2, //y coordinate of first corner
      (x1 + x3) / 2 + (x1 - x2) / 2, //x coordinate of second corner
      (y1 + y3) / 2 + (y1 - y2) / 2, //y coordinate of second corner
      (x1 + x3) / 2, //x coordinate of third corner
      (y1 + y3) / 2  //y coordinate of third corner
    );
  }
}

运行程序后的效果如下:



与运算(AND)法

上面给出的方法只是绘制谢尔宾斯基三角形的众多方式之一。其中一种方式是使用 AND 运算符。将像素的 x 坐标和 y 坐标分别取整,然后对它们进行 AND 运算,若结果为 0,则绘制一种颜色的像素,否则绘制另一种颜色。这样就会出现一个谢尔宾斯基三角形!

对两个整数进行 AND 运算时,将两个整数都视为二进制数,并对每个对应的位执行 AND 运算。位运算 AND 的规则如下:

运算
结果
0 AND 0
0
0 AND 1
0
1 AND 0
0
1 AND 1
1

对整数的每一位都执行此操作,只有当结果整数的二进制表示为 0000000000000000 时,才给该像素赋予另一种颜色。

代码是一个非常简单的双重循环,遍历每个像素并检查 x & y 是否为 0,其中"&"是 C++ 中的二进制 AND 运算符。在 if 条件中使用"x & y"时,只有当 x & y 为 0 时条件才为假,在其他所有情况下都绘制白色像素,因此结果将是白色背景上的黑色谢尔宾斯基三角形。

int main(int argc, char *argv[])
{
  screen(256, 256, 0, "AND");
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    if(x & y) pset(x, y, RGB_White);
  }
  redraw();
  sleep();
  return 0;
}

当屏幕尺寸为 2 的幂次方时效果最佳,否则只能看到三角形的一部分。



随机函数法

另一种绘制谢尔宾斯基三角形近似图的完全不同的方法如下:

1)定义 3 个带坐标的点:a(ax,ay)、b(bx,by)和 c(cx,cy)。这些点将成为大外三角形的顶点。

2)再定义一个点 p(px,py),将其放置在三角形的某个顶点处(例如 px = ax,py = ay)。

3)在位置(px,py)处用铅笔画一个点

4)掷一个理论上的三面骰子(第 0 面、第 1 面和第 2 面),或者在"0"、"1"和"2"之间随机选择一个数。

5)根据掷出的数字,更新点 p 的坐标:
6)回到步骤 2,在 p 的新位置处画点,如此循环直到你不想继续为止

尽管整个过程是随机的,但经过足够多的步骤后,结果会越来越像谢尔宾斯基三角形!

这段代码并不难写,注释会解释一切:

//How much pixels should be randomly chosen and drawn
#define numSteps 10000

int main(int argc, char *argv[])
{
  //create the screen and make it white
  screen(256, 256, 0, "Sierpinski Triangle");
  cls(RGB_White);

  //the 3 corners of the outer triangle
  float ax = 10;
  float ay = h - 10;
  float bx = w - 10;
  float by = h - 10;
  float cx = w / 2;
  float cy = 10;

  //initial coordinates for the point px
  float px = ax;
  float py = ay;

  //do the process numSteps times
  for(int n = 0; n < numSteps; n++)
  {
    //draw the pixel
    pset(int(px), int(py), RGB_Black);

    //pick a random number 0, 1 or 2
    switch(abs(rand() % 3))
    {
      //depending on the number, choose another new coordinate for point p
      case 0:
        px = (px + ax) / 2.0;
        py = (py + ay) / 2.0;
        break;
      case 1:
        px = (px + bx) / 2.0;
        py = (py + by) / 2.0;
        break;
      case 2:
        px = (px + cx) / 2.0;
        py = (py + cy) / 2.0;
        break;
    }
  }

  //redraw, sleep and end the program
  redraw();
  sleep();
  return(0);
}

numSteps 的值越大,绘制的像素就越多。下面这些图像分别展示了 1000、10000 和 100000 步时的结果:

 

通过对公式稍作修改,可以得到其他形状,例如这个:



通过将部分除以 2.0 的运算改为除以 3.0 来获得。

矩形递归法

绘制谢尔宾斯基三角形的另一种方式是使用基于矩形的递归函数。递归过程如下:将每个矩形变换为 L 形:



L 形本身由 3 个矩形组成,这些矩形再次被转换为 L 形,如此循环……



如果持续进行足够长的时间,结果会越来越像谢尔宾斯基三角形。这同样可以用递归函数来实现,与前面给出的递归代码非常相似,但这次绘制的是矩形,每次都需要计算 3 个新矩形的坐标。

//the coordinates are 2 corners of the rectangle, n is the current recursion step
void drawSierpinski(int n, int x1, int y1, int x2, int y2);

//how much recursions maximum (after log_2(screenWidth) recursions the rectangles are smaller than a pixel)
#define maxRecursions 8

int main(int argc, char *argv[])
{
  //create the screen and make it white
  screen(256, 256, 0, "Sierpinski Triangle");

  //start the recursive function
  drawSierpinski(1, 0, 0, w - 1, h - 1);

  //redraw, sleep, etc...
  redraw();
  sleep();
  return(0);
}

void drawSierpinski(int n, int x1, int y1, int x2, int y2)
{
  //draw white rectangle in the upper right part, thereby making a black L
  drawRect((x1 + x2) / 2, y1, x2 - 1, (y1 + y2) / 2 - 1, RGB_White);

  //call itself 3 times again, now for the 3 new rectangles in the L shape
  if(n < maxRecursions)
  {
    drawSierpinski(n + 1, x1, y1, (x1 + x2) / 2, (y1 + y2) / 2);
    drawSierpinski(n + 1, x1, (y1 + y2) / 2, (x1 + x2) / 2, y2);
    drawSierpinski(n + 1, (x1 + x2) / 2, (y1 + y2) / 2, x2, y2);
  }
}

以下是 3、5 和 8 次递归的结果:



为了更好地理解递归,可以试试当 drawSierpinski 函数只调用自身 1 次而不是 3 次时会发生什么(另外 2 次调用被注释掉):此时若 maxRecursions 为 8,则只绘制 8 个矩形,每个都比前一个小:



现在启用两次自身调用,绘制的正方形数量增多了,但形状还不够复杂。此时右侧是 1 个大正方形,其左侧是 2 个半大正方形,再左侧是 4 个 1/4 大小的正方形,再左侧是 8 个 1/8 大小的正方形,以此类推,最左侧一列有 128 个正方形,比右侧大正方形小 128 倍。



最后,如果启用全部 3 次自身调用,就得到了谢尔宾斯基三角形,绘制的正方形数量如此之多,以至于几乎全部都是白色:



还有另一种得到谢尔宾斯基三角形的方法:使用细胞自动机,但这可能会在后续文章中介绍。

谢尔宾斯基地毯

矩形递归法

谢尔宾斯基地毯是另一种分形。手工绘制时,从一个白色正方形开始,然后在中心画一个边长为原正方形 1/3 的黑色正方形:



现在,黑色正方形周围有 8 个白色正方形。在这 8 个正方形中各画一个小 1/8 的黑色正方形,然后在 8*8=64 个新的白色正方形中再次重复:



如此无限重复下去,就得到了谢尔宾斯基地毯!

绘制时,可以使用与矩形递归法绘制谢尔宾斯基三角形非常相似的递归函数,但现在函数需要调用自身 8 次而不是 3 次,并使用不同的坐标。坐标 x1,y1-x2,y2 被划分为 9 个区域,在中心区域用 rect 绘制矩形,其余 8 个区域作为下一层递归调用 drawCarpet 的参数。

//the coordinates are 2 corners of the rectangle, n is the current recursion step
void drawCarpet(int n, float x1, float y1, float x2, float y2);

//how much recursions maximum
#define maxRecursions 6

int main(int argc, char *argv[])
{
  //create the screen and make it white, use powers of 3 for the screen size for best result
  screen(243, 243, 0, "Sierpinski Carpet");
  cls(RGB_White);

  //start the recursive function
  drawCarpet(1, 0, 0, w - 1, h - 1);

  //redraw, sleep, etc...
  redraw();
  sleep();
  return(0);
}

void drawCarpet(int n, float x1, float y1, float x2, float y2)
{
  //draw black rectangle with 1/3th the size in the center of the given coordinates
  drawRect(int((2 * x1 + x2) / 3.0), int((2 * y1 + y2) / 3.0), int((x1 + 2 * x2) / 3.0) - 1, int((y1 + 2 * y2) / 3.0) - 1, RGB_Black);

  //call itself 8 times again, now for the 8 new rectangles around the one that was just drawn
  if(n < maxRecursions)
  {
    drawCarpet(n + 1, x1         , y1         , (2 * x1 + x2) / 3.0, (2 * y1 + y2) / 3.0);
    drawCarpet(n + 1, (2 * x1 + x2) / 3.0, y1         , (x1 + 2 * x2) / 3.0, (2 * y1 + y2) / 3.0);
    drawCarpet(n + 1, (x1 + 2 * x2) / 3.0, y1         , x2         , (2 * y1 + y2) / 3.0);
    drawCarpet(n + 1,  x1        , (2 * y1 + y2) / 3.0, (2 * x1 + x2) / 3.0, (y1 + 2 * y2) / 3.0);
    drawCarpet(n + 1, (x1 + 2 * x2) / 3.0, (2 * y1 + y2) / 3.0, x2         , (y1 + 2 * y2) / 3.0);
    drawCarpet(n + 1, x1         , (y1 + 2 * y2) / 3.0, (2 * x1 + x2) / 3.0,  y2        );
    drawCarpet(n + 1, (2 * x1 + x2) / 3.0, (y1 + 2 * y2) / 3.0, (x1 + 2 * x2) / 3.0,  y2        );
    drawCarpet(n + 1, (x1 + 2 * x2) / 3.0, (y1 + 2 * y2) / 3.0, x2         ,  y2        );
  }
}

以下是 6 次递归步骤的结果:



如果想要更大的图像,可以使用 729*729 像素的分辨率和 7 次递归。

三进制法

还有一种绘制谢尔宾斯基地毯的方法,与绘制谢尔宾斯基三角形的"AND"法类似。即对每个像素进行计算,判断它是否应该着色。

不过这种方法稍微复杂一些,因为需要使用三进制(即以 3 为基数)。三进制数的每一位可以是"0"、"1"或"2"。

该方法的工作原理如下:

将像素的两个坐标都以三进制整数表示。对每一位数字,检查对应位是否不同时为"1"。如果从未出现这种情况,则该点属于地毯。

求一个数的第一位(最右位)三进制数字,将其对 3 取模(对 3 取模的结果如下:0%3=0,1%3=1,2%3=2,3%3=0,4%3=1,5%3=2,6%3=0,7%3=1,以此类推……)。
求第二位三进制数字,先将该数除以 3(整数除法,即去掉小数部分),再对 3 取模。
求第三位三进制数字,先除以 9,再对 3 取模。
求第四位,先除以 27,再对 3 取模,
以此类推……

下面的示例使用 243×243 像素的分辨率,因此只需检查像素坐标的前 5 位三进制数字,因为用 5 位三进制数字可以表示 0 到 242 的所有数字。

"if"中的条件分多行写出,对 5 位数字中的每一位检查 x 坐标和 y 坐标对应位是否不同时为 1。如果条件为真,则在 x,y 处绘制白色像素。

int main(int argc, char *argv[])
{
  screen(243, 243, 0, "Sierpinski Carpet");
  for(int y = 0; y < h; y++)
  for(int x = 0; x < w; x++)
  {
    if
    (
      //Not both the first (rightmost) digits are '1' in base 3
      !(
           (x / 1) % 3 == 1
        && (y / 1) % 3 == 1
      )

      &&

      //Not both the second digits are '1' in base 3
      !(
           (x / 3) % 3 == 1
        && (y / 3) % 3 == 1
      )

      &&

      //Not both the third digits are '1' in base 3
      !(
           (x / 9) % 3 == 1
        && (y / 9) % 3 == 1
      )

      &&

      //Not both the fourth digits are '1' in base 3
      !(
           (x / 27) % 3 == 1
        && (y / 27) % 3 == 1
      )

      &&

      //Not both the fifth digits are '1' in base 3
      !(
           (x / 81) % 3 == 1
        && (y / 81) % 3 == 1
      )
    )
    pset(x, y, RGB_White);

  }
  redraw();
  sleep();
  return(0);
}

结果与前一个程序完全相同,尽管使用了完全不同的方法:



如果想要提高分辨率,每将分辨率扩大三倍,就需要为第 6 位数字添加一个额外的条件(此时除以 243),以此类推……

检查第 5 位数字(从右数第 5 位三进制数字)的条件负责中心的大黑色正方形。
检查第 4 位数字的条件负责中心正方形周围的 8 个较小正方形,
以此类推……

如果去掉第 5 位数字的条件,中心黑色正方形就消失了,取而代之的是:



如果去掉检查第 3 位数字的条件,则所有第三阶的黑色正方形都会消失,而其他的保持不变:



因此,如果将分辨率扩大三倍,现在位于中心的正方形将变成左上角的一个正方形,而新的中心需要一个更大的黑色正方形,这就是为什么需要添加检查第 6 位数字的新条件。


最后编辑于:2004 年

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