Fast simulation of off-resonance artifacts in MRI using FORECAST (Fourier-based Off-REsonanCe Artifact Simulation in the STeady-State)
Frank Zijlstra1, Job G Bouwman1, Ieva Braškutė1, and Peter R Seevinck1

1Image Sciences Institute, UMC Utrecht, Utrecht, Netherlands


We present a fast alternative to Bloch simulation for simulation of off-resonance artifacts in steady-state imaging. By assuming a steady-state, the signal equation can be quickly evaluated by using multiple Fast Fourier Transforms. We show an acceleration factor of over 350 for a 2D simulation of a titanium cylinder phantom, while the differences with Bloch simulation were minor. The speed of the proposed method enables 3D simulations at high resolution and may benefit various applications.


Numerical simulation of MRI is an important investigative tool for MR research. Examples of its application include pulse sequence design, simulation of various artifacts, such as distortions around metal1, and testing of reconstruction algorithms. Most MRI simulators are based on the discrete time solution of the Bloch equations, commonly referred to as Bloch simulations. As simulation of a large number of spins is computationally intensive, Bloch simulations have been accelerated with parallel computing2 and GPU acceleration3. However, the simulation size still scales quadratically (2D) or cubically (3D) with increasing field of view or resolution.

In this abstract we demonstrate an efficient method to accurately simulate off-resonance artifacts with a lower computational complexity than Bloch simulations, as long as some assumptions about the pulse sequence can be made.



In Cartesian 2D gradient echo imaging with perfect encoding gradients, the basic signal equation in the presence of B0 inhomogeneity is given by:

$$s(k_x,k_y) = \int M_\perp(x,y,z,t',t) \cdot e^{-i2\pi(k_x x+k_y y+\gamma\kern-0.4em- \Delta B_0(x,y,z)t'))} dx dy dz$$

Here, $$$M_\perp$$$ is the transverse signal evolution without encoding, $$$t$$$ is the time since the start of the pulse sequence, and $$$t’$$$ is the time since the last excitation (i.e. for the $$$n$$$'th repetition: $$$t' = t - (n-1) {TR}$$$). Note that $$$k_x = \gamma\kern-0.8em- (t' - TE) G_x$$$.

Assume that $$$M_\perp$$$ does not depend on $$$t$$$, i.e. the signal evolution without encoding is the same during every repetition, which primarily requires a repeating pulse sequence with the MR signal in a steady-state. In this case the signal equation is almost exactly the Fourier transform of $$$M_\perp$$$, with the exception of effects depending on $$$t'$$$ (e.g. $$$\Delta B_0$$$). In this situation, the Bloch simulation is closely related to a naïve evaluation of the 2D Discrete Fourier Transform.

We propose to evaluate this signal equation using Fast Fourier Transforms of $$$M_\perp(t')$$$ for each $$$t'$$$ during data sampling (i.e. one FFT per readout sample). Factors that do not depend on $$$t$$$ can be encoded in $$$M_\perp(t')$$$ , for example: selective excitation, T1 weighting, and T2 decay.

Experimental validation:

We acquired a fast 3D gradient spoiled SSFP scan of a titanium cylinder placed in agar gel (matrix 256x256x90, resolution 1 mm, TE/TR 2.1/7.1, flip angle 30°) at 1.5T (Achieva, Philips, Best, The Netherlands).

We simulated a 2D slice through the center of the cylinder with both Bloch simulations and FORECAST. An analytical model of the cylinder was used to create a spin density and $$$\Delta B_0$$$-map4 with 512x512 spins (i.e. 4 spins per voxel). For the Bloch simulation the actual RF pulses were simulated, while FORECAST only used the RF bandwidth to suppress signal from spins that were not excited. For simplicity, no slice excitation was simulated. We also simulated the full 3D scan with 512x512x180 spins (i.e. 8 spins per voxel) using FORECAST.

Finally we compared the computation times of both methods by performing 2D simulations with matrix sizes varying from 4x4 to 256x256.


Figure 1 shows the acquired scan and 2D simulated images. FORECAST shows only minor differences compared with Bloch simulation, and both methods correspond well to the actual scan except for the field induced by the cuboid-shaped agar phantom, which was not simulated. Simulation time was 2478 seconds for Bloch simulations and 6.5 seconds for our method, on a single core of a Xeon E5-1607 CPU. For the 3D simulation, the simulation time using FORECAST was 883 seconds.

Figure 2 shows the computation times of 2D simulations with increasing matrix size; the acceleration is linear with the matrix size.

Discussion and conclusion

We demonstrated FORECAST, a fast alternative to Bloch simulations for the simulation of artifacts resulting from off-resonance in steady-state pulse sequences. For 2D simulation, acceleration factors up to over 350 times were achieved. Numerical results confirm that the acceleration factor increases linearly with the simulation size. Theoretically, for 3D simulation we expect a quadratic acceleration.

In theory, our method generalizes to various pulse sequences, as long as the aforementioned requirements are met. For example, simulation of spin echo and non-Cartesian trajectories is possible. A hybrid approach including traditional Bloch simulation could be used to simulate additional effects, such as the transient state at the start of an acquisition, or RF excitation profiles. Further evaluation with different types of pulse sequences is warranted to test the method and explore its limits.

We foresee a range of applications that may benefit from FORECAST, including estimation of encoding artifacts around metal, QSM, and background field removal algorithms. The unprecedented speed of FORECAST may pave the way for new applications.


No acknowledgement found.


1. Bakker CJG, Bhagwandien R, Moerland MA, Ramos LMP. Simulation of susceptibility artifacts in 2D and 3D fourier transform spin-echo and gradient-echo magnetic resonance imaging. Magnetic Resonance Imaging 1994;12:767–774.

2. Stöcker T, Vahedipour K, Pflugfelder D, Shah NJ. High-performance computing MRI simulations. Magn. Reson. Med. 2010;64:186–193.

3. Liu F, Kijowski R, Block WF. MRiLab: Performing Fast 3D Parallel MRI Numerical Simulation on A Simple PC. ISMRM 21st Annual Meeting 2013, Salt Lake City, Utah, USA, April 20-26, 2013.

4. Bouwman JG, Bakker CJG. Alias subtraction more efficient than conventional zero-padding in the Fourier-based calculation of the susceptibility induced perturbation of the magnetic field in MR. Magn. Reson. Med. 2012;68:621–630.


Figure 1. Cylinder phantom SSFP scan and corresponding simulation using Bloch simulation and FORECAST.

Figure 2. Graphs showing 2D simulation time (A) of Bloch simulation and FORECAST and the corresponding acceleration factor from Bloch simulation to FORECAST (B) for increasing simulation size (from 4x4 to 256x256).

Proc. Intl. Soc. Mag. Reson. Med. 24 (2016)