Iterative RF pulse refinement for magnetic ... - Semantic Scholar
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
41
Iterative RF Pulse Refinement for Magnetic Resonance Imaging Eliot T. Lebsack, Member, IEEE, and Steven M. Wright*, Member, IEEE
Abstract—Selective RF pulses are needed for many applications in magnetic resonance imaging (MRI). The waveform required to produce a desired excitation profile is, to first-order, its Fourier transform. This approximation is most valid for small tip angles and the quality and accuracy of such excitations decreases with increasing tip angle. Since large-tip-angle excitations are required in most types of imaging, a better synthesis technique is necessary. While a variety of analytical and numerical synthesis techniques based on solution of the Bloch equations are available, these techniques fail to consider the effect of the physical scanner hardware and are often accompanied by computational complexity. We present a technique for selective RF pulse refinement which uses real-time feedback techniques in lieu of a solution to the Bloch equations. Physical experiments are conducted to demonstrate the effectiveness of this algorithm and an extension to pulses of 90 is investigated. Index Terms—Feedback control, iterative refinement, magnetic resonance imaging, selective excitations.
I. INTRODUCTION
I
N MAGNETIC resonance imaging (MRI), shaped radio-frequency (RF) pulses are used to selectively encode spins within the sample. Development of such pulses has been a topic of interest to researchers for over 20 years [1] and many methods have been developed. Hoult showed that at small tip angles (30 or less), the characteristic Bloch equations are nearly linear in nature [2]. As a result, assuming linear response in the spectrometer transmitter chain, the Fourier transform (FT) of the desired magnetization profile is a reasonable choice for the RF waveform. Pulses derived in this manner can produce acceptable results at larger tip angles, but as the tip angle approaches and exceeds 90 , the Bloch equations become increasingly nonlinear, with significant distortion in both phase and magnitude. For tomographic imaging applications, such degradation tends to decrease the accuracy of the selected slice. If the candidate pulse is used in a multislice imaging sequence, this distortion may produce unwanted interactions between slices [3]. Recent interest in non-Fourier-based imaging techniques has increased the demand for accurate radio-frequency (RF) excitations, including pulses that produce wavelet excitations [4]. Manuscript received June 7, 2000; revised August 6, 2001 This work was supported by the Texas Higher Education Coordinating Board under Grant 005120315-1997. Asterisk indicates corresponding author. E. T. Lebsack is with the Department of Electrical Engineering, Texas A&M University, College Station, TX 77843-3128 USA. *S. M. Wright is with the Department of Electrical Engineering, M/S 3128, Texas A&M University, College Station, TX 77843-3128 USA (e-mail: [email protected]). Publisher Item Identifier S 0018-9294(02)00090-3.
Researchers have extensively investigated both direct and inverse approaches to RF pulse optimization. Hoult [2] and Mansfield et al. [5] were pioneers in the direct solution of the Bloch equations for shaped pulse design, whose work was verified numerically by Locher [6]. Techniques for analytical inversion of the Bloch equations were also proposed by Caprihan [7] and Silver et al. [8]. Caprihan concluded that numerical methods are the key to effective selective pulse synthesis. Conolly et al. [9] and Murdoch et al. [3] explored optimal control adaptations to RF pulse refinement. This work was reinforced by Ngo and Morris [10], who showed that the Bloch equations respond linearly to large excitations if they are treated as a superposition of perturbations of small tip angle. Shinnar et al. [11], [12] applied digital filtering techniques to obtain improved excitations and Pauly et al. [13] further developed this work to formalize the well-known Shinnar–Le Roux (SLR) method. The SLR technique has proven to be quite popular: Pickup and Popescu [14] extended the SLR technique to produce inversion pulses and Gelman and Wood [15] developed 90 and 180 Haar wavelet pulses as well. Another technique, proposed by Panych and Kyriakos, relies upon numerical Bloch simulation techniques to iteratively correct an existing RF waveform [16]. One limitation of these numerical techniques is that they fail to account for the nonlinear characteristics of the physical magnetic resonance (MR) hardware. The RF and gradient amplifiers are perhaps the most obvious nonlinear components in a typical system. On medical scanners, calibration schemes have been developed to counter such nonlinear behavior. Chan et al.developed a technique of RF pulse prewarping to counteract the effects of the RF amplifier itself [17]. However, there are some components in any system whose behavior cannot easily be calibrated out of the response, in particular the gradient coils and transmit/receive coils. This problem is especially significant on research MR scanners, where such components are often developed in-house and optimized for specific applications. This behavior is illustrated in Fig. 1(a), where a square profile is desired. Using the inverse FT of the square profile as the candidate RF pulse yields an excitation profile which is distorted in a nonlinear fashion by such factors, as shown in Fig. 1(b). Iterative pulse refinement techniques that compensate for such behavior have been proposed for other applications. Christy et al. devised a technique for recursively enhancing spectral features by copying the echo signal due to an RF excitation back into the RF digital-to-analog converter (DAC) buffer [18]. However, that paper does not address the general problem of slice selection refinement. In particular, the motivation for the present paper was that the authors required a method for generating simple, low angle RF pulses for use in wavelet imaging experiments on a noncommercial MR scanner.
0018–9294/02$17.00 © 2002 IEEE
42
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
Fig. 2. Magnetization profile visualization pulse sequence. The shaped RF t is transmitted into the sample in the presence of an encoding pulse , gradient G . The demodulated spin echo signal S t , occurring at t has a spectrum which is weighted by the spectral characteristics of the shaped RF pulse. This sequence is symmetric about the nonselective pulse and the extents of the k -space trajectory are identical on each side of this pulse. This is the sequence used to collect all magnetization profiles in this paper.
B ()
Fig. 1. Example of nonideal response of the system due to both system and Bloch equation nonlinearities. In (a), the ideal profile is specified in a scale consistent with what is to be expected from the system. By exciting the system with the inverse FT of this profile, the profile (b) is acquired, showing clearly the problems attributable to the transfer function of the physical system. The readout gradient strength and digitization dwell time are identical to those values chosen for encoding. The tip angle of the acquired profile is approximately 65 .
Previous researchers have explored feedback techniques using numerical Bloch equation simulation results from a computer to direct the optimization, as well as for controlling other aspects of MRI [19]. This paper presents a new RF pulse synthesis technique that relies upon the response of the physical MR system as part of a real-time feedback method to obtain improved slice-selective RF pulses. These pulses are realized without special consideration of the nonlinear nature of the Bloch equations at large tip angles. This paper also discusses several issues that can affect performance of the technique, including effects due to damping and scaling of the correction signals.
()
=2
Fig. 3. Complete block diagram of the iterative RF pulse refinement algorithm. The MR scanner block yields a magnetization profile z, which is the FT of the acquired spin echo signal. The difference profile z is scaled by the appropriate coefficient, producing the time-domain correction signal t . An optional window is applied, in conjunction with a scalar multiplier k . k is analogous to a damping coefficient in a dynamic system. This signal is added back to the input t . The circular nature of the algorithm is evident.
1M ( )
M ()
1B ( )
B ()
II. METHODS In order to implement the iterative refinement method, it is necessary to obtain an accurate representation of the magnetization profile, both magnitude and phase, at each step of the iteration. The pulse sequence used to both excite the profile and read it out is shown in Fig. 2, as suggested by Devoulon et al. [20]. , is transmitted in the presThe shaped RF waveform, . This creates a spatial transence of an encoding gradient along the axis which reverse magnetization profile . In the second half of the pulse sequence, sembles the FT of is following the 180 pulse, the demodulated echo signal collected in the presence of . The spectrum of is weighted . by the magnetization profile Iterative refinement techniques are implemented with at least one goal in mind. In the case of RF pulse synthesis, one seeks to and the desired magnetization minimize the error between . If these signals are scaled appropriately with profile respect to each other, the error or difference between these two signals is denoted by (1) In a feedback control system, error signals are often used to manipulate the input signal in order to achieve stability or other
desirable characteristics. In this case, the error signal is defined in the spatial domain, while the input signal exists in the time domain. As discussed previously, the Bloch equations are linear for small perturbations, which means that the inverse FT can from be used to obtain a time-domain correction signal . Addition of to from the previous iteration yields a corrected input signal, which is then used to excite the already relaxed spins. By iterating in this manner, an improved waveform can be obtained which compensates for some of the system nonlinearities, including those associated with the nonlinear behavior of the Bloch equations. No numerical solution of the Bloch equations is required; the physical MR system for successive iterations. A block diais used to provide gram of this iterative refinement technique is shown in Fig. 3. The optimization process continues until the algorithm converges upon the desired result. Various methods for determining convergence exist. In the past, researchers have proposed tracking the energy of the difference profile between iterations [9]. The difference energy at the th itera, is more useful if it is normalized with respect tion, to the initial difference energy (2)
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
43
Over successive iterations, the normalized absolute metric will ideally exist on the interval [0 1]. By quantifying the convergence of the algorithm in this manner, it is possible to modify the algorithm so that it can halt optimization once the normalized absolute metric falls below a predefined threshold. This threshold value is chosen with respect to the initial difference energy. The normalized absolute metric of convergence is useful for tracking the progress of the optimization algorithm when large are necessary. This is true when is changes to . However, if resemconsiderably different from with considerable accuracy, only small corrections bles may be required. Since the metric is normalized with respect , establishing a termination threshold that permits to such minor changes may be difficult. If the threshold is set prohibitively low, the algorithm may not terminate itself. Conversely, a large threshold may not produce an acceptable magnetization profile. An alternative metric of convergence may be established to address this shortcoming. By taking into account the change in between iterations, it is possible to define the energy of a relative metric of convergence. This metric describes the rel. ative change in the energy of For this algorithm, the metric of relative convergence is defined, for the th iteration, as (3) This metric may be monitored over the course of several iterations. For this metric, it is also possible to use a threshold to halt optimization progress. However, this threshold is defined with respect to the difference energy of the previous iteration, instead of the initial difference energy. This threshold is set to a value on the interval [0 1]; as the threshold approaches zero, the algorithm will tend to iterate indefinitely. The threshold should be set as small as possible for best results. Obviously, this metric is not self-starting; is undefined and is usually set to a value on the interval [0 1]. During development of this algorithm, several issues were considered. Perhaps the most obvious was the issue of scaling of the various signals used in the iteration procedure. Implementing an iterative design numerically is straightforward, as all scaling steps of the process are well known and repeatable. On the physical system, however, there are a variety of steps in which the signals are scaled. For example, the representais determined by the resolution of the DAC. For a tion of 12-bit DAC, the range is defined on the interval [ 2048, 2047]. might be digitized Upon demodulation, the acquired echo using a pair of 16-bit analog-to-digital converters (ADC), where individual values exist on the interval [ 32 768, 32 767]. The scaling depends on the settings of the various amplifiers and attenuators in the system. Additionally, the FT operator is implemented with different scaling factors on different machines. was While software could account for these factors, chosen as the reference scale for the spatial domain calculations. This selection makes tip angle selection and maintenance be considerations straightforward and it also requires that
scaled appropriately. When is formulated, it is assigned a maximum magnitude of one. On the first iteration of the alis adjusted so that correctly gorithm, the scale of describes the error in the magnetization profile. On subsequent remains unchanged. This arrangeiterations, the scale of ment allows one to select the target tip angle from the first iteration. is rescaled As shown in Fig. 3, the feedback signal is defined as before the inverse FT operation. The norm of the square root of its energy [21]
(4) is similarly defined. If is to be converted between scales, one may simply multiply it by the scalar (5) The time-domain correction signal is simply the inverse FT of this rescaled signal. This transformation denoted by the larger triangle in Fig. 2. A second issue which arose was the effect of discontinuities in the corrections applied to the RF signal. If the desired prohas sharp edges, subtracting can result in similar file discontinuities or spikes, which may result in excessive Gibbs’ . To prevent ringing in the time domain correction signal such nonphysical behavior, a windowing function could be apfollowing the inverse FT operation. Alternatively, plied to or could be windowed. Feedback techniques often include a gain adjustment on the correction signal. Adjustment of the magnitude of the feedback signal is often required to achieve a stable mode of operation. As shown in Fig. 3, the scalar multiplier may be used to maintain stability over successive iterations. As the algorithm converges upon an RF waveform which is nearly ideal, it is sometimes desirable to control the amplitude of the correction to enhance convergence. At small tip angles, where the behavior of the spins may be assumed to be linear, this might not be necessary. However, at large tip angles, perturbations to the input signal may push the algorithm into an unstable mode of operation. Proper adjustment of was found to provide an additional measure of stability to the algorithm, since the Bloch equations are virtually linear for small perturbations from existing magnetization states [10]. As the algorithm approaches convergence, the relative metric waveform is of convergence decreases as a near-optimal attained. This is precisely the behavior required for adaptively was found to provide good stability, adjusting . Setting even under conditions when the behavior of the spins is known to be nonlinear. III. RESULTS A series of experiments explored various aspects of the optimization algorithm. The spin-echo pulse sequence shown in Fig. 2 was used to examine the magnetization profile achieved
44
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
during each iteration. For convenience, this sequence was chosen such that the bandwidth of the excitation is identical to that of the acquired signal. This ensures that a 1:1 relationship exists between the spectrum of the RF waveform and the magnetization profile it excites, as detected by the spectrometer. The authors chose this direct relationship in an effort to simplify implementation and testing of the optimization process. Alternatively, one could employ the same encoding sequence parameters as those used in the intended application. The iterative technique was implemented in an Oxford 2-T horizontal, 31-cm bore superconducting magnet, with a TecMag Libra spectrometer and GE Accustar shielded gradients. The was shimmed using the supplied MacNMR package. static The host computer was a Power Macintosh 7100/66. Phase control of the RF pulse was achieved with a PTS-160 Frequency Synthesizer which was fitted with a digital phase rotation option. The phase rotation option was accessed with digital lines controlled directly from the pulse sequence. Phase modulation of the RF pulse waveform was achieved by directly controlling the phase of the carrier from the PTS. The phantom was a glass test tube measuring 2 23 cm, filled with a 0.5-mM CuSO solution. The phantom was oriented so that its longest dimension lies along the centerline of the RF birdcage and magnet. was defined with a length of The candidate RF waveform 1024 points. Originally, the authors devised the algorithm to obtain accurate excitations for wavelet-encoded MRI. Using 1024 provides sufficient precision for describing all of points for the excitations in our chosen multiresolution wavelet pulse set; the results reported here are derived from wavelet development efforts. The field of view was 6.43 cm, which was selected based and existing wavelet upon the measured homogeneity of imaging pulse sequences. The repetition time (TR) for all pulse sequences was 500 ms and the width of the nonselective inversion pulse was 122 s. When discussing tip angle, it is important to mention the technique used to establish the 90 reference power for a specific pulse. In this paper, the authors determine the reference point by repetitively exciting the phantom, while manually adjusting the applied RF power during the encoding phase until the acquired echo amplitude is maximized. This becomes the 90 reference point for a given RF envelope. For all experiments, the maximum magnitude of the initial profile was selected to be identical to the maximum magniproduced by the unoptimized candidate RF pulse. tude of This choice was based on the desire to maintain tip angle for a given RF pulse. Over time, the algorithm should be able to maintain the initial tip angle, as discussed previously. The first experimental investigation focuses upon the mechanics of the algorithm. In Figs. 4 and 5, the optimization technique is presented in a graphical format, with signals at each step of the process. In Fig. 4(a), the ideal “boxcar” profile is shown, with tip angle specified at roughly 55 . In Fig. 4(b), waveform, simply the FFT of the ideal profile, is the initial plotted with respect to the dynamic range of the DAC. Using the pulse sequence in Fig. 2, the acquired magnetization prousing this RF pulse is acquired as shown in Fig. 4(c). file The difference between the ideal profile and the desired profile, , follows in Fig. 4(d). The correction signal is ob-
Fig. 4. Step-by-step illustration of the optimization algorithm. In (a), the ideal profile is specified; the inverse FT of this profile is used to create (t) (b). This pulse is used to excite the MR system and the resulting (z ) is shown in (c). The resulting error profile is shown in (d).
M
B
tained from the inverse FT of and is plotted in Fig. 5(a). 1.0 and a Hanning In this case, the feedback gain scalar pulse, window was applied. Fig. 5(b) shows the corrected which completes the first iteration. Fig. 5(c) shows the acquired waveform. After magnetization profile due to the updated four additional iterations, the magnetization profile produced by the improved RF waveform is provided in Fig. 5(d). In the next experiment, the optimization was carried out to 20 iterations for tip angles of approximately 30 and 65 . In an effort to study the effect of upon convergence and stability, the 1.0, 0.6, and . The correction algorithm was tested with signal was not windowed. The difference energy is plotted for each of these cases as shown in Fig. 6. In Fig. 6(a), the result for a tip angle of approximately 30 is plotted. In this case, the Bloch equations should be relatively linear, as shown by Hoult [2], and poor RF pulse behavior is likely due to system imperfections. As can be seen, 1.0 will yield the fastest convergence in this case. choosing will produce similarly quick results. However, setting If the tip angle is increased to approximately 65 , as shown in Fig. 6(b), the algorithm diverges after six iterations when 1.0. Reducing the value of or adaptively adjusting it as one nears convergence has the effect of damping on the correction; improved stability is realized as a result. Additional perspective is provided in Figs. 7 and 8. In Fig. 7, 1.0. This correthe selected tip angle is roughly 30 , with
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
45
=
Fig. 7. Graphical progression of the optimization process (k 1.0) for a tip angle of approximately 30 . The unoptimized RF sinc pulse (a) produces the profile (b), which exhibits substantial rounding of the features when compared to the expected boxcar shape. After ten iterations, the optimized RF pulse (c) produces a much-improved magnetization profile (d). Fig. 5. Continued step-by-step illustration of the optimization algorithm. The inverse FT of the difference profile shown in Fig. 4 produces the time domain correction signal t (a). This signal is added to the original RF pulse in Fig. 3, yielding a corrected signal (b). The profile produced by the corrected is shown in (c). After five iterations, the profile shows a marked improvement in shape (d).
1B ( )
B
B
Fig. 6. Convergence trends observed for optimizations performed at approximately 30 and 65 . The normalized energy of the difference profile is plotted for each of these studies. At 30 (a), the optimization converges for all selected values of k . If the tip angle is increased to 65 (b), the results are mixed. If k 1.0, the algorithm proceeds through a local minimum without settling within it. This illustrates the nonlinear behavior of the spin, t is scaled to be roughly 60% as described by the Bloch equation. If of its full value before addition to t , the optimization is well-behaved. If k is set dynamically, as the square root of the relative convergence metric d, stability is also apparent.
=
1B ( ) B ()
sponds to the associated curve in Fig. 7(a). In Fig. 7(b), the magnetization profile due to the uncorrected RF waveform is plotted. There is substantial smoothing due to the system transfer function. In Fig. 7(c), the candidate RF pulse after ten iterations is shown. The sharp features of the desired boxcar shape are obvious. After ten additional iterations, there is minimal change, as shown in Fig. 7(d). 0.6,1.0, and Fig. 8 illustrates optimization results for at a tip angle of approximately 65 . In Fig. 8(b), choosing 0.6 produces a well-defined square profile with negligible peaks in the stopband. If is set to 1.0, as shown in Fig. 8(c), the resulting profile is distorted with substantial ringing in the stopband. The magnitude of this error is also exhibited in the correyields the result sponding curve in Fig. 6(b). Choosing shown in Fig. 6(d), which is nearly identical to that obtained for 0.6, with slightly larger harmonic ripples in the stopband. An application for this algorithm is refinement of excitations of 90 . The unique nature of the apparent magnetization around 90 makes such optimization difficult. Since the Bloch equations are linear for small changes to existing tip angle, a perturbational approach was adopted. This method requires setting . A flowchart detailing this approach is shown in Fig. 9. To begin, a candidate RF pulse is optimized until it has converged sufficiently for a tip angle below 90 . The power is then increased to 90 and the pulse is optimized further until convergence is attained yet again. Convergence in both cases was defined when dropped below 10%. The results from this optimization are shown in Fig. 10. For this experiment, the intermediate tip angle was approximately 60 ; this tip angle was chosen based on the near-linear nature
46
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
Fig. 8. Graphical comparison of the effects of k on the optimization for a tip angle of approximately 65 . The profile in (a) results from the uncorrected RF waveform. If k 0.6, as shown in (b), the resulting profile is very square. If k is chosen to be 1.0 throughout the optimization instead, the resulting profile exhibits highly nonlinear behavior. By choosing the feedback scalar factor d (d), it is possible to produce results identical to (b), without to be k the additional problem of selecting an optimal value for k . Reduction of the feedback scalar k provided the necessary stability for the algorithm to continue to an acceptable solution.
=
=
p
Fig. 10. 90 optimization for a threshold of 10%. The sub-90 ideal profile is shown in (a), followed by the magnetization profile due to the uncorrected RF pulse (b). After seven iterations, the algorithm terminated itself, producing the profile shown in (c). The power is then set to 90 , where subsequent optimization produces the profile (d).
havior was numerically modeled in a Bloch simulator by truncating the lowest five bits of the applied RF pulse in the feedback algorithm, which resulted in similar periodic ripples. In the case of the sinc pulse, the algorithm attempts to compensate for this behavior by increasing the amplitude of the sidelobes, in the same way as the example shown in Fig. 11. As the optimization progresses, the algorithm attempts to overcome the noise threshold and is largely successful, as shown in Figs. 7 and 8. Since this algorithm was developed originally to produce improved wavelet excitations, an illustration of a partial Haar wavelet family is shown in Fig. 12. After three iterations, it is possible to obtain reasonable scaling and wavelet functions. The ripple behavior discussed previously appears here as well.
IV. DISCUSSION Fig. 9. Flowchart detail of the 90 optimization process. The tip angle is first set to a value below 90 and optimization is conducted until d is below a preselected threshold. The power is then increased to a point where the excitation is nearest to 90 and optimization continues for a few more iterations until the threshold is crossed again. The result is a pulse which is correct for a 90 tip angle. The threshold values are identical for each phase of the process. In Part 1, d is set to 1.0 since it is otherwise undefined. In Part 2, d is set to be 0.5.
of the Bloch equations, as well as the proximity in signal amplitude associated with 90 . In Fig. 10(a), the ideal magnetization profile for the sub-90 case is displayed, followed by the magnetization profile due to the uncorrected RF waveform. After seven iterations, fell below 10%, yielding the profile in Fig. 10(c). At this point, the power was increased to 90 and two further optimization iterations yielded the result shown in Fig. 10(d). The spurious harmonic ripples observed in Figs. 7, 8, and 11 were attributed to low dynamic range in the RF DACs. This be-
From the signals shown in Figs. 4 and 5, the simplicity of the algorithm is evident. Simple signal addition and subtraction of signals combined with scaling techniques produces an algorithm that is straightforward. In Fig. 6, the algorithm exhibits significant divergent behavior for tip angles around 65 when is equal to 1.0 [Fig. 6(b)]. which is added This is attributable to the correction signal . Since the Bloch equations are linear for small perturbato tions of tip angle, better waveforms can be obtained by reducing the scale of the correction signals as the tip angle increases. Selecting a smaller, fixed value for is an obvious modification, produces as shown in Fig. 8(b), but allowing to vary as equally effective results, as illustrated by Fig. 8(d). This is an excellent example of linearity in the Bloch equations; as the tip produce angle becomes large, sufficiently small changes to . Such small changes are correspondingly linear changes in
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
47
Fig. 11. Illustration of the threshold compensation behavior of the optimization algorithm. In (a), the unoptimized sinc pulse is shown alongside the resulting slice profile (b). Because the lowest bits of the RF DAC contribute virtually nothing to the output pulse, the optimization algorithm overcomes this loss by emphasizing the amplitude of the sidelobes (c). The resulting slice profile after five iterations of optimization shows the small periodic ripples (d).
necessary for optimization of RF waveforms for tip angles of 90 or 180 . The difference in the rate of convergence due to different is noticeable as well. This is to be expected, since large changes will produce correspondingly dramatic changes to . to The final difference energy was observed to be smallest for both tip angles when was set to 0.6. Clearly, the selection of will influence the quality and extent of convergence in this case and this choice needs to be examined further. This algorithm may be useful in practical applications, such , the as the optimization of excitations of 90 . By setting as it algorithm can make increasingly smaller corrections to approaches convergence. As shown in Fig. 10, the algorithm is capable of producing reasonable excitations of approximately 90 when used in a perturbational scheme. The algorithm has been shown to correct phase errors in the acquired profile. Correction of phase is especially important in non-Fourier-based imaging techniques, such as wavelet-encoded MR. In wavelet-encoded MR, the excitation profiles required by the reconstruction algorithm are phase-sensitive in nature and may require correction [16]. Fig. 12 provides evidence of the utility of this algorithm for such excitations. A subtle trend in Figs. 7 and 8 is the maintenance of tip angle. This is achieved by specifying the ideal profile in terms of the normalized magnetization. After 20 iterations, the peak values of the magnetization profile are nearly identical to the corresponding values of the ideal profile. The reason for tip angle consistency is not as obvious. As a corrective scheme, the algorithm should be allowed to alter
p
Fig. 12. Haar Wavelet optimization example. In all cases, the tip angle is d. After three iterations, the Haar wavelets and approximately 30 , k scaling functions are reasonably reproduced.
=
the magnitude of specific points in . Continuous addition of may result in occasions where the the correction signal to dynamic range of the DAC hardware is exceeded. Therefore, it . This is achieved by is necessary to establish a margin for so that its maximum value is somewhat smaller selecting than the dynamic range of the DAC. For the results shown in Figs. 6–8, the initial magnitude was selected to be 70% of the maximum. Examination of several of the figures reveals a recurring phenomenon. Harmonic ripples would evolve and, in most cases, become negligible after several optimization iterations. As discussed previously, this behavior was determined to originate in the RF transmitter signal path. The algorithm is able to partially compensate for this behavior. It is clear that the algorithm is capable of refining RF waveforms where the small-tip-angle assumptions are not applicable. With some modification, the algorithm can produce improved excitations at 90 . A history-based approach where
48
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
one considers past acquired and difference profiles might make the algorithm more robust for such excitations as well as inversion pulses. However, this will likely require changes in the algorithm to accommodate potentially confusing behavior at the slice edges. V. CONCLUSION The FT properties of the MR spin system provide the basis for rough RF pulse formulations. Methods have been proposed to further refine the RF pulses obtained in this manner, but only two of these methods have considered effects due to the components of the overall imaging system [17], [18]. The physical transfer function, which is affected by magnetic field inhomogeneities and the spectrometer itself, should be taken into account when constructing RF pulses used to excite precise regions within the sample. In this paper, we present a method of RF pulse refinement which iteratively compensates for the apparent physical transfer function of the system. The technique relies on the formulation of an error profile, which is in turn used to correct the excitation RF pulse. We also discuss methods of filtering and scaling which will provide acceptable results when computing correction signals based upon error signals. The optimization tests presented demonstrate the potential effectiveness of such an algorithm. Simple and complex shapes were reproduced with little effort on an experimental MR system and criteria for convergence was introduced. A damping coefficient was introduced which produced favorable results for larger tip angles and dynamic adjustment of this feedback gain was shown to produce beneficial results. In addition, the authors have shown that the method is directly applicable to the problem of realizing acceptable wavelet pulses. The algorithm may be used on host microcomputers and the effectiveness and efficiency of this algorithm lends itself to implementation on a fast imaging workstation for real-time correction of RF pulses. This algorithm provides a simple method of generating highly specific excitation profiles. Unlike analytical and numerical approaches, the iterative feedback approach presented here does not require compensation for nonlinearities in the RF transmitter. The method should be useful in applications where a series of different shaped RF profiles are needed, such as nonFourier based imaging. REFERENCES [1] A. N. Garroway, P. K. Grannell, and P. Mansfield, “Image formation in NMR by a selective irradiative process,” J. Phys. C, vol. 7, pp. L457–L462, 1974. [2] D. I. Hoult, “The solution of the Bloch equations in the presence of a varying field-an approach to selective pulse analysis,” J. Magn. Reson., vol. 35, pp. 69–86, 1979. [3] J. B. Murdoch, A. H. Lent, and M. K. Kritzer, “Computer-optimized narrowband pulses for multislice imaging,” J. Magn. Reson., vol. 74, pp. 226–263, 1987. [4] J. B. Weaver, Y. Xu, D. M. Healy, and J. R. Driscoll, “Wavelet-encoded MR imaging,” Magn. Reson. Med., vol. 24, pp. 275–287, 1992. [5] P. Mansfield, A. A. Maudsley, P. G. Morris, and I. L. Pykett, “Selective pulses in NMR imaging: A reply to criticism,” J. Magn. Reson., vol. 33, pp. 261–274, 1979. [6] P. R. Locher, “Computer simulation of selective excitation in n.m.r. imaging,” Phil. Trans. R. Soc. Lond., vol. B.289, pp. 537–542, 1980.
B
[7] A. Caprihan, “Effect of amplitude modulation on selective excitation in NMR imaging,” IEEE Trans. Med. Imag., vol. MI-4, pp. 169–175, 1983. [8] M. S. Silver, R. I. Joseph, and D. I. Hoult, “Highly selective =2 and pulse generation,” J. Magn. Reson., vol. 59, pp. 347–351, 1984. [9] S. Conolly, D. Nishimura, and A. Macovski, “Optimal control solutions to the magnetic resonance selective excitation problem,” IEEE Trans. Med. Imag., vol. MI-2, pp. 106–115, 1986. [10] J. T. Ngo and P. G. Morris, “General solution to the NMR excitation problem for noninteracting spins,” Magn. Reson. Med., vol. 5, pp. 217–237, 1987. [11] M. Shinnar, L. Bolinger, and J. S. Leigh, “The use of finite impulse response filters in pulse design,” Magn. Reson. Med., vol. 12, pp. 81–87, 1989. , “The synthesis of soft pulses with a specified frequency response,” [12] Magn. Reson. Med., vol. 12, pp. 88–92, 1989. [13] J. Pauly, P. Le Roux, D. Nishimura, and A. Macovski, “Parameter relations for the Shinnar–Le Roux selective excitation pulse design algorithm,” IEEE Trans. Med. Imag., vol. 10, pp. 53–65, Mar. 1991. [14] S. Pickup and M. Popescu, “Efficient design of pulses with trapezoidal magnitude and linear phase response profiles,” Magn. Reson. Med., vol. 38, pp. 137–145, 1997. [15] N. Gelman and M. L. Wood, “Wavelet encoding for improved SNR and retrospective slice adjustment,” Magn. Reson. Med., vol. 39, pp. 383–391, 1998. [16] W. E. Kyriakos and L. P. Panych, Proceedings of the 5th Scientific Meeting of the ISMRM, Vancouver, B.C., Canada, April 1997. [17] F. Chan, J. Pauly, and A. Macovski, “Effects of RF amplifier distortion on selective excitation and their correction by prewarping,” Magn. Reson. Med., vol. 23, pp. 224–238, 1992. [18] P. S. Christy, R. C. Grimm, J. F. Greenleaf, and R. L. Ehman, “Recursive RF excitation,” Magn. Reson. Med., vol. 26, pp. 361–367, 1992. [19] J. L. Schiano, R. L. Magin, and S. M. Wright, “Feedback control of the nuclear magnetization state,” IEEE Trans. Med. Imag., vol. 10, pp. 138–147, Mar. 1991. [20] P. Devoulon, L. Emsley, P. Weber, R. Meuli, M. Decorps, and G. Bodenhausen, “Methods for reconstructing phase sensitive slice profiles in magnetic resonance imaging,” Magn. Reson. Med., vol. 31, pp. 178–83, 1994. [21] L. E. Franks, Signal Theory. Stroudsberg, PA: Dowden & Culver, 1981.
Eliot T. Lebsack (S’97–M’99) was born in Williston, ND, in 1972. He received the B.S. degree in aerospace engineering and the M.S. degree in electrical engineering from Texas A&M University, College Station, in 1995 and 1999, respectively. He was with Lernout & Hauspie Speech Products in Burlington, MA. Since 2000, he has been working in the area of computational electromagnetics at The MITRE Corporation, Bedford, MA. His interests include the application of high-performance computing techniques to electromagnetics (radio-frequency dosimetry, radar cross section (RCS) prediction).
Steven M. Wright (S’78–M81) received the B.S, M.S., and Ph.D. degrees in electrical engineering from the University of Illinois at Urbana-Champaign in 1980, 1981, and 1984, respectively. From 1984 to 1988, he was an Engineer/Scientist for Magnetic Resonance Imaging at Saint Francis Medical Center, Peoria, IL. He joined the faculty at Texas A&M University, College Station, in 1988 and is currently a Professor of Electrical Engineering and Director of the Magnetic Resonance Systems Laboratory. During the summer and fall of 2000, he was a Visiting Professor at the University of Texas M.D. Anderson Cancer Center. His research interests are in the development of instrumentation and techniques for magnetic resonance imaging and spectroscopy and in applied computational electromagnetics. He has authored or co-authored over 130 journal and conference publications and three book chapters Prof. Wright was awarded the Lockheed Fort Worth Company Award for Excellence in Engineering Education in 1993, the James Stone Faculty Fellow Award in 1999 and a TEES Fellow award in 1999.
41
Iterative RF Pulse Refinement for Magnetic Resonance Imaging Eliot T. Lebsack, Member, IEEE, and Steven M. Wright*, Member, IEEE
Abstract—Selective RF pulses are needed for many applications in magnetic resonance imaging (MRI). The waveform required to produce a desired excitation profile is, to first-order, its Fourier transform. This approximation is most valid for small tip angles and the quality and accuracy of such excitations decreases with increasing tip angle. Since large-tip-angle excitations are required in most types of imaging, a better synthesis technique is necessary. While a variety of analytical and numerical synthesis techniques based on solution of the Bloch equations are available, these techniques fail to consider the effect of the physical scanner hardware and are often accompanied by computational complexity. We present a technique for selective RF pulse refinement which uses real-time feedback techniques in lieu of a solution to the Bloch equations. Physical experiments are conducted to demonstrate the effectiveness of this algorithm and an extension to pulses of 90 is investigated. Index Terms—Feedback control, iterative refinement, magnetic resonance imaging, selective excitations.
I. INTRODUCTION
I
N MAGNETIC resonance imaging (MRI), shaped radio-frequency (RF) pulses are used to selectively encode spins within the sample. Development of such pulses has been a topic of interest to researchers for over 20 years [1] and many methods have been developed. Hoult showed that at small tip angles (30 or less), the characteristic Bloch equations are nearly linear in nature [2]. As a result, assuming linear response in the spectrometer transmitter chain, the Fourier transform (FT) of the desired magnetization profile is a reasonable choice for the RF waveform. Pulses derived in this manner can produce acceptable results at larger tip angles, but as the tip angle approaches and exceeds 90 , the Bloch equations become increasingly nonlinear, with significant distortion in both phase and magnitude. For tomographic imaging applications, such degradation tends to decrease the accuracy of the selected slice. If the candidate pulse is used in a multislice imaging sequence, this distortion may produce unwanted interactions between slices [3]. Recent interest in non-Fourier-based imaging techniques has increased the demand for accurate radio-frequency (RF) excitations, including pulses that produce wavelet excitations [4]. Manuscript received June 7, 2000; revised August 6, 2001 This work was supported by the Texas Higher Education Coordinating Board under Grant 005120315-1997. Asterisk indicates corresponding author. E. T. Lebsack is with the Department of Electrical Engineering, Texas A&M University, College Station, TX 77843-3128 USA. *S. M. Wright is with the Department of Electrical Engineering, M/S 3128, Texas A&M University, College Station, TX 77843-3128 USA (e-mail: [email protected]). Publisher Item Identifier S 0018-9294(02)00090-3.
Researchers have extensively investigated both direct and inverse approaches to RF pulse optimization. Hoult [2] and Mansfield et al. [5] were pioneers in the direct solution of the Bloch equations for shaped pulse design, whose work was verified numerically by Locher [6]. Techniques for analytical inversion of the Bloch equations were also proposed by Caprihan [7] and Silver et al. [8]. Caprihan concluded that numerical methods are the key to effective selective pulse synthesis. Conolly et al. [9] and Murdoch et al. [3] explored optimal control adaptations to RF pulse refinement. This work was reinforced by Ngo and Morris [10], who showed that the Bloch equations respond linearly to large excitations if they are treated as a superposition of perturbations of small tip angle. Shinnar et al. [11], [12] applied digital filtering techniques to obtain improved excitations and Pauly et al. [13] further developed this work to formalize the well-known Shinnar–Le Roux (SLR) method. The SLR technique has proven to be quite popular: Pickup and Popescu [14] extended the SLR technique to produce inversion pulses and Gelman and Wood [15] developed 90 and 180 Haar wavelet pulses as well. Another technique, proposed by Panych and Kyriakos, relies upon numerical Bloch simulation techniques to iteratively correct an existing RF waveform [16]. One limitation of these numerical techniques is that they fail to account for the nonlinear characteristics of the physical magnetic resonance (MR) hardware. The RF and gradient amplifiers are perhaps the most obvious nonlinear components in a typical system. On medical scanners, calibration schemes have been developed to counter such nonlinear behavior. Chan et al.developed a technique of RF pulse prewarping to counteract the effects of the RF amplifier itself [17]. However, there are some components in any system whose behavior cannot easily be calibrated out of the response, in particular the gradient coils and transmit/receive coils. This problem is especially significant on research MR scanners, where such components are often developed in-house and optimized for specific applications. This behavior is illustrated in Fig. 1(a), where a square profile is desired. Using the inverse FT of the square profile as the candidate RF pulse yields an excitation profile which is distorted in a nonlinear fashion by such factors, as shown in Fig. 1(b). Iterative pulse refinement techniques that compensate for such behavior have been proposed for other applications. Christy et al. devised a technique for recursively enhancing spectral features by copying the echo signal due to an RF excitation back into the RF digital-to-analog converter (DAC) buffer [18]. However, that paper does not address the general problem of slice selection refinement. In particular, the motivation for the present paper was that the authors required a method for generating simple, low angle RF pulses for use in wavelet imaging experiments on a noncommercial MR scanner.
0018–9294/02$17.00 © 2002 IEEE
42
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
Fig. 2. Magnetization profile visualization pulse sequence. The shaped RF t is transmitted into the sample in the presence of an encoding pulse , gradient G . The demodulated spin echo signal S t , occurring at t has a spectrum which is weighted by the spectral characteristics of the shaped RF pulse. This sequence is symmetric about the nonselective pulse and the extents of the k -space trajectory are identical on each side of this pulse. This is the sequence used to collect all magnetization profiles in this paper.
B ()
Fig. 1. Example of nonideal response of the system due to both system and Bloch equation nonlinearities. In (a), the ideal profile is specified in a scale consistent with what is to be expected from the system. By exciting the system with the inverse FT of this profile, the profile (b) is acquired, showing clearly the problems attributable to the transfer function of the physical system. The readout gradient strength and digitization dwell time are identical to those values chosen for encoding. The tip angle of the acquired profile is approximately 65 .
Previous researchers have explored feedback techniques using numerical Bloch equation simulation results from a computer to direct the optimization, as well as for controlling other aspects of MRI [19]. This paper presents a new RF pulse synthesis technique that relies upon the response of the physical MR system as part of a real-time feedback method to obtain improved slice-selective RF pulses. These pulses are realized without special consideration of the nonlinear nature of the Bloch equations at large tip angles. This paper also discusses several issues that can affect performance of the technique, including effects due to damping and scaling of the correction signals.
()
=2
Fig. 3. Complete block diagram of the iterative RF pulse refinement algorithm. The MR scanner block yields a magnetization profile z, which is the FT of the acquired spin echo signal. The difference profile z is scaled by the appropriate coefficient, producing the time-domain correction signal t . An optional window is applied, in conjunction with a scalar multiplier k . k is analogous to a damping coefficient in a dynamic system. This signal is added back to the input t . The circular nature of the algorithm is evident.
1M ( )
M ()
1B ( )
B ()
II. METHODS In order to implement the iterative refinement method, it is necessary to obtain an accurate representation of the magnetization profile, both magnitude and phase, at each step of the iteration. The pulse sequence used to both excite the profile and read it out is shown in Fig. 2, as suggested by Devoulon et al. [20]. , is transmitted in the presThe shaped RF waveform, . This creates a spatial transence of an encoding gradient along the axis which reverse magnetization profile . In the second half of the pulse sequence, sembles the FT of is following the 180 pulse, the demodulated echo signal collected in the presence of . The spectrum of is weighted . by the magnetization profile Iterative refinement techniques are implemented with at least one goal in mind. In the case of RF pulse synthesis, one seeks to and the desired magnetization minimize the error between . If these signals are scaled appropriately with profile respect to each other, the error or difference between these two signals is denoted by (1) In a feedback control system, error signals are often used to manipulate the input signal in order to achieve stability or other
desirable characteristics. In this case, the error signal is defined in the spatial domain, while the input signal exists in the time domain. As discussed previously, the Bloch equations are linear for small perturbations, which means that the inverse FT can from be used to obtain a time-domain correction signal . Addition of to from the previous iteration yields a corrected input signal, which is then used to excite the already relaxed spins. By iterating in this manner, an improved waveform can be obtained which compensates for some of the system nonlinearities, including those associated with the nonlinear behavior of the Bloch equations. No numerical solution of the Bloch equations is required; the physical MR system for successive iterations. A block diais used to provide gram of this iterative refinement technique is shown in Fig. 3. The optimization process continues until the algorithm converges upon the desired result. Various methods for determining convergence exist. In the past, researchers have proposed tracking the energy of the difference profile between iterations [9]. The difference energy at the th itera, is more useful if it is normalized with respect tion, to the initial difference energy (2)
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
43
Over successive iterations, the normalized absolute metric will ideally exist on the interval [0 1]. By quantifying the convergence of the algorithm in this manner, it is possible to modify the algorithm so that it can halt optimization once the normalized absolute metric falls below a predefined threshold. This threshold value is chosen with respect to the initial difference energy. The normalized absolute metric of convergence is useful for tracking the progress of the optimization algorithm when large are necessary. This is true when is changes to . However, if resemconsiderably different from with considerable accuracy, only small corrections bles may be required. Since the metric is normalized with respect , establishing a termination threshold that permits to such minor changes may be difficult. If the threshold is set prohibitively low, the algorithm may not terminate itself. Conversely, a large threshold may not produce an acceptable magnetization profile. An alternative metric of convergence may be established to address this shortcoming. By taking into account the change in between iterations, it is possible to define the energy of a relative metric of convergence. This metric describes the rel. ative change in the energy of For this algorithm, the metric of relative convergence is defined, for the th iteration, as (3) This metric may be monitored over the course of several iterations. For this metric, it is also possible to use a threshold to halt optimization progress. However, this threshold is defined with respect to the difference energy of the previous iteration, instead of the initial difference energy. This threshold is set to a value on the interval [0 1]; as the threshold approaches zero, the algorithm will tend to iterate indefinitely. The threshold should be set as small as possible for best results. Obviously, this metric is not self-starting; is undefined and is usually set to a value on the interval [0 1]. During development of this algorithm, several issues were considered. Perhaps the most obvious was the issue of scaling of the various signals used in the iteration procedure. Implementing an iterative design numerically is straightforward, as all scaling steps of the process are well known and repeatable. On the physical system, however, there are a variety of steps in which the signals are scaled. For example, the representais determined by the resolution of the DAC. For a tion of 12-bit DAC, the range is defined on the interval [ 2048, 2047]. might be digitized Upon demodulation, the acquired echo using a pair of 16-bit analog-to-digital converters (ADC), where individual values exist on the interval [ 32 768, 32 767]. The scaling depends on the settings of the various amplifiers and attenuators in the system. Additionally, the FT operator is implemented with different scaling factors on different machines. was While software could account for these factors, chosen as the reference scale for the spatial domain calculations. This selection makes tip angle selection and maintenance be considerations straightforward and it also requires that
scaled appropriately. When is formulated, it is assigned a maximum magnitude of one. On the first iteration of the alis adjusted so that correctly gorithm, the scale of describes the error in the magnetization profile. On subsequent remains unchanged. This arrangeiterations, the scale of ment allows one to select the target tip angle from the first iteration. is rescaled As shown in Fig. 3, the feedback signal is defined as before the inverse FT operation. The norm of the square root of its energy [21]
(4) is similarly defined. If is to be converted between scales, one may simply multiply it by the scalar (5) The time-domain correction signal is simply the inverse FT of this rescaled signal. This transformation denoted by the larger triangle in Fig. 2. A second issue which arose was the effect of discontinuities in the corrections applied to the RF signal. If the desired prohas sharp edges, subtracting can result in similar file discontinuities or spikes, which may result in excessive Gibbs’ . To prevent ringing in the time domain correction signal such nonphysical behavior, a windowing function could be apfollowing the inverse FT operation. Alternatively, plied to or could be windowed. Feedback techniques often include a gain adjustment on the correction signal. Adjustment of the magnitude of the feedback signal is often required to achieve a stable mode of operation. As shown in Fig. 3, the scalar multiplier may be used to maintain stability over successive iterations. As the algorithm converges upon an RF waveform which is nearly ideal, it is sometimes desirable to control the amplitude of the correction to enhance convergence. At small tip angles, where the behavior of the spins may be assumed to be linear, this might not be necessary. However, at large tip angles, perturbations to the input signal may push the algorithm into an unstable mode of operation. Proper adjustment of was found to provide an additional measure of stability to the algorithm, since the Bloch equations are virtually linear for small perturbations from existing magnetization states [10]. As the algorithm approaches convergence, the relative metric waveform is of convergence decreases as a near-optimal attained. This is precisely the behavior required for adaptively was found to provide good stability, adjusting . Setting even under conditions when the behavior of the spins is known to be nonlinear. III. RESULTS A series of experiments explored various aspects of the optimization algorithm. The spin-echo pulse sequence shown in Fig. 2 was used to examine the magnetization profile achieved
44
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
during each iteration. For convenience, this sequence was chosen such that the bandwidth of the excitation is identical to that of the acquired signal. This ensures that a 1:1 relationship exists between the spectrum of the RF waveform and the magnetization profile it excites, as detected by the spectrometer. The authors chose this direct relationship in an effort to simplify implementation and testing of the optimization process. Alternatively, one could employ the same encoding sequence parameters as those used in the intended application. The iterative technique was implemented in an Oxford 2-T horizontal, 31-cm bore superconducting magnet, with a TecMag Libra spectrometer and GE Accustar shielded gradients. The was shimmed using the supplied MacNMR package. static The host computer was a Power Macintosh 7100/66. Phase control of the RF pulse was achieved with a PTS-160 Frequency Synthesizer which was fitted with a digital phase rotation option. The phase rotation option was accessed with digital lines controlled directly from the pulse sequence. Phase modulation of the RF pulse waveform was achieved by directly controlling the phase of the carrier from the PTS. The phantom was a glass test tube measuring 2 23 cm, filled with a 0.5-mM CuSO solution. The phantom was oriented so that its longest dimension lies along the centerline of the RF birdcage and magnet. was defined with a length of The candidate RF waveform 1024 points. Originally, the authors devised the algorithm to obtain accurate excitations for wavelet-encoded MRI. Using 1024 provides sufficient precision for describing all of points for the excitations in our chosen multiresolution wavelet pulse set; the results reported here are derived from wavelet development efforts. The field of view was 6.43 cm, which was selected based and existing wavelet upon the measured homogeneity of imaging pulse sequences. The repetition time (TR) for all pulse sequences was 500 ms and the width of the nonselective inversion pulse was 122 s. When discussing tip angle, it is important to mention the technique used to establish the 90 reference power for a specific pulse. In this paper, the authors determine the reference point by repetitively exciting the phantom, while manually adjusting the applied RF power during the encoding phase until the acquired echo amplitude is maximized. This becomes the 90 reference point for a given RF envelope. For all experiments, the maximum magnitude of the initial profile was selected to be identical to the maximum magniproduced by the unoptimized candidate RF pulse. tude of This choice was based on the desire to maintain tip angle for a given RF pulse. Over time, the algorithm should be able to maintain the initial tip angle, as discussed previously. The first experimental investigation focuses upon the mechanics of the algorithm. In Figs. 4 and 5, the optimization technique is presented in a graphical format, with signals at each step of the process. In Fig. 4(a), the ideal “boxcar” profile is shown, with tip angle specified at roughly 55 . In Fig. 4(b), waveform, simply the FFT of the ideal profile, is the initial plotted with respect to the dynamic range of the DAC. Using the pulse sequence in Fig. 2, the acquired magnetization prousing this RF pulse is acquired as shown in Fig. 4(c). file The difference between the ideal profile and the desired profile, , follows in Fig. 4(d). The correction signal is ob-
Fig. 4. Step-by-step illustration of the optimization algorithm. In (a), the ideal profile is specified; the inverse FT of this profile is used to create (t) (b). This pulse is used to excite the MR system and the resulting (z ) is shown in (c). The resulting error profile is shown in (d).
M
B
tained from the inverse FT of and is plotted in Fig. 5(a). 1.0 and a Hanning In this case, the feedback gain scalar pulse, window was applied. Fig. 5(b) shows the corrected which completes the first iteration. Fig. 5(c) shows the acquired waveform. After magnetization profile due to the updated four additional iterations, the magnetization profile produced by the improved RF waveform is provided in Fig. 5(d). In the next experiment, the optimization was carried out to 20 iterations for tip angles of approximately 30 and 65 . In an effort to study the effect of upon convergence and stability, the 1.0, 0.6, and . The correction algorithm was tested with signal was not windowed. The difference energy is plotted for each of these cases as shown in Fig. 6. In Fig. 6(a), the result for a tip angle of approximately 30 is plotted. In this case, the Bloch equations should be relatively linear, as shown by Hoult [2], and poor RF pulse behavior is likely due to system imperfections. As can be seen, 1.0 will yield the fastest convergence in this case. choosing will produce similarly quick results. However, setting If the tip angle is increased to approximately 65 , as shown in Fig. 6(b), the algorithm diverges after six iterations when 1.0. Reducing the value of or adaptively adjusting it as one nears convergence has the effect of damping on the correction; improved stability is realized as a result. Additional perspective is provided in Figs. 7 and 8. In Fig. 7, 1.0. This correthe selected tip angle is roughly 30 , with
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
45
=
Fig. 7. Graphical progression of the optimization process (k 1.0) for a tip angle of approximately 30 . The unoptimized RF sinc pulse (a) produces the profile (b), which exhibits substantial rounding of the features when compared to the expected boxcar shape. After ten iterations, the optimized RF pulse (c) produces a much-improved magnetization profile (d). Fig. 5. Continued step-by-step illustration of the optimization algorithm. The inverse FT of the difference profile shown in Fig. 4 produces the time domain correction signal t (a). This signal is added to the original RF pulse in Fig. 3, yielding a corrected signal (b). The profile produced by the corrected is shown in (c). After five iterations, the profile shows a marked improvement in shape (d).
1B ( )
B
B
Fig. 6. Convergence trends observed for optimizations performed at approximately 30 and 65 . The normalized energy of the difference profile is plotted for each of these studies. At 30 (a), the optimization converges for all selected values of k . If the tip angle is increased to 65 (b), the results are mixed. If k 1.0, the algorithm proceeds through a local minimum without settling within it. This illustrates the nonlinear behavior of the spin, t is scaled to be roughly 60% as described by the Bloch equation. If of its full value before addition to t , the optimization is well-behaved. If k is set dynamically, as the square root of the relative convergence metric d, stability is also apparent.
=
1B ( ) B ()
sponds to the associated curve in Fig. 7(a). In Fig. 7(b), the magnetization profile due to the uncorrected RF waveform is plotted. There is substantial smoothing due to the system transfer function. In Fig. 7(c), the candidate RF pulse after ten iterations is shown. The sharp features of the desired boxcar shape are obvious. After ten additional iterations, there is minimal change, as shown in Fig. 7(d). 0.6,1.0, and Fig. 8 illustrates optimization results for at a tip angle of approximately 65 . In Fig. 8(b), choosing 0.6 produces a well-defined square profile with negligible peaks in the stopband. If is set to 1.0, as shown in Fig. 8(c), the resulting profile is distorted with substantial ringing in the stopband. The magnitude of this error is also exhibited in the correyields the result sponding curve in Fig. 6(b). Choosing shown in Fig. 6(d), which is nearly identical to that obtained for 0.6, with slightly larger harmonic ripples in the stopband. An application for this algorithm is refinement of excitations of 90 . The unique nature of the apparent magnetization around 90 makes such optimization difficult. Since the Bloch equations are linear for small changes to existing tip angle, a perturbational approach was adopted. This method requires setting . A flowchart detailing this approach is shown in Fig. 9. To begin, a candidate RF pulse is optimized until it has converged sufficiently for a tip angle below 90 . The power is then increased to 90 and the pulse is optimized further until convergence is attained yet again. Convergence in both cases was defined when dropped below 10%. The results from this optimization are shown in Fig. 10. For this experiment, the intermediate tip angle was approximately 60 ; this tip angle was chosen based on the near-linear nature
46
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
Fig. 8. Graphical comparison of the effects of k on the optimization for a tip angle of approximately 65 . The profile in (a) results from the uncorrected RF waveform. If k 0.6, as shown in (b), the resulting profile is very square. If k is chosen to be 1.0 throughout the optimization instead, the resulting profile exhibits highly nonlinear behavior. By choosing the feedback scalar factor d (d), it is possible to produce results identical to (b), without to be k the additional problem of selecting an optimal value for k . Reduction of the feedback scalar k provided the necessary stability for the algorithm to continue to an acceptable solution.
=
=
p
Fig. 10. 90 optimization for a threshold of 10%. The sub-90 ideal profile is shown in (a), followed by the magnetization profile due to the uncorrected RF pulse (b). After seven iterations, the algorithm terminated itself, producing the profile shown in (c). The power is then set to 90 , where subsequent optimization produces the profile (d).
havior was numerically modeled in a Bloch simulator by truncating the lowest five bits of the applied RF pulse in the feedback algorithm, which resulted in similar periodic ripples. In the case of the sinc pulse, the algorithm attempts to compensate for this behavior by increasing the amplitude of the sidelobes, in the same way as the example shown in Fig. 11. As the optimization progresses, the algorithm attempts to overcome the noise threshold and is largely successful, as shown in Figs. 7 and 8. Since this algorithm was developed originally to produce improved wavelet excitations, an illustration of a partial Haar wavelet family is shown in Fig. 12. After three iterations, it is possible to obtain reasonable scaling and wavelet functions. The ripple behavior discussed previously appears here as well.
IV. DISCUSSION Fig. 9. Flowchart detail of the 90 optimization process. The tip angle is first set to a value below 90 and optimization is conducted until d is below a preselected threshold. The power is then increased to a point where the excitation is nearest to 90 and optimization continues for a few more iterations until the threshold is crossed again. The result is a pulse which is correct for a 90 tip angle. The threshold values are identical for each phase of the process. In Part 1, d is set to 1.0 since it is otherwise undefined. In Part 2, d is set to be 0.5.
of the Bloch equations, as well as the proximity in signal amplitude associated with 90 . In Fig. 10(a), the ideal magnetization profile for the sub-90 case is displayed, followed by the magnetization profile due to the uncorrected RF waveform. After seven iterations, fell below 10%, yielding the profile in Fig. 10(c). At this point, the power was increased to 90 and two further optimization iterations yielded the result shown in Fig. 10(d). The spurious harmonic ripples observed in Figs. 7, 8, and 11 were attributed to low dynamic range in the RF DACs. This be-
From the signals shown in Figs. 4 and 5, the simplicity of the algorithm is evident. Simple signal addition and subtraction of signals combined with scaling techniques produces an algorithm that is straightforward. In Fig. 6, the algorithm exhibits significant divergent behavior for tip angles around 65 when is equal to 1.0 [Fig. 6(b)]. which is added This is attributable to the correction signal . Since the Bloch equations are linear for small perturbato tions of tip angle, better waveforms can be obtained by reducing the scale of the correction signals as the tip angle increases. Selecting a smaller, fixed value for is an obvious modification, produces as shown in Fig. 8(b), but allowing to vary as equally effective results, as illustrated by Fig. 8(d). This is an excellent example of linearity in the Bloch equations; as the tip produce angle becomes large, sufficiently small changes to . Such small changes are correspondingly linear changes in
LEBSACK AND WRIGHT: ITERATIVE RF PULSE REFINEMENT FOR MRI
47
Fig. 11. Illustration of the threshold compensation behavior of the optimization algorithm. In (a), the unoptimized sinc pulse is shown alongside the resulting slice profile (b). Because the lowest bits of the RF DAC contribute virtually nothing to the output pulse, the optimization algorithm overcomes this loss by emphasizing the amplitude of the sidelobes (c). The resulting slice profile after five iterations of optimization shows the small periodic ripples (d).
necessary for optimization of RF waveforms for tip angles of 90 or 180 . The difference in the rate of convergence due to different is noticeable as well. This is to be expected, since large changes will produce correspondingly dramatic changes to . to The final difference energy was observed to be smallest for both tip angles when was set to 0.6. Clearly, the selection of will influence the quality and extent of convergence in this case and this choice needs to be examined further. This algorithm may be useful in practical applications, such , the as the optimization of excitations of 90 . By setting as it algorithm can make increasingly smaller corrections to approaches convergence. As shown in Fig. 10, the algorithm is capable of producing reasonable excitations of approximately 90 when used in a perturbational scheme. The algorithm has been shown to correct phase errors in the acquired profile. Correction of phase is especially important in non-Fourier-based imaging techniques, such as wavelet-encoded MR. In wavelet-encoded MR, the excitation profiles required by the reconstruction algorithm are phase-sensitive in nature and may require correction [16]. Fig. 12 provides evidence of the utility of this algorithm for such excitations. A subtle trend in Figs. 7 and 8 is the maintenance of tip angle. This is achieved by specifying the ideal profile in terms of the normalized magnetization. After 20 iterations, the peak values of the magnetization profile are nearly identical to the corresponding values of the ideal profile. The reason for tip angle consistency is not as obvious. As a corrective scheme, the algorithm should be allowed to alter
p
Fig. 12. Haar Wavelet optimization example. In all cases, the tip angle is d. After three iterations, the Haar wavelets and approximately 30 , k scaling functions are reasonably reproduced.
=
the magnitude of specific points in . Continuous addition of may result in occasions where the the correction signal to dynamic range of the DAC hardware is exceeded. Therefore, it . This is achieved by is necessary to establish a margin for so that its maximum value is somewhat smaller selecting than the dynamic range of the DAC. For the results shown in Figs. 6–8, the initial magnitude was selected to be 70% of the maximum. Examination of several of the figures reveals a recurring phenomenon. Harmonic ripples would evolve and, in most cases, become negligible after several optimization iterations. As discussed previously, this behavior was determined to originate in the RF transmitter signal path. The algorithm is able to partially compensate for this behavior. It is clear that the algorithm is capable of refining RF waveforms where the small-tip-angle assumptions are not applicable. With some modification, the algorithm can produce improved excitations at 90 . A history-based approach where
48
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 1, JANUARY 2002
one considers past acquired and difference profiles might make the algorithm more robust for such excitations as well as inversion pulses. However, this will likely require changes in the algorithm to accommodate potentially confusing behavior at the slice edges. V. CONCLUSION The FT properties of the MR spin system provide the basis for rough RF pulse formulations. Methods have been proposed to further refine the RF pulses obtained in this manner, but only two of these methods have considered effects due to the components of the overall imaging system [17], [18]. The physical transfer function, which is affected by magnetic field inhomogeneities and the spectrometer itself, should be taken into account when constructing RF pulses used to excite precise regions within the sample. In this paper, we present a method of RF pulse refinement which iteratively compensates for the apparent physical transfer function of the system. The technique relies on the formulation of an error profile, which is in turn used to correct the excitation RF pulse. We also discuss methods of filtering and scaling which will provide acceptable results when computing correction signals based upon error signals. The optimization tests presented demonstrate the potential effectiveness of such an algorithm. Simple and complex shapes were reproduced with little effort on an experimental MR system and criteria for convergence was introduced. A damping coefficient was introduced which produced favorable results for larger tip angles and dynamic adjustment of this feedback gain was shown to produce beneficial results. In addition, the authors have shown that the method is directly applicable to the problem of realizing acceptable wavelet pulses. The algorithm may be used on host microcomputers and the effectiveness and efficiency of this algorithm lends itself to implementation on a fast imaging workstation for real-time correction of RF pulses. This algorithm provides a simple method of generating highly specific excitation profiles. Unlike analytical and numerical approaches, the iterative feedback approach presented here does not require compensation for nonlinearities in the RF transmitter. The method should be useful in applications where a series of different shaped RF profiles are needed, such as nonFourier based imaging. REFERENCES [1] A. N. Garroway, P. K. Grannell, and P. Mansfield, “Image formation in NMR by a selective irradiative process,” J. Phys. C, vol. 7, pp. L457–L462, 1974. [2] D. I. Hoult, “The solution of the Bloch equations in the presence of a varying field-an approach to selective pulse analysis,” J. Magn. Reson., vol. 35, pp. 69–86, 1979. [3] J. B. Murdoch, A. H. Lent, and M. K. Kritzer, “Computer-optimized narrowband pulses for multislice imaging,” J. Magn. Reson., vol. 74, pp. 226–263, 1987. [4] J. B. Weaver, Y. Xu, D. M. Healy, and J. R. Driscoll, “Wavelet-encoded MR imaging,” Magn. Reson. Med., vol. 24, pp. 275–287, 1992. [5] P. Mansfield, A. A. Maudsley, P. G. Morris, and I. L. Pykett, “Selective pulses in NMR imaging: A reply to criticism,” J. Magn. Reson., vol. 33, pp. 261–274, 1979. [6] P. R. Locher, “Computer simulation of selective excitation in n.m.r. imaging,” Phil. Trans. R. Soc. Lond., vol. B.289, pp. 537–542, 1980.
B
[7] A. Caprihan, “Effect of amplitude modulation on selective excitation in NMR imaging,” IEEE Trans. Med. Imag., vol. MI-4, pp. 169–175, 1983. [8] M. S. Silver, R. I. Joseph, and D. I. Hoult, “Highly selective =2 and pulse generation,” J. Magn. Reson., vol. 59, pp. 347–351, 1984. [9] S. Conolly, D. Nishimura, and A. Macovski, “Optimal control solutions to the magnetic resonance selective excitation problem,” IEEE Trans. Med. Imag., vol. MI-2, pp. 106–115, 1986. [10] J. T. Ngo and P. G. Morris, “General solution to the NMR excitation problem for noninteracting spins,” Magn. Reson. Med., vol. 5, pp. 217–237, 1987. [11] M. Shinnar, L. Bolinger, and J. S. Leigh, “The use of finite impulse response filters in pulse design,” Magn. Reson. Med., vol. 12, pp. 81–87, 1989. , “The synthesis of soft pulses with a specified frequency response,” [12] Magn. Reson. Med., vol. 12, pp. 88–92, 1989. [13] J. Pauly, P. Le Roux, D. Nishimura, and A. Macovski, “Parameter relations for the Shinnar–Le Roux selective excitation pulse design algorithm,” IEEE Trans. Med. Imag., vol. 10, pp. 53–65, Mar. 1991. [14] S. Pickup and M. Popescu, “Efficient design of pulses with trapezoidal magnitude and linear phase response profiles,” Magn. Reson. Med., vol. 38, pp. 137–145, 1997. [15] N. Gelman and M. L. Wood, “Wavelet encoding for improved SNR and retrospective slice adjustment,” Magn. Reson. Med., vol. 39, pp. 383–391, 1998. [16] W. E. Kyriakos and L. P. Panych, Proceedings of the 5th Scientific Meeting of the ISMRM, Vancouver, B.C., Canada, April 1997. [17] F. Chan, J. Pauly, and A. Macovski, “Effects of RF amplifier distortion on selective excitation and their correction by prewarping,” Magn. Reson. Med., vol. 23, pp. 224–238, 1992. [18] P. S. Christy, R. C. Grimm, J. F. Greenleaf, and R. L. Ehman, “Recursive RF excitation,” Magn. Reson. Med., vol. 26, pp. 361–367, 1992. [19] J. L. Schiano, R. L. Magin, and S. M. Wright, “Feedback control of the nuclear magnetization state,” IEEE Trans. Med. Imag., vol. 10, pp. 138–147, Mar. 1991. [20] P. Devoulon, L. Emsley, P. Weber, R. Meuli, M. Decorps, and G. Bodenhausen, “Methods for reconstructing phase sensitive slice profiles in magnetic resonance imaging,” Magn. Reson. Med., vol. 31, pp. 178–83, 1994. [21] L. E. Franks, Signal Theory. Stroudsberg, PA: Dowden & Culver, 1981.
Eliot T. Lebsack (S’97–M’99) was born in Williston, ND, in 1972. He received the B.S. degree in aerospace engineering and the M.S. degree in electrical engineering from Texas A&M University, College Station, in 1995 and 1999, respectively. He was with Lernout & Hauspie Speech Products in Burlington, MA. Since 2000, he has been working in the area of computational electromagnetics at The MITRE Corporation, Bedford, MA. His interests include the application of high-performance computing techniques to electromagnetics (radio-frequency dosimetry, radar cross section (RCS) prediction).
Steven M. Wright (S’78–M81) received the B.S, M.S., and Ph.D. degrees in electrical engineering from the University of Illinois at Urbana-Champaign in 1980, 1981, and 1984, respectively. From 1984 to 1988, he was an Engineer/Scientist for Magnetic Resonance Imaging at Saint Francis Medical Center, Peoria, IL. He joined the faculty at Texas A&M University, College Station, in 1988 and is currently a Professor of Electrical Engineering and Director of the Magnetic Resonance Systems Laboratory. During the summer and fall of 2000, he was a Visiting Professor at the University of Texas M.D. Anderson Cancer Center. His research interests are in the development of instrumentation and techniques for magnetic resonance imaging and spectroscopy and in applied computational electromagnetics. He has authored or co-authored over 130 journal and conference publications and three book chapters Prof. Wright was awarded the Lockheed Fort Worth Company Award for Excellence in Engineering Education in 1993, the James Stone Faculty Fellow Award in 1999 and a TEES Fellow award in 1999.
