Fractal objects derive from many interface phenomena, as they arise in, for example, materials science, chemistry and geology. Hence the problem of estimating fractal dimension becomes of both theoretical and practical importance. Existing algorithms implement the standard definitions of fractal dimension directly, but, as we show, often give unreliable results when applied to digitized and quantized data. We present a new algorithm for estimating the fractal dimension of surfaces--the variation method--that is more reliable and robust than the standard ones. It is based on a new definition of fractal dimension particularly suited for graphs of functions. The variation method is validated with both fractional brownian surfaces and Takagi surfaces, two classes of mathematical objects with known fractal dimension, and is shown to give more accurate results than the classical algorithms. Finally, our new algorithm is applied to data from sand-blasted metal surfaces.