1. A method of correcting eddy currents-induced background phase (EC-BP) in magnetic resonance imaging data, comprising: acquiring phase-contrast data from a slice of interest (SOI) using a magnetic resonance imaging (MRI) scanner, with a user-specified table position and SOI position and orientation; the SOI position and orientation, with respect to the magnet coordinate system, being controlled by varying a center frequency of an applied radio-frequency (RF) excitation pulse in combination with a variation of a slice select gradient amplitude and direction; acquiring phase-contrast data from at least one additional slice, the at least one additional slice having the same position and orientation, with respect to the magnet coordinate system, as the SOI and being acquired using the slice select gradient amplitude and direction, the at least one additional slice imaging a physically different location of the subject, achieved by either changing the table position with respect to the position used for the SOI acquisition or by physically moving the subject within the scanner; performing image reconstruction from k-space to image domain for the SOI and the at least one additional slice; subtracting a phase of a reconstructed additional slice from a phase of a reconstructed SOI to extract phase difference maps for the SOI and the at least one additional slice; determining a correction map by: processing the phase difference maps from both the SOI and the at least one additional slice; and estimating EC-BP in accordance with a first relationship: x ^ mSAP = arg ⁢ ⁢ min x ⁢  Bx - z  Ω 2 , wherein xε k×1 represents the k coefficients of a 2D polynomial, wherein B and z are obtained by vertical concatenation of N temporally-averaged slices, such that B [A 1 T |A 2 T |. .. |A N T ] T and z [y 1 T |y 2 T |. .. |y N T ] T , and ∥Bx−z∥ Ω 2 (Bx−z)′Ω −1 (Bx−z), wherein each column of matrix Aε M×k defines the spatial distribution of one of the polynomial terms, wherein vector yε M×1 is the measured temporally-averaged phase map, wherein Ωε NM×NM is a covariance matrix constructed by computing inner products of mean-subtracted temporal profiles across all pixel pairs in z; and correcting the BP of the SOI using the correction map.