1. A magnetic resonance imaging method for the acquisition of a three-dimensional dataset, wherein spatial encoding by three mutually orthogonal magnetic field gradients is performed such that signals are readout under a read-gradient in one spatial direction k 1 with spatial encoding in the other two spatial directions k 2 , k 3 being performed by applying phase encoding gradients in those other two spatial directions prior to signal acquisition, wherein data acquisition is performed in a sequential manner such that, at each acquisition step, signals are acquired under the readout gradient but with different combinations of the two phase encoding gradients, the method comprising the steps of: a) fully sampling a k-space part 1 with a first sampling density, the k-space part 1 being a first subset of a k 2 -k 3 plane, which is symmetric around a k-space center; b) undersampling, with a second, uniform sampling density which is lower than the first sampling density, a k-space part 2, the k-space part 2 being a second subset of the k 2 -k 3 plane, which has higher spatial frequencies than those in k-space part 1 and which is symmetric around the k-space center; c) undersampling, with a third, uniform sampling density which is lower than the second sampling density, a k-space part 3, the k-space part 3 being a third subset in a lower half of the k 2 -k 3 plane, which has higher spatial frequencies than those in part 2; d) undersampling, with a fourth, uniform sampling density which is lower than the third sampling density, a k-space part 4, the k-space part 4 being a fourth subset in an upper half of the k 2 -k 3 plane, which has higher spatial frequencies than those in k-space part 2; e) acquiring no data within a k-space part 5, the k-space part 5 being a fifth subset in the upper half of the k 2 -k 3 plane, which has higher spatial frequencies than those in k-space part 4; and f) reconstructing images by iteratively minimizing a cost function with descent algorithms such as a nonlinear conjugate gradient method, the cost function being a weighted summation of multiple regularization terms, with a phase constraint term R pc being introduced into the cost function, wherein R pc =∥g∥ 1 =∥(x∘e −iP R −|x|)∘W∥ 1 , with x representing an intermediate solution in an iterative minimization; ∥g∥ 1 =Σ k |g k |, with g k being a k th element of matrix g; A∘B representing a Hadamard product of a matrix A and a matrix B; |x| a magnitude of x; P R an estimate of a phase of images to be reconstructed and W a weighting map.