Fast Fourier-based simulation of off-resonance artifacts in steady-state gradient echo MRI applied to metal object localization.

PURPOSE: To accelerate simulation of off-resonance artifacts in steady-state gradient echo MRI by using fast Fourier transforms and demonstrate its applicability to metal object localization. THEORY AND METHODS: By exploiting the repetitive nature of steady-state pulse sequences it is possible to use fast Fourier transforms to calculate the MR signal. Based on this principle, a method for fast simulation of off-resonance artifacts was designed. The method was validated against Bloch simulations and MRI scans. Its clinical relevance was demonstrated by employing it for template matching-based metal object localization, as applied to a titanium cylinder, an oxidized zirconium knee implant, and gold fiducials.
RESULTS: The fast simulations were accurate compared with actual MRI scans of the objects. The differences between the fast simulations and Bloch simulations were minor, while the acceleration scaled linearly with the number of phase-encoding lines. The object localization method accurately localized the various metal objects.
CONCLUSION: The proposed simulation methodology provided accurate 3D simulations of off-resonance artifacts with a lower computational complexity than Bloch simulations. The speed of the simulations opens up possibilities in image reconstructions involving off-resonance phenomena that were previously infeasible due to computational limitations, as demonstrated for metal object localization. Magn Reson Med 78:2035-2041, 2017.

INTRODUCTION

Simulation of MRI has proven successful in characterizing off‐resonance artifacts caused by disturbances of the main magnetic field B0 1, 2, 3. These artifacts include geometric distortion, signal loss, signal pile‐up, and phase dispersion 4. In clinical practice, off‐resonance artifacts are increasingly caused by the presence of metal objects, such as implants, surgical screws, or interventional devices 5. The ability to predict artifacts around metal is an important investigative tool for the development of methods aiming at either imaging around metal objects 6, 7, 8 or localization of metal objects 9, 10, 11.

The most general method for simulating off‐resonance artifacts is Bloch simulation, which uses the discrete time solutions of the Bloch equations to simulate the entire pulse sequence for a lattice of isochromats 12, 13. However, the generality of Bloch simulation comes at the cost of a high computational burden, which increases linearly with the product of the number of simulated isochromats and the number of simulated k‐space coefficients. In simulating a 3D scan, both the number of isochromats and the number of k‐space coefficients increase cubically with increasing scan resolution and/or field of view. Both parallel computation 14 and GPU acceleration 15 have been used to accelerate Bloch simulation. However, these types of acceleration only reduce the simulation time by a constant factor, which means that simulating large 3D scans remains impractical.

In this study we propose a fast alternative to Bloch simulation for simulating off‐resonance artifacts in steady‐state pulse sequences. By assuming the signal is in a steady state, it is possible to simulate the MRI signal with a computational complexity that is lower than that of Bloch simulation, while maintaining some of its generality and without requiring analytical descriptions of the phenomena involved. We validated the speed and accuracy of the fast simulation method by comparing it with Bloch simulations and gradient echo scans of a titanium cylinder in vitro. Furthermore, to demonstrate the clinical relevance of the fast simulation method, we employed it for a metal object localization method, inspired by the work of Wachowicz et al. 10, and tested this method on scans of a titanium cylinder, an oxidized zirconium knee implant, and gold fiducials.

THEORY

Bloch simulation can accurately simulate a large number of phenomena in MR simultaneously, such as radiofrequency (RF) pulse inhomogeneity, nonlinear encoding gradients, off‐resonance effects, and chemical shift 13. Depending on the specific research question, not all of these effects need to be simulated, which can simplify the simulation.

In idealized conditions, MRI encoding establishes a perfect Fourier relationship between the measured signal and the image of the object. For example, in a 2D gradient echo sequence with perfectly rectangular encoding gradients, the demodulated MR signal during the readout interval is given by the following formula:

Here, ρ is the effective spin density at time point t, which, for example, can include the effects of RF excitation, proton density, T1 and decay, and off‐resonance ( ). In a perfect steady state with no motion or other time‐varying factors, the effective spin density without encoding is identical in every repetition of the pulse sequence. Then, for simulation of a discrete isochromat grid with proton density ( ), , and T2 as parameters, the signal formula (Eq. (1)) can be rewritten to the following:

Here, the time point t′ is the time after RF excitation. Because is a function of t′ (i.e., the readout direction in k‐space is traversed over time during one repetition) and is independent of t′, this formula can be split into a form where the phase encoding can be evaluated separately as a 1D fast Fourier transform (FFT) for every and every z:

The complexity of this alternative computation of the signal formula is , where the size of the isochromat grid is , and the size of the simulated k‐space is . For Bloch simulation, the complexity is equivalent to naïve summation: . Theoretically, this corresponds to a speedup in the order of . In 3D simulation, phase encoding in both y and z directions can be calculated with FFTs, resulting in a theoretical speedup in the order of .

METHODS

Fast Simulation

We implemented the FFT‐based fast simulation as described in the Theory section, referred to as FORECAST (Fourier‐based off‐resonance artifact simulation in the steady state) 16. We published the MATLAB (MathWorks, Natick, Massachusetts, USA) source code for FORECAST on the MATLAB File Exchange (http://www.mathworks.com/matlabcentral/fileexchange/56680-mri-simulation-using-forecast--fourier-based-off-resonance-artifact-simulation-in-the-steady-state).

To validate the fast simulation methodology, we created a phantom of a titanium cylinder (length, 86 mm; diameter, 24 mm; magnetic susceptibility, χ = 181 ppm) placed in agar gel. The phantom was scanned with a fast 3D spoiled gradient echo scan with non–slice‐selective excitation using the following parameters: matrix = 256 × 256 × 96; isotropic resolution = 1 mm; echo time/pulse repetition time (TE/TR) = 3.2/6.9 ms; flip angle = 30°; feet–head readout direction; readout bandwidth = ∼57 kHz; excitation bandwidth = ∼6.8 KHz; scan time = 2 min, 13 s) at a field strength of 1.5T (Philips Achieva, Best, Netherlands). To test the accuracy of the simulated artifacts under varying conditions, we rotated the phantom in the coronal plane with angles of approximately 0°, 45°, and 90° relative to . Furthermore, we varied scan parameters that directly influence artifact size and shape: 1) the echo time was doubled to 6.4 ms (TR = 11.2 ms); 2) the readout bandwidth was doubled to ∼114 kHz (TE/TR = 2.1/4.7 ms); and 3) the readout direction was changed to the left‐right axis.

Fast simulations were performed for each of these situations to allow a qualitative comparison of the simulated artifacts and the experimentally obtained artifacts. The parameters of the isochromat grid were chosen as follows: The map was calculated from the known magnetic susceptibility of the titanium cylinder using a fast method for forward field‐shift calculation 17. The proton density was defined to be 1 in the medium surrounding the cylinder, and 0 in the cylinder itself. If at any location the value was outside the RF excitation bandwidth, the effective spin density was set to 0 (i.e., no signal originates from those isochromats). The T2 value of the medium was set to 50 ms. The resolution of the isochromat grid was 0.5 mm.

To measure the speed and accuracy of FORECAST relative to Bloch simulation, we simulated a 2D cross‐section of the cylinder using both simulation methods. The scan matrix size was varied from 16 × 16 to 256 × 256, and the simulation times per matrix size were recorded. The other scan parameters were identical to the 3D experiments. The isochromat grid had a resolution of 0.5 mm and was restricted to a single slice without simulating slice selection. In the Bloch simulation, we included 300 dummy repetitions of the pulse sequence to ensure a steady state at the start of encoding. These repetitions were not included in the total simulation time. Furthermore, we simulated perfect spoiling of the transverse magnetization at the end of each repetition. These timing experiments were run on a single core of an Intel Xeon E5‐1607 CPU.

Object Localization

To demonstrate the accuracy of the simulations in a practical application, we implemented a method for localization of known metal objects inspired by the work of Wachowicz et al. 10. FORECAST was used to generate a library of 3D images of a known object in varying orientations. Each of these images was matched individually to an acquired MRI scan using the phase correlation method 18. The voxel with the highest correlation with any image from the library was chosen as the detected position of the object. The orientation of the object is then known from the orientation that was used to simulate the image with the highest correlation.

We applied this localization method to 3D spoiled gradient echo scans of the titanium cylinder, an oxidized zirconium knee implant (Oxinium, estimated χ = 112 ppm), and gold fiducials (length, 5 mm; diameter, 1 mm; χ = −34 ppm). As described in the previous section, the cylinder was scanned in three orientations: 0°, 45°, and 90°. A phantom containing the knee implant was scanned with the phantom in seven different orientations: one scan approximately aligned with B0, and one clockwise and counterclockwise rotation for each image plane. The scan parameters of the knee implant scans were almost identical to the parameters of the cylinder scans, but with matrix size 128 × 128 × 128 and readout bandwidth ∼29 kHz. One scan of a phantom with three gold fiducials was acquired on a 3T scanner (Philips Ingenia, Best, Netherlands) using the following parameters: matrix = 376 × 292 × 75; isotropic resolution = 1.2 mm; TE/TR = 2.7/4.6 ms; flip angle = 10°; anterior–posterior readout direction; readout bandwidth = ∼25 kHz; excitation bandwidth = ∼7 KHz; scan time = 2 min, 10 s. Each of the fiducials was aligned approximately with one of the main axes. With the same scan parameters we also acquired one in vivo scan of a prostate cancer patient, in whom three gold fiducials were placed in the prostate for position verification. A segmentation of the prostate was used as the region of interest. Informed consent was obtained in accordance to the institutional review board of the University Medical Center Utrecht.

For both the cylinder and the gold fiducials a library of 321 different orientations of the objects was generated, corresponding to an angular resolution of approximately 8°. For the knee implant, a library of 891 orientations was generated in a range of −20° to 20° rotation in each axis, corresponding to an angular resolution of approximately 4°. The isochromat grids of both the cylinder and knee implant were generated at a resolution of 0.5 mm. The isochromat grid for the gold fiducials was generated at a resolution of 0.075 mm in a 25.2 mm field of view around the fiducial.

RESULTS

Figure 1 shows the magnitude, phase, and phase difference of the MRI scans and the FORECAST simulations of the titanium cylinder rotated at angles of approximately 0°, 45°, and 90° (Fig. 1A), and with increased echo time, increased readout bandwidth, and with the readout direction changed from feet‐head to left‐right (Fig. 1B). The 3D FORECAST simulations (matrix size = 256 × 256 × 96) took 1045 s on average on a single core of an Intel Xeon E5‐1607 CPU.

Figure 1

MRI scans (first and third columns) and FORECAST simulations (second and fourth columns) of a titanium cylinder rotated approximately 0°, 45°, and 90° in the coronal plane (A) and with varying scan parameters (B). Echo time increased to 6.2 ms (fourth row), readout bandwidth increased to ∼114 kHz (fifth row), and readout direction changed to left–right (LR) (sixth row). The fifth column shows the difference in phase between the MRI scan and the simulations. The outline of the cylinder is shown in red.

In all of these experiments, the fast simulations of the cylinder appeared very similar to the actual scans. In areas of the scan where intravoxel dephasing occurs, the simulation slightly overestimated the signal voids, which may be attributed to the limited resolution of the isochromat grid, as well as possible imperfections in the size and susceptibility of the cylinder model. The phase difference between the simulations and the actual scans was generally consistent with the presence of merely a low‐order background field, which could have been caused by the shape of the agarose phantom and/or system imperfections such as imperfect shimming. Only very close to the cylinder some local phase differences occurred, in particular in areas with signal pile‐ups. Similar results can be seen for simulations of the knee implant in Supporting Figure S1.

Figure 2 shows the simulated images of both the 2D FORECAST simulation and the 2D Bloch simulation of the cylinder. Some minor differences between these images are visible at the boundary of the cylinder, where the full RF simulation in the Bloch simulation produced a more accurate spectrally selective excitation than the excitation bandwidth threshold used in the fast simulation. In Figure 3A, the simulation times for these images are shown for the FORECAST and Bloch simulations for increasing matrix size. Figure 3B shows the corresponding acceleration factors, which are nearly linear with matrix size, as predicted by the theoretical computational complexity.

Figure 2

Magnitude and phase of a 2D simulation of a titanium cylinder using both FORECAST and Bloch simulations. The third column shows magnitude and phase difference maps. The outline of the cylinder is shown in red.

Figure 3

Simulation times for Bloch and FORECAST simulation for increasing 2D matrix sizes (A) and the corresponding acceleration factor (B).

For the scans with the titanium cylinder rotated approximately 0°, 45°, and 90°, the best matching simulations were exactly those with the cylinder model rotated at 0°, 45°, and 90°. These simulations are depicted in Figure 1A in the second and fourth columns. The actual positions of the cylinder are depicted by red outlines, which show that the signal voids do not accurately represent the location of the cylinder.

Figure 4 shows the localization results for the rotations of the knee implant phantom. The detected locations and orientations of the implant are shown as an overlay on these scans. The detected locations generally showed good correspondence with the MRI scan, with the exception of the clockwise rotation in the sagittal plane, where the localized implant appeared to be slightly rotated in the counterclockwise direction.

Figure 4

Localization of the femoral part of an oxidized zirconium knee implant in MRI scans where the implant was rotated counterclockwise and clockwise in the three orthogonal image planes. The outline of the 3D model of the localized implant is shown in red.

Figure 5 shows the gold fiducial detection results for both the phantom and in vivo experiments. In both cases the phase correlation maps showed high correlation at the locations of the fiducials. Both the locations and the orientations of the detected fiducials were accurate when qualitatively compared to the CT scan, as shown in the maximum intensity projection of the CT scan.

Figure 5

Gold fiducial detection in a phantom (top row) and in vivo (bottom row). From left to right: a transverse slice of the MRI scan showing two of the three fiducials (first column), the corresponding phase correlation map (second column), and a maximum intensity projection (MIP) of the CT scan in each image plane with an overlay of a MIP of the detected fiducials in red (third column). In the in vivo scan and phase correlation map an outline of the segmented prostate is shown in red.

DISCUSSION

In this study, we have presented a theory for fast Fourier–based simulation of off‐resonance artifacts in steady‐state gradient echo MRI. Based on this theory, we implemented the FORECAST simulation method and applied it to the simulation of artifacts induced by metal objects in 3D gradient echo scans. The simulated artifacts were found to be accurate with respect to actual scans of metal objects of varying size and shape. To demonstrate the practical applicability of FORECAST, we applied it in a method for localization of known metal objects in MRI.

Our fast simulation methodology allowed the computation of high‐resolution 3D simulations of off‐resonance artifacts within several minutes on regular computer hardware. FORECAST has a computational complexity that is lower than the complexity of Bloch simulation by a factor of the number of phase‐encoding lines. This lower complexity makes our approach especially useful for simulations of scans with large matrix sizes, such as high‐resolution or 3D scans. In our experiments, we were able to perform 3D FORECAST simulations in 1045 s on a single CPU core, which was faster than a 2D Bloch simulation with the same in‐plane matrix size, which required 2572 s. Because computing discrete Fourier transforms and FFTs is the primary computational load of the FORECAST method, the method can be easily accelerated using parallel computing or GPU acceleration.

In general, the off‐resonance artifacts simulated by FORECAST were similar to the artifacts observed in actual MRI scans, in both the magnitude and the phase of the images. Differences between the simulations and scans can be partly attributed to factors that were not included in the simulation, such as accurate simulation of RF excitation and spoiling mechanisms. These effects can be included in the fast simulation with a hybrid simulation approach, in which Bloch simulation is used to numerically establish the steady state of a sequence, which can subsequently be used as an input to the FORECAST simulation. Additionally, there may have been small imperfections in the simulation models with respect to the real objects, since the size and susceptibility values were based on idealized models. The limited resolution of the models may have also caused imperfect intravoxel dephasing and partial volume effects. Given the favorable computational complexity of our method, the resolution of the simulation model can still be increased to reduce these imperfections while maintaining acceptable simulation times.

The results presented in this study only involved steady‐state gradient echo pulse sequences. However, the simulation methodology is theoretically applicable to other types of steady‐state pulse sequences as well. For example, including refocusing effects in the signal formula (Eq. (3)) would allow simulation of spin echo sequences. Simulation of sequences with non‐Cartesian encoding could be achieved through the nonuniform FFT 19. Although the gridding in the nonuniform FFT would increase simulation time, its computational complexity remains better than that of naïve summation. Another possibility is this simulation of sequences with turbo or EPI factors, although this would increase the simulation time proportionally to the turbo or EPI factor.

Although in this study we applied the fast simulation exclusively to the simulation of off‐resonance artifacts for metal object localization, other research areas may benefit from the speed of the proposed simulation method. For example, FORECAST can be used to quickly generate validation images to investigate the influence of susceptibility and off‐resonance effects on reconstruction and postprocessing algorithms, such as B0 mapping or quantitative susceptibility mapping. Furthermore, the simulations could be used in methods that correct for off‐resonance artifacts. For example, the phase difference images in Figure 1 and Supporting Figure S1 are essentially phase images in which the phase dispersion induced by the metal object is largely corrected. This may improve phase‐sensitive reconstructions that have difficulties near strong field disturbers, such as Dixon water–fat separation 20, temperature mapping, or quantification of metal deposits 21. The speed of the simulations also opens up possibilities for using simulations as part of iterative reconstruction approaches.

As a demonstration of the proposed fast simulation method, we implemented a method that localizes metal objects in MRI scans. The localization results on a cylinder, knee implant, and gold fiducials were mostly accurate to within the limits that could be expected from the angular resolution of the template library and the spatial resolution of the images. Generally, this shows that the simulations were sufficiently accurate for this type of application. Automatically locating known metal objects has a variety of potential applications. First, we have shown the possibility of locating gold fiducials. The same procedure can be applied to other types of metal markers and small metal objects such as brachytherapy seeds and sources 22, 23. Second, locating large implants could be used to evaluate arthroplasty procedures and to improve diagnosis of pathology such as implant loosening. Third, the method could be used to locate interventional devices during interventional procedures, although modifications may be necessary to apply the procedure in real time. Finally, gaining exact knowledge about the position and orientation of any metal object could potentially be used to improve imaging and postprocessing in the vicinity of that object.

In conclusion, FORECAST is a fast alternative to Bloch simulation for 3D simulation of off‐resonance artifacts in gradient echo MRI. The clinical applicability of FORECAST was shown in a metal object localization method, which was able to accurately localize a variety of objects. The low computational complexity of FORECAST makes it ideal as an investigative tool and opens up possibilities in image reconstruction involving off‐resonance phenomena that were previously infeasible due to computational limitations.

Fig. S1. Coronal (first row), sagittal (second row), and transverse (third row) slices of an MRI scan (first and third columns) of the femoral part of an oxidized zirconium knee implant and the best matching simulation as found by the proposed localization method (second and fourth columns). The outline of the localized 3D model of the implant is shown in red.

Click here for additional data file.