Karsten Østergaard Noe – Visual Computing Lab

Oriented Particles in 2D

Karsten Østergaard Noe — Wed, 02 Jan 2013 09:02:47 +0000

We have been asked quite a lot about the Oriented Particles approach to physically based simulation, where particles besides position and velocity also have an ellipsoid shape, an orientation and an angular velocity. So we decided to do a small c++ demo in 2D of how it can be done in practice.

Basically what the demo show is a somewhat simplified version of our Oriented Particles Christmas Card. The demo is a standard GLUT/OpenGL application but it can also be cross-compiled to javascript/WebGL using emscripten, the result of which can be seen here. In the following we will assume that the reader has read the original paper and thus we will just give a brief introduction to the demo code.

In the demo we want to be able to toss this guy around:

Also we want to make it look like the guy has bones in his body. For this purpose we have drawn a “skeleton” and used OpenCV to identify individual bones and fit ellipses to these. The result can be seen below with the ellipses shown in red. Our little OpenCV program outputs code directly which has been inserted into the example code.

Because of some technicalities with emscripten we wanted to avoid loading files from the example, so the shaders have been inlined as strings in c-headers. Also the image of the guy has been saved to a c-header using GIMP and included directly in the example. Source code for the example can be found here.

The class PositionBasedDynamics encapsulates a basic generalized Position Based Dynamics simulation loop, and two constraints, StayAboveLineConstraint and GeneralizedShapeMatchingConstraint, have also been included in the project. In main.cpp particles are attached to nearby particles using implicit (generalized) shape matching constraints in the function CreateObject and the object is inserted to the physics system in the function UploadOrientedParticlesObjectToPhysicsSystem.

For visualization a grid of vertices is created covering the area of the simulated particles. The vertices are triangulated and supplied with texture coordinated and used for rendering. As the particles each have an orientation they each constitute a two-dimensional coordinate system and can thus be used for skinning the grid mesh. For this purpose a parametrisation of vertices is made in the function GenerateGrid2D. Each vertex can be skinned from up to 3 particles and the parametrisation information is stored as 3D texture coordinates where the integer part of each component describes which particle to skin from and the fractional part describes the weighting to use for this particle. Skinning matrices are computed and uploaded to a float RGBA texture and the actual skinning is done in the vertex shader.

To get some motion in the system the little guys bag is moved around the scene. The example is not optimal – for instance the grid mesh could be based on index buffers which would save some calculations. Another problem with the 2D example is that for shape matching involving only 2 particles we sometimes find an optimal rotation that mirrors the particle. A quick fix to remedy this would be to check for mirroring matrices in the 2D matrix class. Also as you might know there are some issues with floating point textures in WebGL – especially when we want to access the texture from vertex shader. So to improve compatability for WebGL/GLES it might be a good idea to do the skinning on the CPU instead of in the vertex program – this was the solution we used when porting the system to iOS.

Hope you’ll have fun with the example!

Oriented Particles in 2D – an example.

Karsten Østergaard Noe — Wed, 02 Jan 2013 09:02:47 +0000

In the demo we want to be able to toss this guy around:

Hope you’ll have fun with the example!

Oriented Particles Christmas Card

Karsten Østergaard Noe — Fri, 07 Dec 2012 11:29:51 +0000

To warm up for Christmas we decided to create a small interactive Christmas card for all of you.

The interactive part is based on a 2D version of the Oriented Particles approach to realtime physics simulation.

The card was implemented using our Alexandra Physics Engine in c++ and cross-compiled to javascript using emscripten.

The card has been reported to work on most systems in Chrome and Firefox but some systems lack the required WebGL capabilities.

You can try if it works on your system here.

There will be a tutorial on the 2D oriented particles later on this blog.

Fun with oriented particles (do knot try this at home).

Karsten Østergaard Noe — Tue, 15 May 2012 08:41:46 +0000

Hi everybody,

Recently we have spent some more time on real-time physics simulation.

One of the things we have been working on can be seen in this video of a real-time volumetric simulation of an elastic torus knot.

The simulation was realized using the ‘oriented particles’ approach to shape matching developed by Matthias Müller and Nuttapong Chentanez.

The implementation will be used in the next version of the Visible Ear Simulator ( ves.cg.alexandra.dk ).

Denmark rendered as a jigsaw puzzle :-)

Karsten Østergaard Noe — Tue, 05 Apr 2011 09:27:43 +0000

Recently the Danish Ministry of the Interior and Health published a lot of key figures about Danish municipalities on a new web page ( http://www.noegletal.dk ). We put a little effort into making a pipeline capable of turning all those numbers into nice graphics. The figure seen above is an example of what our software can do. It is basically a bar diagram where each bar (or ‘brick’) is a municipality and the bar height and color encodes a quantity – here the population figure is shown.

Noe’s tutorial on deforming 3D geometry using RBFs

Karsten Østergaard Noe — Fri, 14 Aug 2009 13:01:50 +0000

In the following I will present a method for deforming three dimensional geometry using a technique relying on radial basis functions (RBFs). These are mathematical functions that take a real number as input argument and return a real number. RBFs can be used for creating a smooth interpolation between values known only at a discrete set of positions. The term radial is used because the input argument given is typically computed as the distance between a fixed position in 3D space and another position at which we would like to evaluate a certain quantity.

The tutorial will give a short introduction to the linear algebra involved. However the source code contains a working implementation of the technique which may be used as a black box.

Source code with Visual Studio 2005 solution can be found here. The code should also compile on other platforms.

Interpolation by radial basis functions

Assume that the value of a scalar valued function $$F : mathbb{R}^3 rightarrow mathbb{R}$$ is known in $$M$$ distinct discrete points $$mathbf{x}_i$$ in three dimensional space. Then RBFs provide a means for creating a smooth interpolation function of $$F$$ in the whole domain of $$mathbb{R}^3$$. This function is written as a sum of $$M$$ evaluations of a radial basis function $$g(r_i) : mathbb{R} rightarrow mathbb{R}$$ where $$r_i$$ is the distance between the point $$mathbf{x} = (x, y, z)$$ to be evaluated and $$mathbf{x}_i$$:

$$!F(mathbf{x}) = sum_{i=1}^M a_i g(||mathbf{x} – mathbf{x}_i||) + c_0 + c_1 x + c_2 y + c_3 z, mathbf{x} = (x,y,z) mathbf{(1)}$$

Here $$a_i$$ are scalar coefficients and the last four terms constitute a first degree polynomial with coefficients $$c_0$$ to $$c_3$$. These terms describe an affine transformation which cannot be realised by the radial basis functions alone. From the $$M$$ known function values $$F( x_i, y_i, z_i ) = F_i$$ we can assemble a system of $$M+4$$ linear equations: $$mathbf{G} mathbf{A} = mathbf{F}$$
where $$mathbf{F} = (F_1, F_2, ldots, F_M, 0, 0, 0, 0)$$, $$mathbf{A} = (a_1, a_2, ldots, a_M, c_0, c_1, c_2, c_3)$$ and $$mathbf{G}$$ is an $$(M+4) times (M+4)$$ matrix :

$$! mathbf{G} = left[begin{array}{cccccccccc}g_{11} & g_{12} & bullet & bullet & bullet & g_{1M} & 1 & x_1 & y_1 & z_1 \ g_{21} & g_{22} & bullet & bullet & bullet & g_{2M} & 1 & x_2 & y_2 & z_2 \ bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet \ bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet \ bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet & bullet \ g_{M1} & g_{M2} & bullet & bullet & bullet & g_{MM} & 1 & x_M & y_M & z_M \ 1 & 1 & bullet & bullet & bullet & 1 & 0 & 0 & 0 & 0 \ x_1 & x_2 & bullet & bullet & bullet & x_M & 0 & 0 & 0 & 0 \ y_1 & y_2 & bullet & bullet & bullet & y_M & 0 & 0 & 0 & 0 \ z_1 & z_2 & bullet & bullet & bullet & z_M & 0 & 0 & 0 & 0end{array}right] $$

Here $$g_{ij} = g(|| mathbf{x}_i – mathbf{x}_j ||)$$. A number of choices for $$g$$ will result in a unique solution of the system. In this tutorial we use the shifted log function:

$$! g(t) = sqrt{log(t^2 + k^2)}, k^2geq 1$$
with $$k = 1$$. Solving the equation system for $$mathbf{A}$$ gives us the coefficients to use in equation $$textbf{(1)}$$ when interpolating between known values.

Interpolating displacements

How can RBF’s be used for deforming geometry? Well assume that the deformation is known for $$M$$ 3D positions $$mathbf{x}_i$$ and that this information is represented by a vector describing 3D displacement $$mathbf{u}_i$$ of the geometry that was positioned at $$mathbf{x}_i$$ in the original, undeformed state. You can think of the $$mathbf{x}_i$$ positions as control points that have been moved to positions $$mathbf{x}_i+mathbf{u}_i$$. The RBF interpolation method can now be used for interpolating these displacements to other positions.

Using the notation $$mathbf{x}_i = (x_i, y_i, z_i)$$ and $$mathbf{u}_i = (u^x_i, u^y_i, u^z_i)$$ three linear systems are set up as above letting the displacements $$u$$ be the quantity we called $$a$$ in the previous section:

$$!mathbf{G} mathbf{A}_x = (u^x_1, u^x_2, ldots, u^x_M, 0, 0, 0, 0)^T$$
$$!mathbf{G} mathbf{A}_y = (u^y_1, u^y_2, ldots, u^y_M, 0, 0, 0, 0)^T$$
$$!mathbf{G} mathbf{A}_z = (u^z_1, u^z_2, ldots, u^z_M, 0, 0, 0, 0)^T$$

where $$mathbf{G}$$ is assembled as described above. Solving for $$mathbf{A}_x$$, $$mathbf{A}_y$$, and $$mathbf{A}_z$$ involves a single matrix inversion and three matrix-vector multiplications and gives us the coefficients for interpolating displacements in all three directions by the expression $$mathbf{(1)}$$

The source code

In the source code accompanying this tutorial you will find the class RBFInterpolator which has an interface like this:

class RBFInterpolator
{
public:
	RBFInterpolator();
	~RBFInterpolator();

	//create an interpolation function f that obeys F_i = f(x_i, y_i, z_i)
	RBFInterpolator(vector x, vector y, vector z, vector F);

	//specify new function values F_i while keeping the same
	void UpdateFunctionValues(vector F);

	//evaluate the interpolation function f at the 3D position (x,y,z)
	real interpolate(float x, float y, float z);

private:
            ...
};

This class implements the interpolation method described above using the newmat matrix library. It is quite easy to use: just fill stl::vectors with the $$x_i$$, $$y_i$$ and $$z_i$$ components of the positions where the value $$F$$ is known and another stl::vector with the $$F_i$$ values. Then pass these vectors to the RBFInterpolator constructor, and it will be ready to interpolate. The $$F$$ value at any position is then evaluated by calling the ‘interpolate’ function. If some of the $$F_i$$ values change at any time, the interpolator can be quickly updated using the ‘UpdateFunctionValues’ method.

We want to deform a triangle surface mesh. These are stored in a class TriangleMesh, and loaded from OBJ files.
In the source code the allocation of stl::vectors describing the control points and the initialisation of RBFInterpolators looks like this:

void loadMeshAndSetupControlPoints()
{
	// open an OBJ file to deform
	string sourceOBJ = "test_dragon.obj";
	undeformedMesh = new TriangleMesh(sourceOBJ);
	deformedMesh = new TriangleMesh(sourceOBJ);

	// we want 11 control points which we place at different vertex positions
	const int numControlPoints = 11;

	const int verticesPerControlPoint = ((int)undeformedMesh->getParticles().size())/numControlPoints;

	for (int i = 0; igetParticles()[i*verticesPerControlPoint].getPos();
		controlPointPosX.push_back(pos[0]);
		controlPointPosY.push_back(pos[1]);
		controlPointPosZ.push_back(pos[2]);
	}

	// allocate vectors for storing displacements
	for (unsigned int i = 0; i
Now all displacements are set to zero vectors – not terribly exciting! To make it a bit more fun we can vary the displacements with time:
// move control points
	for (unsigned int i = 0; i < controlPointPosX.size(); i++ )
	{
		controlPointDisplacementX[i] = displacementMagnitude*cosf(time+i*timeOffset);
		controlPointDisplacementY[i] = displacementMagnitude*sinf(2.0f*(time+i*timeOffset));
		controlPointDisplacementZ[i] = displacementMagnitude*sinf(4.0f*(time+i*timeOffset));
	}

	// update the control points based on the new control point positions
	rbfX.UpdateFunctionValues(controlPointDisplacementX);
	rbfY.UpdateFunctionValues(controlPointDisplacementY);
	rbfZ.UpdateFunctionValues(controlPointDisplacementZ);

	// deform the object to render
	deformObject(deformedMesh, undeformedMesh);

Here the function ‘deformObject’ looks like this:
// Code for deforming the mesh 'initialObject' based on the current interpolation functions (global variables). 
// The deformed vertex positions will be stored in the mesh 'res'
// The triangle connectivity is assumed to be already correct in 'res'  
void deformObject(TriangleMesh* res, TriangleMesh* initialObject)
{
	for (unsigned int i = 0; i < res->getParticles().size(); i++)
	{
		Vector3 oldpos = initialObject->getParticles()[i].getPos();

		Vector3 newpos;
		newpos[0] = oldpos[0] + rbfX.interpolate(oldpos[0], oldpos[1], oldpos[2]);
		newpos[1] = oldpos[1] + rbfY.interpolate(oldpos[0], oldpos[1], oldpos[2]);
		newpos[2] = oldpos[2] + rbfZ.interpolate(oldpos[0], oldpos[1], oldpos[2]);

		res->getParticles()[i].setPos(newpos);
	}
}

That’s it!!! Now I encourage you to download the source code and play with it. Perhaps you can experiment with other radial basis functions? Or make the dragon crawl like a caterpillar? If you code something interesting based on this tutorial send a link to me and we will link to it from this page 
Karsten Noe
I got a mail from Woo Won Kim from Yonsei University in South Korea who has made a head modeling program that can generate 3D human heads from two pictures of the person using code from this RBF tutorial. Check out a video of this here.