Finding shortest path in a time/distance map

I get a distance map output after using a Fast Marching Method. The PDE involved is the Eikonal equation which take the form : $$\begin{cases} c(x).|\nabla u| =1\\ u(x) =\phi(x) \end{cases}$$ where c is a speed function and u the unknown distance/time we compute. Below an example of a 2D distance map (in blue the inital points and in red the further points). Given 2 points in the map I would like to show a simple/shortest path between them. I have tried using a naive algorithm which is not very good since I can have complicated problems with obstacles, multiple sources, different speeds.

find_path(point A, point B)
{
point current_point = B;   // B is ending point
while(current_point is not A)
{
new_current_point is the the neighbor of current_point with the minimum value   // I use a 5 point stencil
current_point = new_current_point
tag every current_point found to show the path
}
}

My data structure is an uni dimensional row major array (I am coding in C) representing distance values in a 2D regular grid. I think gradient descent might be a good solution but I have difficulties finding a simple algorithm. Thank you.

Edit I join an article Application of the fast marching method for outdoor motion planning in robotics (starting page 12) which use a gradient descent in order to compute that. I have already seen that other papers but I cannot see how it works.

Update I am trying to add a backtracking algorithm to my naïve one but I have troubles putting it in applications. I feel like I am doing a brute force and this is not what I want to do. If someone has an idea on how to change my naive algorithm.

• I'm not sure your question makes sense (to me): You compute the distance - presumably given a (nonuniform) metric - from a fixed point to any other point in your domain, and then want to use that distance map as a new metric for computing the distance between two different points? Can you add some context why you are doing this? – Christian Clason Aug 14 '14 at 9:35
• Yes I will edit my post for examples. The distance metric is uniform for every point. The value array can be seen as the shortest time the initial point reaches an other point in the grid. – coincoin Aug 14 '14 at 11:06
• If your distance metric is uniform, then I don't understand the map you show -- the distance map should be uniform as well (in particular, radially symmetric around the initial point). – Christian Clason Aug 14 '14 at 11:09
• The propagation is not uniform. This can be seen as a wave propagation where different speed can apply on the points. Maybe I am wrong about the meaning of distance metric... – coincoin Aug 14 '14 at 11:11
• Yes, the distance metric exactly corresponds to (one over) the wave speed (squared). – Christian Clason Aug 14 '14 at 11:12

If you already have a distance map $T(x)$ that describes the distance (with respect to a given metric $d$) from a given point $x_0$ to any other point $x$ in your domain (i.e., $d(x_0,x) = T(x)$), then the shortest path from $x_0$ to any point $y$ -- called the geodesic between $x_0$ and $y$ -- can be computed by backtracking: Set $X(0) = y$ and then solve (e.g., by explicit Euler) the gradient flow $$\dot X(s) = -\nabla T(X(s)), \qquad s>0,$$ where $\nabla T$ is the gradient of $T$, until for some $S>0$ you have $X(S) = x_0$.

Your algorithm can be seen as a space-discrete variant of this, but it will only ever compute the shortest path from B to the "bluest" point in your distance map -- if that is not A, it will not terminate in general.

(Note that $T$ depends on you initial point $x_0$, and has to be recomputed for each initial point. There are more efficient methods for computing all possible pairs of geodesics at once; this is sometimes referred to all-to-all geodesics, see, e.g., http://www.ias-iss.org/ojs/IAS/article/viewFile/891/789.)

• Thank you Christian for your answer. I think this is the idea to use backtracking technique. I start from y and then I try to find x_0. You say I have to recompute T if I change x_0 but I think it is possible with gradient descent to keep my map and change x_0 and y in the map (see article linked in my edited post). – coincoin Aug 14 '14 at 11:16
• Unless you have an explicit formula for the map $T$ in terms of $x_0$ (which is not what you get from fast marching), then I don't see how you can change $x_0$ without starting over with fast marching. From what I can tell, there's no mention in the linked paper about changing the initial point (which enters the Eikonal equation (4) as a boundary condition, i.e., enforcing $T(x_0)=0$). – Christian Clason Aug 14 '14 at 11:29
• Yes sorry you are right – coincoin Aug 14 '14 at 11:45
• Thanks Christian I took the liberty to join my results with my algorithm above and I am pretty sure it won't work in all cases. I did not implement any backtracking method maybe that is the solution ? – coincoin Aug 14 '14 at 15:20
• No, you are implementing a backtracking if you start with B. Your new pictures don't look unreasonable to me; the path is not smooth, but there's nothing enforcing this. (If I understand it correctly, this is why the reference you linked introduced the matrix $W$.) Anyway, you didn't mention anything about obstacles before; now it's really a different question... – Christian Clason Aug 14 '14 at 16:08

This can be solved using some variant of Dijkstra's algorithm, such as the A* algorithm.

• Thank you moyner but since we have some informations already with the distance map I think it is costly to use over a A* algorithm. – coincoin Aug 14 '14 at 11:04