基于Python turtle库绘制Koch Snowflake

前言

昨日下午一友人问我怎么用C语言绘制Koch Snowflake. 而我从未接触过C语言的图形库, 心想着能搞定吗? 再细细一问, 原来是大一新生C语言的作业, 只需要使用控制台打印即可. 费了几番周折, 总算是写了一个堪堪能用的版本出来. 但是这个Koch曲线倒是引起了我的兴趣. 然而C语言确实是许久未用了, 难免有些生疏, 何况其绘图捉急, 因此选择使用python 的turtle来绘制图像.

何为Koch曲线?

The Koch snowflake (also known as the Koch curve, Koch star, or Koch island) is a fractal curve and one of the earliest fractals) to have been described. It is based on the Koch curve, which appeared in a 1904 paper titled "On a Continuous Curve Without Tangents, Constructible from Elementary Geometry"by the Swedish mathematician Helge von Koch.

引用自:Koch snowflake - Wikipedia

显然, Koch雪花是分形曲线的一种, 它有着有限的面积, 却拥有无限的周长.

下列图像展示了Koch曲线的变化过程:

n1

n2

n3

n4

让我们具体来看一下曲线的生成过程. 绘制伊始, 我们有一条直线段:

line1

取中间三分之一的线段, 向外做一个正三角形.

line2

以此类推, 对每一个线段都取中间长13\frac13的线段, 向外做正三角形. 最终, 我们就得到了一条koch曲线.

显然, 这类问题一般首先就能想到递归法求解. 递归虽然形式简单易懂, 但终究差点意思.

观察图像, 发现Koch曲线实际上只由3种线段组成, 即线段A,B,CA,B,C:

line3在递归法中, 曲线的走向由迭代的深度决定,当递归到最底层的时候,就直接画线段, 不再向外做三角形了.

如果考虑turtle本身的特性, 采用一笔画的方式绘制. 我们可以定义六种基操作, 分别为OP0,OP1,,OP5OP_0,OP_1,\cdots,OP_5向六个方向绘图.

操作方向

定义了基本操作后, 绘制一条Koch曲线就可以表述为:

OP0,OP1,OP5,OP0OP_0,OP_1,OP_5,OP_0

使用矩阵记法, 简写为:

[0150]\left[ \begin{matrix} 0&1&5&0 \end{matrix} \right]

我们将每一个数字都称作操作数, 代表曲线前进的方向和一定的长度.

再拓展为Koch Snowflake, 可以得到一个3×43\times 4的矩阵:

[015045342312]\left[ \begin{matrix} 0&1&5&0\\ 4&5&3&4\\ 2&3&1&2\\ \end{matrix} \right]

它的效果就是:

n2

有了"矩阵"工具, 研究起问题就方便多了, 但是, 怎么才能得到3×n3\times n矩阵呢?

依然观察某一段曲线作为特殊例子, 如果给定一个前进方向, 由于Koch曲线的性质, 曲线会在13\frac13处向外旋转60°, 再向内旋转120°, 最后再23\frac23处旋转回原先的前进方向.

旋转

如果以xx方向作为前进的方向, 那么第一次旋转相对于xx逆时针旋转60°, 第二次旋转相对xx顺时针旋转60°, 最后一段仍然按照xx方向前进. 所以, 对于nn阶Koch Snowflake图形, n+1n+1阶图形就是将每一个线段都向外做正三角形, 也就是将xx方向拆分成上述的四个小部分. 用矩阵的方法来表示,我们可以得到这样一个生成公式:

[x][xx+1x1x]\left[ \begin{matrix} x \end{matrix} \right] \rightarrow \left[ \begin{matrix} x &x+1&x-1&x \end{matrix} \right]

同时, 角度运算是一种类似于时钟的模运算, 时针后退一小时, 相当于时针前十一小时. 同理, 矩阵操作可修正为:

[i][i%6(i+1)%6(i+5)%6i]\left[ \begin{matrix} i \end{matrix} \right] \rightarrow \left[ \begin{matrix} i \% 6 & (i+1)\%6& (i+5) \%6& i \end{matrix} \right]

ii050\sim5的操作数.

举个例子, 一阶Koch曲线生成二阶Koch曲线, 矩阵操作为:

[042][015045342312]\left[ \begin{matrix} 0\\ 4\\ 2\\ \end{matrix} \right] \rightarrow \left[ \begin{matrix} 0&1&5&0\\ 4&5&3&4\\ 2&3&1&2\\ \end{matrix} \right]

以此类推, 如果生成三阶Koch曲线, 那么矩阵规模就已经达到了3×163\times16. nn阶Koch曲线的矩阵大小可达到3×4n3\times 4^n, 呈现指数增长. 因此, 我们绘制的时候八九阶就差不多了.

但是光这样达不到我们的要求, 这每次都需要从1阶开始生成, 和递归没什么区别. 再次观察, 发现, 一阶矩阵其实早已经蕴含在二阶矩阵的第一列. 这并非偶然. 分形曲线的局部可以无限放大, 并且其形状不会发生改变. 这也就意味着n+1n+1阶矩阵是nn阶矩阵的延伸. 这一性质保证了迭代公式是一定存在的.

矩阵生

注意到, 第ii个位置, 如果ii能够被4整除, 那么其操作数继承于上一个矩阵, 也就是自身i//4i//4的位置("//"代表整除) 而对于不能够整除的位置, 其操作数取决于前一个能被4整除的位置. 而n+1n+1阶矩阵又包含了nn阶矩阵, 可以通过已生成的序列得到未生成的序列. 这样 就得到了一个迭代公式:

xi={xi//4i%4=0(xi1+1)%6i%4=1(xi2+5)%6i%4=2xi3i%4=3x_i= \begin{cases} x_{i//4} & i\%4 = 0\\ (x_{i-1}+1)\% 6& i\%4 = 1\\ (x_{i-2}+5)\% 6& i\%4 = 2\\ x_{i-3}& i\%4 = 3\\ \end{cases}

xix_i表示为矩阵的某一列. 将生上述公式转写成代码:

import turtle
import numpy as np
import math as math


def paint(angle, length):
    turtle.seth(angle * 60)
    turtle.fd(length)


def init(size, length):
    a = (pow(3, size - 1)) * length
    b = 2 / 3 * math.sqrt(3) * a
    turtle.setup(a + 50, b + 50)
    if size >= 9:
        turtle.tracer(False)
        # 关闭轨迹会导致曲线最后一段生成出现问题, 不太清楚原因
        # 但是高阶反正已经看不见了, 所以关闭加快速度
    else:
        turtle.speed(0)
        turtle.delay(0)
    turtle.penup()
    turtle.goto(-a / 2, a / 2 / math.sqrt(3))
    turtle.pendown()


if __name__ == '__main__':
    size = 3  # n阶Koch Snowflake
    if size == 0:
        print("N 必须大于零")
        exit()
    length = 2000.0 / math.pow(3, size)
    init(size, length)

    S = np.zeros((3, pow(4, size - 1)))
    S[0][0] = 0
    S[1][0] = 4
    S[2][0] = 2
    for i in range(1, S.shape[1]):
        if i % 4 == 0:
            S[:, i] = S[:, i // 4] % 6
        elif i % 4 == 1:
            S[:, i] = (S[:, i - 1] + 1) % 6
        elif i % 4 == 2:
            S[:, i] = (S[:, i - 2] - 1) % 6
        else:
            S[:, i] = S[:, i - 3]
    for i in range(3):
        for j in range(S.shape[1]):
            paint(S[i][j], length)
    turtle.hideturtle()
    turtle.done()  # 防止Pycharm关闭turtle

运行程序, 就可以得到漂亮的雪花图样.

雪花

最后, 如果想旋转/拉升图像, 只需要对矩阵直接进行变化操作即可. 当然 如果要拉伸图像还得定义另外一个长度矩阵 用于衡量每一小线段的长度. 所以, 使用操作矩阵表示图像能比递归更加容易理解分形的生成过程, 对图像的各种操作也更方便.