Patent Publication Number: US-10324154-B2

Title: Generalized spherical deconvolution in diffusion magnetic resonance imaging

Description:
BACKGROUND 
     In general, magnetic resonance imaging (MRI) examinations are based on the interactions among a primary magnetic field, a radiofrequency (RF) magnetic field, and time varying magnetic gradient fields with gyromagnetic material having nuclear spins within a subject of interest, such as a patient. Certain gyromagnetic materials, such as hydrogen nuclei in water molecules, have characteristic behaviors in response to external magnetic fields. The precession of spins of these nuclei can be influenced by manipulation of the fields to produce RF signals that can be detected, processed, and used to reconstruct a useful image. 
     Techniques have been developed to perform MRI imaging sequences to generate diffusion information that might be indicative of tissue microstructure, and/or to evaluate tissue for various abnormalities, such as neural injuries, neural degenerative disorders, and so on. In diffusion imaging, information is encoded in both image space and diffusion space, the latter typically referred to in the art as q-space. Typically, q-space is utilized to detect the movement of hydrogen nuclei in water molecules. Diffusion encoding gradient pulses are utilized to characterize the three-dimensional movement of the water molecules within and between tissues—more specifically within individual voxels (each spatial location). The movement of the water molecules may be characterized as incoherent motion, which results from diffusion processes, and coherent motion, which results from tissue movement, such as during pulsation. 
     The data that is obtained from q-space acquisitions may be categorized into either diffusion spectrum imaging (DSI), single-shell acquisition (e.g., high-angular diffusion imaging (HARDI)), or multi-shell acquisitions, depending on the sampling technique. However, due to the high dimensionality of the acquisition techniques (e.g., 3D in spatial domain and 3D in q-space), evaluating the diffusion properties of water in vivo may result in extended scan times to acquire acceptable data for image reconstruction. 
     BRIEF DESCRIPTION 
     Certain embodiments commensurate in scope with the originally claimed subject matter are summarized below. These embodiments are not intended to limit the scope of the claimed subject matter, but rather these embodiments are intended only to provide a brief summary of possible forms of the invention. Indeed, the present disclosure may encompass a variety of forms that may be similar to or different from the embodiments set forth below. 
     In one embodiment, a magnetic resonance imaging method includes generating spatially resolved fiber orientation distributions (FODs) from magnetic resonance signals acquired from a patient tissue using a plurality of diffusion encodings, each acquired magnetic resonance signal corresponding to one of the diffusion encodings and being representative of a three-dimensional distribution of displacement of magnetic spins of gyromagnetic nuclei present in each imaging voxel. Generating the spatially resolved FODs includes performing generalized spherical deconvolution using the acquired magnetic resonance signals and a modeled tissue response matrix (TRM) to reconstruct the spatially resolved FODs. The method also includes using the spatially resolved FODs to generate a representation of fibrous tissue within the patient tissue. 
     In a second embodiment, a magnetic resonance imaging (MRI) system includes a primary field magnet that may place gyromagnetic nuclei within a patient into an equilibrium magnetization, a plurality of gradient field coils that may encode positional information into the gyromagnetic nuclei, a radiofrequency (RF) transmit coil that may perturb the gyromagnetic nuclei away from their equilibrium magnetization, a plurality of RF receiving coils configured to receive MR signals from the gyromagnetic nuclei as they relax to their equilibrium magnetization, and control and analysis circuitry programmed to acquire magnetic resonance signals using a plurality of diffusion encodings, each acquired magnetic resonance signal corresponding to one of the diffusion encodings and being representative of a three-dimensional distribution of displacement of magnetic spins of the gyromagnetic nuclei present in each imaging voxel, perform generalized spherical deconvolution using the acquired magnetic resonance signal and a modeled tissue response matrix (TRM) to reconstruct a spatially resolved fiber orientation distribution (FODs), and generate a representation of fibrous tissue within the patient using the spatially resolved FODs. 
     In a third embodiment, one or more tangible, non-transitory, machine-readable media storing instructions executable by one or more processors to cause a magnetic resonance (MR) imaging system to acquire magnetic resonance signals using a plurality of diffusion encodings, each acquired magnetic resonance signal corresponding to one of the diffusion encodings and being representative of a three-dimensional distribution of displacement of magnetic spins of the gyromagnetic nuclei present in each imaging voxel, perform generalized spherical deconvolution using the acquired magnetic resonance signal and a modeled tissue response matrix (TRM) to reconstruct a spatially resolved fiber orientation distribution (FODs); and generate a representation of fibrous tissue within the patient using the spatially resolved FODs. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       These and other features, aspects, and advantages of the present embodiments will become better understood when the following detailed description is read with reference to the accompanying drawings in which like characters represent like parts throughout the drawings, wherein: 
         FIG. 1  is a diagrammatical illustration of an embodiment of a magnetic resonance imaging (MRI) system configured to perform diffusion imaging in accordance with an aspect of the present disclosure; 
         FIG. 2  is a process flow diagram depicting an embodiment of a method for obtaining spatially resolved fiber orientation distributions (FODs) by performing compressed sensing on undersampled diffusion spectrum imaging (DSI) data followed by spherical deconvolution in accordance with an aspect of the present disclosure; 
         FIG. 3  is an embodiment of example fiber crossing plots comparing FODs generated from simulated DSI data according to the method of  FIG. 2  and orientation distribution functions (ODFs) generated by integrating over a set of weighted reconstructions of the diffusion propagator obtained via a Fourier transform of the same simulated DSI data; 
         FIG. 4  depicts example tractography images generated from spatially resolved FODs computed from two different data sets having comparable data acquisition times; one is produced from a single-shell high angular resolution diffusion imaging (HARDI) acquisition and the other one is produced from undersampled DSI data using compressed sensing and spherical deconvolution according to the method of  FIG. 2  in accordance with an aspect of the present disclosure; 
         FIG. 5  is a process flow diagram depicting an embodiment of a method for obtaining spatially resolved FODs and/or quantitative tissue maps from diffusion MRI (dMRI) data based on a generalized spherical deconvolution algorithm; 
         FIG. 6  depicts example fiber crossing plots comparing fiber peaks detected from FODs generated according to the method of  FIG. 5  and from a damped Richardson-Lucy algorithm (dRLA); 
         FIG. 7  is an embodiment of example tractography images generated from spatially resolved FODs generated according to the method of  FIG. 5  and the damped Richardson-Lucy algorithm (dRLA); and 
         FIG. 8  is an embodiment of example quantitative tissue maps computed using the data generated according to the method of  FIG. 5  and using the original data. 
     
    
    
     DETAILED DESCRIPTION 
     One or more specific embodiments will be described below. In an effort to provide a concise description of these embodiments, all features of an actual implementation may not be described in the specification. It should be appreciated that in the development of any such actual implementation, as in any engineering or design project, numerous implementation-specific decisions must be made to achieve the developers&#39; specific goals, such as compliance with system-related and business-related constraints, which may vary from one implementation to another. Moreover, it should be appreciated that such a development effort might be complex and time consuming, but would nevertheless be a routine undertaking of design, fabrication, and manufacture for those of ordinary skill having the benefit of this disclosure. 
     When introducing elements of various embodiments of the present embodiments, the articles “a,” “an,” “the,” and “said” are intended to mean that there are one or more of the elements. The terms “comprising,” “including,” and “having” are intended to be inclusive and mean that there may be additional elements other than the listed elements. Furthermore, any numerical examples in the following discussion are intended to be non-limiting, and thus additional numerical values, ranges, and percentages are within the scope of the disclosed embodiments. 
     Diffusion MRI (dMRI) may be used to evaluate white matter pathways in a patient&#39;s brain (e.g., via fiber tracking also known as tractography) and to generate quantitative tissue maps (e.g., mean diffusivity (MD), fractional anisotropy (FA), mean kurtosis, radial kurtosis, and/or neurite dispersion NODDI)). The dMRI data is generally acquired by repeating a dMRI acquisition with varying diffusion encodings that include directional (e.g., 3D) and weighting (e.g., b-value) information associated with the diffusion properties of the water. The diffusion encodings span a 3D space, also known as q-space. 
     Acquiring a complete data set to fill the q-space (3D space) and generate high resolution images may result in extended scan times. In general, patient movement is limited during MRI acquisitions (e.g., the patient may not be allowed to move). Therefore, the extended scan times may result in patient discomfort, and in certain situations the patient may be unable to remain still during acquisitions. Patient movement during MRI data acquisition may result in poor image quality. Accordingly, it may be desirable to decrease data acquisition times for dMRI techniques to decrease patient discomfort and movement during imaging. A technique for decreasing dMRI acquisition times may be to undersample the q-space data and using mathematical algorithms to transform the undersampled q-space data into a fully sampled data set (e.g., a data set that has not been undersampled during data acquisition), which may be referred to as a reconstructed data set, and reconstruct an image. However, due, in part, to the missing q-space data (e.g., undersampled data), images may include artifacts and have undesirable resolution. In accordance with present embodiments, model-based techniques that combine spherical deconvolution (SD) and compressed sensing are used to denoise and process q-space data (e.g., undersampled q-space data) to generate images similar to images generated from a full q-space data set. In this way, dMRI q-space data acquisition times may be accelerated without affecting image quality. Indeed, the present embodiments improve image resolution and/or quality compared to the conventional dMRI acquisition and processing techniques. 
     The embodiments described herein may be performed by a magnetic resonance imaging (MRI) system, wherein specific imaging routines (e.g., diffusion MRI sequences) are initiated by a user (e.g., a radiologist). Thus, the MRI system may perform data acquisition, data construction, and image synthesis. Accordingly, referring to  FIG. 1 , a magnetic resonance imaging system  10  is illustrated schematically as including a scanner  12 , scanner control circuitry  14 , and system control circuitry  16 . According to the embodiments described herein, the MRI system  10  is generally configured to perform MR imaging, such as imaging sequences for diffusion imaging. System  10  additionally includes remote access and storage systems or devices such as picture archiving and communication systems (PACS)  18 , or other devices such as teleradiology equipment so that data acquired by the system  10  may be accessed on- or off-site. In this way, MR data may be acquired, followed by on- or off-site processing and evaluation. While the MRI system  10  may include any suitable scanner or detector, in the illustrated embodiment, the system  10  includes a full body scanner  12  having a housing  20  through which a bore  22  is formed. A table  24  is moveable into the bore  22  to permit a patient  26  to be positioned therein for imaging selected anatomy within the patient. 
     Scanner  12  includes a series of associated coils for producing controlled magnetic fields for exciting the gyromagnetic material within the anatomy of the subject being imaged. Specifically, a primary magnet coil  28  is provided for generating a primary magnetic field, B 0 , which is generally aligned with the bore  22 . A series of gradient coils  30 ,  32 , and  34  permit controlled magnetic gradient fields to be generated for positional encoding of certain of the gyromagnetic nuclei within the patient  26  during examination sequences. A radio frequency (RF) coil  36  is configured to generate radio frequency pulses for exciting the certain gyromagnetic nuclei within the patient. In addition to the coils that may be local to the scanner  12 , the system  10  also includes a set of receiving coils  38  (e.g., an array of coils) configured for placement proximal (e.g., against) to the patient  26 . As an example, the receiving coils  38  can include cervical/thoracic/lumbar (CTL) coils, head coils, single-sided spine coils, and so forth. Generally, the receiving coils  38  are placed close to or on top of the patient  26  so as to receive the weak RF signals (weak relative to the transmitted pulses generated by the scanner coils) that are generated by certain of the gyromagnetic nuclei within the patient  26  as they return to their relaxed state. 
     The various coils of system  10  are controlled by external circuitry to generate the desired field and pulses, and to read emissions from the gyromagnetic material in a controlled manner. In the illustrated embodiment, a main power supply  40  provides power to the primary field coil  28 . A driver circuit  42  is provides power to pulse the gradient field coils  30 ,  32 , and  34 . Such a circuit may include amplification and control circuitry for supplying current to the coils as defined by digitized pulse sequences output by the scanner control circuit  14 , which in one embodiment may be a diffusion imaging module. Another control circuit  44  is provided for regulating operation of the RF coil  36 . Circuit  44  includes a switching device for alternating between the active and inactive modes of operation, wherein the RF coil  36  transmits and does not transmit signals, respectively. Circuit  44  also includes amplification circuitry configured to generate the RF pulses. Similarly, the receiving coils  38  are connected to switch  46 , which is capable of switching the receiving coils  38  between receiving and non-receiving modes. Thus, the receiving coils  38  resonate with the RF signals produced by relaxing gyromagnetic nuclei from within the patient  26  while in the receiving mode, and they do not resonate with RF energy from the transmitting coils (i.e., coil  36 ) so as to prevent undesirable operation while in the non-receiving mode. Additionally, a receiving circuit  48  is configured to receive the data detected by the receiving coils  38 , and may include one or more multiplexing and/or amplification circuits. 
     It should be noted that while the scanner  12  and the control/amplification circuitry described above are illustrated as being coupled by a single line, that many such lines may occur in an actual instantiation. For example, separate lines may be used for control, data communication, and so on. Further, suitable hardware may be disposed along each type of line for the proper handling of the data. Indeed, various filters, digitizers, and processors may be disposed between the scanner and either or both of the scanner and system control circuitry  14 ,  16 . By way of non-limiting example, certain of the control and analysis circuitry described in detail below, although illustrated as a single unit, includes additional hardware such as image reconstruction hardware configured to perform the data processing techniques described herein. 
     As illustrated, scanner control circuit  14  includes an interface circuit  50 , which outputs signals for driving the gradient field coils and the RF coil and for receiving the data representative of the magnetic resonance signals produced in examination sequences. The interface circuit  50  is coupled to a control and analysis circuit  52 . The control and analysis circuit  52  executes the commands for driving the circuit  42  and circuit  44  based on defined protocols selected via system control circuit  16 . Control and analysis circuit  52  also serves to receive the magnetic resonance signals and performs subsequent processing before transmitting the data to system control circuit  16 . Scanner control circuit  14  also includes one or more memory circuits  54 , which store configuration parameters, pulse sequence descriptions, examination results, and so forth, during operation. As noted above, in one embodiment, the scanner control circuit  14  may correspond to all or a part of a diffusion imaging module. 
     Interface circuit  56  is coupled to the control and analysis circuit  52  for exchanging data between scanner control circuit  14  and system control circuit  16 . In certain embodiments, the control and analysis circuit  52 , while illustrated as a single unit, may include one or more hardware devices. Further, in certain embodiments, system control circuit  16  may be a part of the diffusion imaging module. The system control circuit  16  includes an interface circuit  58 , which receives data from the scanner control circuit  14  and transmits data and commands back to the scanner control circuit  14 . The interface circuit  58  is coupled to a control and analysis circuit  60  which may include a CPU in a multi-purpose or application specific computer or workstation. Control and analysis circuit  60  is coupled to a memory circuit  62  to store programming code for operation of the MRI system  10  and to store the processed image data for later reconstruction, display and transmission. The programming code may execute one or more algorithms that, when executed by a processor, are configured to perform reconstruction of acquired q-space data to obtain both angular and radial information associated with incoherent motion (e.g., diffusion) of the gyromagnetic nuclei throughout the patient&#39;s tissue to generate spatially resolved fiber orientation distributions (FODs), quantitative maps, or other images of the patient&#39;s tissue. The spatially resolved FODs may be used to visualize fiber crossings or to generate fiber tractograms. For example, as discussed in further detail below, the spatially resolved FODs, as well as the quantitative tissue maps, may be generated based on one or more algorithms that employ generalized spherical deconvolution. 
     An additional interface circuit  64  may be provided for exchanging image data, configuration parameters, and so forth with external system components such as remote access and storage devices  18 . Finally, the system control and analysis circuit  60  may include various peripheral devices for facilitating operator interface and for producing hard copies of the reconstructed images. In the illustrated embodiment, these peripherals include a printer  60 , a monitor  62 , and user interface  64  including devices such as a keyboard or a mouse. 
     It should be noted that the MRI system described is merely intended to be exemplary only, and other system types, such as so-called “open” MRI systems may also be used. Similarly, such systems may be rated by the strength of their primary magnet, and any suitably rated system capable of carrying out the data acquisition and processing described below may be employed. Indeed, at least a portion of the disclosed methods may be performed by the system  10  described above with respect to  FIG. 1 . That is, the MRI system  10  may perform the acquisition techniques described herein, and, in some embodiments, the data processing techniques described herein. It should be noted that subsequent to the acquisitions described herein, the system  10  may simply store the acquired data for later access locally and/or remotely, for example in a memory circuit (e.g., memory  62 ). Thus, when accessed locally and/or remotely, the acquired data may be manipulated by one or more processors contained within an application-specific or general-purpose computer. The one or more processors may access the acquired data and execute routines stored on one or more non-transitory, machine readable media collectively storing instructions for performing methods including the diffusion imaging methods described herein. For example, the one or more processors may apply compressed sensing, separately or in combination with spherical deconvolution, to the acquired data. 
     One such method  80  is illustrated as a process flow diagram in  FIG. 2 . In particular, method  80  includes a number of steps configured to process q-space data (e.g., diffusion encoding data, 3D data) to characterize both angular and radial components associated with diffusion characteristics of water within the patient&#39;s tissue being imaged (e.g., nerve fibers). The method  80  begins with performing a diffusion MRI acquisition using undersampled diffusion spectrum imaging (DSI) (block  82 ). The DSI acquisition includes diffusion pulse sequences for acquiring data in the form of a set of signals from a three-dimensional (3D) region of interest (e.g., a tissue of the patient). By way of non-limiting example, the acquisition may include a spin-echo or double-spin-echo echo-planar imaging pulse sequence having diffusion-sensitizing gradient pulses. The pulse sequence produces a series of MR images (after complex Fourier transformation in k-space), with each image representing one point in q-space. The images constitute a complex diffusion-encoded data set containing angular and radial information relating to underlying tissue anisotropy and geometric restrictions in the diffusion process, respectively. The acquired signals at each of the q-space locations represent the three-dimensional (3D) displacement distribution of the spins of the gyromagnetic material (e.g., protons in water) in each voxel. The signals obtained according to block  82  may include signals corresponding to one of the diffusion encodings and may be representative of a three-dimensional distribution of displacement of magnetic spins of the gyromagnetic nuclei (e.g., the proton in the water) present in each imaging voxel. 
     As should be noted, the signals may be acquired using methods other than diffusion spectrum imaging (DSI), which acquires signals on a Cartesian grid. For example, the signals may be acquired using high angular resolution diffusion imaging (HARDI) or q-ball imaging, which represent subsets of q-space imaging with acquisitions on one spherical single shell (e.g., HARDI) or multiple concentric spherical shells instead of a full Cartesian grid. 
     The most comprehensive diffusion information is typically obtained using an acquisition that satisfies the Nyquist criterion. The acquisition according to block  82  will typically be of a high dimensionality—including three dimensions in the spatial domain (e.g., k-space) and three dimensions in the q-space domain, and in certain embodiments, may be performed at the Nyquist rate. As discussed above, the scan time for certain dMRI sequences may be relatively long, which may limit clinical in vivo applications of dMRI. To reduce dMRI scan time, in certain embodiments, less than all of the q-space data may be acquired. In other words, as an alternative acquisition method, the q-space data may be undersampled (e.g., sampled below the Nyquist rate). The acquisition trajectory employed in such embodiments may be any suitable trajectory, such as a non-uniform (e.g., a random or pseudo-random) pattern, a regularly sampled pattern of sparsely sampled locations, and so on. In one embodiment, the pattern may not be evenly distributed. Sampling of the q-space may be within any desired radius of the q-space. In certain embodiments, a center region of the q-space is sampled more densely compared to other non-central regions within the q-space. Examples of undersampling schemes that may be used for q-space data include, but are not limited to, uniform random distribution, Poisson disk, and Gaussian. As should be noted, the number of samples in q-space may be determined by the Nyquist sampling rate and image acquisition time. 
     In embodiments where fewer than the total available q-space locations are sampled, the method  80  may include reconstructing q-space encodings for unsampled q-space locations using compressed sensing reconstruction (block  84 ) so as to generate an ensemble average propagator (EAP) and to produce a set of uniformly distributed q-space data points using the inverse Fourier transform (block  86 ) on the EAP. Due to the random or pseudo-random nature of the undersampling pattern, compressed sensing may enable the reconstruction of missing q-space data. Compressed sensing may recover the missing q-space data by exploiting generic prior knowledge about the q-space data such as signal sparsity or minimal total variation (TV). However, in certain embodiments, parallel imaging reconstruction techniques could potentially be utilized for synthesizing missing spatial data. 
     As will be understood by those of skill in the art, compressed sensing is a technique in which data is both acquired and then reconstructed. In compressed sensing, it is desired that the image have a sparse representation in a known transform domain (such as the wavelet domain) and that aliasing artifacts due to q-space undersampling be incoherent in that transform domain (i.e., noise-like). That is, the data-sampling pattern is chosen so as to reduce coherency in the sparse domain. This incoherence may be achieved by randomly undersampling the q-space encodings. The undersampling of the q-space encodings may result in aliasing, and when the undersampling is random, the aliasing is incoherent and acts as incoherent interference of the sparse transform coefficients. 
     Accordingly, the sampling in accordance with block  82  may be a pattern that is amenable to reconstruction using a compressed sensing algorithm to produce the EAP used for fiber tracking and/or quantitative mapping. While any suitable compressed sensing algorithm may be used for reconstruction (e.g., estimation), in one embodiment, the reconstruction may be performed using a non-linear reconstruction scheme. By way of example, an L1-norm constraint in a cost function minimization method may be used to recover sparse transform coefficients that may be used to compute the EAP. In this way, the EAP may be reconstructed for an undersampled q-space data set. 
     As discussed above, the EAP may be reconstructed using compressed sensing of the q-space data may be retransformed. For example, an inverse Fourier transform may be applied to the EAP (block  86 ). The inverse Fourier transform retransforms the reconstructed EAP into a fully sampled q-space data set, which may be used for fiber tracking and/or quantitative mapping. In the inverse Fourier transform, at least a portion of the reconstructed q-space data (e.g., synthesized data) may be replaced with corresponding acquired q-space data (e.g., the undersampled q-space data generated according to block  82 ), thereby generating a fully sampled q-space data set. As should be noted, in certain embodiments, reconstruction and re-transformation of the undersampled q-space data (e.g., according to blocks  84 ,  86 ) before spherical deconvolution processing may be omitted, as discussed in further detail below with reference to  FIG. 5 . 
     Following re-transformation of the EAP, spherical deconvolution may be applied (block  90 ) to the re-transformed q-space data to produce the spatially resolved fiber orientation distributions (FODs) (block  92 ). For example, in certain embodiments, generating the spatially resolved FODs includes reconstructing the ensemble average propagator (EAP) using compressed sensing, transforming the reconstructed EAP back to a fully sampled DSI data set, and performing spherical deconvolution on the fully sampled DSI data set. Once reconstructed the spatially resolved FODs may be used for visualizing fiber crossings or for fiber tracking (block  98 ). Fiber tracts generated using the reconstructed spatially resolved FODs may have a better resolved fiber peaks compared to tracts generated using orientation distribution functions (ODFs), thereby increasing the accuracy of fiber crossings and fiber tractography generated from undersampled DSI data. 
     Generally, fast Fourier transform is performed in q-space to produce the EAP. The EAP is further processed by integrating over a set of weighted reconstructions of a diffusion propagator to produce the ODF. The spatially resolved ODFs may be used as an input for tractography. 
     However, the angular resolution of the ODFs may be affected by the weighting (e.g., b-value) information associated with the diffusion properties of the water. For example, a low b max  (e.g., less than approximately 3,000 seconds/square millimeter (s/mm 2 )) may result in blurred ODFs, and a high b max  (e.g., greater than approximately 6,000 s/mm 2 ) may have an undesirable signal-to-noise ratio, thereby decreasing the sensitivity of the ODFs for fiber tracking applications. Therefore, to improve the quality of the fiber tracts generated from diffusion spectrum imaging (DSI) q-space data, the reconstructed q-space data may be processed using spherical deconvolution algorithms to generate the spatially resolved FODs that may be used for high resolution fiber tracking rather than using ODFs. 
     The spherical deconvolution model that may be used to reconstruct the spatially resolved FOD is expressed as follows:
 
 S ( b,q )=Σ i   f   i   e   −λ     ∥     b(q     T     u     i     )     2     e   −λ     ⊥     b(1−(q     T     u     i     )     2     )   +v   GM   e   −λ     GM     b   +v   CSF   e   −λ     CSF     b   (EQ.1)
 
 S=[H,h   GM   ,h   CSF   ]*[f;v   GM   ;v   CSF ]  (EQ. 2)
 
where λ | , λ ⊥ , λ GM , and λ CSF  represent diffusivity constants for (anisotropic) white matter (parallel and orthogonal, respectively, to the white matter fiber direction) and for (isotropic) grey matter (GM), and cerebral spinal fluid (CSF), respectively, f i  is the wanted fiber orientation distribution (FOD) in the direction u i , v GM  and v CSF  are the wanted volume fractions of GM and CSF, respectively, b and q are the strength and the direction, respectively, of the diffusion weighting used to acquire the data S. EQ. 2 is the matrix-vector notation for the spherical deconvolution model. EQs. 1 and 2 may be used for multi-shell or Cartesian grid (DSI) acquisition and processing techniques. For single-shell data the model is simplified by removing the isotropic components.
 
       FIG. 3  illustrates reconstructed FODs, according to block  90 , and orientation distribution functions (ODF) for simulated q-space data (Gaussian mixture model, FA=0.8, D ∥ =1E−3 mm 2 /s, 11-cube DSI) at various b max  values, noise levels, and crossing angles. As illustrated, in  FIG. 3 , the fiber crossings generated based on the reconstructed spatially resolved FOD are equally or better defined for b max  values of 3,000-8,000 s/mm 2  compared to ODF, with the exception of fiber crossings at a b max  of 8,000 s/mm 2  (e.g., high b max ) and approximately 10% noise. Accordingly, it is now recognized that by combining the denoising associated with compressed sensing and the use of an FOD generated from spherical deconvolution, higher resolution fiber images may be obtained at a wide range of b values. 
     Moreover, experiments on a healthy volunteer were performed on a 3T GE MR750 clinical MR scanner (GE Healthcare, Milwaukee, Wis., USA) using a 32-Channel Head Coil and acquiring two single shot echo planar imaging (EPI) acquisitions (ASSET factor 2, 2.5 mm isotropic resolution), first an undersampled DSI acquisition (11 cube, R=4, b max =6,000 s/mm 2 ) and second a single-shell HARDI acquisition (150 directions, b=3,000 s/mm 2 ). As illustrated in  FIG. 4 , undersampled diffusion spectrum imaging (DSI) techniques using a combination of compressed sensing, according to block  84 , and spherical deconvolution, according to block  90 , to process the dMRI signals may improve fiber tractography resolution compared to high angular resolution diffusion imaging (HARDI) q-space data processing techniques. For example, images  110 ,  112  illustrate motor-tracts and optic radiation tracts, respectively, in a patient&#39;s brain tissue generated using the method  80 . Similarly, images  116 ,  118  illustrate motor-tracts and optic radiation tracts, respectively, in the patient&#39;s brain tissue generated via single-shell HARDI and spherical deconvolution techniques. As shown by arrows  120 , fiber angles and orientations are better resolved in tract images  110 ,  112  compared to tract images  116 ,  118 . 
     In addition to fiber tracking, the re-transformed q-space data (e.g., synthesized to fully sampled q-space data), generated according to block  86  may also be used to produce additional information associated with the diffusion of the water within the patient&#39;s tissue. For example, in certain embodiments, a higher order model fitting (block  98 ) may be used to fit the re-transformed reconstructed q-space data and generate the quantitative maps (e.g., mean diffusivity (MD), fractional anisotropy (FA), mean kurtosis, radial kurtosis, and/or neurite dispersion NODDI)) of the patient&#39;s tissue. 
     As discussed above, certain embodiments of the present disclosure may use any aspect of compressed sensing for generating spatially resolved FODs (e.g., according to block  84  of method  80 ). For example,  FIG. 5  illustrates a method  130  for analyzing and denoising q-space data to spatially resolve fiber orientation distributions (FODs) using a generalized spherical deconvolution algorithm. Method  130  begins with performing a diffusion MRI (dMRI) acquisition (block  132 ) to generate q-space data. The dMRI acquisition may be done using any suitable dMRI acquisition technique such as a diffusion spectrum imaging (DSI) acquisitions (e.g., Cartesian grid acquisition), a single-shell acquisition (e.g., HARDI), or a multi-shell acquisition. The sampling may be complete (e.g., satisfies the Nyquist criterion) or incomplete (to generate an undersampled data set). 
     Following acquisition of the q-space data, the method  130  applies a generalized spherical deconvolution algorithm (block  140 ) to the q-space generated in accordance with block  132 . The generalized spherical deconvolution algorithm may include a least squares formulation and, in certain embodiments, a compressed sensing regularization term. The compressed sensing regularization term may promote sparsity of the spatially resolved FODs. In certain embodiments, the compressed sensing regularization term may promote sparsity by minimizing the L1-norm associated with the spatially resolved FODs. Additionally, the compressed sensing regularization term may promote smoothness of the spatially resolved FODs. As discussed in further detail below, the smoothness may be promoted by minimization of the total variation (TV) of the spatially resolved FODs. For example, given the measured q-space data d for a particular voxel in a single-shell (HARDI) acquisition, the matrix notation for the standard spherical deconvolution task is as follows:
 
min f   ∥Hf−d∥s.t.f≥ 0  (EQ. 3)
 
where H is the fiber response function matrix and f is the wanted fiber orientation distribution (FOD). A non-negativity constraint is imposed to avoid negative fibers (e.g., physically meaningless information). The number of potential fiber directions (e.g., the size of f) is usually greater than the number of measurements (e.g., the size of d), resulting in an underdetermined system (e.g., the matrix H is non-quadratic). Systemically reducing the number of measured data points (e.g., undersampling d) may decrease determinacy and result in increased blurring of the FOD. However, the application of compressed sensing to the spherical deconvolution task may recover the original resolution (e.g., sharpness) of the FOD. Therefore, in certain embodiments, a compressed sensing (CS) regularizer may be applied to the spherical deconvolution to generate a generalized spherical deconvolution. The matrix notation for the generalized spherical deconvolution may be expressed as shown in EQs. 4.
 
min f Ψ( f ) s.t.∥Hf−d∥&lt;ε,f≥ 0  (EQ. 4)
 
where Ψ is a regularizer (e.g., the total variation (TV) of f, the 1-norm (∥⋅∥ 1 ) of f, or the 1-norm of a transform of f) and ε a noise dependent small scalar.
 
     EQs. 4 and may be solved by using an alternating direction algorithm capable of handling non-negativity constraints or an iterative shrinkage/thresholding algorithm (ISTA), respectively. For example, for Ψ=∥⋅∥ 1 , i.e. if f is assumed to be sparse, the alternating direction algorithm may be applied. 
     For Ψ being the Total Variation (TV) of f, the iterative shrinkage thresholding algorithm (ISTA) approach may be used yielding a smooth plateau of constant value, thereby separating background and distinct peaks in the FOD. Because f is defined on a surface of a single sphere, Ψ is adapted for non-Cartesian sampling. Accordingly, the generalized spherical deconvolution based on EQ. 4 promotes sparsity of f (e.g., fiber direction), and decreases spurious fiber directions compared to conventional spherical deconvolution techniques. 
     In other embodiments, the spherical deconvolution model may be generalized to include multiple diffusivity values for anisotropic components, isotropic components, or a combination thereof within each imaging voxel (e.g., 3D space). As such, the multiplication of the reconstructed FODs with the extended Tissue Response Matrix (TRM) may yield denoised q-space data due, in part, to allowing multiple degrees of freedom via the algorithm. As shown above with regard to EQs. 1 and 2, conventional spherical deconvolution algorithms use a single set of diffusivities (e.g., λ | , λ ⊥ , λ GM , λ CSF ). However, by incorporating multiple diffusivities into the spherical deconvolution algorithm the q-space data may be denoised, thereby increasing image quality for quantitative tissue mapping. The generalized spherical deconvolution algorithm that incorporates multiple diffusivities is expressed as follows in EQs. 5-7:
 
 S ( b,q )=Σ j,i   f   j,i   e   −λ     ∥,j     b(q     T     u     i     )     2     e   −λ     ⊥,j     b(1−(q     T     u     i     )     2     ) +Σ k   v   GM,k   e   −λ     GM,k     b +Σ l   V   CSF,l   e   −λ     CSF,l     b   (EQ. 5)
 
 S=[H   1    . . . H   j   H   GM   H   CSF   ]*[f   1    . . . f   J   ;v   GM   ;v   CSF   ]=:H   gen   f   gen   (EQ. 6)
 
min f     gen     ∥H   gen   *f   gen   −d∥s.t.f≥ 0  (EQ. 7)
 
where the subscripts j, k, and l are subscripts for different diffusivity coefficients for anisotropic and isotropic components, respectively. EQ. 7 may be solved using the compressed sensing regularization discussed above.
 
     Based on the generalized spherical deconvolution algorithms and matrices expressed in EQ. 4, the spatially resolved FODs may be used for fiber tracking (block  96 ) as discussed above. For example, single-shell simulations of fiber crossings were performed using a Gaussian mixture model (FA=0.8, D ∥ =1E−3 mm 2 /s, b=3000 s/mm 2 , 128 directions) with varying crossing angles (e.g., between 30° and 90°) and up to 10% additive Rician noise. The generalized spherical deconvolution algorithm (e.g., EQ. 4, Ψ=∥⋅∥ 1 , Alternating Direction Algorithm) and damped Richardson-Lucy-Algorithm (dRLA) according to EQ. 3 were used to compute FODs. Subsequently, fiber peaks were detected using a simple neighborhood comparison. 
     Experiments on a healthy volunteer were performed on a 3T GE MR750 clinical MR scanner (GE Healthcare, Milwaukee, Wis., USA) using a 32-Channel Head Coil single shot echo planar imaging (EPI) HARDI acquisition (150 directions at b=3,000 s/mm 2 , 1.8 mm isotropic resolution). Data was randomly undersampled with factor R=3 yielding a single-shell data set with 50 directions. Spherical deconvolution was performed on the fully sampled data using damped Richardson-Lucy-Algorithm (dRLA), e.g., EQ. 3, to generate a reference image. In addition, the undersampled q-space data was processed using both dRLA, EQ. 3, and the generalized spherical deconvolution expressed in EQ. 4 with Ψ being the Total Variation using the Iterative Shrinkage/Thresholding Algorithm (ISTA). 
       FIG. 6  illustrates an array of example fiber crossing plots  140  showing detected fiber orientation distribution (FOD) peaks of a cross-fiber simulation study that compared FOD peaks derived from the generalized spherical deconvolution expressed in EQ. 4, Ψ=∥⋅∥ 1 , alternating direction algorithm, and the conventional dRLA expressed in EQ. 3 at different noise levels (e.g., between 0% and 10%). As illustrated in the example fiber crossing plots  140 , fiber peaks generated using conventional dRLA techniques ( 142 ) lead to spurious peaks for noisy data (e.g., noise levels above 0%). However, the peaks computed using the generalized spherical deconvolution algorithm ( 146 ) recovered correct peak counts and directions for various crossing angles independent of the noise level. In addition, the generalized spherical deconvolution derived fiber peaks  146  also show resolved peak crossing at a fiber angle of 30° in the noiseless case (e.g., approximately 0% noise). Note, numbers at the top right of each fiber crossing plot  140  indicate the number of FOD peaks detected, where the number on the left represents the FOD peaks detected using dRLA and the number on the right represents the FOD peaks detected using the method  130 . For example, in the fiber crossing plot  140  for a noise level of 0% and a crossing angle of 30°, the number of FOD peaks detected using dRLA is 1 and the number of FOD peaks detected using the method  130  is 2. 
       FIG. 7  illustrates example fiber tracts using the derived TV-ISTA generalized spherical deconvolution algorithm expressed in EQ. 4, according to block  96 . Besides having higher tract density, images generated based on EQ. 4 (e.g., TV-ISTA images  150 ,  152 ) using undersampled q-space data have similar fiber tracts compared to images generated based on the conventional dRLA algorithm expressed in EQ. 3 (e.g., reference images  154 ,  156 ) using a fully sampled q-space data set. That is, the TV-ISTA generalized spherical deconvolution algorithm in EQ. 4 recovered similar fiber tract information using only one third of the data compared to conventional dRLA processing for motor tracts and fine structures of the SLF branching out to the cerebral cortex. Moreover, as shown in undersampled dRLA images  160 ,  162  generated based on EQ. 3, fiber tract information of undersampled q-space data using conventional dRLA is not clearly resolved. Therefore, data processing q-space data using the disclosed EQs. 4 improves resolution of fiber peaks and fiber tracking compared to conventional dRLA processing. As such, desirable fiber peak and fiber tracking resolution may be achieved using data sets with fewer q-space sampling points, thereby decreasing the overall acquisition time for dMRI. 
     In addition to fiber tracking, the FODs generated using generalized spherical deconvolution according to EQs. 5-7 may be used to denoise q-space data and generate quantitative tissue maps. For instance, the method  130  includes multiplying the spatially resolved FODs generated according to block  140  with an extended tissue response matrix (TRM) (block  168 ), to reconstruct q-space data. In this way, the spatially resolved FODs may be used to go back to the q-space data and/or retransform potentially undersampled q-space data to a fully sampled q-space data set used to compute quantitative tissue maps (block  170 ) based on an output of the multiplication step. Moreover, due the multiple diffusivity components incorporated into the generalized spherical deconvolution algorithm (e.g., EQ. 5-7), the resolution and accuracy of the quantitative maps generated from undersampled q-space data are comparable to quantitative maps generated from a fully sampled q-space data set. Therefore, by using the disclosed generalized spherical deconvolution algorithms (EQ. 5-7) the q-space data may be undersampled and acquisition times may be decreased without decreasing image quality. 
       FIG. 8  illustrates example quantitative tissue maps for mean diffusivity (MD)  172 , fractional anisotropy (FA)  176 , orthogonal kurtosis (orth K)  178 , mean kurtosis (mean K)  180 , and parallel kurtosis (par K)  184  using higher order model fitting on either the original q-space data or the q-space data reconstructed using the generalized spherical deconvolution algorithms expressed in EQs. 5-7 for undersampled and fully sampled q-space data sets. The quantitative maps  172 ,  176 ,  178 ,  180 ,  184  were compared to a ground truth reference image generated from a tenfold averaged fully sampled q-space data. As shown, the quantitative maps  172 ,  176 ,  178 ,  180 ,  184  generated based on the reconstructed q-space data according to block  168  are comparable to the ground truth maps and, have increased image quality compared to maps generated using the original data. As such, the generalized spherical deconvolution algorithms (EQ. 5-7) enable denoising of the q-space data such that quantitative tissue maps may be generated with fewer q-space data points, thereby decreasing dMRI acquisition times. 
     This written description uses examples to disclose embodiments of the invention, including the best mode, and also to enable any person skilled in the art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal languages of the claims.