It is often of interest to infer lower-dimensional structure underlying complex data. As a flexible class of non-linear structures, it is common to focus on Riemannian manifolds. Most existing manifold learning algorithms replace the original data with lower-dimensional coordinates without providing an estimate of the manifold or using the manifold to denoise the original data. This article proposes a new methodology to address these problems, allowing interpolation of the estimated manifold between the fitted data points. The proposed approach is motivated by the novel theoretical properties of local covariance matrices constructed from samples near a manifold. Our results enable us to turn a global manifold reconstruction problem into a local regression problem, allowing for the application of Gaussian processes for probabilistic manifold reconstruction. In addition to the theory justifying our methodology, we provide simulated and real data examples to illustrate the performance.