In the case of Cartesian k-space sampling, Nyquist’s theorem states that the distance Δk _i between adjacent sampling points has to fulfil the condition Δk _i ≤ 1/FOV _i , where FOV is the field of view and i = x, y, z reflects the standard three-dimensional Euclidean basis.

For a variety of non-Cartesian k-space trajectories it appears more convenient to evaluate an upper limit by using Pythagoras’ theorem in 2D or 3D k-space, i.e.

Δ k max = max { d ∈ K } Δ k x , d 2 + Δ k y , d 2 + Δ k z , d 2 , where d is the set of all distances in k-space K between points that are nearest neighbours and do not belong to the same k-space interleave or read-out (in case of multi-shot acquisitions). Since Δk _max ≥Δk _i,d ∀d ∈ K and i = x, y, z, Nyquist’s theorem is indeed fulfilled if min(1/FOV _i ) ≥Δk _max. While e.g. in the case of radial sampling, this evaluation can be restricted to the sampling point of each read-out that is farthest away from the centre of k-space, it appears insufficient for the case of any quasi-random sampling point distribution.

In the scope of this publication it is useful to extend the estimation of a Δk _max even further, since the later introduced distribution of sampling points will not follow any regular or symmetric pattern, as shown for a single sample point C and its six surrounding nearest neighbours P _i with i = 1, …, 6 in Figure 1C for simplification in two dimensions.

In order to calculate a local estimate of the encoded FOV at C, the direct distances d C , P i between C and each P _i are calculated and the mean value of all six distances defines a radius Δr around C. This radius can then define a generalised FOV for C with FOV _C = 1/Δr _C but in the same manner also for every other point in k-space. In the three dimensional case, Δr defines the radius of a sphere, from which an equal FOV along each direction is derived. For this calculation, a sufficient number of nearest neighbours is required, to ensure that Δr is derived from the expression Δ r = Δ x i ̄ 2 + Δ y i ̄ 2 + Δ z i ̄ 2 with Δx _i ≠ 0, Δy _i ≠ 0 and Δz _i ≠ 0, i.e. the radius can define a spherical (3D) volume. The given definition also allows the assignment of an isotropic FOV for various specific regions in k-space as well as for the entire trajectory by averaging Δr for all considered k-space points.

The previous publication [17] has used Jacobi elliptic functions for the construction of a single k-space interleave as the solution to an optimisation problem given by P. Erdös in 2000 [19]. Here, we have generalised the approach by making use of Jacobi theta functions [20], giving rise to the possibility of constructing a multitude of inherently different k-space waveforms. For example, Figure 1A shows two different k-space waveforms of arbitrary length, for which the k _x − and k _y − components are once based on θ ₁ and once based on θ ₄ (k _z = θ ₂). A projection of both waveforms into the k _x − k _y − plane is shown in Figure 1B. All corresponding gradient waveforms are shown in Figure 2.

The waveform presented in the following is meant to prove the imaging concept for a lower limit of read-out durations and is generally based on the concepts and parameters introduced in [17]. The general waveform is generated on the surface of a unit sphere according to the definition: ζ : R 0 + → R 3 , s↦ζ(s) with η ∈ (0, 1) and ζ x s , η = θ 1 s , η 2 ⋅ cos s ⋅ η 2 (1) ζ y s , η = θ 1 s , η 2 ⋅ sin s ⋅ η 2 (2) ζ z s , η = θ 2 s , η 2 (3)

While η is a parameter to adapt the waveform to hardware limitations such as maximum gradient amplitudes and available slew-rates. The combination with sine and cosine terms in the first two components allows for a modifiable change in direction per unit length while the symmetry along the z-direction remains unchanged (Figures 1A,B). The length of the waveform is determined by s and therefore the restriction s ∈ [0, s _max] should be applied, with s _max being sufficiently large according to the desired resolution (extension of k-space).

For this publication, one k-space interleave was constructed with η = 0.5 and by setting the target resolution to 0.85 mm (isotropic). For an extended oversampling of k-space centre α = 1.3 was chosen to facilitate a subsequent CS reconstruction. Each presented trajectory was generated with a maximum gradient strength of 21 mT/m and a slew-rate of 120 T/m/s. Each trajectory was furthermore explicitly optimised by minimising its discrepancy [17], by iterative optimisation of s, η and the angle of rotation around the symmetry axis according to [17] and using an Euclidian arc-length parameterization [21]. An adequate choice of the centre-oversampling parameter α requires information about the energy distribution of the underlying k-space and can therefore not be generally included to the optimisation problem. Nevertheless, defining an interval by experience e.g. α ∈ [1.2, 1.7] allows for another degree of freedom during the optimisation process and is therefore recommended.

The resulting read-out duration for each interleave was 3.52 ms, in order to reach the boundary of the k-space sphere for the defined maximum gradient amplitude and slew-rate limits. Ten interleaves of a final trajectory with 20,000 interleaves are depicted in Figure 1D for illustrative purposes. For this trajectory, we define an undersampling factor of R: = 1 due to an equal scan time as the corresponding Cartesian acquisition at the same resolution and FOV.

Based on the previously mentioned parameters, four sampling point spread functions were calculated according to undersampling factors of R = 1, 8, 12, 16 with respect to the trajectory with the largest number of interleaves, i.e. 20,000. Each trajectory for each undersampling factor was independently generated and optimised with respect to low-coherent sampling properties. All PSF_Ss were obtained by a separate calculation of a Voronoi tessellation [7] in Euclidean space ( R 3 , Euclidean distance) for every generated trajectory. Each PSF_S was obtained by gridding unit k-space data onto a Cartesian grid in combination with a 3D Voronoi tessellation [17, 22] to estimate weights for the necessary density compensation of non-uniformly acquired k-space data. Based on all normalised PSF_Ss, the centre-peak FWHM was determined in order to evaluate relative image sharpness with increasing undersampling factors. Furthermore, the peak/side-lobe ratio PSR = max i ≠ j PSF S i , j PSF S i , i (4)

Of each PSF_S is an appropriate measure for emerging coherences [3].

In order to evaluate the aliasing behaviour, as well as imaging performance, in-vivo head images were acquired using a 3.0 T wholebody MRI system (Achieva 3.0 T, Philips, Best, Netherlands) with an 8-element SENSE Neuro coil (Philips, Best, Netherlands).

Image reconstruction for all 3D ζ-based Spirals was performed as follows: after data acquisition, raw data were exported and processed with MATLAB (MathWorks, Natick, MA, United States). Images were obtained using gridding [23] with an oversampling factor of 1.25, a Kaisser Bessel kernel for interpolation and again in combination with a 3D Voronoi tessellation. Gradient system delays were estimated [24] and used to correct the trajectory before gridding. Further eddy current effects were compensated using a mono-exponential model [25] with a time constant of τ = 39 μs. No post-processing was applied to any of the presented image. Compressed Sensing reconstructions were accomplished using an in-house written CS/SENSE reconstruction based on nonlinear conjugate gradient methods as suggested in [3, 26]. Undersampling was again realized by calculating separate and optimised trajectories with 12,000, 8,000, 4,000 and 1,250 interleaves, leading to undersampling factors R = 1.66, 2.5, 5 and 16 with respect to the initial trajectory with 20,000 interleaves. These undersampling factors were chosen to correspond to scan durations of 90, 60, 30 and 10s. All relevant scan parameters are listed in Table 1.

To provide an experimental assessment of the aliasing properties, we evaluate the noise characteristics of the acquired in-vivo datasets for the 3D ζ-based Spiral, acquired with 12,000, 8,000, 4,000 and 1,250 interleaves each, according to the previous section. A reference dataset with the vendor’s 3D radial Kooshball trajectory was acquired, employing the same spatial resolution and choosing a FOV encompassing the entire head. Image reconstruction for the Kooshball trajectory followed the description given in the previous section for the 3D ζ-based Spiral trajectory, except for the weighting calculation. Kooshball weights were calculated analytically, based on the symmetry of the sampling scheme (spherical shells). The radial dataset was retrospectively undersampled by a random selection of 1/R spokes.

Three regions of interest (ROI) as highlighted in Figure 3D were selected in the background of the reconstructed images, ideally containing no source of any MR signal. Accordingly, the pixel intensities are exclusively governed by artefacts and noise whose characteristics can be analysed by consideration of the associated power spectrum. The calculation of the power spectrum followed [27], generalised to the 3D case. To facilitate the analysis of their characteristics (line shapes), all power spectra were normalised and are presented in arbitrary units. For each acquisition, the power spectra were averaged over all coil elements and evaluated with respect to the three geometrical axis to capture possible similarities (coherences).

According to the definition of the generalised FOV, the presented trajectory with 20,000 interleaves has the following properties: The Nyquist condition Δr _C ≤ 1/FOV = (1/220) mm is fulfilled within a sphere of radius r _N = 0.18 ⋅ k _max. The largest FOV that is stored within a sphere of radius r = 0.01 ⋅ k _max corresponds to 38-times the FOV dimension (220 mm). According to the definition, the Nyquist condition is not fulfilled for points outside a sphere of radius r _N.

Since the number of interleaves for the trajectory with 20,000 interleaves was selected such that the resulting scan duration is approximately equal to a standard Cartesian acquisition (174 s), without considering the actually encoded FOV, it is reasonable to calculate the introduced generalised FOV for this and for all other trajecories, to achieve undersampling factors that truly rely on Nyquist’s condition and not on a relative number of interleaves. Therefore, the trajectory with 20,000 interleaves corresponds to an undersampling factor of R 1 ′ ≈ 3.29 . This value was determined by calculating the mean generalised FOV of all sampling points of each trajectory, stating that Nyquist’s theorem is violated 3.29-times according to the definitions. Accordingly, the equally generated trajectory with 2,500 interleaves led to R 2 ′ ≈ 6.43 , to R 3 ′ ≈ 9.33 for 1,665 interleaves and to R 4 ′ ≈ 12.98 for 1,250 interleaves. The intention of these numbers is merely to classify the presented trajectories than to enforce comparisons to other sampling schemes due to drastic differences in the distribution of points in k-space. For simplification, all images and results that correspond to undersampling factors R 1 ′ , … , R 4 ′ of the 3D ζ-based Spirals are denoted by M ₃, M ₆, M ₉, M ₁₃ according to the mean undersampling factors.

Figure 4 shows all associated sampling point spread functions for these four cases of undersampling in the xy-plane with z = 0.

With increasing undersampling, energies in the PSF_Ss emerge that do not seem to follow any ordered or symmetric pattern. Consequently, all PSF_Ss appear to be governed by a low-coherent distribution of energies with an expected aliasing behaviour that is (in its appearance) vastly similar to an introduction of white noise.

Figures 4E,F shows a logarithmic plot of the centre region of the four PSF_Ss, shown in Figures 4A–D and a cross section of the entire PSF_S in 3b). As the undersampling factor increases, an overall increase in energy can be appreciated in which the side-lobe behaviour shows rarely signs of emerging coherences. Furthermore, image sharpness is preserved for all undersampling factors. The FWHM of the PSF_S centre-peak results to ≈ 2.544   pixel in width (mean) with a maximum deviation of 0.59% between the broadest peak (M ₁₃: 2.549 px) and the narrowest peak (M ₉: 2.534 px) of all undersampling PSF_Ss. All values were obtained in non-logarithmic representation. This finding of retained sharpness is furthermore supported by the in-vivo images, presented in the following section.

Figure 5 shows axial and sagittal slices, acquired with the presented 3D ζ-based Spiral trajectories for different numbers of interleaves (20,000, 8,000 and 1,250). All data was directly gridded and no additional image or data processing was applied before and after gridding. Furthermore, no sensitivity maps were used, in order not to alter the emerging imaging artefacts. Figure 5 also shows a Compressed Sensing reconstruction based on the dataset with 1,250 interleaves using a total variation regularisation [7].

The generalised FOV for 20,000 interleaves is ≈ 67  mm (M ₃, isotropic), ≈ 41  mm for 8,000 interleaves and ≈ 17  mm for 1,250 interleaves (M ₁₃). The reconstruction of a larger FOV, in this case of 220 mm (isotropic) for all datasets results in additional low-coherent aliasing artefacts, which can clearly be appreciated. Despite uncorrected coil sensitivity profiles, all images appear non-degraded by coherent aliasing artefacts, especially visible in the background of all images. As expected from the PSF_S analysis, image sharpness is preserved, and based on the optical impression, equal for all investigated undersampling factors.

The noise analysis for the three cubic regions (ROI) is shown in Figure 3. The normalised power spectra indicate slightly overpronounced DC components for all acquisitions. For the 3D ζ-based Spiral trajectories (20,000, 8,000, 4,000 and 1,250 interleaves), all further spatial frequency components are about equally represented, leading to a widely flat power spectrum in accordance to the behaviour of (bandwidth limited) white noise. As expected, an increasing undersampling factor increases the overall power, but the characteristics (line-shapes) remain widely unchanged. Concerning the Kooshball trajectory, with the analysis being shown in the Supplemental Data (Supplementary Figure S1), all directions show a similar behaviour in terms of overpronounced DC components and a following decline of the power spectrum. But, all directions in each evaluated ROI show clear modulations (wave pattern) in the frequency analysis with again increasing powers as the undersampling factor increases. The latter also introduces obvious changes to the frequency modulations, indicating varying (coherent) aliasing artefacts. While the noise-like characteristics of aliasing artefacts remain unchanged for varying undersampling factors in the case of 3D ζ-based Spiral trajectories, undersampling of the radial Kooshball trajectory introduces streak artefacts.

Calculating the peak/side-lobe ratios for the four PSF _s s shown in Figure 4 yields ratios in non-logarithmic representation of 4.7 ⋅ 10^–3 (M ₃), 5.1 ⋅ 10^–3 (M ₆), 5.2 ⋅ 10^–3 (M ₉), 5.8 ⋅ 10^–3 (M ₁₃), indicating that emerging coherences remain at the same level with increasing power densities in the PSF _s s towards the outer regions as R increases. Compared to Kooshball sampling with the PSF _s s shown in the Supplemental Data (Supplementary Figure S2), the peak/side-lobe ratio for R = 3 is 5.3 ⋅ 10^–3 while it increases to 11.2 ⋅ 10^–3 for R = 6, to 13.8 ⋅ 10^–3 for R = 9 and to 20.6 ⋅ 10^–3 for R = 13, which is in accordance to the visually enhanced coherences in the emerging aliasing artefacts for increasing undersampling factors.