Mathias Brandewinder on .NET, F#, VSTO and Excel development, and quantitative analysis / machine learning.
by Mathias 11. January 2015 18:37

I had the great pleasure to speak at CodeMash this week, and, on my way back, ended up spending a couple of hours at the Atlanta airport waiting for my connecting flight back to the warmer climate of San Francisco – a perfect opportunity for some light-hearted coding fun. A couple of days earlier, I came across this really nice tweet, rendering the results of an L-system:

I ended up looking up L-systems on Wikipedia, and thought this would make for some fun coding exercise. In a nutshell, a L-system is a grammar. It starts with an alphabet of symbols, and a set of rules which govern how each symbol can be transformed into another chain of symbols. By applying these rules to a starting state (the initial axiom), one can evolve it into a succession of states, which can be seen as the growth of an organism. And by mapping each symbol to operations in a logo/turtle like language, each generation can then be rendered as a graphic.

So how could we go about coding this in F#? If you are impatient, you can find the final result as a gist here.

First, I started with representing the core elements of an L-System with a couple of types:

type Symbol = | Sym of char

type State = Symbol list

type Rules = Map<Symbol,State>

type LSystem = 
    { Axiom:State
      Rules:Rules }

A symbol is a char, wrapped in a single-case discriminated union, and a State is simply a list of Symbols. We define the Rules that govern the transformation of Symbols by a Map, which associates a particular Symbol with a State, and an L-System is then an Axiom (the initial State), with a collection of Rules.

Let’s illustrate this on the second example from the Wikipedia page, the Pythagoras tree. Our grammar contains 4 symbols, 0, 1, [ and ], we start with a 0, and we have 2 rules, (1 → 11), and (0 → 1[0]0). This can be encoded in a straightforward manner in our domain, like this:

let lSystem =
    { Axiom = [ Sym('0') ]
      Rules = [ Sym('1'), [ Sym('1'); Sym('1') ]
                Sym('0'), [ Sym('1'); Sym('['); Sym('0'); Sym(']'); Sym('0') ]]
              |> Map.ofList }

Growing the organism by applying the rules is fairly straightforward: given a State, we traverse the list of Symbols, look up for each of them if there is a matching rule, and perform a substitution if it is found, leaving it unchanged otherwise:

(*
Growing from the original axiom
by applying the rules
*)

let applyRules (rs:Rules) (s:Symbol) =
    match (rs.TryFind s) with
    | None -> [s]
    | Some(x) -> x

let evolve (rs:Rules) (s:State) =
    [ for sym in s do yield! (applyRules rs sym) ]

let forward (g:LSystem) =
    let init = g.Axiom
    let gen = evolve g.Rules
    init |> Seq.unfold (fun state -> Some(state, gen state))

// compute nth generation of lSystem
let generation gen lSystem =
    lSystem
    |> forward 
    |> Seq.nth gen
    |> Seq.toList

What does this give us on the Pythagoras Tree?

> lSystem |> generation 1;;
val it : Symbol list = [Sym '1'; Sym '['; Sym '0'; Sym ']'; Sym '0']

Nice and crisp – that part is done. Next up, rendering. The idea here is that for each Symbol in a State, we will perform a substitution with a sequence of instructions, either a Move, drawing a line of a certain length, or a Turn of a certain Angle. We will also have a Stack, where we can Push or Pop the current position of the Turtle, so that we can for instance store the current position and direction on the stack, perform a couple of moves with a Push, and then return to the previous position by a Pop, which will reset the turtle to the previous position. Again, that lends itself to a very natural model:

(*
Modelling the Turtle/Logo instructions
*)

type Length = | Len of float
type Angle = | Deg of float

// override operator later
let add (a1:Angle) (a2:Angle) =
    let d1 = match a1 with Deg(x) -> x
    let d2 = match a2 with Deg(x) -> x
    Deg(d1+d2)

type Inst =
    | Move of Length
    | Turn of Angle
    | Push
    | Pop

let Fwd x = Move(Len(x))
let Lft x = Turn(Deg(x))
let Rgt x = Turn(Deg(-x))

We can now transform our L-system state into a list of instructions, and convert them into a sequence of Operations, in that case Drawing lines between 2 points:

type Pos = { X:float; Y:float; }
type Dir = { L:Length; A:Angle }

type Turtle = { Pos:Pos; Dir:Dir }
type ProgState = { Curr:Turtle; Stack:Turtle list }

let turn angle turtle = 
    let a = turtle.Dir.A |> add angle
    { turtle with Dir = { turtle.Dir with A = a } }

type Translation = Map<Symbol,Inst list>

type Ops = | Draw of Pos * Pos

let pi = System.Math.PI

let line (pos:Pos) (len:Length) (ang:Angle) =
    let l = match len with | Len(l) -> l
    let a = match ang with | Deg(a) -> (a * pi / 180.)
    { X = pos.X + l * cos a ; Y = pos.Y + l * sin a }

let execute (inst:Inst) (state:ProgState) =
    match inst with
    | Push -> None, { state with Stack = state.Curr :: state.Stack }
    | Pop -> 
        let head::tail = state.Stack // assumes more Push than Pop
        None, { state with Curr = head; Stack = tail }
    | Turn(angle) -> 
        None, { state with Curr =  state.Curr |> turn angle }
    | Move(len) -> 
        let startPoint = state.Curr.Pos
        let endPoint = line startPoint len state.Curr.Dir.A
        Some(Draw(startPoint,endPoint)), { state with Curr = { state.Curr with Pos = endPoint } }

let toTurtle (T:Translation) (xs:Symbol list) =

    let startPos = { X = 400.; Y = 400. }
    let startDir = { L = Len(0.); A = Deg(0.) }
    let init = 
        { Curr = { Pos = startPos; Dir = startDir }
          Stack = [] }
    xs 
    |> List.map (fun sym -> T.[sym]) 
    |> List.concat
    |> Seq.scan (fun (op,state) inst -> execute inst state) (None,init)
    |> Seq.map fst
    |> Seq.choose id

We simply map each Symbol to a List of instructions, transform the list of symbols into a list of instructions, and maintain at each step the current position and direction, as well as a Stack (represented as a list) of positions and directions. How does it play out on our Pythagoras Tree? First, we define the mapping from Symbols to Instructions:

let l = 1.
let T = 
    [ Sym('0'), [ Fwd l; ]
      Sym('1'), [ Fwd l; ]
      Sym('['), [ Push; Lft 45.; ]
      Sym(']'), [ Pop; Rgt 45.; ] ]
    |> Map.ofList

… and we simply send that toTurtle, which produces a list of Draw instructions:

> lSystem |> generation 1 |> toTurtle T;;
val it : seq<Ops> =
  seq
    [Draw ({X = 400.0;
            Y = 400.0;},{X = 401.0;
                         Y = 400.0;}); Draw ({X = 401.0;
                                              Y = 400.0;},{X = 401.7071068;
                                                           Y = 400.7071068;});
     Draw ({X = 401.0;
            Y = 400.0;},{X = 401.7071068;
                         Y = 399.2928932;})]

Last step – some pretty pictures. We’ll simply generate a html document, rendering the image using SVG, by creating one SVG line per Draw instruction:

let header = """
<!DOCTYPE html>
<html>
<body>
<svg height="800" width="800">"""

let footer = """
</svg>
</body>
</html>
"""

let toSvg (ops:Ops seq) =
    let asString (op:Ops) = 
        match op with
        | Draw(p1,p2) -> sprintf """<line x1="%f" y1="%f" x2="%f" y2="%f" style="stroke:rgb(0,0,0);stroke-width:1" />""" p1.X p1.Y p2.X p2.Y 

    [ yield header
      for op in ops -> asString op
      yield footer ]
    |> String.concat "\n"

open System.IO

let path = "C:/users/mathias/desktop/lsystem.html"
let save template = File.WriteAllText(path,template)

And we are pretty much done:

> lSystem |> generation 8 |> toTurtle T |> toSvg |> save;;
val it : unit = ()

… which produces the following graphic:

image

Pretty neat! Just for fun, I replicated the Sierpinski Triangle example as well:

let sierpinski () =

    let lSystem =
        { Axiom = [ Sym('A') ]
          Rules = [ Sym('A'), [ Sym('B'); Sym('>'); Sym('A'); Sym('>'); Sym('B') ]
                    Sym('B'), [ Sym('A'); Sym('<'); Sym('B'); Sym('<'); Sym('A') ]]
                  |> Map.ofList }

    let l = 1.
    let T = 
        [ Sym('A'), [ Fwd l; ]
          Sym('B'), [ Fwd l; ]
          Sym('>'), [ Lft 60.; ]
          Sym('<'), [ Rgt 60.; ] ]
        |> Map.ofList

    lSystem 
    |> generation 9
    |> toTurtle T
    |> toSvg 
    |> save

… which results in the following picture:

image

That’s it for tonight! I had a lot of fun coding this (it certainly made the flight less boring), and found the idea of converting code to turtle instructions, with a stack, pretty interesting. Hope you enjoyed it, and if you end up playing with this, share your creations on Twitter and ping me at @brandewinder!

Gist for the whole code here

by Mathias 25. March 2012 10:53

In a previous post, I looked at creating a Sierpinski triangle using F# and WPF. One of the pieces I was not too happy about was the function I used to transform a Triangle into a next generation triangle:

type Point = { X:float; Y:float }
type Triangle = { A:Point; B:Point; C:Point }

let transform (p1, p2, p3) =
   let x1 = p1.X + 0.5 * (p2.X - p1.X) + 0.5 * (p3.X - p1.X)
   let y1 = p1.Y + 0.5 * (p2.Y - p1.Y) + 0.5 * (p3.Y - p1.Y)
   let x2 = p1.X + 1.0 * (p2.X - p1.X) + 0.5 * (p3.X - p1.X)
   let y2 = p1.Y + 1.0 * (p2.Y - p1.Y) + 0.5 * (p3.Y - p1.Y)
   let x3 = p1.X + 0.5 * (p2.X - p1.X) + 1.0 * (p3.X - p1.X)
   let y3 = p1.Y + 0.5 * (p2.Y - p1.Y) + 1.0 * (p3.Y - p1.Y)
   { A = { X = x1; Y = y1 }; B = { X = x2; Y = y2 }; C= { X = x3; Y = y3 }}

Per se, there is nothing wrong with the transform function: it takes 3 points (the triangle corners), and returns a new Triangle. However, what is being “done” to the triangle is not very expressive – and the code looks rather ugly, with clear duplication (the exact same operation is repeated on the X and Y coordinates of every point).

Bringing back blurry memories from past geometry classes, it seems we are missing the notion of a Vector. What we are doing here is taking corner p1 of the Triangle, and adding a linear combinations of the edges p1, p2 and p1, p3 to it, which can be seen as 2 Vectors (p2 – p1) and (p3 – p1). Restated that way, here is what the transform function is really doing:

A –> A + 0.5 x AB + 0.5 x AC

A –> A + 1.0 x AB + 0.5 AC

A –> A + 0.5 x AB + 1.0 x AC 

In graphical form, the first transformation can be represented as follows:

image

In order to achieve this, we need to define a few elements: a Vector, obviously, a way to create a Vector from two Points, to add Vectors, to scale a Vector by a scalar, and to translate a Point by a Vector. Let’s do it:

type Vector = 
   { dX:float; dY:float }
   static member (+) (v1, v2) = { dX = v1.dX + v2.dX; dY = v1.dY + v2.dY }
   static member (*) (sc, v) = { dX = sc * v.dX; dY = sc * v.dY }

type Point = 
   { X:float; Y:float }
   static member (+) (p, v) = { X = p.X + v.dX; Y = p.Y + v.dY }
   static member (-) (p2, p1) = { dX = p2.X - p1.X; dY = p2.Y - p1.Y }

type Triangle = { A:Point; B:Point; C:Point }

Thanks to operators overloading, the transform function can now be re-phrased in a much more palatable way:

let transform (p1:Point, p2, p3) =
   let a = p1 + 0.5 * (p2 - p1) + 0.5 * (p3 - p1)
   let b = p1 + 1.0 * (p2 - p1) + 0.5 * (p3 - p1)
   let c = p1 + 0.5 * (p2 - p1) + 1.0 * (p3 - p1)
   { A = a; B = b; C = c }

… and we are done. The code (posted on fsSnip.net) works exactly as before, but it’s way clearer.

It can also be tweaked more easily now. I got curious about what would happen if slightly different transformations were applied, and the results can be pretty fun. For instance, with a minor modification of the transform function…

let transform (p1:Point, p2, p3) =
   let a = p1 + 0.55 * (p2 - p1) + 0.5 * (p3 - p1)
   let b = p1 + 1.05 * (p2 - p1) + 0.45 * (p3 - p1)
   let c = p1 + 0.5 * (p2 - p1) + 0.95 * (p3 - p1)
   { A = a; B = b; C = c }

… we get the following, bloated “Sierpinski triangle”:

BeerPinski-triangle

Add a bit of transparency, some more tweaks of the linear combinations,

let transform (p1:Point, p2, p3) =
   let a = p1 + 0.3 * (p2 - p1) + 0.6 * (p3 - p1)
   let b = p1 + 0.8 * (p2 - p1) + 0.3 * (p3 - p1)
   let c = p1 + 0.6 * (p2 - p1) + 1.1 * (p3 - p1)
   { A = a; B = b; C = c }

and things get much wilder:

SnowflakePinski-triangle

I don’t think these are really Sierpinski triangles any more, but I had lots of fun playing with this, and figured someone else might enjoy it, too… If you find a nice new combination, post it in the comments!

Source code: fsSnip.net

by Mathias 16. March 2012 08:33

In my last post, I looked into drawing a Sierpinski triangle using F# and WinForms, and noted that the rendering wasn’t too smooth – so I converted it to WPF, to see if the result would be any better, and it is. In the process, I discovered John Liao’s blog, which contains some F# + WPF code examples I found very useful. I posted the code below, as well as on FsSnip. The differences with the WinForms code are minimal, I’ll let the interested reader figure that part out!

One thing I noticed is that the starting point of the Sierpinski sequence is a single triangle – but nothing would prevent a curious user to initialize the sequence with multiple triangles. And while at it, why not use WPF Brush opacity to create semi-transparent triangles, and see how their superposition looks like?

We just change the Brush Color and Opacity, and add a second triangle to the root sequence…

let brush = new SolidColorBrush(Colors.DarkBlue)
brush.Opacity <- 0.6
let renderTriangle = render canvas brush

let triangle = 
    let p1 = { X = 190.0; Y = 170.0 }
    let p2 = { X = 410.0; Y = 210.0}
    let p3 = { X = 220.0; Y = 360.0}
    { A = p1; B = p2; C = p3 }

let triangle2 =
    let p1 = { X = 290.0; Y = 170.0 }
    let p2 = { X = 510.0; Y = 210.0}
    let p3 = { X = 320.0; Y = 360.0}
    { A = p1; B = p2; C = p3 }

let root = seq { yield triangle; yield triangle2 }

… and here we go:

Sierpinski-superposition

Granted, it’s pretty useless, but I thought it looked rather nice!

As an aside, here is something I noted when working in F#: I often end up looking at the code, thinking “can I use this to do something I didn’t think about when I wrote it”? In C#, I tend to think in terms of restrictions: write Components, with a “containment” approach – figure out what the component should do, and enforce safety by constraining the inputs/outputs via an interface. By contrast, because of type inference and the fact that a function doesn’t require an “owner” (it is typically not a member of a class), I find myself less “mentally conditioned”, and instead of a world of IWidgets and ISprockets, I simply see functions that transform elements, and wonder what else they could apply to.

The case we saw here was trivial, but pretty much from the moment I wrote that code, I have been mulling over other extensions. What is the transform function really doing, and what other functions could I replace it with? generateFrom is simply permuting the triangle corners and applying the same transformation – could I generalize this to an arbitrary sequence and write Sierpinski Polygons? Could I even apply it to something that has nothing to do with geometry?

// Requires reference to 
// PresentationCore, PresentationFramework, 
// System.Windows.Presentation, System.Xaml, WindowsBase

open System
open System.Windows
open System.Windows.Media
open System.Windows.Shapes
open System.Windows.Controls

type Point = { X:float; Y:float }
type Triangle = { A:Point; B:Point; C:Point }

let transform (p1, p2, p3) =
   let x1 = p1.X + 0.5 * (p2.X - p1.X) + 0.5 * (p3.X - p1.X)
   let y1 = p1.Y + 0.5 * (p2.Y - p1.Y) + 0.5 * (p3.Y - p1.Y)
   let x2 = p1.X + 1.0 * (p2.X - p1.X) + 0.5 * (p3.X - p1.X)
   let y2 = p1.Y + 1.0 * (p2.Y - p1.Y) + 0.5 * (p3.Y - p1.Y)
   let x3 = p1.X + 0.5 * (p2.X - p1.X) + 1.0 * (p3.X - p1.X)
   let y3 = p1.Y + 0.5 * (p2.Y - p1.Y) + 1.0 * (p3.Y - p1.Y)
   { A = { X = x1; Y = y1 }; B = { X = x2; Y = y2 }; C= { X = x3; Y = y3 }}

let generateFrom triangle = seq {
      yield transform (triangle.A, triangle.B, triangle.C)
      yield transform (triangle.B, triangle.C, triangle.A)
      yield transform (triangle.C, triangle.A, triangle.B)
   }

let nextGeneration triangles =
   Seq.collect generateFrom triangles 
      
let render (target:Canvas) (brush:Brush) triangle =
   let points = new PointCollection()
   points.Add(new System.Windows.Point(triangle.A.X, triangle.A.Y))
   points.Add(new System.Windows.Point(triangle.B.X, triangle.B.Y))
   points.Add(new System.Windows.Point(triangle.C.X, triangle.C.Y))
   let polygon = new Polygon()
   polygon.Points <- points
   polygon.Fill <- brush
   target.Children.Add(polygon) |> ignore
   
let win = new Window()
let canvas = new Canvas()
canvas.Background <- Brushes.White
let brush = new SolidColorBrush(Colors.Black)
brush.Opacity <- 1.0
let renderTriangle = render canvas brush

let triangle = 
    let p1 = { X = 190.0; Y = 170.0 }
    let p2 = { X = 410.0; Y = 210.0}
    let p3 = { X = 220.0; Y = 360.0}
    { A = p1; B = p2; C = p3 }

let root = seq { yield triangle }
let generations = 
   Seq.unfold (fun state -> Some(state, (nextGeneration state))) root
   |> Seq.take 7
Seq.iter (fun gen -> Seq.iter renderTriangle gen) generations

win.Content <- canvas
win.Show()

[<STAThread()>]
do 
   let app =  new Application() in
   app.Run() |> ignore
by Mathias 11. March 2012 11:29

I am midway through the highly-recommended “Real-world functional programming: with examples in C# and F#”, which inspired me to play with graphics using F# and WinForms (hadn’t touched that one in a long, long time), and I figured it would be fun to try generating a Sierpinski Triangle.

The Sierpinski Triangle is generated starting from an initial triangle. 3 half-size copies of the triangle are created and placed outside of the original triangle, each of them having a corner “touching” the middle of one side of the triangle.

Sierpinski-original

That procedure is then repeated for each of the 3 new triangles, creating more and more smaller triangles, which progressively fill in an enclosing triangular shape. The figure below uses the same starting point, stopped after 6 “generations”:

 

Sierpinski-6

The Sierpinski Triangle is an example of a fractal figure, displaying self-similarity: if we were to run the procedure ad infinitum, each part of the Triangle would look like the whole “triangle” itself.

So how could we go about creating this in F#?

I figured this would be a good case for Seq.unfold: given a state (the triangles that have been produced for generation n), and an initial state to start from, provide a function which defines how the next generation of triangles should be produced, defining a “virtual” infinite sequence of triangles – and then use Seq.take to request the number of generations to be plotted.

Here is the entire code I used; you’ll need to add a reference to System.Windows.Forms and System.Drawings for it to run:

open System
open System.Drawing
open System.Windows.Forms

type Point = { X:float32; Y:float32 }
type Triangle = { A:Point; B:Point; C:Point }

let transform (p1, p2, p3) =
   let x1 = p1.X + 0.5f * (p2.X - p1.X) + 0.5f * (p3.X - p1.X)
   let y1 = p1.Y + 0.5f * (p2.Y - p1.Y) + 0.5f * (p3.Y - p1.Y)
   let x2 = p1.X + 1.0f * (p2.X - p1.X) + 0.5f * (p3.X - p1.X)
   let y2 = p1.Y + 1.0f * (p2.Y - p1.Y) + 0.5f * (p3.Y - p1.Y)
   let x3 = p1.X + 0.5f * (p2.X - p1.X) + 1.0f * (p3.X - p1.X)
   let y3 = p1.Y + 0.5f * (p2.Y - p1.Y) + 1.0f * (p3.Y - p1.Y)
   { A = { X = x1; Y = y1 }; B = { X = x2; Y = y2 }; C= { X = x3; Y = y3 }}

let generateFrom triangle = seq {
      yield transform (triangle.A, triangle.B, triangle.C)
      yield transform (triangle.B, triangle.C, triangle.A)
      yield transform (triangle.C, triangle.A, triangle.B)
   }

let nextGeneration triangles =
   Seq.collect generateFrom triangles 
      
let render (target:Graphics) (brush:Brush) triangle =
   let p1 = new PointF(triangle.A.X, triangle.A.Y)
   let p2 = new PointF(triangle.B.X, triangle.B.Y)
   let p3 = new PointF(triangle.C.X, triangle.C.Y)
   let points = List.toArray <| [ p1; p2; p3 ]
   target.FillPolygon(brush, points)
   
let form = new Form(Width=500, Height=500)
let box = new PictureBox(BackColor=Color.White, Dock=DockStyle.Fill)
let image = new Bitmap(500, 500)
let graphics = Graphics.FromImage(image)
let brush = new SolidBrush(Color.FromArgb(0,0,0))
let renderTriangle = render graphics brush

let p1 = { X = 190.0f; Y = 170.0f }
let p2 = { X = 410.0f; Y = 210.0f}
let p3 = { X = 220.0f; Y = 360.0f}
let triangle = { A = p1; B = p2; C = p3 }

let root = seq { yield triangle }
let generations = 
   Seq.unfold (fun state -> Some(state, (nextGeneration state))) root
   |> Seq.take 7
Seq.iter (fun gen -> Seq.iter renderTriangle gen) generations

box.Image <- image
form.Controls.Add(box)

[<STAThread>]
do
   Application.Run(form)

A few comments on the code. I first define 2 types, a Point with 2 float32 coordinates (float32 chosen because that’s what System.Drawing.PointF takes), , and a Triangle defined by 3 Points. The transform function is pretty ugly, and can certainly be cleaned up / prettified. It takes a tuple of 3 Points, and returns the corresponding transformed triangle, shrunk by half and located at the middle of the p1, p2 edge. We can now build up on this with the nextGeneration function, which takes in a Sequence of Triangles (generation n), transforms each of them into 3 new Triangles and uses collect to “flatten” the result into a new Sequence, generation n+1.

The rendering code has been mostly lifted with slight modifications from Chapter 4 of Real-world functional programming. The render function retrieves the 3 points of a Triangle and create a filled Polygon, displayed on a Graphics object which we create in a form.

Running that particular example generates the following Sierpinski triangle; you can play with the coordinates of the root triangle, and the number of generations, to build your own!

As an aside, I was a bit disappointed by the quality of the graphics; beyond 7 or 8 generations, the result gets fairly blurry. I’ll probably give a shot at moving this to XAML, and see if it’s any better.

 

Sierpinski-example

As usual, comments, questions or criticisms are welcome!

Comments

Comment RSS