A better Fast Marching Method?

Question

I am using the Fast Marching Method (FMM) to calculate shortest "distance" (traveltime) from some points.

The way FMM works is: I keep a velocity function in RAM: V(xi,yj,zk). I also keep a priority queue of all points on the front sorted on their V value. I repeatedly propagate the first point in this queue one step outwards. As I do this I remove the used point from the queue and insert those "touched" by the front.

My current implementation have two problems:

I. I must keep the whole cost (slowness) function in RAM. This limits the size of the cost function I can use.

II. I would like it to be even faster.

Any suggestions on how to improve my current implementation? For instance would it be possible to implement this on the GPU?

Answer 1

My current implementation have two problems:

I. I must keep the whole cost (slowness) function in RAM. This limits the size of the cost function I can use.

II. I would like it to be even faster.

The first point is a bit ambiguous. Just computing the cost function and the traveltime for each grid cell and storing them somewhere should not be an issue. The issue is the access pattern of the fast marching method to this data, which is very cache unfriendly and prevents effective parallelization. For a straightforward implementation, the priority queue is the worst offender for the cache. There are cache optimal priority queues, but even then the memory access pattern of the updates is still very irregular.

The second point is probably most related to the parallelization potential of the algorithm. At least for the second point, one is probably forced to give up on fast marching methods and look at fast sweeping methods instead.

The observation behind the fast sweeping methods is that for the updates it is only important in which of $2^n$ "general directions" the characteristics point. Especially any region where this "general direction" is constant can be updated by a sweep. Still, there are many different ways to turn this observation into an algorithm, and which strategy will be most efficient depends a bit on whether your velocity function leads to large regions with constant "general direction", or whether it creates more a sort of labyrinth.