A method for generating subsurface images, comprising: receiving, by a data processing apparatus, a set of seismic data (102a, 102b) associated with a subsurface region, wherein the set of seismic data corresponds to seismic signals received at one or more receiver devices; generating, by the data processing apparatus, source analytic wavefields and receiver analytic wavefields based on the set of seismic data; decomposing the source analytic wavefields and receiver analytic wavefields based on down-going and up-going components of the source analytic wavefields and receiver analytic wavefields; computing (140a, 140b), by the data processing apparatus, directions of propagations for the source analytic wavefields and receiver analytic wavefields based on the decomposed source analytic wavefields and receiver analytic wavefields, wherein the directions of propagations are computed using an optical flow method; computing, by the data processing apparatus, for a plurality of subsurface points, an azimuth angle and a reflection angle for a respective subsurface point based on the directions of propagations, wherein the reflection angle for the respective subsurface point is calculated according to cos θ = k → ⋅ k → r k → k → r, and the azimuth angle for the respective subsurface point is calculated according to cos φ = k → s × k → r ⋅ n → × k → s + k → r k → s × k → r n → × k → s + k → r, where θ is the reflection angle, φ is the azimuth angle, n is a unit vector indicating starting direction of the azimuth angle, k s is a wavefield vector at a source location, and k r is a wavefield vector at a receiver location; generating, for each of the plurality of subsurface points, a weighting function for a respective subsurface point using the azimuth angle and the reflection angle of the respective subsurface point, wherein the weighting function represents illumination energy at the respective subsurface point with respect to the azimuth angle and the reflection angle of the respective subsurface point; and generating a subsurface image using the weighting functions of the plurality of subsurface points. receiving, by a data processing apparatus, a set of seismic data (102a, 102b) associated with a subsurface region, wherein the set of seismic data corresponds to seismic signals received at one or more receiver devices; generating, by the data processing apparatus, source analytic wavefields and receiver analytic wavefields based on the set of seismic data; decomposing the source analytic wavefields and receiver analytic wavefields based on down-going and up-going components of the source analytic wavefields and receiver analytic wavefields; computing (140a, 140b), by the data processing apparatus, directions of propagations for the source analytic wavefields and receiver analytic wavefields based on the decomposed source analytic wavefields and receiver analytic wavefields, wherein the directions of propagations are computed using an optical flow method; computing, by the data processing apparatus, for a plurality of subsurface points, an azimuth angle and a reflection angle for a respective subsurface point based on the directions of propagations, wherein the reflection angle for the respective subsurface point is calculated according to cos θ = k → ⋅ k → r k → k → r, and the azimuth angle for the respective subsurface point is calculated according to cos φ = k → s × k → r ⋅ n → × k → s + k → r k → s × k → r n → × k → s + k → r, where θ is the reflection angle, φ is the azimuth angle, n is a unit vector indicating starting direction of the azimuth angle, k s is a wavefield vector at a source location, and k r is a wavefield vector at a receiver location; generating, for each of the plurality of subsurface points, a weighting function for a respective subsurface point using the azimuth angle and the reflection angle of the respective subsurface point, wherein the weighting function represents illumination energy at the respective subsurface point with respect to the azimuth angle and the reflection angle of the respective subsurface point; and generating a subsurface image using the weighting functions of the plurality of subsurface points.