MR physicists and physicians interested in the evolving landscape of MR image reconstruction which uses system calibrations to relax the constraints of conventional Fourier encoding and the mapping methods that provide the needed calibration data (B0 and B1 fields and non-uniform and spatially changing gradients fields). We update what is possible with an improved knowledge of the imperfections in an MR system and how to measure and capitalize on them.

Traditional image reconstruction methods have always made assumptions about things like the uniformity of B0, B1 and flip angles. But modern computational advances allow images to be reconstructed from a much more general forward model, whereby deviations from the ideal can be measured and then simply included into this forward model. The image is computed by taking the inverse of the linear model. In addition to removing artifacts that would have resulted from a simplistic assumption of Fourier encoding, this strategy can actually make advantageous use of the inhomogenieties. An example is parallel imaging. Long thought a nuisance, in the forward-model formalism, the B₁^- inhomogeneity of small surface coils is put to use to help encode an under-sampled image.

The standard imaging equation describes the signal sampled at each point in time given the “image” which is the transverse spin magnetization m_xy(x,y,z) (weighted by whatever relaxation weighting function). The magnetization vector is amplitude and phase modulated by the receive sensitivity of the jith RF receive coil, B_1j^-₍x,y,z), and phase modulated by the integrated history of the gradient waveforms (given by k for linear gradients). Finally, an additional phase factor, Δϕ(x,y,z) can describe any static B₀ inhomogeneity induced phase shift which accrues with time. A volume integration of these factors is performed by the coil detection process. The imaging equation and its discretized form is shown in Fig.1.

Exchanging t for an index, the signal as a function of time can be thought of as a vector. Similarly the “image”, m_xy(x,y,z) can be written as a long vector of all of the voxels. The imaging equation then becomes a matrix multiplication: S = E m, where m is the image vector a vector with length equal to the number of voxels and S is a vector of length equal to the number of samples * number of coils, as depicted in Fig. 2. The encoding matrix E is generated by the forward model to have the same number of columns as pixels in the image (with the same length as m). Since the unknown vector, m, has as many “unknowns” as there are image voxels, E is (hopefully) a tall, skinny matrix, meaning that there are more equations than unknowns and m can be solved for. See figure below. If E is well conditioned, the image information, m, can be solved with any number of numerical methods, employing any number of tricks such as sparsity constraints and pre-conditioners.

This approach provides powerful insight and motivation for knowing the parameters of the forward model, since the image reconstructed this way should “fix” all of the imperfections, giving an accurate reconstruction of m if the parameters of the forward model are accurately measured. Additionally, this framework also sets the stage for seamlessly incorporating non-traditional encoding methods. For example, it gives an insight into how parallel imaging works. In a single channel, well shimmed “ideal” case, E is just the Fourier terms of exp(-ik . r). But if the fewer k-space points are taken than imaging voxels (i.e. kspace is undersampled), the width of the encoding matrix becomes greater than its height (more unknowns than equations). Then an inverse is not possible. Parallel imaging can be viewed as stepping in to save the situation by providing multiple copies of the encoding matrix and stacking them up so E is now taller than it is wide. Simply making copies of the equations is, of course not enough; the equations must be independent. Fortunately, each “copy” of the Fourier matrix is uniquely modulated by the complex coil sensitivity profile B1,j- . In theory 4-fold under-sampling could be reconstructed with 4 coils (producing a square matrix), but in practice the coil sensitivity profiles are not orthogonal enough and more coils are need to make E well conditioned. The g-factor is a measure of the condition-number of the encoding matrix.

In this picture, the coil sensitivities transform from a nuisance shading of the image to a robust and useful part of the image encoding. Recent work has also turned B0 inhomogenieties in a light-weight magnet into an encoding process by rotating the magnet around the head (1), similarly turning “lemons into lemonade.”

The formalism above outlines one of the most general ways to incorporate additional information into the image reconstruction. But, it requires we have a way of quickly obtaining the B0, B1 or gradient field maps.

B0 field mapping Perhaps the easiest of the “nuisance fields” to map, a relative B0 field map can be obtained by comparing the phase of two gradient echo images, since Δϕ = γΔB ΔTE. The main problem is un-wrapping phase jumps if the phase changes by more than 2π over the image. This is a well-explored area with many tools available such as MATLAB’s “unwrap” (2) or FSL’s FUGUE (3). One strategy is to keep ΔTE short so phase wrapping doesn’t happen, but this reduces the accuracy of the Δϕ measurement. A better way is to acquire data at 3 or more TEs, with ΔTE1 long for good estimations and ΔTE2 short to help predict the phase wraps.(4) Note that other issues arise from areas with fat and water (good fat-sat is the first line of defense), and also areas with no signal.

B0 shimming Not all problems can be taken care of in image reconstruction. For example, signal dropout due to a short T2* and across-voxel de-phasing is signal lost to measurement (except for spin- echo type experiments.) To fix the problem before this irreconcilable loss requires B0 shimming. The last 5 years have seen a revolution of sorts for higher order shimming which uses either higher spatial order spherical harmonic shim sets (5) or higher order matrix arrays (6-9), including matrix B0 shim arrays using the same wire loops employed in the RF receive coils (10, 11). An overview will be given of what is now possible with these methods. Figure 3 above gives an indication that at 3T, a 48 or 64channel matrix shim approach can substantially improve B0 homogeneity.

B1 field mapping The coil sensitivities are needed in the forward model approach described above and in many parallel imaging reconstructions (e.g. SENSE). One of the powers of GRAPPA is the use of relatively easy to obtain auto-calibration data instead of coil sensitivity maps. Multiple advances have been made over the years to get good B1- maps that are un-corrupted by underlying tissue structures. To do this, a fast body coil image is often very useful to help identify regions of coil-fall off from regions simply without spins. It is usually helpful to try to extrapolate the B1- map into regions just outside the body to avoid edge effects. Note that since the B1- are guaranteed by Maxwell’s equations to be spatially smooth at some level, B1- mapping is perfect for Compressed Sensing approaches that enforce smoothness on the final coil maps.

In addition to the receive profile maps (B1-), the excitation profile map (B1+) is needed to correct flip angle variations and for use in parallel transmit acceleration of spatially selective RF excitation pulses. (12, 13) B1+ mapping is further constrained in parallel transmit by the need to do it rapidly because 8-16 coils must be independently mapped. There are good reviews of flip angle mapping and the flip angle dependence of sequences. (14) (15, 16) (17-20) The “double angle method” whereby you compare two gradient echo images, one with flip angle set to α and a second image set to 2α. This works fine in regions where the B1+ is large enough so that the 2α excitation is above ~45° where the sin(2α) function starts to become non-linear, but is useless for small flip angles since it utilizes the non-linear nature of the response. The need to wait for full recovery of the Mz magnetization is also a general problem of many methods, although a clever fix for this has been proposed (21).

Much of the field has adopted a variant of either the “Actual Flip-angle Imaging” (AFI) sequence(22, 23) or the, so called “Bloch-Siegert shift” method(24, 25) in which an off resonance RF pulse imparts phase on the signal dependent on the B1+ amplitude. The “Bloch-Siegert method” should properly be called the “Ramsey method” since the phase imparted by an arbitrary off-resonant RF was first described.(26) Described before the invention of NMR in bulk matter, the Bloch-Siegert shift is a special case where the offset is twice the Larmor frequency. (27) MR Fingerprinting has also recently been used to produce B1+ maps.(28)

Gradient trajectory mapping One of the final corrections that can be put into the forward model is to correct the signal phase for imperfections in the gradient field (although there are easier ways to correct gradient trajectory problems (29)). These can be spatial deviations from linearity (whose effect on phase and k-space location) accumulate over time, or they can be dynamic field effects such as gradient eddy-currents or respiratory induced field changes. If the gradient imperfections don’t change from day-to-day, it is sufficient to map the fields and eddy currents in a phantom or with field probes in a calibration step. (30-32) Recent work has concentrated on dynamically monitoring fields with field-probe measurements during each image. (33-35) This approach also allows respiratory induced field changes to be measured and accounted for.