Click here to Skip to main content
13,007,045 members (62,156 online)
Click here to Skip to main content
Add your own
alternative version


74 bookmarked
Posted 15 Sep 2013

Fractals in theory and practice

, 28 Sep 2013
Rate this:
Please Sign up or sign in to vote.
Describes the theory behind fractals and how to draw fractal images, using L-system, IFS, Mandelbrot, Newton, Burning ship, Nova and Julia set, Markus-Lyapunov fractals, Kronecker products, DLA algorithm and the Sunflower spiral.


I was trying to understand the definitions of fractals and I got an impossible sentence to deal with, as a fractal is religiously defined by a single sentence. The reason for the complexity could be illustrated by the title I have given this article. I must admit that I don't really like sentences that is formulated with the (mis)use of the phrase "theory vs practice", as one of the usages is "There is a difference between theory and practice", but this sentence is in reality a meaningless sentence, as the actual sentence "There is a difference between theory and practice" is actually a theory you have formulated. Thus in a mathematical (logical) way it is meaningless, however it must be said that scientific theories generally have a limited area of divergence, or limited usage, as the equation of kinetic energy 0.5*m*vis only valid for small speeds, as the theory of general relativity

will start to impact the results for speeds roughly 10% of the speed of light, but this is very different from the original statement, or rather an imprecise usage of language.  

This is the second point of writing this article is that normally they usually show you some really nice pictures, and a code of how to generate them, and little or nothing of the mathematics behind them, or "Why are we interested in fractals at all? Surely mathematicians and programmers have better things to do than to look at pretty pictures?"    

Fibonacci and his rabbits

Leonardo from Fibonacci is perhaps the best-known mathematician in the early medieval period, and is best known for the Fibonacci series and bringing Arabic numerals to Europe. The story behind the numerals are a quite complicated one, and his suggestion of using Arabic numerals over roman numbers was not accepted at first. His book, Liber abbaci was published in 1202, but decimal number systems were not in wide use in Europe before the 16th century, but this was mainly to the usage of numbers below 1. This meant that merchants and others could make some additional money, and there you go.   

Another reason for the delayed usage of decimal system was the need for paper, as number system as the roman is quite space saving (at least for most numbers). Paper  was not in wide use before the 14th century in Europe, while it had been used since the 8th century in Baghdad. The numbers that Fibonacci brought with him was the west Arabic numerals, or gubar numerals (gubar is an Arabic word meaning dust, as they were first written on sand covered boards). They had no need for a zero, as the numbers were written in columns, and they thereby didn't need the zero. When they started writing on paper, they put dots over the number to signify the column the number should be in, and it is assumed that the number zero came from this. For the record, our current number system could be said to be based on three main ideas: 

  • The position idea that stems from Babylonian mathematics 
  • The actual numbers, that came from India 
  • The number for the empty space, zero, but the actual source and origin of this is still not completely clear. 

Fibonacci had more rabbits in his hat to pull out, namely his other famous contribution, the Fibonacci series. The problem he wanted to solve was, under ideal conditions, how many rabbits would you have after n months? He assumes that:

  1. They will produce one male and one female offspring
  2. They can reproduce once a month
  3. They can start to reproduce once they are more than one month old
  4. They don't die 

The first month there are two rabbits in the pen, however they can't reproduce yet. The second month there are still the one pair, but now they can start to reproduce. The third month, the original pair will have two offsprings, so now there are two pairs. The fourth month the original pair has yet another pair of rabbits, and there are  the young pair from the previous month, a total of three pairs. And so it continues, on and on. The Fibonacci series, is thus defined by the number series below:   

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

You might want to know what this has to do with fractals at this point. Well, it's all about recursion, which is a fundamental system in which we in mathematics and programming can calculate the nth values of many different series. Take the example of calculation the Factorial of an integer. This could be done in a recursive way (the function calls itself) as shown below:  

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

As it turns out the Fibonacci series could also be programmed in a recursive way, by letting each number in the series detonated by F(n), thus F(1) = 1 F(3) = 2 etc... 

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

And can be programmed using the recursive relation, meaning that the function will call itself, as follows: 

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

Sometimes recursion can lead to nested functions, and this is were Fractals comes from. There is a difference between recursive functions and recursive systems.

However, the numbers can also be related to an exponential spiral growth constructed by a series of squares where each square is related to the corresponding Fibonacci sequence. This spiral is called the Golden spiral, and it is beautifully shown in the video by Cristobal Vila. He sadly got something wrong, as the Nautilus shell is not a golden spiral, but rather a logarithmic spiral. There are all sorts of misconceptions about the Golden ratio all over the web, and indeed in many older textbooks, so it's no wonder that people are confused about it. One of the people that has pointed this out rather eloquently is Professor Keith Devlin who explains the golden ratio and Fibonacci series in his free lecture series at Stanford University. 

The Kahn Academy also have a fun video of the (correct) connections between the Fibonacci golden spiral and pine cones. (You might find some examples in nature that does not exhibit the correct Fibonacci numbers, but this is due to a wrong mutation of a gene as the placing is not optimal). The cones and the golden spirals, can also be seen in the sunflower plant, and could be generated by the use of only one relation, given that n goes to infinity (called the Golden ratio): 

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

This number is very interesting, mainly because it is the number that is the "most irrasional". If you create a continued fraction of any number, the speed of conversion is decidd by how high the number under the fractions are, the hight the number the faster the convergence happens. In the general equation below, this will mean that the higher the numbers a are, the faster it will converge, in other words, it will have more accurate deciumals faster with higher numbers.


So the worst outcome for this series is when all of the a's are equal to 1, and that is exactly what happens at the Golden ratio.

This makes the number desirable to use if you want a unique placemnt of an object on a disk, read; an evenly distributed number of points. This is very well explained by Alexandre Devert in his blog here, were he shows a picture that is created by using rational numbers(left), most irrational numbers(middle), and normal irrational numbers(right):

The top middle one uses the golden ratio, and this spread of points was first shown by Helmut Vogel, and is called Vogel's method. You could also use this in lots of computational methods, as in acoustic raytracing or raytracing of light from a omnidirectional lamp.

I didnt want to use Vogel's method directly, as I wanted it to change the space and size of the circles in the sunflower fractal, to look more Phyllotaxis (natural). The code is still surprisingly simple though:  

Dim Angle, c, R As Double
c = (Math.Sqrt(5) + 1) / 2
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))

The code generates an animated groth pattern seen below:

However none of the subjects in the two videos, or the sunflower generated code, were my first encounter with the Fibonacci series. I was sitting in my math class learning about Pascal's triangle, so I started doodling about, and found that if I summarized the columns in 45 degrees upwards to the right i got the Fibonacci numbers:   

01  0000000000
11  1000000000
11  2100000000
21  3310000000
31  4641000000
51  510105100000
81  6152015610000
131  72135352171000
211  828567056288100
341  936841261268436910
551  104512021025221012045101

This had indeed been discovered many centuries before although the original discoverer is probably lost in time. 


Fractals have severals definitions where one of them are "a rough or fragmented geometric shape that can be subdivided in parts, each of which is (at least approximately) a reduced/size copy of the whole". A more formal way of defining them is by the use of Hausdorff dimensions. If one thinks of them in the way of doubling of the size of several different geometric shapes: 

Object type     Dimension (d)EquationNumber of copys (2^d = N)
Line     1D2*x2^1 = 2
Square     2D2*x+2*y2^2 = 4
Qube     3D2*x+2*y+2*z2^3 = 8

This implies that the number for dimensions is related to the number of copies by the formula Number of copies = 2^d. This relation the formula is called the Hausdorff dimension, and we can use the definition, as we in general could say the number of copies each iteration that is generated in a fractal, and thereby calculate the number d dimension. The first one to produce a curve that did not follow this relation was Guiseppe Peano, who shocked the mathematical community in 1890, when he produced a curve similar to the Hilbert Curve, that had the dimension 2 and not 1 as one would expect from the table above.

The Sierpinski triangle generates three copies per iteration, and it doubles the length in each direction, and thereby on gets the dimension:  

3 = 2^d which gives d = ln(3)/ln(2) = 1.58.... 

A Koch snowflake is a little different, as we start of by three lines that forms a triangle, and for each straight line we split up we get 4 copies of the original straight line, so we can set up the relation:

N = 4 

There is also the case that each line segments is 1/3 of the original line so the scaling factor will 3 and not 2, so:

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

However it is not always possible to calculate the dimension analytically for all different types of fractals, but more examples of dimensions calculation is given on this site.

The calculation, and with argumentation above, one can now define the fractal as the following:

"A fractal is a set for which the Hausdorff dimension strictly exceeds the topological definition."

This definition still sounds quite difficult to grasp, but if you remember my introduction of imprecise definitions, you will not find this very surprising, as mathematicians define words, in general, very precise. This is quite well summarized by Whitehead and Russel in the preference to Principia Mathematica:

"The very abstract simplicity of the ideas of this work defeats language  Language can represent complex ideas more easily. The preposition "a whale is big" represents language at its best, giving terse expression to a complicated fact; while the true analysis of "one is a number" leads, in language, to an intolerable prolixity."

This statement could be explained quite simply as John Derbyshire does in his book Prime Obsession: They weren't kidding. Principia Mathematica takes 345 pages to define the number "1".  This usually means that a  mathematical definition of words usually is deeply rooted in complex theorems, and thus the specifics behind the language can be quite wast. In the definition of the fractal, the missing piece of information is what one means with a topological dimension.  Here it gets rather tricky, as it is quite a complicated derivation, and it involves set theory. However a topological dimension is also called Covering dimension or more formally correct Lebesgue covering dimension, and you can read a more detailed account of the definitions here, but be warned it is not straight forward. 

Mandelbrot and Julia set  

There is no story or summary articles about fractals without the Mandelbrot and Julia set, for one thing it is actually Mandelbrot that coined the phrase Fractal (stemming from greek word "fractus" that means broken or fractured. You could hear Mandelbrot himself explain more about the usage and evolutions of fractals at this video from the TED talks. In general the Mandelbrot and Julia set stems from the convergence of a polynomial system of equations that is expanded on the from: 

Z(n+1) = Zn2 + C  where n = 1, 2, 3, ...  

The difference between the two sets are that the Mandelbrot set is generated by manipulating the Z value in order to generate the picture, while Julia sets manipulates the C value to generate a picture. You can go through this step by step process on this web site to see exactly how it is calculated in the program. It is important to notice that the Mandelbrot and Julia sets are constructed using complex numbers, meaning that you could relate the real part to the x axis and the imaginary part to the y axis, there are also the possibility of using a 3D of the sets, using Hamiltonian rules (the rules that apply when calculating with complex numbers) using quaternion rules instead.

The basic code for creating and iteration through a Mandelbrot set is given below, where one assumse that you are currently standing in point (x0,y0) and that all the points in the given area are looped through:

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

For a Julia set the equation becomes slightly different: 

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

The corresponding picture below shows one of the Julia sets generated by the program:

The code that is given in the two previous snippets isn't very fast, as the Complex.Abs() calls square root, which is in nature a Taylor series, and take up lots of time. It is better to set the ConvergenceTest to 4, and just do ||Z||^2 = x^2 +(i)^2*y^2= x^2-y^2 instead. The colors in the final picture is decided on how many iterations is needed for the ConvergenceTest to leave the ConvergenceCondition

The test of divergence is clearly illustrated by this Wolfram demonstration of a Mandelbrot set, and you'll notice that all the area the is inside the Mandelbrot structure wont leave it no matter how many iterations you perform, while points on the outside use a number of terms before it leaves the area of convergence. 

The colors of the sets are is generated in advance, and for each number of iteration it performs it pics a different color. The coloring is taken from an answer to a StackOverflow question here, however you can replace this with any coloring scheme you desire. 

Newton fractal

There are also another set which is quite well known, that is the Newton Fractal, which uses the Newton-Raphson iteration scheme to generate a drawing, and it involves using the derivative of a function in order to get the value of C (a simple example in C# could be seen here, however it is not quite general enough to allow you to choose the input). I decided to include this functionality as well, however this require some additional coding as well. First off the function needs to be defined, and I decided to start with the function Z<sup>5</sup>-1

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

The second problem is a bit harder, to get the derivative of the function. For a well behaved function the code given below would work well enough, but you can also replace it with the exact, which in this given case will be 5z<sup>4</sup> (it can fail on more advanced functions, however it can be checked using the Cauchy-Riemann equation). 

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

For a more general method of getting the derivative function, the really hard way, one can instead use Cauchy's Integral theorem, which involves the contour integral of a complex function, and in this case we are interested in what's often called Cauchy's differentiation formula. This is so hard that I have yet to find any numerical methods that uses this for fractal generation though. They usually just provide the derivative by analytical means instead. 

Then it is the actual main 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

There are basically two ways of plotting this solution, one as is shown in the code above using the phase of the complex number, and the other one can be seen by commenting out the line that uses Z0.Phase (You can see both here). The phase is used to color what root the Newton-Raphson technique will find, and the other will use the number of terms before an acceptable solution is found. There is also the possibility of coloring the plot by using a combination of limits that are found, to generate more dynamic looking pictures. This is the example of a Newton fractal generated by the program:


It should also be mentioned that a generalization of the Newton fractal can be used to generate so-called Nova fractals, these are pretty simple code to generate once you have the Newton fractal code. Only two modifications are needed, you need to add a constant C, and a relaxation constant called R, and then is just a matter of choosing the "right" values for these two parameters. Some pretty examples can be seen here and here.


Burning ship fractal

The Burning ship fractal is also included, which is just a slight change from the Mandelbrot set, in that the absolute value is taken before squaring the previous Z value.  A zoomed in section of the Burning Ship is shown below:


Markus-Lyapunov fractal

The Markus-Lyapunov fractal is a fractal that is strongly based on the logistic growth pattern (Verhulst model and this formula has been rediscoverd several times). This states that the rate of growth over time (<span class="math">\frac{dN}{dt}</span>) can be written as the equation:

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

In the formula N is the number of people at any given time, K is the carry capasity i.a. how many people it can support, r is the maximum growth rate (well this is actually a simplification, as the number was originally called Malthusian Parameter, but we will ignore this and say that it is only the population growth rate) and t is the time. We can simplegy the equation by replacing N/K with x and rewrite the equation as this:

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

You can solve this equation, either by analytical integration or, as is done in then Markus-Lyapunov fractal, by numerical integration.

You can simplify the logistic growth model by defining a new variable x to represent the portion of the population that’s alive, compared to the total population that the environment could support (and keep alive). So with x = N/K, you get a new differential equation in terms of x. Now we are looking at the rate of change of the population fraction over time. Once x = N/K = 1, the environment can’t support any more members in the population:

Iterated Function System (IFS)    


Iterated Function System (IFS for short) is a series of rendering transforms that is used to generate a fractal. The simplest from to draw using IFS is probably the Sierpinski triangle, and it can be constructed by choosing the tree points that makes up a triangle where all the angels are 60 degrees. One of the points is picked, and than, you will at random pick the halfway point to one of the other two points in the triangle. You now have a new point inside the triangle, and you pick one of the tree outer triangle points, and draw a new point at that halfway distance. If you do this enough times, you will eventually generate the Sierpinski triangle. 

The most commonly know example of IFS is the Barnsley fern that is constructed using 4 different mapping functions, together with a percentage likelihood that a particular transform is used. This is what makes the fern self similar, as a transform will reappear another place at random using a different transformation, and this is explained very well in the second of the Godel, Escher Bach lectures given at MIT.  You will also notice that the Sierpinski triangle is also explained in the same lecture. To see the mapping function, you can color the different values as they are implemented. 

There are several versions of different ferns in the source code, and they are all taken from this web site. They are not completely the same, as the mapping technique is slightly different in the source code. As a side note Ferns are actually one of the earliest plant structures found in nature.

L-systems  and Fractals 


The study of the growth of plants and other organic structures, and especially on how to model the growth processes of cell of plants, lead the biologist Aristid Lindenmayer to develop a new method for modelling growth of these structures in 1968. To see how he thought about creating this relative simple system he thought somewhere along the lines given below.

"In many growth processes of living organisms, especially of plants, regularly repeated appearances of certain multicellular structures are readily noticeable.... In the case of a compound leaf, for instance, some of the lobes (or leaflets), which are parts of a leaf at an advanced stage, have the same shape as the whole leaf has at an earlier stage."

At first he was using his theory to study the development of simple multicellular organisms, but it was subsequently applied to investigate higher plants and plant organs. The actual process could be formulated along these lines; Assume that you have a starting point, either some predetermined axiom (or a starting gene in the DNA), that form the basis of how it will initiate the growth. Further it will also have some rules of how it will expand from the axiom (starting condition) over a discrete period of n (Involving the implementation of the other DNA pattern). These basic principles and ideas were already discovered in the beginning of the 20th century, but was not widely known or used before Noam Chomsky started to use a similar system to describe language syntax in 1957.  The Lindenmayer system (L-system for short) is however different from language in the replacements algorithm, while the Chomsky algorithm is sequential the Lindenmayer system is in parallel. It is this difference that make the L-system so effective for growth patterns in fractals, while the Chomsky algorithm is used nearly only for language syntax.

The first thing that is needed is to generate the growth pattern, by the axiom and the replacements rules. For example set the axiom to F--F--F, this will generate a triangle, provided that the "-" sign stands for 60 degrees. To generate a Koch snowflake you will, in each iteration replace the value F with F+F--F+F (written as F=F+F--F+F), and it gets more and more complicated as you go on. This can be coded in the following way: 

Private Sub Grow()
    Dim Current, [Next] As New StringBuilder

    Dim Found As Boolean = False
    Dim CurrentGeneration As Integer = 0

    While CurrentGeneration < Generations

        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
                    Found = True
                End If
            If Not Found Then
            End If
        CurrentGeneration += 1

    End While
    GrowString = ([Next].ToString)
End Sub   

After the GrowthString is generated, all that is left is to do the actual drawing of the fractal generated. The 2D L-system has 6 different operations, and these are often visualized by a turtle moving in different directions while it is drawing a line while it moves. The basic draw commands are: 

  • "F" - Draw a line   
  • "f" - Move the turtle without drawing 
  • "+" - Add a predetermind angle to the current angle 
  • "-" - Withdraw a predetermined angle to the current angle 
  • "[" - Save the current position to a stack 
  • "]" - Remove the current top point in the stack and sets it as the current position to the turtle

All other symbols are ignored by the turtle (the turtle preserves its current state) and that is very easy to code, assuming that you have created the iterated tree and replaced axioms used in the rewriting system.  

Public Function GetDrawingLines() As List(Of Line)
    'Make the tree from axioms and rules

    'Replace temporary axioms with new ones

    '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
        ElseIf GrowString(i) = "]"c Then
            CurrentTurtle = PointList.Pop
        End If

    Return (result)
End Function

There are a lot of implementations of this algorithm online, but as I mention earlier they are (in my opinion) overly complicated, and might confuse people who want to implement and understand them. It is quite remarkable that these simple rules can generate quite complex structures, and even quite complex trees.  

There is another way of generate a Tree fractal by the use of recursive algorithm instead, as you could draw it by a command similar to:

Private sub DrawBranch(Angle, MaxIterations)
End Sub

This is exactly the approach taken by Mohammad Reza Khosravi in his article DotNet Real Tree, where he allows you many different settings to generate a wide range of different trees. There is one thing that I haven mentioned, and that is how a natural tree will decrease its cross section when it splits into branches. The formula that was proposed by Leonardo Da Vinci, and considered to be a natural truth, were d<sub>1</sub> is the crossection area before the split, and d<sub>2</sub> and d<sub>3</sub> are the cross section area of each of the two new branches, one gets the following relation:  

d<sub>1</sub> = Math.Sqrt(d<sub>2</sub><sup>2</sup> + d<sub>3</sub><sup>2</sup>) 

This relation was proposed in order so that the sap could reach the top of the trees. If the two branches, d<sub>2</sub> and d<sub>3</sub> are equal, the equation suggests that each branch would be d1/Math.Sqrt(2) thinner.

Since L-systems could  generate quite a vide range of natural organisms and plants in nature, people started to wonder what will happen when you would combine the iterating possibilities with a natural selection, meaning choices, like humans could do. One simple scheme of doing this is by a suggestion by Richard Dawkins and an example of this type of algorithm called AI: Dawkins Biomorphs / And Other Evolving Creatures by Sacha Barber.

I should also mention that you could also implement the L-system in three dimensions, which would leave the string generator unchanged, but would add two additional angles. An example of this is the Hilbert curve in 3D as could be seen here. I did not implement it as it would take a lot of coding and wouldn't show any further deep insights into the L-system. There is however a book, written by Lindenmayer et. al. that could be downloaded for free, where a lot of 3D structures appear here.  

To show how little that could differentiate different shapes I have set the axiom to this:


And the exchange rule for each iteration to this:


Given 4 iterations and with the angle equal to 60 degrees we get a Koch like flake: 

Set the angle to 90 degrees instead gives this pattern:

Yet again, set the angle to 120 degrees we get a version of the Sierpinski triangle:

All of these three patterns a quite distinct, but the only thing that has change is the angle. They do however overlap, meaning that they draw the same line multiple times, but it does show how little it takes to create quite different shapes.

Kronecker products and fractals  


To generate a Kronecker fractal all one really needs to do is to create the algorithm for the Kronecker product between two matrices as defined in this Wikipedia article, and this could be done the following way: 

''' <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)
    Return R
End Function 

Now all you need to do is to fill two matrices with values from 0 upto 1 and multiply them together.  To draw the image make 0 white and 1 black and thats all there is to it really.

There are also other possibilities for implementing the drawing and allow numbers outside the range from 0 and 1 but then you would need to change the code accordingly.

DLA - Diffusion Limited Aggregation


DLA is one of the easiest fractals to understand, and to generate it one of then uses the following analogy; A drunken man walking about looking for an inn, and he finds one and just stays there. He starts form an arbitrary point, and moves to one of his nine neighbors at random, and this is repeated until he reaches a neighbor that has an inn.   

The code to generate a starting point for the drunken man is given below:

''' <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>

To move, one needs to choose one arbitrary point that is next to the current point. I have selected to simulate the rectangle that overlaps itselfs, like a sphere like structure. If it wants to go to the right and it currently stands in the rightmost point, it moves to the start on the left hand side etc. The code is given below:

''' <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
                x -= 1
            End If
        Case 2
            'Move to point x;y = 9;9
            If x = 0 Then
                x = a - 1
                x -= 1
            End If

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

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

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

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

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

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

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

This algorithm could be altered in many ways, but the example I have in my, program is the simplest form possible, and many enhancements are possible, as this example shows (code in Java). It is also a process that can simulate many phenomena in nature.  

Coloring sheme 

I have made use of two different coloring schemas in this article, one that generates a distinct and unique color that is different from the surrounding coloring. I have now also added a smooth coloring scheme, that lets you set gradient points of colors, and it simply iterates linear from one color to the next, generating a smooth looking graphics. This is simply done by using this code:

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)))

End Sub

This has however some drawbacks as the coloring is not matched to the given iterations in the Mandelbrot set, thus you must take care so that there is enough colors in the list to all the given iterations in the different fractals.

The colors can be exchanged by simply altering altering the colors in the gradient stop here:

Sub New()

    For i As Integer = 0 To GradientStop.Count - 2
        GetGradients(GradientStop(i).Color, GradientStop(i + 1).Color, 35)

End Sub


As you probably have seen the reason that fractals are so interesting is that it can be used to simulate plants, growth patterns in nature, meaning that from some simple axiom one can generate quite complex forms.

Initial release: 16 September 2013. 


I start off with the most important references here:

  • The MIT "Godel, Echer, Back" Lecture series 
  • "Mathematical tools in computer graphics with C# implementations", written by Alexandre Hardy & Willi-Hans Steeb , World Scientific  
  • The source code from "Mathematical tools in computer graphics with C# implementations"  can be downloaded from here.
  • "Fractals, Chaos, Power Laws" by Manfred Schroeder. This is a very readable book about, with lots of information of the different types of fractals.  

Other useful online resources are:


This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


About the Author

Kenneth Haugland
Norway Norway
I hope that you like the stuff I have created and if you do wish to say thank you then a donation is always appreciated.
You can donate here[^].

You may also be interested in...


Comments and Discussions

QuestionThanks Pin
Member 1172085626-May-15 11:29
memberMember 1172085626-May-15 11:29 
AnswerRe: Thanks Pin
Kenneth Haugland26-May-15 12:10
professionalKenneth Haugland26-May-15 12:10 
GeneralMy vote of 5 Pin
Athari4-Oct-13 4:04
professionalAthari4-Oct-13 4:04 
GeneralRe: My vote of 5 Pin
Kenneth Haugland4-Oct-13 4:06
professionalKenneth Haugland4-Oct-13 4:06 
Questionvote Pin
arussell28-Sep-13 9:22
professionalarussell28-Sep-13 9:22 
AnswerRe: vote Pin
Kenneth Haugland28-Sep-13 9:39
professionalKenneth Haugland28-Sep-13 9:39 
GeneralMy vote of 5 Pin
dyma25-Sep-13 5:13
memberdyma25-Sep-13 5:13 
GeneralRe: My vote of 5 Pin
Kenneth Haugland25-Sep-13 5:17
professionalKenneth Haugland25-Sep-13 5:17 
GeneralMy vote of 5 Pin
Wooters24-Sep-13 17:23
memberWooters24-Sep-13 17:23 
GeneralRe: My vote of 5 Pin
Kenneth Haugland25-Sep-13 0:14
professionalKenneth Haugland25-Sep-13 0:14 
GeneralMy vote of 5 Pin
SeptimusHedgehog 15157620-Sep-13 11:45
memberSeptimusHedgehog 15157620-Sep-13 11:45 
GeneralRe: My vote of 5 Pin
Kenneth Haugland21-Sep-13 0:33
professionalKenneth Haugland21-Sep-13 0:33 
GeneralMy vote of 5 Pin
Southmountain20-Sep-13 5:56
memberSouthmountain20-Sep-13 5:56 
GeneralRe: My vote of 5 Pin
Kenneth Haugland20-Sep-13 5:59
professionalKenneth Haugland20-Sep-13 5:59 
GeneralMy vote of 5 Pin
Paulo Zemek19-Sep-13 18:20
professionalPaulo Zemek19-Sep-13 18:20 
GeneralRe: My vote of 5 Pin
Kenneth Haugland19-Sep-13 19:31
professionalKenneth Haugland19-Sep-13 19:31 
GeneralMy vote of 5 Pin
rspercy6518-Sep-13 13:51
memberrspercy6518-Sep-13 13:51 
GeneralRe: My vote of 5 Pin
Kenneth Haugland18-Sep-13 23:19
professionalKenneth Haugland18-Sep-13 23:19 
GeneralMy vote of 5 Pin
Mr. Xie.G.Gang18-Sep-13 3:58
memberMr. Xie.G.Gang18-Sep-13 3:58 
GeneralRe: My vote of 5 Pin
Kenneth Haugland18-Sep-13 6:18
professionalKenneth Haugland18-Sep-13 6:18 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Permalink | Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.170628.1 | Last Updated 28 Sep 2013
Article Copyright 2013 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2017
Layout: fixed | fluid