Computing device and method for processing seismic traces to produce an image of a subsurface area. The method includes receiving a series of seismic traces related to the subsurface area and recorded by one or more seismic receivers, wherein the one or more seismic sources are originally generated by a source; applying a phase encoding function to the series of seismic traces, at least a portion of said seismic traces comprise signals reflected by geological interfaces of the subsurface area; applying a 3 dimensional (3D) harmonic-source reverse time migration of the series of seismic traces encoded with the phase encoding function; computing a forward wavefield by solving a first wave equation; computing a backward wavefield by solving a second wave equation; and cross-correlating the forward wavefield with the backward wavefield to generate an image of the subsurface.