High-resolution images of three-dimensional (3D) seismic structures are not only of scientific interest, but also of practical importance in predicting strong ground motion after large earthquakes. Given the source and station distributions, resolutions in current regional seismic tomography studies have been limited by two types of simplifying practices: the adoption of high-frequency approximations such as the ray theory and the use of one-dimensional (1D) reference (starting) models. We have developed a new approach to compute accurate finite-frequency 3D Fréchet (sensitivity) kernels of observed travel time and amplitude anomalies relative to 3D reference models. In our approach, we use a fourth-order staggered-grid finite-difference method to model the seismic-wave propagation in 3D media, and the reciprocity property of the Green's tensor to reduce the number of numerical simulations. This approach accounts for the perturbations in compressional- and shear-wave speeds in the same way, leading to a capability of inverting for the shear-wave speed directly from seismic data. The algorithm is readily parallelized to allow for realistic regional high-resolution 3D tomography inversions. We have implemented the algorithm for the Southern California Earthquake Center (scec) Community Velocity Model, sceccvm 3.0, a complex 3D model for Southern California including a number of sedimentary basins. By enabling the inversion of 3D structural perturbations to 3D reference models, our approach provides a practical means of iteratively solving the nonlinear regional tomography problems.