Pixelux Digital Molecular Matter engine - Free license soon!

AMD have just announced there will be a free version of the DMM engine. It sounds like it will be closely integrated with Bullet Physics and provide OpenCL support.

I'm looking forward to this as DMM is the only realtime finite-element (FEM) based physics engine around. I'm interested to see how well it works.

Hopefully I can get a copy soon and some time to integrate it into the Physics Abstraction Layer.

Here is a video of DMM to refresh your memory:

Alpha Beta Filters

(EDIT: See also Kalman filters)

Alpha-Beta filters are fairly well explained on wikipedia, and are generally a very easy first-stop solution before heading to Kalman or particle-filter alternatives. (The next obvious step is to extend this to include acceleration (Alpha-Beta-Gamma), and to limit the error/estimates to a sensible range for your application).

Here is a small example program showing how they work, it generates output similar to this:

Ideal     position: -0.897 -0.443
Mesaured  position: -0.890 -0.421 [diff:0.029]
AlphaBeta position: -0.898 -0.442 [diff:0.001]
Total error if using raw measured: 1.522438
Total error if using a-b filter:   1.059981
Reduction in error: 69% 

C source code follows:

/** A simple alpha-beta filter example by Adrian Boeing 
    www.adrianboeing.com 
 */ 

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

typedef struct { 
    float alpha; //alpha value (effects x, eg pos)
    float beta; //beta value (effects v, eg vel)
    float xk_1; //current x-estimate 
    float vk_1; //current v-estimate 
} AlphaBeta; 

void InitializeAlphaBeta(float x_measured, float alpha, float beta, AlphaBeta* pab) {
    pab->xk_1 = x_measured; 
    pab->vk_1 = 0; 
    pab->alpha = alpha; 
    pab->beta = beta; 
}    

void AlphaBetaFilter(float x_measured, float dt, AlphaBeta* pab) {
    float xk_1 = pab->xk_1; 
    float vk_1 = pab->vk_1; 
    float alpha = pab->alpha; 
    float beta = pab->beta; 
    
    float xk; //current system state (ie: position)
    float vk; //derivative of system state (ie: velocity)
    float rk; //residual error 
     
    //update our (estimated) state 'x' from the system (ie pos = pos + vel (last).dt)
    xk = xk_1 + dt * vk_1; 
    //update (estimated) velocity  
    vk = vk_1; 
    //what is our residual error (mesured - estimated) 
    rk = x_measured - xk;  
    //update our estimates given the residual error. 
    xk = xk + alpha * rk; 
    vk = vk + beta/dt * rk; 
    //finished! 
     
    //now all our "currents" become our "olds" for next time 
    pab->vk_1 = vk; 
    pab->xk_1 = xk; 
} 

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

int main(int argc, char *argv[]) {
    AlphaBeta ab_x; 
    AlphaBeta ab_y; 
    double t; //time 
    double x,y; //ideal x-y coordinates 
    double xm,ym; //measured x-y coordinates 
    double xnoise = 0; //noise we are inserting into our system
    double ynoise = 0; 
    double m_error = 0; //total error (measured vs ideal)
    double ab_error = 0; //total error (ab filter vs ideal)
#define DT 0.1
    //intialize the AB filters 
    InitializeAlphaBeta(1,0.85,0.001,&ab_x); //x position
    InitializeAlphaBeta(0,1.27,0.009,&ab_y); //y position
    srand(0); 

    for (t = 0; t < 4; t+=DT) {
        //our 'true' position (A circle) 
        x = cos(t); 
        y = sin(t); 
        //update our simulated noise & drift 
        xnoise += frand()*0.01;
        ynoise += frand()*0.01;
        //our 'measured' position (has some noise) 
        xm = x + xnoise; 
        ym = y + ynoise; 
        //our 'filtered' position (removes some noise) 
        AlphaBetaFilter(xm,DT, &ab_x); 
        AlphaBetaFilter(ym,DT, &ab_y); 
         
        //print  
        printf("Ideal     position: %6.3f %6.3f\n",x,y); 
        printf("Mesaured  position: %6.3f %6.3f [diff:%.3f]\n",xm,ym,fabs(x-xm) + fabs(y-ym)); 
        printf("AlphaBeta position: %6.3f %6.3f [diff:%.3f]\n",ab_x.xk_1,ab_y.xk_1,fabs(x-ab_x.xk_1) + fabs(y-ab_y.xk_1)); 
         
        //update error sum (for statistics only) 
        m_error += fabs(x-xm) + fabs(y-ym); 
        ab_error += fabs(x-ab_x.xk_1) + fabs(y-ab_y.xk_1); 
    } 
    printf("Total error if using raw measured: %f\n",m_error);
    printf("Total error if using a-b filter:   %f\n",ab_error); 
    printf("Reduction in error: %d%% \n",(int)((ab_error/m_error)*100));
    return 0; 
}

Catchup: Graphics Links Post

Again, a set of interesting links from the last few weeks:

• Simple and Fast Multimedia Library, a SDL replacement with a few more modern concepts. (eg: shaders!)
  
• The OpenCL debugger, gDEBugger CL is in public beta for Windows, Mac OS X and Linux.
  
• ATI Stream Software Development Kit (SDK) v2.01 now supports DirectX / OpenCL interoperability, plus debugging with GDB! A giant leap forward for OpenCL by AMD!
  
• Sander Van Rossen has put together a set of interesting posts on implementing a constructive solid geometry (CSG) system.
  
• OpenRL (raytracing language) has gone beta. I've not had positive experiences with Caustic so far. Wait and see I guess.
  
• Apparently, someone managed a 20x speedup for SQL select queries with the GPU. I guess I need to eat my words. I'm still not convinced that this is really a CPU-intensive issue. (Maybe for scheduling and retrieving flight prices - a dijkstra-like issue)
  
• An interesting article on the history of photoshop, never realised how humble the beginnings were, and how much influence one mans decision to allow optional extras (ie: plug-ins) played on the future of the product.
  
• I never realised this, butAutodesk has a massive price reduction for students, this includes 3ds Max, Maya, and AutoCAD. Wow!
  

And, as per usual, I'll finish off the catchup post with some neat videos. First, one that I think does a great job of explaining the animation process, and then a neat video showing just how far CG has come in film. (I recall speaking to Paul Debevec about Spiderman, and he said the director decided to redo all the actors in all scenes in CG since they looked better and they had more control over facial expressions,etc.)
Cirkus Animation ABC: Stargate studios reel:

Profiling on Mac OSX with Saturn

It is always advisable to profile your code before trying to optimize it. It is often most helpful to role your own (using timers & counters), but often a quick and simple tool will help a lot.

For windows I like the very sleepy profiler. On Mac OSX Apple provide the Saturn profiler (See: Saturn profiler user guide).

To demonstrate we can try compiling a small program:

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

#include <saturn.h> //include the saturn profiler

double func1(double x) { //do some maths
    return sin(x)*cos(x); 
} 

double func2(double x) { //do some more maths
    return pow(x,3); 
} 

int main() { 
    double x=0; 
    double z=0; 
    startSaturn(); //being profiling 
    for (x=0;x<100;x+=0.01) {
        z+=func1(x)+func2(x)*tan(x); 
    } 
    stopSaturn(); //end profiling 
    printf("z:%f\n",z); //make sure compiler doesn't throw our computations away
    return 0; 
} 

Now we need to compile it with profiling support:

g++ x.c -finstrument-functions -lSaturn -m32 -O2

Some common problems include:

• 
  
  
  

  Undefined symbols: ___cyg_profile_func_enter
  

  This comes from the -finstrument-functions option : it requires a special hook for each function - this is provided by saturn, so you must link with it.
  
• 
  
  
  

  ld: warning: in /usr/lib/libSaturn.dylib, file is not of required architecture
  Undefined symbols: _startSaturn
  

  Again, either you forgot -lSaturn, or you have remembered it, but are using a 64 bit OS/chip. Specify '-m32' to force 32bit mode.
  

If it all compiled succesfully, then start Saturn and choose 'Saturn', 'Launch Process'. Then select the executable (eg: a.out) in the dialog box.

Press OK, and Saturn will run and profile your program, and generate a folder with the profiler output data. (eg: Saturn_profile_a.out.000)

You can then select the data file to view the output, and Saturn will display a call graph and the amount of time spent in each function. Now the fun of optimizing can begin. Enjoy!

Intersection of a Convex Hull with a line segment

A ray, or line segment can be represented parametrically as:
S(t) = A + t(B-A)
Where A and B are the endpoints of the segment, and t is the parameter that ranges from –infinity to +infinity for a ray, or just 0..1 for a segment.

A plane can be represented as n.X = d, where n is the plane’s normal, and d is the offset. (Given the plane’s normal, and a single point on the plane, P, we can calculate: d = -n.P)

A convex object can be represented as the area contained within a set of planes. Thus, to find the intersection between a line segment and a convex object, we just need to clip it against all the planes that form the convex object.

First, substitute the line equation into the plane, and solve for t:
i.e.:
n.(A+t(B-A)) = d
n.A + t*n.(B-A) = d
note: using identity ru.sv = rs(u.v), where u,v are vectors, and r,s are scalars
re-arrange to solve for t, the intersection point along the line:

We can determine if the plane faces the segment or not by evaluating the dot product of the plane’s normal, and the line segment’s direction vector.

From this we can determine which point of an intersecting line segment to influence. If the plane is facing the segment direction, then we can clip against the end point, otherwise we can clip against the start point.

As we are testing the intersection against a convex object we can simply keep clipping against each plane and altering the segment endpoints until we have the minimum remaining line length, or the intersection length becomes empty (there is no intersection).

In pseudo-code the entire operation is:

AB = B – A
tFirst = 0
tLast = 0
for all planes:
 denom = N dot AB
 dist = d – N dot A
 t = dist/denom
 if (denom>0 )
  if (t>tFirst) tFirst = t;
 else
  if (ttLast)
  No Intersection

Solving Linear Systems

In school you probably learnt how to solve systems of linear equations with techniques like Gaussian elimination, and Row-Reduced Echelon Form (RREF). However a simple, brute-force way to solve linear systems on a computer is through iteration.

Say we wish to solve the following equations:

4x -  y +  z = 7
 4x - 8y +  z =-21
-2x +  y + 5z = 15

Then we can re-write them as:

We can solve these with Gauss-Seidel iteration just by plugging in the current x,y,z values we calculate from these equations (and begining with an initial estimate.)

Thus, the Gauss-Siedel method in C-code looks something like:

#include <stdio.h> 

int main() { 
//a sparse way of representing the equations 
float eq[3][4];
eq[0][0] = 7/4.0; eq[0][1] = 0; eq[0][2] = 1/4.0; eq[0][3]= -1/4.0;
eq[1][0] = 21/8.0; eq[1][1] = 4/8.0; eq[1][2] = 0; eq[1][3]= 1/8.0;
eq[2][0] = 15/5.0; eq[2][1] = 2/5.0; eq[2][2] = -1/5.0; eq[2][3]= 0;

float x,y,z; 
x=1;y=1;z=2; //initial guess

//10 iterations of gauss-seidel 
for (int i=0;i < 10;i++) {
  x = eq[0][0] + eq[0][2]*y + eq[0][3]*z;
  y = eq[1][0] + eq[1][1]*x + eq[1][3]*z;
  z = eq[2][0] + eq[2][1]*x + eq[2][2]*y;
  printf("%f %f %f\n",x,y,z); 
} 

return 0; 
}

Producing this output:

1.500000 3.625000 2.875000
1.937500 3.953125 2.984375
1.992188 3.994141 2.998047
1.999023 3.999268 2.999756
1.999878 3.999908 2.999969
1.999985 3.999989 2.999996
1.999998 3.999999 3.000000
2.000000 4.000000 3.000000
2.000000 4.000000 3.000000
2.000000 4.000000 3.000000

Converging towards the solution nicely. (Things don't always converge, only when Ax=B, A is diagonally dominant - but that is another story)

Jacobi iteration does not converge as quickly, but is easy to execute in parallel. With Jacobi iteration you simply use the last iterations x,y,z value instead of updating it.

See? Solving systems of equations is easy.

Building MRPT under Windows with MSVC

These are brief instructions on building MRPT under Windows.

1. Download CMake, and install it. (You need a recent version, > 2.6)
   
2. Download OpenCV, and install it. (I used the source package)
   
3. Download wxWidgets, and install it. (Again, I used source.) [Note: wxWidgets is in my opinion one of the worst packages I've ever had to deal with - its build is notoriously unstable]
   
4. Find the wxWidgets 'include/wx/msw/setup.h' file, and enable the use of OpenGL widgets:
   
   

   #define wxUSE_GLCANVAS       1

5. Open (and convert) the VC6 build solution (build/msw: wx.dsw), and perform a batch build (select all, and go!). You will need to do this three or four times. (wx does not have the dependencies correct!)
   
6. Close Visual Studio, and add an environment variable (right click on my computer, properties, advanced). Set wxWidgets_ROOT_DIR to C:\wxWidgets-2.8.10\lib\vc_lib (or your equivalent).
   
   
   
   We are now done with wxWidgets - you may wish to try some samples to make sure it works. (Pick ones that use OpenGL)
   
   
7. Now run CMake with OpenCV. Set the Source and Bin directories to your OpenCV directory, eg:
   
   

   source: C:/OpenCV2.0
   bin: C:/OpenCV2.0
   

   . Press Configure a few times and generate the solution files.
   
8. Compile OpenCV from the generated solution (Should be in C:/OpenCV2.0). Try some of the examples. Many of them require a webcam, so try plugging one to make sure DirectShow works.
   
   
9. Set your visual studio library directory (Projects and Solutions/ VC++ Directories) to include the OpenCV lib directory. Close visual studio.
   

10. Run the MRPT makefile 'win32_rebuild_MSVC9_GUI.bat' (not sure why?), and then start CMake again, clear the cache, and point CMake to MRPT. eg:
    
    
    
    
    

    source: C:/lib/MRPT/mrpt-0.8.0
    binaries: C:/lib/MRPT/mrpt-0.8.0/makefiles/MSVC9
    

11. Build MRPT, all done!
    

Before running the examples you will probably need to copy all openCV dll's from the OpenCV bin directory, and copy all the wx widgets dll's from lib/vc_dll directory into the MRPT bin directories.