The seismic body wave tomography method has been improved and extended to adapt to a general velocity structure with a number of complexly shaped seismic velocity discontinuities (SVDs) and with three-dimensional variations in the velocities in the modeling space. An efficient three dimensional ray tracing algorithm which iteratively uses the pseudobending technique and Snell's law is developed. The large and sparse system of observation equations is solved by using the LSQR algorithm. This method is applied to 18,679 arrival times from 470 shallow and intermediate-depth earthquakes in order to study P and S wave tomographic images beneath northeastern Japan. In addition to the first P and S wave arrivals, clear later arrivals of SP waves converted at the Moho and PS SP waves converted at the upper boundary of the subducted Pacific plate (UBPP) are also used in the inversion. The UBPP, Conrad and Moho are take as three SVDs, and their depth distributions obtained by previous studies are used. High-resolution P and S wave tomographic images down to a depth of 200 km have been determined. Large velocity variations amounting to 6% for P wave and 10% for S wave are revealed in the crust and upper mantle. In the crust low-velocity (low-V) zones exist beneath active volcanoes. In the upper mantle the low-V zones dip toward the west from the volcanic front. A high-velocity (high-V) zone corresponding to the subducted Pacific plate is clearly delineated. Most earthquakes in the lower plane of the double-planed deep seismic zone are found to occur in relatively high-V areas. The obtained tomographic images are also found to explain other seismological observations well. ¿American Geophysical Union 1992