Wednesday, 23 December 2009

Surface plotting redux

Time has been short for the last fortnight, so haven't made as much progress as I'd like on this project. This short post covers one topic: adding labels to the plot from the last post. Here's an example of the result, using a 'pretty' function from F# for Visualization:

An Example Surface

To render text in a Viewport3D, this code follows roughly the same approach described by Eric Sink:

  1. Calculate the extents of the label in 3D space
  2. Create a rectangular "billboard" at the desired location
  3. "Texture Map" the desired text onto the billboard using a VisualBrush and a TextBlock

Positioning the text

To align the center of the label with the grid lines, the width of the text needs to be estimated. This can be done using GlyphTypeface, which provides the advance widths for individual characters:

module Wpf = 
  /// Attempts to retrieve a GlyphTypeface for the specified family,
  /// style and weight. If no matching typeface can be obtained then
  /// this raises an invalid argument exception.
  let GetTypeface (family:FontFamily)(style:FontStyle)(weight:FontWeight)=
    let t = Typeface(family, style, weight, FontStretches.Normal)
        
    let success, typeface = t.TryGetGlyphTypeface()
    if not success then
      invalidArg "family" "Cannot find a GlyphTypeface for the specified font family, style and weight"
            
    typeface

  /// Measures the length of a line of text in the specified typeface, 
  /// assuming a size of 1.0. Callers can scale the resulting width 
  /// according to the size they are using.
  let MeasureLine (text:string) (f:GlyphTypeface) = 
    Seq.fold (fun w (c:char) -> w + f.AdvanceWidths.[f.CharacterToGlyphMap.[int(c)]]) 0.0 text

Creating the billboard

The position of the label in three dimensional space is done using its center point and two vectors - one that points along the line of text and one that points from the bottom to the top of the line. Together the two vectors define the plane on which the text is rendered.

Texture coordinates are assigned to vertices, so in order to allow the text to be viewed the "right way around" from above and underneath we need to create two rectangles on which to display the text, with distinct vertices.

The extents of the rectangle are calculated by moving half of the text width away from the specified center point along the 'along' axis and half of the height on the 'up' axis. Two rectangular billboards are then created and painted using a visual brush containing a text block.

/// Creates a label oriented on a plane defined by the center point of the
/// text and two unit vectors defining its orientation and plane.
let CreateLabel (text:string) (center:Point3D) (up:Vector3D) (along:Vector3D) =           
    let height      = 0.035 // hard-coded font height, 
    let width       = height * Wpf.MeasureLine text this.typeface
    let bottomLeft  = center - width / 2.0 * along - height / 2.0 * up;
    let bottomRight = bottomLeft + (along * width);
    let topLeft     = bottomLeft + (up * height);
    let topRight    = topLeft + (along * width);

    let points = Point3DCollection()
    let rectangle = MeshGeometry3D()
    rectangle.Positions <- points

    // Shorthand for adding a point and associated UV coordinate
    let point p u v =
        points.Add(p)
        rectangle.TextureCoordinates.Add(Point(u, v))

    // Shorthand for creating a triangle from three point indices
    let triangle p0 p1 p2 = 
        rectangle.TriangleIndices.Add(p0)
        rectangle.TriangleIndices.Add(p1)
        rectangle.TriangleIndices.Add(p2)

    // Add the points and texture coordinates for the underside 
    point bottomLeft  0.0 1.0
    point topLeft     0.0 0.0
    point bottomRight 1.0 1.0
    point topRight    1.0 0.0
        
    // Add the points and texture coordinates for the topside
    point bottomLeft  1.0 1.0
    point topLeft     1.0 0.0
    point bottomRight 0.0 1.0
    point topRight    0.0 0.0

    // Define the triangles that make up the two sides of the rectangle
    triangle 0 3 1
    triangle 0 2 3
    triangle 4 5 7
    triangle 4 7 6        

    let block = TextBlock(Run(text))
    block.Foreground <- this.brush
    block.FontFamily <- this.family

    // Hint that the visual brush is not going to change, so that it is
    // not re-rendered on every frame.
    let brush = VisualBrush(block)
    RenderOptions.SetCachingHint(brush, CachingHint.Cache)

    let visual = ModelVisual3D()
    visual.Content <- GeometryModel3D(rectangle, DiffuseMaterial(brush))
    visual        

Performance

Having written the small amount of code described above, I fired up the application to be met with abysmally jittery rendering. Three small changes restored bearable performance:

Enabling caching for the VisualBrush

By default, the contents of the VisualBrush are rendered for every frame. WPF allows the application to hint that the Visual underlying the VisualBrush is not frequently changing and can be cached. If the hint is obeyed, then the Visual is rendered off-screen, reducing the cost to regular texture mapping.

This has an advantage over rendering the text to a bitmap yourself, since WPF will re-render the Visual at different resolutions as its projection onto the screen changes in size to ensure that the appearance doesn't significantly degrade.

To use this hint requires a single line, although the corresponding MSDN page describes how the results can be tuned:

let brush = VisualBrush(textblock)
RenderOptions.SetCachingHint(brush, CachingHint.Cache)

Disabling hit-testing

Leaving hit-testing enabled has a noticable impact on performance. Since it's not used in this application, it can be disabled on the viewport with:

viewport.IsHitTestVisible <- false

Disabling clipping

There isn't a significant amount of content outside the window and the Viewport3D fills the window, so clipping can be disabled for an additional speed up. The Maximize WPF 3D Performance page on MSDN indicates that the Viewport3D's clipping is slow. Disabling it is a single-line change:

viewport.ClipToBounds <- false

Updated code

The updated code can be downloaded from here.

Saturday, 5 December 2009

Plotting a surface using WPF 3D from F#

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 a Viewport3D, 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