We present an approximation algorithm for computing shortest paths in weighted three-dimensional domains. Given a polyhedral domain $$\mathcal D $$ D , consisting of $$n$$ n tetrahedra with positive weights, and a real number $$\varepsilon \in (0,1)$$ ε ∈ ( 0 , 1 ) , our algorithm constructs paths in $$\mathcal D $$ D from a fixed source vertex to all vertices of $$\mathcal D $$ D , the costs of which are at most $$1+\varepsilon $$ 1 + ε times the costs of (weighted) shortest paths, in $$O(\mathcal{C }(\mathcal D )\frac{n}{\varepsilon ^{2.5}}\log \frac{n}{\varepsilon }\log ^3\frac{1}{\varepsilon })$$ O ( C ( D ) n ε 2.5 log n ε log 3 1 ε ) time, where $$\mathcal{C }(\mathcal D )$$ C ( D ) is a geometric parameter related to the aspect ratios of tetrahedra. The efficiency of the proposed algorithm is based on an in-depth study of the local behavior of geodesic paths and additive Voronoi diagrams in weighted three-dimensional domains, which are of independent interest. The paper extends the results of Aleksandrov et al. (J ACM 52(1):25–53, 2005), to three dimensions.