Precision free-form surfaces are widely used in advanced optical and mechanical devices. In order to evaluate the form quality of a free-form surface, it is required to fit the measurement data with the design template and compare the relative deviation between them. A common approach is to minimize the sum of squared differences in the z direction. Its solution is not robust enough and may be biased due to outliers. This paper presents a fitting algorithm which employs the sum of orthogonal distances to evaluate the goodness of fit. The orthogonal projection points are updated simultaneously with the shape and motion parameters. Additionally, the l1 norm is adopted to improve the robustness of the solution. The Monte-Carlo simulation demonstrated that the bias in the fitted intrinsic characteristics of this method is much smaller than the traditional algebraic fitting, whereas the fitted motion parameters have no distinct difference.