  We present a new flexible wavefront propagation algorithm for the boundary
value problem for sub-Riemannian (SR) geodesics in the roto-translation group
$SE(2) = \mathbb{R}^2 \rtimes S^1$ with a metric tensor depending on a smooth
external cost $\mathcal{C}:SE(2) \to [\delta,1]$, $\delta>0$, computed from
image data. The method consists of a first step where a SR-distance map is
computed as a viscosity solution of a Hamilton-Jacobi-Bellman (HJB) system
derived via Pontryagin's Maximum Principle (PMP). Subsequent backward
integration, again relying on PMP, gives the SR-geodesics. For $\mathcal{C}=1$
we show that our method produces the global minimizers. Comparison with exact
solutions shows a remarkable accuracy of the SR-spheres and the SR-geodesics.
We present numerical computations of Maxwell points and cusp points, which we
again verify for the uniform cost case $\mathcal{C}=1$. Regarding image
analysis applications, tracking of elongated structures in retinal and
synthetic images show that our line tracking generically deals with crossings.
We show the benefits of including the sub-Riemannian geometry.
