Keywords

1 Introduction

Diffusion Magnetic Resonance Imaging (dMRI) is the only non-invasive tool to obtain information about the neural architecture within the white matter of the brain in vivo. Several medical applications can be considered such as the 3D fibers reconstruction, the characterization of neuro-degenerative diseases and the surgical planning. The diffusion MRI consists in the study of the water molecules movement. In the case of brain tissues, this diffusion is anisotropic. Therefore, we need to perform models that characterize the diffusion and reproduce the preferred directions.

High Angular Resolution Diffusion Imaging (HARDI) [2] has become one of the most widely used methods to characterize the underlying fiber configurations in the presence of complex fiber crossings from diffusion-weighted MRI data. This technique is based on sampling the signal on a single or many spheres of the q-space. Among the HARDI techniques, there are the QBI (Q-Ball Imaging) method [7] and the spherical deconvolution one [3]. The first approach consists in the definition of a spherical function, the Orientation Distribution Function (ODF), that characterizes the probability of diffusion along a given angular direction. The second technique estimates the Fiber Orientation Distribution function (FOD) in each imaging voxel. It models the diffusion MRI signal as the convolution of the diffusion signal with a single fiber response.

However, although able to identify complex fiber configurations, HARDI techniques require a very large number of DW images to reach a satisfactory angular resolution and to give an adequate characterization of all the features of the diffusion signal. A long measurement time is, therefore, needed. It is problematic for clinical studies involving children and people afflicted with certain diseases.

For this reason, a large variety of methods have been recently proposed in the literature to shorten the acquisition time of the HARDI techniques. Donoho [10] was among the first to address this challenging field. He proposed a novel method qualified by Compressive Sensing (CS) to accurately reconstruct a signal from under sampled measurements acquired below the Shannon-Nyquist rate by using different representations. This method was applied in a large range of domains including image and video compression as well as geophysics and medical imaging. The CS technique was already applied in High Angular Resolution Diffusion Imaging (HARDI) to reduce the acquisition time.

Michailovich et al. [11] used this technique to represent the HARDI signals in the basis of spherical ridgelets by using a relatively small number of representation coefficients. Tristn-Vega et al. [12] proposed to represent the probabilistic Orientation Distribution Function (ODF) in the frame of Spherical Wavelets (SW). Duarte-Carvajalino et al. [4] directly estimated the constant solid angle orientation distribution function (CSA-ODF) from under-sampled multi-shell high angular resolution diffusion imaging (HARDI) datasets. Daducci et al. [5] used the CS method to estimate the FOD function with the spherical deconvolution method from limited samples.

In this work, we propose a novel method to reduce the acquisition time of the HARDI technique. It is based on a geometrical approach to construct the diffusion signal from a reduced number of gradient directions. These missing gradient directions will be obtained according to their neighborhood from a reduced set of data on the sphere of the q-space. The experimentations will be performed on simulated data and on digital phantoms.

Thus, this paper will be structured as follows: We present in the second section a brief recall of the principle of the HARDI acquisition technique. In the third section, we detail all the steps of the proposed approach. We test its accuracy to perform a good estimation of the FOD function using the Constrained Spherical Deconvolution (CSD) method on both simulated and digital phantoms in the fourth section.

2 Brief Recall of the HARDI Acquisition

Diffusion MRI provides a serie of images which is obtained by a direction variation of the encoding gradient. The application of a single pulsed gradient produces one Diffusion Weighted (DW) image that corresponds to one position in the q-space [9]. In HARDI imaging, the N discrete gradient directions are uniformly distributed on a single or many spheres of the q-space as illustrated in Fig. 1 (the case of one sphere). Each red dot on the surface of the unit sphere corresponds to a specific direction of the diffusion gradient. Along each of these gradient directions \( g_{i}(x_{i}, y_{i}, z_{i} ), 1< i < N\), a DW signal S(\( g_{i} \)) is measured.

Hence, the diffusion signal \(S=\{S(g_{1}), S(g_{2}),...S(g_{N})\}\) is obtained by the collection of the magnitudes in all diffusion-weighted images of the brain.

Fig. 1.
figure 1

Explanation of the diffusion-weighted images acquisition procedure in the HARDI technique. (Color figure online)

Full size image

3 Proposed Approach

This work investigates the possibilities of the reconstruction of the diffusion signal with a reduced number of gradient directions. Therefore, the acquisition procedure time is reduced and thus make the HARDI technique clinically feasible. We propose here a novel approach to characterize the complete water diffusion process in the white matter, with a small number of measurements.

Let N denote the full set of uniformly distributed directions on the sphere with their 3D coordinates. The value of N is chosen to be able to ensure a good construction of the DW signal. We consider \((n_{1} \) \(<N)\) a subset of the full set of directions with their DW signal values. It is qualified here by the partial original signal. We intend to estimate the (\( n_{2} = N- n_{1} \)) DW signal values. They correspond to the rest of points from the full set of directions. Figure 2-A illustrates the \( n_{1} \), \( n_{2} \), and the N points on the sphere. The blue dots denote the partial original diffusion data and red ones correspond to those the DW signal values are to be estimated.

Fig. 2.
figure 2

(A) Representation of the partial original signal and the data to be recovered on the sphere. (B) The Delaunay triangulation of the original signal (the \( n_{1} \) set of points). (C) An illustration of a point P (of \( n_{2} \)) and its nearest triangle.

Full size image

We propose, here, a novel geometrical approach in order to reach the real DW signal values of the \( n_{2} \) points. Since the DW signal is a continuous one that was sampled on the sphere of q-space [8], each DW value of a point on the sphere depends on the DW values of its neighboring points. In order to estimate the DW value of each point of the \( n_{2} \) set, its neighboring points from the \( n_{1} \) points (their DW values are known) should be determinated. We construct therefore the Delaunay triangulation of the \( n_{1} \) points. The points of the \( n_{1} \) set composing the nearest triangle to a point of \(n_{2}\) are considered as its neighboring ones. This nearest triangle corresponds to the one that has the minimum Euclidean distance between the considered point of \( n_{2} \) and all the triangles obtained after the Delaunay triangulation. Figure 2-B illustrates the Delaunay triangulation of the set of points \( n_{1} \). Let P be a point of the set \(n_{2}\) and \( T_{P_{1}P_{2}P_{3}}\) its nearest triangle composed by the points \(P_{1}\), \(P_{2}\) and \(P_{3}\) as showed in Fig. 2-C. The DW values of P (that we denote here by DW(P)) depends on the DW values respectively of \(P_{1}\), \(P_{2}\) and \(P_{3}\). We propose that DW(P) will be estimated as a linear combination of \(DW(P_{1})\), \(DW(P_{2})\) and \(DW(P_{3})\) as follows:

$$\begin{aligned} DW(P)=\alpha _{1}DW(P_{1})+\alpha _{2}DW(P_{2})+\alpha _{3}DW(P_{3}) \end{aligned}$$
(1)

with \(\alpha _{1} + \alpha _{2} + \alpha _{3} = 1\).

\(\alpha _{1}\), \(\alpha _{2}\) and \(\alpha _{3}\) are computed in the following sense:

We denote by \( T_{PP_{1}P_{2}}\), \( T_{PP_{1}P_{3}}\) and \( T_{PP_{2}P_{3}}\) the three triangles composed respectively by the points (P, \(P_{1}\) and \(P_{2}\)), (P, \(P_{1}\) and \(P_{3}\)) and (P, \(P_{2}\) and \(P_{3}\)) as shown in Fig. 2-C.

We qualify by \(A_{T_{P_{i}P_{j}P_{k}}}\) the area of a triangle composed by the point \(P_{i}\), \(P_{j}\) and \(P_{k}\). The values \(\{ \alpha _{i}, 1\le i\le 3\} \) characterize respectively the participation of the point \(\{ P_{i}, 1\le i\le 3\}\) in the estimation of DW(P). The more P is close to \(P_{i}\) the more \( \alpha _{i}\) should be large. The two triangles containing the point \(p_{i}\) will be used for the estimation of the value of \(\alpha _{i}\). In the case of the point p too close to \(p_{i}\), the areas of these two triangles are small. We formulate the value of \(\alpha _{i}\) as follows to ensure a better participation of the point \(p_{i}\) in the estimation of \(\alpha _{i}\).

Then:

$$\begin{aligned} \alpha _{i}= 1-A_{p_{i}} \end{aligned}$$
(2)

Where

$$\begin{aligned} A_{p_{i}}= \frac{A_{T_{PP_{i}P_{j}}}+A_{T_{PP_{i}P_{k}}}}{A_{T_{PP_{i}P_{j}}}+A_{T_{PP_{i}P_{k}}}+A_{T_{PP_{j}P_{k}}}} \end{aligned}$$
(3)

\(A_{p_{i}}\) corresponds to the normalized portion of the area of the two triangles containing the point \(P_{i}\).

\(\alpha _{i}\) is equal to (\(1-A_{p_{i}}\)) to ensure that the value of \(\alpha _{i}\) is large when P is too close to \(P_{i}\). Therefore,

$$\begin{aligned} \alpha _{1}= 1- \frac{A_{T_{PP_{1}P_{2}}}+A_{T_{PP_{1}P_{3}}}}{A_{T_{PP_{1}P_{2}}}+A_{T_{PP_{1}P_{3}}}+A_{T_{PP_{2}P_{3}}}} \end{aligned}$$
(4)
$$\begin{aligned} \alpha _{2}= 1- \frac{A_{T_{PP_{2}P_{1}}}+A_{T_{PP_{3}P_{2}}}}{A_{T_{PP_{1}P_{2}}}+A_{T_{PP_{1}P_{3}}}+A_{T_{PP_{2}P_{3}}}} \end{aligned}$$
(5)
$$\begin{aligned} \alpha _{3}= 1- \frac{A_{T_{PP_{3}P_{1}}}+A_{T_{PP_{2}P_{3}}}}{A_{T_{PP_{1}P_{2}}}+A_{T_{PP_{1}P_{3}}}+A_{T_{PP_{2}P_{3}}}} \end{aligned}$$
(6)

4 Experimental Results

In order to evaluate the accuracy of the proposed method to estimate the HARDI signal, we will use the Constrained Spherical Deconvolution (CSD) [6] method. This technique aims to estimate the Fiber Orientation Distribution (FOD) within each voxel directly from the diffusion-weighted (DW) data, using the concept of Spherical Deconvolution (SD) [3]. Tournier et al. [1] prove that a minimum of 45 DW directions is sufficient to fully characterize the DW signal with a maximum harmonic degree of 8 (Lmax = 8) and for a b-value equal to 3000 s/\(\mathrm{{mm}}^{2}\) (The b-value is a parameter that reflects the intensity and the duration of the gradient pulses used to generate the diffusion-weighted images).

We used a full set of directions \(N=65\) uniformly distributed on the sphere to avoid issues with imperfections in the uniformity of the DW gradient directions [1] and we varied the subsets \( n_{1}\) from 26 to 40 in order to evaluate the impact of the number of samples. After recovering the missing data with the proposed method, we obtain a new diffusion signal composed by the \(n_{1}\) directions and the recovered data. Then, we calculated the fiber orientation distribution (FOD) using CSD (Lmax = 8). We reconstruct the FOD function from HARDI data with the novel estimated signal. The obtained FOD with the full original set directions (the N original values of DW signal) is used as a ground truth. In order to test the performance of the proposed method, we generate a two-tensor model simulating two fibers crossing in different angles: 50, 70, and 90. We also evaluate our method on the data acquired from many digital phantoms of fibers crossing.

Results on Simulated Data. We simulate fiber crossing by generating Noise free DW data from the sum of two exponentials characterized by two identical fibers described by a diffusion tensor equal to \([0.3~0.3~1.7]\times 10^{-3}\,\mathrm{{mm}}^{2}/\mathrm{{s}}\) using 65 directions and a b-value equal to 3000 s/\(\mathrm{{mm}}^{2}\)

$$\begin{aligned} S(g_{i})= \exp (-b g_{i}^{t}D_{1}g_{i})+\exp (-b g_{i}^{t}D_{2}g_{i}) \end{aligned}$$
(7)

Where \(D_{1}\) is a diagonal matrix with diagonal entries \( [0.3~0.3~1.7 ]\times 10^{-3}\,\mathrm{{mm}}^{2}/\mathrm{{s}}\) and \(D_{2}\) corresponds to \(D_{1}\) rotated respectively by the angles 90, 70 and 50. For the first part of the experimentation, we calculate the FOD function. For the second one, we extract the local maxima of the FOD to calculate the crossing angle. By increasing the number of the \(n_{1}\) points set from 26 to 40, we show that as few as 34 diffusion directions are sufficient to obtain an accurate FOD estimation. Therefore, we can reduce the number of gradient acquisition by approximately a half and thus the scan time is reduced by 50 %. Table 1 illustrates the obtained crossing angles according to the values of \(n_{1}\). Figure 3 illustrates the obtained result for \(n_{1}= 34\) directions for the angles 90, 70 and 50. From the observation of the obtained results, we can note that the reconstruction of the diffusion signal from under-sampled data using the proposed method yields accurate results compared with the full HARDI acquisition (ground truth).

Table 1. Angular values estimation according to the \(n_{1}\) points set on simulated data.
Full size table
Fig. 3.
figure 3

FOD reconstruction on simulated data from: (a) the full set directions (b) the novel estimated signal for \(n_{1}= 34\).

Full size image

Results on Digital Phantoms. We used the software Phantom\( \alpha \)sFootnote 1 to generate the diffusion weighted images that consist of two fiber bundles crossing at various angles 50, 70 and 90. The data were acquired using 65 gradient directions with b = 3000 \(\mathrm{{s/mm}}^{2}\). Here, we select the region (voxel) corresponding to the intersection zone between the fibers. The first column of Fig. 4 illustrates the generated images and the selected voxels. In order to study the impact of the initial set of \(n_{1}\) points for the estimation of the FOD function, we variate the values of \(n_{1}\) from 26 to 40. Table 2 illustrates the obtained results. Figure 4 illustrates the obtained results using the digital phantoms for \(n_{1}= 34\). These results confirm the accuracy of the proposed method. In fact, the crossing angles values are too close to the ground truth. We show, therefore, with both synthetic diffusion signal and digital phantomas that using only 34 directions is enough to reconstruct the FOD function.

Table 2. Angular values estimation according to the \(n_{1}\) point set of digital phantoms.
Full size table
Fig. 4.
figure 4

FOD reconstruction on digital phantoms from: (b) the full set directions (c) the novel estimated signal for \(n_{1}= 34\).

Full size image

5 Conclusion

In this paper, we have proposed a new solution to accelerate the acquisition time of the HARDI technique. In fact, it requires a big number of gradient directions. We proposed in this work a novel geometrical approach to accurately reconstruct the HARDI signal from reduced number of samples. It is based on recovering non-acquired data obtained according to their neighborhood from a reduced set of data on the sphere of the q-space. The experiments were conducted on both synthetic data and digital phantoms simulating crossing fibers. We computed the FOD function to evaluate the accuracy of the proposed method to estimate the diffusion weighted signal. We have shown that as few as 34 directions are sufficient to obtain accurate FOD estimation.

In future works, we intend to test the accuracy of the proposed method in real data. It will be also interesting to study the robustness of the method under many types of noise.