Coordinate metrology often involves fitting a geometric surface to coordinate data. The increasingly rigorous approaches to uncertainty evaluation being developed across metrology have been reflected in the need to evaluate the uncertainties associated with coordinate data and calculate how these uncertainties are propagated through to uncertainties associated with the parameters describing the fitted surface. Standard least squares algorithms such as orthogonal distance regression (ODR) for finding the best-fit surface implicitly assume that the uncertainties associated with the coordinates are uncorrelated and axis-isotropic, i.e. the uncertainties associated with the x-, y- and z-coordinates are equal, but very few coordinate measuring systems have such uncertainty characteristics. Given a model of the measurement system, it is possible to propagate uncertainties associated with systematic and random effects influencing the measurements to determine an uncertainty matrix associated with the measured coordinates. This uncertainty matrix will generally be full, reflecting correlation amongst all the measured coordinates, but will have a factorization structure that gives a compact representation of the sources of uncertainty. In this paper, we derive surface fitting algorithms that take account of the coordinate uncertainties in order to provide maximum likelihood estimates of the surface parameters. We show that ODR algorithms can be adapted to take into account the type of uncertainty matrices that arise in practice, exploiting the factorization structure to enable very efficient implementations of these algorithms. As a result, a significant class of surface fitting problems can be solved efficiently using standard nonlinear least squares algorithms. We also discuss how the algorithms can be modified to account for probe radius effects. The calculation of the uncertainties associated with the fitted surface parameters is also derived in such a way that systematic effects associated with the measuring system, such as scale errors, are automatically accounted for.