Dynamic contrast-enhanced MRI (DCE-MRI) is an imaging method which is used to study microvascular structure and tissue perfusion. The method has many applications including blood-brain-barrier assessment after acute ischemic stroke Merali2017; Villringer2017 and treatment monitoring of breast cancer Martincich2004; Pickles2005 and glioma Piludu2015. The operation principle of DCE-MRI is to inject a bolus of gadolinium based contrast agent into the blood stream, and acquire a time series of MRI data with a suitable -weighting to obtain a time series of 2D (or 3D) images which exhibit contrast changes induced by concentration changes of the contrast agent in the tissues.
The analysis of the contrast agent dynamics during imaging requires high resolution in space and in time; the high temporal resolution is necessary to measure when the contrast is passing through the artery and it is used for determining the subject specific Arterial Input function (AIF), and the high spatial resolution is necessary to adequately capture boundaries of perfused tissues. In many cases, sufficient time resolution can only be obtained by utilizing an imaging protocol during which only partial k-space sampling can be obtained for each image in the time series. However, acquiring less samples than required by the Nyquist criterion makes the image reconstruction problem ill-posed
(and non-unique), causing artifacts and deterioration of the image quality when conventional reconstruction methods such as inverse Fourier transform or re-gridding are employed.
According to the theory of compressed sensing (CS) Cand`es2006a; Cand`es2006; Cand`es2006b
, images that have a sparse representation can be recovered from undersampled measurements of a linear transform, i.e. sampling rate below Nyquist rate, using appropriate nonlinear reconstruction algorithms and appropriate (random) sampling of the data space. The compressibility of MR images and the fact that MR scanner measures samples of a linear transform of the unknown image (the k-space samples can be mathematically considered as Fourier coefficients) suggest that the idea of CS is applicable to MR imaging, offering thus a potentially significant scan time reduction without sacrificing the image quality. Since the seminal work of Candès, Romberg and Tao in 2006, CS has been extensively applied to MRI. In 2007 CS was applied to MRI inlustig2007_sparseMRI and in 2018 it received FDA approval for clinical use.
The undersampling of the k-space for speeding up the dynamic MRI data acquisition can in principle be done in many different ways. However, for many applications, the low frequency features that are present in the center of the k-space are of importance. In Zong2014, for example, it was demonstrated that center weighted random sampling patterns were preferable to purely random sampling of the k-space within the CS approach. Radial sampling has the advantage that the center of the k-space is sampled densely, even when the sampling (i.e., number of radial spokes) is remarkably reduced.
In this paper, we consider image reconstruction problem of DCE-MRI with sparsely sampled golden angle radial data, where the angle of subsequent spokes is . The number of measurement spokes used for reconstructing a single time frame is chosen to be a Fibonacci number, which was shown to be an optimal choice in Winkelmann2007. This type of sparse sampling between the time frames results in each time frame having different spokes, and eventually the spokes cover the whole k-space.
We also combine the Golden Angle (GA) sampling with concentric squares sampling. This sampling strategy resembles the linogram method Edholm1987 developed for computed tomography imaging, but the angles of subsequent spokes were chosen according to the golden angle method as opposed to the angles being equidistant in in linogram sampling. Unlike the conventional radial sampling pattern with spokes of equal length, the concentric squares sampling strategy also covers the corners of the k-space. The sampling pattern therefore also collects information of the high frequencies in the corners of the k-space, which leads to a reduction of artifacts originating from the lack of sampling in the corners.
To overcome the ill-posedness of our image reconstruction problem, we use a variational framework. Thus we solve the following optimization problem
where denotes the image sequence, the data sequence, the number of time frames, the data-fitting term, the spatial regularization, the temporal regularization and and are the spatial and temporal regularization parameters, respectively. The selection of the regularization parameters is crucial in terms of resulting image quality. There exists several proposed parameter selection methods, but most of them require the computation of a large number of reconstructions with varying parameters. Having to fix two regularization parameters makes the selection even more time consuming.
In this work, we propose to use the S-curve method hamalainen2013; Niinimaki2013; Kolehmainen2012; niinimaki2015 for automatic selection of the spatial regularization parameter . The idea in the S-curve method is to select the regularization parameter so that the reconstruction has a priori defined level of sparsity in the chosen transformation domain. In DCE-MRI, a reliable a priori estimate for the sparsity level can be extracted from an anatomical MRI image which is based on full-sampling of the k-space and is always taken as part of the MRI measurement protocol but is usually used only for visualization purposes. For the selection of the spatial regularization parameter , we employ one time frame of the GA data from the baseline measurement before the contrast agent administration. After fixing the spatial regularization parameter, we compute dynamic reconstructions with several values of parameter and select a suitable temporal regularization parameter manually. Furthermore we study the performance of three different temporal regularization functionals, namely temporal smoothness, temporal total variation and total generalized variation.
The proposed method is evaluated using simulated GA DCE-MRI data from a rat brain phantom. The results are compared to re-gridding approach, which is the most widely used non-iterative algorithm for reconstructing images from non-Cartesian MRI data. Our re-gridding method was developed in IR4M UMR8081, CNRS, Université Paris-Sud using Matlab. This re-gridding approach does not need additional density correction and it was first used in Kusmia2013, see also iddn_codeGG.
2 Image reconstruction in radial DCE-MRI
2.1 Forward problem
The forward problem in 2D MRI can be modelled for most measurement protocols by the Fourier transform
where is the image domain is the unknown image, is the measured data, and denotes the k-space trajectories. In the discrete framework, the Fourier transform is typically approximated with the multidimensional FFT when using cartesian k-space trajectories and with the non-uniform FFT when using non-cartesian k-space trajectories.
In this work, we consider non-uniform k-space trajectories and approximate the Fourier transform by the non-uniform fast Fourier transform (nuFFT) operator FesslerSutton. We discretize our functions as follows; temporal direction is divided into a sequence of
(vectorized) imagesand data vectors , where each and , respectively. The number of data per frame is equal to the number of GA spokes per frame times the number of samples per spoke. The number of image pixels is . Thus, using nuFFT we re-write (2) in a discretized form at time as
is an interpolation matrix between Cartesian k-space and non-cartesian k-space,is the 2D FFT operation and is a scaling matrix.
2.2 Inverse problem of dynamic image reconstruction
The dynamic inverse problem related to the equation (3) is: given measurement time series and the associated k-space trajectories, solve the unknown images . n. To recover from , we define the inverse problem as the optimization problem
where denotes the spatial regularization functional, the spatial regularization parameter, the temporal regularization parameter and the temporal regularization functional.
In this work, we study the applicability of S-curve method for selecting the spatial regularization parameter . Once has been fixed, the temporal regularization parameter is then selected manually by computing estimates with different values of . Our minimization problem is based on -data fidelity term for the measurement model and we use spatial total variation regularization for promoting sparsity of the gradients of each image Rudin1992. Furthermore we study the performance of three different temporal regularization functionals for promoting temporal regularity of the image series. Our spatio-temporal image reconstruction problem thus writes
where the isotropic 2D spatial total variation norm for complex valued image is defined by
where and denoting the real and imaginary parts of respectively and and the discrete forward first differences in horizontal and vertical directions, respectively.
Temporal regularization 1: Temporal smoothness (TS)
The temporal smoothness regularization (hereafter referred as TS) is defined as the norm of forward first differences in time:
This model promotes smooth slowly changing signals, and it has been used in Adluru2007 for radial DCE myocardial perfusion imaging. TS regularization was compared with temporal TV regularization in the same application in Adluru2008.
Temporal regularization 2 : Temporal total variation (TV)
Temporal total variation (hereafter referred as TV) is defined by the norm of the forward first differences in time:
The temporal total variation model promotes sparsity of the time derivative of the pixel signals, being a highly feasible regularization model for reconstruction of piece-wise regular signals which may exhibit large jumps. The smoothed form of temporal total variation was used in Adluru2009 for multislice myocardial perfusion imaging.
Temporal regularization 3: Total generalized variation (TGV)
The total generalized variation model Bredies2010 is a total variation model that is generalized to higher order differences. Here we use the second-order total generalized variation, which in the discrete 1-dimensional form is of the form
where is an auxiliary vector and is a parameter, which balances the first and second order terms, and is set to here.
This functional balances between minimizing the first-order and second-order differences of the signal. The difference with TV regularization is most clear in smooth regions where piecewise linear solutions are favored over the piecewise constant solutions of TV. From hereafter this temporal regularization is referred as TGV.
TGV was first used in MRI as a spatial prior in Knoll2011, and it has also been used in DCE-MRI as a temporal prior in Wang2017, where different temporal priors were compared in cartesian MRI of the breast.
2.3 Regularization parameter selection
Spatial regularization parameter selection
The spatial regularization parameter is selected using the S-curve method, originally proposed in hamalainen2013; Niinimaki2013; Kolehmainen2012; niinimaki2015, but here modified for TV regularization.
Assume that we have an a priori estimate for the total variation norm of the unknown function. In practice we can use an anatomical image of the same slice in order to obtain a reliable estimate for . Such an anatomical image is practically always acquired as part of the DCE MRI acquisition experiment but usually only used for visualization purposes. However, if such an anatomical image was not acquired, we could, in case of GA acquisition, estimate the expected sparsity level from a conventional reconstruction of a long sequence of baseline data taken before the contrast agent injection. Or the anatomical image could be estimated from the entire data-acquisition similarly as the composite image in Mistretta2006.
Now, given the estimate we select the regularization parameter using the S-curve method as follows
Take a sequence of regularization parameters ranging on the interval such that
Compute the corresponding estimates .
With DCE-MRI data, reconstructions are computed as follows; we take the data that correspond to the first time frame, i.e. which has number of elements equal to the number of GA spokes per frame times the number of points per spoke. We reconstruct for given value by
Here it is important that is taken to be so small that the problem is under regularized and the corresponding reconstruction results to a very noisy image with a big TV-norm value and is taken so large that the problem is over regularized and TV norm of reconstruction is very close to zero.
Compute the TV-norms of the recovered estimates .
Fit a smooth interpolation curve to the data and use the interpolated sparsity curve to find the value of for which . For the interpolation we use Matlab’s interp function and we interpolate our original S-curve to a more dense discretization of the regularization parameter .
2.4 Temporal regularization parameter selection
Once the spatial regularization parameter has been fixed using the S-curve, the temporal regularization parameter can be tuned by computing estimates with different values of and selecting a suitable value manually, for example, by visual assessment of the results.
In this work, we compute the results with three different temporal regularization models using simulated measurement data. Since we consider a simulated test case where a ground truth is available, we select an optimal value of for each temporal regularization model by selecting the value of which produces the reconstruction with the smallest root mean square error (RMSE). The RMSE values were calculated separately for three regions; tumor, vascular region and the rest of the image domain. The RMSEs of different ROIs were then used to define a joint RMSE as
where corresponds to pixels in the vascular region, correspond to pixels in the tumour region and corresponds to pixels in rest of the image domain. The RMSE was calculated this way to weigh the small tumour and vascular regions appropriately. In estimating the pharmacokinetic parameters of tissues, obtaining an accurate arterial input function (AIF) is required Tofts1999. The AIF can be obtained via population averaging, however, usage of patient specific AIF produces more accurate estimates of the kinetic parameters Port2001. The AIF is preferably extracted from an artery feeding the tissues of interest, but it can also be estimated from a venous sinus or vein when an artery is not visible Lavini2010. Here the AIF is estimated from the superior sagittal sinus.
3 Materials and methods
A simulated test case modelling DCE-MRI measurements of a glioma in rat brain was created. The rat brain phantom was based on the rat brain atlas in Valdes2011, and scaled to a size of 128x128. The rat brain image was divided into three subdomains of different signal behaviour: vascular region (labelled ’1’ and highlighted with red in figure 1), tumour region (labelled ’2’ and highlighted with blue in figure 1) and the rest of the brain tissue. The vascular signal region corresponds to the location of the superior sagittal sinus.
A time series of 2800 ground truth images was simulated by multiplying the signal of each pixel with the template of the corresponding region and adding that to the baseline value of the pixels. The tumour signal templates were based on an experimental DCE-MRI measurement described in Hanhela2019, where the three different ROIs were identified. Figure 1 shows the signal templates for each of the different tissue regions.
One spoke of GA k-space data was simulated for each of the simulated images, leading to a dynamic experiment with 2800 spokes of k-space data. The time scale of the simulation was set to be similar to the in vivo measurements in Hanhela2019 where the measurement time between consecutive GA spokes was 38.5ms. Gaussian complex noise at 5% of the mean of the absolute values of the signal was added to the simulated k-space signal. The simulated test case was carried out using a k-space trajectory which combines the golden angle and the concentric squares sampling strategies.
In Hanhela2019 it was found that reconstruction of the form (5) performed optimally with segment length111Segment length equals the number of radial spokes per image. The number of elements in the data vector is segment length times number of samples per spoke. of 34 for a similar data set, thus we selected 34 as the segment length for our reconstructions leading to a temporal resolution of s.
All the regularized reconstructions in this work were computed using the Chambolle-Pock primal-dual algorithm Chambolle2011. In the NUFFT implementation of the forward model, the measurements were interpolated into a twice oversampled cartesian grid with min-max Kaiser-Bessel interpolation with a neighbourhood size of 4 FesslerSutton. The regridding reconstructions were computed using a Matlab code developped at Imagerie par résonance magnétique médicale et multi-modalités (IR4M) UMR8081, Université Paris-Sud, France.
We remark that when computing the RMS error (10), the reconstructed time signals of each pixel were linearly interpolated in the temporal direction to match the temporal resolution of the ground truth phantom.
The selection of was carried out using the first 34 spokes (i.e. the first frame ) of the DCE-MRI data and then the selected was used for all the spatio-temporal reconstructions with different temporal regularizations. The rat brain phantom of section 3 was used to compute the a priori level of sparsity, i.e. in our case we computed the TV norm of the first time frame of the dynamic phantom. This resulted in a sparsity level of . The spatial regularization term was selected using the S-curve as described in section 2.3 and resulted in spatial regularization parameter value of . The TV norms of the reconstructions for the S-curve were computed with 15 values of alpha ranging on interval . These 15 values of TV norm were then interpolated using Matlab’s interp function to 405 values. The resulting S-curve for the determination of is presented in figure 2.
In many practical applications, the a priori information, which we use to estimate the value of , may come from a different modality or from acquisition with different pulse sequence than the one used in the dynamic measurement. Therefore, in order to compute meaningful estimate of for the TV-regularized case, the reference image has to be scaled such that it is compatible with the measured dynamic data. This normalization of the reference image can be obtained by
where denotes the reference image, the frame of dynamic data that is used in the S-curve and the respective forward model.
This selection resulted in smallest joint RMSE of 0.085 for TS model which corresponded to , smallest joint RMSE of 0.058 for the TV model corresponding to and smallest joint RMSE of 0.063, corresponding to for the TGV model. As a reference method we selected the regularization parameters and using L-surface method for the case where we use spatial total variation and temporal total variation as our regularization model. The resulted L-surface is presented in figure 3. Application of L-surface on our data resulted parameters and .
Figure 4 shows slices of all reconstructions before, during and after contrast agent administration. Figure 5 shows the reconstructed images with different methods for one time frame. The top row of figure 5 shows the phantom with a red square denoting a domain that is presented as a closeup in figure 6 for all of the reconstructions.
As can be seen from figures 4 - 6, the classical regridding method fails on such high time resolution data as employed here, and thus we leave out the regridded reconstruction from the figures of the temporal evolution of the ROI signals. However the L-surface method seems to work nicely on our data, so we include it in the temporal evolution studies. The temporal evolutions in the vascular domain and in the tumor region are averages of and , respectively.
Figure 7 shows the time signals of the reconstructions in the vascular region () with the different temporal regularization models. Corresponding signals of the reconstructions in the tumor region (i.e. in ) are shown in figure 8. The tumor region is accurately reconstructed with all the methods, with only small differences between the methods, L-surface method being the noisiest. In the vascular region, the methods have some differences with TGV having the best reconstruction quality and TS having the worst reconstruction quality. The TS method shows smoothing at both the maximum and minimum signal levels whereas TGV reconstructs the fast signal change of the vascular region most reliably. L-surface method is again more noisy than the other reconstructions in the vascular region and its performance of reconstructing the vascular signal falls between TS and TV most likely due to the temporal regularization parameter being too small whereas the spatial regularization parameter is on the same order of magnitude as the parameter obtained with S-curve.
Variational regularization based solutions for dynamic MRI problems usually include two regularization parameters, one for the spatial and one for the temporal regularization, that the user has to select. Typically, the selection of both of the parameters is carried out manually based on visual assessment of the reconstructed images. In this work we proposed to use the S-curve method for the automatic selection of the spatial regularization parameter, leaving the time regularization parameter the only free parameter. The S-curve method selects the regularization parameter based on the expected sparsity of unknown images in domain of the regularization functional. Furthermore, the method requires computation of the reconstructions with relatively few values of the parameter, making it computationally efficient. The approach was demonstrated to lead to a feasible choice of the spatial regularization parameter in a simulated DCE-MRI experiment of rat brain. The reconstructions were also computed with three different temporal regularization models with the same fixed spatial regularization parameter, demonstrating the robustness of the approach with different time regularization models.
While we proposed automatic selection of the spatial regularization parameter, the temporal regularization parameter was still selected manually. In this work, we selected the temporal regularization parameter by computing RMS errors with respect to a ground truth. If on the other hand the ground truth is unknown, the temporal regularization parameter could be selected based on visual assessment of reconstructed images, leaving to be the only free parameter. In the future work, we aim to study methods for automatic selection of the time regularization parameter as well. One possibility could be to extend the S-curve to select a parameter which leads to an expected sparsity level in the domain of the temporal regularization. A feasible estimate for the expected level of sparsity in the time direction could potentially be extracted from the changes in the dynamic measurement data.
This work was supported by Jane and Aatos Erkko foundation and the Academy of Finland, Centre of Excellence in Inverse Modelling and imaging (project 312343).