We present the first polynomial time approximation algorithm for computing shortest paths in weighted three-dimensional domains. Given a polyhedral domain $\D$, consisting of $n$ tetrahedra with positive weights, and a real number $\eps\in(0,1)$, our algorithm constructs paths in $\D$ from a fixed source vertex to all vertices of $\D$, whose costs are at most $1+\eps$ times the costs of (weighted) shortest paths, in $O(\C(\D)\frac{n}{\eps^{2.5}}\log\frac{n}{\eps}\log^3\frac{1}{\eps})$ time, where $\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, Maheshwari and Sack [JACM 2005] to three dimensions.