Okay, so now you have your happy eikonal equation solver and you're ready to use it in Continuum Crowds. To be clear you don't actually

*need*any of the stuff here, an eikonal equation solver often simpler than the one given in the previous article possibly combined with RVO works for this kind of flow field pathfinding, but if you want something more like this or this that's what this post is for.

The first thing to note is that, in the previous post, I described $\nabla{T}$ as having one x and y component for each grid node. In the paper, the author actually uses 4 components per node, two x for left and right, and two y for up and down. The method I described averages the x values and y values together, and though it's a little different than the paper, I prefer it. It makes for cleaner implementations, and in practice is basically identical.

The idea behind continuum crowds is that you start with a

*discomfort*and

*height*field, and a set of

*units*. From the positions of these units you compute a

*density*field, and from their velocities you compute a

*velocity*field. Using all these fields you compute a

*speed*field, solve an eioknal equation to find a

*potential*field, then compute the gradient of that solution and multiply by your speed field to find the final velocities of units. You update those units' positions and velocities, and repeat.

From Continuum Crowds:

So first, every field here is the same width and height. Like the eikonal solver from last post, the entire algorithm is running on every grid point every timestep, because that allows for a single computation to be reused for every unit that has the same goal.

The discomfort field $g$ (I have no idea why the choice of 'g') is full of positive floating point values: The higher the value, the less units like to go there. This allows units to prefer things like a paved trail if you give it a lower discomfort than the terrain surrounding it.

The height field $h$ is full of floating point values, and sorta acts like the discomfort field: People prefer to go downhill, aren't affected by no change in slope, and avoid going uphill when possible.

Both of these fields are inputs to Continuum Crowds, and aren't modified by the algorithm - though they can be "dynamic" and changed between timesteps if the caller wants to.

Next, we have the density field $\rho$ (also called the potential field). This field also acts like the discomfort field, helping people avoid places of high density and prefer places of low density. $\bar{\rho}$ is a positive floating-point input parameter between 0.0 and 1.0 (0.7 is fine), where each unit will contribute a minimum of $\bar{\rho}$ to it's own cell, and a maximum of $\bar{\rho}$ to it's neighbor cells. The individual potential field $\rho_i$ for a individual $unit=units_i$ is then created through an interpolating process basically the reverse of the one described in the previous post:

$topLeftX = (int)floor(unit_{posx})$

$topLeftY = (int)floor(unit_{posy})$

$\Delta{x} = x - topLeftX$

$\Delta{y} = y - topLeftY$

$density(topLeftX, topLeftY) = min(1 - \Delta{x}, 1 - \Delta{y})^\lambda$

$density(topLeftX + 1, topLeftY) = min(\Delta{x}, 1 - \Delta{y})^\lambda$

$density(topLeftX, topLeftY + 1) = min(1 - \Delta{x}, \Delta{y})^\lambda$

$density(topLeftX + 1, topLeftY + 1) = min(\Delta{x}, \Delta{y})^\lambda$

Where $\bar{\rho}={1/2}^\lambda$, so solving for $\lambda$ gives $\lambda=log_{1/2}{\bar{\rho}}$

You sum each individual potential field $\rho_i$ to get the potential field $\rho$, IE:

$$\rho= \sum\limits_{i}\rho_i$$

The velocity field $\bar{v}$ can be created at the same time, by multiplying each potential field by that unit's velocity before summing them all together. After summing all the individual velocity fields together the velocity field $\bar{v}$ is also divided by the total potential field $\rho$ at that point. In summary:

$$\bar{v_x}= \frac{\sum\limits_{i}{\rho_i*unit_{velx}}}{\rho}$$

Where $\bar{\rho}={1/2}^\lambda$, so solving for $\lambda$ gives $\lambda=log_{1/2}{\bar{\rho}}$

You sum each individual potential field $\rho_i$ to get the potential field $\rho$, IE:

$$\rho= \sum\limits_{i}\rho_i$$

The velocity field $\bar{v}$ can be created at the same time, by multiplying each potential field by that unit's velocity before summing them all together. After summing all the individual velocity fields together the velocity field $\bar{v}$ is also divided by the total potential field $\rho$ at that point. In summary:

$$\bar{v_x}= \frac{\sum\limits_{i}{\rho_i*unit_{velx}}}{\rho}$$

$$\bar{v_y}= \frac{\sum\limits_{i}{\rho_i*unit_{vely}}}{\rho}$$

Okay, so now we have all those fields down, there's only two left: the potential field $\phi$ and the speed field $f$.

The speed field $f$ is the most painful, so I'll start with that. You can think of this field as being the same as $F$ in the previous post, because it will be used the same way for solving the eikonal equation and moving the units with the normalized final solution. Essentially the idea is we have a

So first, we need to compute the gradient of the height field $\nabla{h}$, which is done in the same way as $\nabla{T}$ in the previous post. This essentially gives us the slope (as an x and y component) for each node. Now we want a speed field that is slow (low values) when units are traveling uphill, and fast (high values) when units are travelling downhill. We'll pretend all units are moving in the positive x and y directions (herby referred to as right and down), and then we can make a small tweak for our result to work in all directions.

Looking at x first, we need to figure out whether a positive gradient values mean we're travelling up or downhill. Because, according to the definition in the previous post, a positive gradient value means that the slope is increasing (as you move to the right), a positive gradient = going uphill. Thus, for higher slope values we want to have lower speeds, and for lower slope values we want to have higher speeds. This is especially true for negative gradient values: that means we're going downhill, so the slope should

Specifically we want a function that returns 1 when the gradient equals 0, returns a value greater than 1 when the gradient is less than 0, and returns a value between 0 and 1 when the gradient is greater than 0. An exponential like $e^{-x}$ works, however it shrinks and grows a little to slowly/quickly. Really we just want a line that passes through 1 at x = 0, and caps off at some minimum speed, I think. I believe this is the reason for using the min and max slope and speed values, and it would probably be nice to have a max speed value, so I'll guess I'll conform. The above logic is also true for y, so for some minimum and maximum speed $f_{min}$ and $f_{max}$, and for some minimum and maximum slope $s_{min}$ and $s_{max}$:

$$f_{T_{x}}(x, y)=f_{min}+\frac{(-\nabla{h_x}(x, y)-s_{min})}{(s_{max}-s_{min})}(f_{max}-f_{min})$$

Where $f_{min}$ and $f_{max}$ would be something like 0.1 and 5.0, and $s_{min}$ and $s_{max}$ would be something like -10.0 and 10.0. This equation works, as when the slope equals $s_{min}$ it will result in $f_{max}$, when the slope equals $s_{max}$ it will result in $f_{min}$, and anything between gives a result between as well.

Next we need to find the flow speed $f_{\bar{v}}$. This uses the velocity field $\bar{v}$ computed above, specifically the velocity field

Finally, given some minimum and maximum densities $\rho_{min}$ and $\rho_{max}$, we'll interpolate between flow speed and topological speed:

$$f_x(x, y) = f_{T_x}(x, y) + \frac{\rho(x + 1, y + 1) - \rho_{min}}{\rho_{max}-\rho_{min}}(f_{\bar{v_x}}(x, y)-f_{T_x}(x, y))$$

$$f_y(x, y) = f_{T_y}(x, y) + \frac{\rho(x + 1, y + 1) - \rho_{min}}{\rho_{max}-\rho_{min}}(f_{\bar{v_y}}(x, y)-f_{T_y}(x, y))$$

The intuition behind this is that when moving into higher densities, we want movement to be dominated by the overall crowd direction there, while at lower crowd densities we want movement to be mostly dominated by the terrain units are moving from. $\rho_{min}$ and $\rho_{max}$ aren't hard limits either: if $\rho(x + 1, y + 1)<\rho_{min}$ then $f(x, y)= f_T(x, y)$, if $\rho(x + 1, y + 1)>\rho_{max}$ then $f(x, y)= f_{\bar{v}}(x, y)$.

If we are moving to, say, $(x - 1, y - 1)$, we would simply look up our $\bar{v_x}$ and $\rho$ values there instead of $(x + 1, y + 1)$. We would also use $\nabla{h}$ instead of $-\nabla{h}$ and $-\bar{v}$ instead of $\bar{v}$, because we're now moving in different directions so uphill and the "same" movement direction are different. You do similar things for $(x - 1, y + 1)$ and $(x + 1, y - 1)$.

This gets back to what I was getting at at the start of this post: At this point, the paper asks you to find $f(x, y)$ for all 4 possible directions ($(x - 1, y - 1)$, $(x + 1, y - 1)$, $(x - 1, y + 1)$, and $(x + 1, y + 1)$), then use them as needed. To understand why this is needed, it's worth talking about how the potential function $\phi$ is used and computed.

Formally, the potential function is written in terms of the eikonal equation:

$$||\nabla{\phi}(x, y)||=C(x, y)$$

Where

$$C(x, y) = \frac{f(x, y)+1+g(x, y)}{f(x, y)}$$

evaluating $f(x, y)$ in the direction of the optimal path to the goal.

What does this mean? In terms of the previous post, $C(x, y)$ is the same as $F(x, y)$. The key here is that you only ever evaluate $F(x, y)$ when you're solving the 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))$$

Look closely at $a$ and $b$ here: The value that we choose tells us what direction we're going to move. Because we always move in the direction where T(x, y) decreases most quickly (because that gets us to the goal fastest), if, say, $T(x + 1, y)$ and $T(x, y - 1)$ are smallest, that means we are moving in direction $(x + 1, y -1)$. Now, if when solving the top quadratic we're using $C(x, y)$ instead of $F(x, y)$, this means we evaluate $f(x, y)$ in the direction $(x + 1, y - 1)$ at that point.

Once you have your solution $T(x, y)$, $\phi(x, y)=T(x, y)$ so you're already done computing $\phi$! All that remains is to compute $\nabla\phi$, normalize it, and multiply by $f(x, y)$ at each unit's position to get the final velocity. You look at the sign of $\nabla\phi_x(x, y)$ and $\nabla\phi_y(x, y)$ to tell whether you need to evaluate $f(x, y)$ in direction $(x - 1, y - 1)$, $(x + 1, y - 1)$, $(x - 1, y + 1)$, or $(x + 1, y + 1)$). Then move all units according to those computed velocities, repeat, and your continuum crowds implementation will be complete :) If you want multiple different "groups," each with different goals, simply have one density and velocity field for all groups, and otherwise compute speed and potential fields independently per group.

Okay, so now we have all those fields down, there's only two left: the potential field $\phi$ and the speed field $f$.

The speed field $f$ is the most painful, so I'll start with that. You can think of this field as being the same as $F$ in the previous post, because it will be used the same way for solving the eikonal equation and moving the units with the normalized final solution. Essentially the idea is we have a

*topological*speed $f_T$, which is how fast we can travel when considering slopes, and a*flow*speed $f_{\bar{v}}$, which is how much our movement is impeded by those around us.**This is the key idea in Continuum Crowds**: units are slowed down if they are moving against the flow, and sped up if they are moving with the flow. It causes the formation of lanes, vortices, and other emergent phenonoma that have been observed in real crowds.So first, we need to compute the gradient of the height field $\nabla{h}$, which is done in the same way as $\nabla{T}$ in the previous post. This essentially gives us the slope (as an x and y component) for each node. Now we want a speed field that is slow (low values) when units are traveling uphill, and fast (high values) when units are travelling downhill. We'll pretend all units are moving in the positive x and y directions (herby referred to as right and down), and then we can make a small tweak for our result to work in all directions.

Looking at x first, we need to figure out whether a positive gradient values mean we're travelling up or downhill. Because, according to the definition in the previous post, a positive gradient value means that the slope is increasing (as you move to the right), a positive gradient = going uphill. Thus, for higher slope values we want to have lower speeds, and for lower slope values we want to have higher speeds. This is especially true for negative gradient values: that means we're going downhill, so the slope should

*speed us up*instead of slow us down. At this point the paper assumes a maximum and minimum slope and a maximum and minimum speed and uses a weird averaging thing accordingly, but I'd prefer to be working with any possible value.Specifically we want a function that returns 1 when the gradient equals 0, returns a value greater than 1 when the gradient is less than 0, and returns a value between 0 and 1 when the gradient is greater than 0. An exponential like $e^{-x}$ works, however it shrinks and grows a little to slowly/quickly. Really we just want a line that passes through 1 at x = 0, and caps off at some minimum speed, I think. I believe this is the reason for using the min and max slope and speed values, and it would probably be nice to have a max speed value, so I'll guess I'll conform. The above logic is also true for y, so for some minimum and maximum speed $f_{min}$ and $f_{max}$, and for some minimum and maximum slope $s_{min}$ and $s_{max}$:

$$f_{T_{x}}(x, y)=f_{min}+\frac{(-\nabla{h_x}(x, y)-s_{min})}{(s_{max}-s_{min})}(f_{max}-f_{min})$$

$$f_{T_{y}}(x, y)=f_{min}+\frac{(-\nabla{h_y}(x, y)-s_{min})}{(s_{max}-s_{min})}(f_{max}-f_{min})$$

Next we need to find the flow speed $f_{\bar{v}}$. This uses the velocity field $\bar{v}$ computed above, specifically the velocity field

*at the point we are moving into*. As described in the paper, "indeed, if not for this offset, a person's speed would be dominated by their own previous speed, an undesirable effect." What does that mean? Since we are assuming everyone is moving right and down, we look at $\bar{v}(x + 1, y + 1)$. If this velocity is going in the same direction as us, we want a high speed, if this velocity is going in the opposite direction as us, we want a low speed. We don't ever want to be pushed backwards though (by having a negative speed), so we'll clamp it at our $f_{min}$ from before.
$$f_{\bar{v}_x}(x, y)=max(f_{min}, \bar{v_x}(x + 1, y + 1))$$

$$f_{\bar{v}_y}(x, y)=max(f_{min}, \bar{v_y}(x + 1, y + 1))$$

$$f_x(x, y) = f_{T_x}(x, y) + \frac{\rho(x + 1, y + 1) - \rho_{min}}{\rho_{max}-\rho_{min}}(f_{\bar{v_x}}(x, y)-f_{T_x}(x, y))$$

$$f_y(x, y) = f_{T_y}(x, y) + \frac{\rho(x + 1, y + 1) - \rho_{min}}{\rho_{max}-\rho_{min}}(f_{\bar{v_y}}(x, y)-f_{T_y}(x, y))$$

The intuition behind this is that when moving into higher densities, we want movement to be dominated by the overall crowd direction there, while at lower crowd densities we want movement to be mostly dominated by the terrain units are moving from. $\rho_{min}$ and $\rho_{max}$ aren't hard limits either: if $\rho(x + 1, y + 1)<\rho_{min}$ then $f(x, y)= f_T(x, y)$, if $\rho(x + 1, y + 1)>\rho_{max}$ then $f(x, y)= f_{\bar{v}}(x, y)$.

If we are moving to, say, $(x - 1, y - 1)$, we would simply look up our $\bar{v_x}$ and $\rho$ values there instead of $(x + 1, y + 1)$. We would also use $\nabla{h}$ instead of $-\nabla{h}$ and $-\bar{v}$ instead of $\bar{v}$, because we're now moving in different directions so uphill and the "same" movement direction are different. You do similar things for $(x - 1, y + 1)$ and $(x + 1, y - 1)$.

This gets back to what I was getting at at the start of this post: At this point, the paper asks you to find $f(x, y)$ for all 4 possible directions ($(x - 1, y - 1)$, $(x + 1, y - 1)$, $(x - 1, y + 1)$, and $(x + 1, y + 1)$), then use them as needed. To understand why this is needed, it's worth talking about how the potential function $\phi$ is used and computed.

Formally, the potential function is written in terms of the eikonal equation:

$$||\nabla{\phi}(x, y)||=C(x, y)$$

Where

$$C(x, y) = \frac{f(x, y)+1+g(x, y)}{f(x, y)}$$

evaluating $f(x, y)$ in the direction of the optimal path to the goal.

What does this mean? In terms of the previous post, $C(x, y)$ is the same as $F(x, y)$. The key here is that you only ever evaluate $F(x, y)$ when you're solving the 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))$$

Look closely at $a$ and $b$ here: The value that we choose tells us what direction we're going to move. Because we always move in the direction where T(x, y) decreases most quickly (because that gets us to the goal fastest), if, say, $T(x + 1, y)$ and $T(x, y - 1)$ are smallest, that means we are moving in direction $(x + 1, y -1)$. Now, if when solving the top quadratic we're using $C(x, y)$ instead of $F(x, y)$, this means we evaluate $f(x, y)$ in the direction $(x + 1, y - 1)$ at that point.

Once you have your solution $T(x, y)$, $\phi(x, y)=T(x, y)$ so you're already done computing $\phi$! All that remains is to compute $\nabla\phi$, normalize it, and multiply by $f(x, y)$ at each unit's position to get the final velocity. You look at the sign of $\nabla\phi_x(x, y)$ and $\nabla\phi_y(x, y)$ to tell whether you need to evaluate $f(x, y)$ in direction $(x - 1, y - 1)$, $(x + 1, y - 1)$, $(x - 1, y + 1)$, or $(x + 1, y + 1)$). Then move all units according to those computed velocities, repeat, and your continuum crowds implementation will be complete :) If you want multiple different "groups," each with different goals, simply have one density and velocity field for all groups, and otherwise compute speed and potential fields independently per group.

Very interesting topic, and well explained !

ReplyDeleteI've noticed some inconsistencies though. In the first part, about calculating T(x,y) it is said that the solution to be considered is Max(solution1, solution2), whereas in the scicomp forum it is the min (I guess it is the min but...).

In the second part, I quote:

if ρ(x+1,y+1)ρmax then f(x,y)=fT(x,y).

Isn't it the opposite ?

My copy paste didn't succeed...

ReplyDeleteif ρ(x+1,y+1)ρmax then f(x,y)=fT(x,y).

Thanks :) The best resource for explaining that is the "Implementation Details of the Fast Marching Methods" in the Fast Marching Method Wikipedia page. I'm mostly referring you to there because, for anyone that is implementing this seriously, the method I described here can be slightly off, and they describe this method here first, and then a better method that's also more complicated to implement. They also describe why both methods work, which I didn't explain here either.

ReplyDeleteThat aside though, if you get two solutions from your quadratic, you want the Max(solution1, solution2), not the min. I fixed my scicomp question accordingly. This is because the current node must be larger than it's neighbors, and the smaller solution won't fulfill this criteria.

On your other question, you're very right, sorry. If our density is lower than pMin we should only use the topological speed, while if it's greater than pMax we should only use the velocity field's speed. That is now fixed.

To be clear, the method I gave does work for, say, pathfinding. It just isn't "perfect" in computing T(x,y), meaning that if you are depending on the exact behavior of this function for an application you need to tweak this method a little to get as close to the exact behavior as possible.