We investigate strategies to numerically integrate closed lines and surfaces that are implicitly defined by level sets (iso-contours) of continuously differentiable toroidally symmetric functions. The “grid-transform” approach transforms quantities given in non-surface-aligned coordinates onto a numerically constructed surface-aligned grid. Here, line and surface integrals, as well as so-called flux-surface averages, can be easily evaluated using high order integration formulas. We compare this method to ones that base on numerical representations of the delta-function.For the grid-transform method we observe high order convergence of line, surface and volume integration. Quantitatively, the errors for line and area integration are several orders of magnitude smaller than previously reported errors for delta-function methods. Furthermore, a delta-function method based on a Gaussian representation shows qualitatively wrong results of surface integrals near O- and X-points. Contrarily, the grid transform method suffers no deterioration near O-points. However, close to X-points we observe reduced first order convergence in volume integral and derivative tests due to the diverging volume element.Finally, we derive a toroidal integration based on toroidal summation and a smoothing kernel that assumes field-alignment of structures between the toroidal planes. The smoothing kernel can be interpreted as a partial flux-surface average. The resulting smoothed toroidal average eliminates unphysical poloidal oscillations that are otherwise present in the simple toroidal average.Our methods can be applied to toroidal and flux-surface averages in simulations of three-dimensional plasma dynamics on non-aligned grids. Further applications include closed line and surface integrals in level set methods. Efficient implementations can be found in the freely available Feltor library.