Single projection driven real-time multi-contrast (SPIDERM) MR imaging using pre-learned spatial subspace and linear transformation

As in Feng et al (2020), a 4D image $I\left({\bf{x}},t\right)$ can be modeled as low-rank using partially separable functions. Using a matrix expression, an image of the following form

$$\begin{eqnarray}{\bf{A}}=\left[\begin{array}{ccc}a\left({{\bf{x}}}_{1},{t}_{1}\right) & \cdots & a\left({{\bf{x}}}_{1},{t}_{{N}_{t}}\right)\\ \vdots & \ddots & \vdots \\ a\left({{\bf{x}}}_{J},{t}_{1}\right) & \cdots & a\left({{\bf{x}}}_{J},{t}_{{N}_{t}}\right)\end{array}\right]\end{eqnarray} \tag{ 1 }$$

can be decomposed as

$$\begin{eqnarray}{\bf{A}}={{\bf{U}}}_{{\bf{x}}}{{\boldsymbol{\Phi }}}_{{\rm{rt}}}=\left[\begin{array}{ccc}{{\boldsymbol{u}}}_{1} & \cdots & {{\boldsymbol{u}}}_{L}\end{array}\right]\left[\begin{array}{ccc}{{\boldsymbol{\phi }}}_{1} & \cdots & {{\boldsymbol{\phi }}}_{{N}_{t}}\end{array}\right],\end{eqnarray} \tag{ 2 }$$

where $J$ is the total voxel number, ${N}_{t}$ is the number of time points (or k-space lines); ${{\bf{U}}}_{{\bf{x}}}\in {{\mathbb{C}}}^{J\times L}$ contains $L$ spatial basis functions, and ${{{\boldsymbol{\Phi }}}_{{\rm{rt}}}\in {\mathbb{C}}}^{L\times {N}_{t}}$ contains real-time temporal weighting functions (depicting relaxation, motion, contrast changes, etc). At a specific time point $t={t}_{s},$ the real-time image ${{\boldsymbol{a}}}_{{t}_{s}}$ is a linear combination of the spatial basis functions, weighted by a vector ${{\boldsymbol{\phi }}}_{{t}_{s}}={\left[\begin{array}{ccc}{\phi }_{1,{t}_{s}} & \cdots & {\phi }_{L,{t}_{s}}\end{array}\right]}^{{\rm{T}}}:$

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{{t}_{s}}=\left[\begin{array}{c}a\left({{\bf{x}}}_{1},{t}_{s}\right)\\ \vdots \\ a\left({{\bf{x}}}_{J},{t}_{s}\right)\end{array}\right]={\phi }_{1,{t}_{s}}{{\boldsymbol{u}}}_{1}+\ldots +{\phi }_{L,{t}_{s}}{{\boldsymbol{u}}}_{L}=\displaystyle \sum _{i=1}^{L}{\phi }_{i,{t}_{s}}{{\boldsymbol{u}}}_{i}.\end{eqnarray} \tag{ 3 }$$

In practice, ${\bf{A}}$ as a whole can be reconstructed by recovering ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ and ${{\bf{U}}}_{{\bf{x}}}$ in a two-step approach (Liang 2007, Pedersen et al 2009, Christodoulou et al 2014, Biswas et al 2015). Typically, temporal weighting functions ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ are first recovered using only 'navigator data' ${{\bf{D}}}_{{\rm{nav}}}\in {{\mathbb{C}}}^{M\times {N}_{{\rm{nav}}}},$ i.e. the central k-space line which is frequently sampled in time, where $M$ is the number of points sampled per k-space line, and ${N}_{{\rm{nav}}}$ is the total number of navigator lines. ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ is often extracted by calculating the singular value decomposition of the frequently sampled navigator data ${{\bf{D}}}_{{\rm{nav}}}$ and selecting the $L$ most significant right singular vectors.

Spatial basis functions ${{\bf{U}}}_{{\bf{x}}}$ are then reconstructed by solving the following problem:

$$\begin{eqnarray}&&{\hat{{\bf{U}}}}_{{\bf{x}}}={{\rm{argmin}}}_{{{\bf{U}}}_{{\bf{x}}}}{\parallel {{\bf{D}}}_{{\rm{im}}}-{\rm{\Omega }}\left({\bf{E}}{{\bf{U}}}_{{\bf{x}}}{{\boldsymbol{\Phi }}}_{{\rm{rt}}}\right)\parallel }_{F}^{2}+\lambda R\left({{\bf{U}}}_{{\bf{x}}}\right),\end{eqnarray} \tag{ 4 }$$

where ${{\bf{D}}}_{{\rm{im}}}\in {{\mathbb{C}}}^{M\times {N}_{{\rm{im}}}}$ denotes the 'imaging data', which is acquired from the entire k-space with sparse sampling schemes, such as randomized Cartesian or golden-angle radial trajectories (${N}_{{\rm{im}}}$ is the total number of imaging lines, ${N}_{{\rm{nav}}}+{N}_{{\rm{im}}}={N}_{t}$), ${\bf{E}}$ is the signal encoding operator, ${\rm{\Omega }}$ is the (k-t)-space undersampling operator, and $R$ is a regularization functional to exploit compressed sensing. Both steps of this reconstruction process are non-causal, and are therefore appropriate for a 'prep' scan but not a real-time, on-the-fly 'live' scan.

Given ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ extracted from the right singular vectors of ${{\bf{D}}}_{{\rm{nav}}},$ there exists a linear transformation ${\bf{T}}$ that maps ${{\bf{D}}}_{{\rm{nav}}}$ to ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}.$ For an individual time point ${t}={t}_{s},$ the navigator data ${{\boldsymbol{d}}}_{{\rm{nav}},{t}_{s}}{\in {\mathbb{C}}}^{M\times 1}$ can therefore be transformed into ${{\boldsymbol{\phi }}}_{{t}_{s}}$ with ${{\boldsymbol{\phi }}}_{{t}_{s}}={\bf{T}}{{\boldsymbol{d}}}_{{\rm{nav}},{t}_{s}}.$ Accordingly, the entire 3D image at $t={t}_{s}$ can be generated with a simple matrix multiplication:

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{{t}_{s}}={{\bf{U}}}_{{\bf{x}}}{{\boldsymbol{\phi }}}_{{t}_{s}}={{\bf{U}}}_{{\bf{x}}}{\bf{T}}{{\boldsymbol{d}}}_{{\rm{nav}},{t}_{s}}.\end{eqnarray} \tag{ 5 }$$

For successive scans using the same sequence, ${{\bf{U}}}_{{\bf{x}}}$ and ${\bf{T}}$ are assumed to remain constant throughout the acquisition process, unless abrupt body motion or unexpectedly introduced contrast mechanisms force the new images outside the range of ${{\bf{U}}}_{{\bf{x}}}.$ Therefore, we developed the SPIDERM technique based on the constant nature of ${{\bf{U}}}_{{\bf{x}}}$ and ${\bf{T}}.$ The 'prep' scan is first applied to learn ${{\bf{U}}}_{{\bf{x}}}$ and ${\bf{T}}.$ Specifically, ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ and ${{\bf{U}}}_{{\bf{x}}}$ are reconstructed in the two-step approach, and ${\bf{T}}$ is calculated as:

$$\begin{eqnarray}&&{\bf{T}}={{\boldsymbol{\Phi }}}_{{\rm{rt}}}{{\bf{D}}}_{{\rm{nav}}}^{+},\end{eqnarray} \tag{ 6 }$$

where ${{\bf{D}}}_{{\rm{nav}}}^{+}$ is the pseudo-inverse of ${{\bf{D}}}_{{\rm{nav}}}.$ Afterwards, the 'live' scan is performed: only ${{\boldsymbol{d}}}_{{\rm{nav}}}$ is acquired, and real-time 3D images can be generated on the fly by applying the matrix multiplication process according to equation (5) (figure 1).

Zoom In Zoom Out Reset image size

In non-steady-state sequences with a periodic signal evolution (e.g. sequence with inversion recovery, saturation recovery, or T2 preparation module, etc), the temporal weighting functions (or the temporal subspace) ${{\boldsymbol{\Phi }}}_{{\rm{rt}}}$ contain the information not only about respiratory motion, but also contrast changes. There is a need in the 'live' scan to separate the motion information and contrast information, so that real-time images can be displayed with a stable contrast weighting, e.g. T1-weighted (T1w), T2-weighted (T2w), or proton-density-weighted (PDw), while maintaining true motion states. This can be achieved by a data-driven image contrast regeneration method, as described below.

With respiratory binning, a multi-bin temporal subspace tensor ${{\rm{\Phi }}\in {\mathbb{C}}}^{L\times {N}_{{\rm{seg}}}\times {N}_{p}}$ can be recovered from the navigator data ${{\bf{D}}}_{{\rm{nav}}}$ with low-rank tensor completion (Christodoulou et al 2018), where ${N}_{{\rm{seg}}}$ is the number of sampling points in a signal evolution cycle, and ${N}_{p}$ is the number of respiratory motion phases.

As the 'live' image comes in, its respiratory phase $p$ is identified with the liver-dome position using nearest-neighbor matching to the 'prep' data, and ${{{\boldsymbol{\Phi }}}_{p}\in {\mathbb{C}}}^{L\times {N}_{{\rm{seg}}}}$ is extracted from ${\rm{\Phi }}$ at phase $p.$ Then the target image contrast that would correspond to time $t={t}_{c}^{({\rm{target}})}$ can be synthesized in the 'live' scan by replacing ${{\boldsymbol{\phi }}}_{{t}_{s}}$ in equation (5) with ${\tilde{{\boldsymbol{\phi }}}}_{{t}_{s}}:$

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{\phi }}}}_{{t}_{s}}\left({t}_{c}^{({\rm{target}})}\right)={{\boldsymbol{\phi }}}_{{t}_{s}}+{\rm{\Delta }}{{\boldsymbol{\phi }}}_{p,{t}_{s}}\left({t}_{c}^{\left({\rm{target}}\right)}\right)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{{\boldsymbol{\phi }}}_{p,{t}_{s}}\left({t}_{c}^{\left({\rm{target}}\right)}\right)={{\boldsymbol{\Phi }}}_{p}\left(:,{t}_{c}^{\left({\rm{target}}\right)}\right)-{{\boldsymbol{\Phi }}}_{p}\left(:,{t}_{c}^{\left({\rm{original}}\right)}\right),\end{eqnarray} \tag{ 8 }$$

where ${t}_{c}^{\left({\rm{original}}\right)}$ refers to the original time point within the signal evolution cycle corresponding to the absolute time point ${t}={t}_{s},$ while ${t}_{c}^{\left({\rm{target}}\right)}$ refers to the time point with targeted contrast of interest, e.g. T1w, T2w, or PDw. This new term subtracts the contribution of the current contrast weighting, ${{\boldsymbol{\Phi }}}_{p}\left(:,{t}_{c}^{\left({\rm{original}}\right)}\right)\in {{\mathbb{C}}}^{L\times 1},$ and replaces it with the desired contrast weighting, ${{\boldsymbol{\Phi }}}_{p}(:,{t}_{c}^{\left({\rm{target}}\right)})\in {{\mathbb{C}}}^{L\times 1}.$

Thus, equation (5) can be adapted as follows for contrast-frozen, motion-maintained real-time imaging:

$$\begin{eqnarray}&&{\tilde{{\boldsymbol{a}}}}_{{t}_{s}}={{\bf{U}}}_{{\bf{x}}}{\tilde{{\boldsymbol{\phi }}}}_{{t}_{s}}={{\bf{U}}}_{{\bf{x}}}\left({\bf{T}}{{\boldsymbol{d}}}_{{\rm{nav}},{t}_{s}}+{\rm{\Delta }}{{\boldsymbol{\phi }}}_{p,{t}_{s}}\left({t}_{c}^{\left({\rm{target}}\right)}\right)\right).\end{eqnarray} \tag{ 9 }$$

Note that the original 'live' image ${{\boldsymbol{a}}}_{{t}_{s}}={{\bf{U}}}_{{\bf{x}}}{\bf{T}}{{\boldsymbol{d}}}_{{\rm{nav}},{t}_{s}}$ is generated directly from the navigator data without binning, as in equation (5). The binning process is only used to generate the contrast update term ${\rm{\Delta }}{{\boldsymbol{\phi }}}_{p,{t}_{s}},$ which ensures better separation of motion and contrast changes.

When using a pulse sequence in which multiple image contrasts present along the signal evolution, SPIDERM is able to generate multi-contrast real-time images using several different ${t}_{c}^{\left({\rm{target}}\right)}$'s. In this work, three different ${t}_{c}^{\left({\rm{target}}\right)}$'s were chosen to represent T1w, T2w, and PDw respectively, as indicated by black arrows in figure 2(b).

Zoom In Zoom Out Reset image size

For simplicity, ${{\boldsymbol{a}}}_{{t}_{s}}$ generated with equation (5) are denoted as contrast-varying (CV-SPIDERM) images, while ${\tilde{{\boldsymbol{a}}}}_{{t}_{s}}$ generated with equation (8) are denoted as contrast-frozen (CF-SPIDERM) images.

Built upon the partial separability model (Liang 2007), the recently proposed MR Multitasking technique (Christodoulou et al 2018) can generate multi-contrast MR images with 3D coverage and high spatiotemporal resolution, without the assistance of external devices for gating or triggering (Hu et al 2020, Wang et al 2020, Han et al 2021). We tested the SPIDERM technique using an abdominal T1/T2 MR Multitasking sequence (figure 2) (Deng et al 2019). A saturation recovery preparation and a T2 preparation are used to generate T1 and T2 contrast weightings, respectively, during each TR (figure 2(b)). A gap of 700 ms is intended to facilitate magnetizations' full recovery and minimize T1 weighting in subsequent PDw and T2w acquisitions. The k-space is continuously sampled with fast low-angle shot (FLASH) readouts using a stack-of-stars acquisition with golden angle ordering in-plane and Gaussian-density randomized ordering in the partition direction (figure 2(a)). The 'imaging data' is interleaved with 'navigator data' (central k-space line along the partition direction, ${k}_{x}={k}_{y}=0$) every 10th readout. The acquired data contains three overlapping dynamics, including respiratory motion, T1 relaxation, and T2 relaxation. General imaging parameters for both phantom and volunteer studies were: axial orientation, TR/TE = 6.0/3.1 ms, flip angle = 5° (following SR preparation) and 10° (following T2 preparation), bandwidth = 762 Hz/pixel, water-excitation for fat suppression, BIREF T2 preparation of 42 ms. The time per scan was 8 min.

The feasibility of SPIDERM was evaluated using an open-source digital phantom in MATLAB (https://github.com/SeiberlichLab/Abdominal_MR_Phantom) (Lo et al 2019). First, k-space data of 60000 readouts (corresponding to a total scan time of 8 min) were simulated using the T1/T2 Multitasking sequence and the sampling pattern shown in figure 2. Additional imaging parameters included: matrix size = 320 × 320 × 40, and voxel size = 1.7 × 1.7 × 6.0 mm³. To demonstrate that our technique does not assume or rely on strictly periodic respiratory cycles, time-varying breathing patterns were simulated with pseudo-randomly interleaved normal (∼4.2 s), long (∼6.3 s) and short (∼3.2 s) respiratory cycles.

The data of 60 000 readouts were viewed as the 'prep' scan, from which the spatial basis ${{\bf{U}}}_{{\bf{x}}},$ the linear transformation matrix ${\bf{T}},$ and the multi-bin temporal subspace tensor ${\rm{\Phi }}$ were generated. Then, 1000 additional time points of navigator data were simulated with variable respiratory motion positions and contrasts. These data were viewed as a portion of the 'live' scan and processed with the SPIDERM technique to generate real-time multi-contrast 3D images.

The in-vivo study was approved by local Institutional Review Board and written informed consent was obtained from all participating subjects. Experiments were performed in eight healthy subjects on a 3.0 T clinical scanner (Biograph mMR, Siemens Healthineers, Erlangen, Germany) equipped with an 18-channel phase array body coil. Additional imaging parameters included: matrix size = 320 × 320 × 52, field-of-view (FOV) = 550 × 550 × 312 mm³, voxel size = 1.7 × 1.7 × 6.0 mm³.

For each subject, two identical T1/T2 MR Multitasking scans were performed successively, serving as a 'prep' scan and a 'live' scan, respectively. Volunteers were instructed to breathe normally during the two scans. Although repetitive acquisition of the same single k-space projection as navigator data is the only essential need in the 'live' scan, we adopted the same sequence as used in the 'prep' scan for the following purposes: (a) the navigator data acquired every 10th readout were used to generate CV-SPIDERM images with equation (5), and CF-SPIDERM images with equation (9), using the proposed SPIDERM technique; (b) the data of the entire scan, including both navigator data and imaging data, were reconstructed retrospectively as in the 'prep' scan, to generate 'reference' real-time contrast-varying (CV-ref) images. Figure 3 illustrates the experimental design of in-vivo studies.

Zoom In Zoom Out Reset image size

Image reconstruction was performed offline in MATLAB 2018a on a Linux workstation equipped with two 2.7-GHz 12-core Intel Xeon CPUs, one NVIDIA Quadro K6000 GPU, and 256 GB RAM.

3.4.1. Digital phantom study

3.4.2.  In-vivo studies

3.4.2.1. Image quality

To assess the image quality of SPIDERM images, CV-SPIDERM and CV-ref images were compared using the following quantitative metrics. The mean values of these metrics among all 6000 time points in the live scan were reported for each subject. The average metrics of each signal evolution cycle were also measured and then plotted as a function of the signal evolution cycle to illustrate their temporal stability over the 8 min scan.

• (1)  
  Normalized root mean square error (NRMSE)

  $$\begin{eqnarray}&&{\rm{NRMSE}}=\sqrt{\frac{{\sum }_{i}{\left({I}_{{\rm{SPIDERM}}}-{I}_{{\rm{ref}}}\right)}^{2}}{{\sum }_{i}{I}_{{\rm{ref}}}^{2}}},\end{eqnarray} \tag{ 10 }$$

  where ${I}_{{\rm{SPIDERM}}}$ and ${I}_{{\rm{ref}}}$ denote the magnitude pixel values of CV-SPIDERM and CV-ref images of the 'live' scan respectively.
• (2)  
  Peak signal-to-noise ratio (PSNR)

  $$\begin{eqnarray}&&{\rm{PSNR}}=-10{\rm{lg}}\left(\frac{{\rm{MSE}}}{{I}_{{\rm{\max }}}^{2}}\right),\end{eqnarray} \tag{ 11 }$$

  where ${I}_{{\rm{\max }}}$ is the maximum magnitude pixel value, and ${\rm{MSE}}$ denotes the mean squared error.
• (3)  
  Structural similarity index (SSIM)

$$\begin{eqnarray}&&{\rm{SSIM}}=\frac{\left(2{\mu }_{{\rm{ref}}}\cdot {\mu }_{{\rm{SPIDERM}}}+{c}_{1}\right)\left(2{\sigma }_{{\rm{ref}},{\rm{SPIDERM}}}+{c}_{2}\right)}{\left({\mu }_{{\rm{ref}}}^{2}+{\mu }_{{\rm{SPIDERM}}}^{2}+{c}_{1}\right)\left({\sigma }_{{\rm{ref}}}^{2}+{\sigma }_{{\rm{SPIDERM}}}^{2}+{c}_{2}\right)},\end{eqnarray} \tag{ 12 }$$

where ${\mu }_{\left(\cdot \right)}$ and ${\sigma }_{\left(\cdot \right)}$ denotes the mean and variance respectively, ${\sigma }_{{\rm{ref}},{\rm{SPIDERM}}}$ is the covariance, ${c}_{1}=0.01,$ ${c}_{2}=0.03.$

3.4.2.2. Accuracy of motion depiction

The respiratory motion-induced displacement of the liver dome was determined in the same coronal view of SPIDERM images (CV-SPIDERM, CF-SPIDERM T1w, and CF-SPIDERM PDw) and CV-ref images. CF-SPIDERM T2w images were not included in this comparison because the banding artifact at the liver dome caused by field inhomogeneity during T2 preparation made it difficult to accurately measure the distance (see Discussion for more details). The first 1000 time points of the 8 min 'live' scan were selected in each volunteer for linear regression analysis.

Figure 4 shows the comparison of true reference, CV-SPIDERM, and CF-SPIDERM images of the digital phantom. Representative time points shown in figures 4(a) and (b) corresponded to end-of-expiration and end-of-inspiration, respectively. The respiratory motion-induced displacement was visually comparable between SPIDERM images and the true reference. CF-SPIDERM T1w, T2w and PDw images showed appropriate contrasts respectively. As shown in figure 5, the displacement of the liver dome measured from CV-SPIDERM images and reference images were strongly correlated in all three different respiratory cycles (normal cycles: slope = 0.90, intercept = 1.48, R² = 0.984; long cycles: slope = 0.90, intercept = 1.43, R² = 0.991; short cycles: slope = 0.88, intercept = 2.35, R² = 0.983).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The pancreas contour analysis showed a Dice similarity coefficient of 0.91, 0.84 and 0.85 for EOE, MED and EOI time points respectively. Mean surface distances were reported as 0.57, 0.95 and 0.87 mm for EOE, MED and EOI respectively.

The average elapsed time from the input of the central k-space line to the generation of real-time contrast-frozen 3D images was approximately 45 ms. Given that the navigator data (one central k-space line) can be acquired in 6 ms, a real-time display latency of 55 ms or less can be reached. This is achieved using MATLAB 2018a on a Linux workstation equipped with two 2.7-GHz 12-core Intel Xeon CPUs.

The NRMSE, PSNR and SSIM values between CV-SPIDERM and CV-ref images for each volunteer are shown in table 1. The average NRMSE, PSNR and SSIM among eight volunteers were 0.141, 30.12 and 0.88 respectively.

Figure 6 shows the time courses of various quantitative image quality metrics during the 8 min live scan. Compared with the first 25 signal evolution cycles, the average NRMSE, PSNR and SSIM changed by +9.1%, −3.2% and −0.9% respectively in the last 25 signal evolution cycles, indicating a slight metric degradation. Abrupt increase in NRMSE or decrease in PSNR/SSIM were visible in some volunteers (such as Volunteer 3), which presumably arose from sudden deep breaths.

Zoom In Zoom Out Reset image size

Figure 7 shows the comparison between reference images and SPIDERM images in in-vivo studies. As shown in figures 7(a) and (b), CV-SPIDERM images had comparable image quality with CV-ref images, and CF-SPIDERM images demonstrated appropriate image contrasts of T1w, T2w and PDw. Banding artifacts caused by imperfect T2 preparation and main field inhomogeneity were visible at the liver dome area, making CF-SPIDERM T2w images' quality suboptimal (see Discussion for more details). A movie showing several respiratory cycles of the same volunteer can be found in supplementary video S1 (available online at stacks.iop.org/PMB/67/135008/mmedia).

Zoom In Zoom Out Reset image size

Figure 8 displays CV-ref images and CF-SPIDERM (T1w, T2w and PDw) images of another volunteer within an arbitrarily selected respiratory cycle. Motion states were consistent between CF-SPIDERM images and CV-ref images, while image contrasts remained stable in CF-SPIDERM images among different motion states.

Zoom In Zoom Out Reset image size

As shown in figure 9, the displacement of the liver dome measured from SPIDERM images and reference images were strongly correlated (CV-SPIDERM images: slope = 0.98, intercept = 1.29, R² = 0.986; CF-SPIDERM T1w: slope = 0.98, intercept = 1.32, R² = 0.983; SPIDERM PDw images: slope = 0.98, intercept = 1.59, R² = 0.983).

Zoom In Zoom Out Reset image size