Summary: Interpolation of scattered data has many applications in different areas. Recently, this problem has gained a lot of interest for CAD applications, in combination with the process of reverse engineering, i.e., the construction of CAD models for existing objects. Until now, no method for scattered data interpolation with a bivariate function has produced surface formats that can be directly integrated into a CAD system. Additionally many of the existing interpolation schemes exhibit undesirable curvature distribution of the reconstructed surface. In this paper we present a method for scattered data interpolation producing tensor--product B-splines with high quality curvature distribution. This method first determines the knot vectors in a way that guarantees the existence of an interpolating B-spline. In a second step the degrees of freedom not specified by the interpolation constraints are automatically set using a data dependent optimization technique. Examples demonstrate the quality of the resulting interpolants w.r.t. curvature distribution and approximation of known surfaces.