Procedures for accurate ray tracing in complex three-dimensional media with interfaces are proposed. The ray tracing equations and the associated paraxial linear equations are solved either by a numerical solver or by an analytical perturbation approach. Interfaces are described with an explicit representation or an implicit representation using B-spline interpolation. For the implicit representation, we exploit two important properties of B-splines, the convex hull and subdivision properties, in order to determine the intersection of the ray with the interface.At the free surface where the recording system is located, a sampling strategy is proposed: limits of branches at caustics, shadow zones, and medium boundaries are detected for a fixed azimuth while the take-off angle is automatically adjusted in order to have a roughly homogeneous spacing between end points of the rays. The same strategy is also possible for a fixed take-off angle. The assumed continuity of the traveltime surface between two adjacent azimuths enables one to obtain the initial condition of a ray arriving at any station located on the portion of surface delimited by these two azimuths. This procedure allows for the classification of rays arriving at a given station as we show on different synthetic examples.