2D transformation matrices baking

While playing with OpenGL, I was facing the situation where the common way to use matrices can quickly lead to CPU overload.

If you know the 3 common transformations in 2D (in homogeneous coordinates):

translation

scale

rotation

you also surely know that you can compose with these simple transformations to build more complex transformations. That is simply done with the use of matrix multiplication.

For example, if you have a sprite that you want to rotate and scale before translating it to have its own center exactly at the center of the screen. You would simply build 4 simple matrices T1, S, R, T2 and do the multiplication to obtain the final matrix M:

M = T2 * R * S * T1

with:
T1 the translation matrix moving your sprite to match his local center with the origin
S a scale matrix
R a rotation matrix
T2 a translation matrix moving your sprite at the center of the screen

That is a very common way to build complex transformations, but when you need to deal with thousands of them, you quickly overload your CPU.

Why ? Because generic matrix 3x3 multiplication involves 27 scalar multiplications and 18 scalar additions, as shown in the code below to compute C = AxB:

C[0][0] = A[0][0] * B[0][0] + A[1][0] * B[0][1] + A[2][0] * B[0][2];
C[1][0] = A[0][0] * B[1][0] + A[1][0] * B[1][1] + A[2][0] * B[1][2];
C[2][0] = A[0][0] * B[2][0] + A[1][0] * B[2][1] + A[2][0] * B[2][2];

C[0][1] = A[0][1] * B[0][0] + A[1][1] * B[0][1] + A[2][1] * B[0][2];
C[1][1] = A[0][1] * B[1][0] + A[1][1] * B[1][1] + A[2][1] * B[1][2];
C[2][1] = A[0][1] * B[2][0] + A[1][1] * B[2][1] + A[2][1] * B[2][2];

C[0][2] = A[0][2] * B[0][0] + A[1][2] * B[0][1] + A[2][2] * B[0][2];
C[1][2] = A[0][2] * B[1][0] + A[1][2] * B[1][1] + A[2][2] * B[1][2];
C[2][2] = A[0][2] * B[2][0] + A[1][2] * B[2][1] + A[2][2] * B[2][2];

Building M involved 3 matrix multiplications (T2 * R * S * T1), so the total for one sprite is 81 scalar multiplications and 53 scalar additions. It is quite a lot when you have that amount of operations each frame for thousand of sprites.


The solution ? Bake your matrix !

That could seem really complicated to set directly a single matrix that is the result of many matrices composition. But in fact that is not, because the 3 common transformation matrices (translate, rotation, scale) contain a lot of 0 and 1.

For example, here is the matrix M baked for T2 * R * S * T1:

with:
t1x, t1y : the 1st translation
sx, sy : the scale factors
a : the rotation angle
t2x, t2y : the 2nd translation

This matrix looks complicated, but it involves only 12 scalar multiplications and 4 scalar additions. That is such an improvement if we compare with the 81 multiplications and 53 additions of the previous method. Putting that in your code instead of explicitly computing the result of T2 * R * S * T1 saves a lot of CPU resources.

Here is my code using float as scalars:

void mat3TRST(float& t2x, float &t2y, float &rot, float &sx, float &sy, float &t1x, float &t1y, Matrix3 &matRes)
{
float cRot, sRot;
cRot = cos(rot);
sRot = sin(rot);
matRes[0][0] = sx*cRot;
matRes[1][0] = sy*-sRot;
matRes[2][0] = t1x*sx*cRot + t1y*sy*-sRot + t2x;
matRes[0][1] = sx*sRot;
matRes[1][1] = sy*cRot;
matRes[2][1] = t1x*sx*sRot + t1y*sy*cRot + t2y;
matRes[0][2] = 0.0f;
matRes[1][2] = 0.0f;
matRes[2][2] = 1.0f;
}


Baking is easy
How did I obtain that result ? I didn't bake the M matrix "by hand". Instead I used WolframAlpha for that.

My input was:

{{1,0,i},{0,1,j},{0,0,1}} . {{e,f,0},{g,h,0},{0,0,1}} . {{c,0,0},{0,d,0},{0,0,1}} . {{1,0,a},{0,1,b},{0,0,1}}

It is my T2 * R * S * T1 matrix composition, but using Wolfram syntax. From that I obtained the resulting matrix, still involving the constants a,b,c,d,e,f,g,h,i,j:

{{c e, d f, a c e + b d f + i}, {c g, d h, a c g + b d h + j}, {0, 0, 1}}

My job was only to copy/paste it in my code and replace the letters with the good variable names !

Bitmap triangulation extraction demo

Simply draw on the surface below. Then press SPACE to extract the triangulation from your picture and embed it into Box2D and see what happen. Press ENTER to clean and play again.

The algo is fast and accurate for almost all pictures. It naturally manages holes and interlocked shapes.

How does it work ? Just few lines of code with the help of Daedalus Lib. Have a look to the wiki page if you need to be conviced !

A playable dungeon map generator

Last days I recycled my dungeon map generator and I plugged it to the engine of my proto 61 Cygni.

You can experiment the result below. Generate a map (with custom parameters or randomly) and then click  the "play it" button.

Just for fun, you can directly catch 4 weapons: the shotgun, the machine gun, the plasmagun and the rocketgun ! You will use these to kill the bots ditributed on the level. Additionally, you can run by pushing quickly the UP key 2 times and dodge with LEFT or RIGHT keys quickly 2 times.

Introducing Daedalus Lib

Project on GitHub

Daedalus is a library I develop. It manages 2D environment modeling and pathfinding. I really hope this library will help designers and developers to invent new gameplays.

When I began to code Daedalus, I had many ideas in mind:

1- focus on 2D
2- fastness and accuracy
3- simplicity of use


1- Why focus only on 2D and not on 3D ? Because constraining to 2D allows to gain simplicity and efficiency. Many great games today are still based on 2D engines involving 2D mechanics, 2D physics and 2D display ; so I hope Daedalus will find his place as a new component for new 2D projects.


2- Fastness and accuracy are reached by using among the best techniques available in the fields of computational geometry : quad-edge structure and fully dynamic Delaunay triangulation. Daedalus algorithms are based on many research publications, among them:

Fully Dynamic Constrained Delaunay Triangulation by Kallmann, Bieri and Thalmann
An improved incremental algorithm for constructing... by Anglada
Efficient Triangulation-Based Pathfinding by Jon Demyen


3- For simplicity, I assumed that the library should work without any pre-generated data. Everything should work in real time : obstacles insertion/motion/deletion and path generation. Also I assumed that the library should be fault tolerant to designer/developer mistakes : obstacles can overlap and can be of any shape : open, convex, concave... At last, the path generation manages non-null size objects in order to avoid any obstacle collision.


Now the basics are exposed, I can show you some demos :

Today Daedalus is coded in Actionscript. But because the library has no dependancies from any other library (I wrote the whole mathematics), it can easily and quickly be translated in any other language. The project is not open-source at this moment, so not available for free download. I plan to use it in some new projects before release it in the open world.

So I am now looking for some collaboration. If you feel interested by Daedalus for your game engine, for a commercial game project or any other application, you can contact me directly at flash dot cedric at google mail service.

A dungeon map generator

I try to find new methods and algorithms to proceduraly build dungeons and mazes. Something that definitevely breaks up with the boring square and rectangular aligned patterns.

Finaly I had interesting results using Daedalus lib, as illustrated in the examples below:

The method I used to generate this map has many interesting properties. These are easier to understand by looking the steps of generation.

Step 1:

We generate a simple Delaunay triangulation. We use an algorithm that iterates through a nxm grid of points and add some of them in the triangulation according to a fixed probability P. With P=1, the triangulation would be a full regulat grid. With P=0, the triangulation would be empty. See a result below with a 20x20 grid and P=0.5 :

Step 2:

Considering the triangluation as a graph, starting from the vertex in the center, we use a custom depth-first algorithm to extract a sub-graph. At depth n, the next node at n+1 is choosen at random among unvisited nodes. Paramaters used are the total nodes count NC, the maximum branch depth BD and the maximum branches count BC. See a result below with NC=24 :

See that sub-graph as the underlying navigation graph of our futur dungeon. Nodes represent rooms and edges represent accessibility bewteen them.

Step 3:

We build the dual of the triangulation. It is similar to the Voronoi diagram but for better results we use the average positions instead of circumcircles centers:

The dual represents the real shape that will be used for the dungeon.

Step 4:

We keep only the dual cells surrounding the previously generated sub-graph at step 2. Then we build a fresh new triangulation from them. Additionally, we dig a door in the edges shared by connected rooms :

Step 5:

Finally we use a chamfer algorithm to add thickness to the walls :

Notice 2 importants properties of the resulting map:

1. It is directly built on a Delaunay triangulation

2: The navigation graph is given

In conclusion, navigation and pathfinding algorithms are ready to use, without any extra cost. It means that AI can navigate efficiently and accurately through the generated map using the pathfinding solution included in Daedalus.

Pathfinding on bitmap triangulation

Because large triangulated maps can be a pain to design by hand, I implemented into Daedalus an algorithm generating optimized triangulations from bitmap images. It is directly inspired by the bitmap segmentation in the Potrace algorithm.

How does-it work ?

Just give a picture made from black and white pixels to the algorithm. For example:

Then the algorithm returns to you a clean and beautiful triangulation:

Neither parameter nor configuration filling are required. Just give the bitmap and get the mesh.

However, black and white pictures are required. But we can consider the use of a threshold to easily convert gray scale and color images.

From that, I experimented with success the use of pathfinding on several triangulations generated from:

- a SNES Mario Kart circuit map

- a Doom map

- a picture of the labyrinth at Grace Cathedral

- a map of Paris created by artist Jazzberry Blue

I used Photoshop to extract properly black and white pictures from the previous maps. Then I just played with the algorithm and enjoyed the fast and accurate pathfinding implemented in Daedalus Lib.

As you can see by yourself:

Cartoons into Daedalus Lib

I recently added a bitmap segmentation algorithm into Daedalus Lib. Inspired by the Potrace algorithm, the main purpose was to use automatic generation of large meshes to check the reliability of  the algorithms inside Daedalus.

Some interesting results:

Mickey

Astérix

Doraemon

Bob

Accurate point in triangle test

Many resources on the web deal with the 2D point in triangle test. This post will summarize the most famous methods to solve it and will show why in some cases they are not accurate and can lead to errors. As a conclusion, we will expose a new algorithm.

At first, we need to remember the problem. We are in 2D. We have a triangle T defined by 3 points p1(x1, y1), p2(x2, y2), p3(x3, y3) and a single point p(x, y). Does this single point p lies inside the triangle T ?

1st method : barycentric coordinate system

Barycentric coordinate allows to express new p coordinates as a linear combination of p1, p2, p3. More precisely, it defines 3 scalars a, b, c such that :

x = a * x1 + b * x2  + c * x3

y = a * y1 + b * y2 + c * y3

a + b + c = 1

The way to compute a, b, c is not difficult :

a = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3)) / ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3))

b = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3)) / ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3))

c = 1 - a - b

Then we just need to apply the interesting following property :

p lies in T if and only if 0 <= a <= 1 and 0 <= b <= 1 and 0 <= c <= 1

Code sample :

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var denominator:Number = ((y2 - y3)*(x1 - x3) + (x3 - x2)*(y1 - y3));
var a:Number = ((y2 - y3)*(x - x3) + (x3 - x2)*(y - y3)) / denominator;
var b:Number = ((y3 - y1)*(x - x3) + (x1 - x3)*(y - y3)) / denominator;
var c:Number = 1 - a - b;

return 0 <= a && a <= 1 && 0 <= b && b <= 1 && 0 <= c && c <= 1;
}

2nd method : parametric equations system

Here the idea is to consider the parametric expressions of the 2 edges [p1, p2] and [p1, p3] in T :

x(t1) = t1*(x2 - x1)

y(t1) = t1*(y2 - y1)

x(t2) = t2*(x3 - x1)

y(t2) = t2*(y3 - y1)

Then express p(x, y) as a linear combination of them :

x = x1 + x(t1) + x(t2)

y = y1 + y(t1) + y(t2)

Solving the system, it gives to us :

t1 = (x*(y3 - y1) + y*(x1 - x3) - x1*y3 + y1*x3) / (x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3)

t2 = (x*(y2 - y1) + y*(x1 - x2) - x1*y2 + y1*x2) / -(x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3)

Finally, we just need to apply the interesting following property :

p lies in T if and only if 0 <= t1 <= 1 and 0 <= t2 <= 1 and t1 + t2 <= 1

Code sample :

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var denominator:Number = (x1*(y2 - y3) + y1*(x3 - x2) + x2*y3 - y2*x3);
var t1:Number = (x*(y3 - y1) + y*(x1 - x3) - x1*y3 + y1*x3) / denominator;
var t2:Number = (x*(y2 - y1) + y*(x1 - x2) - x1*y2 + y1*x2) / -denominator;
var s:Number = t1 + t2;

return 0 <= t1 && t1 <= 1 && 0 <= t2 && t2 <= 1 && s <= 1;
}

3rd method : check sides with dot product

Maybe the most famous method, based on dot product. We assume that p1, p2, p3 are ordered in counterclockwise. Then we can check if p lies at left of the 3 oriented edges [p1, p2], [p2, p3] and [p3, p1].

For that, first we need to consider  the 3 vectors v1, v2 and v3 that are respectively left-orthogonal to [p1, p2], [p2, p3] and [p3, p1] :

v1 = <y2 - y1, -x2 + x1>

v2 = <y3 - y2, -x3 + x2>

v3 = <y1 - y3, -x1 + x3>

Then we get the 3 following vectors :

v1' = <x - x1, y - y1>

v2' = <x - x2, y - y2>

v3' = <x - x3, y - y3>

At last, we compute the 3 dot products :

dot1 = v1 . v1' = (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1)

dot2 = v1 . v2' = (y3 - y2)*(x - x2) + (-x3 + x2)*(y - y2)

dot3 = v3 . v3' = (y1 - y3)*(x - x3) + (-x1 + x3)*(y - y3)

Finally, we can apply the interesting property :

p lies in T if and only if 0 <= dot1 and 0 <= dot2 and 0 <= dot3

Code sample :

function side(x1, y1, x2, y2, x, y:Number):Number
{
return (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1);
}

function pointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var checkSide1:Boolean = side(x1, y1, x2, y2, x, y) >= 0;
var checkSide2:Boolean = side(x2, y2, x3, y3, x, y) >= 0;
var checkSide3:Boolean = side(x3, y3, x1, y1, x, y) >= 0;
return checkSide1 && checkSide2 && checkSide3;
}

These 3 methods are quite good to solve the point in triangle test. Purely mathematically speaking, they must validate any point inside the triangle and even those lying exactly on the boundary (on any edge).

Accuracy problems

Despite the strong mathematical background of our methods, in some cases they can lead to a lack of accuracy because the floating-point number system have limited size and most of the time it deals with approximations. The problem occurs sometimes when the point p should exactly on one triangle's edge ; the approximations lead to fail the test.

Example 1*:

We consider the triangle T defined by 3 points p1(x1, y1), p2(x2, y2), p3(x3, y3) with values :

x1 = 1/10

y1 = 1/9

x2 = 100/8

y2 = 100/3

x3 = 100/4

y3 = 100/9

and a single point p(x, y) lying exactly on the segment [p1, p2] :

x = x1 + (3/7)*(x2 - x1)

y = y1 + (3/7)*(y2 - y1)

If we apply the barycentric method, we get the 3 following values for a, b, c :

a : 0.5714285714285715

b : 0.42857142857142855

c : -5.551115123125783e-17

Because c < 0, the test fails to validate the point inside the triangle. In many applications, this situation is not really a problem because a lack of accuracy is not a tragedy.

But in some situations, it can be really annoying.

Example 2*:

We consider 2 triangles T and T' defined respectively by the points p1(x1, y1), p2(x2, y2), p3(x3, y3) and p1'(x1', y1'), p2'(x2', y2'), p3'(x3', y3'). Values are :

x1 = 1/10

y1 = 1/9

x2 = 100/8

y2 = 100/3

x3 = 100/4

y3 = 100/9

x1' = x1

y1' = y1

x2' = x2

y2' = y2

x3' = -100/8

y3' = 100/6

and a single point p(x, y) lying exactly on the segment [p1, p2] :

x = x1 + (3/7)*(x2 - x1)

y = y1 + (3/7)*(y2 - y1)

The situation is quite simple : 2 non-overlapping triangles sharing one edge and one single point lying on this edge. We know from the previous example that the barycentric method fails to validate p inside the triangle T :

a : 0.5714285714285715

b : 0.42857142857142855

c : -5.551115123125783e-17

But more surprisely, the method fails again when applied to the triangle T' :

a : 0.5714285714285715

b : 0.4285714285714285

c : -1.1102230246251565e-16

Mathematically speaking, the point p belongs to both triangles. In practice, we could tolerate that a lack of accuracy leads to validate only one triangle belonging. But here we face a complete invalidation and the point is detected as outside both triangles.

We could suppose that the barycentric-based method is intrinsically inaccurate and one other method will lead to satisfiable results. But in reality the problem remains.

Example 3*:

We consider the oriented edge E defined by 2 points p1(x1, y1) and p2(x2, y2) with values :

x1 = 1/10

y1 = 1/9

x2 = 100/8

y2 = 100/3

and a single point p(x, y) lying exactly on edge E :

x = x1 + (3/7)*(x2 - x1)

y = y1 + (3/7)*(y2 - y1)

If we apply the check side method with E and p, we get the following dot product :

dot = (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1)

dot = -2.842170943040401e-14

Then if we apply the same method with E' the reversed edge of E, we get the same result :

dot = (y1 - y2)*(x - x2) + (-x1 + x2)*(y - y2)

dot = -2.842170943040401e-14

It means that the check side test gives us 2 contradictory results. The point looks on the right sides of both the edge E and his reversed ; that is of course mathematically impossible. We face the same situation as using barycentric method : if the edge E is shared by 2 non-overlapping triangles T and T', despite the fact that p lies mathematically on E, our method will lead us to tragically conclude that p is outside both triangles T and T'.

Any accurate solution ?

The answer is yes. We can write a complete algorithm leading to a safe point in triangle test by combining many familiar algorithms. The core of the method is to assume a real thickness value for the triangle's edges and vertices ; it contrasts with the original purely mathematical situation where the triangle's edges and vertices have an infinitesimal thickness. We will call it epsilon, because in practice we will keep it very small (near 0.001).

The steps of the algorithm are:

1: use a bounding box as a fast pre-validation test
Because our new algorithm will involve more computations, it is convenient to add a pre-test to reject very quickly most of the false cases.

The bounding box is simply the min/max of the x/y values among the 3 triangle's vertices, slightly inflated by the epsilon value.

Code sample :

const EPSILON:Number = 0.001;

function pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var xMin:Number = Math.min(x1, Math.min(x2, x3)) - EPSILON;
var xMax:Number = Math.max(x1, Math.max(x2, x3)) + EPSILON;
var yMin:Number = Math.min(y1, Math.min(y2, y3)) - EPSILON;
var yMax:Number = Math.max(y1, Math.max(y2, y3)) + EPSILON;

if ( x < xMin || xMax < x || y < yMin || yMax < y )
return false;
else
return true;
}

2: use any method studied prevously (barycentric, parametric or dot product)
If the test is positive, we can trust it and stop the algorithm immediatly. But if the test is negative, maybe we face the situation with the point lying on one triangle's edge : then we need more investigations, involving methods using the epsilon value.

3: use the point to segment distance
For every edge of the triangle, we compute the shortest distance between the edge and the point to evaluate. If the distance is shorter that epsilon, we can validate definitely the test.

In detail, given 3 points  p, p1 and p2, a very tricky use of the dot product allows us to check efficiently the relative position of the orthogonal projection p' of p on the infinite line passing through p1 and p2. If the projection lies between p1 and p2 then we compute the distance p and p'. Otherwise, we compute the distance between p and the nearest among p1 and p2.

Finally, because our algorithm use distances for comparisons only, we will restrict our computations to square distances only (faster because we can omit the square root).

Code sample:

function distanceSquarePointToSegment(x1, y1, x2, y2, x, y:Number):Number
{
var p1_p2_squareLength:Number = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
var dotProduct:Number = ((x - x1)*(x2 - x1) + (y - y1)*(y2 - y1)) / p1_p2_squareLength;
if ( dotProduct < 0 )
{
return (x - x1)*(x - x1) + (y - y1)*(y - y1);
}
else if ( dotProduct <= 1 )
{
var p_p1_squareLength:Number = (x1 - x)*(x1 - x) + (y1 - y)*(y1 - y);
return p_p1_squareLength - dotProduct * dotProduct * p1_p2_squareLength;
}
else
{
return (x - x2)*(x - x2) + (y - y2)*(y - y2);
}
}

Final sample code

Here is the code illustrating the steps we described in the previous section.

const EPSILON:Number = 0.001;
const EPSILON_SQUARE:Number = EPSILON*EPSILON;

function side(x1, y1, x2, y2, x, y:Number):Number
{
return (y2 - y1)*(x - x1) + (-x2 + x1)*(y - y1);
}

function naivePointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var checkSide1:Boolean = side(x1, y1, x2, y2, x, y) >= 0;
var checkSide2:Boolean = side(x2, y2, x3, y3, x, y) >= 0;
var checkSide3:Boolean = side(x3, y3, x1, y1, x, y) >= 0;
return checkSide1 && checkSide2 && checkSide3;
}

function pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
var xMin:Number = Math.min(x1, Math.min(x2, x3)) - EPSILON;
var xMax:Number = Math.max(x1, Math.max(x2, x3)) + EPSILON;
var yMin:Number = Math.min(y1, Math.min(y2, y3)) - EPSILON;
var yMax:Number = Math.max(y1, Math.max(y2, y3)) + EPSILON;

if ( x < xMin || xMax < x || y < yMin || yMax < y )
return false;
else
return true;
}

function distanceSquarePointToSegment(x1, y1, x2, y2, x, y:Number):Number
{
var p1_p2_squareLength:Number = (x2 - x1)*(x2 - x1) + (y2 - y1)*(y2 - y1);
var dotProduct:Number = ((x - x1)*(x2 - x1) + (y - y1)*(y2 - y1)) / p1_p2_squareLength;
if ( dotProduct < 0 )
{
return (x - x1)*(x - x1) + (y - y1)*(y - y1);
}
else if ( dotProduct <= 1 )
{
var p_p1_squareLength:Number = (x1 - x)*(x1 - x) + (y1 - y)*(y1 - y);
return p_p1_squareLength - dotProduct * dotProduct * p1_p2_squareLength;
}
else
{
return (x - x2)*(x - x2) + (y - y2)*(y - y2);
}
}

function accuratePointInTriangle(x1, y1, x2, y2, x3, y3, x, y:Number):Boolean
{
if (! pointInTriangleBoundingBox(x1, y1, x2, y2, x3, y3, x, y))
return false;

if (naivePointInTriangle(x1, y1, x2, y2, x3, y3, x, y))
return true;

if (distanceSquarePointToSegment(x1, y1, x2, y2, x, y) <= EPSILON_SQUARE)
return true;
if (distanceSquarePointToSegment(x2, y2, x3, y3, x, y) <= EPSILON_SQUARE)
return true;
if (distanceSquarePointToSegment(x3, y3, x1, y1, x, y) <= EPSILON_SQUARE)
return true;

return false;
}


I hope you appreciated.

* all the tests are done in Actionscript 3, typing values with the native Number type (IEEE-754 double-precision floating-point number).

Daedalus : group pathfinding and FOV

Today 2 new demos made with Daedalus, showing in action 2 new components: group pathfinding and field of view.

The group pathfinding has many interesting features:
- it can manage any number of entities
- it can manage any size of entities
- generated paths are automatically sampled in order to synchronise the entities and avoid collision between them

We could just reproach the single common path the entities use. Indeed, if the entities have different sizes, it could be legitimate to expect them to reach the goals by different optimized paths and keep synchronisation to avoid collision. But in fact, this behavior is the subject of a new component coming soon, more focused on synchronisation for entities having completely different start and goal positions.


Then the field of view, while looking not so impressive, is a very important component to give awareness to AI entities, allowing them to detect cleverly the other entities. The features are:
- radius, opening angle, position and direction as dynamic properties
- the FOV is broken by constrained edges (other entities can hide behind walls)
- it manages non-null size entites

In the details, the FOV implementation is simply an accurate clipping system. So it can really manages the visibility of an entity through very complex situation with partially overlapping walls.

Demo : fully dynamic Delaunay triangulation

Demo : pathfinding on Delaunay constrained triangulation

A simple planar mesh primitive with half-edge

In a previous article, I explained why and how to implement the half-edge data structure. This article will focus on the implemention of a simple primitive: a 2 polygons rectangular planar mesh. Later we will see how to add some useful tools allowing us to expand the primitive by adding/removing vertices and flipping edges,

The implementation of a simple 2 polygons rectangle primitive could look easy at first glimpse. In reality, with half-edge, it is not so obvious. We need to be very careful and set all the relevant adjacency relations between the 4 vertices, the 10 oriented edges and the 2 faces:

Remember that each edge must references 4 datas: the origin vertex, the opposite edge, the next left edge and the left face:

We immediatly face an annoying problem : the 4 border edges e01, e12, e23 and e30 have neither left face nor next left edge. We could think that leaving null references as datas for these edges is safe, but in fact it would lead to several problems because the missing core datas would break the way we hope to iterate through the mesh. That's why I suggest rather to extend the rectangle in order to obtain a complete closed and consistent mesh:

We just added 2 new edges e13, e31 and 2 new faces f2 (bounded by e12, e23, e31) and f3 (bounded by e30, e01, e13)  This figure could look strange at first glimpse, but in fact it is really intuitive if you imagine our planar mesh lying on a sphere:

Finally, if these 4 new elements are convenient to get a mesh with a complete adjacency structure, they can be cumbersome when you display your mesh on screen or when you iterate in order to navigate from element to element. So feel free to add a new visible property to your classes and set it to false for any object lying outside the rectangle v0, v1, v2, v3. In this way you will be able to simply skip them.

As a conclusion, we will simply give now the complete declarations necessary to implement this simple rectangular primitive:

// (x, y) is any 2d coordinates system
// w is the rectangle width
// h is the rectangle heighy

v0.position = (0, 0)
v0.edge = e01 // or e02 or e03

v0.visible = true

v1.position = (w, 0)
v1.edge = e12 // or e10 or e13

v1.visible = true

v2.position = (w, h)
v2.edge = e23 // or e21 or e20

v2.visible = true

v3.position = (0, h)

v3.edge = e30 // or e32 or e31

v3.visible = true

e01.originVertex = v0
e01.oppositeEdge = e10
e01.nextLeftEdge = e13
e01.leftFace = f3
e01.visible = true

e10.originVertex = v1
e10.oppositeEdge = e01
e10.nextLeftEdge = e02
e10.leftFace = f0
e10.visible = true

e12.originVertex = v1
e12.oppositeEdge = e21
e12.nextLeftEdge = e23
e12.leftFace = f2
e12.visible = true

e21.originVertex = v2
e21.oppositeEdge = e12
e21.nextLeftEdge = e10
e21.leftFace = f0
e21.visible = true

e23.originVertex = v2
e23.oppositeEdge = e32
e23.nextLeftEdge = e31
e23.leftFace = f2
e23.visible = true

e32.originVertex = v3
e32.oppositeEdge = e23
e32.nextLeftEdge = e20
e32.leftFace = f1
e32.visible = true

e30.originVertex = v3
e30.oppositeEdge = e03
e30.nextLeftEdge = e01
e30.leftFace = f3
e30.visible = true

e03.originVertex = v0
e03.oppositeEdge = e30
e03.nextLeftEdge = e32
e03.leftFace = f1
e03.visible = true

e02.originVertex = v0
e02.oppositeEdge = e20
e02.nextLeftEdge = e21
e02.leftFace = f0
e02.visible = true

e20.originVertex = v2
e20.oppositeEdge = e02
e20.nextLeftEdge = e03
e20.leftFace = f1
e20.visible = true

e13.originVertex = v1
e13.oppositeEdge = e31
e13.nextLeftEdge = e30
e13.leftFace = f3
e13.visible = false

e31.originVertex = v3
e31.oppositeEdge = e13
e31.nextLeftEdge = e12
e31.leftFace = f2
e31.visible = false

f0.edge = e10 // or e02 or e21
f0.visible = true

f1.edge = e32 // or e20 or e03
f1.visible = true

f2.edge = e31 // or e12 or e23
f2.visible = false

f3.edge = e13 // or e01 or e30
f3.visible = false

Doubly connected edge list implementation

If you are interested by the half-edge or curious about sophisticated data structures, this post is for you. I know it is very difficult for people new in the field of computational geometry to find relevant informations about implementation of half-edge structure, so I will try to explain the concepts very carefully.

First, what is the motivation behind the use of half-edge ? Maybe you already know data structures like linked-list, tree, or graph, but what is half-edge and why you would need to implement it ?

One answer I would give is the following: you need a half-edge structure when you need to navigate very quickly and easily inside a polygon mesh. Remember that you need polygon mesh not only for 3D, but also for some 2D complex datas representation like Delaunay triangulation (note that in this article we will deal only with 3 sides polygon mesh).

Now the problem is that in many libraries or softwares, the polygon mesh datas are stored inside 2 arrays: one that explicitly stores the vertices positions and one that stores vertices indexes for each triangle. For example, it is the case in Unity3D Mesh class. We can illustrate this by the following diagram:

This data structure is quite simple to understand and to implement, but it is a performance disaster when you need to write an algorithm that needs to navigate and search efficiently inside the mesh. One obvious example of navigation is pathfinding. But why is it not efficient ? Because the adjacency access is not. For example, if you want to find all the vertices connected by an edge to a given vertex, you need to scan the entire triangle array to find the relevant indexes ; so the time you spend to do it increases by the size of your mesh.

But with half-edge you can access adjacency in constant time, because in this structure it becomes the explicit way to keep the datas. For this purpose, the main object you will implement will be not a vertex or a face, but an oriented edge. Think oriented edge as the core node of your data, as we can see in this illustration:

Given an Edge object current edge, a descent half-edge API should give you access to the complete adjacency: Vertex, Edge and Face objects adjacent to current edge. For example, my own interface for half-edge structure as a library contains 4 core classes: Vertex, Edge, Face and Mesh with the following public methods :

public class Vertex

• public get position : Vector // (can be any 2D or 3D coordinates data)
• public get edge : Edge // (any Edge with this Vertex as origin)

public class Edge

• public get originVertex : Vertex
• public get destinationVertex : Vertex
• public get oppositeEdge : Edge
• public get nextLeftEdge : Edge
• public get prevLeftEdge : Edge
• public get nextRightEdge : Edge
• public get prevRightEdge : Edge
• public get rotLeftEdge : Edge
• public get rotRightEdge : Edge
• public get leftFace : Face
• public get rightFace : Face

public class Face

• public get edge : Edge // (any Edge with this face adjacent at left)

public class Mesh

• public get vertices : Vertex []
• public get edges : Edge []
• public get faces : Face []

Equiped with that, you must be able to iterate efficiently in order to navigate through the whole mesh. For example, from any Vertex v, you can iterate to any direct adjacent Vertex by calling a sequence like:

v.edge.destinationVertex // 1st adjacent vertex
v.edge.rotLeftEdge.destinationVertex // 2nd adjacent vertex at left
v.edge.rotLeftEdge.rotLeftEdge.destinationVertex // 3rd adjacent vertex at left

v.edge.rotLeftEdge. ... .rotLeftEdge.destinationVertex // n-th adjacent vertex at left


The same idea works for faces. From any Face f, you can reach the 3 direct adjacent Face just by:

f.edge.rightFace // 1st adjacent face at left
f.edge.nextLeftEdge.rightFace // 2nd adjacent face at left
f.edge.nextLeftEdge.nextLeftEdge.rightFace // 3rd adjacent face at left

When you understand the basic principles of navigation, then maybe you want to expose a compliant API in order to iterate more easily through your mesh. So it can be done very easily by implementing a collection of iterators on top of the half-edge structure. For example, some ideas of iterators you can add to your library:

public class IteratorFromVertexToNeighbourVertices

• public function set fromVertex(value:Vertex)
• public function get next : Vertex

public class IteratorFromFaceToNeighbourFaces

• public function set fromFace(value:Face)
• public function get next : Face

public class IteratorFromEdgeToRotatedEdges

• public function set fromEdge(value:Edge)
• public function get next : Edge

public class IteratorFromVertexToHoldingFaces

• public function set fromVertex(value:Vertex)
• public function get next : Face

Now come back to low level half-edge implementation. We saw that the list of 11 getter methods exposed by the Edge class is quite useful to use. But at first look, it will imply a lot of datas to maintain and in consequence a nightmare to implement. But don't worry, because in reality almost of the adjacency relations between the objects can be deduced from the very small core of datas showed in this illustration:

You can see here that given an Edge object current edge, only 4 datas are required: origin vertex opposite edge, next left edge and left face. As a justification, I will simply give you a possible implementation of the Edge class with the complete deduced adjacency relations:

public class Edge
// core datas:

• private var _originVertex : Vertex
• private var _oppositeEdge : Edge
• private var _nextLeftEdge : Edge
• private var _leftFace : Face

// public interface:

• public get originVertex : Vertex { return _originVertex }
• public get destinationVertex : Vertex { return _oppositeEdge.originVertex }
• public get oppositeEdge : Edge { return _oppositeEdge }
• public get nextLeftEdge : Edge { return _nextLeftEdge }
• public get prevLeftEdge : Edge { return _nextLeftEdge.nextLeftEdge }
• public get nextRightEdge : Edge { return _oppositeEdge.nextLeftEdge.nextLeftEdge.oppositeEdge }
• public get prevRightEdge : Edge { return _oppositeEdge.nextLeftEdge.oppositeEdge }
• public get rotLeftEdge : Edge { return _nextLeftEdge.nextLeftEdge.oppositeEdge }
• public get rotRightEdge : Edge { return _oppositeEdge.nextLeftEdge }
• public get leftFace : Face { return _leftFace  }
• public get rightFace : Face { return _oppositeEdge.leftFace }

This proves that only 4 datas need to be maintained in order to implement a complete half-edge structure.

So now that we have a descent core API, we need some additionnal tools to comfortably build concrete structures. For example, we could code a parser that generates a complete half-edge structured Mesh from a .3ds file (3D Studio Max native export format). An other possibility is to code some simple primitives and tools allowing us to expand them. We choose this second approach and it is the subject of the next articles:

A simple planar mesh primitive with half-edge

Incremental constrained Delaunay triangulation prototype

Delaunay triangulation is a huge subject and it can become very hard when we reach some refinements like constrained triangulation, incremental algorithm, quad-edge data structure... that are globally part of the computational geometry theory.

But there are several applications and they can be very impressive. Of course my mind is focused on artificial intelligence and I think especially of the pathfinding through navigation mesh.

This article is the first of a serie about concepts and algorithms related to computational geometry. It is an illustration of Delaunay triangulation through a demo prototype.

Just click on the board to add new constrained edges:



This demo and the underlying library I wrote are direct implementations of the following article:

Fully Dynamic Constrained Delaunay Triangulation (2003)

The underlying data structure behind the library is the quad-edge:

Quad-edge core implementation explained

A simple planar mesh primitive with quad-edge