Implicitly embedding a surface as a level set of a scalar function φ : ޒ d → ޒ is a powerful technique for computing and manipulating surface geometry. A variety of applications, e.g., level set methods for tracking evolving interfaces, require accurate approximations of minimum distances to or closest points on implicitly defined surfaces. In this paper, we present an efficient method for calculating high-order approximations of closest points on implicit surfaces, applicable to both structured and unstructured meshes in any number of spatial dimensions. In combination with a high-order approximation of φ, the algorithm uses a rapidly converging Newton's method initialised with a guess of the closest point determined by an automatically generated point cloud approximating the surface. In general, the order of accuracy of the algorithm increases with the approximation order of φ. We demonstrate orders of accuracy up to six for smooth problems, while nonsmooth problems reliably reduce to their expected order of accuracy. Accompanying this paper is C++ code that can be used to implement the algorithms in a variety of settings.