Since seismic imaging creates an image of the subsurface structure based on information received from the measured wavefield, it is essential to fully utilize the reflected waves. Full Wavefield Modeling (FWMod) was developed with recursive and iterative up-and-down wavefield propagation, using one-way wave propagation, to model both primary and multiple reflections. Using FWMod as the modeling engine, Full Wavefield Migration (FWM) has been introduced to directly image data including internal multiples, where internal multiple crosstalk is suppressed automatically via an inversion-based data-fitting process. This avoids the need for applying internal multiple removal, which is often challenging. Conventional one-way wave propagators calculated in the wavenumber domain, like the phase shift (PS) operator, have limitations when applied to strongly inhomogeneous media. Even when computing a new operator at each lateral grid point, they still suffer difficulties because the medium is assumed to be locally homogeneous. In the past, matrix eigendecomposition has been proposed as a way to create accurate, local velocity-based one-way propagation operators. In this article, an accurate propagator based on eigendecomposition is incorporated into FWMod and FWM. In the numerical examples, four models with strong lateral velocity variations were used to test the propagator. With a comparison of the conventional FWM based on the PS operator with input data including FWMod and a finite-difference (FD) approach, the numerical examples demonstrated that the proposed method has the potential to significantly enhance image reflectivity, suppress internal multiples, and maintain convergence speed during the least-squares inversion.