WO2014162300A1

WO2014162300A1 - Acceleration of low-rank mri data acquisition

Info

Publication number: WO2014162300A1
Application number: PCT/IB2014/060442
Authority: WO
Inventors: Stephen Smith; Thomas BLUMENSATH; Karla Miller; Mark CHIEW
Original assignee: Isis Innovation Ltd.
Priority date: 2013-04-05
Filing date: 2014-04-04
Publication date: 2014-10-09

Abstract

The present invention provides a method of recovery of an undersampled low-rank Magnetic Resonance Imaging (MRI) dataset. The method comprises the steps of: undersampling the dataset; and employing an Iterative Hard Thresholding and Matrix Shrinkage (IHT+MS) algorithm, on the data set, to directly reconstruct the low-rank data set.

Description

ACCELERATION OF LOW-RANK MRI DATA ACQUISITION

TECHNICAL FIELD

The present invention relates to a method of recovery of low-rank image time-series from under-sampled data that is embedded in a higher dimensional space. In the present context, "low-rank" means that the data is composed of a small number of linearly independent vectors in the presence of noise. Equivalently, it means that there are a small number of non-trivial singular values in a singular value decomposition of the data. The method is applied here to the acquisition of Magnetic Resonance Imaging (MRI) data and particularly to accelerating the data acquisition process. The method is particularly suited to but not limited to functional MRI data. That is, the acquisition and analysis of a time-series of brain images that enables detection of dynamic changes in brain activity.

The method of the present invention has been demonstrated in functional and cardiac MRI, hereinbelow. However, other potential applications include and are not limited to: dynamic angiography; elastography; velocity encoding; and perfusion time-series, dynamic contrast enhanced imaging, diffusion encoding, relaxometry and magnetic resonance spectroscopic imaging.

BACKGROUND

MRI data is commonly acquired in a reciprocal spatial domain (k-space), such that data can be subjected to an inverse Fourier transform to reconstruct an image of the signal variation across space, k-space is most commonly conceived of as a three-dimensional domain corresponding to the three spatial dimensions, discretized as a three-dimensional array of dimensions [N_x, N_y, N_z] . However, one can also consider each voxel independently and represent data as a one-dimensional vector of size N_xxN_yxN_z. Many techniques further collect a time series of MRI data, so that k-space can be extended in the temporal dimension to form a k-t space. This k-t space can be formulated as a 2D Nk x Nt data matrix, representing all spatial data along one dimension and all time-points along the other. In MRI, the time required to fully sample the k domain (to form a complete 3D image of a prescribed resolution and extent) limits the resolution in the temporal domain.

There are a number of existing methods that aim to accelerate MRI data. Classical "parallel imaging" methods reduce the amount of k-space data required to form an image by under- sampling k-space and filling in the missing data by combining data across multiple receive coil channels. A different category of methods, generally termed "compressed sensing", reconstruct images from under-sampled k-space data by imposing the constraint that a small number of coefficients from a chosen basis set can represent the data. However, these methods reconstruct each time point independently (i.e., these methods only consider k), and therefore do not take advantage of any temporal structure in the data. Some k-t methods have been proposed that aim to leverage temporal redundancy in the reconstruction process, however these require the prior specification of (spatial and temporal) basis functions or other structure, which may not be particularly appropriate for the data. One alternative is to use "matrix completion" methods, which do not presume a particular form for the temporal structure, but instead are driven by the structure of the data itself. There are existing k-t matrix completion techniques that have been applied to MRI and biomedical imaging to take advantage of the inter-dependence of spatial and temporal domains. However, these methods have been found to perform poorly on functional MRI data in particular, and in dynamic MRI methods more generally, compared to the methods proposed here.

The method of the present invention is to eliminate or at least to alleviate some of the problems encountered using the prior art methods, and to extend the application of said methods to multiple imaging domains.

SUMMARY

The present invention provides a method of accelerating the data acquisition process for dynamic MRI experiments. According to the present invention there is provided a method of recovery of an undersampled low-rank Magnetic Resonance Imaging (MRI) dataset, said method comprising the steps of (a) undersampling the dataset; and (b) employing an Iterative Hard Thresholding and Matrix Shrinkage (IHT+MS) algorithm on said data set to directly reconstruct said low-rank data set.

The method may comprise an Iterative Hard Thresholding and a Matrix Shrinkage algorithm, wherein the algorithm comprises the following equation:

where xⁿ represents the n iteration of the estimated k-t matrix, y is the sampled data, Φ is the sampling operator that selects measured matrix entries, μ is a step size parameter, and S is the thresholding and shrinkage operator.

Said data set may be a time-series of 2D or 3D images. Said dataset may be acquired using spatial frequency space (k-space) Fourier spatial encoding.

Said dataset may be undersampled pseudorandomly in k-t space by acquiring different k samples at different times in a manner that is incoherent with respect to aliaising due to the undersampling.

Said dataset may be undersampled pseudorandomly in k-v space, where v is a Fourier velocity encoding dimension.

Said dataset may be undersampled pseudorandomly in k-As space, where As is a shear wave propagation delay dimension.

Said dataset may be undersampled pseudorandomly in k-q or k-Δϋ space, where q is a reciprocal water displacement dimension, and AD is a diffusion time dimension.

Said dataset may be undersampled pseudorandomly in k-TE space, where TE is an echo-time dimension.

Approximation of the true missing values in said rank-reduced data matrix may be produced by a Matrix Completion process.

The method may further comprise the direct estimation of a data set post- dimensionality reduction using a post Principal Component Analysis.

Said dataset may be functional MRI (FMRI) data.

Said dataset may be used for independent component analysis of FMRI resting state networks.

Said dataset may be used for task-based FMRI analysis

Said dataset may be cardiac cine MRI data.

Said dataset may be angiography, blood flow or perfusion MRI data, with or without contrast enhancement.

Said dataset may be velocity encoded MRI data.

Said dataset may be elastography MRI data.

Said dataset may be diffusion MRI data.

Said dataset may be multi-echo relaxometry (T2 or T2*) MRI data.

BRIEF DESCRIPTION OF THE DRAWING

Many aspects of the disclosure can be better understood with reference to the following drawings. The components in the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating the principles of the present disclosure. Moreover, in the drawings, like reference numerals designate corresponding parts throughout the several views.

Figure 1 depicts a schematic of the method of the present invention;

Figure 2 illustrates zoomed cardiac cine MRI images of the heart at peak systole for 4/8/ 16x undersampling factors;

Figure 3 depicts time courses of a line through the cardiac images from figure 3 at 16x undersampling;

Figure 4 shows a diagram comparing standard, unaccelerated FMRI data acquisition compared to imaging accelerated with matrix completion for resting state FMRI analysis;

Figure 5 shows an example pseudorandom sampling pattern over k_z, with 8/32 k_z planes sampled every time point;

Figure 6 illustrates a comparison of a portion of 4^th fMRI principal components in the true data, data decimated to match undersampling, IHT+MS and FPCAr data. Results from first 3 principal components are visually indistinguishable;

Figures 7a and 7b are a comparison of Resting State Network (RSN) maps;

Figure 8 plots RSN map spatial correlation coefficients with truth data against smoothness estimates, which are a measure of spatial resolution loss;

Figure 9 shows representative diagrams of k-t undersamping patterns, where

Figure 9(a) illustrates pseudorandom undersampling,

Figure 9 (b) shows sheared grid undersampling, and

Figure 9 (c) shows standard low temporal resolution sampling;

Figure 10 shows example magnitude k-space time-series;

Figure 11 is a z-statistic image of a dual regression output for a single dataset for a regressor corresponding to an occipital (visual) resting state network;

Figure 12 is a z-statistic image of a dual regression output for a single dataset for a regressor corresponding to a bilateral parietal resting state network;

Figure 13 is a boxplot of correlation coefficients from the dual regression for each of the six subjects; and

Figure 14 is an z-statistic image of a dual regression output for a single subject for a regressor corresponding to a left somatosensory resting state network.

DETAILED DESCRIPTION OF THE DRAWINGS As mentioned above, prior art techniques use matrix completion for the recovery of under-sampled matrix data. In contrast to compressed sensing, where acceleration relies on enforcing constraints of sparseness (i.e., a low number of non-zero coefficients) in the image data or some transform domain thereof, the data for matrix completion are not required to be sparse - rather, the matrix needs to be of low rank for successful recovery to occur.

Further, these methods do not require a specified basis set (as is used in compressed sensing) but effectively estimate a basis that gives the best low-rank approximation of the data. These methods start from the observation that an a X b matrix of rank r has r(a + b— r) free parameters, compared to ab for a full rank matrix. Hence, it is possible under certain conditions to recover a low-rank matrix or approximation from a dataset that is not fully sampled. This concept has recently been applied to reconstruction of undersampled k-t MRI data (cardiac cine and dynamic contrast enhanced), where small amounts of coherent motion or signal enhancement in a static background produce data suitable for low-rank matrix recovery.

Figure 1 is a schematic of the method of the present invention and is a new algorithm designed for rank-reduced approximation of undersampled MRI data, based on iterative hard thresholding (IHT). The algorithm of the present invention is IHT coupled with Matrix Shrinkage (MS).

The IHT+MS algorithm can be summarised:

where xⁿ represents the n iteration of the estimated k-t matrix, y is the sampled data, Φ is the sampling operator that selects measured matrix entries, and μ is a step size parameter. S is the hard thresholding and matrix shrinkage operator.

To begin, an initial estimate x° of the full k-t matrix must be produced. In the present implementation, x° is produced from a Singular Value Decomposition (SVD) of the sampled measurement data y and keeping only the first r components. The next step of the algorithm involves applying the sampling operator Φ to x° to produce a sampled estimate Φ(χ°) that is non-zero only for matrix entries that correspond to measurements. Next, a residual difference is calculated between the sampled estimate and the sampled data, y— Φ(χ°). This residual difference is then scaled by a step size parameter μ. This scaled residual term is then added to the estimate to produce a new estimate that is in a sense closer to the true underlying k-t matrix. The final step in the algorithm is the hard thresholding and shrinkage operator S (described below). Application of S to xⁿ + μ(γ— Φ(χ⁰)) produces the next estimate in the iterative reconstruction process.

The shrinkage operator S uses the SVD to find the Γ+1^ώ largest singular value σ_Γ+1, above which all singular values are shrunk as follows;

0 otherwise where a is the i singular value of the matrix, r is the fixed rank cutoff, and τ is a soft- thresholding parameter which shrinks the remaining singular values by a fixed amount. In this example, the threshold is floating such that τ = ca_r+1 where c is a proportionality constant

0 < c < 1.

In the case where c=l, all singular values are shrunk to max (a - σ_Γ+ι, 0). Only the first r singular values survive this shrinkage and thresholding, producing a rank r data estimate after the decomposition is multiplied back with the reduced singular values. To be explicit, in a hard thresholding operation alone, the r largest singular values would be identified and kept, while discarding all the remaining singular values (and associated vectors). No shrinkage of singular values occurs in hard thresholding, meaning that the values σ_χ, σ₂, ... , _r are unmodified. In hard thresholding with matrix shrinkage, the same identification of the r largest singular values occurs, but the numerical value or amplitude of each singular value is shrunk by subtracting from each the numerical value of ca_r+1, the scaled Γ+1^ώ singular value. Figure

1 illustrates the application of this algorithm in the circumstances when c=l.

The algorithm is performed repeatedly until a stopping condition is met. In the present example, the stopping condition is that the iteration continues while the relative Frobenius norm, ||x^l— x^l-1 ||_F > a minimum update criteria or n < a preselected maximum number of iterations.

After the algorithm converges or halts, the final step of the matrix reconstruction is replacement of the originally sampled data back into the matrix estimate (i.e. the corresponding elements in the matrix element are replaced with the orginal, undersampled data).

The algorithm can be summarised formally as follows;

^■ Let the matrix y contain the sampled data, Φ be a sampling matrix with Is in sampled positions, and 0s in unsampled positions, and xⁿ be the n^th matrix estimate

^■ The^■ operator denotes element-wise multiplication ^■ Choose a rank r, step size 0 < μ < 1, and soft-thresholding parameter 0 < c < 1

^■ Initialise matrix estimate x°

^■ While ||xⁿ — x^{n— 1} II F > ^mm update criterion or n < max iterations

^■ z = Φ^■ (y - x')

^■ U∑V^* = ApproximateSVD(xⁿ + μ^■ z)

^■ For j = 1: n

^■ If j < r :∑(j, j) =∑(j, j) - c - σ_Γ+1

^■ Else :∑(j, j) = 0

^■ xⁿ⁺¹ = U∑V^*

^■ ApproximateSVD(x)

^■ QAQ^* = eigendecomposition(x^Tx)

^■ V = Q

^■ ∑ = A^{0 5}

^■ U = min||x - U∑V^* ||₂

u

The algorithm is thus an iterative process of sequentially updating the k-t matrix estimate by adding a scaled residual term computed from the sampled data, and hard thresholding and matrix shrinking that estimate by truncating and reducing the amplitude of the singular values of the estimate. The algorithm converges to a rank r approximation of the true underlying k-t data. Because real MRI data is only ever approximately low rank, the algorithm recovery performance can depend on the decay rate of singular values. While IHT performs well where singular values decay quickly, the IHT+MS algorithm of the present invention improves performance in slow decay regimes such as MRI data.

Accordingly, the method described herein enables faster data acquisition, particularly but not exclusively for fMRI. The k-t space is not fully sampled, which reduces the amount of time necessary for data collection. The "full" k-t space is then reconstructed by estimating its most important components. The reconstructed data is then used to generate one or more images, by a Fourier transform of the reconstructed data or otherwise. It will be apparent that the algorithm described above, and a method of generating an image using the algorithm and captured undersampled data, may be implemented on any suitable computer, whether integrally provided as part of an MRI apparatus or otherwise.

Advantageously, the method does not rely on knowledge of the coil sensitivities of the MRI device used to capture the undersampled data. However, when the coil sensitivities are known, the algorithm can be modified to take advantage of the extra information to reconstruct the undersampled matrix with lower reconstruction errors and at higher acceleration factors. Accordingly, the revised algorithm can be rewritten as;

where Ψ^" and Ψ represent the transforms for the reverse projection into the composite space and the forward projection onto multi-coil space respectively, each of which depend on the respective coil sensitivities. This approach takes advantage of the fact that the number of independent measurements increases in proportion to the number of coils, whereas the number of unknowns remains unchanged because a single 'true' (i.e. non-coil weighted) k-t matrix is estimated. The operator Ψ can be considered as a product of a block diagonal Fourier transform matrix F with a vertical stack of matrices, where each ip_t is a diagonal matrix corresponding the sensitivities of coil /. With / coils, Ψ is an Zm X m matrix, producing a Im X n multi-coil k-t matrix. Conversely, Ψ^" is the SNR optimal coil combination after inverse Fourier transform, represented Ψ- = ((ip_l) * R^~1 (ip_l)) ^~1 (ip_l) * R^~1F* , where * represents the conjugate transpose and R is a coil noise covariance matrix.

The method of the present invention was evaluated for its ability to recover both undersampled k-t FMRI and cardiac cine data. In the cardiac cine data, the IHT+MS algorithm was compared with the incremented rank power factorization (IRPF) and rank-constrained fixed point continuation approximation (FPCAr) algorithms at a reconstruction rank of 20. In the FMRI data, initial comparisons were made between the IHT+MS and FPCAr methods at rank 128, and further evaluation of the FMRI data RSN spatial maps focused on the IHT+MS algorithm only.

The method of the present invention was implemented in MATLAB, although any other suitable programming language with basic mathematical and linear algebra libraries may be used. Short axis cardiac cine data were acquired on a 3T MRI scanner (interpolated 352x512 spatial points, 100 phases over a single cardiac cycle). Retrospective sampling used 8 centre k_y lines and 4/8/ 16x random undersampling factors on the other 504 lines. Sliding window data were regularly undersampled to the same factor for comparison.

Cardiac data reconstruction fidelity was assessed using the relative Frobenius norm error between the true and reconstructed magnitude x-t data and visual inspection of reconstructed images (Figs. 2, 3). In the 4x undersampled cardiac data, the IRPF method produced the lowest error, compared to 1.9% for sliding window. At 8x and 16x, all images show visible artefacts, although qualitatively the IHT+MS images look least affected despite slightly higher error values than the sliding window method.

In Fig. 3, artefacts are least apparent in the IHT+MS images. The IHT+MS images instead appear smoothed or filtered, suggesting a graceful degradation of the spatio-temporal point spread function from the reconstruction algorithm. Approximate reconstruction times for the cardiac data were 150, 90 and 180 min respectively for the IHT+MS, FPCAr and IRPF methods (using a computer with 16x4 core 2.66 GHz processors, and 64 Gb RAM).

Figure 4 illustrates the principles underlying the use of matrix completion for acceleration of resting state FMRI data acquisition. In traditional acquisition strategies, k-t space is fully sampled (I) (use of parallel imaging or compressed sensing acceleration also results in a full k-t matrix). Subsequently, the k-t data are pre-processed using Principal Component Analysis (PC A) (II), and a rank reduction of the dataset is performed by truncating principal components (III). This corresponds to removing some components of the originally sampled data. Finally, this rank-reduced data is input to an Independent Component Analysis (ICA) step for identification of RSNs (IV). In the present invention, matrix completion is used to facilitate undersampling of the acquired data directly, reducing acquisition time (V). The IHT+MS algorithm is then used to directly reconstruct the rank reduced dataset (VI), which is then input to the ICA processing step (VII). In principle, if the data are low-rank, the rank- reduced datasets produced in (III) and (VI) are identical, where in the latter case acquisition time is saved by reducing the proportion of measurements taken, instead of taking a full set of measurements and discarding a proportion, as in the former.

Resting state FMRI data were acquired on a 3T MRI scanner (106x106x32 spatial points at 2x2x2 mm³, 512 time points at TR = 836 ms). This produced a k-t matrix size of 359552x512, where dimensions were limited by algorithm efficiency and do not constitute fundamental reconstruction limits. The present inventors simulated 3D-EPI with 4x undersampling of the k_z dimension, where the 4 centre planes were always sampled, and 4 of the remaining 28 planes were randomly sampled every time point (Fig. 5) for an effective undersampling factor of 4x. The original dataset was decimated by selecting every 4^th time point to provide a comparison to data from a non- accelerated acquisition with similar imaging parameters (designated [1/4 t] data), and a dataset with 25% of the spatial resolution in the z- direction was produced as an additional non-accelerated dataset with equivalent temporal resolution (designated [1/4 kz] data).

Figure 6 shows an example comparison of the 4^th (ordered by variance) spatial and temporal components in the true, [1/4 t], IHT+MS and FPCAr data ([1/4 kz] not shown for clarity) in the FMRI data. The IHT+MS algorithm produced the best estimates of the first 100 FMRI PCA spatial and temporal components. All subsequent results focus on evaluation of the IHT+MS algorithm.

To compare RSNs, a Multivariate Exploratory Linear Optimized Decomposition into Independent Components (MELODIC) probabilistic ICA analysis was run on the original data to produce 100 IC spatial maps. Dual regression was used to estimate equivalent maps for each IC from each reconstructed dataset, and z-scores were null-corrected using mixture modelling. Correlation coefficients were computed between original and reconstructed maps for 10 RSNs, and smoothness values were extracted for all spatial maps by estimating the size of image resolution elements.

Figures 7a and 7b show images from 2 representative RSNs, highlighting the excellent spatial agreement between the original spatial maps and those produced from undersampled data using the IHT+MS algorithm. The axial images of the somatomotor network (Fig. 7a, unthresholded) show nearly identical spatial correspondence between the true and IHT+MS data, whereas the [1/4 t] data show little representation of the network, and the [1/4 kz] data show over representation in that axial plane due to spatial blurring along the z-direction from the resolution decrease. The sagittal images of the fronto-parietal network (Fig. 7b, thresholded at Izl > 3), however, do suggest a slight loss of spatial fidelity in the undersampled z-direction for the IHT+MS dataset as well. Average smoothness values indicated that the IHT+MS data were 29% smoother than the original maps.

Figure 8 evaluates the correlation coefficients of the [1/4 t], [1/4 kz] and IHT+MS reconstructed datasets with the true data RSN maps against relative spatial smoothness estimates in 10 manually identified resting state networks. The smoothness parameter serves as a measure of spatial resolution, with higher smoothness indicating less spatial specificity. The [1/4 t] data show low correlations (mean + standard deviation 0.37+0.04) and relative smoothness values < 1 (0.34+0.03), indicating sparse and poorly represented RSNs. The [1/4 kz] data show higher correlations (0.57+0.03) but greatly increased smoothness (5.29+0.36), suggesting a pronounced loss of spatial resolution. The IHT+MS data have similarly high correlations (0.50+0.05), but with relative smoothness similar to the original data (1.29+0.14), indicating good fidelity RSN recovery with only a moderate loss of spatial resolution.

These results demonstrate the feasibility of matrix completion using the IHT+MS algorithm for accelerating the acquisition of FMRI RSN data. The IHT+MS algorithm shows excellent ability to recover a low rank approximation of undersampled MRI data. Qualitatively, both the spatial and temporal data reconstructed with IHT+MS show the least amount of artefact contamination, although quantitative error estimates can appear high because the algorithm does not preserve overall signal power. This may be irrelevant for applications such as resting state fMRI, in which temporal correlations are more important than signal amplitudes, particularly when variance normalization is performed. Convergence of the IHT+MS is reasonably fast, and efficient SVD approximation methods can be used, particularly when the k-t matrix is highly non-square. Rank-constrained reconstruction algorithms can be sensitive to the choice of rank, and prior information can be useful in selection of an optimal target rank.

The present method, therefore, allows the experimenter to increase the rate of data collection by a significant factor. Here a factor of 4x is demonstrated, meaning that 4 times as many volumetric images can be recorded at the same time, or the same number of images can be measured in ¼ of the time. However, any arbitrary acceleration factor can in principle be used. Absolute limits to the amount of acceleration depend on a number of factors, including the imaging context, the acquisition scheme, underlying rank of the imaging data and the fidelity of the desired reconstruction.

The method described herein was further evaluated in two ways: using retrospectively undersampled data (where the true k-t matrix is known), and in experiments using a prospectively-undersampled 3D EPI acquisition in healthy volunteers at 7T.

Two datasets, including six subjects each, were collected on healthy subjects with informed consent in accordance with local ethics. The first dataset used a simultaneous multi- slice EPI acquisition to enable simulations of retrospectively undersampled data with a known true k-t matrix at high temporal resolution. This data was acquired on a 3T (Siemens Verio) with an acceleration factor of 8 to achieve whole brain coverage with volume TR 836 ms, spatial matrix size of 106x106x64 at 2 mm isotropic resolution, and 1075 time points. The second dataset used prospective undersampling to provide a more realistic comparison of k-t FASTER with a duration-matched fully sampled acquisition. These data were collected on a Siemens 7T using a 32-channel head coil (Nova Medical) with a 3D segmented EPI acquisition modified to allow arbitrary undersampling of the kz dimension. The 3D EPI was acquired using an in-plane acceleration factor of 2 to achieve an appropriate TE for FMRI at 7T. The in-plane acceleration was reconstructed using GRAPPA, independently of and prior to matrix recovery. Acquisition used a spatial matrix size of 100x98x64 at 2mm isotropic resolution, with 72 time points in the fully sampled case and 304 time points for the accelerated data. Each subject was scanned twice with matched 5-minute scan durations in the fully sampled and accelerated cases respectively. Detailed parameters for imaging data are listed in Table 1. Retrospective Prospective undersampling undersampling

Field Strength 3T 7T

TR/TE (ms) 836/40 65/31

Spatial Matrix 106x106x64 100x98x64

Spatial Resolution 2mm Isotropic 2mm Isotropic

Measurements 1075 304

Volume TR (ms) 836 975

Sampling Factor 23.4% - 25% 23.4% - 25%

Rank Constant 128 72

Table 1: List of parameters used in the retrospective and prospectively sampled data.

Retrospective sampling and reconstruction of k-t matrices to was used to compare the method described herein (referred to as k-t FASTER) to 3 additional methods, with knowledge of the true k-t matrix. Table 2 lists the reconstruction algorithms that were compared, along with their associated sampling patterns

Table 2: Reconstruction methods compared, with associated sampling patterns

For k-t FASTER reconstruction, retrospective under- sampling was performed by selecting a pseudorandom subset of kz planes after Fourier transform of the (previously reconstructed) 3D volumes. Undersampling of 23.4% (i.e., ignoring 76.6% of the k-space data) was achieved by sampling the central 8 kz planes and 7 pseudorandomly-selected outer kz planes (in total 15 of 64 planes) at every time point (Fig. 9a). Reconstruction of the data was performed using the IHT+MS algorithm with rank constraint of 128. The shrinkage thresholding parameter c was set to 0.5, and was tuned to produce the best results in the retrospectively sampled data.

A known method, referred to as the "keyhole" method uses data sharing or interpolation to recover missing entries, and a keyhole-type reconstruction of the pseudorandomly undersampled data was performed by interpolating adjacent k-t entries to "fill in" the missing data. This method, denoted here as "k-t INTERP" performs a linear interpolation across time for every k-space location to generate values for the non-sampled matrix entries.

Another algorithm that has been proposed for low-rank k-t acceleration of MRI data: the linearly reconstructed PSF Model, here referred to as "k-t PSF". For this reconstruction, a sheared grid (lattice) pattern was used for sampling as shown in figure 9b. In the k-t PSF approach, the central strip of fully sampled k-space is used to estimate a temporal basis, which are then used to least squares fit a spatial basis using the grid sampled points. The central 8 kz planes were used as training data, combined with sheared grid sampling with an 8x undersampling factor. Rank constraints of 32 and 192 were used, which produced the best results after comparison with a range of candidate ranks (including the 128, the rank used for k-t FASTER reconstruction).

Finally, to examine the effect that temporal dimensionality has on the RSN analysis, the undersampled methods were compared with a standard long-TR fully sampled acquisition of the same total imaging duration, here referred to as "k-t SLOW". To generate this comparison, the truth data were temporally decimated by selecting only every 4th time point (figure 9b), such that a sampling density of 25% of the ground truth k-t matrix was achieved. The k-t SLOW data were constructed by selecting every 4^th volume from the ground truth data.

In the prospectively undersampled data, the k-t FASTER method using pseudorandom undersampling was compared to the slow, fully-sampled 3D EPI acquisition (k-t SLOW), in separate acquisitions from the same subjects. The k-t FASTER reconstruction differs slightly here from the retrospectively sampled case, using a lower rank constraint of 72 because of the reduced total number of time points in the acquisition. Furthermore, 32 k-t matrices were independently reconstructed corresponding to each coil channel. After reconstruction, these were combined using sum-of- squares to emphasize that the reconstruction was not driven by coil sensitivity information. In all other respects, except where explicitly noted, sampling and reconstruction were identical to the retrospectively sampled cases.

Data from the retrospective sampling comparisons were evaluated using dual regression with a set of 82 resting state spatial regressors (representing canonical network spatial maps) derived from high-dimensional group-level ICA of a large, high quality resting FMRI dataset from the Human Connectome Project. For the prospectively sampled data, only the top 64 and 32 RSN map regressors (sorted by variance in the original data) were used, due to the reduced temporal dimensionality in the data. Dual regression is a multiple linear regression framework used here to generate spatial z-statistic maps that reflect the degree of expression in each dataset of the original spatial regressors. The method first performs a spatial regression of each dataset against the canonical maps to extract the time-series corresponding to each map, and then performs a second temporal regression of each dataset against this set of times-series. The output is the specific spatial pattern corresponding to the desired network in that subject (which will, in general, vary from subject to subject). The dual regression z-statistic maps reflect the degree to which each spatial regressor, corresponding to a canonical RSN, is expressed with a unique time-course in the data.

To assess reconstruction fidelity, voxel-by- voxel correlation coefficients were computed for the output z-statistic maps from dual regression against both the corresponding regressor maps and the ground truth output maps to quantify the fidelity of each map in a single summary statistic. This fidelity metric captures the requirement for correct attribution of unique time-courses to corresponding spatial networks during matrix recovery in order to correctly produce a subject-specific version of a given RSN. The z-statistic maps output by the dual regression were null-corrected using a Gaussian and Gamma statistic mixture model

To illustrate the low-rank matrix reconstruction fidelity, Fig. 10 shows an example k- space magnitude time- series from a single representative dataset. That is, the timeseries associated with a single k-space position, or equivalently, a row of the k-t matrix. This point was chosen to show a time-series with median time-series reconstruction error (omitting fully- sampled k-space positions). In Fig. 10, a 100-timepoint zoomed view of the ground truth ('full rank') data is compared to a rank- 128 dimensionality reduced ground truth and k-t FASTER reconstructed time- series, with black circles indicating the sampled points (that is, the data used by k-t FASTER to reconstruct the time-series for this particular k-space location). In this view, many transient features of the k-space time-series can be seen between the sampled points with good overall correspondence. These features would be lost by slower traditional sampling or "keyhole" methods that would simply interpolate between the sampled points. Furthermore, although the transient features are recovered imperfectly, many of the reconstruction artefacts are shared between the k-t FASTER and rank- 128 data (arrows). This highlights the fact that using a rank 128 constraint in the k-t FASTER reproduces (and rejects) many of the same variance components as a rank- 128 decomposition of the ground truth that might be used in RSN analysis. Overall k-t matrix reconstruction fidelity for the various methods was evaluated using the relative Frobenius norm error, defined as:

\\x - y\\_F

err_F{x) = 100 - „ ,, %

\\y\\ F

where x is the estimated matrix, and y is the ground truth. Table 3 reports the err_F values for all 6 retrospectively- sampled datasets respectively, including the err_F obtained by performing a rank- 128 dimensionality reduction of the ground truth data. With this metric, the k-t FASTER method produces the lowest the err_F for all data sets, comparable to the amount of error present in a rank 128 truncated ground truth. The k-t INTERP and k-t PSF (rank 32) methods also produce reasonably low err_F values, due to the err_F metric being normalised by the total signal energy, which is concentrated in the fully-sampled central portion of k-space.

Table: Relative Frobenius norm errors for all 6 datasets across the three different reconstruction methods and the rank- 128 dimensionality reduction of the ground truth data.

Cumulatively, these results indicate that the k-t FASTER method generates significantly improved estimates of non- sampled k-t matrix locations than simple interpolation across time (k-t INTERP), or using explicit rank constrained model fitting (k-t PSF). Paired t- tests between the err_F values for k-t FASTER vs. the other methods, show significant differences (p < 0.05).

Figures 11 and 12 show example z-statistic maps from representative datasets for each method, alongside the ground truth maps, for an occipital visual RSN and bilateral parietal RSN respectively. In each of the examples, characteristics common to each reconstruction method can be observed. The k-t SLOW maps resemble noisier versions of the truth maps, lacking the temporal DOF necessary to capture all the detail in each map. Both k-t INTERP and k-t FASTER produced maps with generally good agreement with the ground-truth data, k- t INTERP generally has lower spatial fidelity than k-t FASTER data, although also appearing slightly less smoothed.

The prospectively-undersampled data were acquired in 5-minute runs, with the fully sampled k-t SLOW data having only 72 time points in that duration (volume TR = 4160 ms), while the k-t FASTER data produced 304 time points (volume TR = 975 ms). The 82 regressors used in the retrospectively sampled data were truncated to 64 for the dual regression analysis to accommodate the reduced temporal DOF in the k-t SLOW version of this data. As an additional test, the set of regressors was further reduced to 32, to further examine the effect of the DOF on these analyses. Figure 13 shows the correlation coefficients for the first 32 regressors (ordered by variance), by performing 2 separate dual regression analyses. The first dual regression used 32 regressors only, while the second dual regression 64 regressors, including the same 32 used in the first analysis. There is a striking increase in correlation in k- t SLOW when halving the number of regressors, suggesting insufficient DOF when using 64 regressors in a dataset with 72 time points. In comparison, k-t FASTER has consistently higher correlations than k-t SLOW, with essentially no change with increased number of regressors. These results demonstrate that, despite identical acquisition durations for the two approaches, the k-t FASTER data contains greater (and more functionally meaningful) temporal DOF than the k-t SLOW acquisition. This provides further evidence that k-t FASTER is not performing a trivial operation (e.g., interpolation) with the sampled points, but is actually identifying meaningful spatial-temporal structure in the data. The relatively low correlations values in figure 13 reflects the fact that these correlations were performed against canonical maps derived from a large set of group-averaged data, which are expected to show substantial mismatch with individual subject maps.

Figure 14 shows an example z-statistic map corresponding to a left somatosensory RSN in a representative dataset, compared to the corresponding spatial regressor map. These maps were chosen to be representative of the mean correlations for each case. The k-t FASTER maps clearly show a greater extent of high (and spatially well localised) z-statistics, most evident in the axial and sagittal views, and are a good reflection of the aggregate correlation differences between the k-t SLOW and k-t FASTER data expressed in the boxplots of Figure 13.

Accordingly, the method described herein is advantageous in that it exploits the spatio- temporal structure of FMRI data by using matrix completion methods to "fill in" the unsampled k-t matrix entries. The approximate low rank structure allows the method to generate the missing points in an intelligent way by preferentially describing high-variance (i.e. signal) components using an iterative SVD-based algorithm. Comparison with other methods shows that the technique described herein has lower errors, captures true temporal features that would be missed by methods that impose temporal smoothing, has higher RSN reproduction fidelity and preserves spatial features.

It should be emphasized that the above-described embodiments are merely examples of possible implementations. Many variations and modifications may be made to the above- described embodiments without departing from the principles of the present disclosure. All such modifications and variations are intended to be included herein within the scope of this disclosure and protected by the following claims.

Claims

1. A method of recovery of an undersampled low-rank Magnetic Resonance Imaging (MRI) dataset, said method comprising the steps of:

a) Undersampling the dataset; and

b) Employing an Iterative Hard Thresholding and Matrix Shrinkage (IHT+MS)

algorithm, on said data set to directly reconstruct said low-rank data set.

2. The method of claim 1, wherein the method comprises an Iterative Hard Thresholding and a Matrix Shrinkage algorithm, wherein the algorithm comprises the following equation:

3. The method of any preceding claim, wherein said data set is a time-series of 2D or 3D images.

4. The method of any preceding claim, wherein said dataset is acquired using spatial frequency space (k-space) Fourier spatial encoding.

5. The method of any preceding claim, wherein said dataset is undersampled

pseudorandomly in k-t space by acquiring different k samples at different times in a manner that is incoherent with respect to aliaising due to the undersampling.

6. The method of claims 1 to 4, wherein said dataset is undersampled pseudorandomly in k-v space, where v is a Fourier velocity encoding dimension.

7. The method of claims 1 to 4, wherein said dataset is undersampled pseudorandomly in k-As space, where As is a shear wave propagation delay dimension.

8. The method of claims 1 to 4, wherein said dataset is undersampled pseudorandomly in k-q or k-Δϋ space, where q is a reciprocal water displacement dimension, and AD is a diffusion time dimension.

9. The method of claims 1 to 4, wherein said dataset is undersampled pseudorandomly in k-TE space, where TE is an echo-time dimension.

10. The method of any preceding claim, wherein approximation of the true missing values in said rank-reduced data matrix is produced by a Matrix Completion process.

11. The method of claim 10, further comprising the direct estimation of a data set post- dimensionality reduction using a post Principal Component Analysis.

12. The method of claims 1 to 5 and 10 to 11, wherein said dataset is functional MRI (FMRI) data.

13. The method of claims 1 to 5 and 10 to 12, wherein said dataset is used for

independent component analysis of FMRI resting state networks.

14. The method of claims 1 to 5 and 10 to 12, wherein said dataset is used for task-based FMRI analysis

15. The method of claims 1 to 5 and 10, wherein said dataset is cardiac cine MRI data.

16. The method of claims 1 to 5, and 10, wherein said dataset is angiography, blood flow or perfusion MRI data, with or without contrast enhancement.

17. The method of claims 1 to 5 and 6, wherein said dataset is velocity encoded MRI data.

18. The method of claims 1 to 5 and 7, wherein said dataset is elastography MRI data.

19. The method of claims 1 to 5 and 8, wherein said dataset is diffusion MRI data.

20. The method of claims 1 to 5 and 9, wherein said dataset is multi-echo relaxometry (T2 or T2*) MRI data.

21. The method of claims 1 to 5 and 9, wherein said dataset is Magnetic Resonance

Spectroscopic imaging data