Written by Oliver
Kreylos, hosted by 2DNow
- The 2D game developer's community!
Path Finding In Complex Maps And The Black Art
Of Linear Algebra
In this article I will describe my musings on
the problem of finding all shortest paths in a map, and how
to apply methods of Linear Algebra to gain a significant speed-up.
(A sample
ANSI C source file implementing the algorithm described
here and solving the problem used as an example on this page
is available for download.)
Let me first state the problem: Given is a map
that describes a virtual world. It doesn't matter if the world
is two- or three-dimensional, but for the sake of simplicity
I will assume it is 2D throughout this article. Thus, the map
is defined by line segments called walls, and the accessible
part of the map consists of one or more general polygons (see
Figure 1, left). If the map is "connected," i.e.,
if every point can be reached from any other point, there is
only one polygon; otherwise there is one for each connected
part of the map. The task is now to answer a large number of
"path queries," that is to report the shortest way
between a pair of points anywhere on the map. This problem often
arises in implementing artificial intelligence for games, where
a character has to go from its current position to one determined
by the AI, without bumping into any walls. To do this the AI
routines have to issue a path query for every character after
any AI decision, this leading to a large number of queries for
the same map.
To minimize the computing time for a single path
query, I trade off pre-processing time with game time and present
a pre-processing algorithm that leads to an O(k) path query
response time, where k is the actual length of the found path.
The first step in this algorithm is to decompose
the accessible polygon(s) into a set of convex polygons (see
Figure 1, right). Convexity ensures that every point inside
a polygon can be reached from any other point inside the same
polygon by walking a straight line, without crossing the polygon's
border. Thus, if a path query's source and destination points
lie both inside the same polygon, reporting a path becomes trivial.

Figure 1: A map and its accesible polygon. Left: General polygon;
right: set of convex polygons.
The next step is to number all polygons from
1 to N, and to find all pairs of directly adjacent polygons,
i.e. those that share a common edge.
What we have after these two steps is a representation
of the map as a set of convex polygons, and a set of polygon
pairs that share a common edge. We now use the following fact
to find all paths between polygons that are not directly connected:
If there exists a path from polygon i to polygon j, and there
also exists a path from polygon j to polygon k, we know there
is a path from polygon i to polygon k, by going from i to j
first and then from j to k.
From an algebraic point of view, the "connectedness"
of polygons defines a relation on the set of polygons. This
relation has the following properties:
- It is reflexive,i.e., every polygon is connected to itself
(trivial).
- It is symmetric,i.e., if there is a way from i to j, there
is also a way from j to i, by going backwards. This is,
of course, just an assumption. In maps where there are one-way
passages like cliffs (you can fall down, but not climb up)
or teleporters as in Doom, this doesn't hold. Luckily, the
following algorithms works for the non-symmetric case without
any changes.
- It is transitive, i.e., if i and j as well as j and k
are connected, so are i and k (see above paragraph for a
proof).
These properties make this relation an "equivalence relation",
but that shall be of no concern to us here.
This algebraic view comes in quite handy, since
the solution to our problem is what is called "transitive
closure" in Algebra, and there is a simple algorithm to
compute it. To do so, we have to write down the "connectedness"
relation as a matrix CM. For N polygons, CM is an N-by-N matrix,
whose entries CM[i][j] contain the length of the shortest path
from polygon i to polygon j, and to which polygon to go next
on that path. Let's define a data type (in C++, because operator
overloading will be of good use later on!):
struct MatrixEntry // Describes one step in the path from polygon i
// to polygon j
{
int pathLength; // Length of path in steps across polygon borders
int nextPolygon; // Number of polygon to go to next
};
We initialize the matrix CM as follows:
- All diagonal entries CM[i][i], which describe the path
from polygon i to polygon i, are set to {0,i}, because we
have to take zero steps and go directly to the destination.
- All entries CM[i][j], where i and j are adjacent, are
set to {1,j}, because we have to take one step and go to
polygon j next. If a shared edge of polygons i and j can
be crossed in both directions, we have to set both entries
CM[i][j] and CM[j][i].
- All other entries are set to {INFINITY,-1}, because the
currently shortest path is infinitely long, i.e., it doesn't
exist, and we don't know where to go next (we define INFINITY
as the largest possible integer, e.g., 65535).
The initial connection matrix CM for the example
map of Figure 1 looks like this (the first number is the path
length, the second is the next polygon; entries containing unknown
paths are left blank):
| |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
| 1 |
0, 1 |
1, 2 |
|
|
|
|
|
|
|
|
|
| 2 |
1, 1 |
0, 2 |
1, 3 |
|
|
|
|
|
|
|
|
| 3 |
|
1, 2 |
0, 3 |
|
|
|
1, 7 |
|
|
1, 10 |
|
| 4 |
|
|
|
0, 4 |
1, 5 |
|
|
|
|
|
1, 11 |
| 5 |
|
|
|
1, 4 |
0, 5 |
|
|
|
|
|
|
| 6 |
|
|
|
|
|
0, 6 |
1, 7 |
|
|
|
1, 11 |
| 7 |
|
|
1, 3 |
|
|
1, 6 |
0, 7 |
1, 8 |
|
|
|
| 8 |
|
|
|
|
|
|
1, 7 |
0, 8 |
1, 9 |
|
|
| 9 |
|
|
|
|
|
|
|
1, 8 |
0, 9 |
|
|
| 10 |
|
|
1, 3 |
|
|
|
|
|
|
0, 10 |
1, 11 |
| 11 |
|
|
|
1, 4 |
|
1, 6 |
|
|
|
1, 10 |
0, 11 |
Now relation theory tells us that in order to
find the transitive closure of the "connected" relation,
i.e., all shortest paths, we just have to multiply this matrix
with itself. CM2 contains all path of length 2, CMn
all paths of length n and CM*=CMoo (oo=infinity)
contains all paths. CM* is called the transitive closure of
CM. Luckily, we don't really have to multiply the matrix an
infinite number of times, because in a set of N polygons the
longest "shortest path" has N steps, so we only have
to calculate CMN (which is the same as CM* in this
case).
There is one little problem left: The entries
of CM are non-numbers, so how do we actually multiply the matrices?
Again, Algebra enters the stage to save the day. We just have
to define meaningful operators "*" (multiplication)
and "+" (addition) and a constant "0" for
the data type MatrixEntry. In this case, multiplication is the
operator that converts a path from i to j and a path from j
to k into a path from i to k (concatenation), and addition is
the operator that selects the shorter of two paths. They work
like this:
MatrixEntry::MatrixEntry(int i) // We define this as a constructor,
// so we can write MatrixEntry me=0
{
if(i==0) // The constant "0" stands for the "unknown path"
{
pathLength=INFINITY;
nextPolygon=-1;
}
else // Only allowed for i==0!
exit(1);
}
MatrixEntry operator*(const MatrixEntry& e1,const MatrixEntry& e2)
{
/* In this function we have to fudge a bit to handle */
/* both undefined and trivial paths correctly: */
MatrixEntry result;
if(e1.pathLength==INFINITY||e2.pathLength==INFINITY)
{
/* If one entry is undefined, the result is undefined, too: */
result.pathLength=INFINITY;
result.nextPolygon=-1;
}
else if(e1.pathLength==0)
{
/* Entry e1 is trivial, so take only e2: */
result=e2;
}
else if(e2.pathLength==0)
{
/* Entry e2 is trivial, so take only e1: */
result=e1;
}
else
{
/* Concatenate e1 and e2 by adding their lengths and going */
/* in direction of entry e1 first: */
result.pathLength=e1.pathLength+e2.pathLength;
result.nextPolygon=e1.nextPolygon;
}
return result;
}
MatrixEntry operator+(const MatrixEntry& e1,const MatrixEntry& e2)
{
if(e1.pathLength<=e2.pathLength) // Return the shorter path
return e1;
else
return e2;
}
If this was a mathematical paper, I would now
have to prove that the data type MatrixEntry, together with
the two operators, forms a "ring," that is a data
type where multiplication and addition behave in the same way
as they do in the prototype ring, the set of integers. Well,
actually, they do. (With one exception: a*b (way from a to b)
is, of course, not equal b*a (way from b to a). But this is
no problem, since multiplication in a ring is not required to
be commutative.) So let's move on.
The algorithm to multiply two matrices of our
type looks exactly like the multiplication algorithm for real
matrices:
MatrixEntry result[N][N],mat1[N][N],mat2[N][N];
for(int i=0;i<N;++i)
for(int j=0;j<N;++j)
{
MatrixEntry sum=0; // Good we defined that constructor!
for(int k=0;k<N;++k)
sum=sum+mat1[i][k]*mat2[k][j]; // "*" and "+" work as usual!
result[i][j]=sum;
}
It may be hard to see, but this algorithm actually
works! If mat1 contains all paths of length l1 (or shorter)
and mat2 contains all paths of length l2 (or shorter), result
will contain all paths of length (l1+l2) (or shorter) afterwards.
This is the magic of Algebra.
After one multiplication, the matrix CM for the
map of Figure 1 contains all paths of length two and looks like
this:
| |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
| 1 |
0, 1 |
1, 2 |
2, 2 |
|
|
|
|
|
|
|
|
| 2 |
1, 1 |
0, 2 |
1, 3 |
|
|
|
2, 3 |
|
|
2, 3 |
|
| 3 |
2, 2 |
1, 2 |
0, 3 |
|
|
2, 7 |
1, 7 |
2, 7 |
|
1, 10 |
2, 10 |
| 4 |
|
|
|
0, 4 |
1, 5 |
2, 11 |
|
|
|
2, 11 |
1, 11 |
| 5 |
|
|
|
1, 4 |
0, 5 |
|
|
|
|
|
2, 4 |
| 6 |
|
|
2, 7 |
2, 11 |
|
0, 6 |
1, 7 |
2, 7 |
|
2, 11 |
1, 11 |
| 7 |
|
2, 3 |
1, 3 |
|
|
1, 6 |
0, 7 |
1, 8 |
2, 8 |
2, 3 |
2, 6 |
| 8 |
|
|
2, 7 |
|
|
2, 7 |
1, 7 |
0, 8 |
1, 9 |
|
|
| 9 |
|
|
|
|
|
|
2, 8 |
1, 8 |
0, 9 |
|
|
| 10 |
|
2, 3 |
1, 3 |
4, 11 |
|
2, 11 |
2, 3 |
|
|
0, 10 |
1, 11 |
| 11 |
|
|
2, 10 |
1, 4 |
2, 4 |
1, 6 |
2, 6 |
|
|
1, 10 |
0, 11 |
So the overall algorithm looks like this:
Initialize CM
Set CMResult = CM
for(int iteration=1;iteration<N;++iteration)
CMResult=CMResult*CM; // Use our matrix multiplication here
CMResult contains all paths up to length N
Let's now analyse how long this iteration takes.
A single matrix multiplication has three nested loops of length
N, and hence takes O(N3) steps. Iterating for N steps
results in O(N4) steps altogether. Initializing the
matrix takes O(N2) steps, as does finding all adjacent
polygons and decomposing the original map polygon into complex
parts. This means our pre-processing algorith has a total runtime
of O(N4) and a total space requirement of O(N2)
to store the matrices CM and CMResult.
And now for the good part: Though this "naïve"
algorithm to find the connection matrix works, it is quite slow.
For 200 polygons (a reasonable number), it takes about 2004=1.6
billion steps, and that can take a while. But, by employing
one more Algebraic method, we can gain a considerable speedup
(factor 26 for N=200)! Here's how: Algebra tells us that taking
the N-th power of anything takes not N, but only ld(N) steps
if it is done correctly (ld: binary logarithm). One can rewrite
the expression CMN as (...((CM2)2)...)2
(m squares), when N=2m. Applying this method to our
algorithm yields:
Initialize CM
for(int iteration=1;iteration<N;iteration*=2)
CM=CM*CM; // Use our matrix multiplication here
CM contains all paths up to length N
If N is not a power of two, this computes a higher
power of CM, but it doesn't matter, since relation theory tells
us that CMT=CMN, if T is bigger than N.
So the application of results from matrix theory to our graphical
problem reduces the overall runtime to O(N3·log(N)),
which is 26 times faster than N4 for N=200.
Now the last missing part is how to use the closure
of the connection matrix when searching for a path in the game.
Let's recall the information that is stored in the matrix. Each
entry CM[i][j] was constructed to show the length of the way
from polygon i to polygon j, and which polygon to go to next.
So, if we want to go from polygon a to polygon b, we look up
CM[a][b] in the matrix and go to the next polygon it shows.
Now we are at polygon c, say, and still have to go to b. Easy.
Looking up CM[c][b], we find the next step to go. We repeat
this until we finally reach polygon b, and we know that this
algorithm will finish after CM[a][b].pathLength steps. In C++,
it looks like:
while(currentPolygon!=destinationPolygon)
{
int nextPolygon=CM[currentPolygon][destinationPolygon];
/* Move the game figure from polygon currentPolygon to nextPolygon */
currentPolygon=nextPolygon;
}
The one thing I didn't address in this paper
is how to move from one polygon to the next. This is easy if
we save the endpoint coordinates of all edges which are shared
by two polygons. Then to go from polygon a to the directly adjacent
polygon b, we could walk from the center of a to the center
of the common edge, and then from there to the center of b.
This results in decent walking behaviour. Another idea, that
leads to shorter walks but is a bit more difficult to implement,
walks from the current position to the center of the first crossed
edge, from there to the center of the second crossed edge, and
finally from the center of the last crossed edge to the destination.
See Figure 2.

Figure 2: Walking through a map. Left: Walking through polygon
centers; right: walking from edge to edge.
Back |
|