High-quality and artifact-free image reconstruction is a perennial task in a plethora of imaging technologies, such as dMRI; a high-fidelity visualization technology with widespread applications in cardiac-cine, dynamic-contrast enhanced and neuroimaging [1], [2]. Due to various physiological constraints, MR scanners are usually slow in data acquisition and unable to keep up with organ motions or fluid flow [1]. Moreover, there is a constant need to speed up the data collection process to make it inexpensive and to cause less patient discomfort. The unpleasant end result is the availability of a limited number of measurements/data, which leads in turn to distorted and aliasing image reconstructions [1], [2].

Naturally, a lot of research focuses on reconstruction algorithms that yield high-fidelity image series given the usually inadequate number of measurements from MR scanners. Recent years have seen a paradigm shift to artificial-intelligence (AI) approaches in image reconstruction, such as convolutional neural networks (CNNs) [3], [4] and generative adversarial networks (GANs) [5], [6]. However, these techniques rely on the availability of large number of training data and a time consuming training phase preceding the reconstruction task.

In contrast to AI approaches, this paper focuses on methods that do not use any training data, but explore and exploit instead latent geometries of the observed data. Among such methods, a mainstream approach is compressed sensing [7]–[11]. For instance, [10] estimates first a temporal basis from image series via singular-value decomposition, followed by the estimation of a spatial basis through a sparsity inducing convex optimization task.

Manifold-learning techniques, identifying intrinsic manifold geometries beneath raw data for data recovery, have found applications in face recognition [12], [13], as well as dMRI [14]. Use of local SVD/PCA in data modeling has gathered also interest in applications not necessarily related to MRI [15]. A graph Laplacian matrix, capturing relations within data clouds, is used as a regularizer in a convex recovery problem in [14]. Motivated by the success of kernel methods [16] in extending linear models to non-linear counterparts, the MRI community has also seen the introduction of kernel methods in MRI data recovery tasks [17]–[19]. Notwithstanding, such kernel-based approaches require explicit and often costly (kernel) preimaging steps, and leave the complex-valued nature of the k-space data underused.

The present study builds on [20] to establish its generalization towards manifold geometries in reproducing kernel Hilbert spaces (RKHSs). To promote a computationally efficient scheme, a small set of landmark points is chosen to describe concicely the data cloud in the input space, and then mapped non-linearly to the RKHS feature space. Driven by the continuity of the mapping (kernel) function, it is hypothesized that the mapped feature points lie on or close to a smooth unknown and low-dimensional manifold. Departing from mainstream manifold-learning approaches, where the graph Laplacian operator is used as a regularizer in an inverse problem, e.g., [14], this study relies heavily on the concept of tangent spaces of smooth manifolds to establish a local and data-adaptive method of revealing affine relations among data points. These affine relations, which take place in the RKHS feature space, amount to non-linear operations in the original input space, yielding thus better data

Fig. 1. (a) The (k,t)-space. “Navigator (pilot) data” comprise the gray-colored area of the (k,t)-space (). (b) k-space with 1-D Cartesian sampling pattern. (b) In the feature space , mapped landmark points are affinely combined to describe possible affine combinations of are represented by the gray-coloured area.

approximations according to kernel-method arguments [16], [21]. Costly preimaging operations to solve for vectors in the original input space are avoided via a bi-linear modeling approach. The framework employs also a dimensionality reduction module to impose a low-rank structure to the data. Furthermore, terms capitalizing on the periodicity that exists often in dMRI are also introduced to achieve high-fidelity image reconstructions. Finally, recently developed convex and non-convex minimization techniques are employed to solve the resultant recovery tasks. It is worth stressing here the points that distinguish the proposed framework from prior art: i) No training data are used; ii) in contrast to mainstream manifold-learning techniques, no graph Laplacian matrix is used to penalize the optimization task; iii) bi-linear terms are used to approximate data; iv) no costly (kernel) preimaging step is necessary to map feature points back to the input space; and v) complex-valued kernel functions are defined to account for the k-space data. Preliminary numerical tests on synthetically generated high-dimensional dMRI phantom data and comparisons against state-of-the-art methodologies underline the rich potential of the method for high-quality data recovery.

An MR image (C is a set of all complex numbers) is observed via at discrete k-space (or frequency domain) locations, spanning an area of where and stand for number of phase and frequency encoding lines, respectively. Measurements Y can be viewed as the two-dimensional (discrete) Fourier transform of X: Y = F(X) [1]. Without any loss of generality, it is assumed that the “low-frequency” part of Y is located around the center of the area. In dMRI, a temporal dimension is added to yield where represents the number of frames collected over time (Fig. 1a), resulting in the augmented (k,t)-space. The (k,t)-space can be seen as the collection of k-space measurements .

In practice, due to physical and physiological phenomena, k-space measurements are missing, causing severe undersampling of the k-space, which in turn results in distorted, artifact-induced reconstructions of [1]. The present framework’s objective is to reconstruct images from limited k-space measurements but free of aliasing effects and distortions. These limited k-space measurements are usually made along predefined trajectories for efficient acquisition of MR data [2]. One such trajectory is the one-dimesnional (1D) Cartesian sampling, exhibited in Fig. 1b, where measurements are recorded only along the “white” lines and the “black” space is assumed to be zero-filled. To this end, the proposed framework considers the availability of these heavily sampled central k-space data, coined navigator data which are usually a small number of phase encoded lines (“gray-coloured region” in Fig. 1a), and the highly (pseudo-randomly) undersampled high frequency region as in [17], [18], [20], [22].

To facilitate processing, (k,t)-space data are vectorised. More specifically, operation stacks one column of below the other to form a vector , where . To avoid notation clutter, F still denotes the two-dimensional (discrete) Fourier transform even when applied to vectorized versions of image frames: . All vectorized k-space frames are gathered in the matrix so that the vectorized image-domain data are . The navigator data of the jth k-space frame (cf. Fig. 1a), are gathered into a vector . All navigator data comprise the matrix .

The navigator data carries useful spatio-temporal information in the (k,t) space. However, tends to be usually large and to promote parsimonious data representations, it is desirable to extract a subset (coined landmark points; ) which concisely describes the data cloud . Stacking as columns in a matrix, is defined. To reveal potential non-linear spatio-temporal dependencies in , embedding into a high dimensional (potentially infinite) complex RKHS (feature space) is pursued. The feature map

is defined via the reproducing kernel [16]. The celebrated kernel trick/ , facilitates processing [16], where is the inner product associated with , and H, which denotes the Hermitian transposition of vectors/matrices, is used here to simplify notations.

The dimensions of the kernel matrix K is contingent on the number of landmark points chosen and therefore, K can still be high dimensional depending on the physical characteristics of the acquired data. To impose a low-rank structure in As. 2, and to meet the restrictions imposed by limited computational resources, dimensionality reduction on K is desirable. To this end, the methodology of [26] is applied:

By As. 2 and 3, the (k,t)-space measurements Y can be modeled as . Upon defining, , , and by the virtue of the linearity of . This establishes the following bi-linear relation between Y and the unknowns . This bi-linear relation holds true also in the image domain: .

As discussed in Sec. II, only a few measurements are available in the (k,t)-space. To account for missing entries, a sampling operator is introduced to emulate sampling trajectories. The operation S(Y) retrospectively under-samples Y retaining the sampled k-space locations and filling the other locations with zeros. Given the positive real-valued parameters , the following inverse problem is formulated:

where denotes the ith column of the identity matrix . In addition to the data fit term T1, (2) contains terms T2 and T4 where the auxiliary variable Z is introduced to capture the periodic process, (e.g., beating heart motion) by imposing sparsity on performs 1D Fourier Transform of the Nfr time profile of every single pixel. The term T3 imposes the sparsity constraint and C2 accounts for the affine constraint on B as in As. 2. Term C1 prevents the unbounded solution of D due to the scaling ambiguity in .

To solve (2), the successive-convex-approximation framework of [28] is employed and is presented in Alg. 1. Convergence to a stationary solution of (2) is guaranteed [28]. Step 8 of Alg. 1 calls for solving convex minimization sub-tasks. Given , for at every iteration of the algorithm, the following estimates are required:

(3a) ˆarg min+ + + s.to (3b)

Tasks in (3) can be viewed as affinely constrained composite convex minimization tasks, hence allowing the use of [27]. The algorithmic and mathematical details of (3) are reserved for an upcoming journal publication.

The proposed framework is validated on the magnetic-resonance extended cardiac-torso (MRXCAT) cine phantom [30]. Further investigation via testing on other datasets is beyond the scope of this paper due to lack of space. The extended cardiac torso (XCAT) framework was used to generate a breath-hold cardiac cine data of spatial size with number of frames. The generated data has a spatial resolution of mmfor a field-of-view (FOV) of mm. The phantom data consists of 15 cardiac cycles and 24 cardiac phases. The data is retrospectively undersampled under 1D Cartesian sampling strategy, as in Fig. 1b, to emulate an undersampling rate [defined by # of acquired voxels)] of 8. Normalised root mean square error (NRMSE), defined as , is used to measure the quality of reconstruction, with X being the fully sampled, high-fidelity, original MR image and the estimate of X obtained from various reconstruction schemes.

The performance of KBiLMDM is tested against various state-of-the-art techniques, namely, the partially separable sparsity aware model (PS-Sparse) [10], smoothness regularization on manifolds (SToRM) [14], kernel low rank (KLR) [18], and the bi-linear modeling of data manifolds (BiLMDM) [20]. KBiLMDM achieves the least NRMSE value of 0.03891 as compared to PS-Sparse’s 0.0556, SToRM’s 0.0503, KLR’s 0.0744, BiLMDM’s 0.0456. KBiLMDM shows less deviation in NRMSE than KLR, SToRM and PS-Sparse while consistently scoring the least NRMSE as compared to BiLMDM, KLR and PS-Sparse across all frames. A further qualitative analysis of the reconstructions, as in Fig. 2 provides evidence that KBiLMDM outperforms the other aforementioned state-of-the-art techniques. KLR and SToRM exhibit severe aliasing effects as marked in Fig. 2, and

Fig. 2. Spatial Results for frame 49 of MRXCAT cardiac cine phantom retrospectively undersampled at acceleration rate of 8. Left to Right: Gold standard (top row); Undersampled Image (bottom row), PS-Sparse (0.0556), SToRM (0.0503), KLR (0.0744), BiLMDM (0.0456), KBiLMDM (0.0389). The numerical values indicate the NRMSE. Top row: Spatial frames generated from aforementioned reconstructed schemes; bottom row: Corressponding error maps.

Fig. 3. Frame-wise NRMSE of PS-Sparse KiBLMDM . The previous numerical values indicate mean and standard deviation of the NRMSE across 360 frames.

the error maps of PS-Sparse provides evidence that reconstructions along the edges are not upto the mark, while there are discrepancies in the contrast as compared to the ground truth. On the other hand, KBiLMDM produces sharper, aliasing-free, distortion-free dMRI time series. A comparison of the error maps between BiLMDM and KBiLMDM shows improvements in the reconstructions due to the introduction of kernels. Details on choice of kernels and comparison with schemes like FRIST-MRI [11] are left for presentation and upcoming journal publications.

A kernel-based framework for reconstructing data on manifolds was proposed and tailored to fit the dMRI-data recovery problem. The proposed methodology exploited simple tangent-space geometries of manifolds and followed classical kernel-approximation arguments to form a bi-linear inverse estimation problem. Departing from state-of-the-art approaches, i) no training data were used; ii) no graph Laplacian matrix was employed to penalize the inverse problem; iii) no costly (kernel) preimaging step was necessary to map feature points back to the input space; and iv) complex-valued kernel functions were defined to account for k-space data. Preliminary numerical tests and qualitative investigation of reconstructed images on retrospectively undersampled synthetic MR images showed the rich potential of the proposed framework against several state-of-the-art techniques. With regards to the road ahead, the incorporation of training data into the framework is currently under investigation. On-going research includes also modeling extensions, further numerical tests on a wide range of synthetic and real dMRI data, which will be reported during the conference.

[1] Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Processing Perspective. IEEE Press, 2000.

[2] ——, “An efficient method for dynamic magnetic resonance imaging,” IEEE Trans. Medical Imag., vol. 13, no. 4, pp. 677–686, 1994.

[3] D. Liang, J. Cheng, Z. Ke, and L. Ying, “Deep magnetic resonance image reconstruction: Inverse problems meet neural networks,” IEEE Signal Processing Magazine, vol. 37, no. 1, pp. 141–151, 2020.

[4] J. Schlemper, J. Caballero, J. V. Hajnal, A. N. Price, and D. Rueckert, “A deep cascade of convolutional neural networks for dynamic MR image reconstruction,” IEEE Trans. Medical Imag., vol. 37, no. 2, pp. 491–503, 2018.

[5] M. Mardani, E. Gong, J. Y. Cheng, S. S. Vasanawala, G. Zaharchuk, L. Xing, and J. M. Pauly, “Deep generative adversarial neural networks for compressive sensing MRI,” IEEE Trans. Medical Imag., vol. 38, no. 1, pp. 167–179, 2018.

[6] K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D. K. Sodickson, T. Pock, and F. Knoll, “Learning a variational network for reconstruction of accelerated MRI data,” Magnetic Resonance in Medicine, vol. 79, no. 6, pp. 3055–3071, 2018.

[7] M. Lustig, J. M. Santos, D. L. Donoho, and J. M. Pauly, “k-t SPARSE: High frame rate dynamic MRI exploiting spatio-temporal sparsity,” in Proc. ISMRM, vol. 2420, 2006.

[8] R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI,” Magnetic Resonance in Medicine, vol. 64, no. 3, pp. 767–776, 2010.

[9] H. Jung, J. C. Ye, and E. Y. Kim, “Improved k-t BLAST and k-t SENSE using FOCUSS,” Physics in Medicine and Biology, vol. 52, no. 11, pp. 3201–3226, 2007.

[10] B. Zhao, J. P. Haldar, A. G. Christodoulou, and Z.-P. Liang, “Image reconstruction from highly undersampled (k,t)-space data with joint partial separability and sparsity constraints,” IEEE Trans. Medical Imag., vol. 31, no. 9, pp. 1809–1820, 2012.

[11] B. Wen, S. Ravishankar, and Y. Bresler, “FRIST—Flipping and rotation invariant sparsifying transform learning and applications,” Inverse Problems, vol. 33, no. 7, p. 074007, 2017.

[12] J. Jiang, R. Hu, Z. Han, Z. Wang, T. Lu, and J. Chen, “Locality-constraint iterative neighbor embedding for face hallucination,” in 2013 IEEE International Conference on Multimedia and Expo (ICME). IEEE, 2013, pp. 1–6.

[13] P. Zhang, X. You, W. Ou, C. P. Chen, and Y.-m. Cheung, “Sparse discriminative multi-manifold embedding for one-sample face identification,” Pattern Recognition, vol. 52, pp. 249–259, 2016.

[14] S. Poddar and M. Jacob, “Dynamic MRI using smoothness regularization on manifolds (SToRM),” IEEE Trans. Medical Imag., vol. 35, no. 4, pp. 1106–1115, 2016.

[15] A. V. Little, M. Maggioni, and L. Rosasco, “Multiscale geometric methods for data sets I: Multiscale SVD, noise and curvature,” Applied and Computational Harmonic Analysis, vol. 43, no. 3, pp. 504–567, 2017.

[16] B. Scholkopf and A. J. Smola, Learning With Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT press, 2001.

[17] S. Poddar and M. Jacob, “Recovery of noisy points on bandlimited surfaces: Kernel methods re-explained,” in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018, pp. 4024–4028.

[18] U. Nakarmi, Y. Wang, J. Lyu, D. Liang, and L. Ying, “A kernel-based low-rank (KLR) model for low-dimensional manifold recovery in highly accelerated dynamic MRI,” IEEE Trans. Medical Imag., vol. 36, no. 11, pp. 2297–2307, 2017.

[19] O. Arif, H. Afzal, H. Abbas, M. F. Amjad, J. Wan, and R. Nawaz, “Accelerated dynamic MRI using kernel-based low rank constraint,” Journal of Medical Systems, vol. 43, no. 8, p. 271, 2019.

[20] G. N. Shetty, K. Slavakis, A. Bose, U. Nakarmi, G. Scutari, and L. Ying, “Bi-linear modeling of data manifolds for dynamic-MRI recovery,” IEEE Trans. Medical Imag., 2019. [Online]. Available: https://doi.org/10.1109/TMI.2019.2934125

[21] M. P. Wand and M. C. Jones, Kernel Smoothing. Chapman and Hall/CRC, 1994.

[22] U. Nakarmi, K. Slavakis, and L. Ying, “MLS: Joint manifold-learning and sparsity-aware framework for highly accelerated dynamic magnetic resonance imaging,” in Proc. ISBI, 2018, pp. 1213–1216.

[23] R. Boloix-Tortosa, F. J. Payán-Somet, E. Arias-de Reyna, and J. J. Murillo-Fuentes, “Complex kernels for proper complex-valued signals: A review,” in 2015 23rd European Signal Processing Conference (EUSIPCO). EURASIP, 2015, pp. 2371–2375.

[24] L. K. Saul and S. T. Roweis, “Think globally, fit locally: Unsupervised learning of low dimensional manifolds,” J. Machine Learning Research, vol. 4, pp. 119–155, 2003.

[25] L. W. Tu, An Introduction to Manifolds. New York: Springer, 2008.

[26] K. Slavakis, G. B. Giannakis, and G. Leus, “Robust sparse embedding and reconstruction via dictionary learning,” in Proc. CISS, Baltimore: USA, Mar. 2013.

[27] K. Slavakis and I. Yamada, “Fejér-monotone hybrid steepest descent method for affinely constrained and composite convex minimization tasks,” Optimization, vol. 67, no. 11, pp. 1963–2001, 2018.

[28] F. Facchinei, G. Scutari, and S. Sagratella, “Parallel selective algorithms for nonconvex big data optimization,” IEEE Trans. Signal Process., vol. 63, no. 7, pp. 1874–1889, 2015.

[29] V. De Silva and J. B. Tenenbaum, “Sparse multidimensional scaling using landmark points,” Stanford University, Tech. Rep., 2004.

[30] L. Wissmann, C. Santelli, W. P. Segars, and S. Kozerke, “MRXCAT: Realistic numerical phantoms for cardiovascular magnetic resonance,” J. Cardiovascular Magnetic Resonance, vol. 16, no. 1, p. 63, 2014.