The inverse gradient problem, finding a scalar field f with a gradient near a given vector field g on some bounded and connected domain Rn , can be solved by means of a Poisson equation with inhomogeneous Neumann boundary conditions. We present an elementary derivation of this partial differential equation and an efficient multigridbased method to numerically compute the inverse gradient on non-rectangular domains. The utility of the method is demonstrated by a range of important medical applications such as phase unwrapping, pressure computation, inverse deformation fields, and fiber bundle tracking.