A magnetic resonance imaging method for the acquisition of a three-dimensional dataset, where spatial encoding by three mutually orthogonal magnetic field gradients is performed such that signal is readout under a read-gradient in one spatial direction k 1 and spatial encoding in the other two spatial directions k 2, k 3 is performed by applying phase encoding gradients in the other two spatial directions prior to signal acquisition and data acquisition is performed in a sequential manner such that at each acquisition step signal is acquired under the said readout gradient but with a different combination of the two phase encoding gradients, characterized in that a first central subset region of the k 2 -k 3 plane referred to as Part 1, which is symmetric around the k-space center, is fully sampled, a second subset region of the k 2 -k 3 plane referred to as Part 2, whose outer edge has higher absolute value spatial frequencies along either of k 2, k 3 or both directions than that of Part 1, the second subset region being symmetric around the k -space center, surrounding the first subset region and being undersampled with uniform density, a third subset region covering the remainder of the lower half of the k 2 -k 3 plane, i.e. of the half with k 2 <0, not being covered by Parts 1 and 2 and being referred to as Part 3, is undersampled with uniform density, a fourth subset region in the upper half of the k 2 -k 3 plane, i.e.in the half with k 2 >=0, referred to as Part 4, whose outer edge in said upper half has higher absolute value spatial frequencies along either of k 2, k 3 or both directions than that of Part 2, said fourth subset region surrounding the second subset region in said upper half and being undersampled with uniform density, a fifth subset region covering the remainder of said upper half of the k 2 -k 3 plane and being referred to as Part 5, which has higher absolute spatial frequencies along either of k 2, k 3 or both directions than Part 4, is completely unacquired; in that the sampling density in Part 2 is lower than that in Part 1, the sampling density in Part 3 is lower than that in Part 2, the sampling density in Part 4 is lower than that in Part 3, no data are acquired in Part 5; and in that images are reconstructed by iteratively minimizing a cost function, i.e. by iteratively solving arg min x f c ( x) with a descent algorithm like in the nonlinear conjugate gradient method, wherein the cost function is a weighted summation of the term ∥U(Fx-v)∥ 2 and multiple regularization terms where y represents the acquired k-space data, F represents the Fourier transform operator, U represents the k-space sampling mask with binary entries, wherein one of said multiple regularization terms is a phase constraint term referred to as R pc of the form R pc = ∥g∥ 1 =∥(x∘e -iP R -|x|)∘W∥, where x represents an intermediate solution that minimizes the cost function f c ( x) in the iterative minimization; ∥g∥ 1 = ∑ k |g k |, g k is the k th element of matrix g and is calculated with the k th element of matrix x,and P R; k represents the integer index of every element of a matrix; A ∘ B represents the Hadamard product of matrix A and matrix B; |x| represents the magnitude of x; P R represents an estimate of the phase of the images to be reconstructed, and P R is taken from an image reconstructed by using only the symmetrically acquired data in Part 1 and Part 2 in said k 2 -k 3 plane; and wherein W represents a weighting map.