Reverse time pre-stack depth migration, which uses the two-way acoustic wave equation, is not a new concept. Conventionally the method has been very computationally intensive and, therefore, has been considered impractical for production 3D depth imaging projects. Here we describe an efficient algorithm that can be used on large scale 3D seismic data. To make it practical and efficient we employ explicit 2nd order in time and high order in space domain finite differences. We also use threading and domain decomposition methods to split the image cube amongst multiple CPU?s, when necessary. Only a few finite differencing layers are communicated between related CPU?s, by message passing during domain decomposition. High order spatial finite differences handle numerical dispersion and allow larger time steps than those possible with the, more conventional, pseudo-spectral method. We will show an overview of the method along with 2D and 3D synthetic and real data examples.