share
1 Introduction
Magnetic Resonance Imaging (MRI) is the most widely used medical imaging modality employed to produce high quality detailed images of both the anatomy and physiology of the human body. The fact that it allows for non-invasive clinical examinations and diagnosis, has played a significant role in its widespread use. MRI however, suffers from very long acquisition times as a result of the fact that MRI data (known as -space or Fourier domain) are acquired sequentially in time, frequently causing patient distress. Additionally, these long acquisition times make it nearly intractable to deliver dynamic treatment such as image-guided targeted radiotherapy in real time.
Truncating the MR imaging time, referred to as accelerating MRI [zbontar2019fastmri], is pivotal not only for improving patient comfort, decreasing scanning time and reducing medical costs, but, it can also make image-guided treatment feasible. Over the past fifteen years accelerating the MRI acquisition has been in the center of many studies yielding advancements such as Parallel Imaging (PI) and Compressed Sensing (CS).
|
|
|
|
In Parallel Imaging multiple radio-frequency receiver coils are utilised to read-out partial -space measurements that depend on the position of each coil [parallel-https://doi.org/10.1002/nbm.1042]. A unique sensitivity map corresponds to each individual receiver coil and can be estimated with different techniques such as SENSE [https://doi.org/10.1002/cmr.a.21244, senseencoding]. An example of a reconstructed brain image with the root-sum-of-squares (RSS) method acquired from multiple coils is illustrated in Figure 1.
Compressed Sensing is a technique used in clinical settings that aims in accelerating the MRI acquisition by acquiring sparse or partial -space measurements [cs-1580791, cs-https://doi.org/10.1002/mrm.21391, cs-lustig-2, cs-https://doi.org/10.1002/cpa.20124, https://doi.org/10.1002/mrm.22428, 1614066]. In CS an acceleration factor is defined (-fold acceleration) as the ratio of the number of the fully sampled -space measurements to the number of measurements of the partially sampled -space. Reconstructing partially observed -space data violates the Nyquist sampling criterion [doi:https://doi.org/10.1002/9781118633953.ch12], therefore leading to aliased and/or blurred reconstructed images, also known as aliasing artifacts which in many cases can be identified but in others can be interpreted as pathology.
In recent years, various Deep Learning approaches have been deployed aiming in reconstructing partially observed MRI data achieving state-of-the-art performance [Aggarwal_2019, hammernik2017learning, zhang2019reducing, Hyun_2018, 8425639, yu2017deep, sriram2020grappanet, schlemper2017deep, han2019kspace]. The majority of these approaches train DL models on retrospectively subsampled measurements by applying rectilinear subsampling schemes (discussed further in Section 2) which even with the use of DL produce imperfect reconstructions.
In this work we investigate the performance of accelerated image reconstruction of a particular DL architecture (RIM) on retrospectively subsampled -space data via two subsampling schemes, rectilinear and radial (Section 3.1). We argue and show that the model trained on data subsampled with the latter scheme, can produce more precise and faithful reconstructions. In summary, the contributions of this paper are the following:
-
[noitemsep]
-
We are the first to employ the RIM architecture for training with radially subsampled data.
-
We show through evaluation that the model we trained utilising radial subsampling produces reconstructions with higher fidelity than the model trained with conventional rectilinear subsampling.
-
We pave the way for other DL approaches to involve prospective or retrospective radial subsampling.
In Section 2, we provide a mathematical formulation of the accelerated MRI reconstruction along with some information about prospective and retrospective subsampling. Section 3 focuses on our experiments implementation and demonstrates some of our results. In Section 4 we conclude with a discussion. Our code is available at https://github.com/directgroup/direct.
2 Background
2.1 Mathematical Formulation
In this section we briefly define the problem of reconstructing accelerated MRI data from multiple coils mathematically. Let denote the true image. For a number of coils, let denote the fully acquired -space measurements and sensitivity maps from coil . Then,
| (1) |
where
denotes the discrete Fast Fourier Transform (FFT) and
the element-wise product. In PI-CS, the RSS reconstruction is commonly used to reconstruct an image from the multiple -space measurements [https://doi.org/10.1002/mrm.1910160203] given by| (2) |
The problem of recovering the true image from accelerated MRI data is an inverse problem with a forward model formulated as
| (3) |
where denote the accelerated measurements, the accumulated noise from coil and a subsampling mask (see 2.2). A solution to the inverse problem can be formulated as a variational problem
| (4) |
where and denote data-fidelity and regularisation terms, respectively [Kaipio:1338003, mathmodelsmri]
. This is equivalent to the Bayesian statistical problem of finding the maximum a posteriori (MAP) estimation
[4407762, putzky2017recurrent, Kaipio:1338003]:| (5) |
where denotes the posterior distribution, and the -space likelihood and prior models, respectively. Let denote the complex conjugate of of . Assuming that
are i.i.d. and derived from a normal distribution with variance
, the negative log-likelihood and its gradient [putzky2017recurrent, LONNING201964] are given by| (6) |
2.2 k-space Subsampling
|
|
|
Prospective Subsampling
In clinical settings several sampling (and subsampling) trajectories have been used since the advent of MRI [kspace_strategies, K-Spaceintheclinic]. The most conventional sampling trajectory is the Cartesian or rectilinear sampling [undersampling_patterns_in_kspace], that is, the -space is measured in a line-by-line scheme on an equidistant rectangular grid (Figure 1(a)). Other non-Cartesian -space trajectories commonly used are the radial and spiral (Figures 1(b) and 1(c)). Since non-Cartesian trajectories are not acquired on an equidistant grid, one of projection reconstruction [https://doi.org/10.1118/1.1677252, https://doi.org/10.1002/mrm.1910190226], regridding [PMID:18243972] or the application of the non-uniform FFT (NUFFT) is required.
Retrospective Subsampling
|
|
To simulate CS subsampling, subsampling masks are retrospectively applied onto the fully sampled -space to produce a subsampled -space hereby referred to as the masked or zero-filled -space. More formally, let denote a binary subsampling mask such that if is masked only if , . This is described by Equation 3. For each subsampling mask an acceleration factor is determined such that . In this paper, we work with data acquired by Cartesian trajectories and thus, for our experiments we generated rectilinear and radial subsampling masks on the Cartesian grid as depicted by Figures 2(a) and 2(b).
3 Experiments
3.1 Implementation
Dataset
The RIM models were trained and evaluated on the Calgary-Campinas public brain multi-coil MRI dataset which was released as part of an accelerated MRI reconstruction challenge [calgary]. The dataset is consisted of 67 three-dimensional raw -space volumes that are T1-weighted, gradient-recalled echo, 1 mm isotropic sagittal and are collected on a Cartesian grid (equidistant). These amount to 17,152 slices of fully sampled -spaces which we randomly split into training (40 volumes), validation (14 volumes) and test (13 volumes) sets. All data were acquired from twelve receiver coils () and the reconstructed images have or pixels.
Subsampling
For the purposes of the experiments in this paper, two different subsampling schemes were used to retrospectively subsample the fully acquired -space data. The first scheme generated rectilinear subsampling masks by first sampling a fraction (10% when or 5% when ) of the central phase encoding columns (low frequencies) and then randomly including others up to a level of acceleration [zbontar2019fastmri]. The second subsampling scheme was radial subsampling in which fully acquired Cartesian measurements were subsampled in a radial fashion to simulate radial acquisition. Radial masks were generated using the CIRcular Cartesian UnderSampling (CIRCUS) method presented in [circus] which fit our settings. Figure 3 illustrates an example of one of each of the subsampling masks we use. Note that the same subsampling mask is applied to all coil data and slices of a single volume sample.
Model Architecture
To build, train and evaluate our models, we exploited our Deep Image Reconstruction (DIRECT) toolkit [DIRECTTOOLKIT]
which uses PyTorch
[NEURIPS2019_9015]. The architecture we opt for is the one of the Recurrent Inference Machine (RIM) as it has been shown to achieve state-of-the art performance in accelerated MRI reconstruction tasks [putzky2017recurrent, putzky2019irim, LONNING201964, beauferris2020multichannel]. At each training iteration RIMs are trained to resemble a gradient descent scheme to find the MAP estimation (Eq. 5) using Eq. 6 by following an iterative scheme of time-steps:| (7) |
where is called the internal state. The update function of the RIM is a sequence of two blocks, each consisted of a convolutional layer (,
kernel sizes) followed by a rectified linear unit (ReLU) activation and a convolutional gated recurrent unit (convGRU)
[shi2015convolutional], followed by a convolution. At time-step , produces an incremental update for and the internal state for time-step .Model Training
 |
 |
For our experiments we trained two RIMs. The first was trained on the radially subsampled data (Radial-RIM) and the second on the rectilinear subsampled data (Rect-RIM) employing the two subsampling schemes presented in 3.1 Subsampling. For both RIMs we set time-steps and used 128 hidden channels for the convolutional layers. In our implementation an initialiser for the internal state was learnt throughout training. We also employed a U-net [ronneberger2015unet] which was trained simultaneously with each RIM that sequentially improved the coil sensitivity maps. The total number of trainable parameters for each RIM amounted to 360,194 parameters and for each U-Net to 484,898 parameters. Models took 2D slices of retrospectively subsampled -space data as input and output 2D reconstructed images. Models were optimised using the Adam optimiser and were trained for approximately 150k iterations with a batch size of 4 slices. A warm-up schedule [goyal2018accurate] was used to linearly increase the learning rate to 0.0001 over 1k warm-up iterations, before decaying with a ratio of 5:1 every 5k iterations. Throughout training, the acceleration factor was arbitrarily set to 5 or 10. Every 500 iterations, the models were evaluated on the 5-fold and 10-fold accelerated validation data. Coil sensitivities were initially calculated using autocalibration data, that is, the center of the subsampled -space [https://doi.org/10.1002/mrm.24751]. Figure 4 illustrates an example of the autocalibration signal (ACS) used to calculate the sensitivity maps corresponding to the masks depicted in 2(a) and 2(b). Figure 5
provides an overview of our model training with radial subsampling. The loss function used for training the models was the sum of the
loss and a loss derived from the SSIM metric [ssim]. Assuming that the RSS reconstructed image is the ground truth image and is the model output, the loss function employed was the following| (8) |
3.2 Results
To evaluate the performance of the models and the quality of the reconstructions, as well as to compare the two distinct subsampling schemes, we used three common metrics; the peak signal-to-noise ratio (PSNR), the structural similarity index meseaure (SSIM) and the visual information fidelity (VIF)
[ssim, zbontar2019fastmri, calgary, 1576816, 5596999]. Note that in inference time, only the last time-step is used as the model prediction. We ran inference on the test data which were retrospectively 5-fold or 10-fold subsampled. For each acceleration and each type of subsampling we selected to present results of the best model checkpoint based on the validation performance on the SSIM metric.The quantitative results are summarised in Table 1
. Specifically, we report the mean metrics along with their respective standard deviations. Some qualitative results are shown by Figures
6 and 7 which illustrate exemplary reconstructions of axial images of brain samples from the test data. In the former are visualised image reconstructions outputted from the models, against the ground truth images and zero-filled reconstructions. The latter depicts 5, 10 accelerated image reconstructions along with error maps from the ground truth images.| Acceleration | Subsampling | Metrics | ||
|---|---|---|---|---|
| PSNR | SSIM | VIF | ||
| 5-fold | Rectilinear | |||
| Radial | 34.836 0.489 | 0.939 0.004 | 0.955 0.006 | |
| 10-fold | Rectilinear | |||
| Radial | 31.144 0.655 | 0.899 0.009 | 0.905 0.017 | |
As shown by Table 1
, Radial-RIM consistently outperformed Rect-RIM for both acceleration factors, as well as, the standard deviations of the evaluation metrics of the Radial-RIM were lower. The superiority in performance of Radial-RIM is also verified by the qualitative results as its predictions better resemble the ground truth images as they are more faithful with significantly less artifacts compared to Rect-RIM’s. The zoomed regions in Figures
6 and 7 show that it’s almost impossible to distinguish which is the true image between the ground truth images and Radial-RIM’s reconstructions. On the other hand, Rect-RIM’s reconstructions of the 10-fold accelerated data exemplified notable artifacts which in fact might be interpreted as pathology. Examining closely the error maps in Figure 7, Radial-RIM’s prediction error maps resemble random noise while in Rect-RIM’s prediction error maps the anatomical structure can be easily inferred as a result of the lower quality reconstructions. Another noteworthy remark of this experiment is the fact that Radial-RIM managed to infer slightly better reconstructions from the 10-fold radially subsampled data than the Rect-RIM did from 5-fold rectilinear subsampled data.
4 Conclusion - Discussion
In this work we have established that training a DL model such as Recurrent Inference Machines on retrospectively radially subsampled data can outperform the traditional approach of rectilinear subsampling. The high incoherence and variable density sampling of the radial subsampling pattern might have contributed to its superior performance. The data we have used in this study were fully acquired by a conventional Cartesian trajectory, therefore, measurements were placed on an equidistant grid of pixel size. In our experiments we retrospectively subsampled the data exploiting a radial subsampling scheme as a way of simulating radial (sub)sampling in real clinical settings. Even though differences exist between the view orderings of the retrospectively and prospectively radially subsampled data, their influence requires further investigation. Additionally, in practice prospective radial (sub)sampling patterns take measurements on a non-Cartesian grid, as depicted by Figure 1(b). As explained in Section 2, this process would require regridding or the application of NUFFT. That is why our results pave the way for incorporating radial subsampling in clinical settings, which combined with Deep Learning can speed up MRI acquisitions and even allow for dynamic treatment.