Game Developers Conference 2012 - Technical summary

GDC2012 is over, and this year there are a huge number of available presentations. You can download the Game Developer Conference 2012 presentations from the GDC vault, Jare / Iguana has also kept a link collection from GDC 2012. I've looked over all the technical publications available and put together this summary post. (Edit: I've updated this to cover some maths, physics, and graphics material I missed on the first pass - thanks Johan & Eric)

I'll start with graphics.

Louis Bavoil / nVidia and Johan Andersson / DICE have a presentation on "Stable SSAO in Battlefield 3 with Selective Temporal Filtering", ambient occlusion is a well established technique now, but they apply a quick way to use past data and the differences in Z buffer states between frames to intelligently reuse the AO results. They also look at filters and optimising blur functions. Similar to established tricks in the realtime raytracing demoscene.

Eban Cook / Naughty Dog presented "Creating Flood Effects in Uncharted 3", a technical artist look at water effects. Unfortunately realtime fluid simulation wasn't used, instead Houdini was used to pre-generate the game content. An overview of the shaders for water, water particles, froth particles, and lighting is given.

Light probe interpolation

Robert Cupisz / Unity discussed light probes, "Light probe interpolation using tetrahedral tessellations", in terms of choosing the appropriate probe and weights using Delaunay Triangulation / Tetrahedrons and Barycentric Coordinates by dividing scenes into convex hulls. Also covers projecting onto the nearest convex hull, covers it all with a fair bit of maths, this would be of interest to physics / collision detection programmers too. There is a collection of nice links and some sample code at the end.

Matthijs De Smedt / Nixxes covers "Deus Ex is in the Details" using DX11 tech. Covers AA (FXAA DLAA MLAA), SSAO, DOF (Gaussian blur), tessellation and soft-shadows.

Colt McAnlis / Google investigates post-compressing DXT textures in his talk "DXT is not enough", trying to out-do zipped DXT's with delta encoding. More info at this blog post or skip it all and download the DXT CRUNCH compressor here.

Matt Swoboda / Sony & Fairlight delves into Signed Distance Fields, a demoscene hot-topic last year, with the talk "Advanced Procedural Rendering in DirectX 11". Investigates converting polgyon mesh data and particle data into signed distance fields. Takes an in-depth look into a optimised marching cubes implementation for a fluid simulation with smooth particle hydrodynamics (SPH), and how to use signed distance fields to do ambient occlusion.

Physically based rendering in realtime 

Yoshiharu Gotanda / tri-Ace research makes a case for physically based rendering with a Blinn-Phong model in the presentation "Practical Physically Based Rendering in Real-time". An indepth look at the BRDF formulation they use.

Wolfgang Engel, Igor Lobanchikov and Timothy Martin / Confetti present "Dynamic Global Illumination from many Lights", just a bunch of pictures, not much information.

Carlos Gonzalez Ochoa / Naughty Dog covers "Water Technology of Uncharted". Covers the shader, animating the normal maps flow, and simulating the ocean water with Gerstner waves, b-spline waves, and wave particles. They go on to look at LOD with "Irregular Geometry Clipmaps" including fixing T joints, and then culling, skylights and underwater fog. Next physics, attaching objects (buoyant), and point queries. Finally, SPU optimization. Quite comprehensive.

Water technology of Uncharted

Ben Hanke / SlantSixGames describes the bone code in "Rigging a Resident Evil". Transforms are described with 9 functions and processed with an optimising compiler, allowing fast retargeting of animations.

Scott Kircher / Volition Inc expands on Inferred Lighting in "Lighting & Simplifying Saints Row: The Third" by looking into lighting for rain, foliage, dynamic decals, and radial ambient occlusion. Then moves on to automated mesh simplification using iterative edge contraction and takes an indepth look at selecting an appropriate error metric.

Nathan Reed / Sucker Punch Productions discusses "Ambient Occlusion Fields and Decals in Infamous 2", going into depth on how to solve the artifacts of this approach.

Marshall Robin / Naughty Dog covers the effects system tools in "Effects Techniques Used in Uncharted 3: Drake’s Deception".

Niklas Smedberg / Epic Games looks at PowerVR GPU processing pipeline and capabilities in "Bringing AAA graphics to mobile platforms" and provides a number of tricks'n'tips on optimising the performance of the mobile GPU, and highlights the cheap operations. In short: AA (fast), hidden surface removal (fast), alpha test (slow), render targets (slow), texture lookups (slow). Takes a more detailed look at material shaders, god rays, and character shadows. All in all, pretend its ~2002, and you'll be right.

Mickael Gilabert / Ubisoft and Nikolay Stefanov / Massive cover the GI system in Far Cry 3 in "Deferred Radiance Transfer Volumes". Light probes get precomputed directional radiance transfer data from a custom raytracer stored using spherical harmonics. Source code for the relighting system is presented, along with optimisations by using volume textures.

John McDonald / nVidia explains CPU/GPU synching for buffers in "Don’t Throw it all Away: Efficient Buffer Management" and provides advice on buffer creation flags.

Bryan Dudash / nVidia suggests using average normals to overcome tesselation issues in "My Tessellation Has Cracks!".

Mastering DX11 with Unity

Renaldas Zioma / Unity and Simon Green / nVidia present "Mastering DirectX 11 with Unity". Starts by looking into Unity's physically based shaders (Oren-Nayar, Cook-Torrance, and energy conservation, then blurry reflections and combining normal maps). Next up, Catmull-Clark Subdivision, tetrahedra light probes (See Robert Cupisz's talk), HBAO, APEX destruction, Hair simulation with guide hairs, Explosions using signed distance fields with noise and color gradients,and finally velocity buffer motion blur.

Tobias Persson / Bitsquid discusses lighting billboards in "Practical Particle Lighting". Looks at normal generation and per-pixel lighting for billboards (including code snippets), applying shadow maps with a domain shader, and shadow casting

Karl Hillesland / AMD investigates realtime Ptex (per-face texture mapping) in "Ptex and Vector Displacement in AMD Demos", and efficient retrieval from the texture atlas, including all MIPs.

Jay McKee / AMD presents the "Technology Behind AMD’s
Leo Demo". He details some of the code behind the forward rendering of 3000 dynamic light sources using a depth pre-pass, light culling (tile-based compute shader to output light list), and light accumulation with materials phase. Basically moves the light management code from CPU to GPU.

Terrain in Battlefield 3

Mattias Widmark / DICE presents "Terrain in Battlefield 3: A modern, complete and scalable system". Begins with an overview of the features for the terrain system (heightfield based, procedurally generated, spline decals, decoration (tree,rock,grass), destruction, water), and presents their quadtree terrain data structure, paying particular attention to LOD. Next, CPU/GPU performance is investigated, and a clip-map based virtual texture system is presented. The large terrain data set is managed by intelligently streaming data to the required detail ('blurriness'), and co-locating data (heighfield/color/mask lumped together, next to the next level of relevant LOD data), which is also compressed (RLE/DXT1). Nodes are then prioritized based on distance, culling, and updates (e.g. destruction). Finally, mesh generation, stitching and tessellation with displacement on the GPU.




Moving on to physics.

Erin Catto covers "Diablo 3 Ragdolls", including representing ragdoll bones, initialising ragdolls from animations, and interacting kinematic and dynamic objects.

François Antoine / Epic talks about Gears of Wars 3 destruction physics in "Pushing for Large Scale Destruction FX" and suggests using particles for dust and debrie, and simplifying meshes for destruction.

Stephen Frye / EA looks at ragdolls in the presentation "Tackling Physics". Highlights aspects of ragdolls that look unrealistic, and suggests adding joint limits and motorized constraints at joints to simulate muscles. Gives two approaches to solving the control problem, first using external forces, second calculating the appropriate torque from world space.

Graham Rhodes / Applied Research Associates presents "Computational Geometry" where he looks at half-edge data structures for triangulating a polygon, splitting a face, splitting an edge, intersection of an edge and a plane and generating a convex hull.

Richard Tonge / nVidia covers "Solving Rigid Body Contacts" and starts with a gentle introduction to rigid body state space and progressively builds a signal-block-diagram of solving a single contact restraint. Then looks at each block in the diagram and deciphers the physics behind it. He then looks at solving multiple contacts, and explains why you can't apply a linear solver to the problem (contacts break), and presents the LCP, and an alternative approach; sequential impulses. He then gives a whirlwind tour of GPU solvers.

Gino van den Bergen / DTECTA presents "Collision Detection", first covering shapes, then configuration space, distance tests, Seperating Axis Tests, and takes a closer look at the GJK algorithm.

Jim Van Verth / Insomniac gives a nice introduction to Navier Stokes in "Fluid Techniques", breaking down the terms for external forces, viscocity, advection and pressure visually. Then looks at three major representations for fluids, grid, particle and surface (wave) based.

Takahiro Harada / AMD examines how heterogeneous compute architectures can achieve large scale dynamic simulations in "Toward A Large Scale Simulation". Begins with an overview of GPU architecture, and GPU rigid body simulation in three key phases: broad-phase, narrow-phase and constraint solving for a system of 128,000 particles and 12,000 convex bodies. He presents a design for overcoming data transfer and minimising synchronisation points whilst dividing the workload between CPU and GPU.

Erwin Coumans / AMD investigates destructive physics in the aptly titled "Destruction". He begins with generating voroni diagrams for shattering geometry and boolean operations, and moves into generating collision shapes with convex decomposition and tetrahedralization. Then moves on to realtime approaches with real-time booleans and breakable constraints and finite element
methods.



Looking at AI.

Bobby Anguelov / IO Interactive, Gabriel Leblanc / Eidos-Montréal and Shawn Harris / Big Huge Games present "Animation-Driven Locomotion For Smoother Navigation". They start with the standard motion graphs and transitioning/blending between animation cycles. Then take an indepth look at footstep planning (IK, Foot sliding) and come up with a system for deciding where steps should be taken to fulfil the navigation goal. They then investigate modifying navigation paths to better fit the animation cycles, and finish by looking into collision avoidance.

Daniel Brewer / Digital Extremes looks at agent perception, reaction, combat chatter, buddy systems and collision avoidance using velocity space Optimal Reciprocal Collision Avoidance in "Building Better Baddies".

Brian Magerko / Georgia Tech covers "How to Teach Game AI from Scratch" including competitions (Mario AI, Google Ants, Poker AI, Starcraft AI).

Dave Mark / Intrinsic and Kevin Dill / Lockheed Martin investigate some examples (snipers, guards) of Utility-Based AI in "Embracing The Dark Art of Mathematical Modeling in Game AI".

Kasper Fauerby / IO Interactive explains "Crowds in Hitman:Absolution" including cell maps, boids, animation and PS3 implementation details. The crowd AI uses a state machine with steering behaviours (pending walk, walk, panic), and behaviour 'zones' with information from the navigation system to select behaviours. Near-optimal Character Animation with Continuous Control was used for animation.

Elan Ruskin / Valve looks at empowering writers and dialog in TF2, Left4Dead, etc, in "Rule Databases for Contextual Dialog and Game Logic". Begins with player triggered lines (extended by environment, memory etc.) and avoiding fill-in-the-blank dialog by using databases. Rules, queries, responses and writers tools are examined next, and ties things off with database query optimisations.

Mike Robbins / Gas Powered Games examines "Neural Networks in Supreme Commander 2", with 34 inputs and 15 output actions and a single hidden layer (98 neurons), with a fitness function composed from 17 inputs trained to control combat platoons.

Ben Sunshine-Hill / Havok investigates LOD for AI in "Perceptually Driven Simulation", and makes a case for using probability of noticing a difference instead of distance as a LOD measure, and presents a market-based "LOD trader" for selecting the appropriate LOD given the constraints on hand.

Moving along to programming and math.

Adisak Pochanayon / Netherrealm covers debugging and timing issues in "Runtime CPU Spike Detection using Manual and Compiler-Automated Instrumentation". First up, manual instrumentation and wrapper functions, Then detours, and automated instrumentation (compiler flags) with an indepth look at the 360. Finally, profiling with threshold functions.

Pete Isensee / Microsoft details how rvalue in C++11 (T&&) can eliminate temporaries in "Faster C++: Move Construction and Perfect Forwarding".

Scott Selfon / Microsoft reviews audio compression technologies in "The State of Ady0 Cmprshn", starting with time-domain compression with PCM (raw, A-Law, U-Law, ADPCM), then frequency-domain compression and discusses the artifacts generated by both, then evaluates the performance of different codecs.

Robin Green / Microsoft and Manny Ko / Dreamworks present "Frames, Quadratures and Global Illumination: New Math for Games". Beings with a review of spherical harmonics, Haar wavelets, and Radial basis functions. Builds up to 'Spherical Needlet' wavelets, by exploring different basis functions ('frames')

Gino van den Bergen / DTECTA presents dual-numbers in "Math for Game Programmers: Dual Numbers", beginning with a look at complex numbers. Automatic differentiation with dual numbers is then described, with code, and examined in curve tangents, directed line geometry (triangle/ray intersections, plucker coordinates, angles), and rigid body transforms/skinning (dual quaternions).

Jim Van Verth / Insomniac explains rotation formats in "Understanding Rotations", including angle (2d) Euler angles, Axis-angle, Matrix (2d/3d), complex (2d) and Quaternion (3d). Interpolation is considered for each case (including slerp).

Eric Lengyel / Terathon presents exterior (Grassmann) algebra in "Fundamentals of Grassmann Algebra". This includes the wedge product, bivectors, trivectors and multivectors. Moves on to cross product transforms, dual-basis 'anti-vectors', regressive 'antiwedge' product, and demonstrates how these can be used in homogeneous and plucker coordinate systems. This leads on to basic intersections (line, plane, point) and distances (point plane, two lines) and finally ray-triangle intersection using bivectors to avoid barycentric coordinates.

Squirrel Eiserloh / TrueThought presents "Interpolation and Splines". Takes us back to basics by looking at averaging and blending, and moves onto interpolation. Begins with quadratic and cubic Bézier curves, then moves into splines and discusses continuity. Cubic Hermite splines are up next, and how to convert between Bézier and Hermite, then Catmull-Rom splines and finishes with the more general Caridnal splines.

John O’Brien / Insomniac covers "Math for Gameplay / AI". Starts with object intersection tests (sphere-sphere, sphere-plane, AABB-AABB, AABB-ray, capsules-capsule, capsule-ray) and projecting onto a plane in a gun turrent AI example. Next up, Bayes' Theorem and conditional probability, followed by fuzzy logic.

The Web up next

Corey Clark and Daniel Montgomery present "Building a Multi-threaded Web-Based Game Engine" covering both client side (WebGL, WebSockets, etc) and server side (NodeJS, Hosting, etc).

Michael Weilbacher / Microsoft looks at server issues in "Dedicated Servers in Gears of War 3".

Michael Goddard "Developing a Javascript Game Engine"
using component based architecture. Takes an indepth look at events/promises and loading content.

Mike Dailly / YoYo investigates packing textures and command list execution for improving performance in "The Voodoo Art of Dynamic WebGL".

Marc O’Morain / Swrve takes a look at a number of issues (including iOS multitouch) in "Building Browser Based Games Using HTML5".

And rounding up everything else

Caruso Daniel explains the "Forza Motorsport Pipeline". Importing assets into the game.

GuayvJean-Francois investigates sound diffraction and absorption in "Real-time sound propagation".

Mike Lewis presents the challenges of multithreading for MMOs "Managing the Masses".

Sean Ahern looks at building better game engine tools in "It stinks and I don't like it"

Clara Fernández-Vara, Jesper Juul, and Noah Wardrip-Fruin make a case that "Game Education Needs Game History"

Chris Jurney presents his idea "Motion Blobs", a fast and crude kinect data "gesture" system, essentially an extension of the typical 2D approach to 3D. Steps are to calculate motion via background subtraction, filtering (open/close), labeling, and then correlation.

Alexander Lucas explains automated testing at Bioware in "The Automation Trap And How Bioware Engineers Quality"

Alex Mejia looks at camera movement in "Saints Row : The Third real time capture tools".

Scott Philips presents "Designing Over the Top SAINTS ROW: THE THIRD Postmortem", and highlights the importance of pre-visualization and playtesting.

Ron Pieket / Insomniac looks at eliminating downtime in "Developing Imperfect Software" via a 'Structured Binary' approach to building engine data by taking advantage of a Data Definition Language.

Benson Russell takes a look at Naughty Dog's approach to polish in "The last 10, going from good to awesome", in essence longer alpha and beta tests.

Luke Muscat takes a look at the lessons learnt while updating Fruit Ninja in "Iterating Design And Fighting Fires: Updating Fruit Ninja And Jetpack Joyride"

Tatyana Dyshlova talks about managing 300+ artists working on 500 car models in "Racing to the Finish"

Summary
Quite a collection this year, but overall seems to be less exciting content than previous years. For graphics, it seems that signed distance fields and physically based rendering is the new theme, AI is still playing catchup and character animation cycles are still a hot topic, following that theme, physics is also looking at characters and ragdolls, with destruction being the hot topic, and the web is focusing on WebGL.

Kalman filters - Part 2

In a previous post I described kalman filters and gave a very simple example of a 1D filter implemented in C. In the real world, having a filter with three or more inputs is common. A standard 3rd order example would be a filter that considers the position, velocity and acceleration of a system.

In this case you would have a 3x3 measurement matrix (H) (with the diagonal set to what you can measure), a 3x3 process noise matrix (Q), a 3x3 measurement noise matrix (R), and a 3x3 state transition matrix. The standard transition matrix would simply contain v+=a*dt, x+=v*dt, but we can use a 2nd order leapfrog integrator to get better accuracy of our position from our acceleration data. So in total we have:

You can either solve this using the standard kalman matrix maths, but you will find that an SVD approach will give you more stable solutions. OpenCV is one library that already has the kalman matrix maths built in. To set up a kalman filter with OpenCV we just need to setup the matrices we need and call the appropriate cvKalman functions.

(ugly) Code snippet follows:

#ifndef KALMAN_H
#define KALMAN_H

/** OpenCV Kalman filter for 1st, 2nd and 3rd order position, velocity and acceleration data.
    (c) Adrian Boeing
    www.adrianboeing.com
 */
#include <cv.h>

typedef enum {
    K_P = 1, //position only
    K_V = 2, //velocity only
    K_A = 4, //acceleration only
    K_PV = 3, //position & velocity
    K_PA = 5, //position & acceleration
    K_VA = 6, //velocity & acceleration
    K_PVA = 7 //position, velocity and acceleration
} Kalman_e;

class KalmanPVA {
public:
    KalmanPVA(int order, Kalman_e type, float initial_pos, float initial_vel, float initial_accel, float dt, float process_noise[3], float measurement_noise[3]) {
        m_order = order;
        //create the OpenCV kalman structures
        kalman = cvCreateKalman(order, 3, 0 );
        z = cvCreateMat(3,1,CV_32FC1);
        
        //construct the transition matrix
        SetDT(dt);
        
        //enable only the measurements we have!
        for (int i=0;i<m_order;i++) {
            cvmSet(kalman->measurement_matrix,i,i,((type&7)>>i)&1);
        }
        //configure the process noise matrix
        cvSetIdentity( kalman->process_noise_cov,    cvRealScalar(1) ); //zero the matrix 
        for (int i=0;i<m_order;i++) {
            cvmSet( kalman->process_noise_cov,i,i, process_noise[i] ); 
        }
        //configure the measurment noise matrix
        cvSetIdentity( kalman->measurement_noise_cov, cvRealScalar(1) ); //zero the matrix
        for (int i=0;i<3;i++) {
            cvmSet(kalman->measurement_noise_cov,i,i,measurement_noise[i]);
        }
        cvSetIdentity( kalman->error_cov_post,cvRealScalar(10) ); //large initial unknown
        
        //setup the initial state
        cvmSet(kalman->state_post,0,0,initial_pos );
        if (m_order>1)
            cvmSet(kalman->state_post,1,0,initial_vel);
        if (m_order>2)
            cvmSet(kalman->state_post,2,0,initial_accel);
        
        
    };
    void SetDT(float dt) {
        //3rd, 2nd, and 1st order transition matrix models
        
        //3rd order uses 2nd order leapfrog integration
        float F3[] = {   1, dt, 0.5*dt*dt,
            0,  1, dt,
            0,  0,  1};
        
        float F2[] = {   1, dt, 
            0,  1};
        
        float F1[] = {1};
        
        switch (m_order) {
            case 3:
                memcpy( kalman->transition_matrix->data.fl, F3, sizeof(F3));
                break;
            case 2:
                memcpy( kalman->transition_matrix->data.fl, F2, sizeof(F2));
                break;
            case 1:
                memcpy( kalman->transition_matrix->data.fl, F1, sizeof(F1));
                break;
        }
    }
    void SetMeasurements(float pos, float vel, float accel) {
        cvmSet(z,0,0,pos);
        //if (z->rows>1)
            cvmSet(z,1,0,vel);
        //if (z->rows>2)
            cvmSet(z,2,0,accel);
    }
    
    void Update(){
        const CvMat* y = cvKalmanPredict( kalman, 0 );
        cvKalmanCorrect( kalman, z );
    }
    float GetPositionEstimate() {
        float v = cvmGet(kalman->state_post,0,0);
        return v;
    }
    float GetVelocityEstimate() {
        if (kalman->state_post->rows>1)
            return cvmGet(kalman->state_post,1,0);
        return 0;
    }
    float GetAccelerationEstimate() {
        if (kalman->state_post->rows>2)
            return cvmGet(kalman->state_post,2,0);
        return 0;
    }
    float GetPositionError() {
        return cvmGet(kalman->error_cov_post,0,0);
    }
    CvKalman* kalman; //openCV Kalman filter
    CvMat* z; //measurement matrix
private:

    void cvKalmanNoObservation() { 
        cvCopy(kalman->error_cov_pre, kalman->error_cov_post); 
        cvCopy(kalman->state_pre, kalman->state_post); 
    } 
    int m_order;
    
};

#endif

A star (A*) pathplanning

Previously I wrote about Dijkstra, this time I'm going to look at heuristic or "AI" style approach the A* algorithm. This is a single-source single-goal pathfinding algorithm. In essence this algorithm "guesses" or estimates the distance that remain to get from a node to the destination node. These estimates guide the direction of the path (i.e. search), with each node that is visited updating the true cost of the path, until a solution is found.

A* considers what to do at each node depending on an overall heuristic score. The total score is a function of the cost of the path that has already been traversed (the "g" cost) and the estimated cost of traveling from the node to the destination (the "h" cost). There are many ways to estimate the "h" cost, typically for pathfinding problems the euclidian distance is used, however it is always better to over-estimate the remaining cost (so you will explore more feasible solutions), and so a common measure used is the manhattan distance.

The A* algorithm works on two lists, a "Closed" list of already explored nodes, and an "Open" list of nodes yet-to-be-explored.

In essence the algorithm is:

• Select a source node to operate with (eg: the start node, or the node with the lowest overall cost)
  
• For each node connected to the source node:
  
  • Calculate the already-traveresed cost
    
  • If the node is on the closed list, and the newly calculated traversed cost exceeds the existing traversed cost for the node, skip the node, otherwise remove it from closed.
    
  • Calculate the estimate cost to find the total cost, and store all these costs with the node.
    
  • Add the node to the Open list
    
• Add the source node to the Closed list
  

In Pseudo-Code (from Bryan Stout):

priorityqueue Open
 list Closed 
s.g = 0 // s is the start node 
s.h = GoalDistEstimate( s ) 
s.f = s.g + s.h 
s.parent = null 
push s on Open 
while Open is not empty 
  pop node n from Open // n has the lowest f 
  if n is a goal node 
    construct path 
    return success 
  for each successor n' of n 
    newg = n.g + cost(n,n') 
    if n' is in Open or Closed, and n'.g < = newg 
      skip 
    n'.parent = n 
    n'.g = newg 
    n'.h = GoalDistEstimate( n' ) 
    n'.f = n'.g + n'.h 
    if n' is in Closed 
      remove it from Closed 
    if n' is not yet in Open 
      push n' on Open 
  push n onto Closed 
return failure // if no path found 

A great interactive demonstration program can be found here: http://www.policyalmanac.org/games/aStarTutorial.htm And some source code for A* is available from the GrinningLizard website (The same author as the fantastic tinyXML library!) Keep in mind that waypoint-based AI navigation for games is a very mid-nineties approach, nav-mesh's are certainly the way to go! For an example of issues with waypoint approaches take a look at this video: Your best bet for navmesh research is Mikko Mononen's blog.

Using Parallel Reduction to Optimize Dijkstra

In a previous blog post I wrote about Dijkstra's algorithm. One optimisation to the naive implementation is to calculate the minimum cost node in parallel - since this can often be the most time consuming operation for large graphs. This is a perfect candidate for a parallel reduction operation.

In serial, if we wished to sum up an array (x) of numbers, the code would be:

int sum = 0;
for(int i=0; i < N; ++i) 
 sum += x[i];

However, in parallel we could speed this up by summing parts of the array and then adding those parts together. For example, to add four numbers a,b,c,d we could first add a+b and c+d in parallel, then sum ab+cd. This is called a 'reduction' operation, and it allows us to complete an otherwise O(N) operation in O(log(N)) steps.

One particularly nifty pattern is the "Butterfly" pattern of adding elements. Best described with an animation:

In code this can be done as:

for (b=n/2;b>0;b/=2) {
 for (i=0;i<n;i++) { 
         t[i] = OP(sum[i],sum[i^b]); 
 } 
 for (i=0;i<n;i++) { 
         sum[i]=t[i]; 
 } 
}

(Where 'OP' can be any serial commutative operation, e.g. add, subtract, etc.)


To see how this works assume we are adding 8 elements (n=8) then for the first step:
b is 0100
and each subsequent operation will be operating on the following elements:

0000^0100:0100
0001^0100:0101
0010^0100:0110
0011^0100:0111
0100^0100:0000
0101^0100:0001
0110^0100:0010
0111^0100:0011

or in decimal:

t[0] = [0] + [4]
t[1] = [1] + [5] 
t[2] = [2] + [6]
t[3] = [3] + [7]
t[4] = [4] + [0]
t[5] = [5] + [1]
t[6] = [6] + [2]
t[7] = [7] + [3]

Now, for the next step:
b is 0010

0000^0010:0010
0001^0010:0011
0010^0010:0000
0011^0010:0001
0100^0010:0110
0101^0010:0111
0110^0010:0100
0111^0010:0101

Again, in decimal:

t[0] = [0] + [2]
t[1] = [1] + [3]
t[2] = [2] + [0]
t[3] = [3] + [1]
t[4] = [4] + [6]
t[5] = [5] + [7]
t[6] = [6] + [4]
t[7] = [7] + [5]

Finally,
b is 0001

0000^0001:0001
0001^0001:0000
0010^0001:0011
0011^0001:0010
0100^0001:0101
0101^0001:0100
0110^0001:0111
0111^0001:0110

And in decimal:

t[0] = [0] + [1]
t[1] = [1] + [0]
t[2] = [2] + [3]
t[3] = [3] + [2]
t[4] = [4] + [5]
t[5] = [5] + [4]
t[6] = [6] + [7]
t[7] = [7] + [6]

Thus, in only 3 parallel operations we have added 8 elements together. We also have the additional bonus of having broadcast the result to all of the eight parallel adding elements - allowing future operations for each processor to use the result, or in the case of a SIMD system, to mask the appropriate processors.

Here is the code in CUDA:

int i = threadIdx.x; // Thread i holds value x_i, one thread per element
 __shared__ int sum[blocksize]; 
sum[i] = x_i; //copy to shared memory 
__syncthreads(); //synchronise 
for(int bit=blocksize/2; bit>0; bit/=2) {
    int t=OP(sum[i],sum[i^bit]);  
    __syncthreads(); 
    sum[i]=t;  
    __syncthreads(); 
}

Nifty!

Dijkstra

Dijkstra's algorithm is a well known algorithm for solving the single-source shortest path problem in a graph. It is a "global" pathfinding algorithm, in that it will find the optimal solution to all the other nodes in the graph (provided the graph has no negative cost edges or cycles). Unfortunately the algorithm is not very easy to run in parallel.

The algorithm requires two lists, a visited list, and a pending list, and an array of the current distances to each node, plus the precursors.

In essence the algorithm begins at a starting node, and all other nodes are added to a 'pending' list. The distances are initialised to 'infinite'. We begin at the starting node that has been specified.

1. From this node we calculate the distance to each neighbouring node. If the distance is less than the current record, we update the distance measurement (cost of getting to current node plus cost to get to neighbour node), and update the precursor to the neighbour node to be the present node.

2. From all nodes, the node with the minimum distance on the pending list is selected as the next node to visit. This new node is then added to the visited list and removed from the pending list..

It's quite simple, but an animation will really help get the point across:

In code, the algorithm is as follows:
(This code stores the graph in the "dist" array, with dist[i][j] representing an edge cost between i and j. If it contains '0' then there is no connection).

The initial loop sets the distances to infinity. The first inner loop finds the minimum cost node, the second inner loop updates the distances based on the current node.
The minimum distance to each node from the starting node is stored in the 'd' array. The visited array is used to record which nodes have already been processed.

void dijkstra(int s) {        
    int i, k, mini;         
    int visited[GRAPHSIZE];         
//initialize all distances to infinite, and set the visited nodes to none.              
    for (i = 1; i <= n; ++i) { 
     d[i] = INFINITY; visited[i] = 0; 
    }         
//set the distance to the starting node 's' to zero 
    d[s] = 0; 
//outer loop - iterate through all nodes 
    for (k = 1; k <= n; ++k) {                
//find the minimum cost node 
//if the node has not been visited, and the minimum cost node is not initialised, 
//then set the minimum cost node to the current node 
//if the distance to the current node is less than the current minimum,
//then update the minimum 
     mini = -1;                 
     for (i = 1; i <= n; ++i)                      
          if (!visited[i]&&((mini==-1)||(d[i]<[mini])))                                        
          mini = i; 
//now we have the current minimum node, set it to visited 
     visited[mini] = 1; //visit node               
//for each node, if there is a valid edge between the current minimum node 
//then calculate the cost to visit this neighbour node via the current node.
//if the cost is less than the current record of minimum costs, update the minimum cost.
     for (i = 1; i <= n; ++i)                               
        if (dist[mini][i]) //update distances                                    
        if (d[mini] + dist[mini][i] < d[i])                                                        
              d[i] = d[mini] + dist[mini][i];         
    } 
}

In future blog posts I'll look at the A* algorithm, a parallel reduction for Dijkstra and some SIMD/GPU/CUDA optimisations for an 8-connected Dijkstra.

Inverting Matrix - SVD (singular value decomposition)

Every once in a while you find yourself needing to solve a set of equations, or invert a matrix, or worse yet, invert a non-square matrix (eg: pseudo-inverse for manipulator inverse kinematics path control (See: Minerva IK control image on right, work I did at TUM) or kalman filtering). This is when a handy mathematical tool 'singular value decomposition' comes in to play.

In essence the SVD decomposes a matrix 'M' into three other matrices with special properties, making them easier to manipulate. The SVD equation is:

Where:

• U is a matrix whose columns are the eigenvectors of the M * M transposed matrix. (AKA left eigenvectors)
  
• S is a diagonal matrix whose diagonal elements are the singular values of M.
  
• V is a matrix whose columns are the eigenvectors of the M transposed * M matrix. (AKA right eigenvectors)
  

Discussing eigenvalues and eigenvectors I'll leave for another time, but they are quite straight-forward to calculate (get the matrix subtract a constant * I , take the determinant, equate to zero and solve, that gets you the eigenvalues!).

Now, how do you use the SVD to invert a matrix? Using a few matrix (inversion) identities (See the Matrix Cook Book, basic identity 1):

So now we just need a way to invert a diagonal matrix, which couldn't be easier. To invert a square diagonal matrix you simple invert each element of the matrix.
What if the matrix isn't square?
Well for that, we can just use the pseudo-inverse: we just invert each element of the matrix, and then transpose it. Done!


A good tutorial on the SVD can be found here by Dr. E. Garcia

Kalman Filters

(EDIT: See also Kalman Filter's Part 2)
I've been looking into a number of sensor fusion techniques, including Particle Filters and Kalman Filters (EKF,UKF, etc). I was looking for something easy to implement, easy to use, and not too computationally intensive.

I've always been a bit sceptical as to the advantages of a Kalman filter over an Alpha-Beta filter, in terms of performance gained vs cost (i.e.: difficulty to implement & use). For most simple tasks, the Alpha-Beta filter seems to work fine. (Similarly ABG vs AB).

Particle filters appeared advantageous to me, as they are relatively intuitive and aren't that hard to use, and easy to extend. The difficulty in implementing them can be overcome by using an existing library,

Kalman filters in my opinion have as a major disadvantage of being difficult to use (you need to construct multiple models), but are relatively easy to implement and computationally efficient. Thus, the use of an existing library doesn't help much.

I spoke to two experts in the field, Julius Ziegler (who worked on the DARPA Grand Challenge Annieway vehicle with me) and Rob Reid (a local, working on visual SLAM with Kalman filters). Both suggested that Kalman filters would be the most appropriate choice. (Reading various papers seems to indicate a merged (E)Kalman & Particle filter approach is the winner)

Wikipedia provides an overview of Kalman filters, but the real problem is in understanding what all the symbols actually mean, and how it works.

Like alpha-beta, Kalman filters are prediction-correction filters. That is, they make a prediction, then correct it to provide the final estimate of the systems state. Alpha-Beta's don't have a natural extension to include control inputs, or system identification techniques. They are also (effectively) constrained to linear position/velocity tracking. Kalman filters can be extended to handle non-linear conditions.

In essence a Kalman filter will (vaguely):

1. Estimate the state (from previous)
   
2. Determine the error
   
3. Predict the next state
   
4. Predict the next error
   
5. Measure the system
   
6. Update the state calculations
   
7. Update the error calculations
   

Bellow are the equations that represent this process, followed by the definitions of each symbol:

x_temp_est = A * x_est(last) + B * u_control_input
P_temp = A * P(last) * A^T + Q
K = P_temp * H^T * (H* P_temp * H^T + R)^-1
x_est = x_temp_est + K * (z_measured - H * x_temp_est) 
P = (I - K*H) * P_temp

• x - is the state we want to estimate (e.g. position)
  
• u - is the input to the control system (if it exists) (e.g. a force)
  
• A - state (model) matrix, how to get from last state, to next state (e.g.: position += velocity *dt). (Some people (eg: engineers) use an alternative notation 'F' for this matrix)
  
• B - is the input matrix (Some people use an alternative notation 'G' for this matrix)
  
• H - is the measurement matrix (what are we measuring?) e.g.: position or velocity, etc.
  
• Q - is the process noise covariance matrix (you just 'know' this)
  
• R - is the measurement noise covariance matrix (likewise, you 'know' this, usually from measurements, data sheets, etc, ie: how noisy are your sensors?)
  
• P - the estimate error covariance matrix (calculated)
  
• K - the kalman gain (calculated)
  

The things you must 'know' are:

• How your system moves from one state to the next (e.g. given a current position, velocity, change in time, what will the next state be?) => This is the 'A' or 'F' matrix
  
• What are you able to measure? (Maybe just acceleration, or just position?) => This is the H matrix
  
• The amount of noise in your measurements and the state-transitions (e.g. the standard deviation of the signal noise, and how 'wrong' your simplified model of the state-transitions are) => These are Q and R matrices
  

Completely confused? Good. Some examples will help clear this whole mess up. First, a simple example where we are observing a noisy-constant. (e.g. measuring the height of a water-tank with a poor quality sensor) We have no control input so we can ignore the B and u terms. We have no state-transition model (it is a constant!) so we can set A to identity (i.e. one). Our first equation becomes:

x_temp_est = x_est(last)

Subsequently the second equation becomes:

P_temp = P(last) + Q

Where Q is our estimation of the variance in our state transition. Since we are measuring a constant, this value should be small. (Remember the real world always has noise, so assuming zero noise is always a bad idea) We measure the state directly (i.e. the noisy value of the constant) so the H matrix is also identity. This means:

K = P_temp * (P_temp + R)^-1
x_est = x_temp_est + K * (z_measured - x_temp_est) 
P = (1- K) * P_temp

Again, we keep R since that is our all-important measurement noise estimate. Now let's see this in action in C:

/** A simple kalman filter example by Adrian Boeing  
 www.adrianboeing.com  
 */  

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 

double frand() { 
    return 2*((rand()/(double)RAND_MAX) - 0.5); 
} 

int main() { 

    //initial values for the kalman filter 
    float x_est_last = 0;
    float P_last = 0; 
    //the noise in the system 
    float Q = 0.022;
    float R = 0.617;
     
    float K; 
    float P; 
    float P_temp; 
    float x_temp_est; 
    float x_est; 
    float z_measured; //the 'noisy' value we measured
    float z_real = 0.5; //the ideal value we wish to measure
     
    srand(0); 
     
    //initialize with a measurement 
    x_est_last = z_real + frand()*0.09;
     
    float sum_error_kalman = 0;
    float sum_error_measure = 0;
     
    for (int i=0;i<30;i++) {
        //do a prediction 
        x_temp_est = x_est_last; 
        P_temp = P_last + Q; 
        //calculate the Kalman gain 
        K = P_temp * (1.0/(P_temp + R));
        //measure 
        z_measured = z_real + frand()*0.09; //the real measurement plus noise
        //correct 
        x_est = x_temp_est + K * (z_measured - x_temp_est);  
        P = (1- K) * P_temp; 
        //we have our new system 
         
        printf("Ideal    position: %6.3f \n",z_real); 
        printf("Mesaured position: %6.3f [diff:%.3f]\n",z_measured,fabs(z_real-z_measured)); 
        printf("Kalman   position: %6.3f [diff:%.3f]\n",x_est,fabs(z_real - x_est)); 
         
        sum_error_kalman += fabs(z_real - x_est); 
        sum_error_measure += fabs(z_real-z_measured); 
         
        //update our last's 
        P_last = P; 
        x_est_last = x_est; 
    } 
     
    printf("Total error if using raw measured:  %f\n",sum_error_measure);
    printf("Total error if using kalman filter: %f\n",sum_error_kalman);
    printf("Reduction in error: %d%% \n",100-(int)((sum_error_kalman/sum_error_measure)*100));
     
     
    return 0; 
}

This should generate output similar to:

Mesaured position:  0.562 [diff:0.062]
Kalman   position:  0.498 [diff:0.002]
Ideal    position:  0.500 
Mesaured position:  0.563 [diff:0.063]
Kalman   position:  0.509 [diff:0.009]
Total error if using raw measured:  1.566003
Total error if using kalman filter: 0.441561
Reduction in error: 72% 

More (non-trivial) examples will follow soon. In the meantime I can suggest the following resources:

• http://www.cs.unc.edu/~tracker/ref/s2001/kalman/index.html
  
• http://www.cs.cmu.edu/~biorobotics/papers/sbp_papers/integrated3/kleeman_kalman_basics.pdf
  
• http://www.acfr.usyd.edu.au/education/multiSensorDataFusion.shtml
  
• http://www.ee.pdx.edu/~greenwd/kalman_overview.pdf
  

Source code/libraries:

• http://kalman.sourceforge.net/
  
• http://bayesclasses.sourceforge.net/
  
• http://www.memsense.com/index.php/Product-Pages/kalman-filter-library.html
  
• http://sites.google.com/site/jordiuavs/Home/kalman_filter_by_Jordi.txt
  
• http://www.orocos.org/bfl
  
• http://www.scipy.org/Cookbook/KalmanFiltering
  

Catchup Post

A short collection of interesting links/articles recently:

• Designing Fast Cross-Platform SIMD Vector Libraries plus source code, an interesting read on rolling your own vector libraries. If your looking for a bigger library you can try WildMagic or the Sony maths libraries.
  
• Compiler Cache helps reduce your compile times.
  
• Howard Yeend has an excellent tutorial on K-means clustering (ie: grouping data together, given you know how many groups you want - always a shame..)
  
• New DirectX SDK is out, not sure how this is exciting..
  
• Apple released the iPad, I have no idea if this is good. But I like the alternative, Linux-friendly iPad.
  
• I stumbled upon Alexandre Mutel's blog on 4k programming
  
• I discovered the IEEE has RSS feeds for most journals.
  
• nVidia have named the GF100 chips: GeForce GTX 480 and GeForce GTX 470.
  
• Anandtech have an
  article on GF100 explaining how the CUDA centric card is also great for gaming.
  
• Facebook released HipHop a PHP to C/gcc conversion tool.
  
• Clang self-hosts! I wonder how long until it becomes a serious g++ competitor.
  
• Symbian OS will be open sourced, perhaps a little too late? Still, Nokia are to be admired.
  
• An interesting gallery of processor cache behaviour, unfortunately I don't think the author is fully qualified on this topic. Still, interesting pictures.
  
• The PS3 hypervisor was hacked
  

Finally, a video of the top Tron AI from the Google AI Challenge:

Catchup post

It's been a while since I updated, reason being it is exams/assignments marking period, and I had two GPU industry projects due (3D groundwater flow fluid simulation and a pathfinding/travelling salesman with profits project). 

The biggest news item was that the ACM decided to start (over)enforcing its rules saying that you can not link to preprint and author pages. Thankfully, it started a call-to-arms and prominent pages like Ke-Sen Huang's SIGGRAPH links have been restored. I wonder how many less public pages have silently slipped away. Frankly, I can't wait until the concept of conferences and journals disappear. My websites have always had far more impact that my publications, and it can't be long until the same can be said universally.

A short update with some interesting things in the last while:

• Obama has made a comment about robotics, it's great to hear the field is getting a lot more visible and more support, but I still feel we have a very long way to go to break the western cultures distrust of robotics.
• Sony/Toshiba/IBM CELL is now officially dead.
• There is a new GPU computing community website, just started up: http://gpucomputing.net/
• nVidia's raytracing technology Optix is finally available, but only for next-gen GPU's. The "RealityServer" was also released, to great disappointment to many.
• Tech report has a more detailed report on nVidias new Fermi based GPU's.
• Caustic claims to be selling raytracing accelerators, but they haven't been able to sell me one yet. Smells like vapourware.
• Chromium OS was released. Not sure who wants this? Fast starts haven't been an issue for me since laptop's could sleep. I guess they are aiming for very low prices.
• Peteris Krumins posted a great blog post on the highlights of MIT's algorithms lectures
• Wolfire did an interesting post on generating vertex weights.
• Felix von Leitner put together a very interesting presentation on compiler optimization tricks. It's good proof that compilers are much better than humans for most optimizations, plenty of new things I learnt from this.
• Finally, I saw a link to this paper:Aggregate Dynamics for Dense Crowd Simulation. I've not had a chance to read it yet, but the results look fantastic and the idea is novel. (Using fluid-simulation techniques for representing crowds.) Video below:

   

Nearest Neighbor

A bunch of links on the nearest-neighbor problem, for higher dimensions:
http://www.mit.edu/~andoni/LSH/
http://www.cs.umd.edu/~mount/ANN/
http://www.cs.sunysb.edu/~algorith/implement/ranger/implement.shtml

And, for good measure a set of comparisons of optical flow algorithms:
http://vision.middlebury.edu/flow/eval/

I'll probably stick to the OpenCV defaults anyway, but its nice to know there are options. It would seem Bruhn et al. is the most accurate..

Talking Piano and AI

Daniel Wedge sent me this interesting link on a talking piano!


A great idea, I wonder if it had been done before..

In other news, the International Joint Conference on Artificial Intelligence archive has been made open to the public covering all the way from 1969 - 2007!

AI space

AI space has a nice set of introductory javascript applets to a number of AI learning techniques. It covers graph searching, Belief and Decision Networks, Decision Trees, Constraint satisfaction problems, Neural Networks and more. It is a more 'traditional' AI site.

I also came across this great SLAM visualization tool for Player. A valuable resource on this is the Robotic Mapping and Exploration book (Springer Tracts in Advanced Robotics).

MAGIC 2010

Do you remember the DARAP Grand Challenge? Well now the Australian DSTO is holding its own challenge, the Multi Autonomous Ground-robotic International Challenge (MAGIC 2010), with a total of $1.6 million US in prize money. The goal is to use a multi-robot team to perform "intelligence, surveillance and reconnaissance mission in a dynamic urban environment."

The MAGIC requires entrants to complete the following tasks:
(i) Accurately and completely explore and map the challenge area;
(ii) correctly locate, classify and recognise all simulated threats; and
(iii) complete all phases within 3.5 hours.

The final event is scheduled to take place during the week of November 8, 2010, somewhere in South Australia.

I'm putting together a team from WA, if you are interested let me know.
The participants conference is in Adelaide next week. See you there!

Paris Game AI

I stumbled across the Paris Game AI conference, there was some interesting stuff.

The most interesting was Mikko's open source automatic navmesh generator. The demo makes the technology look useable enough.

If your not convinced by the navmesh approach yet, the AI blog has an excellent discussion of thebenefits of a navmesh over waypoints. Including some amusing failed videogame pathfinding algorithms.

The most interesting presentations were Coordinating Agents with Behavior Trees for Farcry from Crytech. It's a different approach to behavior management than the standard emergent behavior approach using a blackboard. I'm not convinced its really novel, all they have done is named the approach taken by your standard undergrad robotics students, but they make a good argument for it. The big advantage is that you have far more control over how the behaviors are activated, which is something you really need in games if you have a demanding director.

Another great presentation is Killzone 2 Multiplayer Bots. Not because there is anything amazingly new in it, but just because they cover their entire process, squads, 'strategic reasoning' by annotating special AI nodes with certain behavior, the waypoint network, and influence maps.

While I am on the topic of game AI, the botprize will be having it's first round playoffs at ECU in Perth next week, should be fun!

AdaBoost

There are a couple of algorithms that float around that I always feel I need to know more about. During a number of data mining conferences the top 10 algorithms in data mining were collated.

One that seems to get a lot of attention in computer vision and robotics is Adaptive Boosting (aka AdaBoost). You can use this algorithm to create a strong classifier from a set of weak classifiers. (One way to think of this is that weak classifiers give you a very approximate answer to a classification problem, strong ones give a much much better approximation)

For example in the figures below, you can see how the algorithm builds up a strong classifier from combining multiple weak classifiers.
(from Jan Sochman's AdaBoost talk)











AdaBoost combines these classifiers through a set of weights that it adaptively adjusts.

I can hardly claim to be an expert on the subject, but since I find the mathsy description of it on Wikipedia confusing, here is my summary of it in CS-algorithmish talk:

Start with a training data set of features, labeled with the correct answer. 
ie: data[i], correct_answer[i]
Run each weak classifier with the data set and store its results.
ie:results[x,i] = classifierX(data[i])
Initialize each of the weightings with 1 / (number of entries in the data set).
At this step we are starting of with each classifier having equal weighting.
ie: D[i] = 1/m

For N number of iterations:
 For each classifier:
  For each feature, i: 
   if the output of the classifier does not equal the correct answer, update the error. 
   error+=D[i] 
   (ie: sum of the weights of the misclassified samples)
   if the error is less than the minimum error, mark this as the best classifier
  End For
 Now update the importance of the classifier: 
 alpha[best classifier] = 0.5 * log(1 - minimum error)/minimum error 
 And update all the classifier weightings:
 For each feature,i
  D[i] = *= exp(-alpha[best] * answer[i] * result(best,i)
  z+=D[i] 
 End for
 For each feature, update D[i] = D[i]/z to renormalize.

And.. were done!

In my experiments it seems not to work too well,.. it definitely improves on your weak classifiers a bit, but they have to be pretty close to spot on or able to solve the problem already. I think I'd still prefer using a neural net for classification..