The #1 largest and #1 most cited open-access publisher of neuroscience

Original Research ARTICLE

Front. Integr. Neurosci., 02 April 2013 |

Complete fourier direct magnetic resonance imaging (CFD-MRI) for diffusion MRI

  • Health Research, Arlington Innovation Center, Virginia Polytechnic Institute and State University, Arlington, VA, USA

The foundation for an accurate and unifying Fourier-based theory of diffusion weighted magnetic resonance imaging (DW–MRI) is constructed by carefully re-examining the first principles of DW–MRI signal formation and deriving its mathematical model from scratch. The derivations are specifically obtained for DW–MRI signal by including all of its elements (e.g., imaging gradients) using complex values. Particle methods are utilized in contrast to conventional partial differential equations approach. The signal is shown to be the Fourier transform of the joint distribution of number of the magnetic moments (at a given location at the initial time) and magnetic moment displacement integrals. In effect, the k-space is augmented by three more dimensions, corresponding to the frequency variables dual to displacement integral vectors. The joint distribution function is recovered by applying the Fourier transform to the complete high-dimensional data set. In the process, to obtain a physically meaningful real valued distribution function, phase corrections are applied for the re-establishment of Hermitian symmetry in the signal. Consequently, the method is fully unconstrained and directly presents the distribution of displacement integrals without any assumptions such as symmetry or Markovian property. The joint distribution function is visualized with isosurfaces, which describe the displacement integrals, overlaid on the distribution map of the number of magnetic moments with low mobility. The model provides an accurate description of the molecular motion measurements via DW–MRI. The improvement of the characterization of tissue microstructure leads to a better localization, detection and assessment of biological properties such as white matter integrity. The results are demonstrated on the experimental data obtained from an ex vivo baboon brain.

1. Introduction

Since the conception of mathematical models for the effect of the magnetic moment diffusion in nuclear magnetic resonance (NMR) experiments by Hahn (1950), Carrand Purcell (1954), and Torrey (1956), several methods have been proposed for analysis of diffusion-weighted (DW) magnetic resonance imaging (MRI) signal. These advancements endowed with the non-invasive, in vivo nature of the technique, have propelled the initial utilization of DW imaging measures, e.g., apparent diffusion coefficient in early detection of ischemia (Moseley et al., 1990; Baird and Warach, 1998), to many highly crucial areas in research and clinical imaging: for example in cancer diagnosis (Song et al., 2002; Turkbey et al., 2009; Xu et al., 2009), follow-up on treatment, pre- and post-operative assessment for different organs [e.g., fiber tracking (Conturo et al., 1999; Mori and van Zijl, 2002)] white matter integrity assessment (Budde et al., 2007; Correia et al., 2008) as in monitoring of neurological diseases such as multiple sclerosis (Song et al., 2002) and disorders (Ciccarelli et al., 2008) like schizophrenia (Seal et al., 2008; Voineskos et al., 2010) and Alzheimer's disease (Mielke et al., 2009), as well as neonatal development (McKinstry et al., 2002) and traumatic brain injury (Mac Donald et al., 2011).

In brief, diffusion weighted magnetic resonance imaging (DW–MRI) has become an indispensable and versatile technique playing an important role in several applications by its ability to estimate diffusion. The abundance of DW–MRI models is an indicator of room for improvement as well as the necessity for unification [see Özcan et al. (2012) for a detailed account of the partial differential equation (PDE) based adaptation's implications as well as a thorough mathematical analysis and a description of the background of existing methods].

DW–MRI's aim is to obtain measures and characterization of microstructure by investigating the diffusion process. Several methods and models have been proposed, all originating from the seminal work of Stejskal and Tanner (1965). Therein, under the influence of the additional motion sensitizing magnetic field gradients, the self-diffusion PDE of the magnetic moments is included in the Bloch PDE to model the attenuation in the DW–NMR spectroscopy signal. The result is the estimation of the scalar diffusion coefficient of the entire sample. In a sense, DW–NMR added another dimension, i.e., the magnetic moment motion, to the spectroscopic information even before the introduction of magnetic moment position later by the invention of MR imaging.

Accordingly, herein, DW–MRI model is naturally progressed to a higher dimensional construct that jointly presents magnetic moment position and motion. This is achieved by carefully re-examining the first principles of DW–MRI signal formation using particle methods in the spirit of the work of McCall et al. (1963). The mathematical model constructed in section 2.1 is specifically obtained for DW–MRI signal (rather than DW–NMR) by including all of its elements (e.g., imaging gradients) using complex values without taking the signal's magnitude.

The approach reveals that for an lmr-dimensional MRI slice, the DW–MRI complex valued signal that comes out of the scanner is the (lmr + 3)-dimensional Fourier transform of the joint distribution function of the number of magnetic moments (that are at a given position at the initial time) and their displacement integrals. In other words, the first lmr dimensions correspond to the usual MRI k-space with position information and the remaining three dimensions constitute the frequency space of displacement integrals. The values of imaging and motion sensitizing magnetic field gradient vectors together define in the (lmr + 3)-dimensional Fourier space the sampling points of the joint distribution function's Fourier transform. The distribution function is recovered by taking the Fourier transform of the complete data directly (i.e., without any scaling or use of magnitudes), giving the method its name: Complete Fourier Direct (CFD) MRI (Özcan 2010b).

2. Materials and Methods

2.1. Complete Fourier Direct MRI Signal Formation

The MR signal is generated by the vectorial sum of transverse magnetization of magnetic moments (Haacke et al., 1999):


By neglecting the effect of spin–spin relaxation, the evolution of the transverse magnetization of the ith magnetic moment is described in a standard manner by a rotating magnetization vector according to Bloch equations (see Appendix B):


Here, γ is the gyromagnetic ratio, the transverse magnetization vector, mi, is written in complex number form with mi(t0) denoting the initial magnetization tipped to the transverse plane,

Ωi=t0tG(xi, τ) · xi(τ)dτ(3)

describes the phase (when multiplied by γ) as a function of the magnetic field gradients G(x,t)3, and the position of the magnetic moment is xi3. The time-dependent position, xi, in Equation (3) affects the phase, Ωi, thereby also affecting the total signal in Equation (1).

Any kind of displacement (such as Brownian motion, molecular movement in biological tissue with different medium and obstacles, coherent motion or any combination thereof) is incorporated straight into the model, by modeling the position in a general and direct form herein without any stochastic assumptions [such as Markovian property used in Wedeen et al. (2005)] on the motion


where wi(t)3 represents the displacement of the magnetic moment from its initial position, xi(t0), [i.e., wi(t0) = 0]. The only physical requirement is the continuity of wi(t) since a magnetic moment cannot disappear at a given point and reappear at another.

The signal is calculated using Equations (1– 4) in reverse order. Following the motion described by Equation (4), the phase of the ith magnetic moment in Equation (3) during the digital acquisition period of the two dimensional imaging (lmr = 2) pulsed-gradient spin-echo (PGSE) sequence of Figure 1, is obtained after tedious but routine derivations [see Appendix B for a brief exposition of the derivations and Özcan (2012)] using the definitions of the variables in Figure 1 with G*3 denoting the magnetic field gradient vectors labeled as read-out, ro; phase encode, pe; slice select, ss; and diffusion, D. Omitting routine calculations for trapezoidal shapes for clarity, the derivation is carried out assuming ideal gradient amplifiers providing rectangular shaped gradient pulses. The initial time, t0, is the end time of π/2 radio frequency (RF) pulse when the magnetization is fully tipped to the transversal plane. The resulting phase of the transverse magnetization is a function of time, t, and, the imaging and diffusion gradients (see Appendix B):

Ωi(t, Gpe, Gss, GD)   =((ttacqΔtrw) GroΔtrw Gpe) · xi(t0)Δtrw Gss · xi(t0)(5)   +GD · Wd+((td4td3)(td2td1))GD · xi(t0)+Φπi(6)   +Gro· Wiacq(t)(Gro+Gpe+Gss) · Wirw.(7)

Figure 1. The pulsed-gradient spin-echo (PGSE) pulse sequence and the definition of the variables used in the calculations. SS–EX is for the slice select gradient during the excitation (π/2) pulse, RO for read out, PE for phase encode, SS for is the slice select gradient, Diff marks the motion sensitizing pulses, SS–PI is the slice select gradient during the π-pulse and ACQ stands for digital acquisition period. In practice, the MR pulse sequences implement the rewind (rw) gradients such that the amplifiers are turned on and off at the same times.

The second term in Equation (6) removes the injection of the initial position into the DW signal because of equal pulse duration times, δ = td4td3 = td2td1. The term Φπi describes the systematic phase (see Appendix B1) created by the π-pulse slice select gradient (SS–PI) and in this work it will be automatically taken out by the phase correction algorithm in section 2.2. Equations (6) and (7) incorporate three integrals of the displacement wi(t): (Wdi, Wacqi, Wrwi) corresponding to the displacement integrals for diffusion (d), analog to digital conversion acquisition (acq), and initial rewind (rw) gradient time periods, respectively

Wid=td3td4wi(τ) dτtd1td2wi(τ) dτ,(8)Wiacq(t)=tacqtwi(τ) dτ, Wirw=t0trwwi(τ) dτ.(9)

First term in Equation (5) is the definition of the k-space in regular MRI, kmr2:

kmr(t, Gpe)=(ttacqΔtrw) GroΔtrw Gpe(10)

with the additional term,

φslice=Δtrw Gss · xi(t0),(11)

which is constant because the slice select axis component of xi(t0) is the slice position1.

Without loss of generality, by adopting the imaging coordinate frame defined by the directions of the read-out, phase encode and slice select gradients, Gro = [gro1, 0, 0], Gpe = [0, gpe2, 0], Gss = [0, 0, gss3], time and gpe2 become functions of kmr = [kmr1, kmr2, kmr3] using Equation (10):

t=kmr1/gro1+tacq+Δtrw and gpe2=kmr2/Δtrw.

Accordingly, Wacqi(t) becomes a function of kmr1 and the coefficients of Wrw in Equation (7) are written as a vector which is an affine function of kmr:

krw=[gro1, kmr2/Δtrw, gss3].(12)

Consequently, the phase, Ωi, of Equation (3) is expressed concisely by defining kD = GD:

Ωi=kmr · xi(t0)+kD · Wid+krw · Wirw        +Wi, 1acq(kmr1) gro1+φslice,(13)

reflecting the effect of initial position and displacement integrals on the phase2 of each magnetic moment. Since ϕslice is constant for all i, it is taken out of Equation (13) with the appropriate rotation of the magnetization coordinate frame on a slice by slice basis.

Finally, assuming that all of the magnetic moments have the real valued initial magnetization mi(t0) = m0, Equation (1) can be re-written using Equation (13) to reveal a Fourier relationship,


A more efficient way to evaluate the sum in Equation (14) is first to group the magnetic moments with the same position-displacement properties and then to count the numbers elements in the groups.

Definition 1. The joint position-displacement integral distribution function, Ptotalcfd(x, W), is defined as the number of magnetic moments with the initial position x3 at time t0, possessing the displacement integral values of W=(Wd,Wrw,W1acq)7.

The signal in Equation (14) is calculated by integrating over the whole position-displacement space (absorbing m0 into Ptotalcfd for ease of notation):


Equation (15) is by definition the Fourier transform of Ptotalcfd with non-linearities added by Wacq1.

Among the elements of W, the focus is on the most descriptive MRI observable, Wd. Its distribution is obtained by marginalizing Wrw from Equation (15)

Pcfd(x, Wd)=Pcfdtotal(x, Wd, Wrw) dWrw                  =(kmr, kD)1{Scfd(kmr, kD, 0)}.(17)

However, the affine dependence of krw in Equation (12) makes it impossible to fix krw = 0 and to sample in (kmr, kD, 0) subspace. The following example demonstrates how the affine dependence affects the measurements by using a two dimensional Gaussian function, exp(−(k21 + k22)), with the “undesirable” variable k2 sampled on a line k2 = a k1 + b:


The result is complex valued in comparison to the real valued Fourier transform of exp(−k21) which can be obtained by setting a = b = 0.

2.2. Phase Corrections for the Estimation of Pcfd

In addition to inherent affine dependence and non-linearities, different experimental factors, noise, hardware imperfections etc., affect the DW–MR signal adversely. CFD–MRI addresses these issues by adopting a pivotal, physically meaningful standpoint originating from the following Fourier transform property (Bracewell, 2000):

(real valued function) ⇆ ℱ ⇆ (Hermitian symmetric function).

Accordingly CFD–MRI reconstruction is based on the following:


Furthermore, an immediate implication of the transform property and the Hermitian symmetry of Scfd is that theoretically, taking its magnitude before Fourier transforming will result in a symmetric real valued distribution under ideal conditions. As noted above, in practice the experimental Scfd is never Hermitian symmetric resulting in an asymmetric magnitude. Consequently, the magnitude's Fourier transform used in existing methods (Callaghan, 1991; Wedeen et al., 2005), results in a complex valued (Hermitian symmetric) distribution function. The difficulty of a physical interpretation forced those methods to take the magnitude of the transform as well to obtain a real valued function.

Herein, in order to obtain a real valued Pcfd, the re-establishment of Hermitian symmetry in the signal during the computation of the inverse transform for Equation (17), is realized by phase corrections. The strategy is similar in principle to the correction of the kmr-space center's (echo time) shift during the read-out period of acquisition in MR imaging. The resulting linear phase shift in the physical read-out axis uniformly and systematically appears in all of the data. The shifts in both phase-encode (e.g., due to sample shaking) and read-out directions are corrected by first determining from the Fourier transform in kmr,

Ikmrcomplex((xro, xpe), kD)=kmr1{Scfd(kmr, kD, krw)},(18)

the angle (∠ Ikcomplexmr) from the signal regions along the center lines of each physical direction at kD = 0,

rro(xro)Ikmrcomplex((xro, 0), 0), rpe(xpe)Ikmrcomplex((0, xpe), 0).(19)

The phase corrections are then applied systematically at each value of kD (see Figure 3):


The Fourier transform in the remaining variables,

Pcfd(x, Wd)=kD1{Ikmr(x, kD)},(21)

is evaluated sequentially in each kD-dimension with the aim of re-establishing the Hermitian property, Ikmr(x, −kD) = Ik*mr(x, kD), using the following steps.

CFD Phase correction algorithm:

  1. Pick a pixel at location xc, preferably near the center of the image where tissue or a good signal area is located.
  2. Starting from the first direction, l = 1 of the kD space calculate the phase on the line passing through the origin (i.e., [kD1, 0, 0], [0, kD2, 0], [0, 0, kD3], respectively for l = 1, 2, 3), e.g., ∠Ikmr(xc,(kD1, 0, 0)).
  3. Investigate the plot of the phase versus kDl. Pick as many as possible consecutive values of kDl near 0 without sudden changes to assure high signal to noise value.
  4. Construct a polynomial of degree m (with m less than the number of points) that approximates the phase at the selected points. The polynomial's constant term must be set to be 0 to guarantee that Ikmr(xc, 0) remains unchanged. For example, for the first direction, at selected values of kD1 the polynomial looks like
    Ikmr(xc,(kD1,0,0))rD1(kD1)am(kD1)m  +am1(kD1)m1++a1kD1(22)
    as demonstrated in Figure 2.
  5. Apply the same phase correction systematically to the entirety of the data along the lth direction at each of the other dimensions, at all of the pixel locations. For example in the first direction, kD1, update Ikmr to be equal to yes for all kD2, kD3 and x.
  6. Repeat steps 2–5 for the remaining directions.

Figure 2. Özcan (2012): Top row: The Nyquist plots of uncorrected (left) and corrected (right) Scfd obtained from the experimental data described in section 2.4. The plots show data acquired at each diffusion gradient value kD1 on the complex plane. Bottom row: The magnitude and phase plots of the data. Uncorrected data (bottom row left, second column) exhibit a linear phase shift around 0 frequency, indicative of coherent motion. After the phase corrections obtained using the polynomial 0.266kD1 estimated from the points kD1 = ‒6, 0, 6, 12 Gauss/cm, the magnitude is unchanged but the signal's imaginary part is smaller for the corrected values visible by the difference between the vertical axis spans of Nyquist plots and the phase plots.

The algorithm transfers the signal to the real channel by preserving its energy as the phase corrected spin-echo image without diffusion gradients, Ikmr(x, 0). The distribution of the magnetic moments with low mobility in all three directions [i.e., Pcfd(x, 0)] shows the result of the transfer in Figure 3 for the sample described in section 2.4.


Figure 3. On the left, the fixed baboon brain images acquired in the anatomical-coronal plane. On the top row, the real and imaginary parts of Ikmr(x, 0) and on the bottom row, Pcfd(x, 0) are displayed. The imaginary part is approximately 10% of the real part in both cases. On the right, full representation of Pcfd with isosurfaces (P¯cfd=0.17, see section 3.1) around CC and EC junction. Starting from left bottom going clockwise, the sample pixels are from cerebro-spinal fluid (CSF), CC, white matter (WM) and CC, and EC junction, respectively (see also Figure 5).

In Pcfd(x, 0), areas with high level of organization inducing low mobility, such as the corpus callosum (CC), the external capsule (EC), the mid-brain and the pons, appear brighter. The image is not an anisotropy map, e.g., mineral oil would appear brighter than water due to a smaller diffusion coefficient despite both liquids being isotropic. Spin-echo image is more blurred because it is a low pass filtered version of Pcfd:

Ikmrcomplex(x, kD)=kmr1{Scfd(kmr, kD, krw)}Ikmr(x, 0)                          =2πγPcfdtotal(x, W) dW(23)

(see Appendix C).

CFD phase correction algorithm outperformed the fitting of the phase values up to the fourth degree multinomials in 3. The reasons behind this outcome, which will provide information about DW–MR signal artifacts, as well as inclusion of different functions for corrections will be investigated in the future.

2.3. CFD-MRI Sampling and Windowing

Whereas the standard MRI field of view (FOV) calculations (Haacke et al., 1999) are used for kmr-space, the infinite bandwidth in kD-space due to Pcfd's finite support in Wd-space (originating from finite length displacements) falls beyond the reach of the gradient hardware's limits for small diffusion gradient duration and separation times (δ and Δ, respectively in Figure 1). Even with a powerful gradient system, a large magnitude of kD causes substantial signal uncertainties due to an increasing performance deterioration as the power requirements push the hardware to its limitations.

With such a hardware constraint, in order to reduce ripple effects caused by truncation, Pcfd's bandwidth (i.e., Scfd's support) is shrunk by increasing δ and Δ causing the dispersion (covariance) of the displacement integral Wd (and therefore Pcfd's support) to increase. This is directly visible in the special case of Brownian Motion characterized by the diffusion tensor D in Figure 4. Pcfd and Scfd are zero mean Gaussians with covariances, respectively equal to (see Appendix A)


Figure 4. Effects of fitting the bandwidth within the field of view (FOV) in the Fourier space on the left, with the covariance of the Gaussians given as D^=(δ2(Δδ/3)D)1. The sampling is done in the fixed interval, [‒5, 5], of the Fourier Space followed by reconstruction on the right using the discrete Fourier transform. Theoretically this is equivalent to convolution with the sinc function (see Brigham, 1988) in physical space. The ripple effects created by sinc lobes diminish on the left as the Gaussian falls into the FOV with increasing δ, Δ but constant D.

E[Wd(Wd)T]=bt D and (bt D)1 where btδ2(Δδ/3)(24)

(see Özcan (2009, 2010a) because the Fourier transform of a Gaussian with a covariance matrix D^ is also a Gaussian with covariance D^1:

{exp((Wd)TD^1 Wd)}~exp(kDT D^ kD).(25)

The procedure is graphically displayed in Figure 4 also emphasizing the effect of ripples on the small values of Pcfd which are especially important in revealing microstructure as explained in section 3.1.

The second sampling criterion is an appropriate sampling rate i.e., sufficient number of points in kD-space to prevent aliasing artifacts on Pcfd. This is constrained by the time available for acquisitions as each point in kD-space requires the scan time of the entire kmr-space.

2.4. Experimental Setup and Analysis Methods

A fixed baboon brain immersed in 4% paraformaldehyde was used for the experiments. The primate was prematurely delivered on the 125th day and sacrificed on the 59th day after delivery. All animal husbandry, handling, and procedures were performed at the Southwest Foundation for Biomedical Research, San Antonio, Texas. Animal handling and ethics were approved to conform to American Association for Accreditation of Laboratory Animal Care (AAALAC) guidelines. Further details of the preparation are described in Kroenke et al. (2005).

The experiments were carried out on a 4.7 Tesla MR scanner (Varian NMR Systems, Palo Alto, CA, USA) with a 15 cm inner diameter gradient system, 45 Gauss/cm maximum gradient strength and 0.2 ms rise time using a cylindrical quadrature birdcage coil (Varian NMR Systems, Palo Alto, CA, USA) with 63 mm inner diameter. CFD–MRI data were obtained using the standard pulsed-gradient spin-echo multi-slice sequence. The kmr-space was sampled to result in images of 128 × 128 pixels with a FOV 64 × 64 mm2 and 0.5 mm slice thickness. The kD-space was sampled in a uniformly spaced Cartesian grid in a cube [−30, 30 Gauss/cm]3 with 6 Gauss/cm sampling intervals at each dimension resulting in 11× 11 × 11 voxels. The repetition time TR = 1 s, echo time TE = 56.5 ms, diffusion pulse separation time Δ = 30 ms and diffusion pulse duration δ = 15 ms were used.

The data were transferred to a two quad core 2.3 GHz Intel Xeon® cpu and 8 GB memory Dell Precision Workstation 490 running Windows XP® 64-bit operating system. The DWI data were placed in a 5-dimensional array in the computer memory and the discrete Fourier transform (DFT) was computed along with the phase corrections. In-house Matlab® (Mathworks, Natick, MA, USA) programs were used for all of the computations and to display the graphics and maps.

3. Results

3.1. Visualization of the CFD Distribution

The joint distribution's high-dimensionality [e.g., two dimensions for position in regular MR images (lmr = 2), plus three more for displacement integrals] creates a visualization challenge which is addressed herein by using Pcfd(x, 0) as the background image. Furthermore, the isosurface3 of normalized Pcfd,

P¯cfd(x, Wd)=Pcfd(x, Wd)Pcfd(x, 0)

is overlayed on the pixel at location x, as in Figure 3. For the sake of an objective assessment, the isosurfaces are defined using a common level value c (0 < c ≤ 1),

{Wd3:P¯cfd(x, Wd)=c},

over all locations. The key point is the choice of an appropriate c-value that will reveal the outskirts of Pcfd corresponding to the small number of “scout” magnetic moments that travel further away thereby portraying the microstructure of the environment. In summary,

  1. Too high values do not provide enough structural information (see first rows in Figure 5).
  2. The appropriately informative value depends on the properties of the motion (thus of the microstructure) at a given location (compare columns of Figure 5, right side).
  3. Too low values force the isosurfaces to become extremely noisy (see last row of Figure 5).

As the motion in highly organized tissue is less dispersed (i.e., a smaller support for Pcfd which implies a larger support for Scfd thereby causing bigger truncation effects), increasing the diffusion gradient times δ and Δ in section 2.3 will create almost “flat” displacement integral distributions at an isotropic medium like CSF. In this case, the small valued distribution [caused by constant integral value ∫ Pcfd(x, Wd) dWd = number of particles] is susceptible to noise, creating noisy isosurfaces of Figure 5 for all level values. In contrast, for the experiments conducted with an isotropic (water) phantom at much lower δ and Δ values, the isosurfaces were spheres for a wide range of c-values (not shown).


Figure 5. On the left, one dimensional graphical representation of the choice of c-value. The drawing below the horizontal axis displays the structure of the sample in the infinitesimal volume element [‒dx, dx] as seen by the magnetic moments from their initial position. The right side of the sample has less restrictive properties as on the boundary of tissue with a liquid, such as CC and cerebrospinal fluid (CSF). The noiseless isosurfaces consist of two points shown as dots on the graph. Low c-values correspond to the magnetic moments with longer travel paths providing more structural information than high c-values. However, too low values create noisy and disconnected isosurfaces represented with more than two points on the drawing. On the right, isosurfaces from different pixels in the baboon brain marked in Figure 3 demonstrate the effect of c-value on the information content from top = insufficient to bottom = noisy (see the text for the CSF column).

Figure 6 presents different isosurfaces that elucidate the tissue structure on the pons and the mid-brain of the fixed baboon brain sample described in section 2.4. The tracts on the left and right side of the mid-brain are visible with the ellipsoidal looking isosurfaces. Five isosurfaces from the same row of pixels marked on Figure 6 are displayed on Figure 6 corresponding to two different c-values. The green isosurfaces with larger level values are smoother and less informative than the red ones with smaller c-values. Different viewpoints at each row of Figure 6 emphasize that the isosurfaces are 3-D objects. The figure demonstrates one of the challenges of presentation: the displacement integral values must be considered in 3 to grasp the complete information offered by CFD–MRI.


Figure 6. The isosurfaces (P¯cfd=0.14) on the pons and the mid-brain. Each of the boxes indicate an isosurface presented, respectively on the right. Each column presents the same isosurface from different view angles. The dots on top of the frames are placed for orientative purposes. The surfaces are not necessarily ellipsoids and they have mostly an asymmetric structure. The outer red surface is the set P¯cfd=0.14 and the green surface is P¯cfd=0.21. The red surface envelops the green one.

Overall, the isosurfaces are not constrained to given forms like Gaussians, spherical harmonics or to any expansions. In fact, they are typically not even symmetric. They are structureless, general and direct.

Isosurface visualizations constitute only one method to present the high dimensional information obtained from CFD–MRI. Another example is the dimension reduction by means of computing the so called orientation distribution function (ODF) (Wedeen et al., 2005) obtained from radial integrals. However, for CFD–MRI the ODF raises the concern of inadequately presenting the “outskirts” of the Pcfd because the dispersion of the outgoing rays shown in Figure 7 jeopardizes the inclusion of the values further away from the origin (see also Figure 6).


Figure 7. An example of a set of rays from the origin along which the distribution function is integrated to obtain orientation distribution function. As the rays disperse with increasing distance from the origin the points describing the displacement of a smaller number of particles contributes less and less to the numerical integration due to sparse sampling on both surfaces. The encapsulation of the isosurfaces on the right with larger level values by the smaller ones in Figure 6 shows the information that would be missed with numerical radial integration. The utilization of isosurfaces is more informative as discussed in Figure 5.

New methods, additional elaborate schemes such as color coding for representation of three dimensional functions aimed at displaying relevant microstructural information of CFD are left for future studies.

4. Discussion

4.1. Comparison with Existing Methods

From a fundamental point of view, guided by the microstructure that surrounds them, the molecules are displaced due to thermal energy whether they are in the scanner or not. All existing DW–MR methods are designed with the same goal in mind: the reconstruction of the propagator4 that describes the displacement of the magnetic moment from the DW–MR signal.

However, as CFD–MRI demonstrates, from a systems science perspective the MRI scanner acts as a time-delay linear system with the input wi [displacement in Equation (4)], and the output Wi [displacement integral in Equation (8)], in Figure 8. The parameters Δ and δ define the delay and filter parameters. Special attention is paid in CFD–MRI to isolate the Fourier variable, kD = GD from these parameters in contrast to the q-space variable (Callaghan, 1991): q = (2π)−1 γ δ GD and the b-value of DTI: b = γ2GD22 δ2 (Δ − δ/3) (see Appendix A).


Figure 8. The DW–MRI signal from a signals and systems perspective (Özcan, 2010b) describing the MRI observables of the diffusion phenomenon. The input is the magnetic moment displacement wi and the output is the displacement integral, Wdi, defined in Equation (8).

The inverse problem of obtaining the propagator from the distribution of the displacement integrals is singular because of the one-to-many relationship between the displacement and its integral:

Wid=td3td4wi(τ) dτtd1td2wi(τ) dτ      wi(t),(26)

because all of the displacements with the same low frequency content in time are mapped to the same displacement integral value. This statistical accumulation prevents the determination of the propagator in a general environment from the distribution of displacement integrals5).

Existing methods' attempts to estimate the propagator relies on the narrow pulse approximation by assuming negligible pulse duration (δ = td2td1 = td4td3 in Figure 1) compared to pulse separation time, Δ, i.e., (δ << Δ) ⇒ (δ → 0) specifically in (Δ−δ/3) (Callaghan, 1991, p. 342). Under this short integration time assumption, in Wedeen et al. (2005) it is further argued that the approximation, Wdi ≈ (xi(td3) − xi(td1)) δ = (wi(td3) − wi(td1)) δ, is plausible. Although by the intermediate value theorem and the sample path continuity of the Brownian motion (Shiryaev, 1995), the values of each integral in Wdi are attained at a time point within the respective integration intervals, [td1, td2] and [td3, td4], the nowhere-differentiability of the sample paths (Shiryaev, 1995) implies that the intermediate time points satisfying the equality are not fixed as td1 and td3, but are themselves random variables. In consequence, without the inference of the displacements at fixed time points the propagator cannot be reconstructed.

Moreover, elaborate derivations carried in Wedeen et al. (2005) to model the propagator as a conditional probability, p(x(td3)|x(td1)), describing a Markovian process raise concerns specially in environments such as biological tissue since the particle's past collisions with microstructure guides its future displacements. In fact, while it violates the conditions of Wiener process (see Appendix A), this displacement memory provides the inference of the microstructure by way of affecting the displacements and consequently their displacement integral distributions. Accordingly, in CFD–MRI is indifferent to memory properties by modeling Pcfd as a joint distribution function of random variables instead of the conditional probability of a stochastic process.

A summary of CFD–MRI's detailed comparison with existing methods (Özcan, 2011; Özcan et al., 2012) is presented below for completeness of exposure. Namely, there exits two avenues for the path from DW–NMR to DW–MRI in the literature:

1. Model matching methods initiated by diffusion tensor imaging (DTI) (Basser et al., 1994; Mattiello et al., 1994) and further expanded with high angular resolution DW imaging (HARDI) (Frank, 2001), composite hindered and restricted model of diffusion (CHARMED) (Assaf and Basser, 2005), diffusion orientation transform (DOT) (Özarslan et al., 2006), two versions of the generalized DTI (GDTI) (Özarslan and Mareci, 2003; Liu et al., 2004), and diffusional kurtosis imaging (DKI) (Jensen et al., 2005).

2. Spectral methods originating from Callaghan's q-space (Callaghan et al., 1988) followed by the diffusion spectrum imaging (DSI) (Wedeen et al., 2005) and Q-ball imaging (Tuch, 2004).

With the exception of the GDTI presented in Liu et al. (2004) [see also the discussion in Özarslan et al., (2009)], all of the DW–MRI methods estimate symmetric quantities. The model matching methods project the data onto symmetric structures, such as ellipsoids in DTI or spherical harmonics in HARDI. The spectral methods use the magnitude of the signal in the Fourier transform resulting in symmetric functions (see section 2.2). It is difficult to imagine that molecular motion in a biological environment populated with different types of fluids, barriers and tissue would be symmetric at any given location, e.g., at the fiber junctions. Symmetry or lack of it ought to be determined by the data free of any constraints imposed by the model as in the implementation of CFD–MRI's unconstrained structure.

The Fourier relationship between signal and joint distribution function provides a complete understanding of model matching methods. The methods start by applying DFT to the data in the first lmr (imaging) dimensions. Thus, in the analysis of model matching methods the first lmr independent variables are the physical location. The three remaining untransformed variables are the independent variables of the Fourier reciprocal of displacement integral space, i.e., they are in the Fourier domain. The goal of model matching methods is to fit the preferred model to the displacement distribution function's Fourier transform, sampled at the (vector) values defined by the diffusion sensitizing gradients. In the case of Brownian motion this mixed variable (physical and frequency variables) approach is applicable because the diffusion coefficient D can be directly estimated from the Fourier domain by Equations (24) and (25). The mixed space, which works well for DW–NMR signal peak attenuation, is translated to DW–MRI at each position x by

|Ikmrcomplex(x, kD)|=|Ikmrcomplex(x, 0)|exp(H(kD))(27)

where Ikcomplexmr is given in Equation (18) and the function H defines the model, e.g., the quadratic form of DTI, spherical harmonics of HARDI or higher tensor expansions of GDTI [see Özcan et al. (2012) for a detailed exposition]. In a sense, these methods' aim could be summarized as expanding the Fourier portion of the mixed signal space. The basic example with a single term in the expansion is DTI for which:

H(kD)=γ2 bt kDT D kD(28)

where the calculation of bt from the PDE approach in Özcan (2010a) and covariance of displacement integrals in Appendix A resulted in the same value: δ2 (Δ − δ/3). In CFD parlance, by the Fourier relationship between Gaussians in Equation (25), the diffusion quadratic form, D, is estimated in kD-space, without recourse to a Fourier transform because of its direct appearance in Equation (28) rather than its inverse, D−1, in the Gaussian of motion space. The coefficient matrix defined by carefully selected vectors in kD-space that satisfy the invertibility conditions (Özcan, 2005) is used for the estimation of D in the linear algebraic framework of symmetric matrices (Papadakis et al., 1999; Özcan 2010a) [also refer to Özcan (2010a) and Özcan (2012) for the correspondence with the B-matrix formulation of Basser et al., (1994)].

The magnitude-based Fourier relationship presented in q-space methodology (Callaghan et al., 1988) is the origin of spectral methods. In Callaghan's book (Callaghan, 1991), parallel to the historical development of DW models, the theory is first developed for NMR experiments [see Callaghan, (1991, Chap. 6)] using polarized neutron scattering analogy. However, the translation from NMR to MRI is presented [see asserting without proof that the imaging and displacement portions of the signal are separable [see Callaghan, (1991, Chap. 8, pp. 440)]. The derivations of section 2.1 demonstrate that this is not the case.

In addition, by the affine dependence of krw on kmr CFD derivations show that the inseparability partially accounts for the non-Hermitian nature of the Scfd. Taking the magnitude of the DW–MRI signal, as in the case of DSI (Wedeen et al., 2005), does not count as a phase correction. Figure 9 demonstrates that by preserving Hermitian property, CFD–MRI captures correctly the crossing fibers at the junction of the CC and EC.


Figure 9. The comparison of Pcfd (Özcan, 2011) (first row) with the diffusion spectrum imaging orientation distribution function (DSI–ODF) (second row). Both functions are calculated from the same data on the right junction of the corpus callosum (CC) and the external capsule (EC), specifically from the pixels marked on the Figure 3. The isosurface (P¯cfd=0.17) captures the asymmetric structure of the fiber crossings while the ODF is constrained to be symmetric for all of the pixels. Note that in CSF, ODF detects structure which is not present in reality as indicated by CFD.

4.2. Conclusion and Future Studies

In the biomedical imaging modalities' grand aim of biomarker capability establishment, the discovery path for CFD–MRI passes through the distribution function:


With Pcfd in the middle, both sides of the path present themselves with important challenges.

First and foremost, in DW–MRI, the displacements without reference to initial positions [see Equation (6)] prevent the inference of microstructure position. For example, the distribution function of the biological phantom (Özcan et al., 2011) constructed with two crossing rat trigeminal nerve fibers is always in the form of two crossing bars across the origin regardless of the nerves' position as long as their relative angle is kept the same. Also in the same phantom, the agar gel (isotropic component) appears as a sphere around the origin of Pcfd domain without the possibility of identifying its location. As the distributions from various types of microstructural components accumulate around the origin, the discrimination level of overlaps, more prominent with increasing biological tissue complexity, directly defines the sensitivity and specificity for microstructure changing pathologies. The important goal is the assessment of the theoretical aspects of the distribution function in order to understand whether it can detect in a timely manner, e.g., before significant disease progression, those changes. The determination of biophysical conditions behind the asymmetry (see also Özarslan et al., 2008; Özarslan, 2009) in the distribution functions is also part of the same goal.

However, the absence of analytical descriptions for Pcfd even in simple environments requires the investigations to be conducted with numerical simulations (Özcan et al., 2011) of particle motion within carefully designed geometries (Landman et al., 2010) and locally variable diffusivities. Along with numerical phantoms mimicking biological ones (Özcan et al., 2011), histopathological information is also being used for interpretation and validation (Budde and Frank, 2012; Budde and Annese, 2013). Additionally varying the time parameters δ and Δ will exploit the filtering effects caused by the displacement integral that will determine whether further information extraction is possible by expanding data acquisition with an appropriate set of parameter values.

On the other hand, on the discovery path's initiation by CFD–MRI signal formation, the re-establishment of Hermitian symmetry requires, in addition to the theoretical reasons presented herein, the analysis and quantification of Hermitian disruptive artifacts and systematic conditions in real data. Constructed by initially experimenting with elementary phantoms (e.g., water and mineral oil), this signal model expansion is necessary for the development of more elaborate systematic phase corrections, possibly by utilizing complex analysis theory. Specifically, accurate estimation of the pertinent Fourier transform of Pcfd from real data points in a clinical setup is targeted by the adaptation of CFD–MRI to fast sequences6, such as echo planar imaging (EPI) prone to Eddy current artifacts. The model will be expanded up to the point of reaching only minimal incremental improvements with new phase correction algorithms. Thereafter, relying on the residuals' content, which is free from displacement effects consequent to the application of system-wide uniform phase corrections7, more effective biomarker construction would be possible by the inclusion of extra information such as tissue susceptibility (Liu et al., 2013). Likewise, on a larger scale CFD–MRI's general aim is to improve outcomes of multimodal imaging, e.g., in prostate cancer strategies (Turkbey and Choyke, 2012).

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


This study was supported, in part, by the National Cancer Institute funded Washington University Small Animal Imaging Resource (U24-CA83060), the NIH/NINDS grant Biomarkers and Pathogenesis of MS (P01-NS059560), and the US Army grant, NeuroPerformance Imaging (W23RYX1089N603). Thanks to J. Quirk, T. Conturo, S.-K. Song, K. Uḡurbil, G. Sapiro, K. Wong, N. Tsekos, and S. K. Mun for valuable discussions, the reviewers for their valuable comments and C. D. Kroenke for the sample preparation.

The manuscript is dedicated to Ayla and Halil Özcan.


Assaf, Y., and Basser, P. (2005). Composite hindered and restricted model of diffusion (CHARMED) MR imaging of the human brain. Neuroimage 27, 48–58.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Baird, A. E., and Warach, S. (1998). Magnetic resonance imaging of acute stroke. J. Cereb. Blood Flow Metab. 18, 583–609.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Basser, P., Mattiello, J., and Lebihan, D. (1994). Estimation of the effective self-diffusion tensor from the NMR spin echo. J. Magn. Reson. Ser. B 103, 247–254.

Pubmed Abstract | Pubmed Full Text

Bernstein, M. A., King, K. F., and Zhou, X. J. (2004). Handbook of MRI pulse sequences. Amsterdam; Boston: Academic Press.

Bracewell, R. N. (2000). The Fourier Transform and its Applications, 3rd Edn. Boston, MA: McGraw Hill.

Brigham, E. O. (1988). The Fast Fourier Transform and its Applications. Englewood Cliffs, NJ: Prentice Hall.

Budde, M. D., and Annese, J. (2013). Quantifying anisotropy and fiber orientation in human brain histological sections. Front. Integr. Neurosci. 7:3. doi: 10.3389/fnint.2013.00003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Budde, M. D., and Frank, J. A. (2012). Examining brain microstructure using structure tensor analysis of histological sections. Neuroimage 63, 1–10.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Budde, M. D., Kim, J. H., Liang, H.-F., Schmidt, R. E., Russell, J. H., Cross, A. H., et al. (2007). Toward accurate diagnosis of white matter pathology using diffusion tensor imaging. Magn. Reson. Med. 57, 688–695.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Callaghan, P. T. (1991). Principles of Nuclear Magnetic Resonance Microscopy. Oxford; New York: Oxford University Press.

Callaghan, P. T., Eccles, C. D., and Xia, Y. (1988). NMR microscopy of dynamic displacements: k-space and q-space imaging. J. Phys. E Sci. Instrum. 21, 820.

Carr, H. Y., and Purcell, E. M. (1954). Effects of diffusion on free precession in nuclear magnetic resonance experiments. Phys. Rev. 94, 630–638.

Ciccarelli, O., Catani, M., Johansen-Berg, H., Clark, C., and Thompson, A. (2008). Diffusion-based tractography in neurological disorders: concepts, applications, and future developments. Lancet Neurol. 7, 715–727.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Conturo, T. E., Lori, N. F., Cull, T. S., Akbudak, E., Snyder, A. Z., Shimony, J. S., et al. (1999). Tracking neuronal fiber pathways in the living human brain. Proc. Natl. Acad. Sci. U.S.A. 96, 10422–10427.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Correia, S., Lee, S. Y., Voorn, T., Tate, D. F., Paul, R. H., Zhang, S., et al. (2008). Quantitative tractography metrics of white matter integrity in diffusion-tensor MRI. Neuroimage 42, 568–581.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cory, D. G., and Garroway, A. N. (1990). Measurement of translational displacement probabilities by NMR: an indicator of compartmentation. Magn. Reson. Med. 14, 435–444.

Pubmed Abstract | Pubmed Full Text

Doob, J. L. (1990). Stochastic Processes. New York, NY: Wiley-Interscience.

Frank, L. R. (2001). Anisotropy in high angular resolution diffusion-weighted MRI. Magn. Reson. Med. 45, 935–1141.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Haacke, E. M., Brown, R. W., Thompson, M. R., and Venkatesan, R. (1999). Magnetic Resonance Imaging: Physical Principles and Sequence Design. New York, NY: John Wiley and Sons.

Hahn, E. L. (1950). Spin echoes. Phys. Rev. 80, 580–594.

Hinshaw, W. S., and Lent, A. H. (1983). An introduction to NMR imaging: from the bloch equation to the imaging equation. Proc. IEEE 71, 338–350.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jensen, J., Helpern, J., Ramani, A., Lu, H., and Kaczynski, K. (2005). Diffusional kurtosis imaging: the quantification of non-gaussian water diffusion by means of magnetic resonance imaging. Magn. Reson. Med. 53, 1432–1440.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Kroenke, C. D., Bretthorst, G. L, Inder, T. E., and Neil, J. J. (2005). Diffusion MR imaging characteristics of the developing primate brain. Neuroimage 25, 1205–1213.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Landman, B. A., Bogovic, J. A., Wan, H., ElShahaby, F. E. Z., Bazin, P.-L., and Prince, J. L. (2012). Resolution of crossing fibers with constrained compressed sensing using diffusion tensor mri. Neuroimage 59, 2175–2186.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Landman, B. A., Farrell, J. A. D., Smith, S. A., Reich, D. S., Calabresi, P. A., and van Zijl, P. C. M. (2010). Complex geometric models of diffusion and relaxation in healthy and damaged white matter. NMR Biomed. 23, 152–162.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, C., Bammer, R., Acar, B., and Moseley, M. E. (2004). Characterizing non-gaussian diffusion by using generalized diffusion tensors. Magn. Reson. Med. 51, 924–937.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Liu, C., Murphy, N., and Li, W. (2013). Probing white-matter microstructure with higher-order diffusion tensors and susceptibility tensor mri. Front. Integr. Neurosci. 7:11. doi: 10.3389/fnint.2013.00011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mac Donald, C. L., Johnson, A. M., Cooper, D., Nelson, E. C., Werner, N. J., Shimony, J. S., et al. (2011). Detection of blast-related traumatic brain injury in U.S. military personnel. N. Eng. J. Med. 364, 2091–2100.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mattiello, J., Basser, P., and Lebihan, D. (1994). Analytical expressions for the b matrix in NMR diffusion imaging and spectroscopy. J. Magn. Reson. Ser. A 108, 131–141.

McCall, D. W., Douglass, D. C., and Anderson, E. W. (1963). Self-diffusion studies by means of nuclear magnetic resonance spin-echo techniques. Berichte der Bunsengesellschaft für Physikalische Chemie 67, 336–340.

McKinstry, R. C., Mathur, A., Miller, J. H., Özcan, A., Snyder, A. Z., Schefft, G. L., et al. (2002). Radial organization of developing preterm human cerebral cortex revealed by non-invasive water diffusion anisotropy mri. Cereb. Cortex 12, 1237–1243.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mielke, M., Kozauer, N., Chan, K., George, M., Toroney, J., Zerrate, M., et al. (2009). Regionally-specific diffusion tensor imaging in mild cognitive impairment and Alzheimer's disease. Neuroimage 46, 47–55.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mitra, P. P., and Halperin, B. I. (1995). Effects of finite gradient-pulse widths in pulsed-field-gradient diffusion measurements. J. Magn. Reson. Ser. A 113, 94–101.

Mori, S., and van Zijl, P. C. M. (2002). Fiber tracking: principles and strategies – a technical review. NMR Biomed. 15, 468–480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Moseley, M. E., Cohen, Y., Mintorovitch, J., Chileuitt, L., Shimizu, H., Kucharczyk, J., et al. (1990). Early detection of regional cerebral ischemia in cats: comparison of diffusion- and T2-weighted MRI and spectroscopy. Magn. Reson. Med. 14, 330–346.

Pubmed Abstract | Pubmed Full Text

Özarslan, E. (2009). Compartment shape anisotropy (CSA) revealed by double pulsed field gradient MR. J. Magn. Reson. 199, 56–67.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Koay, C. G., and Basser, P. J. (2009). Remarks on q-space MR propagator in partially restricted, axially-symmetric, and isotropic environments. Magn. Reson. Imaging 27, 834–844.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., and Mareci, T. H. (2003). Generalized diffusion tensor imaging and analytical relationships between diffusion tensor imaging and high angular resolution diffusion imaging. Magn. Reson. Med. 50, 955–965.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Nevo, U., and Basser, P. J. (2008). Anisotropy induced by macroscopic boundaries: surface-normal mapping using diffusion-weighted imaging. Biophys. J. 94, 2809–2818.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özarslan, E., Shepherd, T. M., Vemuri, B. C., Blackband, S. J., and Mareci, T. H. (2006). Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT). Neuroimage 31, 1086–1103.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2005). (Mathematical) necessary conditions for the selection of gradient vectors in DTI. J. Magn. Reson. 172, 238–241.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2009). “Theoretical and experimental analysis of imaging gradients in DTI,” in Proceedings of the 31st Annual International Conference of the IEEE EMB Society (Minneapolis, MN), 2703–2706.

Özcan, A. (2010a). Characterization of imaging gradients in diffusion tensor imaging. J. Magn. Reson. 207, 24–33.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2010b). “A new model for diffusion weighted MRI: complete fourier direct MRI,” in Proceedings of the 32st Annual International Conference of the IEEE EMB Society (Buenos Aires), 2710–2713.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2011). “Comparison of the complete fourier direct MRI with existing diffusion weighted MRI methods,” in Proceedings of the 2011 IEEE International Symposium on Biomedical Imaging (Chicago, IL), 931–934.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A. (2012). “Complete fourier direct magnetic resonance imaging (CFD–MRI) equations for diffusion MRI signal,” in Proceedings of the 15th International Conference on Medical Image Computing and Computer Assisted Intervention (Nice).

Özcan, A., Quirk, J., Wang, Y., Wang, Q., Sun, P., Spees, W., et al. (2011). “The validation of complete fourier direct MR method for diffusion MRI via biological and numerical phantoms,” in Engineering in Medicine and Biology Society, EMBC, 2011 Annual International Conference of the IEEE (Boston, MA), 3756–3759.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Özcan, A., Wong, K. H., Larson-Prior, L., Cho, Z.-H., and Mun, S. K. (2012). Background and mathematical analysis of diffusion mri methods. Int. J. Imaging Syst. Technol. 22, 44–52.

Papadakis, N. A., Xing, D., Huang, C. L.-H., Hall, L. D., and Carpenter, T. A. (1999). A comparative study of acquisition schemes for diffusion tensor imaging using MRI. J. Magn. Reson. 137, 67–82.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Rugh, W. J. (1995). Linear System Theory, 2nd Edn. New York, NY: Prentice Hall.

Seal, M. L., Yücel, M., Fornito, A., Wood, S. J., Harrison, B. J., Walterfang, M., et al. (2008). Abnormal white matter microstructure in schizophrenia: a voxelwise analysis of axial and radial diffusivity. Schizophr. Res. 101, 106–110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Shiryaev, A. N. (1995). Probability, 2nd Edn. New York, NY: Springer Verlag.

Song, S.-K., Qu, Z., Garabedian, E. M., Gordon, J. I., Milbrandt, J., and Ackerman, J. J. H. (2002). Improved magnetic resonance imaging detection of prostate cancer in a transgenic mouse model. Cancer Res. 62, 1555–1558.

Pubmed Abstract | Pubmed Full Text

Song, S.-K., Sun, S.-W., Ramsbottom, M. J., Chang, C., Russell, J., and Cross, A. H. (2002). Dysmyelination revealed through MRI as increased radial (but unchanged axial) diffusion of water. Neuroimage 17, 1429–1436.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stejskal, E. O., and Tanner, J. (1965). Spin diffusion measurements: spin echoes in the presence of a time-dependent field. J. Chem. Phys. 42, 288–292.

Torrey, H. C. (1956). Bloch equations with diffusion terms. Phys. Rev. 104, 563–565.

Tuch, D. S. (2004). Q-ball imaging. Magn. Reson. Med. 52, 1358–1372.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turkbey, B., Albert, P. S., Kurdziel, K., and Choyke, P. L. (2009). Imaging localized prostate cancer: current approaches and new developments. Am. J. Roentgenol. 192, 1471–1480.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turkbey, B., and Choyke, P. L. (2012). Multiparametric mri and prostate cancer diagnosis and risk stratification. Curr. Opin. Urol. 22, 310–315.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Voineskos, A. N., Lobaugh, N. J., Bouix, S., Rajji, T. K., Miranda, D., Kennedy, J. L., et al. (2010). Diffusion tensor tractography findings in schizophrenia across the adult lifespan. Brain 133, 1494–1504.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wedeen, V. J., Hagmann, P., Tseng, W.-Y. I., Reese, T. G., and Weisskoff, R. M. (2005). Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging. Magn. Reson. Med. 54, 1377–1385.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Xu, J., Humphrey, P. A., Kibel, A. S., Snyder, A. Z., Narra, V. R., Ackerman, J. J., et al. (2009). Magnetic resonance diffusion characteristics of histologically defined prostate cancer in humans. Magn. Reson. Med. 61, 842–850.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text


wi: ith magnetic moment's displacement from its initial position xi(t0).

Wd: The difference of magnetic moment displacement integrals during the periods diffusion (d) gradients are turned on. Its distribution is the main quantity of interest providing information about microstructure. Joined with the other displacement integrals during acquisition and rewind times, Wacq, Wrw, respectively, it forms the total displacement integral: W=(Wd,Wrw,W1acq)7.

Pcfd(x, Wd): CFD joint distribution function of the number of magnetic moments with initial position x possessing displacement integral values of Wd. It is obtained by marginalizing Ptotalcfd(x, W).

Scfd(kmr, kD, krw): DW–MRI signal that comes out of the scanner which is equal to the Fourier transform of Ptotalcfd. Theoretically it must be Hermitian symmetric since Ptotalcfd is real valued.

(kmr, kD): kcfd-space variables corresponding to imaging gradients (mr) and diffusion (D) gradients, respectively.

krw: Rewinding (rw) frequency variable affine dependent on kmr. It cannot be sampled at 0 because of its dependence. For Ptotalcfd marginalization, phase corrections are applied to estimate Scfd at krw = 0.

Ikcomplexmr(x, kD): Complex valued images obtained by taking the Fourier transform of Scfd with respect to imaging frequency kmr.

Ikmr(x, kD): Image domain phase corrected Ikcomplexmr(x, kD).


A. Second order Statistics for Displacement Integrals of Self Diffusion

In this section, the covariance of displacement integrals is derived for the special case of particles executing Brownian motion. The following preliminaries are necessary to calculate the second order statistics for the displacement integrals of Equation (9). It is assumed for ease of notation below that tmtm + 1.

The Wiener process (Shiryaev, 1995), w, which describes the diffusion in an isotropic homogenous medium, is sample path continuous, has w(0) = 0 and independent increments (i.e., Markovian property) with a normal distribution, meaning w(t)w(s)~(0,(ts)D) where t > s ≥ 0 and D > 0. Using these properties, the covariance of the displacement is:

E[w(t) wT(s)]=E[(w(t)w(s)+w(s)) wT(s)]                     =E[w(s)wT(s)]=min(t, s) D.(29)

It is straightforward to show that the mean for the displacement integral processes is 0. The covariance in the case of a single interval is calculated as


and for a non-overlapping interval:

t3t4t1t2E[wi(τ) wiT(s)] dτ ds=t3t4t1t2min(τ, s)D dτ ds=12(t4t3)(t22t12)D.(31)

Finally, using the formulas above, the time points of Figure 1, [t1 t2 t3 t4] = [td1 td2 td3 td4], and these substitutions,


the following is obtained for Wdi:

E[Wid(Wid)T]=[((t2t1)2(2t1+t2)+(t4t3)2                           (2t3+t4))/3(t4t3)(t22t12)]D                      =δ2(Δδ/3)D.(32)

The scalar factor, δ2 (Δ − δ/3), defined as bt-value in Özcan (2009, 2010a), is a function of time directly related to the measured quantity, the displacement integrals. It is completely detached from the measurement parameters such as diffusion gradient strength, in contrast to the b-value used in the literature, b = γ2GD22 bt. The dependence of the derivations on Markovian property restricts the validity of the calculations to isotropic samples where the diffusion is characterized by the diffusion coefficient, D, in the distribution of magnetic moment displacement w. For an isotropic homogenous medium, D is estimated from the displacement integral's (Wd) covariance bt D = δ2 (Δ − δ/3)D using DW–MRI acquisitions.

The rigorous treatment of the theory of the stochastic processes in Doob, (1990, p. 62) demonstrates that the sample paths of stochastic processes are Lebesgue integrable under certain conditions and these integrals are well defined random variables. Although a rigorous mathematical analysis, which proves that these conditions are met, is out of the scope of this manuscript, it is important to note that Equation (23) does not prove that the displacement integrals are Gaussian random variables. In this case, the central limit theorem does not apply because the displacement integrals of Equations (8, 9), are not the sums of independent variables since in the integral approximating sum

t0t1wi(τ) dτlimNk=1Nwi(t0+k dτ) dτ      (N, dτ0, N dτ=t1t0),(33)

wi(τ), rather than independent increments, wi(t) − wi(s), is present. A variant of central limit theorem for sums of dependent variables satisfying certain conditions (Shiryaev, 1995, p. 541) might provide the theoretical validation of this highly theoretical open problem.

B. Derivation of the Phase Equation

B.1. Bloch equations and their solutions

The time evolution of the magnetization, m¯i(t)3, that generates the MR signal is modeled by the Bloch equations:

dm¯i(t)dt=γ m¯i(t)×B(xi, t)[1T2i0001T2i0001T1i]m¯i(t)              +1T1i[00m3i](34)

with T1i and T2i denoting spin–lattice and spin–spin relaxation constants, respectively. In the MR scanner the magnetic field as a function of time and position is given by

B(x, t)=[00B0]+B1m(t)+[00G(x, t) · x(t)](35)

where B0 is the static (strong) magnetic field, Bm1(t) is the radio frequency (RF) pulse modulated at the precession frequency of B0 and G(x,t)3 describes the magnetic field gradients. The expression of the magnetic field simplifies further in the rotating frame:

B(x, t)=B1(t)+[00G(x, t)  x(t)].(36)

When the magnetic field gradients are turned on without the RF pulse the system matrix becomes

m¯i×[00ω(xi(t))]=[0ω(xi(t))0ω(xi(t))00000]m¯i=A(t) m¯i.(37)

with ω(xi(t)) = G(xi, t) · xi(t). As the time-dependent system matrix A(t) commutes with its integral ∫ A(τ)dτ, when the relaxation terms are negligible the solution of Equation (34) is obtained [see Rugh, (1995, pg. 59)] as:

m¯i(t)=exp(γt0tA(τ)dτ) m¯i(t0)(38)

by calculating the matrix exponential using Ω(xi(t), t0) = ∫tt0 ω(xi(τ)) dτ = ∫tt0 G(xi,τ)· xi(τ) :

exp(γt0tA(τ)dτ)=[cos(γΩ(xi(t), t0))sin(γΩ(xi(t), t0))0sin(γΩ(xi(t), t0))cos(γΩ(xi(t), t0))0001].

Equation (38) defines a system with a rotating transverse magnetization component that can be written in a complex form yes. The corresponding first order differential equation is


with the solution


For the portion of the pulse sequence involving the RF pulses, in general (Bernstein et al., 2004) the magnetic moments are assumed to be immobile [i.e., wi(t) ≡ 0 in Equation (4)] thereby defining a constant system matrix that has a standard matrix exponential that solves the differential equation. Herein, however, the effect of the displacements specifically during the π-pulse results in a time varying system matrix. Without loss of generality, the π-pulse is applied in the direction8. [1, 0, 0]T as in Hinshaw and Lent (1983), along with the slice select gradient, Gπ = [0, 0, gπ3], active only during the RF pulse. The RF pulse modulated at the same frequency corresponding to the center of the slice, B0 + gπ3 xc3, becomes a constant vector, B1(t) = [B1a, 0, 0]T, in Equation (36). The system of equations in this rotating frame is

mi×[B1a0ωπ(xi(t))]=[0ωπ(xi(t))0ωπ(xi(t))0B1a0B1a0]mi=A(t) mi(41)

with9 ωπ(xi(t)) = Gπ · (xi(t0) + wi(t) − xc) = gπ3 (x3i(t0) + w3i(t) − xc3). The procedure would be to solve Equation (41) and then express the solution in the initial rotating frame to obtain a phase component, Φπi, that would be a function of the displacement integral and the initial position from the slice center.

However, the new system matrix, A(t), does not commute with its integral preventing the calculation of a matrix exponential formulating a convenient analytical solution as in Equation (40). With the rigorous analysis describing the exact manifestation of the RF pulse left to future studies, herein due to the large magnitude of the RF field in comparison to the magnetic field gradient, its effect will be approximated as a sign reversal of the phase. In this work, the phase change induced by the magnetic field gradient is considered as a factor that is automatically corrected by the phase correction algorithm of section 2.2.

On the other hand, the effect of the displacement induced by the RF slice select gradient (SS–EX in Figure 1) during the excitation period is systematic. As both the RF pulse and the magnetic field gradient remain fixed, statistically speaking, the phase of the total transverse magnetization is affected by the displacements in the same amount at each acquisition prior to the tipping. The position and displacement encoding occurs thereafter during the remainder of the pulse sequence (see also footnote 1 for the phase correction along the slice select direction within the slab) leaving the π/2-pulse out of the derivation of the phase equations below.

B2. Phase equation

For clarity of exposition, the calculations are carried out with the assumption of ideal gradient amplifiers creating rectangular shaped gradient pulses in the time domain and perfect linearity in physical space. In short, during the time interval when the amplifiers are powered on G(x,t)=G3.

The derivation of Equations (5– 7) provided in Özcan (2012) is briefly presented in this section for reference purposes. Using the definitions of the variables in Figure 1, the evolution of the phase is described as:

1. First, imaging gradients for read out rewinding, Gro3, phase encoding, Gpe3 and slice select Gss3 on the time interval [t0, trw] are turned on. After these gradients are applied the phase and the magnetization become

Ωrw=t0trw(Gro(τ)+Gpe(τ)+Gss(τ)) ·xi(τ) dτ      =(Gro+Gpe+Gss) · (Δtrwxi(t0)+t0trwwi(τ) dτ)      =Δtrw(Gro+Gpe+Gss) · xi(t0)         +(Gro+Gpe+Gss) · t0trwwi(τ) dτ.(42)


2. The RF π-pulse between the diffusion gradient pulses, GD3, and the accompanying slice select gradient provide theoretical sign reversal of the phase. The slice select gradient's encoding of magnetic moment motion into the signal is expressed by Φπi (see Appendix B1). Since xi(tdk) = xi(t0) + wi(tdk), k = 1, …, 4; at t = td4

ΩD=Φπi+td3td4GD · xi(τ) dτtd1td2GD · xi(τ) dτ(44)     =Φπi+GD · ((td4td3)xi(t0)+td3td4wi(τ) dτ)         GD · ((td2td1)xi(t0)+td1td2wi(τ) dτ)     =Φπi+GD · (td3td4wi(τ) dτtd1td2wi(τ) dτ)         +((td4td3)(td2td1))GD · xi(t0).(45)

If the diffusion gradient times are equal, i.e., td4td3 = td2td1 = δ, then the last term in Equation (6) is equal to zero, erasing the influence of initial position from motion encoding part of the signal. This is the insight gained by using the formulation with displacement integrals, ∫ wi(τ) dτ, rather than the center of mass (COM) of random walk, ∫ x(τ) dτ, introduced in Mitra and Halperin (1995).

The sign reversal affects also previously accumulated phase:


3. The last part of the sequence is where the data are collected:

Ωro(t)=tacqtGro· (xi(t0)+wi(τ)) dτ=(ttacq) Gro               · xi(t0)+tacqtGro · wi(τ) dτ,(47)

leading to


and Equations (5– 7).

C. Total number of Particles from CFD-MRI

The fundamental Fourier relationship of Equation (16), Scfd(kmr, kD, krw) = ℱ{Ptotalcfd}(kmr, kD, krw) establishes the relationship between the standard MR image space and the higher dimensional CFD space. By Equation (16)


and by Equation (18), Ikcomplexmr(x, kD, krw) = ℱ−1kmr{Scfd(kmr, kD,krw)}, the image obtained without the diffusion encoding magnetic field gradients is



=Pcfdtotal(x,W)δ(γ2π(xx))dxdW=2πγPcfdtotal(x,W)dW~total number of particles at x.(52)(53)(54)


  1. ^The refocusing lobe of the slice select gradient after the RF pulse denoted by SS in Figure 1 adjusts ϕslice to be constant throughout the slice in the slice select direction.
  2. ^The unit for kmr is magnetic field×timedistance but for kD and krw it is magnetic fielddistance. They are both in accordance with the units of the position and the displacement integrals in Equation (13). After kcfd is multiplied by the gyromagnetic ratio γ, with unit (magnetic field × time)−1, the product becomes unitless in Equation (14).
  3. ^Another approach in the literature is to present the value of the function over a sphere (Tuch, 2004; Wedeen et al., 2005). However, uniqueness is lost in this representation as demonstrated by these functions: f(x, y, z) = x2 + y2 + z2, h(x, y, z) = (f(x, y, z))2, which have the same value on the unit sphere but different isosurfaces.
  4. ^Also named nuclear spin self-correlation function (Callaghan et al., 1988) and translation probability (Cory and Garroway, 1990).
  5. ^One exception is the basic case of isotropic diffusion analyzed in Appendix A where the diffusion coefficient that describes the Gaussian propagator can be recovered using the adjustment factor for the displacement integral covariance namely bt in Equation (24) and the b-value in DTI formulation.
  6. ^Acquisition time is shortened with reduced sampling schemes aimed at specific goals, e.g., compressed sensing for tractography (Landman et al., 2012) at the expense of some information loss.
  7. ^Instead of pixel by pixel corrections that would completely eliminate the imaginary part in Figure 3.
  8. ^The choice of direction does not matter (Haacke et al., 1999, p. 387).
  9. ^The absence of a rewinding slice select gradient, in contrast to π/2-pulse in Figure 1, requires the addition of the slice center, xc, in the formulation as there will be a phase shift along the slice select direction in the slab [see Equation (11) and footnote 1].

Keywords: magnetic resonance imaging, diffusion weighted imaging, fourier transform

Citation: Özcan A (2013) Complete fourier direct magnetic resonance imaging (CFD-MRI) for diffusion MRI. Front. Integr. Neurosci. 7:18. doi: 10.3389/fnint.2013.00018

Received: 01 November 2012; Paper pending published: 29 January 2013;
Accepted: 08 March 2013; Published online: 02 April 2013.

Edited by:

Elizabeth B. Torres, Rutgers University, USA

Reviewed by:

Bennett A. Landman, Vanderbilt University, USA
Evren Ozarslan, Brigham and Women's Hospital, USA

Copyright © 2013 Özcan. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in other forums, provided the original authors and source are credited and subject to any copyright notices concerning any third-party graphics etc.

*Correspondence: Alpay Özcan, Health Research, Arlington Innovation Center, Virginia Polytechnic Institute and State University, 900 N. Glebe Road, Arlington, VA 22203, USA. e-mail: