If by "2D surface z = f(x, y)" you mean a curve in the plane, that is, the set of solutions in $x$ and $y$ to this equation for a particular value of $z$, which is how I interpret the question in any event:
Here's a quick description of a procedure for finding the curves $f(x, y) = z = \text{constant}$, first, assume one knows a "sufficiently nice", i.e., differentiable, functional form for $f(x, y)$. Then taking $\nabla f(x, y) = ((\partial f/\partial x)(x, y), (\partial f /\partial y)(x, y))$ one obtains the gradient vector field of $f(x, y)$. Then define the vector field $\vec T(x, y)$ via
$\vec T(x, y) = J \nabla f(x, y), \tag{1}$
where
$J = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}. \tag{2}$
It is easy to see that
$\nabla f \cdot \vec T = \nabla f \cdot J \nabla f = 0, \tag{3}$
since $J$ is a skew-symmetric matrix; in fact, $J$ is an orthogonal matrix which represents a rotation by $\pi/2$ in the counter-clockwise direction. In any event, (3) shows that $\vec T$ is in fact tangent to the level sets of $f(x, y)$. Thus computing the integral curves of $\vec T$ is tantamount to computing the level sets of $f(x, y)$; if we pick some point $(x_0, y_0)$ then $f(x, y) = f(x_0, y_0)$ for any $(x, y)$ on the trajectry through $(x_0, y_0)$.
Any number of methods may be used to calculate the trajectories of $\vec T$; one concern will be to keep the timestep $h$ sufficiently small to obtain adequate resolution; similar concerns apply to the spacing of the initial points $(x_0, y_0)$;
a tighter spacing of initial values and a smaller timestep will generally lead to a higher resolution at the cost, of course, of greater run times to complete the process.
Of course, if a functional form for $f(x, y)$ isn't readily available, e.g., you are working with a set of measured data instead of a given expression for $f$, then you would have to construct $\nabla f(x, y)$ from information extracted from the data sets; this is where algorithms such as marching squares or "triangulation and slicing" make their appearance. Either technique is a structured way of searching the given data set with an eye to (approximately) constructing the level sets of $f(x, y)$, but they can also be used to extract the normal field $\nabla f$ from the raw input. If this is done in a sufficiently high resolution manner, then we still obtain a vector function $\nabla f(x, y)$ which may be used as above to find and integrate the vector field $\vec T$. But in the process of so doing, one will have effectively obtained the level sets of $f$, making the integration step I suggested above redundant.
But, there is another way to do things: if the data for $f(x, y)$ lies on a grid, one can in principle directly compute $\nabla f(x, y)$ via differencing, setting for example
$\dfrac{\partial f}{\partial x}(i, j) = \dfrac{f(i + 1, j) - f(i - 1, j)}{2h}, \tag{4}$
where $h$ is the step-size in the $x$ direction, with an analogous formula for $\dfrac{\partial f}{\partial y}(i, j)$; there are also other differencing schemes available. Then one can use the vector array $\nabla f(i, j)$ to perform the procedures previously outlined. If the data set is triangulated, one can interpolate to a grid and so obtain an array of the form $\nabla f(i, j)$, and proceed from there.
It is a challenge to estimate the complexity or run-time of a code addressing this problem constructed from choices between various possible subsidiary algorithms such as the ones mentioned here. It might be more efficient to work directly from marching squares or sliced triangulations. But either of these approaches require a degree of interpolation which itself can be time-consuming. On the other hand, calculating straight from an array of the form $f(i, j)$ to obtain $\nabla f(i, j)$ may in fact be less work. I don't have a ready answer at this point; perhaps someone with a firmer command of these topics could chime in . . .
I worked at an avant garde computer graphics firm in a previous incarnation, and ideas such as those presented here often came up in the coffee room. Ultimately, the gridifiers/triangulators/polygonizers won the day in the market place, largely due to the (relative) ease with which such algorithms may be mapped to, and hence implemented in, hardware. But such notions may yet again surface and find there place in the world of commerce. Meanwhile, it's a pleasure to reminisce and find inspiration in novel (to yours truly at least) approaches to such problems.
Hope this helps. Cheers,
and as always,
Fiat Lux!!!
