Analytic heuristics for a fast DSC-MRI

TL;DR

This paper introduces a deterministic heuristic method for fast DSC-MRI reconstruction, offering a computationally efficient alternative to existing compressed sensing techniques, with demonstrated effectiveness on real and synthetic noisy data.

Contribution

It presents a novel deterministic heuristic approach for DSC-MRI reconstruction, supported by mathematical analysis and comparative results.

Findings

Heuristic performs well on real images.

Effective on noisy artificial phantoms.

Offers computational advantages over compressed sensing.

Abstract

In this paper we propose a deterministic approach for the reconstruction of Dynamic Susceptibility Contrast magnetic resonance imaging data and compare it with the compressed sensing solution existing in the literature for the same problem. Our study is based on the mathematical analysis of the problem, which is computationally intractable because of its non polynomial complexity, but suggests simple heuristics that perform quite well. We give results on real images and on artificial phantoms with added noise.

Tables

Table 1. Table 1: Relative percentage errors of the real data set reconstructed by the algorithms with different number of measures.

Number of measures (on the total):	10%	20%	33%	50%
1	1.78	1.20	0.77	0.45
2	2.10	1.66	1.37	0.99
3	1.79	1.22	0.77	0.45
4	1.91	1.44	0.95	0.82
IHT	4.18	2.40	1.41	0.63
LCAMP	3.06	2.16	1.99	4.04

Table 2. Table 2: Relative percent errors of the simulated dataset reconstructed by the algorithms with different number of measures.

Number of measures (on the total):	10%	20%	33%	50%
1	1.24	1.08	0.91	0.70
2	1.23	1.00	0.90	0.79
3	1.24	1.07	0.91	0.70
4	1.19	1.02	0.89	0.76
IHT	10.36	5.96	3.46	3.60
LCAMP	2.52	1.98	2.88	-

Table 3. Table 3: Relative percent errors of the simulated dataset (with noise) reconstructed by the algorithms with different number of measures.

Number of measures (on the total):	10%	20%	33%	50%
1	1.80	1.56	1.32	1.01
2	1.79	1.45	1.30	1.15
3	1.80	1.56	1.32	1.01
4	1.74	1.48	1.31	1.14
IHT	12.67	5.82	3.51	3.34
LCAMP	3.59	2.87	3.58	-

Table 4. Table 4: Relative percent errors of the real sequence reconstructed by the modified algorithm 1 with different values of the parameter a 𝑎 a .

a	1.0	0.9	0.8	0.7	0.6	0.5	0.4	0.3	0.2	0.1	0.0
err	1.78	1.74	1.71	1.82	1.95	2.01	2.02	2.29	2.26	2.16	2.07

Equations

$$F=\left[\begin{array}[]{c}F_{J}\\ F_{R}\end{array}\right],\quad f=\left(\begin{array}[]{c}f_{J}\\ f_{R}\end{array}\right),\quad y=\left(\begin{array}[]{c}y_{I}\\ y_{S}\end{array}\right).$$

$$\begin{array}[]{ccc}WF^{T}\left(\begin{array}[]{c}f_{J}\\ 0\end{array}\right)&=&WF^{T}\left(f-\left(\begin{array}[]{c}0\\ f_{R}\end{array}\right)\right)\\ &=&y-W\left[F_{J}^{T}F_{R}^{T}\right]\left(\begin{array}[]{c}0\\ f_{R}\end{array}\right)\\ &=&y-W\left(\begin{array}[]{c}0\\ F_{R}^{T}f_{R}\end{array}\right).\end{array}$$

$$W\left(\begin{array}[]{c}0\\ F_{R}^{T}f_{R}\end{array}\right)$$

$$\left[\begin{array}[]{cc}W_{IJ}&W_{IR}\\ W_{SJ}&W_{SR}\end{array}\right],$$

$$e_{rec}=W_{IR}{F}_R^T{f}_R.$$

$$({\mathcal{W}}_{\psi }\circ {\mathcal{F}}^{-1}f)(j,k)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{+\infty }({\mathcal{W}}_{\psi }{e}^{i\omega \cdot })(j,k)f(\omega )d\omega ,j,k\in \mathbb{Z},$$

$$(\mathcal{F}\circ {\mathcal{W}}_{\psi }^{-1}y)(\omega )=\sum _{j,k=-\infty }^{+\infty }{y}_{j,k}{\hat{\psi }}_{j,k}(\omega ),\omega \in \mathbb{R}.$$

$$({\mathcal{W}}_{\psi }\circ {\mathcal{F}}^{-1}f)(j,k)=\frac{1}{\sqrt{2\pi }}({\mathcal{W}}_{\psi }{e}^{i\bar{\omega }\cdot })(j,k),j,k\in \mathbb{Z},$$

$$(\mathcal{F}\circ {\mathcal{W}}_{\psi }^{-1}y)(\omega )={\hat{\psi }}_{\bar{j},\bar{k}}(\omega )=2^{-\bar{j}/2}{e}^{-i{2}^{1-\bar{j}}\pi \omega \bar{k}}\hat{\psi }(2^{-\bar{j}}\omega ),$$

$$\begin{array}[]{ccl}(\mathcal{F}\circ\mathcal{W}_{\psi}^{-1}y)(\omega)&=&e^{-i2^{1-j_{1}}\pi\omega k_{1}}[2^{-j_{1}/2}\hat{\psi}(2^{-j_{1}}\omega)+\\ &&+2^{-j_{2}/2}\hat{\psi}(2^{-j_{2}}\omega)e^{-i2\pi\omega(2^{-j_{2}}k_{2}-2^{-j_{1}}k_{1})}].\end{array}$$

$${\bar{x}}^{(0)}=\bar{x},{\bar{x}}^{(k+1)}=a{\bar{x}}^{(k)}+(1-a)\tilde{x}$$

$$(\mathcal{F}\circ {\mathcal{W}}^{-1}y)(\omega )=2^{-\bar{j}/2}{e}^{-i{2}^{1-\bar{j}}\pi \omega \bar{k}}\hat{\psi }(2^{-\bar{j}}\omega ),\text{if}y={\delta }_{j,\bar{j}}\cdot {\delta }_{k,\bar{k}}.$$

Full text

††journal: Applied Numerical Mathematics

1 Introduction

The magnetic resonance images formation process (MRI) is based on radio-frequency energy absorption and emission of the hydrogen atoms nucleus: by varying the intensity of two orthogonal magnetic fields on a particular region, we can obtain an echo signal on a wide frequency spectrum. The physical image in the spatial coordinates, which represents the spatial distribution of hydrogen atoms density, is reconstructed by applying the bi-dimensional inverse Fourier transform of the echo frequency signal [8].

In this paper we consider in particular the Dynamic Susceptibility Contrast MRI (DSC-MRI) [12, 13]: a temporal sequence of images is reconstructed while a para-magnetic contrast liquid flows into the vessels of the human brain. This experiment is fundamental for the analysis of the emo-dynamics of the region. This sequence of images, whose purpose is to highlight the flow path of the contrast liquid, exhibits a strong spatiotemporal correlation, that can be exploited to predict which measurements will bring mostly the new information contained in the next frames.

This is crucial because in principle, for each frame, in order to reconstruct the whole spatial image we need to measure the complete spectrum in the frequency domain, but in practice, it is highly recommended a rapid-MRI scheme involving the smallest possible number of in-vivo measurements (i.e. a frequency under-sampling), so the time to sample each frame is reduced. In general, the sampling speed is an important issue in many applications of the MRI, e.g. some medical exams need several minutes to be done, and this can represent a limit to this kind of exams. Sampling speed is limited by constraints that can be physical (e.g. sampling frequency) and physiological (e.g. bloodstream). The focus of many current researches is to study methods to reduce the number of measurement samples needed for each frame without degrading the image quality [9]. For the DSC-MRI, the frequency under-sampling of single images does not imply automatically the decreasing of the duration of the whole clinic procedure (the flow of the contrast liquid has inevitable physical speed limits), but it can be exploited to make more frequent space or time acquisitions: the second option increases the time resolution and allows the analysis of fast dynamics not observed at this moment.

When the frequency domain is under-sampled, the Nyquist criterion is in general violated and the Fourier reconstruction shows aliasing artifacts. In DSC-MRI, however, by exploiting the redundancy in MR images, it is possible to define methods for rapid reconstructions that produce minimal artifacts in the recovered images. Several techniques can be used to reduce the number of measures and at the same time keep the alias in artifacts small: some methods generates incoherent or less visual apparent artifacts at the cost of a reduced signal-noise ratio (SNR), other methods exploit the redundancy of the signal to recover, in its domain or in other domains [9].

A novel set of techniques exploits the fact that, although the image to recover can be dense, its representation in other domains can be sparse. This is the case of medical images, that many transforms can make sparse, e.g. the wavelet transform [3, 7]. The possibility of reducing the memory space needed to store the information of an image has set up the question if it could be possible even to undersample the measured signal (in our case the frequency spectrum) in order to recover the image, apparently violating the Nyquist criterion. This approach goes under the name of Compressed Sensing [2, 4, 5, 6] and from the last recent years there’s a growing literature exploting this approach on several applications, and also the MRI [9, 16].

In the literature, notable results in sparse signals recovery have been achieved by Compressed Sensing (CS) techniques, and also for MRI applications [9, 16]. However, these methods are based on statistical properties rather than deterministic ones. In particular they do not make any assumption on the possible relations between sampled and recovered signals. This lack of information limits performances; a recently proposed solution is introducing deterministic constraints in the CS procedure [16]: the space-time correlation of the sequence of images allows to define locations contraints, i.e. the prediction of where the image will evolve mostly in the next frames. The difficult task is to choose the measurement points to form the minimal set that best describes this evolution. In CS this is made with a random approach.

In this paper we want to study analytically the compound Fourier + wavelet transform, involved in the processes of reconstruction and sparsification of MR images, in order to define a deterministic technique for a rapid-MRI. By studying the wavelet transform, known for producing sparse representations for natural images, we exploit deterministically the relations between the wavelet sparse representation of the recovered signal (MR image) and the frequency samples. We then compared image reconstructions implemented with deterministic and CS methods.

The article structure is the following: in Section 2 we define the problem of recovering a signal with sparse wavelet representation by undersampling its frequency spectrum; in Section 3 we study the involved transforms (two particular compositions of the Fourier and wavelet transforms); in Section 4 we describe four different heuristic approaches for solving the problem; in Section 5 we recall two well-known iterative algorithms for sparse recovery; in Section 6 we describe the datasets used for the simulations; in Section 7 we compare the results of all the implemented algorithms, while in Section 8 we propose some possible directions we can explore to obtain better results.

2 The problem

Let $x\in\mathbb{R}^{N}$ be the signal to recover in the spatial coordinates (for simplicity we suppose in this section that the signal be one-dimensional), and $N\gg 0$. Let $y\in\mathbb{R}^{N},y=Wx$ be the vector of its wavelet coefficients, where $W$ is the Discrete Wavelet Transform (DWT). Let $y$ be $n$-sparse, i.e. $y_{i}=0\hskip 5.69046pt\forall i\in I,\ I\subseteq\{1,..,N\},\ |I|=n,\ n\ll N,$ (except an arrangement of the terms). Let $f\in\mathbb{C}^{N},f=Fx$ be the vector of the Fourier coefficients of $x$, where $F$ is the Discrete Fourier Transform (DFT).

We want to find an optimal $m$-undersampling of $f$ for recovering $y$, i.e. given the number of measurements $m$, with $m\geq n$ and possibly also $m\ll N$, we want to determine the set of indexes $J\subseteq\{1,..,N\},|J|=m,$ such that the recovered signal \hat{y}=WF^{T}\left(\begin{array}[]{c}f_{J}\\ 0\end{array}\right), where $f_{J}=(f_{h})_{h\in J}$, leads to a minimal error with respect to the original signal $y$, measured with the norm $||\hat{y}-y||_{2}$.

Note that between the wavelet and Fourier coefficients it holds the relation $FW^{T}y=f$. Now, let $J$ a frequency undersampling, it holds $F_{J}W^{T}y=f_{J}$, where $F_{J}$ denotes the submatrix of the rows of $F$ that correspond to the non-zero coefficients $f_{J}$ of $f$. Let us use the notation

[TABLE]

By applying the recover $WF^{T}$ to the zero-filled $f_{J}$, we have

[TABLE]

So we need to find the set of indexes $R$ that makes minimal the norm of the non-zero coefficients of the vector W\left(\begin{array}[]{c}0\\ F_{R}^{T}f_{R}\end{array}\right). Being $y$ sparse, it becomes relevant only to minimize the $n$ coefficients at indexes $I$ of the error array

[TABLE]

that is, if we decompose $W$ as

[TABLE]

it means to minimize the norm of the vector

[TABLE]

The problem is in theory solvable combinatorially choosing the appropriate set $R$ (or $J$) of frequencies to disregard (or to measure). The number of possible choices for $J$, as $N$ and $m$ increase, becomes however too big to be solvable with brute-force search. As an example, the real images we analyzed had $N=128^{2}=16384$ pixel, and if we wanted to perform the 25% of the total possible measurements, the possible choices for $R$ or $J$ would be \left(\begin{array}[]{c}N\\ m\end{array}\right)=\left(\begin{array}[]{c}16384\\ 4096\end{array}\right)\approx 10^{4000}.

The combinatorial approach is therefore intractable, so we defined other algorithms that could obtain good results to the same problem nevertheless.

3 The compound transforms $\mathcal{W}\circ\mathcal{F}^{-1}$ and $\mathcal{F}\circ\mathcal{W}^{-1}$

The analytic study of the problem is based on studying the composition of two transforms: the Inverse Fourier Transform $\mathcal{F}^{-1}$ that reconstructs the image in the spatial domain from the measured frequencies, and the Wavelet Transform $\mathcal{W}$ that sparsifies the image in its wavelet coefficients. In particular we focused on their compositions $\mathcal{W}\circ\mathcal{F}^{-1}$ and $\mathcal{F}\circ\mathcal{W}^{-1}$, as $\mathcal{W}$ composed to the inverse of $\mathcal{F}$ shows the error in the wavelet domain of every frequency disregarded by the undersampling, and $\mathcal{F}$ composed to the inverse of $\mathcal{W}$ shows the frequency spectrum of the sparse wavelet coefficients to sample.

In order to keep the notation as clean and simple as possible, we considered $\mathcal{F}$ as the usual continue Fourier transform, while $\mathcal{W}$ as the wavelet transform defined by a discrete orthonormal basis of $L^{2}(\mathbb{R})$, since using the continuous wavelet transform would have requested to introduce some notions not relevant to our purposes.

Let $\psi\in L^{2}(\mathbb{R})$ be a discrete orthonormal mother wavelet, i.e. let $\{\psi_{j,k}:\psi_{j,k}(z):=$ $2^{j/2}\psi(2^{j}z-k),j,k\in\mathbb{Z}\}$ be an orthonormal basis of $L^{2}(\mathbb{R})$.

Then it is easy to prove that for all $f:\mathbb{R}\rightarrow\mathbb{C}$ we have:

[TABLE]

and vice versa for all $y:\mathbb{Z}^{2}\rightarrow\mathbb{R}$ we have:

[TABLE]

In particular, if we consider a single non-zero Fourier coefficient $f=\delta(\omega-\bar{\omega})$, then we have:

[TABLE]

that are wavelet representations of single pure sinusoids; so low frequencies produce mainly low resolution components in the error $e_{rec}$ (1), while high frequencies produce mainly high resolution components in the error (figure 2).

Vice versa, if $y=\delta_{j,\bar{j}}\cdot\delta_{k,\bar{k}}$ then:

[TABLE]

so high resolutions have a frequency spectrum scattered on high frequencies and low resolutions have a frequency spectrum peaked in low frequencies. Moreover, fixed a resolution $j$ but varying the position $k$, the frequency spectrum is the same in modulus (figure 3).

When we transform more then one non-zero wavelet coefficients, the modulus of the frequency spectrum shows the typical phenomenon of wave interference (figure 4): if $y=\delta_{j,j_{1}}\cdot\delta_{k,k_{1}}+\delta_{j,j_{2}}\cdot\delta_{k,k_{2}}$, then

[TABLE]

4 Heuristic algorithms for determining subsets of measurements

We present four algorithms that determine good subsets $J$ of measurements (undersampling). They all implement some possible heuristics, since the optimal solution cannot be computed analytically.

The algorithms receive as input some informations about the spatial vector $x$ to recover, at least the indexes of its non-zero wavelet coefficients $I$, and at most also its values. As output the algorithms compute a set $J$ of frequency indexes to be sampled among the whole set $\{1,...,N\}$.

In a DSC-MRI exam, the time-series frames $x^{(t)}$ to acquire , $t=1,...,T$, are strongly spatio-temporal correlated. So we give as input to the algorithm the data relative to the image $\bar{x}=\operatorname{mean}_{t=1,...,{\tau}}{x^{(t)}}$ calculated as the mean of the first $\tau\ll T$ frames recovered in a complete way (i.e. with as much as possible frequencies measured). The set $J$ determined are then used to sample every other frame to acquire.

We also considered a time-variable mask in order to evolve the set $J$ in accordance with the evolution of the physical process. In particular, we update $\bar{x}^{(k)}$ as follows:

[TABLE]

where $\tilde{x}$ is the last reconstructed frame. The results in section 7 show that this is meaningful in the DSC-MRI context.

Moreover, the non-measured frequencies of a frame (elements of $R=\{1,...,N\}\setminus J$), instead of being set to zero, are set to their value in the spectrum of the image mean of the first frames.

4.1 Algorithm 1

Let us calculate the $m\in\mathbb{N}$ frequencies with maximum modulus of the mean recovered image $\bar{x}\in\mathbb{R}^{N}$ (a priori $F\bar{x}\in\mathbb{C}^{N}$):

This algorithm computes the set $J$ of frequencies to measure by thresholding the ones with maximum modulus (figure 5).

Such an under-sampling $J$ obtains an approximation \hat{y}=WF^{T}\left(\begin{array}[]{c}f_{J}\\ 0\end{array}\right) for $y=W\bar{x}$ such that the errors $|f_{k}|WF^{T}_{\cdot,k}$ for every single excluded frequency $k\in R$ are minimal, since the frequencies $f_{k}$ set to zero are the ones with minimum modulus. This is true for excluding single frequencies but does not hold necessarily when considering the exclusion of groups of frequencies, so the error could not be globally optimal, but it is locally optimal, in the sense that given $\mu+1$ non-zero frequencies, if you want measure $\mu$ of them (excluding just one of them), the frequency to be discarded is precisely the one with the minimum modulus.

4.2 Algorithm 2

Suppose to know the distribution of the wavelet non-zero coefficients of $y=W\bar{x}$ among the different resolution orders, then we can calculate the frequencies that are expected to have the maximum modulus and add them to the set of measurements $J$.

Let us exploit the analysis of the transforms $\mathcal{W}$ and $\mathcal{F}$. We found the following formula for the frequency spectrum of a single wavelet coefficient:

[TABLE]

So it is clear that wavelet coefficients of the same resolution order $\bar{j}$ have the same frequency spectrum in modulus. The same formula highlights that varying the order of resolution the frequency spectrum is subjected to dilatation and contraction effects: the maximum modulus frequency of a frequency spectrum of a wavelet coefficient of order $j$ has modulus proportional to $2^{-j/2}$, while the support of the frequency spectrum dilates with the factor $2^{j}$. These observations are useful to define sets of measures for every single wavelet coefficients: given a single wavelet coefficient $y=\delta_{j,\bar{j}}\cdot\delta_{k,\bar{k}}$, and fixed $\bar{j},\bar{k}$, the cardinality $m^{\prime}=|J^{\prime}|$ of the set $J^{\prime}=\{i:|(FW^{T}y)_{i}|>threshold\}$ of frequencies with modulus greater then a fixed threshold is proportional to $2^{\bar{j}}\cdot 2^{-\bar{j}/2}=2^{\bar{j}/2}$.

Suppose that the wavelet representation $y=W\bar{x}$ of the signal $\bar{x}$ has only a set $I,\hskip 5.69046pt|I|=n$ of non-zero coefficients, and let $n_{j},j=0,...,ord_{max}$ the numbers of the non-zero coefficients in every subsequence of different resolution (let the minimum order trend subsequence have index [math], as the maximum order fluctuation sequence has index $ord_{max}$). We can now exploit the facts just described to compute the set $J$: for every resolution order of the wavelet representation $y$ we add to $J$ a group of frequencies with maximum modulus such that their cardinality is proportional to the number of non-zero wavelet coefficients, to the size of the spectrum support, and to the maximum modulus of the spectrum (figure 6).

The algorithm is then:

4.3 Algorithm 3

The analysis of the transforms $\mathcal{W}$ and $\mathcal{F}$ has however pointed out an important event that the previous algorithm does not consider: wavelet representations in which there are more than one non-zero coefficient exhibit frequency spectra with interference-like effects. So we defined a more refined algorithm that computes better measurements sets in respect of this occurrence. The idea is to iteratively compute frequency spectra of an increasing number of wavelet coefficients, in which we can select a new frequency to be added to the set $J$.

Let $y=W\bar{x}$ be the wavelet representation of the signal $\bar{x}$. Let $I=\{i:y_{i}\neq 0\}$ be the set of indices of the non-zero wavelet coefficients, $|I|=n$. Let $I$ be sorted with respect to the descending values $|y_{i}|,\hskip 5.69046pti\in I$. Let $y_{(l)}=\sum_{i\in I,i<l}y_{i}e_{i}$ (i.e. $y$ with every element smaller then the $l$ elements of maximum modulus set to zero). Then for every $l=1,...,n$ we calculate the spectrum $f=FW^{T}y_{(l)}$, and we add $\arg\max|f|$ to the set $J$; while this frequency is already an element of $J$, we set it to zero and add $\arg\max|f|$ to $J$ again (figure 7).

The algorithm is the following:

4.4 Algorithm 4

Suppose that the positions of the non-zero wavelet coefficients of the signal $\bar{x}$ are known. Let us find the frequencies that “influence more” such coefficients: with this expression we mean the frequencies that have a wavelet spectrum (by the transform $\mathcal{W}\circ\mathcal{F}^{-1}$) with values over a fixed threshold in the same positions of the non-zero coefficients of $y=W\bar{x}$.

Let $I=\{i\in\{1,..,N\}:y_{i}\neq 0\}$ the positions of the non-zero coefficients of $y$, then the frequencies $f_{h}$ we want to measure are the ones that transformed by $WF^{T}$ have at least one of the components in $I$ over a fixed threshold. Since the $h$-th frequency corresponds to the $h$-th row of $FW^{T}$ and the $i$-th component of $y$ to the $i$-th column of $FW^{T}$, the problem reduces to compute the indices $h\in\{1,..,N\}$ such that $\max(WF^{T}e_{h})_{i\in I}>s_{1}$, with $s_{1}$ a fixed threshold (figure 8).

Dually, we discard the frequencies that would influence more the zero wavelet coefficients of the image. We compute the frequency spectrum of the zero wavelet coefficient set to the unitary value (but they could also be weighted differently) and we exclude the frequencies whose modulus exceeds a fixed value. Since $S=\{i\in\{1,..,N\}:y_{i}=0\}$, then we discard the frequencies correspondent to the indices $h\in\{1,..,N\}$ such that $(FW^{T}y_{S})_{h}>s_{2}$, with $s_{2}$ a fixed value (figure 9).

The algorithm is as follows:

5 Iterative algorithms

In order to compare the results of the previous algorithms with other sparse recovery methods already studied in literature, we considered two iterative algorithms from the compressed sensing theory: the Iterative Hard Thresholding (IHT) algorithm [1] and the Location Constraint Approximate Message Passing (LCAMP) algorithm [16]. Both algorithms use $F_{J}W^{T}$ as the sensing matrix and an undersampling set $J\subset\{1,...,N\}$ that has a casual distribution concentrated in low frequency.

The IHT algorithm:

where $H_{k}(z)$ is the hard thresholding operator on $z$ that set all the $z$ components to zero except the $k$ components with the maximum modulus.

The LCAMP algorithm:

The stop criterion in both algorithms can consist in a limit of the iteration number or a convergence criterion for $\hat{x}$.

6 Data

6.1 Real data

A series of $T2^{*}$ weighted volumes ($128\times 128\times 20$ voxels with in plane resolution of $1.8\times 1.8$ mm and slice thickness of 4 mm, $TR=1500$ ms, $TE=40$ ms) were acquired from a patient affected by artero-venous malformation, after a bolus of gadolinium (Gd-TPA).

6.2 Simulated data

We simulated a set of a DSC-MRI sequence, using a $256\times 256$ pixels Shepp-Logan phantom [14], commonly used in the simulations in the field of MRI. To the phantom we added white noise to decrease the signal to noise ratio (SNR) to $15$ dB. Moreover, the simulation of the local transit of contrast agent was obtained by adding to specific regions a monovariate gamma function $y_{gad}(t)$ with fixed parameters with a secondary gamma function $y_{rec}(t)$ representing the recirculation [15]; to the gamma function a log-normally distributed noise was added.

7 Comparing the reconstructions

We ran all the above algorithms on three sets of time-series images: a sequence of images from a real medical DSC-MRI exam on the encephalic region of a arterio-venous malformation (AVM) patient, and the Shepp-Logan phantom imitating an encephalic region with the addition of localized flow of a contrast agent, with and without noise.

For every sequence we show the reconstructions of one particular frame made by all the algorithms described in the previous sections and a table of the relative percent errors (2-norm) of any single frame averaged over the whole sequence. The error of a particular reconstructed (under-sampled) frame from the original frame is computed only on the central encephalic region (the background is discarded).

7.1 Real data

The algorithms 1 and 3 give the best reconstructions with any number of measurements we tried. The algorithms 2 and 4 performs better than the CS algorithms while the number of measures is lower than the 33% of the total number of possible measures. Note that between the errors corresponding to the CS algorithms and the algorithms 1 or 3 there is a factor $1.5$/$2$, a little less then the factor $2$ between the necessary number of random measures and the sparsity of the signal to recover in the estimates in the CS theory [4]. As we saw, the algorithms 1 and 3 do not produce the optimal set $J$, maybe such an optimal set would achieve a factor 2. LCAMP algorithm obtain good results only with a few measures, the results degrade as their number is increased: this problem is intrinsic in the AMP algorithm: it recovers a signal correctly only when the quantity $|\mbox{\rm supp }{y}|/\dim{y}$ is lesser then a certain increasing function of $\dim{y}/|J|$ [10].

The recovered images in the figure 10 reflect the error table: note how the deterministic algorithms perform better then the CS algorithms recovering the higher resolution details.

7.2 Phantom without noise

The deterministic algorithms give the best reconstructions then the CS algorithms with any number of measures we tried. All of them seem to obtain similar results. Between the two CS algorithms LCAMP has the better results at a lower number of measures.

The recovered images in the figure 11 reflect the error table. Note how the first four reconstructions are close to the original frame while the number of measures are only the 20% of all the possible measures.

7.3 Phantom with noise

The deterministic algorithms give the best reconstructions then the CS algorithms with any number of measures we tried. All of them seem to obtain similar results. Between the two CS algorithms, LCAMP has the best results at a lower number of measures. The errors are globally worse then the phantom without noise.

The recovered images in the figure 12 are very similar to the previous ones, the added noise is not sufficiently intense to be visible.

7.4 Adaptive mask

As we mentioned in section 4, we tried to identify better measure sets $J$ by adapting the image that the algorithms receive as input with the last reconstructions.

We modified the algorithm 1 in a way that it updates the input mask $\bar{x}$ for every instant $t>\tau$ according to the formula (2).

The table 4 shows the relative percent errors, in the $2$-norm, between reconstructed and original images, averaged over the whole sequence, in the case of real images with 10% of measures.

As $a$ decreases, the error lightly decreases too, then it increases again. The graph in figure 13 shows the error time evolution varying the parameter $a$. Note that the most “prompt” algorithm ($a=0.0-0.2$) presents a greater global error, but it is sensibly more accurate in the peak corresponding to the moment at which the contrast fluid flows in the region, with an error reduction of about the 8%. Then one can improve the tracking of the dynamic phase of the test by updating the mask $\bar{x}$. It is worth noting that from a diagnostic point of view, the late portion of the examination (approx. from frame 50 onward), is usually discarded from the analysis. Algorithms that are able to provide better performance in the first part of dynamic phase (rising edge, peak and washout) are preferable.

8 Conclusions

The adoption of techniques deriving from the mathematical analysis of the DSC-MRI problem enabled us to develop algorithms for the reconstruction of MRI images that perform sensibly better than the random based approaches cited in the literature.

• 2.

The results are sensitive to the error function adopted, e.g. the best images are not necessarily the most accurate in the $2$-norm. Improvements in this direction must be obtained, in our opinion, with a joint work with the physicians.

• 3.

We have also noticed that these results can still be improved by introducing a dynamic model of the problem, that could be improved with an estimator of the direction and speed of the contrast propagation in the time-spatial domain rather than only in the time domain, see e.g. [11] for an example. With this aim, it should be relevant to exploit the fluid dynamics of the contrast liquid to sharpen the prediction.

Acknowledgments

The first author is supported by the project ex-60% of the University of Padua, named “Metodi numerici per problemi inversi in meccanica computazionale”.