There are many definitions of the fractal dimension of an object, including box dimension. Bouligand-Minkowski dimension, and intersection dimension. Although they are all equivalent in the continuous domain, when discretized and applied to digitized data, they differ substantially. We show that the standard implementations of these definitions on self-afline curves with known fractal dimension (Weierstrass-Mandelbrot, Kiesswetter, fractional Brownian motion) yield results with errors ranging up to 5 or 10%. An analysis of the source of these errors led to a new algorithm in 1-D. called the variation method, which yielded accurate results. The variation method uses the notion of e-variation to measure the amplitude of the one-dimensional function in an e-neighborhood. The order of growth of the integral of the e-variation, as E tends toward zero, is directly related to the fractal dimension. In this paper, we extend the variation method to higher dimensions and show that, in the limit, it is equivalent to the classical box counting method. The result is an algorithm for reliably estimating the fractal dimension of surfaces or, more generally, graphs of functions of several variables. The algorithm is tested on surfaces with known fractal dimension and is applied to the study of rough surfaces.© (1987) COPYRIGHT SPIE--The International Society for Optical Engineering. Downloading of the abstract is permitted for personal use only.