TWIGZ Path Finding

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:

  1. It is reflexive,i.e., every polygon is connected to itself (trivial).
  2. 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.
  3. 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

© 2004 Shattered Realm Productions