Let [Formula: see text] be a [Formula: see text]-smooth compact submanifold of dimension [Formula: see text]. Assume that the volume of [Formula: see text] is at most [Formula: see text] and the reach (i.e. the normal injectivity radius) of [Formula: see text] is greater than [Formula: see text]. Moreover, let [Formula: see text] be a probability measure on [Formula: see text] whose density on [Formula: see text] is a strictly positive Lipschitz-smooth function. Let [Formula: see text], [Formula: see text] be [Formula: see text] independent random samples from distribution [Formula: see text]. Also, let [Formula: see text], [Formula: see text] be independent random samples from a Gaussian random variable in [Formula: see text] having covariance [Formula: see text], where [Formula: see text] is less than a certain specified function of [Formula: see text] and [Formula: see text]. We assume that we are given the data points [Formula: see text] [Formula: see text], modeling random points of [Formula: see text] with measurement noise. We develop an algorithm which produces from these data, with high probability, a [Formula: see text] dimensional submanifold [Formula: see text] whose Hausdorff distance to [Formula: see text] is less than [Formula: see text] for [Formula: see text] and whose reach is greater than [Formula: see text] with universal constants [Formula: see text]. The number [Formula: see text] of random samples required depends almost linearly on [Formula: see text], polynomially on [Formula: see text] and exponentially on [Formula: see text].