65.9K
CodeProject 正在变化。 阅读更多。
Home

分形的理论与实践

Kenneth Haugland
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.96/5 (41投票s)

2013年9月16日

CPOL

27分钟阅读

viewsIcon

99849

downloadIcon

4339

描述了分形背后的理论以及如何使用L系统、IFS、Mandelbrot、Newton、Burning Ship、Nova和Julia集、Markus-Lyapunov分形、Kronecker乘积、DLA算法和向日葵螺旋来绘制分形图像。

介绍 

我试图理解分形的定义，但我遇到了一个难以处理的句子，因为分形被严格地定义为一个单一的句子。这种复杂性的原因可以用我给这篇文章的标题来说明。我必须承认，我不太喜欢用“理论与实践”这个短语（误）用的句子，因为其中一种用法是“理论与实践之间存在差异”，但这个句子实际上是一个毫无意义的句子，因为实际的句子“理论与实践之间存在差异”实际上是**你**提出的一个理论。因此，从数学（逻辑）角度来看它是毫无意义的，但必须指出，科学理论通常具有有限的**发散**区域或有限的用途，因为动能方程 0.5*m*v2只适用于小速度，而广义相对论

在速度达到光速约10%时就会开始影响结果，但这与最初的说法，或者说语言的不精确使用，有很大不同。

撰写本文的第二点是，通常他们只会给你展示一些非常漂亮的图片，以及如何生成它们的代码，而很少或根本不涉及其背后的数学原理，或者“我们为什么要对分形感兴趣呢？数学家和程序员难道没有比看漂亮图片更重要的事情要做吗？”

斐波那契和他的兔子

莱昂纳多·斐波那契可能是中世纪早期最著名的数学家，他最著名的贡献是斐波那契数列以及将阿拉伯数字引入欧洲。关于这些数字的故事相当复杂，他建议使用阿拉伯数字而非罗马数字的建议最初并未被接受。他的著作《算盘书》(Liber abbaci)于1202年出版，但十进制数字系统在16世纪之前并未在欧洲广泛使用，这主要是因为1以下的数字的使用。这意味着商人和其他人可以赚取一些额外的钱，就是这样。

十进制系统延迟使用的另一个原因是需要纸张，因为罗马数字系统相当节省空间（至少对于大多数数字而言）。纸张在欧洲直到14世纪才广泛使用，而在巴格达自8世纪以来就一直在使用。斐波那契带来的数字是西阿拉伯数字，或古巴数字（gubar是一个阿拉伯词，意为尘土，因为它们最初写在铺满沙子的板上）。它们不需要零，因为数字是写在列中的，因此不需要零。当他们开始在纸上书写时，他们在数字上方放置点以表示数字应在的列，据推测，数字零就是由此产生的。顺便说一下，我们当前的数字系统可以说基于三个主要思想：

  • 源自巴比伦数学的位置思想
  • 实际数字，来自印度
  • 空位数字，零，但其确切来源和起源仍不完全清楚。

斐波那契还有另一个著名的贡献，即斐波那契数列。他想要解决的问题是，在理想条件下，n个月后你会有多少只兔子？他假设

  1. 它们会生出一公一雌的后代
  2. 它们每月可以繁殖一次
  3. 它们满一个月后就可以开始繁殖
  4. 它们不会死亡

第一个月围栏里有两只兔子，但它们还不能繁殖。第二个月仍然只有一对，但现在它们可以开始繁殖了。第三个月，原始的那一对将生下两只后代，所以现在有两对。第四个月，原始的那一对又生下了一对兔子，加上上个月的那对小兔子，总共有三对。就这样，一直持续下去。斐波那契数列因此由下面的数列定义：

1, 1, 2, 3, 5, 8, 13, 21, 34, ...    

你可能想知道这与分形有什么关系。好吧，这完全与递归有关，递归是数学和编程中计算许多不同数列的第n个值的基本系统。以计算整数阶乘为例。这可以通过递归方式（函数调用自身）完成，如下所示：

    Private Function Factorial(ByVal v As Integer) As Integer
        If v > 1 Then
            Return v * Factorial(v - 1)
        Else
            Return 1
        End If
    End Function

事实证明，斐波那契数列也可以通过递归方式编程，将数列中的每个数字表示为F(n)，因此F(1) = 1，F(3) = 2 等等...

F(n) = F(n-1) + F(n-2)

并且可以使用递归关系进行编程，这意味着函数将按如下方式调用自身：

    Private Function Fibonacci(ByVal n As Integer) As Integer
        If n > 1 Then
            Return Fibonacci(n - 1) + Fibonacci(n - 2)
        Else
            Return 1
        End If
    End Function

有时递归会导致嵌套函数，这就是分形的由来。递归函数和递归系统之间存在差异。

然而，这些数字也可以与由一系列正方形构成的指数螺旋增长相关联，其中每个正方形都与相应的斐波那契数列相关。这个螺旋被称为黄金螺旋，它在Cristobal Vila的视频中 beautifully 地展示出来。遗憾的是，他犯了一个错误，因为鹦鹉螺壳不是黄金螺旋，而是一个对数螺旋。关于黄金比例的各种误解遍布网络，甚至在许多旧教科书中也存在，因此人们对此感到困惑也就不足为奇了。其中一位相当雄辩地指出这一点的人是基思·德夫林教授，他在斯坦福大学的免费讲座系列中解释了黄金比例和斐波那契数列。

可汗学院也有一个关于斐波那契黄金螺旋和松果之间（正确）联系的有趣视频。（你可能会在自然界中发现一些不符合正确斐波那契数列的例子，但这是由于基因错误突变造成的，因为排列并非最佳）。松果和黄金螺旋也可以在向日葵植物中看到，并且可以通过仅使用一个关系生成，前提是 n 趋于无穷大（称为黄金比例）：

$\phi = \lim_{n \to \infty} \frac{F(n)}{F(n-1)} = \frac{5^{\frac{1}{2}} +1}[2}$

这个数字非常有趣，主要是因为它是一个“最不理性”的数字。如果你创建任何数字的连分数，收敛速度取决于分数下方数字的大小，数字越高，收敛速度越快。在下面的通用方程中，这意味着 a 的数字越高，它将收敛得越快，换句话说，它将更快地获得更准确的小数。

$a_0+\cfrac{1}{a_1+\cfrac{1}{a_2+\cfrac{1}{a_3+\cdots}}}$

所以这个级数最糟糕的结果是当所有的a都等于1时，而这正是黄金比例发生的情况。

这使得当你想在一个磁盘上唯一地放置一个物体时，即均匀分布的点数时，这个数字变得很有用。Alexandre Devert在他的博客这里很好地解释了这一点，他在博客中展示了一张使用有理数（左）、最无理数（中）和普通无理数（右）创建的图片。

顶部中间的那个使用了黄金比例，这种点的分布方式最初是由赫尔穆特·沃格尔展示的，被称为沃格尔方法。你也可以在许多计算方法中使用它，例如声学射线追踪或全向灯光的射线追踪。

我不想直接使用Vogel的方法，因为我想改变向日葵分形中圆形的空间和大小，使其看起来更叶序（自然）。代码仍然出奇地简单：

        Dim Angle, c, R As Double
        c = (Math.Sqrt(5) + 1) / 2
        SunFlowerCanvas.Children.Clear()
        Dim x, y As Double
        Dim numberofseeds As Integer = 3000
        SunFlowerCanvas.Background = Brushes.SaddleBrown
        For i As Integer = 0 To numberofseeds
            R = i ^ c / (numberofseeds / 2)
            Angle = 2 * Math.PI * c * i
            x = R * Math.Sin(Angle) + 300
            y = R * Math.Cos(Angle) + 300
            DrawEllipse(New Point(x, y), i / (numberofseeds / 10))
            DoEvents(Me.Dispatcher)
        Next

代码生成如下所示的动画生长模式

然而，两个视频中的主题，或者向日葵生成的代码，都不是我第一次接触斐波那契数列。我当时在数学课上学习帕斯卡三角形，于是我开始乱涂乱画，发现如果我将列以45度角向右上方求和，我就得到了斐波那契数：

斐波那契 帕斯卡
三角形		                  
0 1     0 0 0 0 0 0 0 0 0 0
1 1     1 0 0 0 0 0 0 0 0 0
1 1     2 1 0 0 0 0 0 0 0 0
2 1     3 3 1 0 0 0 0 0 0 0
3 1     4 6 4 1 0 0 0 0 0 0
5 1     5 10 10 5 1 0 0 0 0 0
8 1     6 15 20 15 6 1 0 0 0 0
13 1     7 21 35 35 21 7 1 0 0 0
21 1     8 28 56 70 56 28 8 1 0 0
34 1     9 36 84 126 126 84 36 9 1 0
55 1     10 45 120 210 252 210 120 45 10 1

这在几个世纪前就已经被发现了，尽管最初的发现者可能已湮没于历史长河。

背景

分形有多种定义，其中之一是“一种粗糙或破碎的几何形状，可以细分为多个部分，每个部分（至少近似地）是整体的缩小/尺寸副本”。更正式的定义方式是使用豪斯多夫维数。如果从几何形状尺寸加倍的角度来思考它们：

物体类型           维度 (d) 方程 副本数量 (2^d = N)
Line           一维 2*x 2^1 = 2
平方           二维 2*x+2*y 2^2 = 4
立方体           3D 2*x+2*y+2*z 2^3 = 8

这意味着维数与副本数通过公式“副本数 = 2^d”相关。这个关系公式被称为豪斯多夫维数，我们可以使用这个定义，因为我们通常可以说分形中每次迭代生成的副本数，从而计算出维数 d。第一个产生不遵循此关系的曲线是朱塞佩·皮亚诺，他在1890年震惊了数学界，当时他产生了一条类似于希尔伯特曲线的曲线，其维数为 2 而非像上表所预期的 1。

谢尔宾斯基三角形每次迭代生成三个副本，并且在每个方向上长度加倍，因此得到维度：

3 = 2^d，因此 d = ln(3)/ln(2) = 1.58....

科赫雪花有点不同，因为我们从三条线组成一个三角形开始，并且对于每条直线，我们都会将其分成4个原始直线的副本，所以我们可以建立关系

N = 4

还有一种情况是，每个线段是原始线段的1/3，所以缩放因子是3而不是2，因此

d = ln(4)/ln(3) = 1.26...

然而，并非所有不同类型的分形都能通过解析方式计算其维度，但此网站提供了更多维度计算的例子。

经过上述计算和论证，现在可以将分形定义如下：

“分形是一个其豪斯多夫维数严格超过拓扑维数的集合。”

这个定义听起来仍然很难理解，但如果你还记得我介绍不精确定义时所说的，你不会觉得这很奇怪，因为数学家通常对词语的定义非常精确。怀特海和罗素在《数学原理》的前言中对此总结得很好：

“这项工作思想的抽象简洁性本身就击败了语言。语言更容易表达复杂的思想。‘鲸鱼很大’这句话代表了语言的最佳状态，简洁地表达了一个复杂的事实；而‘一是一个数’的真实分析，在语言上则导致了无法忍受的冗长。”

这句话可以被约翰·德比郡在他的著作《素数迷恋》中简单地解释为：_他们不是开玩笑。《数学原理》用了345页来定义数字“1”。_这通常意味着一个词的数学定义通常深植于复杂的定理中，因此语言背后的细节可能相当庞大。在分形的定义中，缺失的信息是拓扑维度的含义。这里变得相当棘手，因为它是一个相当复杂的推导，并且涉及集合论。然而，拓扑维度也称为覆盖维度，或者更正式地称为勒贝格覆盖维度，你可以此处阅读更详细的定义说明，但请注意，它并非一目了然。

曼德尔布罗特集合和朱利亚集合

任何关于分形的故事或总结文章都离不开曼德尔布罗特集合和朱利亚集合，原因之一是曼德尔布罗特本人创造了“分形”这个词（源自希腊词“fractus”，意为破碎或裂开）。你可以在TED演讲的这个视频中听到曼德尔布罗特本人更多地解释分形的使用和演变。通常，曼德尔布罗特集合和朱利亚集合源于一个多项式方程系统的收敛，该系统以以下形式扩展：

Z(n+1) = Zn2 + C，其中 n = 1, 2, 3, ...

这两个集合的区别在于，Mandelbrot集合是通过操纵Z值来生成图像的，而Julia集合则通过操纵C值来生成图像。你可以在这个网站上逐步了解如何在程序中精确计算它。重要的是要L注意到Mandelbrot和Julia集合是使用复数构建的，这意味着你可以将实部与x轴相关联，将虚部与y轴相关联，也可以使用四元数规则（计算复数时应用的规则）而不是汉密尔顿规则来使用这些集合的3D表示。

创建和迭代Mandelbrot集合的基本代码如下，其中假设你当前位于点(x0,y0)，并且遍历给定区域中的所有点

                Dim Z As New Complex(x0, y0)
                Dim C As New Complex(x0, y0)

                Dim MaxIterations As Integer = 200
                Dim CurrentIteration As Integer = 0
                Dim ConvergenceCondition As Double = 2.0
                Dim ConvergeneTest As Double = 0


                While CurrentIteration < MaxIterations AndAlso ConvergeneTest < ConvergenceCondition
                    Z = Complex.Pow(Z, 2) + C
                    ConvergeneTest = Complex.Abs(Z)
                    CurrentIteration += 1
                End While

对于朱利亚集，方程略有不同：

                Dim Z As New Complex(x0, y0)
                Dim C As New Complex
                C = YourComplexVariable

                Dim MaxIterations As Integer = 200
                Dim CurrentIteration As Integer = 0
                Dim ConvergenceCondition As Double = 2.0
                Dim ConvergeneTest As Double = 0


                While CurrentIteration < MaxIterations AndAlso ConvergeneTest < ConvergenceCondition
                    Z = Complex.Pow(Z, 2) + C
                    ConvergeneTest = Complex.Abs(Z)
                    CurrentIteration += 1
                End While

下面的相应图片显示了程序生成的一个朱利亚集

前面两个代码片段中给出的代码速度不是很快，因为`Complex.Abs()`调用了平方根，这本质上是泰勒级数，会占用大量时间。最好将`ConvergenceTest`设置为4，然后只执行`||Z||^2 = x^2 +(i)^2*y^2= x^2-y^2`。最终图片中的颜色取决于`ConvergenceTest`离开`ConvergenceCondition`所需的迭代次数。

这个Wolfram演示清晰地说明了Mandelbrot集合的散度测试，你会注意到Mandelbrot结构内部的所有区域无论你执行多少次迭代都不会离开它，而外部的点则需要一定数量的项才能离开收敛区域。

集合的颜色是预先生成的，并且每次迭代都会选择不同的颜色。着色取自StackOverflow上此处的一个问题的答案，但你可以用任何你想要的着色方案替换它。

牛顿分形

还有一个相当著名的集合，那就是牛顿分形，它使用牛顿-拉夫森迭代方案来生成图形，并且涉及使用函数的导数来获取C值（一个C#的简单例子可以在这里看到，但它不够通用，无法让你选择输入）。我决定也包含这个功能，但这需要一些额外的编码。首先需要定义函数，我决定从函数`Z5-1`开始：

    Private Function NewtonFractalFuntion(ByVal CC As Complex) As Complex
        Return Complex.Pow(CC, 5) - 1
    End Function

第二个问题有点难，就是获取函数的导数。对于行为良好的函数，下面给出的代码就足够了，但你也可以用精确的导数替换它，在这种情况下，精确的导数将是`5z4`（它在更高级的函数上可能会失败，但可以使用柯西-黎曼方程进行检查）。

    Private Function DerivateNewtonFractalFuntion(ByVal CC As Complex) As Complex
        Dim Delta As New Complex(10 ^ -10, 10 ^ -10)
        Return (NewtonFractalFuntion(CC + Delta) - NewtonFractalFuntion(CC)) / (Delta)
    End Function

为了获得更通用的求导函数方法，也就是非常困难的方法，可以使用柯西积分定理，它涉及到复函数的围线积分，在这种情况下，我们感兴趣的是通常所说的柯西微分公式。这太难了，以至于我还没有找到任何用于分形生成的数值方法。它们通常只是通过解析方法提供导数。

接下来是实际的主函数：

            Dim a, b, kmax As Integer
            a = 400
            b = 400
            kmax = 50

            Dim image(a, b) As Integer
            Dim x0, y0, dx, dy As Double

            dx = (xmax - xmin) / (a - 1)
            dy = (ymax - ymin) / (b - 1)

            Dim Z, Z0 As New Complex
            Dim k As Integer = 0

            For nx As Integer = 0 To a - 1
                For ny As Integer = 0 To b - 1

                    k = 0
                    'Find the current location
                    x0 = xmin + nx * dx
                    y0 = ymin + ny * dy

                    Z = New Complex(x0, y0)

                    While (k <= kmax)
                        Z0 = Z - (NewtonFractalFuntion(Z) / DerivateNewtonFractalFuntion(Z))
                        k += 1
                        If Complex.Abs(Z0 - Z) < 10 ^ -6 Then
                            k = (Z0.Phase) + 8
                            Exit While
                        End If
                        Z = Z0
                    End While
                    image(ny, nx) = k
                Next
            Next 

绘制此解决方案基本上有两种方法，一种如上述代码所示，使用复数的相位，另一种可以通过注释掉使用`Z0.Phase`的行来看到（你可以在这里看到两种方法）。相位用于为牛顿-拉夫逊技术找到的根着色，另一种方法则使用找到可接受解决方案之前的项数。还可以通过使用找到的极限组合来为图着色，以生成更具动态感的图片。这是程序生成的牛顿分形示例：

 

还应提及，牛顿分形的推广可用于生成所谓的新星分形，一旦拥有牛顿分形代码，生成这些分形的代码就相当简单。只需要进行两处修改，需要添加一个常数C和一个松弛常数R，然后只需选择这两个参数的“正确”值即可。一些漂亮的例子可以在此处和此处看到。

 

燃烧船分形

燃烧船分形（Burning Ship fractal）也被包含在内，它只是对曼德尔布罗特集合进行了微小的改动，即在对前一个Z值进行平方之前取其绝对值。下面显示了燃烧船的一个放大区域

 

马库斯-李雅普诺夫分形

马库斯-李雅普诺夫分形是一种强烈基于逻辑增长模式（费尔胡斯特模型，这个公式曾多次被重新发现）的分形。它指出增长率随时间变化（\frac{dN}{dt}）可以写成以下方程：

$\frac{dN}{dt} = \frac{r N ( K - N )}{K}$

公式中，N是任意给定时间的人口数量，K是承载能力，即能支持多少人口，r是最大增长率（这实际上是一种简化，因为这个数字最初被称为马尔萨斯参数，但我们将忽略这一点，只说它是人口增长率），t是时间。我们可以将方程N/K替换为x，然后将方程改写为：

$\frac{dx}{dt} = r x ( 1 - x )$

你可以通过解析积分或数值积分（如马库斯-李雅普诺夫分形中所示）来解这个方程。

你可以通过定义一个新变量 x 来简化逻辑增长模型，它表示当前存活的人口比例，与环境能支持（并保持存活）的总人口相比。所以，当 x = N/K 时，你得到了一个关于 x 的新微分方程。现在我们关注的是人口比例随时间的变化率。一旦 x = N/K = 1，环境就不能再支持更多的人口成员了。

迭代函数系统 (IFS)

 

迭代函数系统（简称IFS）是一系列用于生成分形的渲染变换。最简单的使用IFS绘制的图形可能是谢尔宾斯基三角形，它可以通过选择构成一个所有角度都为60度的三角形的三个点来构造。选择其中一个点，然后随机选择到三角形另外两个点之一的中点。现在你在三角形内部有了一个新点，你选择三个外部三角形点之一，并在该中点距离处绘制一个新点。如果你这样做足够多次，你最终将生成谢尔宾斯基三角形。

IFS最著名的例子是Barnsley蕨，它由4个不同的映射函数构建而成，并结合了特定变换被使用的百分比可能性。这就是蕨类植物自相似的原因，因为一个变换会随机地以不同的变换形式在另一个地方重新出现，麻省理工学院的哥德尔、埃舍尔、巴赫系列讲座的第二部分对此解释得非常好。你还会注意到，谢尔宾斯基三角形也在同一个讲座中得到解释。要查看映射函数，你可以根据它们的实现方式对不同的值进行着色。

源代码中包含了几种不同版本的蕨类植物，它们都取自这个网站。它们并非完全相同，因为源代码中的映射技术略有不同。顺便说一句，蕨类植物实际上是自然界中发现的最早的植物结构之一。

L-系统和分形

 

对植物和其他有机结构生长以及如何模拟植物细胞生长过程的研究，促使生物学家阿里斯蒂德·林登迈耶在1968年开发了一种模拟这些结构生长的新方法。为了理解他是如何思考创建这个相对简单的系统的，他大概是沿着以下思路进行的。

“在许多生物体，特别是植物的生长过程中，某些多细胞结构的规律性重复出现是显而易见的……例如，在复合叶片的情况下，一些叶片裂片（或小叶），作为成熟叶片的一部分，具有与整个叶片早期阶段相同的形状。”

起初，他用他的理论研究简单的多细胞生物的发展，但后来被应用于调查高等植物和植物器官。实际过程可以这样表述：假设你有一个起点，无论是预定的公理（或DNA中的起始基因），它构成了生长的基础。此外，它还将有一些关于如何从公理（起始条件）在离散的n个时期内扩展的规则（涉及其他DNA模式的实现）。这些基本原理和思想在20世纪初就已经被发现，但直到诺姆·乔姆斯基在1957年开始使用类似的系统来描述语言语法之前，并没有被广泛知晓或使用。然而，林登迈耶系统（简称L-系统）在替换算法上与语言不同，乔姆斯基算法是顺序的，而林登迈耶系统是并行的。正是这种差异使得L-系统在分形中的生长模式方面如此有效，而乔姆斯基算法几乎只用于语言语法。

首先需要通过公理和替换规则生成生长模式。例如，将公理设置为 F--F--F，这将生成一个三角形，前提是“-"符号代表 60 度。要生成科赫雪花，你将在每次迭代中将值 F 替换为 F+F--F+F（写为 F=F+F--F+F），并且随着迭代的进行，它会变得越来越复杂。这可以通过以下方式编码：

    Private Sub Grow()
        Dim Current, [Next] As New StringBuilder
        Current.Append(Axiom)
        [Next].Append(Axiom)
 
        Dim Found As Boolean = False
        Dim CurrentGeneration As Integer = 0
 
        While CurrentGeneration < Generations
            Current.Clear()
            Current.Append([Next])
            [Next].Clear()
 
            For i As Integer = 0 To Current.Length - 1
                Found = False
                For j As Integer = 0 To Rules.Count - 1
                    If Current(i) = Rules(j).Rule Then
                        [Next].Append(Rules(j).Pattern)
                        Found = True
                    End If
                Next
                If Not Found Then
                    [Next].Append(Current(i))
                End If
            Next
            CurrentGeneration += 1
 
        End While
        GrowString = ([Next].ToString)
    End Sub   

生成 GrowthString 后，剩下的就是实际绘制生成的分形。2D L-系统有 6 种不同的操作，这些操作通常通过乌龟在不同方向移动时绘制线条来可视化。基本的绘制命令是：

  • "F" - 绘制一条线
  • "f" - 不绘制线条地移动乌龟
  • "+" - 在当前角度上增加一个预定角度
  • "-" - 在当前角度上减去一个预定角度
  • "[" - 将当前位置保存到堆栈中
  • "]" - 移除堆栈中最顶部的点并将其设置为乌龟的当前位置

乌龟会忽略所有其他符号（乌龟会保持其当前状态），这很容易编写代码，前提是你已经创建了迭代树并替换了重写系统中使用的公理。

    Public Function GetDrawingLines() As List(Of Line)
        'Make the tree from axioms and rules
        Grow()
 
        'Replace temporary axioms with new ones
        ReplaceExchangeAxiomsInGrowString()
 
        'Set start position for the turtle
        Dim CurrentTurtle As New LindemayerCurrentPoint
        CurrentTurtle.Posision = New Point(200, 200)
        CurrentTurtle.CurrentAngle = -90
        CurrentTurtle.LineLength = LineLength
 
        'Create a temporary next point
        Dim NextTurtle As New LindemayerCurrentPoint
 
        'Create a stack to allow branching
        Dim PointList As New Stack(Of LindemayerCurrentPoint)
 
        'The drawn lines
        Dim result As New List(Of Line)
 
        For i As Integer = 0 To GrowString.Length - 1
            If GrowString(i) = "F"c Then
                'Calculate the next point
                NextTurtle = GetNextPoint(CurrentTurtle)
 
                'Draw the curve
                result.Add(GetLine(CurrentTurtle, NextTurtle))
 
                'Move the turtle
                CurrentTurtle = NextTurtle.Clone
            ElseIf GrowString(i) = "f"c Then
                'Move the turtle without drawing a line
                CurrentTurtle = GetNextPoint(CurrentTurtle)
            ElseIf GrowString(i) = "+"c Then
                CurrentTurtle.CurrentAngle += Angle
            ElseIf GrowString(i) = "-"c Then
                CurrentTurtle.CurrentAngle -= Angle
            ElseIf GrowString(i) = "["c Then
                PointList.Push(CurrentTurtle.Clone)
            ElseIf GrowString(i) = "]"c Then
                CurrentTurtle = PointList.Pop
            End If
        Next
 
        Return (result)
    End Function

网上有很多这个算法的实现，但在我看来（恕我直言）它们过于复杂，可能会让想实现和理解它们的人感到困惑。这些简单的规则能生成如此复杂的结构，甚至是相当复杂的树，这真是令人惊叹。

还有另一种生成树分形的方法，即使用递归算法，你可以通过类似于以下命令来绘制它

Private sub DrawBranch(Angle, MaxIterations)
 
... 
 
... 
 
DrawBranch(Angle+IncrementAngle,MaxIterations-1)
DrawBranch(Angle-IncrementAngle,MaxIterations-1) 
 
End Sub

穆罕默德·礼萨·霍斯拉维在他的文章DotNet Real Tree中正是采用了这种方法，他允许你设置许多不同的参数来生成各种不同的树。我还有一点没有提到，那就是自然树在分叉时如何减小其横截面。莱昂纳多·达·芬奇提出的公式被认为是自然真理，其中d1是分叉前的横截面积，d2和d3分别是两个新分叉的横截面积，得到以下关系：

d1 = Math.Sqrt(d22 + d32) 

提出此关系是为了让树液能够到达树顶。如果两个分支，d2和d3相等，则该方程表明每个分支的粗细将为d1/Math.Sqrt(2)。

由于L-系统可以生成自然界中种类繁多的天然生物和植物，人们开始思考，当迭代可能性与自然选择（即像人类一样做出选择）结合时会发生什么。理查德·道金斯提出了一个简单的方案来实现这一点，萨沙·巴伯的文章AI：道金斯生物形态/和其他进化生物就展示了这种算法的一个例子。

我还要提到，你也可以在三维空间中实现L-系统，这将使字符串生成器保持不变，但会增加两个额外的角度。一个例子是三维希尔伯特曲线，可以在此处看到。我没有实现它，因为它会占用大量编码工作，而且不会显示出L-系统的任何进一步深入的见解。然而，有一本由林登迈耶等人撰写的书，可以免费下载，其中出现了很多三维结构，网址是此处。

为了说明不同形状的差异有多小，我将公理设置为：

-F-F-F-F-F-F

并且每次迭代的交换规则为：

F=F+F-F-F+F

给定4次迭代，角度为60度，我们得到一个类似科赫雪花的图形：

将角度改为90度，则得到以下图案

再次，将角度设置为120度，我们得到一个谢尔宾斯基三角形的版本

这三种图案都非常独特，但唯一改变的只是角度。然而，它们确实有重叠，这意味着它们多次绘制同一条线，但这确实表明创建截然不同的形状只需很小的变化。

克罗内克积与分形

 

要生成克罗内克分形，真正需要做的就是创建维基百科文章中定义的两个矩阵之间的克罗内克积算法，这可以通过以下方式完成：

    ''' <summary>
    ''' From "Mathematicl tools in comuter graphics with C# implementations", by Hardy & Steeb, World Scientific
    ''' </summary>
    ''' <param name="M" />Matrix x
    ''' <param name="N" />MAtrix B
    ''' <returns />
    ''' <remarks />
    Private Function Kronecker(ByVal M(,) As Double, ByVal N(,) As Double) As Double(,)
        Dim R(M.GetLength(0) * N.GetLength(0), M.GetLength(1) * N.GetLength(1)) As Double
 
        For i As Integer = 0 To M.GetLength(0) - 1
            For j As Integer = 0 To M.GetLength(1) - 1
                For k As Integer = 0 To N.GetLength(0) - 1
                    For l As Integer = 0 To N.GetLength(1) - 1
                        R(i * N.GetLength(0) + k, j * N.GetLength(1) + l) = M(i, j) * N(k, l)
                    Next
                Next
            Next
        Next
        Return R
    End Function 

现在你所需要做的就是用0到1之间的值填充两个矩阵，并将它们相乘。要绘制图像，将0设为白色，1设为黑色，仅此而已。

还有其他可能性可以实现绘图并允许超出0到1范围的数字，但那样你需要相应地更改代码。

DLA - 扩散限制聚合

 

DLA是最容易理解的分形之一，生成它的方法之一是使用以下类比：一个醉汉四处寻找客栈，他找到一个后就待在那里。他从一个任意点开始，随机移动到他的九个邻居之一，这个过程重复进行，直到他到达一个有客栈的邻居。

生成醉汉起始点的代码如下

''' <summary>
''' Select a starting point on one of the 4 edges, and select one pixel
''' </summary>
''' <param name="width" />
''' <param name="height" />
''' <returns>
''' <remarks>
Private Function GenerateStartPoint(ByVal width As Integer, ByVal height As Integer) As Point
    Dim st, x, y As Integer
    st = Rand.Next(1, 5)
    Select Case st
        Case 1
            x = 0
            y = Rand.Next(0, height - 1)
        Case 2
            x = Rand.Next(0, width - 1)
            y = 0
        Case 3
            x = width - 1
            y = Rand.Next(0, height - 1)
        Case 4
            x = Rand.Next(0, width - 1)
            y = height - 1
    End Select
    Return New Point(x, y)
End Function</remarks>

要移动，需要选择一个紧邻当前点的任意点。我选择了模拟自身重叠的矩形，就像球体结构一样。如果它想向右移动，并且当前位于最右边的点，它就会移动到左侧的起点等等。代码如下所示

''' <summary>
''' This will assume that the rectangle is connected in a sphere, when it reaches maximum it will start on 0 again for both directions
''' </summary>
''' <param name="CurrentPoint" />The current point location
''' <param name="height" />The height of the surounding rectangle
''' <param name="width" />The height of the surounding rectangle
''' <returns>A random move to one of its neighbor points</returns>
''' <remarks>
Private Function MoveToNeighbor(ByVal CurrentPoint As Point, ByVal height As Integer, ByVal width As Integer) As Point
    Dim x, y, a, b As Integer
    x = CurrentPoint.X
    y = CurrentPoint.Y
    a = width
    b = height

    Dim seed As Integer = Rand.Next(1, 9)
    ' If one assumes the starting point x;y is 10;10 youll see the
    ' that it will move to the following neighboring point
    Select Case seed
        Case 1
            'Move to point x;y = 9;10
            If x = 0 Then
                x = a - 1
            Else
                x -= 1
            End If
        Case 2
            'Move to point x;y = 9;9
            If x = 0 Then
                x = a - 1
            Else
                x -= 1
            End If

            If y = 0 Then
                y = b - 1
            Else
                y -= 1
            End If

        Case 3
            'Move to point x;y = 10,9
            If y = 0 Then
                y = b - 1
            Else
                y -= 1
            End If

        Case 4
            'Move to point x;y = 11;10
            If x = a - 1 Then
                x = 0
            Else
                x += 1
            End If
        Case 5
            'Move to point x;y = 11;11
            If x = a - 1 Then
                x = 0
            Else
                x += 1
            End If

            If y = b - 1 Then
                y = 0
            Else
                y += 1
            End If
        Case 6
            'Move to point x;y = 10;11
            If y = b - 1 Then
                y = 0
            Else
                y += 1
            End If
        Case 7
            'Move to point x;y = 11;9
            If x = a - 1 Then
                x = 0
            Else
                x += 1
            End If

            If y = 0 Then
                y = b - 1
            Else
                y -= 1
            End If
        Case 8
            'Move to point x;y = 9;11
            If x = 0 Then
                x = a - 1
            Else
                x -= 1
            End If

            If y = b - 1 Then
                y = 0
            Else
                y += 1
            End If
    End Select

    Return New Point(x, y)
End Function</remarks>

这个算法可以用多种方式修改，但我程序中的示例是最简单的形式，并且可以进行许多增强，如这个示例（Java代码）所示。它还可以模拟自然界中的许多现象。

着色方案

我在本文中使用了两种不同的着色方案，一种生成与周围着色不同的独特颜色。我现在还添加了一种平滑着色方案，它允许您设置颜色的渐变点，并且它只是从一种颜色线性迭代到下一种颜色，生成平滑的图形。这只需使用以下代码即可完成

Public Sub GetGradients(ByVal starts As Color, ByVal ends As Color, ByVal steps As Integer)
    Dim stepA As Integer = Math.Truncate((CInt(ends.A) - CInt(starts.A)) / (CInt(steps) - 1))
    Dim stepR As Integer = Math.Truncate((CInt(ends.R) - CInt(starts.R)) / (CInt(steps) - 1))
    Dim stepG As Integer = Math.Truncate((CInt(ends.G) - CInt(starts.G)) / (CInt(steps) - 1))
    Dim stepB As Integer = Math.Truncate((CInt(ends.B) - CInt(starts.B)) / (CInt(steps) - 1))

    For i As Integer = 0 To steps - 1

        ColorList.Add(Color.FromArgb(CInt(starts.A) + CInt(stepA) * CInt(i),
                        CInt(starts.R) + CInt(stepR) * CInt(i),
                        CInt(starts.G) + CInt(stepG) * CInt(i),
                        CInt(starts.B) + CInt(stepB) * CInt(i)))
    Next

End Sub

然而，这有一些缺点，因为着色与Mandelbrot集中给定的迭代次数不匹配，因此您必须注意列表中有足够的颜色来覆盖不同分形中所有给定的迭代次数。

只需在此处更改渐变停止点中的颜色即可替换颜色。

Sub New()

    GradientStop.Add(Brushes.Red)
    GradientStop.Add(Brushes.Yellow)
    GradientStop.Add(Brushes.Green)
    GradientStop.Add(Brushes.Magenta)
    GradientStop.Add(Brushes.Blue)
    GradientStop.Add(Brushes.Purple)
    GradientStop.Add(Brushes.Gold)
    For i As Integer = 0 To GradientStop.Count - 2
        GetGradients(GradientStop(i).Color, GradientStop(i + 1).Color, 35)
    Next

End Sub

历史

正如你可能已经看到的，分形之所以如此有趣，是因为它能够模拟植物、自然界中的生长模式，这意味着从一些简单的公理中可以生成相当复杂的形状。

首次发布：2013年9月16日。

参考文献

我将从这里最重要的参考文献开始

  • 麻省理工学院“哥德尔、埃舍尔、巴赫”系列讲座 
  • Alexandre Hardy & Willi-Hans Steeb 撰写的《带有 C# 实现的计算机图形学数学工具》，World Scientific
  • 《带有C#实现的计算机图形学数学工具》的源代码可以从此处下载。
  • 曼弗雷德·施罗德的《分形、混沌、幂律》。这是一本非常易读的书，其中包含大量关于不同类型分形的信息。

其他有用的在线资源包括：

© . All rights reserved.