IMAGING

The invention relates to a method for performing parallel magnetic resonance imaging (pMRI) of a body, including parallel, dynamical (time-resolved) magnetic resonance imaging, such as functional MRI (fMRI).

Reducing global acquisition time is of main interest in medical magnetic resonance imaging (MRI), and even more when dynamic imaging such as fMRI is concerned. Actually, a short acquisition time allows improving the spatial/temporal resolution of acquired fMRI data, which leads to a more efficient statistical analysis. In addition, by reducing the global imaging time, some additional artifacts caused by the patient motion can be avoided. For this reason, parallel imaging systems have been developed: multiple receiver surface coils with complementary sensitivity profiles located around the underlying object or body are employed to simultaneously collect in the frequency domain (i.e. the so-called k-space), data sampled at a rate R times lower than the Nyquist sampling rate along at least one spatial direction, i.e. the phase encoding one; R is usually called the "reduction factor". Therefore, the total acquisition time is R times shorter than with conventional non parallel imaging. A reconstruction step is then performed to build a full Field of View (FOV) image by unfolding the undersampled "elementary" images acquired by the individual receivers. This reconstruction is a challenging task because of the low Signal to Noise Ratio (SNR) in parallel MRI (pMRI) caused by aliasing artifacts related to the undersampling rate, those caused by noise during the acquisition process and also the presence of errors in the estimation of coil sensitivity maps.

The Simultaneous Acquisition of Spatial Harmonics (SMASH) [Sodickson et al., 1997] was the first reconstruction method, operating in the k-space domain. It uses a linear combination of pre-estimated coil sensitivity maps to generate the missing phase encoding steps.

Some other k-space based reconstruction techniques have also been proposed like GRAPPA (Generalized Autocalibrating Partially Parallel Acquisitions) [Griswold et al., 2002], and SENSE (Sensitivity Encoding) [Pruessmann et al., 1999]. SENSE is a two-step procedure relying first on a reconstruction of reduced FOV images and second on a spatial unfolding technique, which amounts to a weighted least squares estimation. This technique requires a precise estimation of coil sensitivity maps using a reference scan (usually a 2D Gradient-Echo (GRE)). It is presently the most frequently employed pMRI technique, applied in particular to brain and cardiac imaging.

For a general overview of reconstruction methods in pMRI see [Hoge et al., 2005].

SENSE is often supposed to achieve an exact reconstruction in the case of noiseless data and perfect coil sensitivity maps knowledge, which is also true for all above mentioned methods. However, in practice, the presence of noise in the data and inaccuracies in the estimation of coil sensitivity maps are unavoidable and make the reconstruction problem ΊΙΙ- conditioned.

As image reconstruction is an ill-posed inverse problem, regularization techniques are commonly applied to better estimate the full FOV image. Most of these techniques operate in the image domain; in particular, this is the case for Tikhonov regularization [Ying et al., 2004], which uses a quadratic penalty term either to promote smoothness constraints or to account for the squared difference between the reconstructed image and an a priori reference image. Despite the use of regularization, however, high reduction factors (exceeding a value of R=2) are generally considered unfeasible when low magnetic field intensities (up to 1.5 Tesla) are used, since the reconstructed images are affected by severe aliasing artifacts.

In [Chaari et al.2008] and [Chaari et al. 2009] the present inventors have described a method of performing regularized image reconstruction in parallel MRI, using a wavelet-based regularization scheme, allowing to increase the reduction factor R.

The present invention aims at providing several improvements of said method, including extending it to dynamical imaging (e.g. fMRI) and making it fully or partially auto-calibrated (or "unsupervised"). An object of the present invention is then a method of parallel magnetic resonance imaging of a body, comprising:

acquiring a set of elementary magnetic resonance images of said body from respective receiving antennas having known or estimated sensibility maps and noise covariance matrices, said elementary images being under-sampled in k-space; and

performing regularized reconstruction of a magnetic resonance image of said body;

wherein said step of performing regularized reconstruction of a magnetic resonance image is carried out in a discrete frame space by minimizing a cost function comprising:

an error term, representative of a likelihood of a reconstructed image, given said acquired elementary images; and

a frame penalty term, representative of a deviation between an actual statistical distribution of frame coefficients of said reconstructed image and an a priori distribution of said coefficients;

said a priori distributions of the frame coefficients of the reconstructed image being estimated on the basis of an auxiliary magnetic resonance image of said body.

According to different embodiments of said method:

Said error term can be representative of a neg-log- likelihood of said reconstructed image, given said acquired elementary images.

Said step of performing regularized reconstruction of a magnetic resonance image of said body can be carried out by maximizing, in said frame space, a full posterior distribution of a set of frame coefficients defining an image of the body, given said acquired elementary magnetic resonance images and said a priori distribution of the frame coefficients.

Said auxiliary magnetic resonance image of said body can be reconstructed from said acquired elementary magnetic resonance images. More particularly, said auxiliary magnetic resonance image of said body can be reconstructed from said acquired elementary magnetic resonance images using a SENSitivity Encoding - SENSE - algorithm. Even more particularly, said auxiliary magnetic resonance image of said body can be reconstructed from said acquired elementary magnetic resonance images using an algorithm chosen between: an unregularized SENSE algorithm; a SENSE algorithm regularized in image space; and a SENSE algorithm regularized in k-space.

A generalized Gauss-Laplace a priori statistical distribution of said frame coefficients can assumed, and parameters of said distribution can be estimated on the basis of said auxiliary magnetic resonance image of said body, using a maximum-likelihood or a posterior mean estimator.

Said error term can be a quadratic mean error term.

Said acquired elementary images can be three- dimensional images, and said step of performing regularized reconstruction of a magnetic resonance image can be carried out in a discrete three- dimensional frame space. More particularly, said acquired three-dimensional elementary images can be obtained by stacking bi-dimensional elementary images of slices of the object to be imaged.

Said step of performing regularized reconstruction of a magnetic resonance image can be based on a redundant wavelet frame representation. Alternatively, said step of performing regularized reconstruction of a magnetic resonance can be based on a non-redundant wavelet representation.

Said cost function can also comprise at least one spatial domain penalty term chosen between: a total variation norm of the reconstructed image; and a convex constraint.

Another object of the present invention is a method of performing dynamical and parallel magnetic resonance imaging of a body, comprising:

acquiring a set of time series of elementary magnetic resonance images of said body from respective receiving antennas having known or estimated sensibility maps and noise covariance matrices, said elementary images being under-sampled in k-space; and performing regularized reconstruction of a time series of magnetic resonance images of said body;

wherein said step of performing regularized reconstruction of a time series of elementary magnetic resonance images is carried out by minimizing a cost function comprising:

an error term, representative of a likelihood of each reconstructed image, given corresponding acquired elementary images; and a temporal penalty term, representative of a pixel-by-pixel or voxel-by-voxel difference between consecutive image of the series.

According to different embodiments of said method:

- Said temporal penalty term can be based on an edge- preserving function, more particularly a convex edge-preserving function and even more particularly on an L _p norm with p> 1 and preferably 1<p<1.5

Said temporal penalty term can be given by the sum of a first partial temporal penalty term and a second partial temporal penalty term, wherein: the first partial temporal penalty term is representative of pixel-by- pixel or voxel-by-voxel differences between each even-numbered image of the series and a preceding odd-numbered image; and the second partial temporal penalty term is representative of pixel-by-pixel or voxel-by-voxel differences between each odd-numbered image of the series and a preceding even- numbered image; said cost function being minimized by using proximity operators for said first and second partial temporal penalty terms.

Said step of performing regularized reconstruction of a magnetic resonance image can be carried out in a discrete frame space, and said cost function can also comprise a frame penalty term, representative of a deviation between statistical distributions of frame coefficients of each reconstructed image and an a priori distributions of said coefficients; said a priori distributions of the frame coefficients of the reconstructed images being estimated on the basis of an auxiliary magnetic resonance image of said body. In this case, said elementary images can be three-dimensional images, and said discrete frame space can be a discrete three-dimensional frame space. The method can further comprise a step of automatically determining, using a maximum-likelihood estimator; a weighting parameter of said temporal penalty term. More precisely, the method can comprise estimating said weighting parameter of the temporal penalty term for each pixel or voxel, or set of neighboring pixels or voxels, of the image to be reconstructed. Moreover, said temporal penalty term can be based on an L _p norm and said weighting parameter of the temporal penalty term and the parameter p can be jointly determined using said maximum-likelihood estimator. In particular, said weighting parameter of the temporal penalty term and the parameter p can be jointly determined using said maximum-likelihood estimator under the constraint p>1.

Said error term can depend on geometrical parameters defining a rigid transformation of each of said elementary magnetic resonance images with respect to an elementary magnetic resonance image taken as a reference, and wherein said step of performing regularized reconstruction of a time series of elementary magnetic resonance images is carried out by minimizing said function also with respect to said geometrical parameters.

Said elementary images can be acquired by echoplanar imaging.

- Said elementary images can be under-sampled with a reduction factor higher or equal to 4.

Additional features and advantages of the present invention will become apparent from the subsequent description, taken in conjunction with the accompanying drawings, which show:

- Figure 1 , empirical histograms of real and imaginary parts of wavelet coefficients of a magnetic resonance image of a human brain;

Figures 2 and 3, nine anatomical axial slices of a human brain obtained using the Tikhonov-regularized SENSE pMRI reconstruction method known from prior art, with R=2 and R=4 respectively;

- Figures 4, an anatomical axial slice of the same human brain obtained using the TV-regularized SENSE pMRI reconstruction method known from prior art, with R=2 (left panel) and R=4 (right panel); Figures 5 and 6, nine anatomical axial slices of the same human brain obtained using an autocalibrated 2D wavelet transform-based regularization scheme (UWR-SENSE) according to an embodiment of the invention, with R=2 and R=4 respectively;

- Figures 7 and 8, nine anatomical axial slices of the same human brain obtained using a constrained, autocalibrated 2D wavelet transform-based regularization scheme (CWR-SENSE) according to another embodiment of the invention, with R=2 and R=4 respectively;

Figure 9, nine anatomical axial slices of the same human brain obtained using an autocalibrated combined wavelet-total variation regularization scheme according to another embodiment of the invention, with R=4 and using image decomposition on a non-redundant orthonormal wavelet basis;

Figure 10, nine anatomical axial slices of the same human brain obtained using an autocalibrated combined wavelet-total variation regularization scheme according to another embodiment of the invention, with R=4 and using image decomposition on a redundant wavelet frame constituted by the union of two orthonormal bases;

Figure 11 , three anatomical axial slices of the same human brain obtained using: conventional non-parellel MRI (top row), mSENSE parallel MRI reconstruction with R=4 (middle row) and an autocalibrated 3D wavelet transform-based regularization scheme (3D-UWR- SENSE) according to another embodiment of the invention (bottom row);

Figure 12, three anatomical sagittal slices of the same human brain obtained using: conventional non-parellel MRI (top row), mSENSE parallel MRI reconstruction with R=4 (middle row) and an autocalibrated 3D wavelet transform-based regularization scheme (3D-UWR- SENSE) according to another embodiment of the invention (bottom row);

Figure 13, three anatomical axial slices of the same human brain (upper row) and magnified details thereof (bottom row) obtained using: TV regularization (left); combined wavelet-total variation regularization scheme according to an embodiment of the invention using image decomposition on a non-redundant wavelet basis (center) and combined wavelet-total variation regularization scheme according to an embodiment of the invention using image decomposition on a redundant wavelet frame constituted by the union of two orthonormal bases (right);

- Figure 14, Axial, Coronal and Sagittal slices of a human brain, acquired using Ecoplanar (EPI) fMRI and reconstructed using the mSENSE method, known from prior art, and an autocalibrated 2D wavelet transform-based regularization scheme (4D-UWR-SENSE) according to an embodiment of the invention;

- Figure 15, subject-level student-t maps, superimposed to anatomical imaging, of the aC-aS contrast detected using EPI fMRI; data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE methods, respectively, with R=2 (top of the figure) and R=4 (bottom of the figure); sagittal, coronal and axial views are displayed;

- Figure 16, subject-level student-t maps of the aC-aS contrast for two subjects, reconstructed at R=2 using the mSENSE, UWR- SENSE and 4D-UWR-SENSE methods, respectively; sagittal, coronal and axial views are displayed;

Figure 17, subject-level student-t maps, superimposed to anatomical imaging, of the Lc-Rc contrast detected using EPI fMRI; data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE methods, respectively, with R=2 (top of the figure) and R=4 (bottom of the figure); sagittal, coronal and axial views are displayed;

Figure 18, subject-level student-t maps of the Lc-Rc contrast for two subjects, reconstructed at R=4 using the mSENSE, UWR- SENSE and 4D-UWR-SENSE methods, respectively; sagittal, coronal and axial views are displayed;

Figure 19, group-level student-t maps for the aC-aS contrast, where data have been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSE for R = 2 and R = 4; sagittal, coronal and axial views are displayed; and Figure 20, group-level student-t maps for the Lc-Rc contrast, where data have been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSE for R = 2 and R = 4; sagittal, coronal and axial views are displayed.

Before describing the invention in detail it will be necessary to recall some basic facts about pMRI (and, in particular, SENSE).

MRI is an imaging technique that can proceed either in 2D (two-dimensional) or directly in 3D (three-dimensional) depending on the involved RF pulse design. In the 2D case, a volume is covered using adjacent slices. PMRI may adapt to both situations since this method makes the k- space scanning faster, whatever its definition (i.e. in 2 or 3D). For the sake of simplicity, we only present the 2D case in the following.

An array of C coils is employed to measure the spin density p into the object under investigation (e.g. a brain, in the head of a patient). Image acquisition is based on specific imaging sequences; in exemplary embodiments of the invention, anatomical MRI is based on the 3D MPRAGE sequence while functional MRI involves is performed using 2D echoplanar imaging (EPI). The following description will be focused on the 2D case; then, the signal d _c received by each coil I (1<c<C) is the Fourier transform of the desired 2D field p weighted by the coil sensitivity profile s _c , evaluated at some locations k = (k _y , k _x ) ^T in the k-space (the apex ^T means transposition). The received signal is therefore defined by the sampling scheme:

■rj

where ' °C is a realization of an Additive White Gaussian Noise (AWGN) and r ⁼ (y. ^X ) ^T is the spatial position in the image domain. For the sake of simplicity, a Cartesian coordinate system is generally adopted. In its simplest form, SENSE imaging amounts to solving a one-dimensional inversion problem due to the separability of the Fourier transform. Let Ay = Y/R be the sampling period where Y is the size of the field of view (FOV) along the phase encoding direction, let y be the position in the image domain along the same direction, x the position in the image domain along the frequency encoding direction and R < L the reduction factor. A 2D inverse Fourier transform allows recovering the measured signal in the spatial domain. By accounting for the undersampling of the k-space by R, (1) can be re-expressed in the following matrix form: d(r)— S(r)p(r) 4- n(r)

(2) where:

In equation (2), (n(r)) _r is a sequence of circular zero-mean Gaussian complex-valued random vectors. These noise vectors are i.i.d. (independent, identically distributed) and spatially independent with covariance matrix Ψ of size C * C. In practice, Ψ is estimated by acquiring data from all coils without radio frequency pulses, and its generic entry Ψ(θι,θ2) corresponding to the covariance between the two coils Ci and C2 is given by:

= ^ ^ d^ fc ^ fcx), V(c _{1 > C2} ) G {1, . . . , C} ² (4). where (□) ^* stands for the complex conjugate. Note that a statistical model can be assumed for matrix Ψ such as the statistical independence between coils or a nearest-neighbors statistical dependence model. In the fir se, matrix

Ψ becomes diagonal and the coil-dependent variances can be estimated using the Mean Absolute Deviation (MAD) technique described by [Donoho, 1995], which relies on an additive Gaussian white noise (AGWN) assumption. In the simplest case, Ψ can be simply be taken equal to the identity matrix. In the presence of non-zero off-diagonal terms, the empirical covariance given in Eq (4) is used to estimate them.

Therefore, the reconstruction step consists in inverting (2) and recovering from d(r) at spatial positions r = (y, x) ^T . Note that the data (dc)i≤i≤c and the unknown image are complex-valued, although | | is only considered for visualization purposes.

A simple reconstruction method also called the SENSE approach [Pruessmann et al., 1999], is based on the minimization of the Weighted Least Squares (WLS) criterion.

The objective is to find a vector at each spatial location r such that: wxs(* = arg min JwLs( (t ¹ ))

^Pir) (5) where Hi d{r) - S(r)p(r) % _{( > ( )H stands for †he} transposed complex conjugate, (.)* stands for the pseudo-inverse and defines a norm on As discussed above, this inverse problem is generally ill- posed and requires regularization, e.g. Tikhonov regularization. The regularization process typically consists in computing p _PW Ls( ^r ) as the minimizer of the following Penalized Weighted Least Squares (PWLS) criterion:

JPWLSW' ^' ) ) =

where IR is the R-dimensional identity matrix. The regularization parameter K>0 ensures a balance between the closeness to the data and the penalty term, which controls the deviation from a given reference vector p _r (r). The solution PPWLS admits the following closed-form expression:

PPWLS(^) = Pr(^)+

(S ^H (r)* ^_1 S(r) + K/ _R ) ^"1 S ^H (r) *- ¹ (rf(r) - S(r)p _T (r))

(7).

Note that the accuracy of the solution depends on the reference vector and the choice of the regularization parameter κ.

Tikhonov regularization is known to introduce blurring in the image.

To overcome this limitation, "edge-preserving" penalty terms have been proposed, which are applied in the image domain and make the regularization more efficient by limiting blurring effects and preserving the image boundaries. However, the introduction of these terms may lead to a non-differentiable optimization problem which is not always easy to solve numerically.

As mentioned above, in [Chaari et al.2008] and [Chaari et al.2009] the present inventors have proposed a new regularization scheme; based on wavelet transforms (WT); which will be summarized here.

Indeed, in SENSE-based reconstructed images, well spatially- localized artifacts appear as distorted curves with either very high or very low intensity, and the WT has been recognized as a powerful tool that enables a good space and frequency localization of useful information. A general introduction to wavelets and wavelet transforms is provided by [Pesquet-Popescu, Pesquet].

In what follows, T stands for the WT operator. It corresponds to a discrete decomposition onto a separable 2D M-band wavelet basis performed over j _max resolution levels. The objective image p of size Y *X can be viewed as an element of the Euclidean space with K = Y ^χ X endowed with the standard inner product < · | · > and norm || · ||.

Let (e _/ J _f <i(<Kbe the considered discrete wavelet basis of the space C ^K . The wavelet decomposition operator T is defined as the linear operator:

T: C ^K → C ^K

p f→ ( {p I e _k ) )i<k<K

The adjoint operator T ^* serving for reconstruction purposes is then defined as the bijective linear operator:

T ^* : C ^K → C ^K

K

fc=l

The resulting wavelet coefficient field of a target image function p is defined by ζ = (^a,k)i<k≤Kjmax , ^oj,k)i≤j≤jmax,i≤k≤Kj) where Kj = KM ^~2j is the number of wavelet coefficients in a given subband at resolution j (by assuming that Y and X are multiple of M ^jmax ) and the coefficients have been reindexed in such a way that denotes an approximation coefficient at resolution level j _max and ζ _0, ^ denotes a detail coefficient at resolution level j and orientation o e 0 = {0 M— 1} ² \{(0,0)}.

In the dyadic case (M = 2), there are three orientations corresponding to the horizontal, vertical or diagonal directions. When an orthonormal wavelet basis is considered, the adjoint operator T ^* reduces to the inverse WT operator T ^"1 and the operator norm ||T|| of T is equal to 1.

An estimate of the target image p is generated through the reconstruction wavelet operator T ^* . Let ξ be the unknown wavelet coefficients such that = Τ*ξ. The aim is to build an estimate ζ of the vector of coefficients from the observations d. This is based on a Bayesian approach relying on suitable priors on the wavelet coefficients.

Given the observation model in Eq. (2.) and the assumptions regarding the noise (i.i.d. circular Gaussian with zero-mean and between-coil correlation matrix Ψ), the likelihood function factorizes over pixels lying in the Y x X FOV: p{d j T ^* 0 = H p(d{r) J p(r)) o JJ exp(-J Ls(p(r)))

reiI....,V7R}x{ 1....JO Γ€{Ϊ,-.. ₍ }7Κ1Χ{1,-,Α ^' >

ae exp (-J _L (T ^* ))

(8) where p = Τ ^* ζ and

re{i,..,y} x{i,...,X>

(9)

Let f be the prior probability density function (pdf) of the image in the wavelet domain. It will be assumed here that the real and imaginary parts of the wavelet coefficients are independent. It will also be assumed that the real (resp. imaginary) parts of the wavelet coefficients are independent and identically distributed (i.i.d.) in each subband. Their statistical characteristics may however vary between two distinct subbands. Furthermore, by looking at the empirical distributions of the real and imaginary parts of the considered wavelet coefficients, the present inventors have noticed that their empirical histograms are well-fitted by a "Generalized Gauss-Laplace" (GGL) distribution, which presents a single mode and whose shape varies between the Gaussian and Laplacian densities. The corresponding pdf reads:

" ^{' ' ;} V 2ΪΓ crfc( -4 ) where a e E¾ ₊ and β e M ^* ₊ are hyper-parameters to be estimated. Figure 1 illustrates the empirical histograms of real and imaginary parts of the horizontal detail subband at the first resolution level using the dyadic (M = 2) wavelet decomposition with Daubechies filters of length 8. This figure shows also that the adopted GGL distribution better fits the empirical histogram than a Generalized Gaussian (GG) pdf.

At the coarsest resolution level j _max , the distributions of both the real and imaginary parts of the approximation coefficients can be modeled by a Gaussian distribution since they belong to a low frequency subband.

Due to its familiarity and simplicity, the MAP (Maximum A

Posteriori) estimator will be used for the estimation purpose, although different choices are possible. Based on the prior and the likelihood given hereabove, the MAP estimator is computed by maximizing the full posterior distribution or minimizing its negative log-likelihood:

( ^AP = _a rg niax( lii/><d | T ^* < ) 4- In /(ζ) )

e ^£K (11) or equivalent^ by minimizing the following criterion:

ζ MAP _{= niin} j _WT (<;)

ee*

= arg min(JLC ) + JviQ)

ieC* (12) with:

∑ < CoJfc) +∑∑∑ *

fc)l ² + a¾|Ira(C.j.fc)l

.Willi

j.k ,

(13). Hereabove, Re( ) and lm( ) (or ^Re and ^lm ) stand for the real and imaginary parts, respectively. The prior parameters (hyperparameters) ao,r(ao,j ^Re , a ₀ , ^m )e (M* ₊ ) ² , Mo,j" ( ) ² and aoj=(a _0, j ^Re , a _0i f ^m )≡ (R+) ² are unknown and need to be estimated. The estimation of hyperparameters is an important part of the invention and it will be discussed extensively later.

While JL is differentiate with a Lipschitz-continuous gradient, JP is not differentiate. Therefore, although the penalty function Jwr is convex, the optimization procedure cannot rely on conventional convex optimization techniques like the pseudo-conjugate gradient.

Optimization can be performed using a generalized form of the iterative optimization procedure developed in [Daubechies et al., 2004; Chaux et al., 2007], which is based on the Forward-Backward (FB) algorithm and on the concept of proximity operator, generalized to the case of functions of a complex variable.

For the function

Φ : C ^K →}— DO. +oo]

x (-»· ^¾ (Ηί;(¾) ) + ^1η ' (lm(r)} where φ ^Κβ and φ ^ΙηΊ are functions in r ₀ ( ^{a K} ) and Re(x) (resp. Im(x)) is the vector of the real parts (resp. imaginary parts) of the components of x E C the proximity operator is defined as

ρΐΌ φ : C ^K → C ^K

x +-> + ·φΓ0Χφϊιη(Ιιη(:Γ)) .

By extending the algorithm in [Chaux et al., 2007] to the complex case, a minimizer of Jwr can be iteratively computed according to the algorithm below (Algorithm 1), where the gradient of JL is first calculated, and then the frame coefficients are updated. It can be noticed that λ _η and γ _η correspond to relaxation and step-size parameters, respectively. Algorithm 1:

Let (7n)n>o and (λ„) _η> ο be sequences of positive reals.

1: Set « = 0 aad £ > 0. Initialize and set jW ~ Λντ(ζ ^<η) ).

¾ repeat

3: Reconstruct the image by setting pl ⁿ — Τ*ζ( ^η >.

4: Compute the image u'"- ¹ such that:

Vr€ {1. Y/R} x {I A ^' },

u ⁱⁿ⁾ (r) = 25 ^Η (τ)Φ ¹ (Sfr)/>< ⁿ >(r) - d(r)),

where the vector «'"'{*') is defined from u ^in> in the same way as (r) is defined from 7J (see Eq. (2.17».

5: Determine the wavelet coefficients v' ⁿ⁾ = Xu ^1' "' = (t» _et (t>oj)o*0 KJ^ _^ ) of u' ^T1 & Update the approximation coefficients of the reconstructed image pl ⁿ⁺ >

v* e { i, ... , },dT ^1} = a + ^(Ρ^Φ. «& ^} - J> - <&>) .

7: Update the detail coefficients of the reconstructed image p ⁽ " ⁺¹ >:

Vo e O. Vj e {1, ... , },V*€ {!,..., ft ^' ,},

» + 1> _ ·(«) . ..^ ₍ An) _ In) v ,,{«) \

8: Compute '" ' ^< - JwT« ^{(n+1 !} )- 9: n <- n + 1

10 _; until - < _ε ,7<· ^» -υ

11: return = Τ' ^η

For every rc [l y} x{l, ... ,A j _{|et Θγ≥0 be tne} maximum eigenvalue of the Hermitian positive semi-definite matrix S"(v)1>- ^l S{r) _{and |et} i = max _r e _{1 ..... _r) x _{ i A} β, > 0 _Jo guarantee the convergence of Algorithm 1, the step-size and relaxation parameters are subject to the following conditions:

(i) inf _n>0 7„ > 0 and sup _n>0 „ < ο ψ,

(ii) inf _n >o λ _η > 0 and sup _ri>0 λ _η < 1.

As discussed in [Chaari et al.2009], the optimization problem (11), (12) can be modified to include additional (preferably convex) constraints, e.g. upper and lower bounds of the image intensity value.

The present inventors have found that using wavelet regularization, as discussed above, allows preserving the image details while smoothing reconstruction artifacts but may introduce some irregularity in homogeneous area of the image. On the other hand, other (non wavelet- based) regularizatlon schemes are known to be adapted to regularize smooth regions, but at the expense of an oversmoothing of image details. This is particular the case of "Total Variation" TV regularizatlon, described inter alia by [Raj et al., 2007] and [Block et al., 2007].

A first improvement of the method of [Chaari et al. 2008] and

[Chaari et al. 2009], constituting an aspect of the present invention, consists in combining WT and TV in a joint regularizatlon framework in order to take advantage of their different properties, allowing them to alleviate the drawbacks of each other. A joint Wavelet Transform-Total Variation (WTTV) regularizatlon can reduce to the optimization of the following penalized criterion: WT-TV = T* arg min WT-TV(C)

= T* arg min J _L (T* ) + K _{l t} 7p(0 + /¾||T*C IITV

< (15) where J _L and J _P have bee, defined above, κι > 0 and κ ₂ > 0 are regularizatlon parameters, and T ^* is the wavelet adjoint operator. The discrete total variation p e C ^{x A*} . .

norm of an image is given by:

v x y=i x= i (16) where for every , V·) is the horizontally smoothed gradient operator defined by i(0) = hp(y + l..r + 1) - p{y + 1) + / (]/ + l. a ^* ) - ρ(ν, ΐ))-.< _ν <γ,ί< _χ <χ _{Ι Λ} ^ 2 - - - - (17) For the sake of presentation, it has been assumed that the image p is periodic or equivalently that the image boundaries are toroidal.

Minimizing the optimality criterion ( 5) is much more difficult than minimizing the space-only criteria given by (1 1 ) - (13), because more than two terms are involved, while the Forward-Backward (FB) algorithm discussed above only applies to the minimization of criteria comprising two terms. Iterative minimization of a non-differentiable convex function comprising more than two terms can be performed using the so-called Parallel ProXimal Algorithm (PPXA) [Combettes and Pesquet, 2008], which also requires calculating the proximity operator of each of three involved terms.

The difficulty stems from the calculation of the proximity operator of . To circumvent this difficulty, the TV penalization in equation (15) can be split into four terms as [Combettes and Pesquet, 2008]. The TV penalty term therefore reads:

3

, , (ff, r)€{0, l} ^a

where for every

Y/'J Λ/2

ft¾a+r<p) = Z v ^/ l(Vjp)(2y + q.2x + τ + \(V _x p ^T Y(2y + g, 2x + r)|-' i

For every q and r in {0, 1}, let ^r be the decimation operator defined by

I . fW x2X , ^YxX V = (ί „.3.)ι<„<2ν',ΐ<αΕ<2 ^,→ (^2v-Hi.2r+r ) ί<ν<2ΥΛ <*<2X · and U _q+2r be the following operator: U _q +2r : C ^y - ^V → C ^Fx

₀ p (f) X Λ

where for every ^μ , the operators V ₀ and V ₂ are defined by

V ₀ (p) = ^ + 1. Λ· + 1) + + 1) + p<y + 1. x) ÷ ρ($. χ))ι< _ν≤ γ _Λ≤9 <χ and

1

V ₂ (p) = -(p(y + 1. x + 1) - p(y. x + 1 ) - p(y + t, x) + p(t ,x)) i<y<v,i<i<x

Let also h be the function defined on Υβ λ ^' /2

h(p) =∑∑ > l(Vip)(|f _f Ϊ + β + |(V ₁ ^T ) ^T _{(?/ +} > ₇ 2,Ϊ)| ²

B=l χ=1

It turns out from the equations above that

Vi e (0, 1, 2,3}, tVi = h Consequently, the optimization problem in equation (15) can be rewritten as:

T ^* ∞gmm[J _L (T ) + «IJP(C) + κ ₂ fVi(T*<)]

ζ to

3

T ^* argmm[J _L (T ) + m J _P (Q + ί¾ΓΛ(¾(ΓΟ)

The regularized WT-TV reconstruction approach is summarized in Algorithm 2 where the PPXA algorithm [Combettes and Pesquet, 2008] is used to minimize the optimality criterion in equation (18) which is made up of J = 6 convex functions.

Algorithm 2

ch that∑Lo-* = « _j ⁿ⁾ )i<i<e e (C ^K )® for every i e {1.....6}. Set also ε > 0,

1: repeat =

11: ( ^(n÷1 - =( ⁿ⁾ + A _n (Pl"> - C- ⁿ >).

12: Compute J< ⁿ⁺¹ > = JwT-rviC*" ^f '-')-

13: Ji - il - 1.

14; until IJ< ⁿ > -J< ⁿ i>| < sj^ i

is: return = T ' ' In the method of [Chaari et al. 2008] and [Chaari et al. 2009], three-dimensional (3D) imaging is performed by stacking regularized two- dimensional (2D) images. A further improvement of the method, constituting another aspect of the present invention, consists in performing direct, regularized 3D image reconstruction. In the 3D case, the image reconstruction problem can still be written: d(r) = S(r)p(r) + n(r) like in equation (12); however, r = (y, x, z) is now the three-dimensional spatial position, ze{1 , . . . ,Z} being the position along the third direction (slice selection). Moreover, the penalty term Jp will depend on the distribution of the wavelet coefficient of a 3D dyadic wavelet transform.

A further improvement of the method, constituting another aspect of the present invention, consists in performing four-dimension (4D), spatio-temporal regularization in dynamic MRI. This embodiment of the invention can be applied e.g. to fMRI of the brain, wherein the whole brain volume has to be acquired several times, yielding a 4D dataset. In conventional fMRI, the 3D images are supposed independent although they belong to the same fMRI session. However, in practice, the 3D temporal images are somehow dependent since they belong to the same fMRI session involving the same experimental paradigm. The BOLD (Blood Oxygen Level Dependent) time-series and the acquisition noise are in fact correlated in time in a brain fMRI session. For this reason, taking into account temporal dependencies between 3D images helps to increase the SNR through the acquired images, and therefore enhances the reliability of the statistical analysis in fMRI. However, since in dynamic MRI the imaged object geometry generally changes during the acquisition, joining the reconstruction process to the temporal regularization is very difficult.

To deal with a 4D reconstruction of a time series of N _r three- dimensional images (N _r is usually even), the observation model of equations (12) and (19) will be rewritten as follows: d _t (r) = S(r)p _t {r) + n _t {r)

(20) where te{1 N _r } is the acquisition time. Using a dyadic 3D wavelet operator T, the coefficients will be reindexed so that = ( i.«ij)o®,i<j<j _{m with} ® = {0, 1} ³ \ {(0.0,0)}

Accounting for an additional temporal l _p regularization term, reconstruction of the 4D "volume" (i.e. of the time series of 3D images or volumes) is performed through the minimization of the following optimality criterion: ζ = argmm JsT«)

= argmjn J ∑ |R(H - 8{τ (Γ ^ι ){τ)\ . +∑ J _P ( ^l )

I I re{l _i ... _: V7 i}x{L....A ^' |x{l,... _! Z} t:r.A t=2

(21) where ζ= (ζ ¹ , ζ ² , . . . , ζ ^ΝΓ )Τ, κ>0 is a regularization parameter and J _P is defined as in (13).

In equation (21),

t=l r€{l....,V7R}x{l Xjx {1 Z] _{js the term} » representative of a likelihood of each reconstructed image, given ||Γ· '- T ^" 'US- corresponding acquired elementary images; ⁽ 2 js a temporal penalty term, representative of a pixel-by-pixel or voxel-by-voxel

∑ <C ^l )

difference between consecutive images of the series and ^{1 1} is a wavelet penalty term used for spatial regularization. It will be understood that a different (e.g. non wavelet-based) spatial regularization term could also be used. The adjoint wavelet operator T ^* is then applied to each component ζ' of ζ to obtain the reconstructed 3D image p ^l at the acquisition time t by taking into account the time dependencies with the other acquired images.

It will be noted that the temporal penalty term is not based on wavelets, but is representative of a pixel-by-pixel or voxel-by-voxel difference between consecutive images of the series. More precisely, in equation (21) the temporal penalty term is based on an L _p norm; the present inventors have found that p has to be greater or equal to one, and that it should preferably satisfy: 1<p<1.5. Other form of the temporal penalty term could be used, in particular based on (preferably convex) edge-preserving functions of a pixel- by-pixel or voxel-by-voxel difference between consecutive images of the series.

Minimizing the optimality criterion (21) is much more difficult than minimizing the space-only criteria given by (11) - (13), because more than two terms are involved, while the Forward-Backward (FB) algorithm discussed above only applies to the minimization of criteria comprising two terms. Iterative minimization of a non-differentiable convex function comprising more than two terms can be performed using the already-cited Parallel ProXimal Algorithm (PPXA) [Combettes and Pesquet, 2008], which also requires calculating the proximity operator of each of three involved terms. This task is quite simple for the two first terms of equation (21) since they are separable with respect to the time variable t and the spatial position r. However, this is not the case for the time penalization term (third term in Eq. (21), which makes non-trivial the calculation of the corresponding proximity operator. It is then proposed to rewrite the optimality criteria in J _S T by decomposing the time penalization into two terms which are separable with respect to the time variable t (a given acquisition time t is either involved in the first or in the second term), and for which the proximity operators are easy to calculate. More specifically, the temporal penalty term J _T is expressed as the sum of a first partial temporal penalty term J _T ¹ and a second partial temporal penalty term Jj , wherein: the first partial temporal penalty term is representative of pixel-by-pixel or voxel-by-voxel differences between each even-numbered image of the series and a preceding odd-numbered image; and the second partial temporal penalty term is representative of pixel-by-pixel or voxel-by-voxel differences between each odd-numbered image of the series and a preceding even-numbered image:

+ Ji( ) + JfiO

iV _r /2- l

JT<0 = « ∑ ||7-< ^{¾+ 1} - T ^* < ^w ||5

1=1 (22)

Since Jj ¹ and J7- ² are separable with respect to the time variable t, the corresponding proximity operator can easily be calculated based on the proximity operator of each of the involved terms. Let us consider the following function

: C x C* ^' — R.

(C ^f ,C* ^"l > κ|| ·ζ ^£ - Ί~ζ'- \% = φ ο Η( {ζ ^ι χ<- ¹ ))

^{1 'ρ} and H is a linear operator defined by

H : C* x C ^K — € ^K

( a , b) t- a— b

Its associated adjoint operator H ^* is therefore given by

H ^* :€ ^K — € ^K x€ ^K

a s→ (a.—a).

The proximity operator of Φ is then given by: ρΓθΧψ = prox, _¾%/ = Id + i-J ^* (prox< _¼ - Id) o II

1 (23) The resulting algorithm (Algorithm 3) for the minimization of the space-time optimality criterion JST is given below.

Algorithm 3

l: repeat ,

12: if f 1 then

13: Set P ^:f~ - "> = EL ^ ^r,) -

14: end if

22: en or

23: C ⁱⁿ⁺¹⁾ = << ⁿ⁾ + λ„( ί ^{η )} - ^inS ).

24: Compute J-<"+ ^l = JsT« ^(n+l) ).

25: 91 4- + 1 .

26: until < cj ^{(n", }} .

27: Set < = < ^{ " _' .

28: return pi = Τ*ζ' for every ί e { 1 N _r }.

It will be understood that space-time regularization can also be applied to the case of a time series of bi-dimensional images, even if this is less common. Moreover, the spatial penalty term may also include a total variation component and/or constraints, as in the case of purely spatial 2D or 3D regularization. It is worth noticing that in a 4D fMRI dataset, and since the imaged body volume is acquired many times, said body (i.e. a brain) may slightly move between successive scans. The inherent motion artifacts may be dramatic for the fMRI analysis. In this context, and since a standard fMRI study involves a motion correction step through applying a rigid transformation, another extension of the estimation method is proposed. This method accounts for the eventual motion artifacts by involving the rigid transformation parameters 6r in the estimation step. The transformation is first applied to the standard SENSE data which will be used for the estimation step. An extension of the above described methods is also proposed performing the estimation step not on the standard SENSE reconstructed data, but on a quadratically regularized (QR) version of SENSE using Tikhonov regularization with a smoothing matrix ( i.e. QR-SENSE). The QR- SENSE appears more efficient for improving the robustness of the estimation step than introducing an additional filtering step on the SENSE solution and remains computationally efficient since it admits a closed form expression. In order to further consider the motion artifacts occurring in fMRI, the proposed 4D reconstruction method (Algorithm 3) can also be extended to account for the rigid transformation parameters during the reconstruction process. The optimality criterion in equation (22) is therefore rewritten as follows :

JST«)

(22bis) where after applying the motion correction transformation and assuming that Ν _Γ /2

£=1 Γ6{1,...,ν7Λ} ^χ {1,... } ^χ {1, .-,Ζ}

and

= ∑ «r |(r ⁺¹ )(r) - (r<f)(r)|' f =1 re{ \,...,Y/R) x {l,... ,X } x {l,...,Z)

Since ¹ (resp. ^τ ) is separable with respect to the time variable t, its proximity operator can easily be calculated based on the proximity operator of each of the involved terms in the sums above.

The rigid transformation parameters 5r can be found by minimizing JST(0 , expressed by equation (22bis), also with respect to these parameters. Alternatively, said rigid transformation parameters can be estimated separately, on the basis of a conventional SENSE reconstruction.

All the regularization methods which have been described here (wavelet-based 2D spatial regularization, wavelet-based 3D spatial regularization, hybrid wavelet-based/total variation 2D and 3D spatial regularization, 3D and 4D space-time regularization) involve numerical parameters which have to be set, either manually or automatically using a suitable estimation algorithm. These parameters include the prior parameters, of the statistical distribution of the wavelet coefficients, or "hyperparameters": Moj ^,r « ² and 0 _o, j=(a _o,j ^Re , a _0,j ^lm )e (¾ ₊ ) ² . They can also include the temporal regularization parameter κ, the so-called "shape parameter", i.e. the value of "p" defining the L _p norm || || _p used in temporal regularization (see equation 21), and the parameters Κ ₂ determining the relative weight of the wavelet and total variation penalty terms in hybrid regularization.

The already-cited documents [Chaari et al.2008] and [Chaari et al. 2009] do not teach how to determine the relevant parameters. A further improvement of the method described by these documents, constituting another aspect of the present invention, consists in providing automatic calibration algorithms, allowing a fully or partly non-supervised wavelet-based regularization. Two different frameworks are considered to achieve this goal depending on the type of decomposition. In case of orthonormal wavelet decomposition, the hyperparameters can be accurately estimated using a maximum likelihood procedure (see below). In contrast, when resorting to redundant frame decomposition, a stochastic sampling procedure seems more efficient. In the latter case, the hyperparameters are estimated using the Minimum Mean Square Error (MMSE) or equivalently the posterior mean estimator (see below). The posterior mean estimator can also be used with an orthonormal wavelet decomposition, as this is a special case of frame.

First of all, the determination of the hyper-parameters defining the prior statistical distribution of the wavelet coefficients used for 2D or 3D spatial regularization will be considered.

Θ = A first possibility for determining the hyper-parameters

(μ, σ, (« , ¾) _oe0)1 ^^ ) _{js t0 maxjmjze the} integrated likelihood according to the following equation:

0 = ar ma (d: 0) = argmax / p(d | Τ ^* ζ) /(£; Θ) άζ,

θ θ J

(24) Maximizing Eq. (24) is a missing data problem since ζ is unknown. It requires integrating out the sought image decomposition ζ and iterating between image reconstruction and hyper-parameter estimation using the intensive EM algorithm described by [Dempster et al., 1977]. In order to alleviate the computational burden, it is advantageous to proceed differently by assuming that a reference or "auxiliary" full FOV image ^ is available, and so is its wavelet decomposition ^ ^— ^ ^ P. In practice, the auxiliary image ^ can be obtained using 1 D-SENSE reconstruction at the same R value, either un regularized or quadratically-regularized in a conventional way (e.g. using Tikhonov regularization) in the image space or in the k-space.

Then, the maximum likelihood (ML) estimation procedure consists of assuming that this auxiliary image is a realization of the full prior distribution and thus in fitting Θ directly on it:

In statistics, this solution is referred to as the complete data maximum likelihood as opposed to the above mentioned missing data maximum likelihood estimator. This procedure can be decomposed in two independent steps, the first one involving the setting of the Gaussian prior parameters (μ,σ) attached to the approximation coefficients and the second one being related to the estimation of the GGL prior parameters

( ^α α βο, ο£® < <ϊτ ^η * _χ f _{rom tne} corresponding detail coefficients

( O.j ) OGO , 1 < J < Jfmax _ On the one hand, ML estimators ^ ^t " ^{> σ} ^ are explicitly given by the empirical mean and standard deviation:

For each resolution level j and orientation o, and

O°.J} a arroe o ecsttiimmaa+toerdl f frrnomm ^ ^05, ^ as follows: (5¾,0¾) = argm^ /(Re« _oJ ); α, 9)

= arg max log / (R ε(ζ _{0 j fc} ) ; a, 0)

arg max ^Q |R*(Coj.Jt) R«(C ₀ ,j.jt)

2 3 Io _g ( ) ₊ ¾ lo _g (erfc(- -)) }

(26)

- , . . ¾ and ¾™ , , .

The hyperparameters °>·? °^ are estimated in the same way, by replacing Re( ) by Im(-) in equation (26).

This two-dimensional minimization problem does not admit a closed form solution. Hence, the ML parameters are computed using a numerical optimization method, such as direct search method (eg, Rosenbrock's method, see [Bertsekas, 1995, chapl , p159-165]). Alternative solutions based on Monte Carlo methods or the Stein principle can also be thought of, at the expense of an increased computational burden.

When temporal regularization is considered, two additional parameters have to be determined: the temporal regularization parameter κ and the shape parameter p (see equations 21 and 22).

Indeed, in a subject-level analysis, only data recorded for one subject is processed. For this reason, manually fixing the regularization parameter κ and the shape parameter p may be suitable regarding to the reasonable amount of data, even if this makes the proposed method not completely automatic and remains somehow user-dependent. When a group- level analysis is considered, the processed data becomes much larger, and since many subjects are considered, setting the temporal regularization parameter manually for all subjects at once is definitely sub-optimal due to the between-subject variability. This variability may be significant either from an anatomical or functional viewpoint. Manually setting this parameter for each subject apart would be possible, but cumbersome and it would make the proposed method more user-dependent regarding to the involved number of subjects in a group-level analysis (typically 15 subjects involved).

Even in a group-level analysis, the estimation is performed for each subject apart. Moreover, a parameter value is preferably estimated for each voxel of the 3D volume, which means that each voxel is processed independently of the others. Nonetheless, for the sake of computational efficiency, the same shape parameter p is considered for the whole volume. Otherwise stated, instead of determining a unique temporal regularization parameter, κ, a respective parameter κ _Γ is determined for each voxel.

Again, a maximum likelihood estimation of said parameters is performed. Let Hr) ... p ^Nr r) be the 3D images forming the 4D dataset, and let define V _r as:

V _r = (p ² (r) - p r), p r) - p ² (r) p ^2Nr (r) - p^ ))

Then, the joint estimation of the parameters p and κ _Γ in the ML sense is performed by minimizing the following criterion:

L p, oc -log(K _r ) + K _r (27) preferably under the constraint p>1.

Obviously, a given voxel is supposed to more likely behave similarly to its neighbors, as well as for anatomical and functional perspective. An alternative method consists of estimating a parameter value per functional region of the brain (i.e. a group of neighboring voxels having similar functional roles). This estimation is based on an a priori classification step of the brain in order to identify the different functional areas. This extension makes the estimation process more robust since more samples are available at a given estimation step. Indeed, the temporal signals related to all the voxels belonging to the same functional region are gathered in the same vector. If we denote by S, the size of a functional region i, the Maximum Likelihood estimation will be made on a signal involving S, times more samples. The resulting estimation is therefore more robust to atypical observations (outliers), which may be occur in an fMRI session. The image reconstruction method of [Chaari et al. 2008] and [Chaari et al. 2009] is based on the decomposition of the observed volume on a wavelet basis. A further improvement, constituting another aspect of the present invention, consists in generalizing said method to the case where decomposition is performed on a frame.

Frames generalize the concepts of wavelets, curvelets, bandlets, etc. A frame is a family of functions on which a signal (or image, which is simply a multi-dimensional signal) can be decomposed. In general, the resulting representation of the signal is redundant; a frame is called a basis when it is not redundant. Frames are generally used to help capturing more geometrical details of images under investigation. For this reason, using frames in pMRI reconstruction is helpful in order to track reconstruction artefacts'.

More precisely, let us consider real-valued digital signals of length L as elements of the Euclidean space Unendowed with the usual scalar product and norm denoted as <· I ·> and 11- 11, respectively.

Let K be an integer greater than or equal to L. A family of vectors (^ ^* ) ΐ<Α:<Λ ^' _{jn tne} fj _n jte-dimensional space R ^K is a frame when there exists a constant μ in ]0, +∞[ such that

K

(28)

If the inequality (28) becomes an equality, is called a i&k ^l <k<K _t jg^ _t f _rame The bounded linear frame analysis operator F and the adjoint synthesis frame operator F ^* are defined as

F M ^L →E ^K : «y|ej »i<Jt<Ar

A' r=l

(29) Note that F is injective whereas F* is surjective. When P =F ^* , then ^ k/l<k<K _{js an or} †h ₀ normal basis. A simple example of a redundant frame is the union of orthonormal bases. In this case, the frame is tight with μ=Μ and thus F ^* F=M\, where I is the identity operator.

In the following, two frame estimation problems of increasing complexity will be addressed, the first one corresponding to a denoising formulation and the second one to an extension involving a degradation operator such as the sensitivity matrix S in parallel MRI reconstruction. The denoising problem is addressed first and its application to pMRI matches the situation in which a noisy reference image, denoted as y in what follows, (e.g. such as the SENSE reconstruction) is available to estimate the hyperparameters of the frame representation. The second problem will be discussed after presenting Algorithm 4 and would correspond to estimating the hyperparameters of the frame decomposition directly from the reduced FOV images in the pMRI context.

„ _ J?* < _j > _L fi

An observed signal ^& can be written according to its frame representation (FR) involving coefficients 2 ^» as follows:

y = F*x + n ₍₃₀₎ where n is the error between the observed signal and its FR F x. This error is modeled by imposing that belongs to the closed convex set

C ₆ = {x€ U ^K I N(y - F*x) < 6} where ^ ^ [ i Οθ[ j _{s some} error bound and N can be any norm on M ^L .

In signal/image recovery problems, n is an additive noise that corrupts the measured data. The following developments will be focused on the case of a bounded observation error modeled by uniform noise. By adopting a probabilistic approach, y and x are assumed to be realizations of random vectors Y and X. In this context, the goal is to characterize the probability distribution of X|Y by considering some parametric probabilistic model and by estimating the associated hyperparameters. The estimation problem is much more difficult than in the case of the decomposition on a non- redundant basis, e.g. a wavelet basis, because F is not bijective.

In a Bayesian framework, it is necessary to define prior distributions for the frame coefficients. For instance, this prior may be chosen so as to promote the sparsity of the representation. In the following f(x,0), denotes the probability density function (pdf) of the frame coefficients that depends on an unknown hyperparameter vector Θ and†(0), is the a priori pdf for the hyperparameter vector Θ. In compliance with the constraint Ο _δ , n is assumed to be uniformly distributed on the ball

Βο, _δ — {α ζ. R ^L I N(a) < 6} _{and the unifom pdf on} the closed convex ball:

Bp-xfi = {y R ^L \ N{y - F*x) < 6} The hyperparameter vector 0 is also considered a realization of the random vector Θ, and the conditional pdf of (Χ,Θ) given Yean be written as:

For the sake of simplicity, it will be assumed that frame coefficients are a priori independent with marginal GG (Generalized Gaussian) distributions even if, as discussed above, a GGL distribution is better suited for applications to pMRI This leads to the following frame coefficient prior:

where ^L ¾ ^> ^ " (with ^{7 7} > ) are the scale and shape parameters associated with X _k , which is the component of the frame coefficient vector and Γ(·) is the Gamma function.

By introducing ^{ft k} , the frame prior can be rewritten as:

The distribution of a frame coefficient generally differs from one coefficient to another. However, some frame coefficients can have very similar distributions, which can be defined by the same hyperparameters. As a consequence, it is proposed to split the frame coefficients into G different roups. The g ^th group will be parameterized by a unique hyperparameter prior can be

(34) where the summation covers the index set S _g of the elements of the g ^th group containing n _g elements and e _g -=-(9^..e _G ). E.g., each group can correspond to a given wavelet subband. A coarser classification may be mad when using multiscale frame representations by considering that all the frame coefficients at a given resolution level belong to a same group.

The hierarchical Bayesian model for the frame decomposition is completed by the following improper hyperprior: where 1 _Α (ξ) is the function defined on Ac lT¾ by 1 _Α (ξ)=1 if ξεΑ and 1 _Α (ξ)=0 otherwise.

The motivations for using this kind of prior are the following

• the interval [0,3] covers all possible values of β _g encountered in practical applications, and there is no additional information about the parameter .

• The prior for the parameter γ ₉ is a Jeffrey's distribution that reflects the absence of knowledge about this parameter.

The resulting posterior distribution is therefore given by:

(36) The Bayesian estimators [e.g., the maximum a posteriori

(MAP) or minimum mean square error (MMSE) estimators also known as the posterior mean estimator] associated with the posterior distribution (35) have no simple closed-form expression.

In what follows, a stochastic procedure is proposed to compute the MMSE estimator, which relies on hybrid Markov Chain Monte Carlo (MCMC) algorithms. The idea is to use a suitable algorithm to generate samples distributed according to the posterior distribution (35). After convergence, the generated samples are used to compute the MMSE estimates of the unknown model parameter and hyperparameter vectors x and Θ, respectively. In the pMRI context, x is the sought frame representation of the reference image.

When the considered frame is the union of orthonormal M bases and N() is the Euclidean norm, the well-known Gibbs sampler (GS) can be used [Geman, 1984], which iteratively generates samples distributed according to conditional distributions associated with the target distribution. More precisely, the basic GS iteratively generates samples distributed according to f(x\0,y) and f(G x, y) to simulate realizations of the full posterior. f(x, 0 | y).

Straightforward calculations yield the following conditional distribution:

This conditional distribution is a product of GG distributions truncated on C ₅ . Actually, sampling according to this truncated distribution is not always easy to perform since the adjoint frame operator F_ is usually of large dimension. However, two alternative sampling strategies are detailed in what follows.

Naive sampling proceeds by sampling according to independent GG distributions

and then accepting the proposed candidate x only if ^ ^ F x) < S This method can be used for any frame decomposition and any norm. However, it can be quite inefficient because of a very low acceptance ratio, especially when δ takes small values.

The Gibbs sampler is designed to sample more efficiently from the conditional distribution in Eq. (36) when the considered frame is the union of M orthonormal bases and N( ) is the Euclidean norm. In this case, the analysis frame operator and the corresponding adjoint can be written as respectively, where Vme{1 , . . . ,M}, F _m is the decomposition operator onto the m ^th orthonormal basis such as F ^* _m F _m =

F _m F _m = Id. In what follows, every ^{LtJ m} with K=ML is decomposed as ^{X =} · ^{* * ?} Λ/ _wh ere ^x m€ 1 _{VmG{1 M}.}

The GS for the generation of frame coefficients draws vectors according to the conditional distribution f{ ^x n \ -n , %/·> $) under the constraint^ ^— ^ ^x ) — ^ where x-„ is the reduced size vector of dimension ^L built from x by removing the n ^th vector x _n . If N( ) is the Euclidean norm, then Vne{1 , . . . ,M}:

M

N (y -∑ FZ _lXrn ) < 6

m= 1

M

^ II f *(F _n y - T F _n F^x _m ) \\ < δ

= 1

II F _n y -∑ F _n F _n x _m - x _n ||< δ (since Vz e M ^L , || F*z \\= \\ z ||)

(38) n = F _n y ^ ^ F _n F _m x _m

where m≠n

To sample each x _n , it is proposed to use an MH step whose proposal distribution is supported on the ball > defined by:

B _Cn = {a e R ^L I iV(a - c ) < 5} _(3g)

Random generation from a pdf q ₅ defined on will now be discussed.

First of all, it will be considered how to sample vectors in the unit ? _p ball (p e]0,+∞]) of M ^L .ln the special case p=+∞, this can be easily performed by sampling independently along each space coordinate according to a distribution on the interval [-1 , 1]. The parameter "p" considered here should not be confused with the shape parameter used in temporal regularization, see equations 21 and 22.

Thus, this section focuses on the more difficult problem associated with a finite value of p. In the following, || · || _p denotes the t _p norm.

Let A = [Ai, . . . ,A _L ] ^T be the random vector of i.i.d. components which have the following GG(p ^1/p , p) pdf:

Let U = [Ui, . . . ,Uf ] ^T = A/||A|| _P . Then, it can be shown ([Gupta et al., 1997]) that the random vector U is uniformly distributed on the surface of the f _p unit sphere of K ^L and the joint pdf of U-i, . . . ,U _L -i is

wh ere D _p = {(«! , ..., « _£ ,_,)€ R ^L ' ¹ 1 < 1} .

T / '

The uniform distribution on the unit £ _p sphere of R ^' will be denoted by ^u ( ^L ', p) . Therefore U = \U ¾,] ^T - {L p).

For every e {1 , . . . ,L'-1}, it can be shown ([Gupta et al., 1997]) that the pdf of V = [Ui, U _L ] ^T is given by:

— ^{UL) 1}

(40)

In particular, if P ^{fc 1} · ^* and L'=L+p, equation (40) provides the uniform distribution on the unit t _p sphere of * ^" . Sampling from any distribution q _n on the t _p ball of radius ^ ^ Ή ^ is straightforwardly deduced by scaling V. Having a closed form expression of this pdf is important to be able to calculate the acceptance ratio of the MH move. To take into account the value of x _n ⁽ ' ^~1) obtained at the previous iteration (7 - 1), it may however be preferable to choose a proposal distribution supported on a restricted ball of radius e ]0, δ[ containing x„ ° ^~1> . This strategy, similar to the random walk MH, results in a better exploration of regions associated with large values of the conditional distribution More precisely, it is proposed to

J3* (t- i ) .

choose a proposal distribution defined on ^mn M where and P is the projection onto the ball ^Ο,δ-η defined as:

This choice of the center of the ball guarantees that

¾(*- *)„ ^c ¾ ,i

Moreover, any point of B ^c n≠ can be reached after consecutive draws in ^x >*7. The radius η has to be adjusted to ensure a good exploration of B «' ^d . In practice, it may also be interesting to fix a small enough value of η (compared with δ) so as to improve the acceptance ratio.

Instead of sampling the hyper-parameter vector Θ according to ί(θ\ ^χ , y) j _t j _s proposed to iteratively sample according to I¾ > ²¹ >^ and 9 * * s ' ^*■ ' J straightforward calculations yield the following results:

Consequently, -^ ^j ?| ¾ ^{? '■} is the pdf of the inverse gamma distribution ^{9 J} that is easy to sample.

Conversely, it is more difficult to sample according to the truncated pdf fiPg g ^{i x i} y) j _n j _s j _s achieved by using an MH move whose proposal

A. g I ng ^} j _{s a} Gaussian distribution truncated on the interval [0, 3] with standard deviation ^σ βα = 0.05 [Dobigeon et al, 2007], [Robert, 1995]. The mode of this distribution is the value of the parameter '' 9 at the previous iteration (i-1). The resulting method is the hybrid GS summarized in Algorithm 4.

If one wants to integrate altogether the reconstruction and hyperparameter estimation steps in the same Bayesian framework, observation model in Eq. 30 can be extended to y=SF ^* x+n, where S is the sensitivity matrix defined in Eq. (3). The inference has then to be reconducted accordingly.

Algorithm 4

l: Initialize with some θ ^{0) = {B ) _{x 9≤(:} = f* ^)\< ₉ <o ami ar"" e C¾.

and set i— 1.

..: repeat

:¾·. Sampling a?:

!: for rt - 1 to M do

Compute ci,° = ^' „(ϊ/ ∑,,„>„ ^F ', '< ")

and x„ = i {x„ „ ; f c„ .

<>: Simulate x„ ^rl na follows:

• Generat j,° ~ !/,)!,, - xj,' where </,, i-s defined on Z¾

• Com ute the ratio

and accept the proposed candidate with the probability min{ 1 , r( x '„''', xi' ¹ ' )}.

7: end for

* _: Sampling Θ:

!>: for <y = 1 to (V do

1 ^; Generale ^" ° - 10 ( _· ^,. \*ΪΥ ^: ' ^' ")·

] ] : Simulate us fulUnvs:

• General e f,, \ J.

• Compute t he rat io

,. η I) . _ -^Λ hj■ *· vh ,. > )

^' » " ^' '-l ' ,·. ,1» I I I» .- . ·. ,7,1)) ! »i> I t-, and accept the proposed candidate with the probability ininf 1. r( ·¾", · ' ' ' )}·

3 : end for

y.i: Set ?— i + i.

3 I: until Convergence

Although this algorithm is intuitive and simple to implement, it must be pointed out that it was derived under the restrictive assumption that the considered frame is the union of M orthonormal bases. When these assumptions do not hold, another algorithm, discussed below, allows sampling frame coefficients and the related hyperparameters by exploiting algebraic properties of frames. Indeed, as a direct generation of samples according to f(x\&, y) is generally impossible, an alternative method is proposed which replaces the Gibbs move by an MH move. This MH move aims at sampling globally a candidate x according to a proposal distribution. This candidate is accepted or rejected with the standard MH acceptance ratio. The efficiency of the MH move strongly depends on the choice of the proposal distribution for x.

We denote as x ⁽ⁱ⁾ be the i ^th accepted sample of the algorithm and q(x \ x°~ ¹⁾ ) the proposal that is used to generate a candidate at iteration / ^' . The main difficulty for choosing q(x | x° ^~1) ) stems from the fact that it must guarantee that x e C ₅ while yielding a tractable expression of q(x ^<M) I x)/q(x \ x ^{( 1)} ).

For this reason, it is proposed to exploit the algebraic properties of frame representations. More precisely, any frame coefficient vector can be decomposed as x = XH + XHX, where XH and x _H ± are realizations of random vectors taking their values in H = Ran(F) and H ^x , = [Ran(F)] ^x = Null(F*), respectively. It is recalled that the range of F is RanfF) = R ^K i/ » ^L . F _S , = x) _{and tne} null space of F* is

NuilCF" ) = {x e H ^K |F"x = 0}

The proposal distribution used here allows generating samples x _H e H and XH±, e H ¹ . More precisely, the following separable form of the proposal pdf will be considered: g(x I x ^l ) = q (x _H

where ^{X/i IJ} > ^X H ^{E W} _AND ^ ~ H + ¾I . | _n other words, independent sampling of XH and xniwill be performed.

If the decomposition X=XH + XHX is considered, sampling x in

C ₅ is equivalent to sampling λ e where - {A e P lNd/ - F ^* FX) < 6} , _{ndeedi h = Ρλ where} λ e M^and, since ^'W ^{1 e} N !l( ^j P ^* ) ^, i- x ^* = ^* F X Sampling λ in _ca n be easily achieved, e.g., by generating u from a distribution on the ball f< ^d and by taking λ = (F*F) ^"1 w. To make the sampling of XH at iteration / more efficient, taking into account the sampled value at the previous iteration

=™<'-» - F(F-F)-W- _{may be interest|ng}

Similarly to random walk generation techniques, u is generated in ^D ^ ^~ - where η ε]0, δ[ and = Ρ{^ ^~1) - V) + V . This allows drawing a vector u such that ^x n = (f F) ' ^{u e} ^ and

N(u - u ^{{i "} ) < 2η

The generation of u can then be performed according to equation (40) provided that N( ) is an t _p norm with p€ [1 ,+∞].

Once that ^XH = FX e H n C$ ( _wn j _c |-, ensures that x is in

C ₅ ), has been simulated, x _H ± has to be sampled as an element of H ¹ . Since If = F*x 4- n = F ^* xii + n _t there is no information in y about χ _Η χ. As a consequence, it is proposed to sample XH by drawing z according to the Gaussian distribution and by projecting z onto H ¹ , i.e.,

tt _H ± = I— F(F ^* F) ^{~ l} F ^* j _s the orthogonal projection operator onto H ¹

It can then be shown that the expression of the proposal pdf q(x ^{i' - ^] 1 ¾) - P(u - y) - y) This expression remains valid in the degenerate case when K

= L (yielding x _H ±= 0). Finally, it is important to note that, if q _n is the uniform distribution on the ball B ₀ , _η , the above ratio reduces to 1 , which simplifies the computation of the MH acceptance ratio. The final algorithm is summarized in Algorithm 5. It is noteworthy that the sampling of the hyper-parameter vector is performed as for the hybrid GS of Algorithm 4. Algorithm 5 l: Initialize with some Θ ^Μ = = 4° ⁵ )' ₉ <(; and ii ^l9 >€ B _{y S} .

Set x ⁽⁰⁾ = F{F ^* F)- ^l u.W and i = 1.

2 repeat

3: Sampling

Compute? ύ P(u ⁽ '-" -y) + y.

Generate u ⁰ ' wliM'c £/,, is defined on

Compute x ¹ , - F(F'F) 'tt ⁽ .

Generate «<·> ^. (x ⁽ ' ^{"_ l)} ,c»£l ^' ).

Compute = Π _Η . z' and = x ^l ÷ .

J& ^,) ≠ ⁱ - ^l .y) _q!! (ul-" ~ Ρ{ΰ '> y) y)

Compute the ratio Hi' ¹ ' , x ¹ ' ¹ ')

/<*<· "!0<'- ". y)

and accept the proposed candidates and with probability

min[l,r( ⁽ .-c ⁽ '- ^,) )}.

Sampling Θ:

for <)— 1 to (7 do

Generate - *Q (^ΤΪ ·Σ^ , ΐ Ύ ^Γ °)·

Simulate .^ί ' s follows:

Generate .¾" - ' ")

Compute the ratio

ami accept the proposed candidate with the probability iuin{ I . r{/ij _; ''. " ')} end for

Set i *- i 4 I .

until Convergence

The hierarchical Bayesian model of equations (31 - 35) can be extended to include an additional term in the prior depending on the Total Variation (TV) of the image to be reconstructed. Like redundant frame representations, using TV priors leads also to a hyper-parameter estimation problem.

Assuming an exponential shape, the new prior can be expressed as:

(47) where ^ ^{= 1} ' ' ^{" K} ^ is the new hyper-parameter vector with K > 0, H II _TV is the TV semi-norm and Ζ(θ) is a normalization constant. The new hierarchical Bayesian model for the frame decomposition is completed by the following improper hyperprior:

G G

/(fl) = Ζ(θ)φ{κ) f[ φ(θ ₉ ) = Ζ(Θ)ψ(κ) f] {ψ(Ύ _β )ψ(β ₉ )}

9=1 g=i

(48) where K _MAX is a positive real to be fixed (in a preferred embodiment, K _max =10).

It can be shown that Ζ(θ) is uniformly bounded with respect to Yg and therefore that the hyperprior f(0) in Eq. (47) has a stable asymptotic behaviour when γ ₉ →+∞.

The resulting new posterior distribution is therefore given by: x exp (-K||F* ||Tv) 1[0, _Κπι „](κ) (49).

The Bayesian estimators associated with the posterior distribution in Eq. (48) still have no simple closed-form expression. For this reason, the same sampling strategy discussed above will be applied, and the frame coefficients will be sampled as in Algorithm 5. However, for the hyper- parameters vector, straightforward calculations show that the posterior distribution for the hyper-parameters γ ₉ , β ₉ and κ will be expressed as: κ, as, y) oc g ^ ¹ ex ( -— |¾| ΐ*+ (¾)

^% *e¾ (50)

and

/(κ|7ΐ, ·· · ,7σ,£ι... · , ϊσ, ^χ ,!/) oc exp(- ||F ^* «i| _T v) («) ( ₅₂ ) respectively.

Consequently, ^^l^' ' ³ ^^ j _s the pdf of the inverse gamma , ^β . ^a . . . Sampling y _g will therefore be performed exactly as in Algorithm 4. Conversely, it is more difficult to sample according to /(Aba. and /( ^K l7t, · · ·, m - · ,0G,x,y). This task is achieved by using two MH moves whose proposal distributions l ⁾ ) and g(« I κ ^{(ι 1} ) _{gre Gaussjan} distributions truncated on the intervals [0, 3] and [0, Kmax] with standard deviations ap _g — 0.05 and σ _κ — 0.01, respectively. These standard deviation values have been fixed based on our practical observations. The resulting method to sample according to the posterior distribution in Eq. (49) is summarized in Algorithm 6.

Algorithm 6

1. Initialize with some 0 ^(o = ((0<, ^O) )i< ₉ c « ⁽⁰ >) = ((vf, Λ< _* <«. κ ⁽⁰⁾ )

and u< ^0> € Set * ⁽⁰⁾ - (F ^* / ^" i ^{_i ( )} and i = 1.

2: repeat

Sampling ¾:

• Compute ti ⁽ ' = (u ⁽ ' '> - y) + y.

» Generate u*'' -> — u ⁽ ' ^~ ") where </,, is defined on B _n _, _r

• Compute = F{F ^* F) ^~ u ' .

• Generate z ^{l ~ λ ^Γ (χ- ^ι) ,σ ΐ).

• Compute χ^', = Π _Η . and a: ^(l) = x ^' }' - x' , _. .

• Compute t he rat 10 r ( x ' ' , x ^{■ ■} ) = r— 7? ; :

»»|^'- ^ι) ,ν)¾(ΰ ⁽ " ώ ^{(ί u} ) and accept the proposed candidates w ¹ '' and ^(,¾ with probability ηηη{1,Γ(χ ^ί: ",¾: ^,, - ¹ ')}.

• Generate - ii ' ·

•

,,ί,,' ^{- 1} ')}. j<: end For

0: Simulate κ ⁽ '* as follows:

• Generate </(κ | ⁽ ' ⁿ ).

• Compute t e ratio

,,.;.,,·, _K ! , .... ) , _

and accept the proposed candidate with the probability min{ 1. r{ , to Set i ÷- # ~ 1.

II until Convergence

The technical results of different embodiments of the inventive methods will now be discussed with reference to figure 2 - 20.

The anatomical results reported in Figures 2-10 and 13 were obtained in the following conditions. Experiments have been conducted on real data sets comprising 256x256x14 Gradient-Echo (GE) anatomical with 0.93x0.93x8 mm3 spatial resolution. GE anatomical images were acquired with TE/TR = 10/500 ms and BW=31.25 kHz . Note also that these images have been acquired using acceleration factors R=2 and R=4 on a Signa 1.5 Tesla GE Healthcare scanner with an eight-channel head coil. Interestingly, the scanning time of anatomical data lasted 5~mn in non-parallel imaging, while acquisition duration was decreased to 3mn10s and 2mn20s in parallel imaging with R=2 and R=4, respectively.

Figures 2 and 3 show 9 anatomical slices of a human brain obtained using SENSE with Tikhonov regularization. The Tikhonov parameter has been fixed manually. The reduction factor is R=2 for figure 2, and R=4 for figure 3. Despite the regularization, some aliasing artifacts are visible in the form of curved lines, particularly on figure 3. Some oversmoothing of the images is also visible.

Figure 4 has been obtained, in the same experimental conditions, using TV regularization (left panel: R=2; right panel: R=4 - a single slice is shown) with a manually tuned regularization parameter κ. Quite strong aliasing artifacts and some "staircase" defects are visible; these defects can be mitigated by increasing the value of κ, at the expense of the informational content of the images.

Figures 5 and 6 have been obtained, in the same experimental conditions, using an autocalibrated 2D wavelet transform-based regularization scheme, as described above; this method will be called 2D- UWR-SENSE, where "UWR" stands for "Unconstrained Wavelet Regularization". Again, R=2 for figure 5 and R=4 for figure 6.

More precisely, dyadic (M = 2) Symmlet orthonormal wavelet bases associated with filters of length 8 were used over j _max = 3 resolution levels. Regarding the wavelet coefficients, the prior of equation (10) has been employed.

The related hyper-parameters (a couple of hyper-parameters is fitted for real/imaginary parts of each subband, i.e. each approximation/detail coefficients at each resolution level and orientation) were estimated using the Bayesian approach described above. Full FOV image reconstruction was then performed using wavelet-based regularization.

The smoothing effects observed in figures 2 and 3 with Tikhonov regularization no longer exist in the WT regularized images of figures 5 and 6, where a quite accurate reconstruction is performed within the brain mask without introducing the staircase effects observed with TV regularization (see figure 4).

A further improvement of the image quality can be obtained by incorporating an additional constraint in the method described hereabove in order to better regularize artifact regions. The ensuing algorithm is called CWR-SENSE, where "CWR" stands for "Constrained Wavelet Regularization".

This method implies imposing local lower and upper bounds on the image intensity values in artifact areas, regardless of their shape and/or location. These bounds define the nonempty closed convex set: C = {p€ C ^K I Vr e {1. . . - , Y/R} x {1, . . . . X}. p(r)€ C _r ] _(g3) where the constraint introduced on the range values at position re2 {1 , . . . , Y/R} * {1 , . . . ,X} is modeled by:

C = {ξ e C I Re( e lm( e ll ^m )

(54) jrRe _ fjRe [ Re 1 irlrn _ frlni rim

with L½iin,r÷ ^A nia ,rJ _anc | ½ — L½iin _t ' ¹ n>ax, ^j . The optimality criterion becomes then:

c ^* = (c e ^K I T ^* c C)

where ^{~ 1 1 "} and ^l c is the indicator function of the closed convex set C defined by: - f 0 if < C*

[ -oc otherwise. (55) A morphological gradient [Serra, 1982] can used to detect artifact regions (very low/high transitions in the gradient image) on which the additional convex constraints are applied. The upper and lower bounds that define the convex sets C _r are spatially varying since they depend on r. They can thus be computed using a morphological opening and closing operations applied to the basic-SENSE reconstructed image (auxiliary, or reference, image) in order to discard very low and high intensities.

Figures 7 and 8 illustrate the results of the CWR-SENSE method, for R=2 and R=4 respectively. It can be seen that the surviving artifacts in figures 5 and 6 have now been removed due to the anisotropic smoothing using the additional convex constraint.

From a quantitative point of view, significant improvements were achieved by the UWR/CWR-SENSE algorithms in comparison with basic-SENSE and Tikhonov reconstructions. The following table 1 reports the signal-noise-ratio (SNR) values in dB corresponding to the basic-SENSE, Tikhonov regularization and the proposed UWR/CWR-SENSE techniques for the illustrated slices of the anatomical brain volume shown in figures 3, 6 and 8 (R=4). On average, the gain drew from the proposed constrained regularization strategy amounts to 1.08 dB and a better visual quality.

Table 1

The influence of the choice of the wavelet basis has also been studied. More particularly, four different bases have been considered: dyadic Symmlet 8, dyadic Daubechies 8, dyadic Haar and Meyer with M = 4 bands [Chaux et al., 2006b]. The first three bases give quite similar results (apart from some blocking effects which only occur with the Haar basis), the dyadic Symmlet 8 leading to a slightly higher SNR. On the contrary, the Meyer 4- band wavelet basis leads to a significantly lower SNR.

Figure 9 shows the results obtained by combined wavelet-total variation regularization (algorithm 2) using a dyadic (M=2) Symmlet orthonormal wavelet basis associated with filters of length 8 over j _ma x=3 resolution levels, with a reduction factor R=4. Figure 10 show the results obtained by combined wavelet-total variation regularization using a redundant wavelet frame constituted by a Union of 2 Orthonormal Bases (U2OB) with Symmlet 8 and Symmlet 4 filters. The hyper-parameters of the wavelet prior and the TV regularization parameter were estimated using the approach described above. These figures show that the reconstructed images present better regularity than the ones reconstructed using the UWR-SENSE algorithm (figure 6). However, from a visual viewpoint, using redundant WT do not necessarily lead to better reconstruction quality. From a quantitative viewpoint, SNR values in dB are provided in the table 2 below. Comparisons with SNR values in table 1 confirm the usefulness of combining WT and TV in a joint regularization framework. This table shows that a slight improvement of 0.02 dB is obtained compared to the UWR-SENSE algorithm when using Symmlet wavelets. However, more significant improvement (0.24 dB) is reached when using the U2OB redundant wavelet frame. It turns out then that even if the reconstruction performance obtained using an orthonormal wavelet basis and a redundant wavelet frame seems equivalent, SNR values indicate that using redundant frames is fruitful.

Figure 13 shows the reconstructed images using TV regularization (left), 3D-UWR-SENSE with a wavelet frame constituted by a union of two orthonormal wavelet basis (middle) and a hybrid frame-TV regularization method (right). The top row shows whole images, the bottom row a magnified detail. The two orthonormal bases whose union constitutes the wavelet frame used for reconstruction (middle and right columns) are a Daubechies basis of length 4 and a shifted Daubechies basis of length 8. Three resolution levels have been used, which means that G = 20 groups of wavelet coefficients are considered. For the frame representations, the hyperparameters have been chosen using the hybrid MCMC algorithm described above.

Table 2

The anatomical results reported in Figures 11 and 12. were obtained in the following conditions. The anatomical MRI scan was performed on a 3T Tim Trio Siemens scanner using a 3D T1 -weighted MP-RAGE pulse sequence and a matrix array head coil consisting of 32 receive channels. The scan parameter were chosen as follows: slice orientation= sagittal (Right- >Left), slice thickness=1. mm; TE=2.98ms, TR=2300ms; Tl=900ms, Flip Angle: 9°, BW=61 kHz, single shot. The field of View was 256x240x176 mm3 with a matrix size of 256x240x160 corresponding to anisotropic resolution of 1x1x1.1 mm3. The scanning time was TA_conv=9min14s using conventional MRI i.e. without acceleration (neither pMRI nor partial Fourier). The use of a 6/8 partial Fourier scheme enabled to decrease the scanning time to 7min46s. The parallel MRI data were collected without partial Fourier using an acceleration factor R=2 or R=4. The respective scanning times for these exams were 5min03s and 2min59s and not TA_conv/R as expected. The reason for this lies in the Siemens k-space sampling strategy. The manufacturer adopts a full sampling scheme for the 24 central lines inducing an actual acceleration factor lower than the prescribed one (1.83 instead of 2 and 3.09 instead of 4).

Figure 11-12 illustrate different MRI reconstruction algorithms from single subject T1 -weighted MRI data acquired at 3 Tesla (Siemens Tim Trio) using a matrix array 32-channel head coil. Figure 11 and 12 show axial and coronal views, respectively for the same subject. Top row in these figures shows the ground truth i.e. the reconstruction performed from full k-space acquisition. Then, the mSENSE (middle row) and 3D-UWR-SENSE (bottom row) algorithms are compared in a noisy context ie for R=4. For the latter algorithm the results have been obtained for J=2 resolution levels and Daubechies wavelets (non-redundant). Also, the hyperparameters of the wavelet representation have been estimated using a complete-data maximum likelihood procedure, as detailed above in the manuscript. MSENSE reconstruction artefacts appearing as white matter spots on the grey matter's boundaries are indicated by white circles in both views outlines. Clearly, the 3D-UWR-SENSE algorithm outperforms the mSENSE technique since these artefacts do not appear at the same location. Also, the elliptical-shaped artefacts in the centre of the axial view are strongly filtered by the inventive wavelet-based algorithm. Coronal views confirm the substantial noise reduction in the brain stem and subcortical regions using the proposed invention.

The results obtained using 2D wavelet bases are not shown here, as they are qualitatively close to the ones reported in figures 11 and 12. However, the 3D wavelet-based reconstruction quantitatively outperforms the 2D-wavelet-based reconstruction in terms of Signal to Noise Ratio (SNR). Note that the 3D reconstruction generated a SNR improvement of 1.3 dB with respect to 2D technique.

Until here, only results regarding the case of static, anatomical imaging have been discussed. However, as discussed above, the present invention also applies to the space-time regularization of 4D fMRI image series, obtained e.g. by echoplanar imaging (EPI). The corresponding results are illustrated by figures 14 - 20. For validation purpose, fMRI data were acquired on a 3 T Siemens Trio magnet using a Gradient-Echo EPI (GE-EPI) sequence (TE = 30 ms, TR = 2400 ms, slice thickness = 3 mm, transversal orientation, FOV = 192 mm ² ) during a cognitive localizer [Pinel et al, 2007] protocol. This experiment has been designed to map auditory, visual and motor brain functions as well as higher cognitive tasks such as number processing and language comprehension (listening and reading). It consisted of a single session of N _r = 128 scans. The paradigm was a fast event-related design comprising sixty auditory, visual and motor stimuli, defined in ten experimental conditions (auditory and visual sentences, auditory and visual calculations, left/right auditory and visual clicks, horizontal and vertical checkerboards). An L = 32 channel coil was used to enable parallel imaging.

For each of 15 subjects, fMRI data were collected at the 2 x 2 mm ² spatial in-plane resolution using different reduction factors (R = 2 or R = 4). Based on the raw data files delivered by the scanner, reduced FOV EPI images were reconstructed using two specific treatments:

i) k-space regridding to account for the non-uniform k-space sampling during readout gradient ramp, which occurs in fast MRI sequences like GE-EPI;

ii) Nyquist ghosting correction to remove the odd-even echo inconsistencies during k-space acquisition of EPI images.

Then, the 4D-UWR-SENSE and the already discussed 3D- UWR-SENSE (simply called "UWR-SENSE" for short) algorithms have been utilized in a final step to reconstruct the full FOV EPI images and compared to the mSENSE solution (mSENSE is the unregularized SENSE algorithm implemented by Siemens scanners).

Figure 14 compares the two pMRI reconstruction algorithms to illustrate on axial, coronal and sagittal slices how the mSENSE reconstruction artifacts have been removed using the 4D-UWR-SENSE approach. The mSENSE reconstructed images actually present large artifacts located both at the center and boundaries of the brain in sensory and cognitive regions (temporal lobes, frontal and motor cortices); on the figure, the artifacts are outlined by ellipses. This results in SNR loss and thus may have a dramatic impact for activation detection in these brain regions.

Irrespective of the reconstruction pipeline, the full FOV fMRI images were then preprocessed using the SPM5 software (http://www.fil.ion.ucl.ac.uk/spm/software/spm5/): preprocessing involved realignment, correction for motion and differences in slice acquisition time, spatial normalization, and smoothing with an isotropic Gaussian kernel of 4mm full-width at half-maximum. Anatomical normalization to MNI space was performed by coregistration of the functional images with the anatomical T1 scan acquired with the thirty two channel-head coil. Parameters for the normalization to MNI space were estimated by normalizing this scan to the T1 MNI template provided by SPM5, and were subsequently applied to all functional images.

For performing subject-level analysis, a General Linear Model (GLM) was constructed to capture stimulus-related BOLD response. The design matrix relied on ten experimental conditions and thus made up of twenty-one regressors corresponding to stick functions convolved with the canonical Hemodynamic Response Function (HRF) and its first temporal derivative, the last regressor modeling the baseline. This GLM was then fitted to the same acquired images and reconstructed using the Siemens reconstructor, UWR-SENSE and 4D-UWR-SENSE.

Contrast estimate images for motor responses and higher cognitive functions (computation, language) were subjected to further analyses at the subject and group levels. These two contrasts are complementary since the expected activations lie in different brain regions and thus can be differentially corrupted by reconstruction artifacts.

More precisely, were studied:

• the Auditory computation vs. Auditory sentence (aC-aS) contrast which is supposed to elicit evoked activity in the frontal and parietal lobes, since solving mental arithmetic task involves working memory and more specifically the intra-parietal sulcus; • the Left click vs. Right click (Lc-Rc) contrast for which we expect evoked activity in the right motor cortex (precentral gyrus, middle frontal gyrus). Indeed, the Lc-Rc contrast defines a compound comparison which involves two motor stimuli which are presented either in the visual or auditory modality. This comparison aims therefore at detecting lateralization effect in the motor cortex.

These two contrasts were chosen because they summarized well different situations (large vs small activation clusters, distributed vs focal activation pattern, bilateral vs unilateral activity) that occurred for this paradigm when looking at sensory areas (visual, auditory, motor) or regions involved in higher cognitive functions (reading, calculation).

Figure 15 shows subject-level student-t maps superimposed to anatomical MRI for the aC-aS contrast. Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively, with R=2 (top of the figure) and R=4 (bottom of the figure). The neurological convention (left is left) is adopted. The cross shows the maximum activation peak.

For the most significant slice and R = 2, all pMRI reconstruction algorithms succeed in finding evoked activity in the left parietal and frontal cortices, more precisely in the inferior parietal lobule and middle frontal gyrus. However, for R = 4 only UWR-SENSE and 4D-UWR-SENSE - and preferentially the latter - enable to retrieve reliable frontal activity elicited by mental calculation, which is lost by the mSENSE algorithm. From a quantitative viewpoint, the proposed 4D-UWR-SENSE algorithm finds larger clusters whose local maxima are more significant than the ones obtained using mSENSE and UWR-SENSE, as reported in Table 3. Concerning the most significant cluster for R = 2, the peak positions remain stable whatever the reconstruction algorithm. However, examining their significance level, one can first measure the benefits of wavelet-based regularization when comparing UWR-SENSE with mSENSE results and then additional positive effects of temporal regularization and 3D wavelet decomposition when looking at the 4D-UWR- SENSE results. These benefits are also demonstrated for R = 4. Table 3 shows the significant statistical results at the subject-level for the aC-aS contrast (corrected for multiple comparisons at the significance level of a = 0.05, which means that the null hypothesis is rejected if the p-value is smaller or equal to a).

The 'p-value" of the statistical significance test should not be confused with the "shape parameter" used in temporal regularization, see equations 21 and 22, or with the parameter "p" used to illustrate the MCMC algorithm.

Table 3

Figure 16 illustrates between-subject variability of detected activation for the aC-aS contrast at R = 2. Indeed, when comparing subject- level student-t maps reconstructed using the different pipelines (R = 2), it can be observed that the mSENSE algorithm fails to detect any activation cluster in the expected regions for the second subject. In contrast, the 4D-UWR- SENSE method retrieves more coherent activity while not exactly at the same position as for the first subject.

Figure 17 shows subject-level student-t maps superimposed to anatomical MRI for the Lc-Rc contrast. Data have been reconstructed using the mSENSE, UWR-SENSE and 4D-UWR-SENSE, respectively.

It can be seen that all reconstruction methods enable to retrieve expected activation in the right precentral gyrus. However, when looking more carefully at the statistical results (see Table 4), the UWR-SENSE and more preferentially the 4D-UWR-SENSE algorithms retrieve an additional cluster in the right middle frontal gyrus. On data acquired with R = 4, the same Lc-Rc contrast elicits similar activations, i.e. in the same region. As it can be seen on the bottom of the figure, this activity is enhanced when pMRI reconstruction is performed according with the methods of the invention.

Quantitative results in Table 4 confirm numerically what can be observed in the figure: larger clusters with higher local t-scores are detected using the 4D-UWR-SENSE algorithm, both for R = 2 and R = 4. More precisely, table 4 shows significant statistical results at the subject-level for the Lc-Rc contrast (corrected for multiple comparisons at the significance level of a = 0.05, which means that the null hypothesis is rejected if the p-value is smaller or equal to a).

Table 4

Figure 18 reports on the robustness of the proposed pMRI pipeline to the between-subject variability for this motor contrast. Since sensory functions are expected to generate larger BOLD effects (higher SNR) and appears more stable, our comparison takes place at R = 4. Two subject- level student-t maps reconstructed using the different pMRI algorithms are compared. For the second subject, one can observe that the mSENSE algorithm fails to detect any activation cluster in the right motor cortex. In contrast, the 4D-UWR-SENSE method retrieves more coherent activity for this second subject in the expected region.

To summarize, on these two contrasts the 4D-UWR-SENSE algorithm always outperforms the alternative reconstruction methods in terms of statistical significance (number of clusters, cluster extent, peak values,...) but also in terms of robustness.

Due to between-subject anatomical and functional variability, group-level analysis is necessary in order to derive robust and reproducible conclusions at the population level. For this validation, random effect analyses (RFX) involving fifteen healthy subjects have been conducted on the contrast maps we previously investigated at the subject level. More precisely, one- sample Student-t test was performed on the subject-level contrast images (e.g., Lc-Rc, aC-aS,. images) using SPM5.

Figure 19 shows group-level student-t maps for the aC-aS contrast where data have been reconstructed using the mSENSE, UWRSENSE and 4D-UWR-SENSE for R = 2 and R = 4. Neurological convention has been used. Arrows indicate the global maximum activation peak.

These maps illustrate that irrespective of the reconstruction method larger and more significant activations are found on datasets acquired with R = 2 given the better SNR. Second, for R = 2, visual inspection confirms that only the 4D-UWR-SENSE algorithm allows to retrieve significant bilateral activations in the parietal cortices (see axial MIP slices) in addition to larger cluster extent and a gain in significance level for the stable clusters across the different reconstructors. Similar conclusions can be drawn when looking at the bottom of the figure, for R = 4. Complementary results are available in Table 5 for R = 2 and R = 4 and numerically confirms this visual comparison:

• Whatever the reconstruction method in use, the statistical performance is much more significant using R = 2, especially at the cluster level since the cluster extent decreases by one order of magnitude.

· Voxel and cluster-level results are enhanced using the 4D-

UWR-SENSE approach instead of the mSENSE or UWR-SENSE.

Table 5 - Significant statistical results at the group-level for the aC-aS contrast (corrected for multiple comparisons at p = 0.05).

Figure 20 reports similar group-level MIP results for R = 2 and

R = 4 concerning the Lc-Rc contrast. It is shown that whatever the acceleration factor R in use, our pipeline enables to detect a much more spatially extended activation area in the motor cortex. This visual inspection is quantitatively confirmed in Table 6 when comparing the detected clusters using the 4D-UWR-SENSE approach with those found by mSENSE, again irrespective of R. Finally, the 4D-UWR-SENSE algorithm outperforms the UWR-SENSE one, which corroborates the benefits of the proposed spatio- temporal regularization scheme. Table 6 - Significant statistical results at the group-level for the Lc-Rc contrast (corrected for multiple comparisons at p = 0.05).

References

[Bertsekas, 1995] Bertsekas, D. P. (1995). Nonlinear programming, Second Edition. Athena Scientific, Belmont, USA. In particular, pages 159-165.

[Block et al., 2007] Block, K. T., Uecker, M., and Frahm, J. (2007). Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint. Magnetic Resonance in Medicine, 56(7): 1086-1098.

[Chaari et al.2008]: L. Chaari, J.-C. Pesquet, A. Benazza- Benyahia, P. Ciuciu, Autocalibrated Parallel MRI Reconstruction in the Wavelet Domain, in: IEEE International Symposium on Biomedical Imaging, Paris, France, 2008, pp. 756-759. [Chaari et al. 2009]: L. Chaari, J.-C. Pesquet, Ph. Ciuciu, and A. Benazza-Benyahia An Iterative Method for Parallel MRI SENSE-based Reconstruction in the Wavelet Domain,. arXiv:0909.0368v1 [math.OC]

[Chaux et al., 2007] : Chaux, C, Combettes, P., Pesquet, J.- C, and Wajs, V. (2007). A variational formulation for frame-based inverse problems. Inverse Problems, 23(4):1495-1518.

[Combettes and Pesquet, 2008]: Combettes, P. L. and Pesquet, J.-C. (2008). A proximal decomposition method for solving convex variational inverse problems. Inverse Problems, 24(6):27.

[Daubechies et al., 2004] : I. Daubechies, M. Defrise, C. DeMol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications on Pure and Applied Mathematics 57 (11) (2004) 1413-1457.

[Dempster, 1997]: Dempster, A., Laird, A., and Rubin, D. (1977). Maximum likelihood from incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society, Series B, 39:1-38.

[Dobigeon et al, 2007]: Dobigeon, N. and Tourneret, J.-Y. (2007). Truncated multivariate Gaussian distribution on a simplex. Technical report, University of Toulouse.

[Donoho, 1995]: D. Donoho. De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41 (3):613-627, 1995.

[Geman, 1984]: S. Geman and D. Geman, "Stochastic relaxation, Gibbs distribution and the Bayesian restoration of image," IEEE Trans. Pattern Anal. Mach. Intell., vol. 6, pp. 721-741 , 1984. [Griswold et al., 2002] : M. A. Griswold, P. M. Jakob, R. M. Heidemann, M. Nittka, V. Jellus, J. Wang, B. Kiefer, A. Haase, Generalized autocalibrating partially parallel acquisitions GRAPPA, Magnetic Resonance in Medicine 47 (6) (2002) 1202-1210.

[Gupta et al., 1997]: Gupta, A. K. and Song, D. (1997). Journal of Statistical Planning and Inference. Volume 60, Pages 241-260. [Hoge et al., 2005] : W. S. Hoge, D. H. Brooks, B. Madore, W.

E. Kyriakos, A tour of accelerated parallel MR imaging from a linear systems perspective, Concepts in Magnetic Resonance 27A (1) (2005) 17-37.

[Pesquet-Popescu, Pesquet] Ondelettes et applications, Techniques de I'lngenieur, traite Telecoms, TE 5 215

[Pinel et al., 2007]: P. Pinel, B. Thirion, S. Meriaux, A. Jobert, J. Serres, D. Le Bihan, J.-B. Poline, and S. Dehaene, "Fast reproducible identification and large-scale databasing of individual functional cognitive networks" BMC Neurosci., vol. 8, no. 1 , pp. 91 , Oct. 2007.

[Pruessmann et al., 1999] : K. P. Pruessmann, M. Weiger, M. B. Scheidegger, P. Boesiger, SENSE : sensitivity encoding for fast MRI, Magnetic Resonance in Medicine 42 (5) (1999) 952-962.

[Raj et al., 2007] : Raj, A., Singh, G., Zabih, R., Kressler, B., Wang, Y., Schuff, N., and Weiner, M. (2007). Bayesian parallel imaging with edge-preserving priors. Magnetic Resonance in Medicine, 57(1):8-21 . [Robert, 1995]: "Simulation of truncated normal variables,"

Statistics and Computing, vol. 5, pp. 121-125, 1995. [Serra, 1982] Serra, J. (1982). Image Analysis and Mathematical Morphology. Academic Press, London.

[Sodickson et al., 1997] : Sodickson, D. K. and Manning, W. J. (1997). Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magnetic Resonance in Medicine, 38(4):591-603.

[Ying et al., 2004] : L. Ying, D. Xu, Z.-P. Liang, On Tikhonov Regularization for image reconstruction in parallel MRI, in: IEEE Engineering in Medicine and Biology Society, San Francisco, USA, 2004, pp. 1056-1059.