Method for performing diffusion weighted magnetic resonance measurements
11525880 · 2022-12-13
Assignee
Inventors
Cpc classification
G01R33/5611
PHYSICS
G01R33/50
PHYSICS
A61B5/055
HUMAN NECESSITIES
International classification
G01R33/565
PHYSICS
G01R33/561
PHYSICS
G01R33/50
PHYSICS
Abstract
Disclosed is a method for generating a time-dependent magnetic field gradient in diffusion weighted magnetic resonance imaging G(t)=[G.sub.x(t)G.sub.y(t)G.sub.z(t)].sup.T, which is asymmetric in time with respect to a refocusing pulse, by meeting one or more of the requirements: A=∫.sub.0.sup.TEh(t)G(t)G(t).sup.Tdt is zero, where TE is an echo time and h(t) is a function of time which is positive during an interval prior to the refocusing pulse and negative during a time interval after the refocusing pulse); minimize A or m=(Tr[AA]).sup.1/2 where A=∫.sub.P1G(t)G(t).sup.Tdt−∫.sub.P2G(t)G(t).sup.Tdt where P1 and P2 represent time intervals prior to and subsequent to the refocusing pulse; m is smaller than a threshold value. an attenuation factor
due to T2* relaxation is one. Signal attenuation due to concomitant field gradients, regardless of the shape or orientation of the diffusion encoding b-tensor and the location of signal is hereby minimized.
Claims
1. A method for diffusion weighted magnetic resonance imaging, comprising: generating by a gradient coil of a magnetic resonance imaging scanner a time-dependent magnetic field gradient G(t)=[G.sub.x(t)G.sub.y(t)G.sub.z(t)].sup.T, wherein the gradient G is asymmetric in time with respect to a refocusing pulse and wherein the gradient G is such that A=∫.sub.0.sup.TEh(t)G(t)G(t).sup.Tdt is zerois zero, where TE is an echo time and h(t) is a function of time which is positive during an interval prior to the refocusing pulse and negative during a time interval after the refocusing pulse.
2. A method according to claim 1, wherein h(t) is a function of time which is positive for 0<t<TE/2 and negative for TE/2<t<TE.
3. A method according to claim 1, wherein the magnetic field gradient G forms part of a spin-echo sequence.
4. A method according to claim 1, wherein the gradient G is configured such that a diffusion encoding tensor representation b has at least two non-zero eigenvalues, where: b=∫.sub.0.sup.TEq(t)q(t).sup.Tdt, and q(t)=γ∫.sub.0.sup.tG(t′)dt′ is a dephasing vector.
5. A method according to claim 1, wherein the gradient G is adapted to cause isotropic diffusion encoding.
6. A method for diffusion weighted magnetic resonance imaging, comprising: generating by a gradient coil of a magnetic resonance imaging scanner a time-dependent magnetic field gradient G=[G.sub.x(t) G.sub.y(t) G.sub.z(t)].sup.T, wherein the gradient G is asymmetric in time with respect to a refocusing pulse and wherein the gradient G is such that an attenuation factor
7. A method according to claim 6, wherein the gradient G further is such that an attenuation factor AF.sub.s=sinc(k.sub.s.Math.ST) is one, where k.sub.s−k.Math.n.sub.s, where n.sub.s, is a normal direction of a slice plane and where ST is a slice thickness.
8. A method according to claim 6, wherein the magnetic field gradient G forms part of a spin-echo sequence.
9. A method according to claim 6, wherein the gradient G is configured such that a diffusion encoding tensor representation b has at least two non-zero eigenvalues, where: b=∫.sub.0.sup.TEq(t)q(t).sup.Tdt, and q(t)=γ∫.sub.0.sup.tG(t′)dt′ is a dephasing vector.
10. A method according to claim 6, wherein the gradient G is adapted to cause isotropic diffusion encoding.
11. A method for diffusion weighted magnetic resonance imaging, comprising: generating by a gradient coil of a magnetic resonance imaging scanner a time-dependent magnetic field gradient G(t)=[G.sub.x(t) G.sub.y(t) G.sub.z(t)].sup.T, wherein the gradient G is asymmetric in time with respect to a refocusing pulse and wherein the gradient G is such that m=(Tr[AA]).sup.1/2 is smaller than a threshold value, where A=∫.sub.P1G(t)G(t).sup.Tdt−∫.sub.P2G(t)G(t).sup.Tdt where P1 and P2 represent time intervals prior to and subsequent to the refocusing pulse.
12. A method according to claim 11, wherein the gradient G is such that m is zero.
13. A method according claim 11, wherein the gradient G is configured such that a diffusion encoding tensor representation b has at least two non-zero eigenvalues, where: b=∫.sub.0.sup.TEq(t)q(t).sup.Tdt, and q(t)=γ∫.sub.0.sup.tG(t′)dt′ is a dephasing vector.
14. A method according to claim 11, wherein the gradient G is adapted to cause isotropic diffusion encoding.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
(1) The above, as well as additional objects, features and advantages of the present inventive concept, will be better understood through the following illustrative and non-limiting detailed description of preferred embodiments of the present inventive concept, with reference to the appended drawings. In the drawings like reference numerals will be used for like elements unless stated otherwise.
(2)
(3)
(4)
DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
(5) To facilitate understanding of the present inventive concept, a discussion of some theoretical concepts will now be provided.
(6) Theory
(7) The Maxwell equations show that linear field gradients produced by the gradient coils of an MRI scanner are accompanied by additional fields (concomitant fields) that are spatially dependent (Norris and Hutchison, 1990; Bernstein et al, 1998; Zhou et al., 1998; Meier et al, 2008). Their contribution to the field can be approximated by an expansion of the magnetic field in the Cartesian coordinates x, y, and z with origin at the isocenter
(8)
Where B.sub.0 is the main magnetic field, G.sub.x x, G.sub.y y, and G.sub.z z are the time-dependent field gradients created by the gradient coil and the other terms are the so-called Maxwell terms.
(9) The deviation from the ‘intended’ gradients may cause undesired signal attenuation, which is a problem particularly for diffusion MRI (Baron et al., 2012). The undesired attenuation is caused by two effects: a through-plane phase dispersion, and a delay of the echo formation. Both effects can be predicted from the residual gradient moment k, which for a spin-echo sequence is given by
(10)
where t.sub.0 and t.sub.1 are the times from the excitation to the 180 pulse for k.sub.pre, and from the 180 pulse to the echo time for k.sub.post. In a spin-echo sequence the moment can be nulled if g(t) is identical for the time periods before and after the 180 pulse, because then k.sub.pre=k.sub.port. In other cases, additional analysis is required.
(11) We can compute the effective gradient g(t) from concomitant gradients from the spatial derivatives of Eq. 1, according to
g(t)=G(t)+Mr Eq. 4
where G(t)=[G.sub.x(t) G.sub.y(t) G.sub.z(t)].sup.T is the desired time-varying gradient produced by the gradient coil, r=[x y z].sup.T is the location in space, and
(12)
Through-plane phase dispersion is caused by a non-zero moment along the slice encoding direction, given by
k.sub.s=k.Math.n.sub.s Eq. 6
where n.sub.s is a normal vector of the slice plane. The resulting attenuation is defined by the slice selection profile. Here, we assume that it is a rectangular function, and thus the attenuation factor (AF) is given by a sinc function (Baron et al., 2012)
AF.sub.s=sinc(k.sub.s.Math.ST), Eq. 7
where ST is the slice thickness.
(13) As realized by the inventors, there is a second potential effect that causes a signal attenuation, not described before. Here, a delay in the echo time is given from the residual gradient moment along the phase encoding direction (assuming an EPI readout is used), given by
k.sub.p=k.Math.n.sub.p Eq. 8
This gradient moment is translated to an echo-displacement time t (which alternatively may be labelled t′) according to
(14)
where Δk is the distance between two acquired k-space lines, given by Δk=n/FOV, where n is the parallel imaging factor and FOV is the field of view, and Δt is the time between the acquisitions of the central echo of two consecutive k-space lines. The echo-displacement can yield additional attenuation due to T2* relaxation according to
(15)
(16) How to Minimize Effects of Concomitant Gradients
(17) As realized by the inventors, both k.sub.s and k.sub.p must be minimized in order to get rid of unwanted signal attenuation. A sufficient (but not necessary) condition for this to be achieved is the fulfilment of the following condition:
A=∫.sub.0.sup.TEh(t)G(t)G(t).sup.Tdt=0 Eq. 11
where TE is the echo time, h(t) is a function that is positive when 0<t<TE/2 and negative when TE/2<t<TE to account for effects of the 180 pulse in a spin echo experiment, and G(t)G(t).sup.T is the outer product of G(t) with itself, i.e., a matrix of size 3 by 3. More generally, h(t) may be a function which is positive during an interval prior to a refocusing pulse and negative during a time interval after a refocusing pulse. The diagonal elements of A will contain the integrals of G.sub.x.sup.2, G.sub.y.sup.2 and G.sub.z.sup.2 terms after accounting for the 180 pulse, and the off-diagonal elements will contain the G G terms. If all of these terms are nulled, then the average of M from excitation to echo formation will be nulled, and thus k will be nulled so that no unwanted signal attenuation should appear. The condition that A=0 can be realised as an additional condition to the waveform optimization procedure described in Sjölund et al (2015).
(18) An important consequence of the problem formulation in Eq. 11 is that any gradient waveform that fulfils the condition can be rotated and rescaled freely, because rotating and rescaling by the waveform by R yields A′=R A R.sup.T=0 when A=0. It will thus guarantee a minimization of concomitant field artefacts regardless of the tensor shape.
(19) Hereafter follows an example implementation of a wave form optimizer taking advantage of the above derivations.
(20) Methods: Implementation of the concomitant-gradient constraint in a waveform optimizer
(21) A waveform optimizer may be implemented using the NOW toolbox, which may be downloaded from https://gifthub.com. (with commit ID 6c69462). The following modifications may be made to produce the desired gradient waveforms:
(22) 1. Optimization Constraint
(23) This constraint is used to minimize the components of the outer product of the gradient waveform with itself. Additionally, effects of the refocusing pulse is taken into account by changing the sign of the integral of the second part. The two intervals, i.e., the gradient waveform before and after the refocusing are denoted by the indices pre and post, and the axes along which the gradient are applied are denoted x, y and z. To minimize the magnitude of each component, the sum of squares is considered in the optimization, down to a tolerance level denoted tolMaxwell. The new part of the optimizer therefore takes on the form of the constraint (const), according to
const=(sum(g(pre,x).*g(pre,x))−sum(g(post,x).*g(post,x))).{circumflex over ( )}2+(sum(g(pre,x).*g(pre,y))−sum(g(post,x).*g(post,y))).{circumflex over ( )}2+(sum(g(pre,x).*g(pre,z))−sum(g(post,x).*g(post,z))).{circumflex over ( )}2+(sum(g(pre,y).*g(pre,y))−sum(g(post,y).*g(post,y))).{circumflex over ( )}2+(sum(g(pre,y).*g(pre,z))−sum(g(post,y).*g(post,z))).{circumflex over ( )}2+(sum(g(pre,z).*g(pre,z))−sum(g(post,z).*g(post,z))). {circumflex over ( )}2−tolMaxwell
where g is the time dependent gradient waveform along three orthogonal axes.
(24) An implementation of the modification (in Matlab code) createConstraintGradientFunction.m of the NOW framework is shown below.
(25) createConstraintGradientFunction.m
(26) TABLE-US-00001 function createConstraintGradientFunction(N,useMaxNorm) q = sym(‘q’,[3*N,1]); targetTensor = sym(‘targetTensor’,[3,3]); syms s B gMax tolIsotropy integralConstraint c1 c2 tolMaxwell real integrationMatrix = eye(N); integrationMatrix(1,1) = 0.5; integrationMatrix(N,N) = 0.5; B = (reshape(q,[N 3])‘*integrationMatrix*reshape(q,[N 3]))*1/(N−1); c1 = norm(B−s*targetTensor,‘fro’){circumflex over ( )}2−(s*tolIsotropy){circumflex over ( )}2; % Center difference, shifted forward by half a step firstDerivativeMatrix = −diag(ones(N,1))+diag(ones(N−1,1),1); firstDerivativeMatrix = firstDerivativeMatrix(1:end−1,:); firstDerivativeMatrix = sparse(firstDerivativeMatrix); g = firstDerivativeMatrix*reshape(q,[N 3]); %No need to include the zeros at start and end %Power constraint: constrains the integral of g(t){circumflex over ( )}2. Since the gradients %already represent the average of each time interval the trapezoid rule %should not be used c3 = g(:,1)‘*g(:,1)−integralConstraint; c4 = g(:,2)‘*g(:,2)−integralConstraint; c5 = g(:,3)‘*g(:,3)−integralConstraint; hlf = ceil(N/2); pre = 1:hlf; post = (hlf+1):(N−1); x=1; y=2; z=3; % Maxwell compensated constraint c6 = (sum(g(pre,x).*g(pre,x)) − sum(g(post,x).*g(post,x))).{circumflex over ( )}2 + ... (sum(g(pre,x).*g(pre,y)) − sum(g(post,x).*g(post,y))),{circumflex over ( )}2 + ... (sum(g(pre,x).*g(pre,z)) − sum(g(post,x).*g(post,z))),{circumflex over ( )}2 + ... (sum(g(pre,y).*g(pre,y)) − sum(g(post,y).*g(post,y))).{circumflex over ( )}2 + ... (sum(g(pre,y).*g(pre,z)) − sum(g(post,y).*g(post,z))).{circumflex over ( )}2 + ... (sum(g(pre,z).*g(pre,z)) − sum(g(post,z).*g(post,z))).{circumflex over ( )}2 − ... tolMaxwell; if useMaxNorm == false c2 = (sum(g.{circumflex over ( )}2,2)−gMax{circumflex over ( )}2)‘;%Nonlinear inequality constraint <= 0 c = [c1 c2 c3 c4 c5 c6]; gradc = jacobian(c,[q;s]).‘; % transpose to put in correct form else c = [c1 c3 c4 c5 c6]; gradc = jacobian(c,[q;s]).‘; % transpose to put in correct form end if useMaxNorm fileName = [‘private/nonlcon’ num2str(N) ‘pointsMaxNorm’]; else fileName = [‘private/nonlcon’ num2str(N) ‘points2Norm’]; end matlabFunction(c,[ ],gradc,[ ],‘file’,fileName,‘vars’,{[q;s],tolIsotropy,gMax, integralConstraint,targetTensor, tolMaxwell}); %The [ ] are for the inequality constraints
(27) 2. Tolerance Inclusion
(28) To engage the constraint, the minimization function needs to be provided a tolerance. The call to fmincon is therefore expanded with problem.tolMaxwell, according to
(29) [x,fval,exitflag,output,lambda,grad]=fmincon(@(x) objFun(x), x0, A,b,Aeq,beq,[ ],[ ],@(x) feval(nonlconFileName, x, problem.tolIsotropy, problem.gMaxConstraint, problem.integralConstraint, problem.targetTensor, problem.tolMaxwell), options);
(30) Results
(31) The additional integral constraint was implemented in the waveform optimizer. Four waveforms were generated, including the use of Euclidean or Max norms for the gradient amplitude, or with and without the concomitant-gradient constraint. Waveforms are illustrated in
(32) TABLE-US-00002 For FIG. 1: sum(G.sup.T * G)= 7.5e+03 3.7e+04 3.5e+04 3.7e+04 8.3e+04 4.3e+04 3.5e+04 4.3e+04 −4.4e+04 For FIG. 2: sum(G.sup.T * G)= 3.5e+04 1.5e+05 1.7e+05 1.5e+05 3.7e+04 4.3e+04 1.7e+05 4.3e+04 4.3e+04 For FIG. 3: sum(G.sup.T * G)= 0.24 −0.096 −0.58 −0.096 0.12 0.34 −0.58 0.34 0.53 For FIG. 4: sum(G.sup.T * G)= 0.078 0.41 0.28 0.41 0.41 −0.022 0.28 −0.022 0.053
(33) The results demonstrate that it is possible to obtain efficient gradient waveform for isotropic diffusion encoding without causing undesired accumulation of additional gradients due to concomitant fields.
(34) Alternative Waveform Optimization Procedure
(35) Asymmetric magnetic field gradient waveforms G(t)=[G.sub.x(t)G.sub.y(t)G.sub.z(t)].sup.T enabling the effects of concomitant fields to be reduced to a negligible level may also be designed using that the ‘Maxwell index’ m−(Tr[AA]).sup.1/2 should be smaller than a predetermined threshold, i.e. within a tolerance, as an additional condition to the waveform optimization procedure. Such a constraint may be engaged in the NOW framework by providing a tolerance, similar to what is outlined above for A=0.
DESCRIPTION OF EMBODIMENTS
(36) In accordance with the present inventive concept there is provided a method for diffusion magnetic resonance imaging (dMRI). The method may be performed using an MRI scanner 1, schematically illustrated in
(37) The encoding block may further comprise a radio frequency (RF) sequence configured to provide relaxation encoding (i.e. attenuation due to longitudinal and/or transverse relaxation). The RF sequence may be generated by a RF transmission system 14 of the scanner 1. The RF sequence and the magnetic gradient G may be generated and applied to the sample or measurement object “S”, which for example may be a biological sample including water, such as brain tissue or biopsy samples of (suspensions) of any organs cell. More generally, the sample includes a nuclear spin system whose properties may be measured by nuclear magnetic resonance techniques. As indicated in
(38) The MRI scanner may comprise a controller 24 for controlling the operation of the scanner 1, in particular the magnet 10, the gradient coil 12, the RF transmission and receiving systems 14, 16 and the signal acquisition, etc. The controller 24 may be implemented on one or more processors of the scanner 1 wherein control data for generating the magnetic gradient and RF sequence may be implemented using software instructions which may be stored on a computer readable media (e.g. on a non-transitory computer readable storage medium) and be executed by the one or more processors of the device. The software instructions may for example be stored in a program/control section of a memory of the device, to which the one or more processors of the device has access. It is however also possible that, instead of using software instructions, the controller method may be implemented in the form of dedicated circuitry of the device/computer such as in one or more integrated circuits, in one or more application-specific integrated circuits (ASICs) or field-programmable gate arrays (FPGAs), to name a few examples.
(39) According to the method, the time-dependent magnetic field gradient G(t) generated by the gradient coil 12 may be configured to present asymmetry about a 180° refocusing pulse of the RF sequence. The gradient G may further be configured/designed in such a manner that A=∫.sub.0.sup.TEh(t)G(t)G(t).sup.Tdt is zero, where TE is an echo time and h(t) is a function of time which is positive during an interval prior to the refocusing pulse and negative during a time interval after the refocusing pulse.
(40) The gradient G may also be defined such that an attenuation factor
(41)
due to T2* relaxation is one, where t is an echo displacement time according to
(42)
where k is a residual gradient moment and n.sub.p defines the phase encoding direction. Δk is a distance between two acquired k-space lines given by Δk=n/FOV, where n is a parallel imaging factor and FOV is a field of view. Δt is a time between acquisitions of a central echo of two consecutive k-space lines.
(43) The gradient G may further be configured such that an attenuation factor AF.sub.s=sinc(k.sub.s.Math.ST) is one, where k.sub.s=k.Math.n.sub.s where n.sub.s is a normal direction of a slice plane and where ST is a slice thickness.
(44) The gradient G may also be defined such that m−(Tr[AA]).sup.1/2 is zero or at least smaller than a threshold value. The threshold value may be a predetermined threshold value set such that the effects of concomitant fields may be reduced to a negligible level.
(45) A gradient waveform G(t)=[G.sub.x(t)G.sub.y(t)G.sub.z(t)].sup.T meeting one or more of these conditions may for instance be determined using a waveform optimizer, such as the modified waveform optimizer described above. More generally, an asymmetric gradient waveform G(t) may be designed by performing a waveform optimization procedure configured to determine waveform components [G.sub.x(t)G.sub.y(t)G.sub.z(t)].sup.T of the magnetic field gradient G. The optimization procedure may be performed under one of the following constraints:
minimize A=∫.sub.0.sup.TEh(t)G(t)G(t).sup.Tdt, where TE is an echo time and h(t) is a function of time which is positive during an interval prior to a time instant of a refocusing pulse and negative during a time interval after the time instant of the refocusing pulse, or
minimize m=(Tr[AA]).sup.1/2, where A=∫.sub.p1G(t)G(t).sup.Tdt−∫.sub.p2G(t)G(t).sup.Tdt where P1 and P2 represent time intervals prior to and subsequent to a time instant of a refocusing pulse. The waveform optimizer may be configured with a minimization tolerance. In other words, the waveform optimizer may be configured to minimize A or m to be within a predetermined tolerance (i.e. within a tolerance from zero), which may be determined in accordance with equipment sensitivity and measurement requirements. The waveform optimizer may be implemented as a computer program product comprising a non-transitory computer-readable storage medium with instructions adapted to carry out the method for determining the wave form components of the gradient G when executed by a processor (of the computer). The term “non-transitory” as used herein may cover any tangible storage medium such as non-volatile memory, such as ROM, (intended to provide persistent storage of data) as well as volatile memory, such as RAM, (not intended to store data permanently). “Non-transitory” should hence be considered a limitation of the storage medium per se rather than a limitation on data storage persistency.
(46) Once the gradient waveform G(t) has been determined, data defining the encoding sequence and the gradient waveform may be supplied to the scanner 1 wherein the dMRI measurement may be performed.
(47) The inventive manner of designing the gradient waveform G is compatible with planar, prolate as well as spherical b-tensor encoding. For instance, the encoding gradient G may be of such a form that a diffusion encoding tensor representation b has at least two non-zero eigenvalues. As is known in the art, the tensor b is given by: b=∫.sub.0.sup.TEq(t)q(t).sup.Tdt, where q(t)=γ∫.sub.0.sup.tG(t′)dt′ is a dephasing vector. The encoding gradient G may be configured such that the tensor b has exactly three non-zero eigenvalues, thereby providing 3D diffusion encoding. Both anisotropic and isotropic 3D diffusion encoding is possible, where anistropic 3D encoding implies that the eigenvalues of the tensor b are not all equal and isotropic 3D encoding implies that the eigenvalues of the tensor b are all equal.
(48) In the above the inventive concept has mainly been described with reference to a limited number of examples. However, as is readily appreciated by a person skilled in the art, other examples than the ones disclosed above are equally possible within the scope of the inventive concept, as defined by the appended claims.
REFERENCES
(49) Here follows a listing of the references referred to in the above: Baron, C. A., Lebel, R. M., Wilman, A. H. & Beaulieu, C. 2012. The effect of concomitant gradient fields on diffusion tensor imaging. Magnetic Resonance in Medicine, 68, 1190-201. Bernstein, M. A., Zhou, X. J., Polzin, J. A., King, K. F., Ganin, A., Pelc, N. J., & Glover, G. H. (1998). Concomitant gradient terms in phase contrast MR: analysis and correction. Magn Reson Med, 39(2), 300-308. Lasic̆, S., Szczepankiewicz, F., Eriksson, S., Nilsson, M. & Topgaard, D. 2014. Microanisotropy imaging: quantification of microscopic diffusion anisotropy and orientational order parameter by diffusion MRI with magic-angle spinning of the q-vector. Frontiers in Physics, 2, 11. Meier, C., Zwanger, M., Feiweier, T., & Porter, D. (2008). Concomitant field terms for asymmetric gradient coils: consequences for diffusion, flow, and echo-planar imaging. Magn Reson Med, 60(1), 128-134. http://doi.org/10.1002/mrm.21615 Nilsson, M., van Westen, D., Ståhlberg, F., Sundgren, P. C., & Lätt, J. (2013). The role of tissue microstructure and water exchange in biophysical modelling of diffusion in white matter. Magnetic Resonance Materials in Physics, Biology and Medicine, 26(4), 345-370.
(50) Norris, D. G. & Hutchison, J. M. S. 1990. Concomitant magnetic field gradients and their effects on imaging at low magnetic field strengths. Magnetic Resonance Imaging, 8, 33-37. Sjölund, J., Szczepankiewicz, F., Nilsson, M., Topgaard, D., Westin, C. F. & Knutsson, H. 2015. Constrained optimization of gradient waveforms for generalized diffusion encoding. Journal of Magnetic Resonance, 261, 157-168. Szczepankiewicz, F., Van Westen, D., Englund, E., Westin, C. F., Stahlberg, F., Latt, J., Sundgren, P. C. & Nilsson, M. 2016. The link between diffusion MRI and tumor heterogeneity: Mapping cell eccentricity and density by diffusional variance decomposition (DIVIDE). Neuroimage, 142, 522-532. Westin, C.-F., Knutsson, H., Pasternak, O., Szczepankiewicz, F., Ozarslan, E., van Westen, D., et al. (2016). Q-space trajectory imaging for multidimensional diffusion MRI of the human brain. NeuroImage, 135(C), 345-362. Zhou, X. J., Du, Y. P., Bernstein, M. A., Reynolds, H. G., Maier, J. K. & Polzin, J. A. 1998. Concomitant magnetic-field-induced artifacts in axial echo planar imaging. Magnetic Resonance in Medicine, 39, 596-605.