Data-Driven Regularization Parameter Selection in Dynamic MRI

Abstract

In dynamic MRI, sufficient temporal resolution can often only be obtained using imaging protocols which produce undersampled data for each image in the time series. This has led to the popularity of compressed sensing (CS) based reconstructions. One problem in CS approaches is determining the regularization parameters, which control the balance between data fidelity and regularization. We propose a data-driven approach for the total variation regularization parameter selection, where reconstructions yield expected sparsity levels in the regularization domains. The expected sparsity levels are obtained from the measurement data for temporal regularization and from a reference image for spatial regularization. Two formulations are proposed. Simultaneous search for a parameter pair yielding expected sparsity in both domains (S-surface), and a sequential parameter selection using the S-curve method (Sequential S-curve). The approaches are evaluated using simulated and experimental DCE-MRI. In the simulated test case, both methods produce a parameter pair and reconstruction that is close to the root mean square error (RMSE) optimal pair and reconstruction. In the experimental test case, the methods produce almost equal parameter selection, and the reconstructions are of high perceived quality. Both methods lead to a highly feasible selection of the regularization parameters in both test cases while the sequential method is computationally more efficient.

1. Introduction

Dynamic magnetic resonance imaging is used to study hemodynamics, microvascular structure and function by contrast agent or stimulus-related changes in a time series of MR images. The high spatial and temporal resolution of the dynamic image series is often required for accurate analysis of the contrast agent or stimulus dynamics. In many cases, sufficient time resolution can only be obtained by utilizing an imaging protocol that produces undersampled data for each image in the time series. This, however, has the complication that reconstructing undersampled datasets with conventional MR image reconstruction methods, such as re-gridding [1], lead to noisy image series with poor spatial resolution.

Recently, the compressed sensing (CS) framework has led to significant advances in MRI with undersampled data. The theory of CS states that a target signal or image that is sparse in some basis, which is also incoherent with the measurement basis, can be perfectly reconstructed from undersampled data with a high probability [2,3,4]. Compressed sensing-based approaches have been developed for numerous applications in both static and dynamic MRI, see for example the review [5].

Provided that the temporal resolution of the MR image series is high enough, one can expect high redundancy in the image series in the sense that the image intensity changes between successive image frames are small and occur only in some parts of the image domain. Therefore, a compressed sensing approach based on the sparsity of the time derivative of the images is warranted.

The basic structure in the CS approaches to dynamic MRI is to reconstruct the whole time series of images simultaneously using an appropriate joint reconstruction formulation where a temporal regularization functional is employed for coupling the undersampled data across the time series of images. The most popular approach has been to use total variation (TV) regularization to promote sparsity of the spatial and temporal derivatives of the images, leading to a formulation [6,7]

$$\widehat{u}=\underset{u}{arg\; min}\left\{D(u,m)+\alpha {\mathrm{TV}}_{\mathrm{S}}\left(u\right)+\beta {\mathrm{TV}}_{\mathrm{T}}\left(u\right)\right\}$$

(1)

where $u=\{u^1,u^2,\dots ,u^T\}$ denotes the time series of unknown images, $m=\{m^1,m^2,\dots ,m^T\}$ is the MRI data time series, T is the number of time frames, D is the data fidelity term, ${\mathrm{TV}}_{\mathrm{S}}$ is the spatial total variation regularization functional, ${\mathrm{TV}}_{\mathrm{T}}$ is the temporal TV regularization functional and $\alpha$ and $\beta$ are the spatial and temporal regularization parameters, respectively.

The selection of the regularization parameters is crucial in terms of resulting image quality. Classical parameter selection methods, such as the Morozov discrepancy principle [8], exist, but these have been designed for the selection of a single regularization parameter and derived for a restricted class of problems. Alternatively, the parameter selection can be done by the L-curve method [9], Monte-Carlo SURE [10,11,12] or by subjective visual assessment of the reconstructions.

In this work, we propose a data-driven approach for the selection of the regularization parameters $\alpha$ and $\beta$ in (1). The proposed idea is to select a parameter pair that produces an a priori expected level of sparsity in both the domain of the spatial TV regularization and the temporal TV regularization.

As a part of the dynamic MRI measurement protocol, fully sampled cartesian measurements of the target are typically taken before the dynamic measurements. The anatomical MRI reconstruction is often used only for visualization, but it can also be used to obtain an a priori estimate for the spatial sparsity of the target. In case an anatomical reference image is not acquired, the expected spatial sparsity could also be estimated from a static reconstruction of a long sequence of the baseline of the dynamic data or from an anatomical atlas image.

An a priori estimate for the expected sparsity in the temporal regularization domain can be extracted from the dynamic MRI data. More specifically, the DC component of the k-space data gives information about the integral intensity of the unknown image at each frame. This information can then be used to approximate the expected sparsity level of the time derivative of the images when the image changes are based on contrast changes instead of tissue movement. The approach leads to a two-dimensional search from a set of reconstructions over a grid of parameter values $(\alpha ,\beta )$.

The idea in the proposed method is a 2D extension of the S-curve method, which was originally proposed for the selection of regularization parameter in sparsity promoting regularization in [13] and later applied to parameter selection in static x-ray tomography in [14,15,16]. The S-curve method was also applied to spatial TV regularization parameter selection in DCE-MRI in [17].

Since the proposed 2D parameter search (S-surface) can require a large number of reconstructions over a grid of values of $(\alpha ,\beta )$, we also consider a second method (sequential S-curve) where the parameters $\alpha$ and $\beta$ are selected separately by applying the S-curve method sequentially to first the temporal regularization parameter $\beta$ and then the spatial regularization parameter $\alpha$.

The proposed parameter selection approaches are evaluated using simulated and experimental DCE-MRI data from a rat brain study. We compare the proposed approaches to the L-curve [9] and Monte-Carlo SURE [10,11,12] parameter selection methods and in the simulated case, also to the true target.

In DCE-MRI, a bolus of gadolinium-based contrast agent is injected into the blood circulation and a time series of MRI data with an appropriate $T_1$-weighting is measured 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. DCE-MRI is used, for example, in treatment monitoring of breast cancer [18,19] and glioma [20], and to detect perfusion abnormalities in brain diseases [21,22].

2. Theory

2.1. Forward Model

The discrete MRI measurement model for a single image with cartesian measurements is of the form

$$m=\mathcal{F}u+e,$$

(2)

where $m\in {\mathbb{C}}^M$ is a complex-valued measurement vector, where M is the number of k-space measurement points, $\mathcal{F}$ is the discrete Fourier transform, $u\in {\mathbb{C}}^N$ is a complex-valued (vectorized) image, where N is the number of pixels in the image, and e models measurement noise. In the case of a non-cartesian k-space sampling trajectory, the Fourier transform can be approximated with the non-uniform fast Fourier transform (NUFFT) operation [23].

When NUFFT is used as the forward model, the measurement model can be written as

$$m=P\mathcal{F}Su+e,$$

(3)

where P is an interpolation and sampling matrix from the cartesian k-space to the non-cartesian k-space trajectory and S is a scaling matrix. Hereafter we denote $A:=P\mathcal{F}S$.

When considering dynamic MRI with complementary k-space sampling, where different (undersampled) trajectories of the k-space are measured at different time points, the forward model changes depending on the time point. The forward model can then be written as

$$m^t=P^t\mathcal{F}{S}^t{u}^t+e^t=A^t{u}^t+e^t,$$

(4)

where the integer-valued superscript t denotes the time index of the measurement and image series, and $m^t$ is the vector of k-space data for a single image in the time series, which is dependent on the selected data segmentation length.

2.2. Joint Reconstruction Formulation of the Dynamic Inverse Problem

The CS based joint image reconstruction formulation we consider is

$$\widehat{u}=\underset{u=u^1,u^2,\dots ,u^T}{arg\; min}\left\{\sum _{t=1}^T\left[\parallel {A^t{u}^t-m^t\parallel }_2^2+\alpha \parallel {{\nabla }_{\mathrm{S}}{u}^t\parallel }_1\right]+\beta \parallel {{\nabla }_{\mathrm{T}}u\parallel }_1\right\},$$

(5)

where T is the number of image frames in the problem, ${\nabla }_{\mathrm{S}}$ is the discrete spatial gradient operator, ${\nabla }_{\mathrm{T}}$ is the discrete temporal gradient operator, and $\alpha$ and $\beta$ are regularization parameters for spatial and temporal regularizations respectively.

Here, we use the isotropic form of 2D spatial TV, where the total variation functional for a complex valued vectorized image $u^t$ is defined as

$$\parallel {{\nabla }_{\mathrm{S}}{u}^t\parallel }_1=\sum _{k=1}^N{\left({(Re\left({\mathrm{D}}_{\mathrm{x}}^k{u}^t\right))}^2+{(Re\left({\mathrm{D}}_{\mathrm{y}}^k{u}^t\right))}^2+{(Im\left({\mathrm{D}}_{\mathrm{x}}^k{u}^t\right))}^2+{(Im\left({\mathrm{D}}_{\mathrm{y}}^k{u}^t\right))}^2\right)}^{1/2},$$

(6)

where Re and Im denote taking the real and the imaginary parts of the complex valued image, k denotes the spatial index in the 2D images, and $D_x^k$ and $D_y^k$ are discrete forward differences in the horizontal and vertical image directions of the k’th pixel defined as

$${\mathrm{D}}_{\mathrm{x}}^k{u}^t=\{\begin{array}{cc}\begin{array}{c}-u_k^t+u_{k+n}^t,\hfill \end{array}\hfill & \mathrm{if} k\le N-n\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}$$

(7)

$${\mathrm{D}}_{\mathrm{y}}^k{u}^t=\{\begin{array}{cc}\begin{array}{c}-u_k^t+u_{k+1}^t,\hfill \end{array}\hfill & \mathrm{if} k\left(\mathrm{mod}n\right)\ne 0\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}$$

(8)

where n is the number of rows and columns in the image that is assumed to be square.

Similarly, the temporal TV is defined by

$$\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1=\sum _{t=1}^T\sum _{k=1}^N\sqrt{{(Re\left({\mathrm{D}}_{\mathrm{T}}^t{u}_k\right))}^2+{(Im\left({\mathrm{D}}_{\mathrm{T}}^t{u}_k\right))}^2},$$

(9)

where $u_k=u_k^1,\dots ,u_k^T$ and ${\mathrm{D}}_{\mathrm{T}}^t$ is the discrete forward difference in the temporal direction of the t’th image defined as

$${\mathrm{D}}_{\mathrm{T}}^t{u}_k=\{\begin{array}{cc}-u_k^t+u_k^{t+1},\hfill & \mathrm{if} t\ne T\hfill \\ 0,\hfill & \mathrm{otherwise}.\end{array}$$

(10)

2.3. Proposed Parameter Selection Methods

The S-curve was originally proposed for the selection of a regularization parameter in wavelet-based sparsity promoting regularization in [13]. The idea in the S-curve method is to select the regularization parameter such that the reconstructed image has an a priori expected sparsity level in the regularization domain.

The CS formulation of the dynamic MRI problem in (5) includes two regularization parameters, $\alpha$ and $\beta$, making the parameter selection a two-dimensional problem. Here, the idea of the S-curve method is extended to the selection of two regularization parameters, which we propose to select by finding the parameter pair $(\widehat{\alpha },\widehat{\beta })$ that produces the expected sparsity level in both the spatial and temporal regularization domains, with a simultaneous or sequential selection of the parameters.

2.3.1. Selection of the a Priori Sparsity Estimates

Assume that we have an a priori estimate

$${\widehat{S}}_{\mathrm{S}}={\mathrm{TV}}_{\mathrm{S}}\left(u^{\mathrm{ref}}\right):=\parallel {{\nabla }_{\mathrm{S}}{u}^{\mathrm{ref}}\parallel }_1,$$

(11)

for the spatial total variation norm of a single unknown image $u^t$, obtained from a reference image $u^{\mathrm{ref}}$. The reference image can be for example an anatomical image from fully sampled measurements, a conventional reconstruction of a long sequence of the baseline of the dynamic data or an anatomical atlas image.

Remark that in cases where the expected sparsity level ${\widehat{S}}_{\mathrm{S}}$ is estimated from a reference image $u^{\mathrm{ref}}$ that has been acquired with different contrast parameters or comes from a different imaging modality, the reference image may need to be normalized to the scale of the dynamic images. This normalization of the reference image can be done, for example, by

$$u^{\mathrm{ref}}\leftarrow \frac{\parallel {m^1\parallel }_2}{\parallel {A^1{u}^{\mathrm{ref}}\parallel }_2}{u}^{\mathrm{ref}},$$

(12)

where $m^1$ is the first baseline frame of dynamic MRI data and $A^1$ is the respective forward model.

When the temporal changes in the image are based on (mostly) unidirectional contrast changes, an estimate for the sparsity level ${\widehat{S}}_{\mathrm{T}}$ of the temporal gradient can be obtained from the dynamic data using the zero frequency (DC) components of the k-space data. The sparsity estimate is based on the property that the DC component of the Fourier transform of a function f equals the total intensity of the image, i.e.,

$$\widehat{f}\left(0\right)={\int }_{\Omega }f\left(r\right)\mathrm{d}r.$$

(13)

Therefore, for the total intensity difference between two images $(w,z)$ we have

$$|{\int }_{\Omega }w\left(r\right)\mathrm{d}r-{\int }_{\Omega }z\left(r\right)\mathrm{d}r=|\widehat{w}\left(0\right)-\widehat{z}\left(0\right)|.$$

(14)

Thus, an estimate for the temporal total variation

$${\mathrm{TV}}_{\mathrm{T}}\left(u\right)=\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1$$

(15)

of the unknown image sequence u can be obtained from the measurement data by

$${\widehat{S}}_{\mathrm{T}}\approx \sum _{t=1}^{T-1}|m_{(k=0)}^{t+1}-m_{(k=0)}^t|,$$

(16)

where the subscript $(k=0)$ refers to one of the zero frequency components of the k-space data vector $m^t$.

Note that under ideal noiseless measurement data, the estimate ${\widehat{S}}_{\mathrm{T}}$ would be a lower bound for the temporal TV of the image sequence. The difference between the zero-frequency k-space coefficients of two images would be equal to the temporal TV of the images when there is no noise or tissue movement and the contrast changes between the two images are completely unidirectional. When the contrast changes between the two images are of different signs in different parts of the target or the image changes are caused by motion, the difference of the zero-frequency k-space coefficients would underestimate the temporal TV between the images. In presence of measurement noise or forward model errors, the estimate (16) may also be larger than the true temporal TV of the unknown image sequence.

2.3.2. Simultaneous Selection of the Parameters (S-Surface)

Given the estimates ${\widehat{S}}_{\mathrm{T}}$ and ${\widehat{S}}_{\mathrm{S}}$ of the a priori expected temporal and spatial sparsity, we select the regularization parameters $(\alpha ,\beta )$ as follows:

(1)

Take a grid of regularization parameters $\alpha$ and $\beta$ ranging on the interval $(0,\infty )$ such that

$$0<{\beta }^{\left(1\right)}<{\beta }^{\left(2\right)}<\cdots <{\beta }^{\left(P\right)}<\infty$$

and

$$0<{\alpha }^{\left(1\right)}<{\alpha }^{\left(2\right)}<\cdots <{\alpha }^{\left(L\right)}<\infty .$$

(2)

Compute the corresponding estimates $\widehat{u}({\alpha }^{(\ell )},{\beta }^{\left(p\right)}),\ell =1,\dots ,L,p=1,\dots ,P$. The reconstructions $\widehat{u}({\alpha }^{(\ell )},{\beta }^{\left(p\right)})$ are computed by

$$\widehat{u}=\underset{u=u^1,u^2,\dots ,u^T}{arg\; min}\left\{\sum _{t=1}^T\left[\parallel {A^t{u}^t-m^t\parallel }_2^2+{\alpha }^{(\ell )}\parallel {{\nabla }_{\mathrm{S}}{u}^t\parallel }_1\right]+{\beta }^{\left(p\right)}\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1\right\}.$$

(17)

Here, it is crucial that ${\beta }^{\left(1\right)}$ is selected so small that the reconstructions have a larger temporal TV than the expected temporal sparsity level, and ${\alpha }^{\left(1\right)}$ is so small that the reconstructions have a larger spatial TV than the expected spatial sparsity level. Similarly, the values ${\beta }^{\left(P\right)}$ and ${\alpha }^{\left(L\right)}$ are selected so large that the problem becomes over-regularized and the respective TV values of the reconstructions are small.

(3)

Compute the temporal TV of the image series ${\mathrm{TV}}_{\mathrm{T}}\left(\widehat{u}(\alpha ,\beta )\right)$ and the spatial TV of one time frame of the image series ${\mathrm{TV}}_{\mathrm{S}}\left({\widehat{u}}^t(\alpha ,\beta )\right)$ for each grid point $({\alpha }^{(\ell )},{\beta }^{\left(p\right)})$.

(4)

Evaluate the value of the merit functional

$$\Psi (\alpha ,\beta )=\frac{|{\mathrm{TV}}_{\mathrm{T}}\left(\widehat{u}(\alpha ,\beta )\right)-{\widehat{S}}_{\mathrm{T}}|}{2{\widehat{S}}_{\mathrm{T}}}+\frac{|{\mathrm{TV}}_{\mathrm{S}}\left({\widehat{u}}^t(\alpha ,\beta )\right)-{\widehat{S}}_{\mathrm{S}}|}{2{\widehat{S}}_{\mathrm{S}}}$$

(18)

for each grid point.

(5)

Select the parameter pair $(\widehat{\alpha },\widehat{\beta })$ which minimizes $\Psi (\alpha ,\beta )$.

2.3.3. Sequential Selection of the Parameters (Sequential S-Curve)

The two-dimensional selection of the parameters $(\alpha ,\beta )$ in Section 2.3.2 requires computing a large number of reconstructions of the form of (5) over a 2D grid of regularization parameter values, making it computationally expensive, especially in large 4D problems. Therefore, to alleviate the computational burden, we also consider an approach where we split the parameter search into two 1D problems.

We employ the S-curve first for the selection of the temporal regularization parameter $\beta$ without any spatial regularization (i.e., $\alpha =0$ in (5)), and then we employ the S-curve for the selection of the spatial regularization parameter $\alpha$ using the value of $\beta$ found in the first step. We select the parameter $\beta$ first since the reconstruction from undersampled measurements is more dependent on the temporal regularization than the spatial regularization.

Selection of $\mathit{\beta }$: Given the estimate ${\widehat{S}}_{\mathrm{T}}$ of the a priori temporal sparsity level, we first select the temporal regularization parameter $\beta$ using the S-curve as follows:

(1)

Take a sequence of the temporal regularization parameters $\beta$ ranging on the interval $(0,\infty )$ such that

$$0<{\beta }^{\left(1\right)}<{\beta }^{\left(2\right)}<\cdots <{\beta }^{\left(P\right)}<\infty .$$

(2)

Compute the corresponding estimates $\widehat{u}\left({\beta }^{\left(1\right)}\right),\dots ,\widehat{u}\left({\beta }^{\left(P\right)}\right)$ according to

$$\widehat{u}\left({\beta }^{\left(p\right)}\right)=\underset{u=u^1,u^2,\dots ,u^T}{arg\; min}\left\{\sum _{t=1}^T\parallel {A^t{u}^t-m^t\parallel }^2+{\beta }^{\left(p\right)}\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1\right\}.$$

(19)

Here too, ${\beta }^{\left(1\right)}$ has to be so small that the problem is under-regularized and the corresponding reconstruction $\widehat{u}\left({\beta }^{\left(1\right)}\right)$ is a noisy image sequence with a temporal TV larger than ${\widehat{S}}_{\mathrm{T}}$ and ${\beta }^{\left(P\right)}$ is so large that the problem is over-regularized and the temporal TV of the reconstruction is close to zero.

(3)

Compute the temporal TV of the recovered estimates $\widehat{u}\left({\beta }^{\left(p\right)}\right),p=1,\dots ,P$.

(4)

Fit a smooth interpolation curve to the data $\{{\beta }^{\left(p\right)},{\mathrm{TV}}_{\mathrm{T}}\left(\widehat{u}\left({\beta }^{\left(p\right)}\right)\right),p=1,\dots ,P\}$ and use the interpolated sparsity curve to find the value $\widehat{\beta }$ for which ${\mathrm{TV}}_{\mathrm{T}}\left(\widehat{\beta }\right)={\widehat{S}}_{\mathrm{T}}$.

Selection of $\mathit{\alpha }$: Given the value $\widehat{\beta }$ and the a priori expected spatial sparsity level ${\widehat{S}}_{\mathrm{S}}$, the spatial regularization parameter $\alpha$ is selected according to the following procedure:

(1)

Take a sequence of the spatial regularization parameters $\alpha$ ranging on the interval $(0,\infty )$ such that

$$0<{\alpha }^{\left(1\right)}<{\alpha }^{\left(2\right)}<\cdots <{\alpha }^{\left(L\right)}<\infty .$$

(2)

Compute the corresponding estimates $\widehat{u}\left({\alpha }^{\left(1\right)}\right),\dots ,\widehat{u}\left({\alpha }^{\left(L\right)}\right)$ by

$$\widehat{u}({\alpha }^{(\ell )},\widehat{\beta })=\underset{u=u^1,u^2,\dots ,u^T}{arg\; min}\left\{\sum _{t=1}^T\left[\parallel {A^t{u}^t-m^t\parallel }_2^2+{\alpha }^{(\ell )}\parallel {{\nabla }_{\mathrm{S}}{u}^t\parallel }_1\right]+\widehat{\beta }\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1\right\}.$$

(20)

Again, ${\alpha }^{\left(1\right)}$ needs to be so small that a frame in the reconstruction $\widehat{u}\left({\alpha }^{\left(1\right)}\right)$ results in a spatial TV larger than ${\widehat{S}}_{\mathrm{S}}$ and ${\alpha }^{\left(L\right)}$ needs to be so large that the spatial TV of a frame $\widehat{u}\left({\alpha }^{\left(L\right)}\right)$ is very small.

(3)

Compute the spatial TV of a frame of the recovered estimates $\widehat{u}({\alpha }^{(\ell )},\widehat{\beta }),\ell =1,\dots ,L$, i.e., ${\mathrm{TV}}_{\mathrm{S}}\left({\widehat{u}}^t({\alpha }^{\ell },\widehat{\beta })\right)$.

(4)

Fit a smooth interpolation curve to the data $\{{\alpha }^{(\ell )},{\mathrm{TV}}_{\mathrm{S}}\left({\widehat{u}}^t({\alpha }^{(\ell )},\widehat{\beta })\right),\ell =1,\dots ,L\}$, and use the interpolated sparsity curve to find the value $\widehat{\alpha }$ for which ${\mathrm{TV}}_{\mathrm{S}}(\widehat{\alpha },\widehat{\beta })={\widehat{S}}_{\mathrm{S}}$.

The final reconstruction for analysis and diagnostic tasks can then be computed using the parameter pair $(\widehat{\alpha },\widehat{\beta })$.

2.4. Comparison Methods

2.4.1. L-Curve Parameter Selection

The L-curve parameter selection [9] was computed as a reference for the proposed sparsity-based methods. Similarly to the Sequential S-curve method, the L-curve was implemented sequentially such that the temporal regularization parameter $\beta$ was selected first, and then the spatial regularization parameter $\alpha$ was selected using a fixed $\beta$.

Selection of regularization parameters: To compute the ensemble of reconstructions with different $\beta$:s or $\alpha$:s, steps (1) and (2) of the respective S-curve parameter selection in Section 2.3.3 were first computed. After this the L-curve method was used as follows:

(3)

Set $\rho ={log}_{10}{\sum }_{t=1}^T\parallel {A^t{u}^t-m^t\parallel }_2^2$, $\eta ={log}_{10}\parallel R\left(u\right)\parallel$, where, $R\left(u\right)$ is either the sum of the spatial total variations of all the frames ${\sum }_{t=1}^T\parallel {{\nabla }_{\mathrm{S}}{u}^t\parallel }_1$ or the temporal total variation $\parallel {{\nabla }_{\mathrm{T}}u\parallel }_1$.

(4)

Interpolate the L-curve ($\rho$,$\eta$) to a dense grid using spline interpolation.

(5)

Calculate the maximum curvature point of the L-curve using

$$\lambda =argmax\kappa \left(\lambda \right)=argmax\frac{{\rho }^{\prime }{\eta }^{″}-{\rho }^{″}{\eta }^{\prime }}{{({\left({\rho }^{\prime }\right)}^2+{\left({\eta }^{\prime }\right)}^2)}^{3/2}},$$

(21)

where $\lambda$ is either $\alpha$ or $\beta$.

2.4.2. Monte-Carlo SURE Parameter Selection

Monte-Carlo SURE (MC-SURE) parameter selection [10,11,12] was also computed as a reference. MC-SURE estimates a weighted mean squared error using the principles of Stein’s unbiased risk estimate [24] by perturbing the data vector with a complex-valued random vector and analyzing the response to the perturbation. Similar to both the Sequential S-curve method and the L-curve method, the MC-SURE parameter selection was first done for the temporal regularization parameter $\beta$ and then the spatial regularization parameter $\alpha$ with a fixed $\beta$.

Selection of regularization parameters: The full MC-SURE function is of the form

$$\mathrm{MC}-\mathrm{SURE}\left(\lambda \right)=M^{-1}\parallel {A{u}_{\lambda }\left(m\right)-m\parallel }_W^2-M^{-1}\mathrm{tr}\{\Omega W\}+2M^{-1}Re\left\{\mathrm{tr}\{\Gamma J_{u_{\lambda }\left(m\right)}\}\right\},$$

(22)

where $\lambda$ is the parameter pair $(\alpha ,\beta )$, $u_{\lambda }\left(m\right)$ is the reconstruction computed with the regularization parameter pair $\lambda$ and measurements m, W is a weighting matrix, $\Omega$ is the noise covariance matrix, $\Gamma =\Omega WA$ and M is the number of measurement points. Further, the second trace is approximated by

$$\mathrm{tr}\{\Gamma J_{u_{\lambda }\left(m\right)}\}\approx {ϵ}^{-1}{b}^{\prime }{\Lambda }^{-1}\Gamma \rho (u_{\lambda },m,\Lambda b,ϵ),$$

(23)

where $ϵ$ is a data perturbation multiplier, $\Lambda$ is a matrix allowing different amounts of perturbation for different elements of m, b is a complex valued random vector and $\rho$ is the perturbation between the reconstructions from the perturbed data and the original data

$$\rho (u_{\lambda },m,\Lambda b,ϵ)=u_{\lambda }(m+ϵ\Lambda b)-u_{\lambda }\left(m\right).$$

(24)

Here, as in [11], we set $b=(b_{Re}+i{b}_{Im})/\sqrt{2}$, where $b_{Re}$ and $b_{Im}$ are independent binary random vectors of −1 and 1 with equal probability. In addition, we set $\Lambda =I$ and $W=I$. Since we are minimizing (22) w.r.t. $\lambda =(0,\beta )$ or $\lambda =(\alpha ,\widehat{\beta })$, where $\widehat{\beta }$ is fixed, we can in addition leave out the second term, which is the same for all $\lambda$, and the multiplier $M^{-1}$. The approximated and shortened form of the MC-SURE function is thus

$$\mathrm{MC}-\mathrm{SURE}\left(\lambda \right)\approx \parallel {A{u}_{\lambda }\left(m\right)-m\parallel }^2+2Re\{{ϵ}^{-1}{b}^{\prime }\Omega A(u_{\lambda }(m+ϵb)-u_{\lambda }\left(m\right))\}.$$

(25)

In practice, MC-SURE is implemented as follows:

(1)

Compute the reconstructions (19) with the original data and the perturbed data for parameter pairs $\lambda =(\alpha =0,{\beta }^{\left(p\right)})$ with $p=1,\dots ,P$.

(2)

Calculate MC-SURE according to (25) and choose the parameter $\widehat{\beta }$ corresponding to the minimum MC-SURE.

(3)

Compute the reconstructions (20) with the original data and the perturbed data for parameter pairs $\lambda =({\alpha }^{(\ell )},\widehat{\beta })$ with $\ell =1,\dots ,L$.

(4)

Calculate MC-SURE according to (25) and choose the parameter $\widehat{\alpha }$ corresponding to the minimum MC-SURE.

The parameter selection by MC-SURE thus requires computing a total of $2(P+L)$ reconstructions.

3. Test Case Setting

3.1. Simulated Test Case

A simulation modeling DCE-MRI of a glioma in a rat brain was generated. The rat brain image used as the basis for the simulation was the rat brain atlas image [25], and the image was scaled to be of size 128 × 128. The rat brain tissue area was split into three subdomains with varying signal behavior: (1) vascular region corresponding to the location of the superior sagittal sinus (highlighted with blue and labeled ‘1’ in Figure 1), (2) generated tumor region (highlighted with red and labeled ‘2’ in Figure 1) and (3) the rest of the brain tissue. The superior sagittal sinus can be used as a proxy for estimating the arterial input function (AIF) in DCE-MRI of the brain [26]. The AIF is used in the estimation of pharmacokinetic parameters in DCE-MRI [27], and the vascular region is thus of special interest.

Figure 1 shows the signal templates generated for each of the three regions of interest. We created 2800 ground truth images based on the signal templates by multiplying the signal of each pixel in the original image with the signal template of the corresponding region and adding the result of the multiplication to the original value of the pixel. The three template signals were based on an experimental DCE-MRI measurement, which is described briefly in Section 3.2 and also in [28,29], from which the three different regions of interest (ROIs) were identified. The same simulated test case has been used previously in [17,28].

The simulated test case was carried out using a k-space trajectory which combines the golden angle (GA) [30] and the concentric squares sampling strategies [31,32].

For each of the 2800 ground truth images, one spoke of k-space data was generated. The repetition time of the experimental measurements was 38.5 ms, and the signal dynamics of the simulation were set to reflect this. Finally, complex Gaussian noise at 2%, 5% and 10% of the mean of the absolute values of the signal without noise was added to the simulated k-space signal to create data realizations with different noise levels.

Error Metric for the Simulation

Root mean square error (RMSE) values were calculated for the three regions of interest: (1) vascular region, (2) tumor region and (3) the rest of the image. For the calculation of the RMSE, the reconstructed signals of each pixel were linearly interpolated in the temporal direction to match with the temporal resolution of the ground truth phantom (i.e., segment length of one spoke). After the interpolation, RMSE to the ground truth phantom was calculated for all regions separately. After the RMSEs for the three ROIs were computed, a joint RMSE was computed by taking the norm of the separate RMSEs. This was done to weigh the ROIs with small support (vascular, tumor) appropriately in the error metric. The error metric was selected independent of the parameter selection, and the same metric was used in [17,28]. The error metric is described in more detail in [28].

3.2. In Vivo Measurements from a Rat Glioma Model

The animal experiment was approved by the Animal Health Welfare and Ethics Committee of the University of Eastern Finland. DCE-MRI data were acquired from a female Wistar rat with a glioma model [33]. The experiment was performed 10 days post-implantation of the glioma cells into the rat brain. During the imaging, the animal was anesthetized and kept in a fixed position in a holder that was inserted into the magnet. A needle was placed into the animal’s tail vein for the injection of the contrast agent.

The DCE data were collected using a 9.4 T horizontal magnet interfaced to an Agilent imaging console and a volume coil transmit/quadrature surface coil receive pair (Rapid Biomed, Rimpar, Germany). The DCE-MRI data were measured with radial golden angle sampling using a gradient-echo based radial pulse sequence with field-of view 32 mm × 32 mm, slice thickness 1.5 mm, repetition time 38.5 ms, echo time 9 ms, flip angle 30 degrees and 128 points in each spoke. 610 spokes were collected in sequential order, after which the next spoke would differ by less than 0.1 degrees from the first spoke, so the cycle of 610 spokes was repeated to simplify the sequence. This cycle was repeated for a total of 25 times, which leads to a sequence of 15,250 spokes of data for a measurement duration of nearly 10 min. The measurement time for a single cycle of 610 spokes was 610 × 38.5 ms = 23.46 s. One minute after the beginning of the dynamic sequence, Gadovist (1 mmol/kg) was injected i.v. over a period of 3 s. From the begining of the measurements, 7320 spokes of data were used for testing the proposed parameter selection approach in the experimental data results section.

Anatomical reference images were also scanned from the same slice before and after the dynamic imaging using a gradient-echo pulse sequence with similar parameters with the dynamic data sequence but using a full Cartesian sampling of 128 × 128 points of k-space data. Anatomical images reconstructed from the full Cartesian data with spatial regularization from before and after the experiment are shown as a reference to the dynamical reconstructions in Figure 2.

The dynamical experiment was also used as a basis for creating the three signal templates shown in Figure 1 for the simulated test case. For the simulation, three regions were identified from the experimental data reconstructions: vascular region (superior sagittal sinus), glioma region, and the rest of the brain tissue.

3.3. Computation

The Chambolle-Pock primal-dual algorithm [34,35] was used to solve the minimization problem of (5). The ratio of the primal and dual step sizes was varied according to the regularization coefficients such that the primal step size was smaller (and accordingly the dual step size was larger) for larger regularization parameters. Asymmetrical primal and dual step sizes in the algorithm have been shown to lead to faster convergence in some cases in both linear [36] and non-linear [37] problems, and this was observed here as well. The operator norms of the forward problem, and the spatial and temporal total variation terms were all scaled to 1.

4. Results

4.1. Simulated Data

For both the simulated and experimental test case, we use a segment length of 34, i.e., each image frame is computed from 34 spokes of GA data. This was found to be an optimal segment length for a similar dynamic simulation in [28].

For the simulation, we use the parameter pair and the reconstruction that yields the lowest joint RMSE as references of the best possible parameter pair and reconstruction. We refer to this parameter selection method as MinRMSE and compare the S-surface, Sequential S-curve, L-curve and MC-SURE methods to the MinRMSE. The regularization parameters selected with the different methods for the simulation with 5% noise are shown in

Table 1. The spatial ($\widehat{\alpha }$) and temporal ($\widehat{\beta }$) regularization parameters for the simulated test case with 5% noise level obtained with the S-surface, Sequential S-curve, L-curve and Monte-Carlo SURE (MC-SURE) methods and the MinRMSE parameters used as reference. The table shows that the parameters obtained with the proposed methods (S-surface and Sequential S-curve) and MC-SURE are close to the optimal (MinRMSE) whereas the parameters obtained with the L-curve-method are further from the optimal.

	MinRMSE	S-Surface	Seq. S-Curve	L-Curve	MC-SURE
$\widehat{\alpha }$	0.0001	0.00032	0.00019	8.8 × ${10}^{-5}$	0.0001
$\widehat{\beta }$	0.0032	0.0032	0.0046	2.7 × ${10}^{-5}$	0.001

.

The reference temporal sparsity for the S-surface and Sequential S-curve methods was obtained by taking the spoke closest to vertical from each segment of 34 spokes, and calculating the total intensity changes using the DC signal of those spokes. This selection was done due to the DC signal of the real measurements having some dependence on the angle of the measurement spoke. The reference spatial sparsity level used was the spatial total variation of the first frame of the ground truth phantom. Note here, that since the temporal regularization parameter was first selected, the selection of the spatial regularization parameter value is not as sensitive as the temporal regularization parameter selection.

For MC-SURE parameter selection, $ϵ$ was set to ${10}^{-3}$ and the noise covariance $\Omega$ was set to the scalar noise level ${\sigma }^2$, which was obtained from the sample variance of both ends of all the measurement spokes.

Figure 3, Figure 4 and Figure 5 show the results using the simulated rat brain DCE data. Figure 3 shows the joint sparsity contour used to select the S-surface parameter pair as well as the parameter selection curves for the Sequential S-curve, L-curve and MC-SURE methods in the simulated case with 5% noise.

Figure 4 shows the joint RMSE contours of the reconstructions calculated on a wide grid of the spatial and temporal regularization parameters for the simulations with noise levels 2%, 5% and 10%. Figure 4 also shows the joint RMSE values of the reconstructions with the parameters chosen with the three different methods. The joint RMS error values for the Sequential S-curve, S-surface and MC-SURE reconstructions are similar to each other with the noise levels 2% and 5% and close to those of the optimal reconstruction (MinRMSE), and the regularization parameters obtained with these methods are also close to the MinRMSE parameters. In the simulation with 10% noise level, the S-surface performs worse than Sequential S-curve and MC-SURE in the joint RMSE measure. The L-curve method performs worse than the other methods with all three noise levels.

Figure 5 shows single frames of the reconstructions of the data with 5% noise after contrast agent injection as well as single-pixel signals from the tumor and vascular areas as well as healthy tissue on the right side of the cortex and the RMSEs of the three regions. The MinRMSE, S-surface and Sequential S-curve reconstructions (with 5% noise level) in Figure 4 appear mostly similar, whereas the L-curve and MC-SURE reconstructions are slightly noisier. Figure 5 also shows the RMS errors of the three different regions. In the region-based RMS error measures, the Sequential S-curve and S-surface yield the best accuracy in the tumor and tissue signals, while MC-SURE has the best accuracy in the vascular region.

4.2. Experimental Data

Figure 6 and Figure 7 show the results for the experimental data. Figure 6 shows the joint sparsity grid of reconstructions with a large grid of spatial and temporal regularization parameters where the S-surface, Sequential S-curve, L-curve and MC-SURE parameter pairs are marked. The figure also shows all the parameter selection curves used to select both $\alpha$ and $\beta$ with the Sequential S-curve, L-curve and MC-SURE methods. The parameters obtained with the S-surface and Sequential S-curve methods are almost the same, namely, the Sequential S-curve parameters are $\alpha =0.00034$, $\beta =0.0097$, and the S-surface parameters are $\alpha =0.00032$, $\beta =0.01$. The parameters obtained with the L-curve method are $\alpha =0.00025$, $\beta =0.00024$ and the parameters obtained with the MC-SURE method are $\alpha =0.0001$, $\beta =0.001$.

For MC-SURE parameter selection, similar to the simulation, $ϵ$ was set to ${10}^{-3}$ and the noise covariance $\Omega$ was set to the scalar noise level ${\sigma }^2$, which was obtained from the sample variance of both ends of all the measurement spokes.

Figure 7 shows the cartesian reference reconstruction with spatial regularization from before contrast agent injection and single-frame images from the dynamic reconstructions using parameters obtained with the different methods from after the contrast agent injection. The reference image is the same as in Figure 2 with different contrast. Figure 7 also contains single-pixel signals from the tumor and vascular regions. The reconstructions are mostly visually similar, and the single-pixel signals of the Sequential S-curve and S-surface reconstructions show less variation than those of the L-curve and MC-SURE reconstructions which yield smaller values of the temporal regularization parameter $\beta$ than the sparsity-based selections.

5. Discussion

Compressed sensing-based models for dynamic MRI problems usually include regularization for sparsity in both the spatial and temporal domain. Therefore these models typically include two regularization parameters, one for the spatial and one for the temporal regularization, that the user has to select. Often, the selection of both of the parameters is carried out manually based on visual assessment of the reconstructed images.

In this work, we investigated a two-dimensional parameter selection problem and proposed two-parameter selection methods using a priori image sparsity estimates obtained from the k-space data for temporal sparsity and from a reference image for spatial sparsity. The evaluations of the methods were carried out using both simulated and experimental golden angle DCE-MRI data. The two-parameter selection methods proposed were the S-surface and Sequential S-curve methods. The proposed methods were compared with the L-curve and MC-SURE parameter selection methods. In the simulated test case, the parameter selection methods were also compared with a parameter pair that gives the smallest joint RMS reconstruction error with respect to the true target images.

Both proposed parameter estimation methods produced a parameter pair that is close to the parameter pair with the minimum joint RMSE in the simulated test cases, and the RMS reconstruction errors with these parameter selections were only slightly larger than the minimum RMSE, especially with the Sequential S-curve method. With 10% noise, the S-surface and Sequential S-curve methods both yielded a bit too high spatial regularization parameter leading to a larger than optimal joint RMS error. In the simulated test cases, the L-curve yielded the worst parameter selections with all the different noise levels. The MC-SURE parameter selection performance in the joint RMSE measure in the simulated test case was slightly worse than with the Sequential S-curve, but slightly better than the S-surface.

In the experimental case, there is no minimum RMSE solution or other known optimal reference solution available, but the parameter pairs obtained with the proposed methods are close to a manually selected parameter pair. All the reconstructions are visually similar with the L-curve and MC-SURE reconstructions being slightly noisier, and having more signal fluctuation.

While the S-surface method is based on a 2D search which requires reconstructions over a grid of $P\times L$ parameter pairs $(\alpha ,\beta )$, the Sequential S-curve method is based on two 1D searches, and therefore it needs only $P+L+1$ reconstructions making it computationally more efficient than the S-surface method. The Sequential S-curve method is also computationally more efficient than the MC-SURE method, which requires the computation of two reconstructions for each parameter pair, i.e., requiring $2(P+L)$ reconstructions in total. This suggests that the Sequential S-curve selection could be a more favorable choice, especially in large-scale 4D problems, as the methods produced similar parameter selection and reconstruction accuracy.

As the expected temporal sparsity level ${\widehat{S}}_{\mathrm{T}}$ for the selection of the temporal regularization parameter is estimated from differences of the zero-frequency components of the measured k-space data between consecutive time steps, the estimate provides an approximation for the temporal total variation of the unknown image series. In the simulated test case with 5% noise, the discrepancy between the true temporal TV of the ground truth images and the estimate ${\widehat{S}}_{\mathrm{T}}$ obtained from the noisy k-space data using (16) was approximately $22\%$ due to the noise level in the simulated k-space data. Despite this discrepancy, the proposed parameter selections yielded highly feasible regularization parameters.

However, we remark that in cases where the images between consecutive time steps exhibit changes due to motion or large opposite signed contrast changes that would lead to only small changes in the zero frequency coefficients, the data based sparsity estimate could be clearly smaller than the actual temporal TV, potentially leading to a less optimal parameter choice. The widely used reconstruction with temporal TV regularization (19) and also the proposed parameter selection approaches are based on the assumption that the target does not move during the experiment. In preclinical small animal studies, this is usually a feasible assumption as the animal specimen is typically in anesthesia and kept in a fixed location. However, if motion should occur, for example in clinical brain imaging, then the measured k-space data could be corrected for rigid motion retrospectively before the parameter selection and reconstruction using e.g., [38].

6. Conclusions

In this work, we proposed to use two sparsity-based methods based on a priori sparsity estimates for the automatic selection of the regularization parameters for the TV regularized CS problem. The S-surface method selects the regularization parameters based on the expected sparsity of unknown images in the two domains of regularization simultaneously. We also adopted a faster method called the Sequential S-curve, where the regularization parameters were selected one at a time, which was shown to produce similar reconstruction accuracy and parameter selection. The Sequential S-curve method requires the computation of fewer reconstructions making it computationally more efficient.

Both of the approaches were demonstrated to lead to a highly feasible choice of the temporal and spatial regularization parameters in both the simulated and the experimental DCE-MRI experiment of a rat brain.