We present a novel approach for the direct computation of integral surfaces in general vector fields. As opposed to previous work, which we analyze in detail, our approach is based on a separation of integral surface computation into two stages: surface approximation and generation of a graphical representation. This allows us to overcome several limitations of previous techniques. We first describe an algorithm for surface integration that approximates a series of timelines using iterative refinement and computes a skeleton of the integral surface. In a second step, we generate a well-conditioned triangulation. The presented approach allows a highly accurate treatment of very large time-varying vector fields in an efficient, streaming fashion. We examine the properties of the presented methods on several example datasets and perform a numerical study of its correctness and accuracy. Finally, we examine some visualization aspects of integral surfaces.