ELSEVIER

Contents lists available at ScienceDirect

Magnetic Resonance Imaging

journal homepage: www.mrijournal.com

CrossMark

Original contribution

Intra voxel analysis in magnetic resonance imaging

Michele Ambrosanioa, Fabio Baselice3-*, Giampaolo Ferraiolib, Flavia Lentic, Vito Pascazioa

aDipartimento di Ingegneria, Université di Napoli Parthenope, Italy b Dipartimento di Scienze e Tecnologie, Université di Napoli Parthenope, Italy cEUMETSAT, Darmstadt, Germany

ARTICLE INFO

ABSTRACT

Article history: Received 18 April 2016 Accepted 12 November 2016 Available online xxxx

Keywords:

Magnetic resonance imaging Sparse Bayesian learning Intra voxel analysis Statistical signal processing

A technique for analyzing the composition of each voxel, in the magnetic resonance imaging (MRI) framework, is presented. By combining different acquisitions, a novel methodology, called intra voxel analysis (IVA), for the detection of multiple tissues and the estimation of their spin-spin relaxation times is proposed. The methodology exploits the sparse Bayesian learning (SBL) approach in order to solve a highly underde-termined problem imposing the solution sparsity. IVA, developed for spin echo imaging sequence, can be easily extended to any acquisition scheme. For validating the approach, simulated and real data sets are considered. Monte Carlo simulations have been implemented for evaluating the performances of IVA compared to methods existing in literature. Two clinical datasets acquired with a 3T scanner have been considered for validating the approach. With respect to other approaches presented in literature, IVA has proved to be more effective in the voxel composition analysis, in particular in the case of few acquired images. Results are interesting and very promising: IVA is expected to have a remarkable impact on the research community and on the diagnostic field.

© 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND

license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

Magnetic resonance imaging (MRI) is a very useful technique for generating images of body slices, with particular reference to soft tissues. The nuclear magnetic resonance (NMR) effects are exploited in order to let hydrogen nuclei transmit a radio frequency signal. By collecting and properly processing such signals, classical MR images are obtained [1].

In order to spatially differentiate the contributions from different regions, magnetic field gradients are employed. In particular, phase and frequency encoding schemes are adopted. As a consequence, the imaged object slice is divided into voxels, i.e. volume elements that produce each pixel of the final image. The size of the voxel is an indicator of the system spatial resolution and, together with the dimension of the field of view (FOV), define the number of pixels of the obtained images. In order to achieve a finer resolution, i.e. to be able to identify very small and close objects or structures, more frequency and phase encoding steps have to be adopted. While frequency encoding can be improved by exploiting more powerful gradients, phase encoding steps cannot be increased limitless, as this leads to high acquisition times. A way to limit the problem drives

* Corresponding author.

into the implementation of parallel acquisition schemes, such as SENSE [2] or GRAPPA [3].

Although modern acquisition systems are able to reach a submillimeter spatial resolution, very small anatomical elements hardly become visible. Signals coming from the hydrogen nuclei of such structures are coherently mixed with signals coming from other anatomical elements of the rest of the voxel. The possibility of discriminating small structures starting from the coherent sum of the recorded signals (the acquired image) would be very helpful for diagnostic purposes. To guarantee this possibility, a new formulation of the problem is required: instead of further improving the spatial resolution of the system, a methodology could be developed for analyzing the signal related to each voxel in order to find its composition. The approach is similar to spectroscopy, in which a volume is analyzed looking for the presence of different chemical structures. Anyway, differently from spectroscopy which is characterized by a coarse resolution (in the order of cm3), the proposed methodology should be able to work with images acquired at full resolution (in the order of mm3).

Techniques for detecting and estimating signals coming from different contributions but received in overlapped manner have been developed for a wide range of research fields, such as radar domain [4] and microwave imaging [5]. The main difficulty consists in handling highly underdetermined mathematical problems [6-8]. With respect to MRI field, some approaches for detecting the

http://dx.doi.org/10.1016/j.mri.2016.11.009

0730-725X/ © 2016 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

presence of different tissues and estimating their T2 relaxation times have been proposed, such as non-negative least square (NNLS) or f1 norm minimization [9,10]. These methodologies have not been broadly adopted due to the high number of required acquisitions (tens to hundreds of images). This is mainly a mixed effect of the high underdetermination of considered problems and of the high correlation between the acquisition matrix rows. The number of acquisitions can be largely reduced if a strong a priori information is available. For example, in the estimation of Myelin water fraction(MWF) the knowledge of relaxation times of involved components allows a reduction of required data [11].

Within this manuscript we propose an approach based on sparse Bayesian learning(SBL) technique, which has demonstrated to be more capable of finding sparse solution than conventional approaches (convex f1 norm minimization), in particular in the case of ill-conditioned problems, as the one we are referring to [12,13]. More in detail, classical f1 norm minimization performs poorly when the restricted isometry is violated, while SBL is better suited in such kind of problems [14,15].

The proposed technique, called intra voxel analysis (IVA), aims at detecting and differentiating the different signals that are merged within each pixel of the obtained image. Specifically, for each voxel the detection of then number of the involved contributions is performed. Jointly with this step, the T2 relaxation times of each detected tissue are estimated. With the proposed methodology, reliable results can be achieved with less than 10 acquisitions.

In the following, after the presentation of the methodology, results on both simulated and real datasets are reported. Finally, conclusion are drawn.

2. Material and method

In the following paragraphs the acquisition model and the proposed IVA approach are presented. After introducing the acquired signal as the coherent summation of multiple tissues contributions, the procedure for the detection of the number of the involved contributions and the estimation of their spin-spin relaxation time is derived.

2.1. Acquisition model

In MRI, signals are acquired in the so called k-space and corrupted by an additive white circular complex Gaussian noise [16]. After the application of inverse 2D Fourier transform, the image in the space domain is formed [17,18]. Due to the linearity of Fourier transform, in space domain noise samples are still distributed as circular complex Gaussian variables. Generally, only the amplitude of the acquired data, corrupted by Rice distributed noise [16,19], is considered in clinical diagnosis.

Several imaging sequences have been proposed. The choice primary depends on the application we are interested in. In the framework of T2 relaxation time estimation, spin echo (SE) imaging sequence has been widely adopted [20,21]. By considering the i-th row and j-th column voxel within the s-th slice in the noise free case, the signal amplitude gj can be written as

different images, each one with a specific echo time TE(n), n = 1,..., N. For each voxel, a data vector is collected:

gi,j,s = kij,s exp (-TT-)

where T2,y,s is the spin-spin relaxation time of the considered voxel, kij,s is a coefficient which takes into account proton density and spin-lattice relaxation time and TE is the echo time of the imaging sequence. In the following, subscripts i, j, and s will be neglected in order to simplify the mathematical treatment. Let us acquire N

g(n) = k exp ( -

Spin-spin relaxation time estimation consists in evaluating T2 by exploiting signals g(n) acquired with different TE(n) [19,20,22-24].

Let us focus on a voxel containing P different tissues, each one characterized by its own spin-spin relaxation time T2(p) and intensity k(p). The acquisition model of Eq. (2) can be recast in order to include all the contributions, leading to

g(n)=£ «M - Tgj).

2.2. Intra voxel analysis

The purpose of IVA consists in detecting the number of tissues P inside each voxel and in estimating their spin-spin relaxation times T2(p). This allows us to identify the presence of elements smaller than the resolution cell. For this purpose, we have to define an effective way for describing the composition of any imaged voxel. Let us consider two tissues present in the considered voxel, with parameters k(1) = 1, T2(1) = 50 ms and k(2) = 1.5, T2(2) = 85 ms, respectively. First, let us define a vector T of possible T2 values, whose M elements are uniformly sampled within the interval [T2MIN, T2MAX]. For example, we could assume T = [1,2, ...,150] ms. Moreover, we define x as the vector of k related to each T2 value contained in T. Clearly, x has the same length of T. In the considered case, the vector x assumes value k(1) in position T2(1) and value k(2) in position T2(2), and all other values equal to zero. Therefore, x(50) = 1 and x(85) = 1.5. In other words, each non-zero element of x describes a different tissue, as its index is related to the relaxation time, and its value to k. In a more general way, x = [k[T2MIN]... k[T2MAX]]T, where k[T2] is the coefficient related to a specific spin-spin relaxation time.

Fig. 1. Example of x vector describing a voxel containing two tissues with k(1) : T2(1) = 50 ms and k(2) = 1.5, T2(2) = 85 ms. Note the sparsity of the x vector.

In Fig. 1, the vector x related to the considered case is reported. It is worthy to underline that, by analyzing such a vector, we are at a glance able to detect presence of multiple tissues within the considered voxel and their spin-spin relaxation times. Since all but P elements of the vector are equal zero, x can be considered sparse.

In order to write the acquisition model of Eq. (3) in a vectorial form, the echo times vector TE = [TE(1)TE(2)... TE(N)]T and the data vector g = [g(1)g(2).. .g(N)]T are defined, obtaining

g = Ax

where the (n, m) element of transformation matrix A is defined as

A(n, m) = exp I -

(-T^ï ; n V T(m)/ ;

= 1, ••• ;N, m = 1, ••• , M

It can be noticed that matrix A has a number of rows equal to the number of acquired images, while the columns number depends on the length of vector T, i.e. it is related to T2min, T2MAX and to the considered discretization step. Due to the high correlation between rows, which are exponential curves decreasing with different rates, matrix A is ill-conditioned. The high underdetermination of the problem implies that the solution is not unique. In order to force solution uniqueness, we have to impose some a priori information on the signal to be reconstructed. Thus, to invert the model of Eq. (4), the maximum a posteriori (MAP) approach is implemented, which consists in finding

xmap = arg maxp(xjg)

where p(x|g) is the a posteriori distribution of the unknown vector x elements given the acquired data g. Such distribution can be fac-torized as the product of the likelihood function p(g|x) and the so called a prior distribution p(x), modifying the maximization problem of Eq. (6) to

xmap = arg min - logp(gjx) - logp(x)

If independent zero-mean Gaussian model is assumed for noise samples, multi-variate likelihood of the observations g is given by

p(gjx) = (2n)- N s-N

exp{-II g - Ax||2}.

where s is the noise standard deviation.

Several approaches belong to MAP class of problems, the main difference is given by the choice of the a priori model. The most natural a priori distribution is the Gaussian function:

p(x) ~ exp

If the hypothesis of non-negative x values is done, the nonnegative least square (NNLS) approach is obtained:

arg min ^ || g - Ax||2+ || x||

x>0 s2N

Such method has been widely adopted [9-11], but it is well-known to produce non-sparse solution. We are interested in achieving high accuracy in T2 estimation, i.e. large M, and meanwhile in working with a limited number of acquired data, i.e. small N, thus NNLS is not well suited, since it requires a high amount of data.

The x vector is expected to have few non-zero coefficients (few tissues within each voxel), i.e. it is expected to be sparse. Thus some a priori model p(x) for encouraging sparsity should be adopted. One of the most adopted a priori distribution for promoting sparsity is the generalized Gaussian one [25]:

p(x) ~ exp

where a assumes values within the [0,1] interval.

When a = 1, the classical l1 -norm minimization is obtained [26]. It performs poorly when the restricted isometry property is violated: the global minimum of the cost function does not coincide with the sparsest solution of Eq. (4) [14]. When a e [0,1), the cost function has several local minima and for a tends to 0, the minimization of Eq. (7) tends to an NP problem [27].

We propose to adopt for inverting Eq. (4) an approach called sparse Bayesian learning (SBL). This method allows automatic a priori model selection maximizing the marginal likelihood.

In details, we choose as the prior distribution a parametrized Gaussian:

M _i x2 p(x, c) = n (2pyi) 1 exp(- 2Ç- ), i=1 c

where y = [y1,..., yM] is a vector of hyperparameters. It is the key of the model sparsity.

The hyperparameters can be estimated from the data. To find y, it is enough to maximize the marginalized pdf p(g, y, s2), that is equivalent to minimizing the effective SBL cost function:

L = log jSg j + gT S-1g,

where Sg = ct2I + ATA and r = diag(y). Once hyperparameters values have been estimated, the solution x is given by

xSBL = (ATA + s2r-1) UTg

More details about the algorithm and the role of the hyperparameters can be found in [15,28]. In particular, in [15] it is proven that the SBL algorithm produce a solution with higher sparsity than the l1 one and that the local minima of the SBL cost function are achieved at sparse solutions.

3. Simulation and experiments

In order to validate the proposed technique, tests on both simulated and real datasets have been carried out. The performances of IVA have been compared to classical NNLS approach. In the following, we considered the signal to noise ratio (SNR) expressed in dB and computed according to

SNR = 10log10

N=1jg.j2 ^

All the processing has been done with Matlab® software on a Linux Debian based Dell® Optiplex 990 workstation.

A comment on the discretization step of vector T, i.e. about the length of vector x, is worthy to be done. Although we found it not to be a critical aspect of the methodology, an indication about its correct value range can be found by looking at the Cramer-Rao lower bounds

(CRLB) of T2 estimator, i.e. by looking at the expected performances of the optimal spin-spin relaxation time estimator. Evaluations of CRLB for different parameters, such as SNR and number of acquisitions N, can be found in [19] and [23]. From the statistical estimation theory, a good discretization step should be smaller than CRLB in order not to negatively affect estimation accuracy. From the results reported in [23], in the case of N < 10, CRLB are approximately 5 ms for the T2 estimator, thus a smaller discretization step should be used. In the following, a discretization step equal to 2 ms is considered.

3.1. Simulated data

The Monte Carlo method has been adopted for evaluating the performances of proposed IVA. In order to provide a reference, the NNLS technique, which is a well established method for carrying out the multiexponential decomposition, has been considered. Estimations varying both the number of acquired signals and the SNR have been conducted, each one composed of 2000 simulations. A single voxel containing two tissues with k = [0.5,1] and T2 = [60,160] ms has been considered. The unknown vector x has been constructed considering T2 relaxation times between 8 ms and 250 ms, with a sampling step equal to 2 ms. Both IVA and NNLS methods have been tested for a number of spin echo acquisitions N between 4 and 20, with echo times TE equally spaced in the [20, 300] ms interval, and considering SNR values between 20 and 60 dB. For each Monte Carlo

run, the maximum value of the estimated vector x is detected. A threshold equal to 1 /20 of the maximum is applied to x: the values higher than the threshold are identified as contributions and their positions in terms of T2 value are stored. This criterion is adopted for both IVA and NNLS. The threshold value, i.e. 1/20 of the maximum, has empirically proven to be a good choice. Results reported in Figs. 2 and 3 are described and discussed in the following sections.

3.2. Real data

IVA methodology has also been tested on two real datasets. The available images were acquired by using a Philips Achieva 3.0T MRI scanner and refer to a 49 years old female patient, whose informed consent was obtained.

The first dataset is related to a cerebral edema in the occipital region and is composed of 5 spin echo images acquired with TE = [82.45,130,180,250,350]T ms and Repetition Time fixed to 3400 ms. The image with TE = 82.45 ms is reported in Fig. 4a, while in Fig. 4b the image with contrast enhancement is shown. The latter has been acquired few days later. We restricted the analysis to a Region Of Interest (ROI) of the image, highlighted by the red square in Fig. 4.

For any image pixel, x vector has been estimated. Again, we searched for components with T2 in the [10, 250] ms range, with a discretization step equal to 2 ms. Results are reported in Figs. 5-8.

A second case study is reported. The dataset refers to an occipital necrotic lesion. In this case, 4 SE acquisitions were available, with

Fig. 2. Histograms of the Monte Carlo simulation results in the case of different acquisition number N and SNR combinations. Blue curves are related to the IVA solutions, while red curves are for the NNLS approach. The case study takes into account two tissues with spin-spin relaxation times equal to 60 ms and 160 ms. The two populations, centered at 60 and 160, should be clearly separated. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Mean value (line marked with asterisks) and 95% confidence interval (dotted lines) for IVA (red lines) and NNLS (blue lines). Reference values are reported in black. The images report the results in the case of tissue with T2 = 60 ms and SNR = 20 dB (a), tissue with T2=60 ms and SNR = 40 dB (b), tissue with T2 = 160 ms and SNR = 20 dB (c), and tissue with T2 = 160 ms and SNR = 40 dB (d). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

TE = [82.45,130,180,250]T ms, respectively. The image acquired with TE = 82.45 ms is reported in Fig. 9a, while in Fig. 9b the image with contrast enhancement is shown. Again, we restricted IVA

analysis to a region of the image, highlighted by red square in Fig. 9, where the lesion can be easily seen. Results are reported in Figs. 10 and 11.

Fig. 4. First real case study: SE image (a) and IR with contrast enhancement (b). Red square highlights the region of interest. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. Selection of two points (red and green) for IVA application (a). IVA results for selected green (b) and red (c) points. For both of them, a single tissue has been detected with spin-spin relaxation times equal to 56 and 144ms, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4. Results

4.1. Simulated data

In Fig. 2, the histograms of the estimated spin-spin relaxation times are reported for both IVA (blue curve) and NNLS (red curve). In order to better appreciate the performances of IVA and NNLS, the estimated values with the 95% confidence interval for the T2 estimation of the two tissues (T2 = 60 ms and T2 = 160 ms) in the case of SNR = 20 dB and SNR = 40 dB have been computed and are reported in Fig. 3.

4.2. Real data

First, in order to validate the discrimination capability of the methodology, the detection has been focused on two points, marked with green and red dots in Fig. 5a. IVA results are reported in panels b and c inFig. 5. For both voxels a single tissue has been identified, with a spin-spin relaxation time equal to 56 ms for the green dot and 144 ms for the red one. Starting from the signals related to the two tissues, the signal coming from the macro-voxel containing the two tissues has been generated as the coherent sum of the two voxels. IVA has been applied to this macro-voxel, providing results of Fig. 6. It is evident that IVA is reliable in detecting the presence of the two involved tissues and their relaxation times.

We extended the analysis to the whole ROI. First, the presence of multiple contributions within each voxel is detected by using IVA. In Fig. 7, pixels where IVA detected the presence of more than one tissue are highlighted in red.

IVA and NNLS estimations for the considered ROI are reported in the first and in the second row of Fig. 8. In particular, the left column shows the estimation of T2 value for the first component, while the right column refers to the second one. In order to validate the estimated spin-spin relaxation time values, a Least Square fitting approach by adopting the bi-exponential model has been implemented, and its results are reported in panels e and f in Fig. 8 [29].

Clearly, this is not exactly a fair comparison as IVA and NNLS should be able to both detect the number of components and estimate their T2, while LS fitting approach based on bi-exponential model a priori assumes the presence of two components in each pixel and estimates their spin-spin relaxation times, i.e. has no detection step. In other words, the LS fitting approach requires the knowledge of the involved tissues number, while IVA and NNLS not. The lack of the detection step implies the impossibility of masking no signal areas, thus the gray color mask (no detected tissue) of IVA and NNLS results cannot be applied. That said, the estimated T2 can be used as reference in order to validate IVA and NNLS results.

Moving to the second real dataset, regions where a second contribution has been identified are shown in red, overlapped to the original image in Fig. 10. Estimated spin-spin relaxation times for

1 -1-T-

0.8 0.6 0.4 0.20- --

0 50 100 150 200 250

Fig. 6. Combination of signals related to two voxels, marked with a blue area (a), and IVA result (b). As expected, the two involved tissues are both detected. Relative errors in T2 estimations are equal to 7% and 11% for left and right peak, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version ofthis article.)

both detected components are reported in Fig. 11. Again, for comparison purposes, results achieved via NNLS approach are reported in panels c and d in Fig. 11, while the estimation via LS fitting method in panels e and f in Fig. 11.

5. Discussion

5.1. Simulated data

From Fig. 2, the effectiveness of IVA in estimating two separate distributions, related to the T2 = 60 ms and to the T2 = 160 ms contributions, is evident. On the other hand, NNLS tends to merge them, in particular in the case of low SNR. More in detail, some additional comments can be done:

• In the case of high SNR (> 40 dB) and few acquisitions (N = 4), NLLS approach fails in retrieving the correct spin-spin relaxation time values, while IVA is effective in finding the correct T2 values and in separating the two distributions.

• Both techniques show good results in the case of high SNR (> 40 dB) and several acquisitions (N > 10).

Fig. 7. IVA multiple contribution detection. Pixels with the presence of more than one tissue have been highlighted in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

• In the case of mid to low SNR (<40 dB), NNLS often finds a tissue with T2 within the [100,150] ms range, which is not correct.

• Only in the case of low SNR/small N and high SNR/big N combinations, NNLS and IVA results are comparable. In all other cases, IVA histograms are closer to the ideal ones than NNLS.

Moving to the 95% confidence intervals for the T2 estimation (Fig. 3), with respect to the T2 = 60 ms tissue (first line), IVA is globally characterized by a smaller bias and shows a narrower 95% confidence interval. Moreover, the performance improvement adding further acquisitions is clearly visible in the case of SNR = 40 dB. In the case of the T2 = 160 ms tissue results, some aspects can be highlighted. Looking at the mean estimated value, IVA is better in the case of SNR = 20 dB, while NNLS shows a lower bias in the case of SNR = 40 dB. Looking at the 95% confidence intervals, they are very similar in the case of SNR = 20 dB, while the IVA one is sensibly narrower in the SNR = 40 dB case. By considering jointly both findings, the T2 estimator provided by IVA is more reliable than NNLS in all cases.

5.2. Real data

By looking at the images reported in the first row of Fig. 8, IVA results seem reliable since, as expected, most of multiple contributions are within the edema region. Some other red spots are present, mainly where tissues edges are located.

Besides the detection of the contribution number, IVA is able to estimate the spin-spin relaxation time of each of them (Fig. 8a and b). By focusing on the edema area in the center of the picture, we can easily note that the first component is characterized by a spinspin relaxation time of about 30 ms, which is slightly lower than the surroundings and with respect to white matter T2 values reported in literature [30]. The same behavior can be noticed in Fig. 8b, where a T2 value of about 180 ms is measured while a slightly higher value is reported in [30] for blood. A possible explanation to this could be related to the blood accumulation and stagnation. Of course, a deeper biological investigation is required, but it is beyond the scope of the paper. The NNLS approach has been also applied to this dataset. Results are reported in panels c and d in Fig. 8 for the two components, with the same color maps adopted for IVA. In this case, the detected voxels containing a second component (Fig. 8d) are

Fig. 8. Estimation results, first dataset. The images of the left column show the first detected tissue within each voxel for IVA (a), NNLS (c) and bi-exponential model LS fitting (e), and the color provides the T2 estimated values. The images of the right column show points where a second tissue has been detected for the three approaches, again T2 values are color coded. Points marked in gray mean no tissue detected. In the case of bi-exponential model LS fitting (e and f), no point has been masked as no detection step is considered. Note that the images of the two tissue components (left and right) have different color maps in order to improve the readability. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

isolated, and the solution is not fully convincing. The first component (Fig. 8c) is characterized by a much higher T2 value with respect to IVA solution, while the second is almost not existing.

By focusing on panels a and e in Fig. 8, i.e. the estimation of component with lower T2, LS fitting provides generally lower T2 values with respect to IVA. Clearly, as the true T2 values are not known, it is hard to state whether the bias should be attributed to IVA or LS solution. However, even if the bias is ascribed to IVA, it would not be a dramatic issue, since it is limited in value, as confirmed by Monte Carlo simulations of Fig. 3, and stable across the image. By inspecting the edema area, where both techniques detect two contributions, the LS finds a first component characterized by a very low spin-spin relaxation time (about 20 ms) (Fig. 8e), which does not match with the existing literature [30]. On the other hand, IVA provides for that

component a more consistent value of about 30 ms. The NNLS result in the case of lower T2 tissue reported in Fig. 8c shows completely different values in most of the regions, i.e. the estimation cannot be considered effective. Moving to the comparison of results related to the higher T2 component, the edema region is the most interesting area. Results of IVA and LS fitting are very similar, while NNLS solution is totally unsatisfactory, as very few components are detected and a lot of false positive are evident.

Considering the second real dataset (Fig. 9), results reported in Fig. 11 appear less accurate, probably due to the reduced number of available images (4 instead of 5). With respect to the lesion area, an interesting result can be highlighted: the presence of the second component is detected by IVA in the top and in the lower parts (Fig. 11b), while in the middle region only one contribution is

(a) (b)

Fig. 9. Second real case study: SE image (a) and IR with contrast enhancement (b). Red square highlights the region of interest. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

revealed. This could be related to the difference between active and necrotic regions within the lesion. This result is in accordance with other studies on the same dataset [17,31]. Concerning the comparison with the other considered techniques (second and third rows of Fig. 11), results are in accordance with the previous case study: NNLS completely fails in detecting the second tissue, while IVA and LS, neglecting a small bias, are consistent.

6. Conclusions

Within this manuscript, a novel methodology based on sparse Bayesian learning is presented in the magnetic resonance imaging field. The approach, named intra voxel analysis (IVA), is capable of evaluating the presence of multiple tissues, or in general of different contributions, within each voxel. IVA is also able to estimate the spin-spin relaxation time of each involved tissue, allowing their discrimination after the detection step. With respect to existing approaches, such as NNLS, IVA has shown better performances in discriminating tissues, due to its high effectiveness in handling highly

underdetermined inversion problems. More in details, IVA showed good accuracy both in detecting the correct number of involved tissues and in their T2 estimation, without any a priori knowledge, with interesting performances also in the case of a limited number of acquired images (less than 10). At the present, IVA has been developed for spin echo imaging sequence, but it can easily extended to other sequences. Monte Carlo simulations quantitatively validated the methodology and compared it to other approaches. First results on clinical data have been presented. Further studies will be focused on the extension of IVA to the estimation of other MR intrinsic parameters, such as spin-lattice relaxation time, and on the validation of the methodology on a wider range of clinical datasets.

Acknowledgments

The authors would like to thank Rocchina Caivano and Aldo Cammarota, Department of Radiology, IRCSS CROB, Rionero in Vulture, Italy, for supplying real datasets and for the interesting discussions about obtained results.

(a) (b)

Fig. 10. IVA multiple contribution detection. Pixels with the presence of more than one tissue have been highlighted in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11. Estimation results, second dataset. The images of the left column show the first detected tissue within each voxel for IVA (a), NNLS (c) and bi-exponential model LS fitting (e), and the color provides the T2 estimated values. The images of the right column show points where a second tissue has been detected for the three approaches, again T2 values are color coded. Points marked in gray mean no tissue detected. In the case of bi-exponential model LS fitting (e and f), no point has been masked as no detection step is considered. Note that the images of the two tissue components (left and right) have different color maps in order to improve the readability. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

References

[1] Slichter CP. Principles of magnetic resonance. 3rd ed., New York: SpringerVerlag Berlin; 1990.

[2] Pruessmann KP, Weiger M, Schdeidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med 1999;42(5):952-62.

[3] Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med 2002;47(6):1202-10.

[4] Baselice F, Budillon A, Ferraioli G, Pascazio V. Layover solution in SAR imaging: a statistical approach. IEEE Geosci Remote Sens Lett 2009;6(3):577-81.

[5] Autieri R, Ferraiuolo G, Pascazio V. Bayesian regularization in nonlinear imaging: reconstructions from experimental data in nonlinearized microwave tomography. IEEE Trans Geosci Remote Sens 2011;49(2):801-13.

[6] Romberg J. Imaging via compressive sampling. iEeE Signal Process Mag 2008;25(2):14-20.

[7] Budillon A, Evangelista A, Schirinzi G. Three-dimensional SAR focusing from multipass signals using compressive sampling. IEEE Trans Geosci Remote Sens 2011;49(1):488-99.

[8] Autieri R, Pascazio V, D'Urso M. Compressive sensing image reconstruction in aspect-limited data. Adv Radar Remote Sens (TyWRRS), 2012 Tyrrhenian Workshop on 2012;304-8.

[9] Whittall KP, MacKay AL. Quantitative interpretation of {NMR} relaxation data. J Magn Reson (1969) 1989;84(1):134-52.

[10] Istratov AndreiA, Vyvenko OlegF. Exponential analysis in physical phenomena. RevSci Instrum 1999;70(2):1233-57.

[11] Bjork M, Zachariah D, Kullberg J, Stoica P. A multicomponent T2 relaxom-etry algorithm for myelin water imaging of the brain. Magn Reson Med 2016;75(1):390-402.

[12] Li Y, Yu ZL, Bi N, Xu Y, Gu Z, Amari SI. Sparse representation for brain signal processing: a tutorial on methods and applications. IEEE Signal Process Mag 2014;31(3):96-106.

[13] Zhang Z, Rao B. Extension of SBL algorithms for the recovery of block sparse signals with intra-block correlation. IEEE Trans Signal Process 2013;61(8):2009-15.

[14] Wu W, Nagarajan S, Chen Z. Bayesian machine learning: EEG/MEG signal processing measurements. IEEE Signal Process Mag 2016;33(1):14-36.

[15] Wipf DP, Rao BD. Sparse bayesian learning for basis selection. IEEE Trans Signal Process 2004;52(8):2153-64.

[16] den Dekker A, Sijbers J. Data distributions in magnetic resonance images: a review. Phys Med 2014; (0).

[17] Baselice F, Ferraioli G, Pascazio V. A Bayesian approach for relaxation times estimation in MRI. Magn Reson Imaging 2016;34(3):312-25.

[18] Eggers H, Knopp T, Potts D. Field inhomogeneity correction based on grid-ding reconstruction for magnetic resonance imaging. IEEE Trans Med Imaging 2007;26(3):374-84.

[19] Baselice F, Ferraioli G, Pascazio V. Relaxation time estimation from complex magnetic resonance images. Sensors 2010;10(4):3611-25.

[20] Bonny JM, Zanca M, Boire JY, Veyre A. T2 maximum likelihood estimation from multiple spin-echo amplitude images. Magn Reson Med 1996;36287-93.

[21] Sijbers J, den Dekker A, Raman E, Van Dyck D. Parameter estimation from amplitude MR images. IntJ Imag Syst Technol 1999;10109-14.

[22] Liu J, Nieminen A, Koenig JL. Calculation of T1, T2 and proton spin density in nuclear magnetic resonance imaging. J Magn Reson 1989;8595-110.

[23] Baselice F, Ferraioli G, Grassia A, Pascazio V. Optimal configuration for relaxation times estimation in complex spin echo imaging. Sensors 2014;14(2):2182-98.

[24] Baselice F, Ferraioli G, Pascazio V. A novel statistical approach for brain MR images segmentation based on relaxation times. Biomed Res Int 2015;201513.

[25] NadarajahS. A generalized normal distribution. J Appl Stat 2005;32(7):685-94.

[26] Fornasier M, Rauhut H. Compressive Ssensing. In: O, Scherzer, ed. Handbook of mathematical methods in imaging. Springer.2011. p. 187-228.

[27] Natarajan BK. Sparse approximate solutions to linear systems. SIAM J Comput 1995;24(2):227-34.

[28] Tipping ME. Sparse Bayesian learning and the relevance vector machine. J Mach Learn Res 2001;1211-44.

[29] Huang C, Galons JP, Graff CG, Clarkson EW, Bilgin A, Kalb B. Correcting partial volume effects in biexponential T2 estimation of small lesions. Magn Reson Med 2014; n/a-n/a.

[30] Stanisz GJ, Odrobina EE, Pun J, Escaravage M, Graham SJ, Bronskill MJ. T1, T2 relaxation and magnetization transfer in tissue at 3T. Magn Reson Med 2005;54(3):507-12.

[31] Baselice F, Caivano R, Cammarota A, Ferraioli G, Pascazio V. T1 and T2 estimation in complex domain: first results on clinical data. Concepts Magn Reson Part A 2014;43(5):166-76.