One of the features that drew me to F# was the ease of interoperability with code written in other CLR based languages, allowing me to easily re-use existing code with F# as a concise and expressive glue language to write the core of the application – something I’ve used IronPython for in the past.
Continuing with the option pricing theme, I’m going to put together an application to visualize the sensitivity of the option price to the various inputs in the pricing. Windows Presentation Foundation seemed to be an easy way of achieving this, so this post walks through the creation of a small F# application to plot a function of two variables as a 3D surface.
Getting started
Being familiar with 3D graphics, but never having used WPF from F# before (or indeed any other language) it took me a little while to get my bearings. After creating a new F# project the following steps were necessary before I could use get going:
- Add references to WindowsBase, PresentationCore, PresentationFramework and System.Xml
- Downloaded 3D Tools from and added a reference to that library
3DTools provides two pieces of functionality that will be used in this application – a trackball that allows the camera to be rotated by clicking and dragging in the viewport, and a mechanism by which to draw lines of constant screen-space width (ScreenSpaceLines3D).
Constructing an approximation of the surface
The API for displaying a surface will accept a set of points on the X-axis, a second set of points on the Z-axis and a function y = f(x,z) to calculate the Y-axis displacement for each (X,Z) point. For each (Xi, Zi) two triangles will be created in the square made up of (Xi, Zi), (Xi+1,Zi), (Xi+1,Zi+1) and (Xi,Zi+1), e.g.
The following code constructs a mesh as described above:
/// Evaluates the surface defined by f at each point in the cartesian
/// product of xs and zs. xs and zs are required to be ordered in ascending
/// numerical order.
let EvaluateSurface (xs:float array) (zs:float array) (f: float -> float -> float) =
if zs.Length < 2 then invalidArg "zs" "Cannot construct surface with fewer than two points on the z-axis"
if xs.Length < 2 then invalidArg "xs" "Cannot construct surface with fewer than two points on the x-axis"
// Find the smallest and largest value on each axis; the members of
// each axis array must be passed in ascending order.
let xMin = xs.[0]
let xMax = xs.[xs.Length-1]
let xRange = xMax - xMin
let zMin = zs.[0]
let zMax = zs.[zs.Length-1]
let zRange = zMax - zMin
// Evaluate the function f for each pairing of points in xs and ys. Can
// be done with List.iter but the imperative form generates simpler,
// more efficient IL
let points = Point3DCollection()
let mutable yMin = Double.MaxValue
let mutable yMax = Double.MinValue
for x in xs do
for z in zs do
let y = f x z
if y < yMin then yMin <- y elif y > yMax then yMax <- y
points.Add(Point3D(x, y, z))
// Reposition the points in object space so (yMax + yMin)/2 is
// positioned at y=0, and similarly for the X and Y axes.
let yRange = yMax - yMin
let translation = Vector3D( -(xMax+xMin)/2.0, -(yMax+yMin)/2.0, -(zMax+zMin)/2.0 )
let scale = Vector3D( 1.0 / xRange, 1.0 / yRange, 1.0 / zRange )
let texture = PointCollection()
for i in 0..points.Count-1 do
let p = points.[i]
points.[i] <- Point3D( (p.X + translation.X) * scale.X,
(p.Y + translation.Y) * scale.Y,
(p.Z + translation.Z) * scale.Z )
texture.Add(Point((p.Y-yMin)/yRange, 0.0))
// We now have |xs| * |zs| points. We need to form a triangle mesh
// over the top. For each index i, we construct two triangles from
// points (i, i+1, i+|xs|) and (i+1, i+|xs|, i+|xs|+1)
let numSquares = xs.Length * (zs.Length-1)
let triangles = Int32Collection(numSquares*6)
for i in 0..zs.Length-2 do
let row = i * xs.Length
for j in 0..xs.Length-2 do
// For convenience of notation, give the indices of the four
// corners of the square labels.
let bl = row+j
let br = bl+1
let tl = bl+xs.Length
let tr = tl+1
// Construct two triangles (bl, tl, br) and (tl, tr, br)
triangles.Add(bl)
triangles.Add(tl)
triangles.Add(br)
triangles.Add(tl)
triangles.Add(tr)
triangles.Add(br)
let mesh = MeshGeometry3D()
mesh.Positions <- points
mesh.TextureCoordinates <- texture
mesh.TriangleIndices <- triangles
mesh
Adding contour lines
To make it easy to see the relative height of different points on the surface and to highlight the shape of the surface, it will be coloured according to its displacement on the Y-axis. An easy way of doing this is by assigning UV texture coordinates to each point based on the ratio of the Y coordinate to the Y-axis range - (Y-YMIN)/(YMAX-YMIN) – as seen in the code above. A GradientBrush can then be used to specify the colour for each Y value. Contour lines can be added by adding thin bands of a single colour to the gradient brush. The following code will construct such a brush:
/// Creates a gradient brush containing equally-sized stripes of each
/// colour in the list, separated by a thin gray contour line. If the
/// blend parameter is true then each band will blend smoothly from one
/// colour to the next. If the blend parameter is false then each band
/// will be a solid colour.
let CreateGradientBrush (stripes:int list) blend =
/// Utility method to convert from an integer xxRRGGBB representation
/// to a WPF Color
let inline ColorFromRRGGBB hex =
let r = byte((hex >>> 16) &&& 0xff)
let g = byte((hex >>> 8) &&& 0xff)
let b = byte(hex &&& 0xff)
Color.FromRgb(r, g, b)
let gradient = LinearGradientBrush()
gradient.StartPoint <- Point(0.0, 0.5)
gradient.EndPoint <- Point(1.0, 0.5)
// A small change in texture space used to create a visually sharp
// edge between two adjacent colours that shouldn't appear blended
// into each other.
let epsilon = 0.0000001
let lineWidth = 0.003
let lineColour = Color.FromRgb(0x66uy, 0x66uy, 0x66uy)
let limit = float(stripes.Length)
for i in 0..stripes.Length-2 do
let c = ColorFromRRGGBB stripes.[i]
let c' = if blend then ColorFromRRGGBB stripes.[i+1] else c
let bottom = float(i) / limit
let top = (float(i+1) / limit) - epsilon
let lineMin = top - lineWidth
let colorMax = lineMin - epsilon
gradient.GradientStops.Add(GradientStop(c, bottom))
gradient.GradientStops.Add(GradientStop(c', colorMax))
gradient.GradientStops.Add(GradientStop(lineColour, lineMin))
gradient.GradientStops.Add(GradientStop(lineColour, top))
// Handle the final colour separately so that there isn't a gray
// contour line at the top - it just finishes in the final colour.
let numStops = gradient.GradientStops.Count
let lastOffset = gradient.GradientStops.[numStops-1].Offset
let topColour = ColorFromRRGGBB stripes.[stripes.Length-1]
gradient.GradientStops.Add(GradientStop(topColour, lastOffset+epsilon))
gradient.GradientStops.Add(GradientStop(topColour, 1.0))
gradient
Rendering the mesh
Having constructed the model and the texture that we are going to apply, then there are not many more steps required to render it on the screen:
1. Apply the gradient to the mesh
To apply the gradient to the mesh, we need to create a GeometryModel3D that
binds together the mesh with the material. In order for both sides of the mesh to be rendered
it is necessary to set both the front and the back material:
let mesh = EvaluateSurface xs ys f let gradient = CreateGradientBrush [0x62a4ce;0x79aed0;0x90b9d3;0xa7c3d6;0xbeced9] false let model = GeometryModel3D() model.Geometry <- mesh model.Material <- DiffuseMaterial(gradient) model.BackMaterial <- model.Material
2. Define the scene
In order to see the surface, it is necessary to add some lighting to the scene. To apply a uniform light to all faces an AmbientLight can be used. The model and the light are bound together in a single Model3DGroup (because we don't need to move or transform them separately) and then wrapped in a ModelVisual3D for display.
let light = AmbientLight() light.Color <- Colors.White let group = Model3DGroup() group.Children.Add( model ) group.Children.Add( light ) let visual = ModelVisual3D() visual.Content <- group
3. Create the user interface
The final steps are to create aViewport3D, add the ModelVisual3D to it, add a camera and place it in a Window to be displayed. Wrapping the Viewport3D in an instance of the 3DTools TrackballDecorator is all that is required to allow the user to rotate the camera around the origin:
// Create a camera through which to view the model let camera = PerspectiveCamera() camera.FieldOfView <- 60.0 camera.Position <- Point3D(-1., 1., 1.5) camera.LookDirection <- Vector3D(1., -1. -1.5) // Create a viewport to display the visual let viewport = Viewport3D() viewport.Camera <- camera viewport.Children.Add( visual ) let trackball = TrackballDecorator() trackball.Content <- viewport let window = Window() window.Width <- 800.0 window.Height <- 800.0 window.Title <- "Surface plot" window.Content <- trackball window.Show() Application().Run() |> ignore
The results
You can download the code here and do with it as you wish. If you make any enhancements it would be good if you could let me know. Here are a couple of examples of the output:
open System
open Psi.Charting
open System.Windows
open Utils
[<EntryPoint>]
[<STAThread>]
let main args =
// Display two windows - both plotting the same function but over
// different ranges and using different colour schemes.
let f = (fun x y -> Math.Cos(x) + Math.Cos(y))
let p1 = plot_surface (linspace -1.0 1.0 40) (linspace -1.0 1.0 40) f
let p2 = plot_surface (linspace -4.5 4.5 100) (linspace -4.5 4.5 100) f
p2.SetColourScheme(AlternativeColourScheme, true)
Application().Run() |> ignore
0
Plot 1: Showing cos(x)+cos(y) over the range -1..1 on both axes
Plot 2: Showing cos(x)+cos(y) over the range -4.5 to 4.5 on both axes
No comments:
Post a Comment