We present a novel way to efficiently compute Riemannian geodesic distance over a two- or three-dimensional domain. It is based on a previously presented method for computation of geodesic distances on surface meshes. Our method is adapted for rectangular grids, equipped with a variable anisotropic metric tensor. Processing and visualization of such tensor fields is common in certain applications, for instance structure tensor fields in image analysis and diffusion tensor fields in medical imaging. The included benchmark study shows that our method provides significantly better results in anisotropic regions in 2-D and 3-D and is faster than current stat-of-the-art solvers in 2-D grids. Additionally, our method is straightforward to code; the test implementation is less than 150 lines of C++ code. The paper is an extension of a previously presented conference paper and includes new sections on 3-D grids in particular.