Receiver functions provide an indispensable tool for producing discontinuity images of the crust and upper mantle from teleseismic earthquake arrivals. However, image quality can be significantly compromised by instabilities inherent in the deconvolution process and by irregularities in data coverage driven by source and/or receiver geometries. We present receiver function estimation and prestack migration techniques which both reduce receiver function deconvolution instability and produce regularized, multimode receiver function images. The receiver function estimation technique exploits frequency--wave number filtering methods using frequency--pseudowave number domain filtering of gathers to enhance arrivals with consistent moveout characteristics. Synthetic tests confirm robust recovery of features in the presence of high noise levels. To process regularized receiver functions into accurate images of the subsurface seismic scattering potential, we employ a regularized Kirchhoff migration that migrates both direct and reverberated P-to-S converted modes to their correct subsurface locations. The subsequent images can be stacked to mute mismigrated modes. Synthetic studies and model resolution tests demonstrate that this methodology is well suited to irregularly spaced stations and uneven data coverage and does not require the recording of a single earthquake at multiple stations.