How does Flow Field pathfinding work?

Question

Supreme Commander 2 has something called flow field path finding. How does it work? Is there some article available I can read up on how it works?

Answer 1

I wrote flow fields for sup com 2, and I wrote an article explaining the details. It can be found in the upcoming book "Game AI Pro: Collected Wisdom of Game AI Professionals".

Also, I recently did a video stream talking about flow fields for Planetary Annihilation. I show some debug views and explain how it works at a high level. http://youtu.be/5Qyl7h7D1Q8?t=24m24s

Hope this helps

Answer 2

I've been looking for this term as well, and this paper is the only major one I could find that references flow fields directly: http://www.aaai.org/Papers/AIIDE/2008/AIIDE08-031.pdf

This approach involves each pathfinding agent being influenced by a global vector field, and in turn influencing that field with their resulting path. You still need some basic object avoidance code to kick things off, so is only really applicable for swarms crossing paths, rather than individual agents.

However, the SupComm2 guys mention research from the University of Washington directly, and this is the most applicable paper that I could find from that institution: http://grail.cs.washington.edu/projects/crowd-flows/

This approach seems more promising but I need to read more about it.

Answer 3

The official specs of the pathfinding algorithm in Supreme Commander 2 refer to the Crowd Continuum Study at University of Washington.

There are several papers and demonstrations at that link.

Answer 4

I'll include the description of how to solve an eioknal equation here, because it's the messy bit to implementing flow field path finding. See Continuum Crowds for an example of how to use this solution for path finding. Basically, the idea is that you are computing the time a wavefront spreading out from initial points reaches every other point in a grid, e.g.:

(Image from Implementation Details of the Fast Marching Methods, which I highly recommend reading)

It's sorta like pathfinding starting from a set of goal points, and computing the amount of time T(x, y) it would take to reach the nearest goal point from any point on a grid if you were following the optimal path and moving at speed F(x, y) (the F(x,y) values you can chose and are part of the map, and something like terrain roughness could make them lower). This is nice because you are solving for every point on a map, so every agent with the same goal can use this same result, and you only have to do one run of the algorithm (per timestep) for all of them. To get these T(x, y) values, there are three options I know of: The Fast Marching Method, Fast Sweeping Method, and Fast Iterative Method.

The Fast Marching Method works fairly well for small grids, takes time O(nlog(n)), where n is the number of points in the grid, and is best described in the "Implementation Details of the Fast Marching Methods" link at the bottom of the wiki page. However, for larger grids (say in 3D), this log(n) factor becomes prohibitive, so that's where the Fast Sweeping Method and Fast Iterative Method come in. More importantly, the Fast Marching Method uses a binary heap, which means it's not straightforward to parallelize, and has poor cache coherency. The Fast Sweeping Method has the potential to perform better, yet is slowed down by any direction changes that could happen, say, in the case of a winding maze. The Fast Iterative Method is a pretty good compromise between these two, working pretty well in winding mazes, being somewhat amenable to parallelization, and having good cache coherency. It also is slower in smaller maps, however, so I'd probably recommend implementing the Fast Marching Method first, then if it's too slow working towards implementing the Fast Iterative Method in something like OpenCL using Scans.

I'll describe the Fast Marching Method in more detail here, as the part about how to solve the quadratic is equally applicable to the Fast Iterative Method. In short, the Fast Marching Method starts by setting all beginning nodes as frozen, with T(x, y)=0. Then every node adjacent (4 neighbor neighborhood in 2D and 6 neighbor neighborhood in 3D) to these nodes is marked as being in the "narrow band," and their T(x, y) is computed using a quadratic:

$$(T(x, y) - a)^2 + (T(x, y) - b)^2 = (1/F(x, y))^2$$

Where

$$a = Min(T(x + 1, y), T(x - 1, y))$$

$$b = Min(T(x, y + 1), T(x, y - 1))$$

only considering values for T(...) of neighboring frozen nodes. Thus, for example, if the node above and the node below aren't frozen, the (T(x, y) - b)^2 term drops out of the equation. This can be solved by hand, or with wolfram, and you use the greater of the two solutions. In 3D it's very similar. The reason for this quadratic is because it makes the solution look smoother than simply spreading out linearly from the goal, also for math consistency reasons I don't fully understand.

When implementing this there were times where this quadratic had no solution, however, so I asked about this on the scicomp stack exchange, and got a very good answer that's best explained there: basically in that case you just need to solve the quadratic differently by falling back to a "lower dimensional wavefront." You also need to do this if the solution you get is lower than any frozen node's T(x, y) values nearby - because then the quadratic solution is invalid because if you use it the wave is somehow travelling back in time.

You add each of these narrow band nodes to a binary heap, sorted by their newly computed T(x, y) values. Then the node at the top is removed, made a frozen node, and any new narrow band nodes are added to the heap. It's important to note that this means narrow band nodes can be added to the heap more than once, which can be addressed by simply ignoring any frozen nodes that are at the top of the heap - since we've already processed them with a more accurate T(x, y) value before.

This process is repeated until all nodes are set to frozen (or inaccessable), at which case you're done solving the eikonal equation. This is also cool because you can easliy make boundaries by either never allowing boundary nodes to be added to the "narrow band" so they never become frozen, or by giving them very low F(x, y) values so agents prefer to go around them. Just be warned that the falling back to a "lower dimensional wavefront" (as described in the scicomp stack exchange link above) is needed if your F(x, y) values change too quickly, as they would with arbitrary boundaries.

Finally, you need to be able to compute

$$\nabla{T}$$

which is a vector with an x and y component (I'll call them dx and dy) that describes the rate of change in their respective directions. These are defined as:

$$dx(x, y) = Average(T(x + 1, y) - T(x, y), T(x, y) - T(x - 1, y) )$$

$$dy(x, y) = Average(T(x, y + 1) - T(x, y), T(x, y) - T(x, y - 1))$$

which continuum crowds also has you normalize (divide each dx(x, y) and dy(x, y) by $$\sqrt{dx(x, y)^2 + dy(x, y)^2}$$ unless they're both zero). If a neighboring node doesn't exist (say if it's outside the boundary or within an obstacle), you can just use the other neighboring node's value and not compute any average, and if neither neighbor node exists it probably makes sense to use some default value like 0. This gives us vectors which point in the direction of the optimal path to the goal, and you can move units along these optimal paths by moving them in this direction at a velocity equal to F(x, y).

Since we're working with a discrete grid, however, all units will be between four points. To find the value of a function at any point we need to interpolate between the four surrounding values:

float interpolateBetweenValues(float x, float y, float[,] array)
{
    int topLeftX = (int)Math.Floor(x);
    int topLeftY = (int)Math.Floor(y);

    if (topLeftX < 0 || topLeftX + 1 >= array.GetLength(0) ||
        topLeftY < 0 || topLeftY + 1 >= array.GetLength(1))
        throw new Exception("Out of bounds");

    float xAmountRight = x - topLeftX;
    float xAmountLeft = 1.0f - xAmountRight;
    float yAmountBottom = y - topLeftY;
    float yAmountTop = 1.0f - yAmountBottom;

    float averagedXTop = array[topLeftX, topLeftY] * xAmountLeft + array[topLeftX + 1, topLeftY] * xAmountRight;
    float averagedXBottom = array[topLeftX, topLeftY+1] * xAmountLeft + array[topLeftX + 1, topLeftY+1] * xAmountRight;

    float averagedYTotal = averagedXTop * yAmountTop + averagedXBottom * yAmountBottom;

    return averagedYTotal;
}