Solving the quadratic in the Fast Marching Method

Question

The Fast Marching Method is a way of solving the Eikonal Equation on a discrete grid, essentially just computing a wavefront speading out from initial points, IE:

The idea is that we want to compute the time at which a wavefront, moving across a field with a velocity function $F(x, y)$, reaches each point in a grid, IE:

This is done by, as described in the "Implementation Details of the Fast Marching Methods" link of the Wikipedia page, setting all beginning nodes as frozen, with T=0. Then every node adjacent (4 neighbor neighborhood) 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.

Using these newly computed $T(x, y)$ values, each narrow band node is added to a stack, then the one with the lowest $(T(x, y))$ is removed, set as "frozen," and this process is repeated until all accessible nodes are "frozen." Because this also means some "narrow band" nodes are recomputed with more information as more of their neighbors become frozen, when this re-computation occurs they will simply be added to the stack again, and if a previously "narrow band" node is removed from the stack that is now frozen it is just ignored.

This process typically works very well, and is widely used in many fields.

However, when I tried implementing this myself, I solved the quadratic as follows:

$$firstSolution = \frac{\sqrt{-c^2(a^2c^2-2abc^2+b^2c^2-2)}+ac^2+bc^2)}{2c^2}$$

or

$$secondSolution = \frac{-\sqrt{-c^2(a^2c^2-2abc^2+b^2c^2-2)}+ac^2+bc^2)}{2c^2}$$

where

$$c = F(x, y)$$

and

$$T(x, y) = Max(firstSolution, secondSolution)$$

Yet, there are some cases where this method fails, for example, with $a = 53.7432442$, $b = 53.8721733$, and $c = 11.18034$, square root insides equal $-9.730039,$ meaning there is no (real) solution.

I thought this may just be because I was changing the value of f(x, y) too quickly (we're approximating a continuous function, and too sharp changes make it not differentiable), but it seemed to happen even when I was trying to change it as smoothly as possible.

Is there a reason for this, or a better way to approximate these T(x, y) values? Or is this just an inherent problem with this method?

The article linked to in the Wikipedia page suggested second order approximations, but I don't think this solves the problem. I'm sorry if this is the wrong place to post such a question as well (being somewhat implementation oriented), feel free to migrate as needed.

Answer 1

Physically, this problem arises when the two input nodes are "impossibly far apart" according to the planar wavefront approximation. Specifically, if two nodes separated by distance $\Delta x \sqrt{2}$ on the grid have $T$ values differing by more than $\Delta x \sqrt{2}$, there is no solution plane that fits both distance values.

Even if the quadratic does have a solution, it might be the case that two inputs $T_0 < T_1$ result in an output value $T_2$ with $T_0 < T_2 < T_1$. This is physically invalid, as it corresponds to information flowing backwards through time.

The solution in both cases is to fall back to a lower dimension wave front, ignoring the value that is too large, and corresponding to an axis aligned wavefront approximation. This can be implemented branch free. If we assume $\Delta x = 1$ and that all distance values are positive, the example code is

double step(double px, double py) {
  small_sort(px,py);
  py = min(py,px+1);
  return .5*(px+py+sqrt(2-sqr(px-py)));
}

double step(double px, double py, double pz) {
  small_sort(px,py,pz);
  py = min(py,px+1);
  pz = min(pz,.5*(px+py+sqrt(2-sqr(px-py))));
  return 1./3*(px+py+pz+sqrt(3-sqr(px-py)-sqr(py-pz)-sqr(pz-px)));
}

Here small_sort sorts two or three numbers, and each min operation implements falling back from 3D to 2D or 2D to 1D planar wavefront approximations.