F# Flocking (boids)

A couple of years ago I was blogging about some simple boid behaviour I had implemented in C# (here). I came across this paper which had some interesting ideas I wanted to try out: "Self-organised complex aerial displays of thousands of starlings: a model" by H. Hildenbrandt, C. Carere, C-K. Hemelrijk.

As I'm getting more into F#, I thought this would be a simple reasonable sized application to create to get a feel for how F# works on a slightly larger project (though definitely not large-scale).

Before taking a look at the code, I'll quickly describe which parts of the paper I implemented.

The video generated is available on http://www.youtube.com/watch?v=eil5K7Ir3i8

I haven't implemented all of the ideas in the paper, so my simple application doesn't exhibit the same realistic flocking behaviour, but it does have more realistic behaviour than my earlier efforts. The interesting behaviours were the desire to not stray too far from the roost, the attempt to maintain a steady cruise speed, and to use the nearest seven topological neighbours for cohesion and alignment. I implemented the simulation in 2D for simplicity.

I hardcoded a perception radius instead of implementing the continual adjustment of the perception radius for seven topological neighbours. I also have totally omitted the flight dynamic simulation (no gravity, no banked turns). The paper discusses that the banked turns gives a further air of realism to the simulation.  

The coding

This was a real pleasure to implement in F# - a lot less time was thinking about application design, classes and their interactions, the logic in the paper was easy to transfer directly to the keyboard. F# is very compact, so the core logic can be seen on just a couple of screens.

I didn't feel the lack of intellisense too much when developing the core algorithm, I did miss the ability to refactor functions and navigate to usages. I especially missed intellisense when developing the WPF visualisation part of the app, when interacting with .NET objects; I did miss the ability for Visual Studio to automatically add the required namespaces. I must have being spoilt by Visual Studio and Resharper for too long!

The actual WPF application wasn’t such a good fit for F# – there’s no generated code behind file for the XAML, and I feel that using F# would be painful in a WPF or Silverlight application (but just for the view, it should be OK for the ViewModels down).

I implemented a F# Vector type, which can be specified with units of measure (or none). I used units of measure throughout – this was really powerful, and did let me quickly find a few bugs in the implementation.

(NOTE: I still haven’t found a code websnippet tool I’m happy with – you need to click in each of the regions below and scroll down and right to see the whole code). Alternately, the zipped up solution can be downloaded from http://www.taumuon.co.uk/jabuka/FSharpFlock.zip 

Vector3D.fs:

module Vector3D

type Vector3D<[<Measure>] 'u>(x : float<'u>, y : float<'u>, z : float<'u>) =
static member Zero() = Vector3D<_>(0.0<_>, 0.0<_>, 0.0<_>)

member v.X = x
member v.Y = y
member v.Z = z

static member (+) (lhs:Vector3D<_>, rhs:Vector3D<_>) =
Vector3D(lhs.X + rhs.X, lhs.Y + rhs.Y, lhs.Z + rhs.Z)

static member (-) (lhs:Vector3D<_>, rhs:Vector3D<_>) =
Vector3D(lhs.X - rhs.X, lhs.Y - rhs.Y, lhs.Z - rhs.Z)

static member (*) (v:Vector3D<_>, a:float<_>) =
Vector3D(v.X * a, v.Y * a, v.Z * a)

static member (*) (a:float<_>, v:Vector3D<_>) =
Vector3D(a * v.X, a * v.Y, a * v.Z)

static member (/) (v:Vector3D<_>, a) =
Vector3D(v.X / a, v.Y / a, v.Z / a)

member v.DotProduct(rhs:Vector3D<_>) = (v.X * rhs.X) + (v.Y * rhs.Y) + (v.Z * rhs.Z)

member v.magnitude = sqrt(v.DotProduct(v)) * 1.0<_>

member lhs.CrossProduct(rhs:Vector3D<_>) =
Vector3D((lhs.Y * rhs.Z - lhs.Z * rhs.Y) * 1.0<_>,
(-lhs.X * rhs.Z + lhs.Z * rhs.X) * 1.0<_>,
(lhs.X * rhs.Y - lhs.Y * rhs.X) * 1.0<_>)

member v.normalise =
let magnitude = float v.magnitude
Vector3D<_>((v.X / magnitude), (v.Y / magnitude), (v.Z / magnitude))

let sumVectors(vectors : Vector3D<_>[]) =
let initial = Vector3D<_>(0.0<_>, 0.0<_>, 0.0<_>)
Array.fold (+) initial vectors



BoidUtils.fs:



module BoidUtils

open Microsoft.FSharp.Math
open Vector3D
open SI

let radiusRoost = 150.0<m>
let hardRadius = 2.0<m> // 0.2<m>
let mass = 0.08<kg>
let timeStep = 0.005<s>
let relaxationTime = 0.05<s>
let cruiseSpeed = 20.0<m/s>
let horizontalRoostWeighting = 0.01<N/m>
let weightingAlignment = 0.5<kg * s^-2>
let weightingCohesion = 1.0<kg s^-2>
let weightingSeparation = 2.0
let perceptionRadius = 50.0<m>

type BodyAxes =
{ Forward:Vector3D<1>;
Side:Vector3D<1>;
Up:Vector3D<1> }

type Boid =
{ Position:Vector3D<m>;
Speed:float<m / s>;
Orientation:BodyAxes; }

// All parameterless functions are evaluated once, just on module opening, so pass random in.
let InitialiseRandomPosition(rand:System.Random) =
Vector3D<m>((300.0 * (rand.NextDouble()-0.5)) * 1.0<m>,
((300.0 * (rand.NextDouble()-0.5)) * 1.0<m>),
0.0<m>)

let InitialiseRandomVelocity(rand:System.Random) =
Vector3D<m/s>((100.0 * (-0.5 + rand.NextDouble()) * 1.0<m/s>),
(100.0 * (-0.5 + rand.NextDouble())) * 1.0<m/s>,
0.0<m/s>)

let InitialiseRandomOrientation(rand:System.Random) =
{Forward=Vector3D(0.0, 1.0, 0.0);
Side=Vector3D(1.0, 0.0, 0.0);
Up=Vector3D(0.0, 0.0, 1.0)}

let setOrientation(oldOrientation:BodyAxes, velocity:Vector3D<m/s>) =
let normalisedVelocity = velocity.normalise
let y = normalisedVelocity.CrossProduct(Vector3D<m / s>(0.0<m/s>, 0.0<m/s>, 1.0<m/s>))
{oldOrientation with Forward=normalisedVelocity * 1.0<m^-1 s>; Side=y*1.0<m^-2 s^2>}

let calculateCruiseSpeedForce (boid:Boid) =
(mass / relaxationTime) * (cruiseSpeed - boid.Speed) * boid.Orientation.Forward

let calculateRoostForce (boid:Boid) =
let horizontalPosition = Vector3D(boid.Position.X, boid.Position.Y, 0.0<_>)
let distanceFromOrigin = horizontalPosition.magnitude
match (distanceFromOrigin) with
| _ when distanceFromOrigin < radiusRoost -> Vector3D<N>(0.0<N>, 0.0<N>, 0.0<N>)
| _ -> let normalRoostingArea = horizontalPosition.normalise
let d = boid.Orientation.Forward.DotProduct normalRoostingArea
let distanceFromRoost = distanceFromOrigin - radiusRoost
let orientationRoostDotProduct = boid.Orientation.Side.DotProduct normalRoostingArea
let weightingFactor = match (orientationRoostDotProduct) with
| _ when orientationRoostDotProduct > 0.0<m> -> -1.0
| _ -> 1.0
weightingFactor * (radiusRoost * horizontalRoostWeighting * (0.5 + (0.5<m^-1> * d)) * (boid.Orientation.Side))

let findDistanceBetweenBoids(boid:Boid, other:Boid) =
(boid.Position - other.Position).magnitude

let findNearestNeighbours(boid:Boid, boids:Boid[]) =
let sortedByDistance = boids |> Array.sortBy(fun other -> findDistanceBetweenBoids(boid, other))
Array.sub sortedByDistance 0 7

let findAverageForwardDirectionDifference(boid:Boid, boids:Boid[]) =
let differences = boids |> Array.map (fun i -> 1.0<m> * (i.Orientation.Forward - boid.Orientation.Forward))
let sumDifferences = sumVectors(differences)
(1.0 / (float sumDifferences.magnitude)) * sumDifferences

let calculateAlignmentForce(boid:Boid, nearest:Boid[]) =
let averageDifference = findAverageForwardDirectionDifference(boid, nearest)
weightingAlignment * averageDifference

let findAveragePosition(boid:Boid, boids:Boid[]) =
let positions = boids |> Array.map (fun i -> i.Position)
let sumPositions = sumVectors(positions)
(1.0 / float boids.Length) * sumPositions

let findNeighboursInRadius(boid:Boid, boids:Boid[], radius:float<m>) =
boids |> Array.filter(fun other -> other <> boid && findDistanceBetweenBoids(boid, other) <= radius)

let calculateCentrality(boid:Boid, boids:Boid[]) =
let separations = boids |> Array.map(fun i -> (i.Position - boid.Position).normalise)
let sumSeparations = sumVectors(separations)
let count = boids.Length
match (count) with
| 0 -> 1.0
| _ -> (1.0 / float count) * (sumSeparations.magnitude / 1.0<m>)

let calculateCohesionForce(boid:Boid, nearest:Boid[], boidsInPerceptionRadius:Boid[]) =
let boidsOutsideHardRadius = nearest |> Array.filter(fun i -> abs ((boid.Position - i.Position).magnitude) > hardRadius)
let centrality = calculateCentrality(boid, boidsInPerceptionRadius)
let averagePosition = findAveragePosition(boid, nearest)
centrality * weightingCohesion * (averagePosition - boid.Position)

let calculateSeparationForce(boid:Boid, boidsInPerceptionRadius:Boid[]) =
let nearest = boidsInPerceptionRadius
let separations = nearest |> Array.map(fun i -> i.Position - boid.Position)
let sigma = 1.8
let forcesToNeighbours = separations |> Array.map(fun i ->
let magnitude = i.magnitude
let multiplier =
match (magnitude) with
| _ when magnitude < hardRadius -> 1.0
| _ -> System.Math.Exp(-((magnitude - hardRadius)*(magnitude - hardRadius)/1.0<m^2>) / (sigma * sigma))
multiplier * magnitude * (i.normalise) * 1.0<kg * m^-1 * s^-2>)
let sumForces = sumVectors(forcesToNeighbours)
match (nearest.Length) with
| _ when (nearest.Length) = 0 -> Vector3D<N>.Zero()
| _ -> (-weightingSeparation / float nearest.Length) * sumForces

let calculateSocialForce(boid:Boid, boids:Boid[]) =
let nearest = findNearestNeighbours(boid, boids)
let boidsInPerceptionRadius = findNeighboursInRadius(boid, boids, perceptionRadius)
calculateAlignmentForce(boid, nearest)
+ calculateCohesionForce(boid, nearest, boidsInPerceptionRadius)
+ calculateSeparationForce(boid, boidsInPerceptionRadius)

let calculateForce (boid:Boid, boids:Boid[]) =
(boid |> calculateRoostForce)
+ (boid |> calculateCruiseSpeedForce)
+ (calculateSocialForce(boid, boids))

let iterateBoid (boid:Boid, boids:Boid[]) =
let originalPosition = boid.Position
let originalVelocity = boid.Speed * boid.Orientation.Forward
let force = calculateForce(boid, boids)
let acceleration = force/mass
let velocity = originalVelocity + (acceleration * timeStep)
let position = originalPosition + (velocity * timeStep)
let newOrientation = setOrientation(boid.Orientation, velocity)
{Position=position;Speed=velocity.magnitude;Orientation=newOrientation}




 

Program.fs:

module Boids

open SI
open Vector3D
open BoidUtils

open System
open System.IO
open System.Windows
open System.Windows.Threading
open System.Windows.Controls
open System.Windows.Shapes
open System.Windows.Media
open System.Windows.Media.Imaging

let window = Application.LoadComponent(new System.Uri("/FSharpFlock;component/MainWindow.xaml",
System.UriKind.Relative)) :?> Window

let rand = new System.Random()
let viewWidth = 480.0
let viewHeight = 360.0

let mutable frameCount = 0
let mutable boids =
[|for i in 0 .. 300 ->
let position = InitialiseRandomPosition(rand)
let velocity = InitialiseRandomVelocity(rand)
let tempOrientation = {Forward=Vector3D(0.0, 1.0, 0.0);
Side=Vector3D(1.0, 0.0, 0.0);
Up=Vector3D(0.0, 0.0, 1.0)}
let orientation = setOrientation(tempOrientation, velocity)
{Position = position; Speed=velocity.magnitude; Orientation=orientation}
|];;

let updateBoids(boids) =
boids |> Array.map (fun boid -> iterateBoid(boid, boids))

let (?) (fe:FrameworkElement) name : 'T =
fe.FindName(name) :?> 'T

let GetRotationForBoid(boid:Boid) =
let forward = boid.Orientation.Forward
let angleRadians = Math.Atan2(forward.X, forward.Y)
let angleDegrees = angleRadians * 180.0 / Math.PI
let rotateTransform = new Media.RotateTransform(angleDegrees, 0.0, 0.0)
rotateTransform

let color = Media.Color.FromArgb(64uy, 255uy, 0uy, 0uy)

let saveFrame(canvas : Canvas) =
let size = new Size(viewWidth, viewHeight)
canvas.Measure(size)
canvas.Arrange(new Rect(size))
let renderTargetBitmap = new RenderTargetBitmap(int viewWidth, int viewHeight, 96.0, 96.0, PixelFormats.Pbgra32)
let sourceBrush = new VisualBrush(canvas)
renderTargetBitmap.Render(canvas)

let bitmapFrame = BitmapFrame.Create(renderTargetBitmap)
let jpegEncoder = new JpegBitmapEncoder()
jpegEncoder.Frames.Add(bitmapFrame)
let filename = System.String.Format("C:\Anim\{0:0000}.jpg", frameCount)
use stream = new FileStream(filename, FileMode.CreateNew)
jpegEncoder.Save(stream)

let createBoidGraphics(boid:Boid) =
let obj : Polygon = new Polygon()
obj.Points <- new Media.PointCollection( [|new Point(-10.0, -5.0);
new Point(0.0, 0.0);
new Point(-10.0, 5.0);
new Point(-5.0, 0.0)|] )
obj.RenderTransform <- GetRotationForBoid(boid)
obj.Fill <- new Media.SolidColorBrush(color)
obj.Stroke <- Media.Brushes.Black
obj.StrokeThickness <- 1.0
obj

let drawBoids() =
let win : Window = window
let canvas : Canvas = win?canvas
canvas.Children.Clear()
for i = 0 to boids.Length - 1 do
let graphicalBoid = createBoidGraphics(boids.[i])
let unitlessPosition = boids.[i].Position * 1.0<m^-1>
System.Windows.Controls.Canvas.SetTop(graphicalBoid, (viewHeight / 2.0) + (unitlessPosition.X))
System.Windows.Controls.Canvas.SetLeft(graphicalBoid, (viewWidth / 2.0) + (unitlessPosition.Y))
canvas.Children.Add(graphicalBoid) |> ignore
// saveFrame(canvas)
frameCount <- frameCount + 1

let timer = new DispatcherTimer();
timer.Tick.Add(fun evArgs ->
[0..5]
|> Seq.iter(fun x -> boids <- boids |> updateBoids)
drawBoids()
)
timer.Interval <- TimeSpan.FromMilliseconds(1.0);

let setup(win: Window) =
win.Loaded.Add(fun _ -> timer.Start())

[<STAThread>]
[<EntryPoint>]
let main args =
setup window
(new Application()).Run(window) |> ignore
0


MainWindow.xaml:



<Window xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
Title="F# Flock" Height="360" Width="480">
<Canvas Name="canvas" HorizontalAlignment="Stretch" VerticalAlignment="Stretch"
Background="LightBlue"/>
</Window>



The rest of the code I hope is simple to read and understand. Coding this up I had some of the most fun programming for a while. With F#, I found I was able to concentrate on the core of the algorithm and very quickly see results.



For this type of project, very algorithmic, F# was a perfect fit. I'd like to get a feel for how to implement F# on a large scale project, and to see how it feels to expose the functionality via an object oriented layer.

Labels: , ,