This paper describes a fast algorithm for scattered data interpolation and approximation. Multilevel B-splines are introduced to compute a C-2-continuous surface through a set of irregularly spaced points. The algorithm makes use of a coarse-to-fine hierarchy of control lattices to generate a sequence of bicubic B-spline functions whose sum approaches the desired interpolation function. Large performance gains are realized by using B-spline refinement to reduce the sum of these functions into one equivalent B-spline function. Experimental results demonstrate that high-fidelity reconstruction is possible from a selected set of sparse and irregular samples.