- Email: [email protected]
Stochastic NMR imaging
JOURNAL
OF MAGNETIC
RESONANCE
60.37-45
(1984)
Stochastic NMR Imaging B. BLUMICH Makromolekulare
Chemie II, Universitiit Bayreuth, Postfach 3008, 8580 Bayreuth, Federal Republic of Germany
Received February 24, 1984; revised April 20, 1984 NMR imaging is treated in terms of white-noise analysis of linear systems. A stochastic rf input signal excites a broad frequency region at low input power. One possible realization of stochastic three-dimensional imaging with random field gradient modulation is supported by a numerical simulation. 0 1984 Academic Req Inc.
The imaging experiment (1) will be treated in terms of system analysis. The underlying concept is that of a black box which is tested with an input signal s(t) and responds with an output signal y(t). From a comparison of input and output signals the features of the black box are sought to be revealed. This concept is termed correlation(2). It is successfully applied in many areas of science, including chemistry (3) and NMR spectroscopy (4, 5). In NMR imaging the feature to be revealed is the spin density distribution of the sample as a function of space coordinates. To achieve this goal it is necessary to apply time-dependent magnetic field gradients across the sample in addition to a radiofrequency input signal. Thus there are four variables in an imaging experiment, the shape of the radiofrequency input signal, and the modulations of magnetic x, y, and z gradients. Consequently the imaging experiment may be described in terms of a four-input,time-invariantsystem(2). Alternatively, however, it can be considered a one-input,time-dependent system,where the time dependence is imposed by the magnetic field gradients and the one input is the rf signal. The latter point of view is chosen below. More specifically, the NMR imaging experiment will be described in terms of a linear timedependent system. Stochastic excitation denotes a noise-modulated rf input signal, which is supplied for the entire length of the imaging experiment. Radiofrequency noise is conveniently generated by a sequence of rf pulses with random flip angles (5). The pulses are separated by the sampling interval At. For each input pulse one value of the response is sampled. The magnetic field gradients are varied independently of the rf input. A particular scheme with random field gradient modulation is worked out below. The stochastic rf excitation delivers its excitation energy to the system within a characteristic relaxation time Ti, while pulsed excitation concentrates the same energy into a single pulse of length 6t. Therefore stochastic excitation requires less power by a factor T,/& than does pulsed exit&ion. The anticipated power saving of a few orders of magnitude justifies elaboration of the concept of stochastic NMR imaging. 37
0022-23&l/84
$3.00
Copy&M 0 1984 by Academic FVem, Inc. All rights of reproduction in any form i-estrvcd.
38
B. BLUMICH
The spin density is retrieved from cross-correlation (2) of input and output records. The longer these records, the better will be the achievable signal-to-noise ratio of the spin density. Since no rigid pulse timing scheme must be followed, every additional input/output data pair will increase the signal-to-noise ratio of the entire image, provided the excitation values are linearily independent. If the spin density is not stationary, such as in living systems, the data records to be used in the cross-correlation can be restricted to those, during the acquisition of which the spin density was approximately in the same state (6). THE NMR IMAGING
SYSTEM
To gain insight into stochastic NMR imaging, ,the linear response function for arbitrary rf excitation and field gradient modulation will be investigated. It is obtained from the impulse response function h(t) of a magnetic spin system in the presence of time-dependent magnetic field gradients in the rotating frame (7)
with h(t < 0) = 0. The spin density distribution denotes the space coordinates
is designated by p(r), where r
The transverse relaxation time is denoted by T2 and y is the magneto-gyric The time-dependent magnetic field gradient is given by g(f) = VB#)
=
g&? i 1 gy(t’> . gzw
ratio.
PI
The gradient function g(f) is integrated in Eq. [l] from the radiofrequency pulse at time t’ = 0 to time t’ = t. Thus the integral depends on the values of the field gradient modulation after the pulse. Assuming for the moment an uninterrupted random modulation then it is obvious, that the experimentally determined shape of the impulse response depends on the particular time at which the pulse is applied. For convenience, this time has been set to zero in Eq. [l]. The formulation [l] of the impulse response thus demonstrates, that the imaging system can indeed be considered a timedependent system where the time dependence is imposed by the field gradient modulation. For general rf input the linear system response y(t) is the convolution of the impulse response h(t) with the rf input signal s(t). IJsing the abbreviation
YsI f-7
[31
STOCHASTIC
NMR
39
IMAGING
in [ 11, the system response at time t is given by m YW =
s0
sss
p(x , y, z)e-T/T2eixu(t,T)eiyv(t,r)errw(r.r) &dydzs(t
- 7)dTe
141
This equation describes the linear response of the imaging system. It can be verified by choosing a delta pulse as input signal, s(t) = 6(t). In this case the impulse response, Eq. [l], results. RETRIEVAL
OF
THE
SPIN
DENSITY
In stochastic system analysis the input signal usually is white noise with zero mean so that the autocorrelation of the real-valued input is proportional to a delta function frnW f 07s(t)s(t - 0)dt = &(a). s
PI
The proportionality constant pz is the second moment or the power of the input. For time-invariant systems the linear impulse response function can be retrieved from cross-correlation of input and output signals according to FT + ,’ y(t)s(t - a)dt = &z(a). s
161
The cross-correlation is a time average of the product of the system response with the delayed input signal. The result is the linear impulse response as a function of the delay u (2). In case of time-dependent systems this concept cannot be applied. However, for the imaging system [4] the goal is different. Not the complete time-dependent impulse response function is wanted, but only its time-invariant part, the spin density. Since the time-dependent part is known, the system response y(t) can be correlated with the product of all the timedependent factors in the response integral [41 P(t - u) = e-ir’u(r,r)e-iv’u(r,u)e-iz’w(t,o)S(t _ u)e i71 The primed coordinates as well as u become variables in the cross-correlation. Its evaluation yields the spin density as a function of the cross-correlation variables: T
lim -!
~-rnT
so Y(Ola - uw
=(2r))r2SSS~x,y,Z)Sgme-‘/T*~x-x?6Cv-y’)6(z-z? X 6(~ - ujdrdxdydz = (2rj3p2p(xI, y’, z?e+‘*.
[8b]
[W
40
B. BLiiMICH
In Eq. [8a] the cross-correlation has been written down explicitly but the order of integration has been changed. The step from [8a] to [8b] involves the evaluation of the innermost integral. The fourfold product of delta functions results from Eq. [5], the definition of the delta function 1 * ~‘w-“o’& G s -03
= (j(x - xg)
]91
and a postulated orthogonality of the four function,s in the product of Eq. [7]. The step from Eq. [8b] to [8c] is the evaluation of the delta functions. For the orthogonality to be satisfied three conditions must be met: (1) The input signal s(t) must have zero mean:
(2) The time integrals over the magnetic field gradients must always be different from zero: u(t, a), v(t, a), w(t, a) # 0
1111
for all t and some given u. (3) The difference of the local phase angles in different directions x, y, z should never vanish for all primed and unprimed values of x, y, and z xu(t, a) - yu(t, a) f 0, yu(t, u) - zw(t, a) # 0, zw(t, a) - xu(t, a) # 0,
1121
for all times t, the u given, and all x, y, z. Condition 1 is trivial. Condition 2 can be fulfilled by modulating a field gradient such that no gradient component ever changes sign. The last condition is the most restricting one. One way to observe it is by using significantly different field gradients in each direction x, y, and z. A different gradient scheme is worked out below for the case where the spin density shall be recovered at discrete coordinates x’, y’, z’. The final result of the cross-correlation [8] is the spin density as a function of space coordinates, scaled by the input power, LQ, and a negative exponential involving the transverse relaxation time T2 and the cross-correlation delay u. Though, formally, any delay can be admitted, it is obvious, that when choosing u too big, the amplitude of the result will be too small for practical needs, while choosing u too small, the orthogonality condition 3 will be hard to fulfil. In fact, for u close to zero, magnetic field gradients will hardly affect the outcome of the cross-correlation at all and no spatial resolution can be achieved. Thus u should be set to some value close to Tz. Also it should be noted that the input power, c(~, cannot be made arbitrarily large in order to increase the system response, since eventually the system will become nonlinear and the response will become damped. However, there exists an input
STOCHASTIC
NMR IMAGING
41
power level which yields a response maximum and consequently optimum signalto-noise ratio in the acquisition of the response. In addition to scaling the system response the power level can be used to select the degree of saturation of the spins and in this way introduce T, contrast in the spin density. Since all T, effects are nonlinear, the T, contrast cannot be deduced from the treated linear model. A PARTICULAR
Discretization
SCHEME
of the Spin Density
For numerical data processing the orthogonality conditions should be satisfied with discrete functions for excitation s, field gradient integrals u, u, w and space coordinates x’, y’, z’. The excitation s can be realized with either a white sequence of pulses or continuous noise with zero mean. The particular distribution function is of no influence for evaluation of a l&a&y driven system. Contrary to continuous excitation, pulsed excitation enables time sharing to prevent receiver overload. In particular a binary sequence of flip angles will be sufhcient. Such sequences are easy to generate with feedback shift registers (5, 8). The second orthogonal&y condition [l l] is fulfilled by choosing distribution functions for the field gradients g,, g,,, g, with nonzero mean. For example, a symmetrical distribution function will do for each gradient component gi such that for the mean values i = x, y, z. ii # 0, ]131 In consideration of the third orthogonality condition [ 121 it will be assumed for now that the measured response derives from a spin density, which, contrary to reality, is defined only on a discrete set of equidistant coordinate points x, y, z. This set of points shall be the same as the one on which the spin density p(x’, y’, z’) will be reconstructed by cross-correlation. A continuous spin density of the sample will be considered later. To avoid overlap of values, an interleaving set of discrete coordinates is constructed by setting x’ = (nx + t&, Y’ = by + +& z’ = (n, + c&z,
[141
where ni is an integer, 0 < ti < 1 for i = x, y, z, and the dimension of a discrete spin density element or picture cell is specified by (2. Similarly a discrete range of values is introduced for the gradient integrals
24= (nu+ &, 2,= b” + C”k, w = (n, + %k,
1151
where ni is an integer, 0 < 9 < 1 for i = U, u, w, and the product of the minimum field gradient step and the gradient switching interval At is denoted by g. Without
42
B. BLijMICH
loss of generality the field gradient switching interval will be assumed to be identical to the sampling interval; hence the same symbol. The unit of g is tad/m. The value of g measures the minimum precession phase per picture cell which is acquired during one gradient step of length At. For a symmetrical distribution of gradient numbers ni, the mean field gradients are gi = cigjdt. Evaluation of one of equations [ 121 provides insight into a choice for tx, E,,, t, and c,,, cv, E,+Z x’u - y’u = [(nx + &z,
+ EJ - (ny + E&0(” + f”)] *a ‘g.
[I61
The orthogonality condition states, that the right-hand side must be different from zero for all values of n,, n,,, n, and n,, n,, n,. Siting fx = cv, eY = E”, and E, = t, one obtains x’u - y’u = [n,n, + (n, + n&
+ E: - ?&A,,- (n, + rz”)Q - c;] * a ‘g.
[I71
Similar expressions result for the other two equations in [ 121. With 6x =
t,
= l/7,
ty = t, = 317, cz = t, = 517,
[I81
it is easily verified, that the third orthogonality condition is fulfilled by the interleaving coordinate and gradient integral values according to Eqs. [ 141, [ 151, and [18].
Numerical Simulation Bearing in mind that the validity of the third orthogonality condition is based on the hypothesis of a discrete spin density p(x, y, z) in the sample, the results derived so far have been tested by a computer simulation. The linear response has been calculated according to Eqs. [4], [ 141, [ 151, and [18]. The spin density was defined on a 3 X 3 X 3 mesh (Fig. 1) with values reflecting the proton densities of living tissues (9). The transverse relaxation time was T2 = 0.1 s (10). To keep the computation times at an acceptable level, the product of spatial resolution and precession phase was set to ug = 2.5 rad. The sampling interval was 0.01 s. For each sampling interval new field gradient values and a new excitation value were derived from a 20 stage binary feedback shift register (5, 8). The excitation could assume two values s = [-l/2, +1/2] with equal probability, while the field gradient numbers could assume four values n,, n,,, n, = [-3/2, -l/2, +1/2, +3/2] with probabilities :l/8, 3/8, 3/8, and l/8. The field gradient numbers were obtained by shifting the feedback register thrice and adding the binary outputs of each shift. Simulation of 65,000 response values required 437 min on a VAX 780 computer. Execution of the cross-correlation [8] needed 39 min. The cross-correlation was executed for a 5 X 5 X 5 point mesh of the spin density p(x’, y’, z’) and a time shift u = T2. The result is shown in Fig. 2 in terms of slices
STOCHASTIC
43
NMR IMAGING
+a
5a 7 i
--a72 a
RG. I. The spin density p(x, y, z) used in the simulation. It is defined as a 3 X 3 X 3 array of cubic picture cells with dimension u.
showing y’z’ planes for different values of x’. The numbers denote the retrieved spin density values in percent. The agreement with the input values, Fig. 1, is satisfactory. The C-shaped elements with high density are clearly identified and the average density of these elements correctly decreases with increasing x’. The background noise results from the finite sample record length. It can be reduced in several ways; for instance, by (1) optimization of the phase angle step ug, (2) evaluation of longer sample records, (3) averaging of spin density values for different values of the crosscorrelation delay 6, and (4) optimization of gradient modulation and excitation statistics.
3 3 2 3 -1 9 -12 7 4 -6 -7 0 -1 1 -7 ,I=-fi
-4 -1 -5 0 -3 0 12 -1 2 3
7
Q
0 2 -6 -9 4 10 77 74 a2 1 -9 79 13 11 7 2 a6 94 aa -6 0 7 4 4 0 *C-$a 7
2
-1 5 -7 0 -2 7 5 0 6 -2 0 0 0 -2 7
1 -2 o-2 2 0 2 4 0 -7
Y
I
j&+150
7
7-11 -3 62 4 a4 9 77 4 -3
4 6 0 71 80 5 12 lo -5 73 77 -a -7 3 6 XL+ 10 7
-5 4 5 -5 0 8 59 76 70 -2 -17 73 12 12 -a 4 69 62 77 -5 6-5 0 0 5 Xk++J
FIG. 2. Slices of the spin density recovered by cross-correlation of 65,OtlO input/output data pairs. Depicted are y’z’ planes for five different vahtes of x’. The reconstructed image volume covered 5 X 5 X 5 picture cells.
44
B. BLiiMICH
The numerical simulation proves that the above scheme with interleaving coordinate and gradient values is applicable to three-dimensional imaging, provided that the spin density of the sample is given as an army of discrete points.
ContinuousSampleDensity For samples with continuous spin density p(x, y, z) the results obtained so far must be corrected. The orthogonality condition 3 will not be fulfdled anymore, since there always will be some values of the continuous coordinates x, y, z which match the discrete coordinates x’, y’, z’ such that for instance x’u - yu = 0, y’u - zw = 0,
or
z’w- xu = 0,
x’u - zw = 0, y’u - .ru = 0,
z’w- yu = 0.
1191
In these cases ghost signals of p(y, z, x) or p(z, x, y) will arise at coordinates x’, y’, z’ in addition to the wanted signal of p(x, y, z). Noting that the gradients g,, g,, g, are time-dependent random variables, it can be concluded from the central limit theorem (2), that the gradient integrals U, u, w approximate Gaussian variables. Thus the ghost signals will be smeared out with Gaussian probability, while the wanted signals will add coherently in the crosscorrelation. The variance of the ghost signals will be the larger, the more independent gradient values are summed up in the integrals (2). Consequently, the image quality is expected to gain with the speed of gradient modulation. The ghost signals may be broadened even Wther by shifting the origin of the field gradients through the sample during data acquisition. SUMMARY
The linear response of a NMR imaging system has been formulated and conditions were derived under which the spin density can be calculated from crosscorrelation of input and response signals. For continuous input signals such as white noise, radiofiequency power demands are considerably less than for sparsely pulsed input. Due to the random nature of the excitation no deterministic pulse timing scheme must be followed, and every single data point contributes toward the quality of the entire image. A particular three-dimensional imaging scheme has been outlined, which in addition to stochastic radiofrequency input utilizes random modulation of magnetic field gradients. Violation of the third orthogonality condition by this scheme leads to a washed out background of the reconstructed spin density. Other stochastic imaging schemes can be devised. To optimize the image quality, different excitation statistics and gradient modulation functions need to be investigated. The resolution of the final image is determined by the step width of the crosscorrelation variables x’, y’, and z’. Similarly the TZ contrast is determined by the particular value chosen for the cross-correlation delay u. Consequently images with different resolution and different TZ contrast can be computed from one and the same set of raw data.
STOCHASTIC
NMR IMAGING
45
The cross-correlation can be viewed in terms of a transformation of the timedependent system response into three-dimensional x’-y?-z’ space. For practical applications of the method it will be worthwhile to search for a fast transformation algorithm in connection with deterministic or pseudo-random radiofrequency and field gradient modulation in order to replace the rather time-consuming crosscorrelation used for reconstruction of the spin density. An alternative route to digital data processing is the analog computation of the cross-correlation. It is possible, at least in principle, to build one hardware crosscorrelation channel per picture element. An array of such channels can be combined to form a two-dimensional window which would enable real-time imaging. The only time delay for image construction in this case would be the integration time of the cross-correlation. ACKNOWLEDGMENTS The first part of this work was done at the University of New Brunswick in Fredericton, Canada, with support of a fellowship by the DAAD (Deutscher Akademischer Austauschdienst) under the research committee of NATO. Continuation of the project was made possible by a research stipend of the DFG (Deutsche ForschungsgemeinschaB) at the University of Bayreuth, Germany. Stimulating discussions with Professor R. Kaiser are gratefully acknowledged. REFERENCES 1. (a) P. MANSFIELD AND P. G. MORRIS, “Advances in Magnetic Resonances,” Suppl. 2, “NMR Imaging in Biomedicine,” Academic Press, New York, 1982; (b) L. KAUFMAN, L. E. CROOKS, AND A. R. MARGULIS, eds., “Nuclear Magnetic Resonance Imaging in Medicine,” Igaku-Shoin, New York, 1981; (c) P. A. BOTTOMLEY, Rev. Sci. Instrum. 53, 1319 (1983). 2. J. S. BENDAT AND A. G. PIERSOL, “Random Data: Analysis and Measurement Procedures,” WileyInterscience, New York, 1971. 3. (a) G. HORLIK AND G. M. HIEFTIE, Contemp. Top. Anal. Clin. Chem. 3, 153 (1978); (b) G. M. HIEFTJE,Am. Lab. 13, 76 (1981). 4. (a) R. R. ERNST, J. Mugs. Reson. 3, 10 (1970); (b) R. KAISER, J. Mugs. Reson. 3, 28 (1970); (c) J. DADOK AND R. F. SPRECHER,J. Mugs. Reson. 13,243 (1974); (d) B. BLUMICH AND D. ZIESSOW, J. Chem. Phys. 78, 1059 (1983). 5. B. BLUMICH AND D. Zmssow, J. Mugn. Reson. 52,42, (1983). 6. D. I. HOULT, private communication. 7. L. F. FEINER AND P. R. LOCHER, Appl. Phys. 22,257 (1980). 8. G. A. KORN, “Random Process, Simulation and Measurement,” p. 83 If, McGraw-Hill, New York, 1966. 9. P. MANSFIELD AND I. L. PYKETT, J. Mugn. Reson. 29,355 (1978). 10. J. D. DE CERTAINES, L. BENOIST, F. DARCEL-MENAULT, M. CHATEL, AND A. M. BERNARD, Bull. Magn. Reson. 5, 267 (1983).
OF MAGNETIC
RESONANCE
60.37-45
(1984)
Stochastic NMR Imaging B. BLUMICH Makromolekulare
Chemie II, Universitiit Bayreuth, Postfach 3008, 8580 Bayreuth, Federal Republic of Germany
Received February 24, 1984; revised April 20, 1984 NMR imaging is treated in terms of white-noise analysis of linear systems. A stochastic rf input signal excites a broad frequency region at low input power. One possible realization of stochastic three-dimensional imaging with random field gradient modulation is supported by a numerical simulation. 0 1984 Academic Req Inc.
The imaging experiment (1) will be treated in terms of system analysis. The underlying concept is that of a black box which is tested with an input signal s(t) and responds with an output signal y(t). From a comparison of input and output signals the features of the black box are sought to be revealed. This concept is termed correlation(2). It is successfully applied in many areas of science, including chemistry (3) and NMR spectroscopy (4, 5). In NMR imaging the feature to be revealed is the spin density distribution of the sample as a function of space coordinates. To achieve this goal it is necessary to apply time-dependent magnetic field gradients across the sample in addition to a radiofrequency input signal. Thus there are four variables in an imaging experiment, the shape of the radiofrequency input signal, and the modulations of magnetic x, y, and z gradients. Consequently the imaging experiment may be described in terms of a four-input,time-invariantsystem(2). Alternatively, however, it can be considered a one-input,time-dependent system,where the time dependence is imposed by the magnetic field gradients and the one input is the rf signal. The latter point of view is chosen below. More specifically, the NMR imaging experiment will be described in terms of a linear timedependent system. Stochastic excitation denotes a noise-modulated rf input signal, which is supplied for the entire length of the imaging experiment. Radiofrequency noise is conveniently generated by a sequence of rf pulses with random flip angles (5). The pulses are separated by the sampling interval At. For each input pulse one value of the response is sampled. The magnetic field gradients are varied independently of the rf input. A particular scheme with random field gradient modulation is worked out below. The stochastic rf excitation delivers its excitation energy to the system within a characteristic relaxation time Ti, while pulsed excitation concentrates the same energy into a single pulse of length 6t. Therefore stochastic excitation requires less power by a factor T,/& than does pulsed exit&ion. The anticipated power saving of a few orders of magnitude justifies elaboration of the concept of stochastic NMR imaging. 37
0022-23&l/84
$3.00
Copy&M 0 1984 by Academic FVem, Inc. All rights of reproduction in any form i-estrvcd.
38
B. BLUMICH
The spin density is retrieved from cross-correlation (2) of input and output records. The longer these records, the better will be the achievable signal-to-noise ratio of the spin density. Since no rigid pulse timing scheme must be followed, every additional input/output data pair will increase the signal-to-noise ratio of the entire image, provided the excitation values are linearily independent. If the spin density is not stationary, such as in living systems, the data records to be used in the cross-correlation can be restricted to those, during the acquisition of which the spin density was approximately in the same state (6). THE NMR IMAGING
SYSTEM
To gain insight into stochastic NMR imaging, ,the linear response function for arbitrary rf excitation and field gradient modulation will be investigated. It is obtained from the impulse response function h(t) of a magnetic spin system in the presence of time-dependent magnetic field gradients in the rotating frame (7)
with h(t < 0) = 0. The spin density distribution denotes the space coordinates
is designated by p(r), where r
The transverse relaxation time is denoted by T2 and y is the magneto-gyric The time-dependent magnetic field gradient is given by g(f) = VB#)
=
g&? i 1 gy(t’> . gzw
ratio.
PI
The gradient function g(f) is integrated in Eq. [l] from the radiofrequency pulse at time t’ = 0 to time t’ = t. Thus the integral depends on the values of the field gradient modulation after the pulse. Assuming for the moment an uninterrupted random modulation then it is obvious, that the experimentally determined shape of the impulse response depends on the particular time at which the pulse is applied. For convenience, this time has been set to zero in Eq. [l]. The formulation [l] of the impulse response thus demonstrates, that the imaging system can indeed be considered a timedependent system where the time dependence is imposed by the field gradient modulation. For general rf input the linear system response y(t) is the convolution of the impulse response h(t) with the rf input signal s(t). IJsing the abbreviation
YsI f-7
[31
STOCHASTIC
NMR
39
IMAGING
in [ 11, the system response at time t is given by m YW =
s0
sss
p(x , y, z)e-T/T2eixu(t,T)eiyv(t,r)errw(r.r) &dydzs(t
- 7)dTe
141
This equation describes the linear response of the imaging system. It can be verified by choosing a delta pulse as input signal, s(t) = 6(t). In this case the impulse response, Eq. [l], results. RETRIEVAL
OF
THE
SPIN
DENSITY
In stochastic system analysis the input signal usually is white noise with zero mean so that the autocorrelation of the real-valued input is proportional to a delta function frnW f 07s(t)s(t - 0)dt = &(a). s
PI
The proportionality constant pz is the second moment or the power of the input. For time-invariant systems the linear impulse response function can be retrieved from cross-correlation of input and output signals according to FT + ,’ y(t)s(t - a)dt = &z(a). s
161
The cross-correlation is a time average of the product of the system response with the delayed input signal. The result is the linear impulse response as a function of the delay u (2). In case of time-dependent systems this concept cannot be applied. However, for the imaging system [4] the goal is different. Not the complete time-dependent impulse response function is wanted, but only its time-invariant part, the spin density. Since the time-dependent part is known, the system response y(t) can be correlated with the product of all the timedependent factors in the response integral [41 P(t - u) = e-ir’u(r,r)e-iv’u(r,u)e-iz’w(t,o)S(t _ u)e i71 The primed coordinates as well as u become variables in the cross-correlation. Its evaluation yields the spin density as a function of the cross-correlation variables: T
lim -!
~-rnT
so Y(Ola - uw
=(2r))r2SSS~x,y,Z)Sgme-‘/T*~x-x?6Cv-y’)6(z-z? X 6(~ - ujdrdxdydz = (2rj3p2p(xI, y’, z?e+‘*.
[8b]
[W
40
B. BLiiMICH
In Eq. [8a] the cross-correlation has been written down explicitly but the order of integration has been changed. The step from [8a] to [8b] involves the evaluation of the innermost integral. The fourfold product of delta functions results from Eq. [5], the definition of the delta function 1 * ~‘w-“o’& G s -03
= (j(x - xg)
]91
and a postulated orthogonality of the four function,s in the product of Eq. [7]. The step from Eq. [8b] to [8c] is the evaluation of the delta functions. For the orthogonality to be satisfied three conditions must be met: (1) The input signal s(t) must have zero mean:
(2) The time integrals over the magnetic field gradients must always be different from zero: u(t, a), v(t, a), w(t, a) # 0
1111
for all t and some given u. (3) The difference of the local phase angles in different directions x, y, z should never vanish for all primed and unprimed values of x, y, and z xu(t, a) - yu(t, a) f 0, yu(t, u) - zw(t, a) # 0, zw(t, a) - xu(t, a) # 0,
1121
for all times t, the u given, and all x, y, z. Condition 1 is trivial. Condition 2 can be fulfilled by modulating a field gradient such that no gradient component ever changes sign. The last condition is the most restricting one. One way to observe it is by using significantly different field gradients in each direction x, y, and z. A different gradient scheme is worked out below for the case where the spin density shall be recovered at discrete coordinates x’, y’, z’. The final result of the cross-correlation [8] is the spin density as a function of space coordinates, scaled by the input power, LQ, and a negative exponential involving the transverse relaxation time T2 and the cross-correlation delay u. Though, formally, any delay can be admitted, it is obvious, that when choosing u too big, the amplitude of the result will be too small for practical needs, while choosing u too small, the orthogonality condition 3 will be hard to fulfil. In fact, for u close to zero, magnetic field gradients will hardly affect the outcome of the cross-correlation at all and no spatial resolution can be achieved. Thus u should be set to some value close to Tz. Also it should be noted that the input power, c(~, cannot be made arbitrarily large in order to increase the system response, since eventually the system will become nonlinear and the response will become damped. However, there exists an input
STOCHASTIC
NMR IMAGING
41
power level which yields a response maximum and consequently optimum signalto-noise ratio in the acquisition of the response. In addition to scaling the system response the power level can be used to select the degree of saturation of the spins and in this way introduce T, contrast in the spin density. Since all T, effects are nonlinear, the T, contrast cannot be deduced from the treated linear model. A PARTICULAR
Discretization
SCHEME
of the Spin Density
For numerical data processing the orthogonality conditions should be satisfied with discrete functions for excitation s, field gradient integrals u, u, w and space coordinates x’, y’, z’. The excitation s can be realized with either a white sequence of pulses or continuous noise with zero mean. The particular distribution function is of no influence for evaluation of a l&a&y driven system. Contrary to continuous excitation, pulsed excitation enables time sharing to prevent receiver overload. In particular a binary sequence of flip angles will be sufhcient. Such sequences are easy to generate with feedback shift registers (5, 8). The second orthogonal&y condition [l l] is fulfilled by choosing distribution functions for the field gradients g,, g,,, g, with nonzero mean. For example, a symmetrical distribution function will do for each gradient component gi such that for the mean values i = x, y, z. ii # 0, ]131 In consideration of the third orthogonality condition [ 121 it will be assumed for now that the measured response derives from a spin density, which, contrary to reality, is defined only on a discrete set of equidistant coordinate points x, y, z. This set of points shall be the same as the one on which the spin density p(x’, y’, z’) will be reconstructed by cross-correlation. A continuous spin density of the sample will be considered later. To avoid overlap of values, an interleaving set of discrete coordinates is constructed by setting x’ = (nx + t&, Y’ = by + +& z’ = (n, + c&z,
[141
where ni is an integer, 0 < ti < 1 for i = x, y, z, and the dimension of a discrete spin density element or picture cell is specified by (2. Similarly a discrete range of values is introduced for the gradient integrals
24= (nu+ &, 2,= b” + C”k, w = (n, + %k,
1151
where ni is an integer, 0 < 9 < 1 for i = U, u, w, and the product of the minimum field gradient step and the gradient switching interval At is denoted by g. Without
42
B. BLijMICH
loss of generality the field gradient switching interval will be assumed to be identical to the sampling interval; hence the same symbol. The unit of g is tad/m. The value of g measures the minimum precession phase per picture cell which is acquired during one gradient step of length At. For a symmetrical distribution of gradient numbers ni, the mean field gradients are gi = cigjdt. Evaluation of one of equations [ 121 provides insight into a choice for tx, E,,, t, and c,,, cv, E,+Z x’u - y’u = [(nx + &z,
+ EJ - (ny + E&0(” + f”)] *a ‘g.
[I61
The orthogonality condition states, that the right-hand side must be different from zero for all values of n,, n,,, n, and n,, n,, n,. Siting fx = cv, eY = E”, and E, = t, one obtains x’u - y’u = [n,n, + (n, + n&
+ E: - ?&A,,- (n, + rz”)Q - c;] * a ‘g.
[I71
Similar expressions result for the other two equations in [ 121. With 6x =
t,
= l/7,
ty = t, = 317, cz = t, = 517,
[I81
it is easily verified, that the third orthogonality condition is fulfilled by the interleaving coordinate and gradient integral values according to Eqs. [ 141, [ 151, and [18].
Numerical Simulation Bearing in mind that the validity of the third orthogonality condition is based on the hypothesis of a discrete spin density p(x, y, z) in the sample, the results derived so far have been tested by a computer simulation. The linear response has been calculated according to Eqs. [4], [ 141, [ 151, and [18]. The spin density was defined on a 3 X 3 X 3 mesh (Fig. 1) with values reflecting the proton densities of living tissues (9). The transverse relaxation time was T2 = 0.1 s (10). To keep the computation times at an acceptable level, the product of spatial resolution and precession phase was set to ug = 2.5 rad. The sampling interval was 0.01 s. For each sampling interval new field gradient values and a new excitation value were derived from a 20 stage binary feedback shift register (5, 8). The excitation could assume two values s = [-l/2, +1/2] with equal probability, while the field gradient numbers could assume four values n,, n,,, n, = [-3/2, -l/2, +1/2, +3/2] with probabilities :l/8, 3/8, 3/8, and l/8. The field gradient numbers were obtained by shifting the feedback register thrice and adding the binary outputs of each shift. Simulation of 65,000 response values required 437 min on a VAX 780 computer. Execution of the cross-correlation [8] needed 39 min. The cross-correlation was executed for a 5 X 5 X 5 point mesh of the spin density p(x’, y’, z’) and a time shift u = T2. The result is shown in Fig. 2 in terms of slices
STOCHASTIC
43
NMR IMAGING
+a
5a 7 i
--a72 a
RG. I. The spin density p(x, y, z) used in the simulation. It is defined as a 3 X 3 X 3 array of cubic picture cells with dimension u.
showing y’z’ planes for different values of x’. The numbers denote the retrieved spin density values in percent. The agreement with the input values, Fig. 1, is satisfactory. The C-shaped elements with high density are clearly identified and the average density of these elements correctly decreases with increasing x’. The background noise results from the finite sample record length. It can be reduced in several ways; for instance, by (1) optimization of the phase angle step ug, (2) evaluation of longer sample records, (3) averaging of spin density values for different values of the crosscorrelation delay 6, and (4) optimization of gradient modulation and excitation statistics.
3 3 2 3 -1 9 -12 7 4 -6 -7 0 -1 1 -7 ,I=-fi
-4 -1 -5 0 -3 0 12 -1 2 3
7
Q
0 2 -6 -9 4 10 77 74 a2 1 -9 79 13 11 7 2 a6 94 aa -6 0 7 4 4 0 *C-$a 7
2
-1 5 -7 0 -2 7 5 0 6 -2 0 0 0 -2 7
1 -2 o-2 2 0 2 4 0 -7
Y
I
j&+150
7
7-11 -3 62 4 a4 9 77 4 -3
4 6 0 71 80 5 12 lo -5 73 77 -a -7 3 6 XL+ 10 7
-5 4 5 -5 0 8 59 76 70 -2 -17 73 12 12 -a 4 69 62 77 -5 6-5 0 0 5 Xk++J
FIG. 2. Slices of the spin density recovered by cross-correlation of 65,OtlO input/output data pairs. Depicted are y’z’ planes for five different vahtes of x’. The reconstructed image volume covered 5 X 5 X 5 picture cells.
44
B. BLiiMICH
The numerical simulation proves that the above scheme with interleaving coordinate and gradient values is applicable to three-dimensional imaging, provided that the spin density of the sample is given as an army of discrete points.
ContinuousSampleDensity For samples with continuous spin density p(x, y, z) the results obtained so far must be corrected. The orthogonality condition 3 will not be fulfdled anymore, since there always will be some values of the continuous coordinates x, y, z which match the discrete coordinates x’, y’, z’ such that for instance x’u - yu = 0, y’u - zw = 0,
or
z’w- xu = 0,
x’u - zw = 0, y’u - .ru = 0,
z’w- yu = 0.
1191
In these cases ghost signals of p(y, z, x) or p(z, x, y) will arise at coordinates x’, y’, z’ in addition to the wanted signal of p(x, y, z). Noting that the gradients g,, g,, g, are time-dependent random variables, it can be concluded from the central limit theorem (2), that the gradient integrals U, u, w approximate Gaussian variables. Thus the ghost signals will be smeared out with Gaussian probability, while the wanted signals will add coherently in the crosscorrelation. The variance of the ghost signals will be the larger, the more independent gradient values are summed up in the integrals (2). Consequently, the image quality is expected to gain with the speed of gradient modulation. The ghost signals may be broadened even Wther by shifting the origin of the field gradients through the sample during data acquisition. SUMMARY
The linear response of a NMR imaging system has been formulated and conditions were derived under which the spin density can be calculated from crosscorrelation of input and response signals. For continuous input signals such as white noise, radiofiequency power demands are considerably less than for sparsely pulsed input. Due to the random nature of the excitation no deterministic pulse timing scheme must be followed, and every single data point contributes toward the quality of the entire image. A particular three-dimensional imaging scheme has been outlined, which in addition to stochastic radiofrequency input utilizes random modulation of magnetic field gradients. Violation of the third orthogonality condition by this scheme leads to a washed out background of the reconstructed spin density. Other stochastic imaging schemes can be devised. To optimize the image quality, different excitation statistics and gradient modulation functions need to be investigated. The resolution of the final image is determined by the step width of the crosscorrelation variables x’, y’, and z’. Similarly the TZ contrast is determined by the particular value chosen for the cross-correlation delay u. Consequently images with different resolution and different TZ contrast can be computed from one and the same set of raw data.
STOCHASTIC
NMR IMAGING
45
The cross-correlation can be viewed in terms of a transformation of the timedependent system response into three-dimensional x’-y?-z’ space. For practical applications of the method it will be worthwhile to search for a fast transformation algorithm in connection with deterministic or pseudo-random radiofrequency and field gradient modulation in order to replace the rather time-consuming crosscorrelation used for reconstruction of the spin density. An alternative route to digital data processing is the analog computation of the cross-correlation. It is possible, at least in principle, to build one hardware crosscorrelation channel per picture element. An array of such channels can be combined to form a two-dimensional window which would enable real-time imaging. The only time delay for image construction in this case would be the integration time of the cross-correlation. ACKNOWLEDGMENTS The first part of this work was done at the University of New Brunswick in Fredericton, Canada, with support of a fellowship by the DAAD (Deutscher Akademischer Austauschdienst) under the research committee of NATO. Continuation of the project was made possible by a research stipend of the DFG (Deutsche ForschungsgemeinschaB) at the University of Bayreuth, Germany. Stimulating discussions with Professor R. Kaiser are gratefully acknowledged. REFERENCES 1. (a) P. MANSFIELD AND P. G. MORRIS, “Advances in Magnetic Resonances,” Suppl. 2, “NMR Imaging in Biomedicine,” Academic Press, New York, 1982; (b) L. KAUFMAN, L. E. CROOKS, AND A. R. MARGULIS, eds., “Nuclear Magnetic Resonance Imaging in Medicine,” Igaku-Shoin, New York, 1981; (c) P. A. BOTTOMLEY, Rev. Sci. Instrum. 53, 1319 (1983). 2. J. S. BENDAT AND A. G. PIERSOL, “Random Data: Analysis and Measurement Procedures,” WileyInterscience, New York, 1971. 3. (a) G. HORLIK AND G. M. HIEFTIE, Contemp. Top. Anal. Clin. Chem. 3, 153 (1978); (b) G. M. HIEFTJE,Am. Lab. 13, 76 (1981). 4. (a) R. R. ERNST, J. Mugs. Reson. 3, 10 (1970); (b) R. KAISER, J. Mugs. Reson. 3, 28 (1970); (c) J. DADOK AND R. F. SPRECHER,J. Mugs. Reson. 13,243 (1974); (d) B. BLUMICH AND D. ZIESSOW, J. Chem. Phys. 78, 1059 (1983). 5. B. BLUMICH AND D. Zmssow, J. Mugn. Reson. 52,42, (1983). 6. D. I. HOULT, private communication. 7. L. F. FEINER AND P. R. LOCHER, Appl. Phys. 22,257 (1980). 8. G. A. KORN, “Random Process, Simulation and Measurement,” p. 83 If, McGraw-Hill, New York, 1966. 9. P. MANSFIELD AND I. L. PYKETT, J. Mugn. Reson. 29,355 (1978). 10. J. D. DE CERTAINES, L. BENOIST, F. DARCEL-MENAULT, M. CHATEL, AND A. M. BERNARD, Bull. Magn. Reson. 5, 267 (1983).