This paper algorithmically solves the problem of finding an approximating isotopic simplicial surface for a compact regular level set \([f=0]\) of a \(C^2\) function \(f\) defined in space. This topic is highly relevant for computing and visualization, unfortunately the authors mention that they have not yet implemented their algorithm. The basic setting is a suitable simplicial decomposition \(T\) of the domain of interest and \(f\)'s associated piecewise-linear interpolant \(\widehat f\). The main theorem is as follows: Supposing \(f\), \(T\) are generic in a well-specified way, consider a subcomplex \(W\leq T\) such that \(f\neq 0\) in \(\partial W\), neither \(f\) nor \(\widehat f\) has critical points in \(W\), \(W\) collapses to \([\widehat f=0]\), and such that \(f,\widehat f\) have the same index in each component of \(W^c\). Then \([f=0]\) and \([\widehat f=0]\) are isotopic in \(W\), with a bound on their Hausdorff distance. It is further shown how to construct \(W\) from given \(f\) and \(T\) (`watershed' construction). The algorithm itself proceeds iteratively, and is initialized by using a hierarchical octree decomposition of the domain into boxes such that for each box either \(f\neq 0\) or \(\nabla f\neq 0\). For that purpose it is assumed that \(f\) is accessible via interval arithmetic and questions regarding nonvanishing of \(f,\nabla f\) can be answered. Proofs rely on the smooth and stratified Morse theories as well as on an appropriate discrete concept of critical points and Morse theory.