An illustrated comparison of processing methods for MR phase imaging and QSM: combining array coil signals and phase unwrapping

Phase imaging benefits from strong susceptibility effects at very high field and the high signal‐to‐noise ratio (SNR) afforded by multi‐channel coils. Combining the information from coils is not trivial, however, as the phase that originates in local field effects (the source of interesting contrast) is modified by the inhomogeneous sensitivity of each coil. This has historically been addressed by referencing individual coil sensitivities to that of a volume coil, but alternative approaches are required for ultra‐high field systems in which no such coil is available. An additional challenge in phase imaging is that the phase that develops up to the echo time is “wrapped” into a range of 2π radians. Phase wraps need to be removed in order to reveal the underlying phase distribution of interest. Beginning with a coil combination using a homogeneous reference volume coil – the Roemer approach – which can be applied at 3 T and lower field strengths, we review alternative methods for combining single‐echo and multi‐echo phase images where no such reference coil is available. These are applied to high‐resolution data acquired at 7 T and their effectiveness assessed via an index of agreement between phase values over channels and the contrast‐to‐noise ratio in combined images. The virtual receiver coil and COMPOSER approaches were both found to be computationally efficient and effective. The main features of spatial and temporal phase unwrapping methods are reviewed, placing particular emphasis on recent developments in temporal phase unwrapping and Laplacian approaches. The features and performance of these are illustrated in application to simulated and high‐resolution in vivo data. Temporal unwrapping was the fastest of the methods tested and the Laplacian the most robust in images with low SNR. © 2016 The Authors. NMR in Biomedicine published by John Wiley & Sons Ltd.


INTRODUCTION
Although phase information is frequently disregarded in MRI, every MR measurement yields phase data as the inherent counterpart to the magnitude information arising from the Fourier transform of the complex-valued resonance signal. The phase of T 2 *-weighted gradient echo scans contains a particular richness of anatomical and biophysical information, as has been outlined in other contributions to this special edition. This review focuses on preprocessing steps required to generate combined, artifact-free phase images from multi-channel coils. Its counterpart, by Schweser et al. (1), considers approaches to the removal of the background field. These steps precede the generation of quantitative susceptibility maps (QSMs) (2,3) and susceptibilityweighted images (SWIs) (4,5) and are labeled "phase processing" in the overview of a typical pipeline for the generation of QSMs and SWIs shown in Figure 1.
The phase differs from the magnitude in some important aspects, which means that it requires dedicated methods to allow the information it contains to be faithfully identified. Outside the area of the region of the image occupied by the object, for instance, noise voxels have the same range of values as voxels in the object itself, meaning that the object is inherently more difficult to discern. The phase image is also subject to wraps when the phase, which has evolved from excitation to the echo time, falls outside a range of 2π rad. The phase measured with each coil in an RF array is also subject to the sensitivity of that coil, making it problematic to combine signals from a number of coils and preserve phase information arising from susceptibilities in the object.
Many of the topics covered here have been addressed in a number of other recent reviews (3,6,7). Our aim in this article and the related review by Schweser et al. in this special issue is to complement those by providing a practical, illustrated guide to phase processing. By applying the most relevant and effective methods for phase combination, unwrapping and background field removal to simulated and ultra-high-field (UHF) in vivo data, we hope to afford the reader with clear guidelines as to the most suitable phase processing choices for a wide range of QSM and other phase contrast applications.

COMBINING MULTI-CHANNEL PHASE IMAGES
An assembly of RF coils, each with its own amplifier and receiver channel and designed for simultaneous reception of signal, is called a 'phased array' (8). Phased array coils provide higher signal-to-noise ratio (SNR) than volume coils (8), allow the acquisition to be accelerated using parallel imaging (9)(10)(11)(12) and enable control over patterns of transmit RF (B 1 + ) via parallel transmit excitation (13,14). These features are particularly beneficial at ultra-high static magnetic field (7 T and above) due to inhomogeneous B 0 and B 1 fields.
Complex signals from each coil in a phased array can be combined at a number of stages in the QSM processing. At one extreme, signal can be acquired and processed, from k-space signal to QSM, entirely separately for each channel, and an average QSM over all channels calculated (15). This circumvents the problem of combining raw coil raw phases, which have disparate values due to inhomogeneity in the coil sensitivities, but is computationally demanding, particularly as the number of coils in arrays increases. At the other extreme, coil signals can be combined at the beginning of the receive chain, so that only the combined signals are digitized and recorded (16). This reduces data storage and processing requirements but can lead to destructive interference between the complex coil signals. The approaches in the following sections all combine phase wileyonlinelibrary.com/journal/nbm information from signals which are acquired and reconstructed into separate complex images for each channel. This is the earliest stage in the processing pipeline at which it is possible to apply spatially varying (i.e. non-constant) corrections for the sensitivity of each coil.

The phase measured with a single RF coil
The phase measured with a single RF coil comprises both an echo-time-dependent contribution and a time-independent contribution. The time-dependent contribution arises from local changes to the static magnetic field, ΔB 0 , which originate, at least in part, from tissues with different magnetic susceptibilities. The phase that accrues due to these effects in a measurement with K echoes, at TE k , is where φ 0 r → ð Þ is the echo-time-independent phase offset, sometimes referred to as the transceive phase, γ the gyromagnetic ratio and r → the spatial location vector. Equation [1] neglects measurement noise and second-order, non-linear phase evolution, which manifests if tissues contain multiple water compartments with different resonant frequencies (17)(18)(19). The term φ 0 r → ð Þis generally taken to include all contributions to the phase, including coil sensitivity, which would be measured at a nominal echo time of TE = 0 ms.
The phase that has developed by time TE k , φ r → ; TE k , is 'wrapped' into the inherent encoding range of 2π. The limits of the range are given by (φ L , φ L + 2π], where φ L , the lower limit, can be chosen at will but is usually ascribed the value -π or 0. At a given position and TE k , the measured, or wrapped, phase Φ is given by where mod denotes the modulo operation and the straight brackets indicate rounding down. Wraps are phase isocontours, which should form closed loops within the object or begin and terminate at the object boundary. Wraps which terminate within the object are often referred to as 'open-ended fringelines' (see, e.g., the phase image in Figure 2, with 'No Correction' at the arrow position). Open-ended fringelines are characterized by the presence of a 'residue', or loop of 2 × 2 voxels at the object-end termination of the fringeline for which the sum of the wrapped phase differences is not zero. Open-ended fringelines cannot be unwrapped by adding an integer value of 2π to a circumscribed region, as the residue poses an unwrapping paradox analogous to the Penrose stairs (made famous by Escher's drawings) (20). Residues occur where the phase evolution between neighboring voxels is greater than π (which is more likely to occur at low resolution and where there is a large ΔB 0 T E product) or where noise dominates (21,22). In this latter case, the signal magnitude is close to zero. In a single coil, this can occur because the sensitivity of the coil is low, either due to the distance from the object, because the transverse component of the receive field is zero at the given point or due to coupling with other coils. In combined phase images, the combination of complex data from many coils may sum to zero magnitude (and arbitrary phase) due to interference between the signals. The reader is referred to References (21) and (22) for more detailed explanations of this phenomenon.
The wraps that affect the measured phase frustrate the determination of ΔB 0 r → and also φ 0 r → in some contexts. Approaches to recovering the underlying phase by unwrapping are described in the next section.

The phase measured in each element of a RF array coil
In the context of processing phased array data, it is important to recognize that the phase component that accrues from ΔB 0 is common to the signal in all coils in an array, whereas the phase offset, φ j 0 r → , is different for each coil in the array. That is, for the set of J array coils, Equation [1] has to be rewritten as Figure 2. Comparison of phase combination outcomes for the methods under consideration. The absolute value of the phase-matched complex sum for Echo 2 'Magnitude of phase-matched, combined' (as in Equation [5] in Reference (55)) is shown for each method (other than for the SVD method, where the difference image is shown for Echo 2 À Echo 1), as well as the phase matching quality index (Q), and the combined phase image (all sagittal). The histogram of the mean Q over all echoes and all voxels in the brain is shown for each method (red line), with the Q mean results of the Roemer shown for reference (black line) (note logarithmic scale on the vertical axis).
The phase offset, φ 0 r → , comprises contributions that are common to all coil signals and contributions that are different for each coil in the array. Contributions that are common to all coils include B 1 + phase, eddy current effects, Maxwell terms and a phase gradient in the readout direction caused by mistiming between the gradient and the acquisition (23). In contrast to these effects, the signal measured by each coil is modified by a complex-valued function that reflects the object distance from the coil and the RF wavelength (24): the receive sensitivity or B 1 À . This, and a constant term reflecting the length of the receiver chain, are different for each coil in the array, and need particular attention as they lead to destructive interference in a complex sum of coil signals. The spatial variation of the receive sensitivity of each coil is determined by the interaction with the signals from other coils (coupling) and the wavelength of RF radiation in tissue, which at high field is of the order of size of the brain (circa 30 cm at 3 T and 12 cm at 7 T (25)). Given the conductivity and permittivity of the object, this could, in principle, be modeled for each subject on the basis of the position relative to the coils during measurement adjustment, and used to address the problem of phase combination, although in practice calculations with large arrays are complex and demanding (26). The sign of eddy current and readout gradient phase is dependent on whether signal is acquired during positive or negative readout gradient lobes, leading some authors to specify this contribution in a dedicated term φ l(k) , where for bipolar acquisitions l(k) = 1 when k is odd and l(k) = 2 when k is even (27). Equation [3] and further descriptions here apply to monopolar acquisitions and all even or all odd echoes in a bipolar acquisition.
In the context of QSM, φ 0 r → is considered to contain no interesting information and is treated as a nuisance effect. In fact, the phase distribution of the B 1 field allows electrical properties such as conductivity and permittivity to be quantified, an approach known as electric property tomography (EPT) (28). While EPT is beyond the scope of this review, readers interested in EPT may find the methods, which allow φ 0 r → to be determined explicitly, to be of interest. It is clear from Equation [3] that subtracting the phase offset, , eliminates the coil-dependent term. We will refer to this process, which is the basis of combining phase information from a number of coils, as 'matching' the coil phases.
The remaining, susceptibility-related phase, ϑ j r → ; TE k , is the same for all coils. Equation [3] can be solved uniquely for multi-echo measurements (see later in this section), but is underdetermined for single-echo acquisitions (K = 1). Solutions for this case require additional information from a scan with a reference coil or involve assumptions about the smoothness of φ 0 j r → . These approaches are also described later in this section.

Common operations on phase
The methods in the following sections involve some operations on the phase which it is worthwhile to define here explicitly, to provide readers new to the field with a starting point for implementing phase combination methods. In phase imaging one often needs to calculate the difference between two wrapped phase measurements of φ, Φ n and Φ m .
Subtracting the two, Φ m À Φ n , propagates wraps present in both Φ n and Φ m into the difference. The complex difference, Φ m À n , on the other hand, is only wrapped where the difference between the (unknown) unwrapped phases, φ m À φ n , is outside the range (φ L À π, φ L + π]. The complex difference can be formulated as (a) the angle of the difference between the two wrapped phase values in complex exponential form, where ∠ denotes the four-quadrant tangent inverse (which is usually called atan2 in computer languages), or (b) a simple subtraction of the second angle from the first, with 0, 2π or À2π being added to the difference, depending on the relative sizes of the two angles: : [5] The second expression (Equation [5]) is particularly useful to avoid calculations involving complex numbers, e.g. if the available computer memory is limited.
The maximum and minimum effective frequencies in the rotating frame, Δf MR , max and Δf MR , min , which are not wrapped in the difference image Φ m À n , are given by (Δf MR , max À Δf MR , min ) × (TE m À TE n ) < 1. The significance of φ L in Equations. [2]- [4], which is usually omitted from the literature, is that it may be chosen judiciously to extend the region of the image over which Φ m À n is wrap free. This is pertinent in calculating the difference between phase images acquired at two echo times in the brain, for instance, where (despite shimming) the ΔB 0 values superior to Z ≅ À 20 in MNI space (29) are not distributed equally about zero, but skewed to either positive or negative values, depending on the vendor's convention (30). For example, effective frequencies in the rotating frame of À50 Hz to +50 Hz lead to a complex phase difference without wraps for 10 ms if φ L = 0. Choosing φ L = π/2 would change this range to À25 Hz to +75 Hz, which may be more appropriate for the range of frequencies encountered. Practically, this may allow a wrap-free phase difference image to be generated for a particular region and echo time difference if typical values of ΔB 0 are known.
Another common operation in phase processing is smoothing, either for noise reduction or as a step in high-pass filtering (31).
Smoothing of wrapped phase needs to be performed in the complex representation to avoid blurring wraps: i.e., the real and imaginary components should be smoothed independently. The smoothed phase Φ s can be expressed as Φ s = ∠ (S(M sin (Φ)) + iS(M cos (Φ))), where M is the magnitude and S is an (unspecified) smoothing function.
Phase combination solutions where there is a reference coil: the Roemer/SENSE approach If the MR system has a coil with homogeneous receive sensitivity, such as a body volume coil, the sensitivities of each element of the phased array can be measured with respect to the sensitivity of this (8). A similar approach is used with sensitivity encoding for fast MRI (SENSE) reconstruction, which likewise generates phase and magnitude images with close-to-optimum SNR (9,32). Rather than matching the phases of the coil signals using the phase offset, φ 0 j r → , the complex signals in a voxel are generally combined to a complex image value P, given by P = λp T R À1 b, where p is a vector of complex-valued signal values of all coils in a certain voxel, the vector b contains the complex-valued coil sensitivities in the voxel, λ is a scaling factor and R is the noise correlation matrix (8). For a number of coils J, R is a J × J matrix in which the diagonal elements are the noise levels for each coil and the off-diagonal elements are the correlations between coils x and y (where x and y run from 1 to J) (33). R can be determined from a prescan with no RF excitation. The scaling factor λ is often set to unity, but can be used to correct for inhomogeneities in the sensitivity of the reference volume coil. While this approach is useful at clinical field strengths, UHF MR scanners usually do not have a body volume coil with which to perform the reference measurement needed for the Roemer approach. Some transmit coils are engineered to be able to receive signal for this purpose (i.e. to operate as transceive coils), but these generally do not have very homogeneous sensitivity at UHF. This does not detract from the quality of phase matching, but introduces the inhomogeneous B 1 À phase of the transceiver coil, in addition to the already present B 1 + phase, into the combined image. This can affect the QSM, if not removed (34,35). More problematic is the fact that for many coils, particularly the forthcoming generation of parallel transmit coils, it may be costly or impractical to engineer the transmit array to receive signal due to the need to include transmit-receive switches and preamplifiers. Receiver channels also need to be allocated to these coils, which may limit the number of elements available to the receive-only array. The additional electronics may also disrupt the other elements in the coil, particularly if the array includes B 0 shimming coils (36). These considerations, and the sensitivity of this method to motion between the acquisition of sensitivity maps and the scan of interest, provide the motivation for the following phase combination methods, none of which requires a volume coil reference measurement.

Reference-coil-free solutions for single-echo acquisitions
Frequency filtering (homodyne filtering and unwrap-and-filter) Spatial frequency filtering approaches to phase combination were developed for susceptibility-weighted imaging (SWI). The increasing use of UHF, multi-channel coils and methodological developments in the field have allowed a move from SWI to phase imaging (37), phase value quantification (38)(39)(40) and quantitative susceptibility mapping (QSM) (2,41,42). This has been accompanied by a shift in focus from small veins to ironcontaining gray matter (GM) structures, including those of several millimeters in size, such as the putamen and pallidum (43). The modified contrast generated by frequency filtering for SWI confounds QSM, as QSM interprets phase quantitatively in the ill-posed inverse problem. Frequency filtering methods are included in this review to cover the historical development of phase combination and for completeness.
In vivo phase images comprise a wide range of spatial frequencies. The macroscopic changes to the local magnetic field arising from the interfaces between tissues with large susceptibility differences such as bone, GM, cerebro-spinal fluid and air, particularly in inferior and frontal regions, are dominated by low spatial frequencies, as are phase offsets. Venous vessels -the original focus of interest in venography (4,5) and SWI (44)give rise to high spatial frequencies. This offered the attractive prospect of being able to suppress both phase offsets and macroscopic susceptibility effects by high-pass filtering the data.
The simplest means to achieve this is to divide the complex image data for each channel by a smoothed version of the same data (45). This homodyne filtering not only removes phase offsets but also reduces the effective phase range, removing wraps in the process. While this is quite effective in superior brain regions at low and intermediate fields, it often fails in regions with very strong field gradients, such as regions close to the sinus cavities, where the high-pass filtering is not sufficient to reduce the filtered phase to the range (0, 2π]. The process needs to remove all wraps, as any residual wraps take the form of open-ended fringelines, which cannot be subsequently unwrapped by conventional means (see the next section). This constraint leads to the need to apply relatively strong filtering at high field strength, where ΔB 0 TE is generally higher, leading to heavily modified contrast and often still leaving residual wraps in frontal and ventral regions, making the resulting phase images unusable for QSM. As unsatisfactory as this solution is at very high field, it is still frequently used in commercial MR systems. An improvement to the method is offered by the unwrap-and-filter approach, introduced for single-channel data at 3 T by Rauscher et al. (31) and later applied to combining array coil UHF data by Koopmans et al. (46). Here, phase images from each channel are spatially unwrapped then high-pass filtered before being combined using a magnitude-weighted mean. The removal of wraps allows greatly reduced filtering to be used: only that needed to remove phase offsets. The need to unwrap phase images is a drawback, however, as this process is computationally demanding and prone to errors in the presence of phase noise, which is particularly problematic with small surface coils with inhomogeneous sensitivity profiles. Applying Laplacian unwrapping rather than pathfollowing spatial unwrapping (see next section) could prove advantageous (15), but with the unwrap-and-filter method having been superseded by others described in this section this seems not to have been widely pursued.

Scalar phase matching
A simple approach to estimating the channel-dependent phase offsets, φ 0 j r → , to match the phase images, is to approximate them to channel-dependent constants (47,48). This has been variously called scalar phase matching (SPM), multi-channel phase combination using constant phase offsets (MCPC-C) or the Hammond method. SPM estimates the offset value for each channel from the mean phase value in a matching region of interest, ROI M . The es- where À φ j ROI M ð Þis the mean phase over ROI M . The choice of ROI M is problematic. There must be sufficient signal in all coils in ROI M to estimate φ 0 j at that position. This requires there to be a point of mutual detection. ROI M should be close to the center of the object to minimize the error in estimating the phase offsets, Hammond et al. chose ROI M to be at the center of the image (48). Schäfer and Turner noted the problems that arise when there is no signal at the center of the image (50) and suggested that ROI M be centered on the voxel in which the product of the single-channel magnitude images over coils is a maximum (50). This ensures that there is signal in all coils at this position, but ROI M may be far from the object center.
The fundamental problem with SPM, however, is that the approximation that the phase offsets, φ 0 j r → , are constant over the image is not well adhered to at UHF. The error in this approx- shorter wavelength at higher field, and distance from ROI M . The regional loss of SNR in combined phase images may only be apparent with a larger number of coils (the original work was with an eight-channel coil (48)), and when a quantitative index of phase matching quality is employed (such as Q, introduced later in this section). SPM can break down completely, though, and lead to an open-ended fringeline at a position that depends on the RF wavelength in tissue (i.e. the field strength), the object size and the number and arrangement of coils. Our experience is that this happens circa 5 cm inferior to ROI M with a 32-channel coil at 7 T. As will be seen in the next section, wherever SPM generates a combined image without open-ended fringelines, a virtual coil method can improve on this.

Virtual coil methods
In the absence of a body coil, a virtual reference image (VRI) can be generated from the array coil data itself. In virtual coil methods, the phase of each coil in the array is referenced to that of the VRI, similar to the body-coil image in the Roemer technique.
The virtual body coil (VBC) (16,51) and virtual reference coil (VRC) (52) methods differ in how the virtual image is generated. In the VBC approach, the VRI is the dominant singular vector in a singular value decomposition (SVD) of the array coil data. Originally conceived as a means of reducing the volume of data needing to be reconstructed (16), it was extended to use the VBC as a reference for phase combination (51). The robustness of the VBC approach with different coil designs and field strengths remains to be investigated. In the VRC method, the VRI is generated using SPM.
In both the VBC and VRC approaches, the phase differences between each channel and the VRI are complex smoothed (see earlier in this section) and subtracted from the respective channels to phase match the signals. It might be supposed that any choice of VRI would yield a perfectly matched result. This is not the case, and interesting insights into the phase combination problem are provided by consideration of the choice of the VRI. If the VRI is taken to be zero in all voxels, for instance, the difference between the signal from each channel and that of the VRI is simply the signal from each channel itself, and the VRC approach is equivalent to homodyne filteringthe limitations of which have been discussed above. The sophistication in the VRC approach lies in the fact that the VRI, despite having low SNR away from the center (see earlier in this section), contains all the phase contrast common to all channels, so the difference between each channel and the VRI is small (generally not wrapped) and smooth.
The practical appeal of the VRC approach is that neither a reference coil measurement nor other reference scans are needed, and that deficiencies in the phase matching in the SPM method are remedied in the second step (phase matching to the VRI), as long as the phase in the VRI is well defined. The VRC fails, however, when there are open-ended fringelines in the VRI. That is, it breaks down wherever the SPM method itself breaks down. As such, while it yields phase images with higher SNR than SPMparticularly at a distance from ROI Mit is not more robust than SPM.
More sophisticated solutions for generating a VRI are required for large objects at 7 T, such as the abdomen, and for neuroimaging at higher field strengths. One possibility would be to modify the SPM approach so that, rather than using a set of phase correction constants that yield perfect phase matching only at ROI M , a set of J complex scalar weights could be determined that would yield homogeneous (but imperfect) phase matching throughout the object. This is the receive equivalent of RF shimming in parallel transmission (13,53).
It should be noted that, while phase matching can be excellent with virtual coil methods, the phase corrections applied to each channel are not equal to φ 0 j r → . The fact that the VRI is inhomogeneous means that arbitrary phase fluctuations are introduced into both sides of Equation [3], and the combined phase is subject to the same phase variation, which generally varies slowly in space. In conclusion, the virtual coil approach can always be applied when SPM can be applied, but provides better phase matching. As such, SPM should only be used to generate a VRI (subject to the provisos above) and not the final combined phase image.
Multi-channel phase combination using 3D phase offsets derived from a (separate) dual echo scan (MCPC-3D-II) In the absence of wraps, the phases of two images acquired at TE 1 and TE 2 can be used to calculate the phase offsets of each channel j (see Equation [3] and Reference (54)): Although this approach uses temporal evolution of the phase, we include it under combination approaches for single-echo phase imaging because, in the MCPC-3D-II variant of this method, φ 0 r → are determined from a separate (low-resolution) dual-echo acquisition but applied to a single-echo scan. Combined phase images with this method are free of The main drawbacks with this approach are vulnerability to motion between acquisitions, the need to coregister (or upscale) the low-resolution maps of φ 0 r → with the data to be combined and the need to spatially unwrap phase images. This final consideration makes it problematic to implement on the scanner console and renders it sensitive to errors in low-SNR regions.
Combining phased array data using offsets from a short echo-time reference (COMPOSER) The simple observation underlying coil combination with this wileyonlinelibrary.com/journal/nbm to zero as TE k tends to zero, at which point φ j r → ; TE k approximates to φ 0 j r → . Whereas MCPC-3D-II extrapolates the temporal evolution of the phase to TE = 0 to identify phase offsets, COM-POSER measures an approximation to them directly (55). It is practical to acquire the short-echo-time reference phase data, ; TE; SHORT , in a fast, low-resolution scan with a short-TE sequence such as vTE (56), UTE (57), PETRA (58) or a conventional GE scan with short TE (~1-2 ms; using a short excitation pulse, high bandwidth and asymmetric echo), and apply this to highresolution single-echo or multi-echo data. The steps in COM-POSER are coregistration or upscaling of φ ref;j r → ; TE; SER to the high-resolution scan, and subtraction from the high-resolution data of interest.
COMPOSER requires reference data, and, like MCPC-3D, needs this to be coregistered to the scans of interest. Without the need to unwrap phase data it has significant advantages over MCPC-3D in matching quality, robustness and computational load, and should supplant MCPC-3D in all applications.

Reference-coil-free solutions for multi-echo acquisitions
When phase data is acquired at multiple echo times, the temporal evolution of the signal over the echoes can be used to retrieve the phase offset. Methods found in literature to compute phase images can be classified into two groups: those that reconstruct phase difference images and those that reconstruct phase images for each echo.

Phase difference methods
One approach to eliminating the channel-dependent phase offsets, φ 0 j r → , is to calculate a phase difference for each coil using the phase from one or more pairs of echoes acquired at two different echo times TE n and TE m (59-61).
Complex signals I j r can be combined over channels in the Hermitian inner product (HiP) to generate estimations of the phase evolution between echoes n and m, Θ d r → ; n; m : where κ d j is a weight that accounts for varying noise in different coil channels (61). This method can be performed without the requirement of unwrapping the phase image for each channel, but the SNR is reduced because the weights are retrieved from the images themselves and contain noise. Smoothing the weights decreases the noise as long as the coil sensitivities are not corrupted (62).
When data is acquired with multiple echoes, the SVD method (see Appendix A1) can be applied to the data from all echoes to calculate the coil sensitivities, resulting in an optimum SNR for both the magnitude and phase (63). This approach has been shown to be effective for spectroscopy data (64). The pixel-bypixel SVD factorization of the channel-echo time matrix combines the data from the different coils. Hereby, the first singular value is the maximal coherently constructed signal from all channels and echoes and the eigenvectors contain the coil sensitivity estimations as well as the complex signal, S, of the acquired echoes. The phase of the complex data has an arbitrary offset due to the pixel-by-pixel nature of the method, but not the phase differences between the first echo and each subsequent echo.
Combined field maps are calculated using the unwrapped phase difference between each different echo acquired at TE k and the first echo TE 1 . The final field map is obtained using weighting factors as in Reference (65), where M d is the absolute value of the signal evolution between Echo 1 and echo k. This results in the field map: (see Appendix A2 for a derivation). Note that the weighting factors for phase difference images differ from those for phase images (see later in this section).

Phase imaging methods
Instead of cancelling the phase offsets, φ 0 j r → , by using the phase difference, MCPD-3D-I calculates them for each channel using a multi-echo scan (54), as described earlier in this section. This has the advantage over the HiP approach that φ 0 j r → can be smoothed before being subtracted from the phase at each echo time, leading to a higher SNR result. The computational complexity may be mitigated either by using a low-resolution multi-echo scan to calculate the phase offsets (applying MCPC-3D-II, described earlier in this section, to the multi-echo data) or by downsampling the high-resolution data for the phase offset calculation step, but it is nonetheless subject to the shortcomings of phase unwrapping. In the MAGPI approach (maximum ambiguity distance for phase imaging) (66), the measured phase from Equation [3] is extended to [12] where φ jk noise denotes the additive noise for each coil j and echo time TE k and φ jk wrap denotes 2π phase wraps. This method uses the likelihood function, resulting in an increase of the SNR in the case of a three-echo measurement. The reconstruction of the corrected phase image is performed in three steps. In the first step, an estimation of the most likely tissue-based frequency, which describes the phase difference between echoes, is calculated and removed from the original data. The angle difference is assumed to be Θ d r; n; m ð Þ¼ where the phase wrapping is assigned the integer R, which forces the angle difference to be in the range ( À π , π]. The remaining data is associated with random noise, φ j n , as well as the phase offset, φ 0 j r → , which are separated in the second step. Finally, the most likely tissue frequency that explains the three echoes is calculated. This method estimates the underlying phase without phase unwrapping or denoising and outperforms previous methods for measurements with low SNR, with the drawback of being computationally intensive. In summary, the phase difference method using the HiP approach can be used to combine the phase images from multiple channels. This can be performed very quickly without the need for phase unwrapping, although the SNR of the combined phase image is reduced by the voxel-by-voxel subtraction. Using the 3D correction of the coil sensitivities improves the reconstructed phase image at the cost of computationally intensive phase unwrapping for each channel and the two echoes. The whole dataset is used for the calculation of the complex signal evaluation for the SVD method. Both phase unwrapping and phase offset are incorporated in the maximal likelihood (MAGPI) calculation, which performs better at low SNR.

Combining matched phase signals over channels and echoes
In the previous sections we have looked at ways to estimate or approximate phase offsets. The individual phase images can be phase matched by subtraction of the phase offset contribution: The matched phases from each coil must then be combined over channels. Optimal SNR is achieved by including both coil sensitivities and consideration of noise correlation between signals from different channels, as described in the spatially matched filter approach (8,67). In combining phases, many authors omit the effects of noise correlation and simply sum the phase values in complex exponential form, with weights κ j r → ; TE k : TE k e ÀTE k = T 2 Ã ; [14] which weights each echo by its SNR (68), as used previously in multi-echo EPI (69). The weighting factors for phase difference methods are given earlier in this section.

Assessing the quality of combined phase images
In the general absence of a ground truth, assessing the quality of combined phase images is not trivial. In the literature, the absence of open-ended field lines or the fact that phase images can be unwrapped is often mistakenly seen as being indicative of effective matching and combination of individual channel phase images. However, this takes no account of the loss of SNR that occurs when phase matching is imperfect but has not led to complete destructive loss of signal. We have suggested that the quality of phase matching be represented by the metric Q (54).
When the phases of the individual signal vectors are in good agreement (i.e., they are matched), the length of the complex sum of the signals (numerator) is equal to the sum of the length of the individual signal vectors (denominator), i.e. Q = 100. Q should only be evaluated in voxels containing signal. In background, or noise voxels, the complex sum (the numerator) tends to zero for large j (a desirable property), whereas noise is additive in the magnitude sum (the denominator).

A quantitative comparison of phase combination methods
To allow the phase combination methods described in earlier sections to be compared quantitatively, we applied them to a human in vivo multi-echo data set acquired with a 7 T Siemens MAGNETOM (Siemens Healthcare, Erlangen, Germany) scanner and a 32-channel Nova Medical coil. This was a triple-echo gradient-echo acquisition with a matrix size of 512 × 512 × 208 (voxels of 0.4 × 0.4 × 0.7 mm 3 ), monopolar readout and TE = [8.0, 14.0, 21.0] ms, T R = 26 ms. The bandwidth was 250 Hz/pixel, TA = 10 min 17 s. Phase images were reconstructed according to Equation [13], with the magnitude data as weights, κ j . Noise correlation between coils was not included in the reconstruction. The phase data were combined with no correction for phase offsets ("no correction", to illustrate the phase combination problem) and with the phase imaging methods Roemer, SPM, VRC, MCPC-3D-II and COMPOSER and the phase difference methods HiP and SVD (see earlier in this section).
The effectiveness of each method was assessed via the quality of phase matching between channels, Q, and the contrast-tonoise ratio (CNR) between GM and white matter (WM). Exact correspondence of phase values to the (scaled) B 0 fieldmap was not considered essential, as the additional slow phase variations introduced by some methods may be attenuated by later steps in QSM processing (3).
Maps of Q were calculated for each echo (for illustration), and the mean of Q for each voxel inside the brain over all echoes, Q mean , for quantification. Histograms of Q mean and the median value of Q mean over the brain were calculated for each method. GM-WM CNR was also assessed in six ROIs. The GM ROIs were thin strips of cortex; the corresponding WM ROIs were similarly sized and immediately neighboring. The six regions were distributed throughout the brain. Two-sided paired t tests were calculated between the GM-WM CNR measured with the Roemer method (which was treated as the reference) and each of the other approaches.
To assess the effect of phase combination on backgroundcorrected phase and QSMs, background correction with V-SHARP was performed with a maximum kernel size of 9 mm and high-pass-based regularization of V-SHARP with a cut-off frequency of 0.0089 mm À1 (parameters determined to be optimal S. D. ROBINSON ET AL. wileyonlinelibrary.com/journal/nbm in the review by Özbay et al. (70) in this special issue). Phase images were resliced to isotropic resolution before applying V-SHARP using linear interpolation in the spatial domain. The same internal region for background correction was used for all techniques, defined based on the Roemer-reconstructed magnitude using FSL BET. Susceptibility maps were reconstructed with a simple k-space-based algorithm (71).
Phase matching was high throughout the brain with the Roemer method (median Q mean of 98.7%), the reference method.
The median value of Q mean was 18.7% without removal of φ 0 r → ("No Correction"), with an ostensibly arbitrary distribution of regions of destructive interference between coil signals (visible as low Q values (blue) and open-ended fringeline at arrow). This was improved with SPM (median Q mean of 50.0%), where the phase matching was perfect at ROI M , at the center of the brain (black square in Q map). The directional dependence of the coils on Q is apparent, with higher values at the top of the brain, where coils point predominantly in the same direction (headfoot), than at the same distance from ROI M near the bottom of the brain, where they are arranged around the object (see values close to zero in the cerebellum and noise in the phase at arrow), leading to destructive interference. Phase matching was high and uniform with VRC (median Q mean of 98.8%). With MCPC-3D-II, Q values were close to those of VRC (median Q mean of 97.8%), although isolated weaknesses from unwrapping errors were apparent at arrow positions. Q values were high and uniform with COMPOSER (98.6%) and SVD (99.2%). Within the brain, HiP phase values (not shown) were very similar to those with SVD (typically within 0.1%). Q was high in both the image and the background in SVD, as the method maximizes Q in every voxel. In summary, phase matching was as effective with VRC, COMPOSER, HiP and SVD as with Roemer (see histograms).
The GM-WM CNR was not significantly different between Roemer and the other phase imaging methods: SPM, VRC, MCPC-3D and COMPOSER. The GM-WM CNR was 41% less on average, however, with the two phase difference methods HiP (2.27 ± 0.75) and SVD (2.28 ± 0.75) than with Roemer (3.87 ± 1.06) (a significant difference; p = 0.002). The echo times used in this comparison are not ideal for HiP and SVD, however. For the phase difference methods, the first echo time should be as short as possible and the last echo time close to the T 2 * of the tissues of interest (see Equation [11]).
Background-corrected phase images and susceptibility maps from Echo 2 of each phase combination method are shown in Figure 3. Isolated patches of noise (at Arrows 1 and 2) are present in the No Correction and, to a lesser extent, SPM results. MCPC-3D-II results are subject to errors in frontal areas (see Arrow 3), which arise from inclusion in the internal region of background correction of inferior regions in which there are isolated abrupt phase changes, where unwrapping failed (see arrows in Figure 2). If included in the mask, these cause background field estimation errors. These pervade up to 4 cm superior to the source of the problem, in this case. This effect could be avoided by excluding regions of low Q from the mask. SVD results are uniformly noisier than results with the Romer, VRC and COMPOSER methods.

Practical considerations
The most suitable phase combination method depends on whether a volume reference coil is available, whether multi-echo measurements are made and whether it is, for a particular application, important to remove all sources of phase other than those arising from susceptibility effects.
The requirements, features and shortcomings of the phase combination approaches described here are summarized in Table 1. In brief, Roemer remains the method of choice if a homogeneous reference coil (e.g. body coil) is available. The Roemer method can also be used if only an inhomogeneous transmit coil can be used to receive signal (e.g. head coil at UHF), although the transmit B1 contributions to the combined phase are not generally harmonic, so are not removed by background correction unless this entails regularization (35,70).
VRC and COMPOSER proved to be the most effective methods for single-echo acquisitions where there is no volume reference coil. VRC needs no reference measurement, but (like the SPM method on which it depends) fails where the VRC image yields open-ended fringelines. COMPOSER is applicable to all coil designs and field strengths but requires a reference scan.
Despite the fact that multi-echo acquisitions contain the information to identify φ 0 j r → ; TE k and ΔB r → uniquely and without recourse to reference measurements, the methods considered here -HiP, SVD and MAGPIall have significant shortcomings. HiP and SVD have reduced CNR, particularly if the echo times are not optimized for phase difference calculation, and MAGPI is computationally intensive. We see the need for computationally undemanding solutions to the multi-echo problem that yield combined images without reducing CNR. The single-echo methods VRC and COMPOSER can, of course, be applied to multi-echo data.

Outline of the problem and requirements
Wraps may be present along any of the four dimensions of a multi-echo data set. The problem of unwrapping is ill posed, as many values of true phase lead to a wrapped value (i.e., it is surjective but non-injective). In practice, assumptions about spatial or temporal continuity (or both) can be drawn upon to attempt to solve the problem. The success of any method is dependent, however, on how rapidly phase changes between voxels or time points and the level of noise. Spatial unwrapping methods (see later in this section) draw on the fact that phase generally changes slowly from voxel to voxel unless a wrap has occurred. Spatial approaches can be subdivided into path-following methods and 'Laplacian' unwrapping (see later in this section). Path-following methods determine a wrap to have occurred if the phase change between two voxels along a path is greater than π. The value of +π or Àπ is added to all further voxels along the path (depending on whether the phase jump was positive or negative). Unless errors occur, this restores the exact value of the underlying phase. Laplacian unwrapping attempts to identify the unwrapped phase whose local derivatives are most similar to the derivatives of the wrapped phase. This removes discontinuities in the phase butunless an exact solution is implemented (see later in this section and Appendix A3)it modifies all phase values in the image in a spatially slowly varying way, introducing some background field suppression (71). This does not propagate into the background-corrected field (71), so does not affect QSM results, but may prove problematic if the phase is being used quantitatively without background elimination (e.g. to determine φ 0 r → (55)). Finally, phase change between echoes can be used to perform a temporal unwrapping of the data (see later in this section) in a 1D unwrapping approach.

Path-following methods
An assumption underlying spatial phase unwrapping is that a wrap has occurred between two neighboring voxels if the difference between them is greater than π. Algorithms can apply this assumption in a number of ways in one, two or three dimensions, unwrapping the phase along paths guided by the quality of the information that voxels contain (72) and the need to avoid 'branch cut' lines (73). Branch cut lines are imposed barriers to paths. They connect problematic path-dependent 'residues', such as occur in 2 × 2 voxels containing the termination of an open-ended fringe line. The result of unwrapping the phase along any path that does not cross such a branch cut line is independent of the path. Unwrapping in a higher number of dimensions reduces the sensitivity to noise but increases the complexity of the calculation. Both 2D and 3D algorithms are in use in MRI. The number, complexity and range of these makes it impossible to cover them in depth in this review; the reader is instead referred to Reference (22) for a comprehensive overview. Instead, we try to  (75), BEST-PATH (76) and PHUN (72). We quantitatively compare the performance of these algorithms with Laplacian and temporal methods in unwrapping highly wrapped and noisy distributions (see later in this section).

Laplacian phase unwrapping
Laplacian unwrapping is distinct from other spatial and temporal unwrapping methods because, while it preserves much phase variation in the image, it does not generally yield quantitative phase values. In the light of its quite recent application in MRI (from the field of interferometry), we dedicate additional space to a comprehensive and rigorous explanation of this method to explain its function and features but also explore the potential to extend it to yield quantitative results. The principal idea behind Laplacian phase unwrapping (as well as its variants) is the observation that, assuming that the unwrapped phase data φ unwrap is sufficiently differentiable, some differential operator applied to the unwrapped phase can be obtained from the wrapped phase φ wrap (77). As the name suggests, this differential operator is the Laplacian, i.e.
as exp(iφ unwrap ) = exp (iφ wrap ), we have, with ψ = exp (iφ wrap ), that Δφ unwrap ¼I ψ À1 Δψ À Á ; [16] where I is the imaginary component. The idea is now to solve the corresponding Poisson equation, i.e. to obtain φ unwrap from Δφ unwrap . However, at this stage, φ unwrap is only determined up to harmonic phase maps, i.e. those whose Laplacian vanishes, so additional information is needed to obtain a solution for the unwrapped phase. One way to incorporate such information is, for instance, to assume that φ unwrap is periodic. This may be achieved by mirror extension or, equivalently, by assuming homogeneous Neumann conditions, i.e. ∂φ unwrap ∂ν ¼ 0 on the boundary. Then, φ unwrap is uniquely determined up to a global constant, which can, for instance, be determined by fixing φ unwrap at a single point or fixing the integral ∫φ unwrap . Alternatively, some applications such as QSM do not require boundary conditions, as the unwrapped phase only has to be determined up to harmonic components, which will in any case be removed in the subsequent background elimination step (see the review article by Schweser et al. in this issue (1) and Reference (70)).
On a rectangle or cuboid, one can express the solution of the Poisson equation with the help of Fourier cosine series. Denoting by F the Fourier analysis operator and by F À1 its inverse, the Fourier synthesis operator, we compute, with the help of Equation [16], the Laplacian as where ω j j 2 ¼ ω 2 1 þ ω 2 2 þ ::: þ ω 2 n denotes the squared Euclidean norm of the nD Fourier indices ω = (ω 1, ω 2, . . ., ω n ). The solution of the Poisson equation may then be computed via where, for ω = 0, no division should be performed (which fixes ∫φ unwrap ). While this approach leads to exact phase unwrapping in the continuous setting, its practical application is limited by the fact that we are dealing with discrete data, usually on a uniform grid. Discretization is therefore necessary. In view of Equations [17] and [18], this means replacing Fand F À1 by discrete counterparts, as proposed in Reference (77). Doing so either involves the (normalized) discrete cosine transform (DCT) of Type II and its inverse, which is the DCT of Type III (78), i.e., F = DCT II and F À1 = DCT III . Alternatively, an equivalent operation such as mirror symmetrization and the fast Fourier transform can be utilized (79). However, Equations [17] and [18] fail in the discrete setting, so applying them only yields an approximation φ apprx , and, in particular, the difference between φ apprx and φ wrap will in general not correspond to integer multiples of 2π, a property that is satisfied by the exact difference φ unwrap À φ wrap . If one is interested in quantitatively exact results, one has to resort to heuristics. One such approach would be to round the approximation φ apprx À φ wrap to multiples of 2π and add φ wrap , which gives an approximation to φ unwrap that is exact where the rounding is exact. Of course, this may or may not be the case, and the approach works particularly well if the error φ unwrap À φ apprx was already small. A more sophisticated method is to perform the following congruence operation.
Let D j = the number of discontinuities in φ apprx j end.
• Find j min for which D is smallest.
There are several possibilities to detect discrete discontinuities; for instance, one can say that a discontinuity occurs if the values in two neighboring pixels differ at least by 2π. For noiseless data, the above-mentioned single-step congruence operation with h = 0 is usually sufficient. However, testing N global shifts h 1 , … , h N may increase the accuracy of the congruence operation in the case of noisy data.
Alternative representations to Equation [16] that are equivalent in the continuous setting give rise to different discrete variants of Equation [17]. For instance, in Reference (80), computing Δφ unwrap is proposed according to where a Á b denotes the scalar product between two vectors. Together with Equation [18], it constitutes another phaseunwrapping approach of Laplace type. The two variants have in common that they involve only a few deterministic steps with well-understood mathematical properties. One can show, for instance, that DCT-based discrete versions of Equations [17], [18], [19] are Lipschitz continuous. This means in particular that these computations are stable with respect to noise, in the sense that noise is amplified up to a finite factor, which turns out to be moderate in practice. Besides Fourier methods, the approximation of Δφ unwrap can also be performed only in the spatial domain by employing finite-difference techniques. One advantage of these approaches is that the solution step for the discrete Poisson equation can be evaluated in a quantitatively exact manner, i.e., φ unwrap may be obtained exactly if one is able to compute its finite-difference Laplacian Δ h φ unwrap . Here, again, one has generally only an approximation Δ h φ apprx available. However, the Poisson inversion process is stable, so the result φ apprx is an approximation of φ unwrap and the error can be controlled. Furthermore, by employing heuristics, it is possible to compute Δ h φ unwrap for wrapped data with moderate phase wraps, leading to exact results. We refer to Appendix A4 for details.
In practice, applying unweighted finite-difference techniques (as in Appendix A4) leaves substantial errors, especially in the ventral brain. These are mainly caused by noisy phase data in the background where the magnitude is small. Introducing binary weights w (1, tissue; 0, background) and considering the weighted finite-difference operator div h (w∇ h ) instead of the Laplacian (A4, Equation [21]) allow this noisy background phase to be masked out, eliminating the main error source. Such a modification still leads to a linear equation, which can then be represented by a sparse matrix, calling for iterative solution strategies. A particularly well-suited procedure is the preconditioned conjugate gradient (PCG) method with DCT-based solution of the unweighted Poisson problem (Equation [23] in Appendix A4) used as preconditioner in each iteration step. The pseudocode for this method is presented in Appendix A3, and a MATLAB implementation is also available (81). We refer the reader to Reference (22) for an extended description. This PCG algorithm is typically characterized by a rapid convergence with welldefined termination conditions (such as a user-defined maximum number of iterations and convergence factor). The results obtained with the weighted Laplacian approach by virtue of the PCG method tend to be much closer to the exact solution than those for the unweighted approach. Here, applying the previously mentioned congruence operation also yields further improvements.
In summary, Laplacian phase unwrapping techniques constitute a fast and easy-to-implement class of methods that are robust to noise. This robustness generally comes at the cost of obtaining only approximations of the exact phase. The latter issue can, to a certain extent, be resolved by finetuning the discretization, introducing appropriate weighting and fast linear solution techniques as well as employing heuristics and congruence operations. Generally, the methods discussed can be recommended for data sets with moderate to low SNR where it is not crucial to obtain exact quantitative values.

Temporal phase unwrapping
Temporal unwrapping uses the evolution of the phase over measurements at a number of echo times to remove wraps that are caused by time-dependent effects, i.e. due to ΔB 0 r → . It therefore requires the information available only in multi-echo acquisitions, and, until recent extensions, has only been capable of removing wraps that occur over a time greater than the interecho period, and only generates unwrapped phase difference images, not unwrapped phase images. If no a priori information is assumed about the sign of the phase change between two echoes (φ L = À π in Equations [4], [5]), then the range of the phase difference that is wrap free is (Àπ, +π]. This limit corresponds to a difference between echo times of (2 Δf MR ) À1 s. With a frequency range of interest of ±50 Hz/T in the head, the inter-echo time needs to be less than 10 ms T, or circa 1.4 ms at 7 T and 3.3 ms at 3 T. This poses an undesirable constraint on the acquisition, and one which may be at or beyond achievable gradient switching rates and physiological dB/dT limits. This, and the fact that phase difference images have poor noise properties compared with phase images (see the previous section, and Equation [20] in Appendix A2), have seen temporal phase unwrapping largely treated as obsolete for UHF applications. Two approaches have recently been presented that overcome both of these restrictions, however.
In UMPIRE (Unwrapping Multi-echo Phase Images with Irregular Echo spacings) (82) and MAGPI (83), the authors show that the limitation on the inter-echo time described above can be overcome by using more than two unequally spaced echoes in a multi-echo train. In UMPIRE, a short additional delay δTE is inserted between one pair of echoes (e.g. between Echoes 2 and 3). The phase evolution in this period, Φ 3 À 2 À Φ 2 À 1 , contains no wraps if δTE, which can be chosen at will, satisfies δTE ¼ wrap-free estimate of Δf MR , which is used to identify and remove wraps in Φ 3 À 2 and Φ 2 À 1 and between all further pairs of echoes. Any phase offset is then identified and removed before, in a final step, unwrapping the phase images themselves. In MAGPI, the authors focus on the optimum spacing of echoes to disambiguate the estimation of Δf MR , which is identified as the closest L2 distance between the probability distributions representing the possible field map values for each pair of echoes (83). In general, temporal phase unwrapping has as advantages that it is unaffectedin both time and accuracyby the complexity of the wrapped topography, and as disadvantages that it is more prone to noise than spatial methods (as it is based on phase differences) and that it is dependent on linear phase evolution, which can be disrupted by tissue characteristics (17), uncompensated flow or idiosyncrasies in the reconstruction.

An illustrated comparison of phase unwrapping methods
In order to illustrate the features of the spatial, Laplacian and temporal phase unwrapping methods described in the previous section, these were applied to simulated data and in vivo multiecho phase data acquired at 7 T.
Two types of phase distribution were simulated. The first was a complex phase topography with no noise, the second a simple 3D Gaussian distribution with varying amounts of noise. These were designed to test the methods' ability to resolve complex patterns and to assess unwrapping accuracy in a low-SNR context respectively. Each phase distribution was modeled as having a purely linear evolution and simulated at three echo times (TE = 5, 10, 16 msthe uneven echo time to allow temporal unwrapping with UMPIRE). More details about the simulated data are given in Reference (82). The complex topographies shown in the upper part of Figure 4 correspond to "Simulated data 2" with complexity level 4 in Reference (82).
The complex topography comprises slices with differing numbers of wraps and complexity. A sample slice is shown at the top left of wileyonlinelibrary.com/journal/nbm isotropic voxels, TE = [7.4, 14.5, 20.8] ms, that is, with an additional delay of 1 ms between the second and third echoes.
Of the spatial methods, PRELUDE (74) was the most accurate in unwrapping both simulated and in vivo data, but was slow where the topographies were noisy. Cusack's method (75) was prone to errors for high noise and suffered isolated errors in in vivo data (at arrow position). BEST PATH (76) was relatively robust. PHUN (72) failed to unwrap most slices in the complex shape and some regions of the in vivo data and was prone to the most errors at high noise. Unweighted LAPLACE unwrapping (84) was no better than other methods in complex topographies, but was reliable in high-noise and in vivo data. The HiP failed where the limit (2 Δf MR ) À1 was exceeded (at arrows). UMPIRE performed poorly in high noise but had no errors in complex topographies and was free from errors in in vivo data above MNI z = À14.
The amplitude of residual phase ambiguities in the weighted Laplace solution was assessed by comparison with PRELUDE, which, without apparent errors in the slices evaluated, was taken to yield an exact solution. Results for the third echo (with TE = 20.8 ms) in two ventral slices are shown in Figure 5. The residual phase errors were assessed for the unweighted Laplace solution, and the weighted Laplace solution without and with the applied congruence operation. Results for two ventral slices are shown in Figure 5. Histograms of phase differences over the whole volume for each of these cases are also presented. Residuals for the unweighted Laplace solution are generally slowly varying. They increase in spatial frequency in the inferior slices, however, where values of ±6 rad are not untypical. Those residues were substantially reduced after six steps of the PCG algorithm for the unweighted Laplace approach, which yields an asymmetric distribution of residual phase errors, with the values mostly within [À5; 1] rad. The congruence algorithm left no residues in most voxels, indicating an exact solution. The exception was a small number of voxels in noisy inferior regions close to the brain boundary, where the residual errors in the weighted Laplace solution were large and led to a wrap in the congruent solution (at the arrow position). Those residuals differ by multiples of 2π from the exact solution, illustrated in a logarithmic plot inset in the top right corner of the histogram.
Any assessment of the performance of a phase unwrapping method is dependent on the features of the test object, as well as how the method is implemented and the choice of  wileyonlinelibrary.com/journal/nbm parameters such as thresholds, seed voxels etc. The reader is referred to References (82), (85) and (86) for other comparative evaluations in a broad range of contexts.

CONCLUSION
The phase combination methods compared in this review differ greatly in performance and applicability. Despite a widespread understanding that the Roemer method cannot be used at UHF because of the absence of a body coil, we find that effective phase matching can be achieved with this approach using an inhomogeneous transmit-receive coil. The transmit phase is transferred to the combined image, however, andas a non-harmonic phase variationis not generally removed by background correction. The most effective alternative solutions, applicable to cases in which the transmit element is not engineered to receive signal, were the VRC method and COMPOSER. Multi-echo approaches were found to match the phase well but be computationally expensive or reduce GM-WM CNR.
A comparison of phase unwrapping methods showed many path-following spatial methods to be fragile in the presence of rapid phase fluctuations. Temporal unwrapping, while fast, was more prone to errors in low SNR. Laplacian unwrapping was both fast and effective when applied to simulated data with low SNR and in vivo data acquired at 7 T. Center for Advancing Translational Sciences of the National Institutes of Health under award Number UL1TR001412 (the content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH). B. Dymerska was additionally supported by a DOC fellowship from the Austrian Academy of Science.

APPENDIX A1: SINGULAR VALUE DECOMPOSITION
Let A be an n × m complex matrix, σ a non-negative scalar, u a nonzero vector of size m and v a non-zero vector of size n, such that Av ¼ σu: In this case, u and v are called the left-singular and rightsingular vectors, respectively, for the singular value σ.
The matrix product The number of non-zero singular values is the rank r of the matrix A.
U and V are unitary matrices (U * U = U U * = I, with I the identity matrix) and Σ is a diagonal matrix.
For an SVD for an n × m matrix A of rank r, the following applies: the r diagonal elements of Σ are the square roots of the r eigenvalues of A * A the columns of V are the normalized eigenvectors, v i , of A * A to the eigenvalue σ i and satisfy Av i = σ i u i the columns of U are the normalized eigenvectors, u i , of AA * to the eigenvalue σ i and satisfy A * u i = σ i v i .
With these properties, the SVD can be interpreted as a transformation of correlated variables (matrix A) into a set of uncorrelated ones (the eigenvectors for different eigenvalues are orthogonal to each other).

APPENDIX A2: DERIVATION OF OPTIMUM WEIGHTING FACTORS IN COMBINING PHASE DIFFERENCE IMAGES
The noise of a phase image acquired at an echo time TE k , φ k , is given by The noise of the phase difference, Θ k1 , between φ k and φ 1 can be calculated by error propagation: The noise of a frequency map, ω k1 ¼ Θ k1 =TE k À TE 1 can then be expressed by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi e 2 TE k ð Þ=T2 Ã þ e 2 TE1 ð Þ=T2 Ã p : . [20] APPENDIX A3: PRECONDITIONED CONJUGATE GRADIENT ALGORITHM FOR WEIGHTED LAPLACIAN UNWRAPPING The main task for weighted Laplacian unwrapping is to solve the finite-difference equation L h φ apprx = b wrap , where the forward operator is given by L h = À div h (w∇ h ) with binary weight w, b wrap = À div h (w∇ h φ apprx ) and ∇ h φ apprx is given by Equation [24]. Here, the operator L h can be represented by a sparse matrix. The PCG solution to this problem can be described as a pseudocode in the following steps. Initialize the iteration step: j = 1, the convergence factor: ϵ j = 1 = 1, and the solution matrix: φ apprx 1 ¼ 0. Compute the weighted Laplacian of the wrapped phase b wrap (as described above), set r 1 = b wrap .
Let t j ¼ ∑ k;l r j À Á k;l φ precon j k;l , where k , l represent indices in 2D space.
If j = 1, then Calculate the weighted Laplacian of p j , i.e. q j = L h p j . Let α j ¼ tj ∑ k:l p j ð Þ k;l Á q j ð Þ k;l .
Let the solution at step j be: φ apprx jþ1 ¼ φ apprx j þ α j p j . Let r j + 1 = r j À α j q j .
Replace j with j + 1. end.

APPENDIX A4: FINITE-DIFFERENCE TECHNIQUES FOR LAPLACIAN PHASE UNWRAPPING
Besides Fourier methods, the approximation of Δφ unwrap can also be performed only in the spatial domain. Here, the most common way of approximating the Laplacian is by finite-difference stencils, which read, for the 2D case and resolution h x > 0 and h y > 0, as x φ kÀ1;l À 2φ k;l þ φ kþ1;l À Á þ 1 h 2 y φ k;lÀ1 À 2φ k;l þ φ k;lþ1 À Á [21] with a straightforward adaptation for higher dimensions (86). Substitution of this into Equation [16] leads to the discrete Poisson problem for the approximate unwrapped phase φ apprx according to where the boundary condition also has to be understood in a discrete finite-difference sense. Again, on rectangular domains, Equation [22] is solvable with transform techniques that diagonalize the discrete Laplace operator. It turns out that that can again be achieved with the DCT of Type II, because of the homogeneous Neumann conditions. This leads to the solution formula Þþ4 sin 2 ω 2 π= 2M ð Þ ð Þ ! ; [23] where (N, M) denotes the size of the data and, as before, no division should be performed for ω = 0. It is worth noting that here, in contrast to the discrete version of Equation [18], the solution step for the discrete Poisson equation is numerically exact (87) and the only source of error is the fact that in general Δ h φ unwrap ≠ Δ h φ apprx . Nonetheless, the solution operator transforms local errors into global errors, meaning that, even in the case where the above inequality holds in only one pixel, we have φ unwrap ≠ φ apprx in every pixel; that is, the solution will in general be exact nowhere, as is also the case for classical Fourier approaches. Nevertheless, one can say that, when the discrete version of φ unwrap originates from a continuously defined smooth function, the error to the exact solution can be estimated accord- , where C > 0 is a constant within a reasonable range. As solving the discrete Poisson equation is stable with respect to the resolution (88), it also holds that φ unwrap k;l À φ apprx k;l ⩽C′ h 2 x þ h 2 y for some C ′ > 0 not depending on h x and h y , i.e., the pixelwise error can be controlled and vanishes with a quadratic rate in the resolution. The latter is, however, only of theoretical interest since in practice phase data cannot be measured with arbitrary resolution. Furthermore, the spatial approximation of Δ h φ unwrap makes it possible to obtain quantitatively exact results in some cases using Laplacian methods. The general strategy here is to replace the computation of Δ h φ apprx in Equation [22] with the aim that it matches or at least gives a better approximation to the unknown Δ h φ unwrap . An idea in this direction is founded on the observation that, in the continuous setting, one can already obtain the gradient of φ unwrap from φ wrap via ∇φ unwrap = À iψ À1 ∇ψ (see also Equation [19]). As before, this does not transfer to discrete data, however, suggesting computation of an approximation of the discrete exact gradient ∇ h φ unwrap from φ wrap . Indeed, assuming the lowest possible phase change between two pixels, for instance, one arrives at ∇ h φ apprx ð Þ k;l ¼ h À1 x arg ψ kþ1;l =ψ k;l À Á h À1 y arg ψ k;lþ1 =ψ k;l À Á ! [24] where arg denotes the argument of a complex number, i.e., the angle that lies here in the interval [ À π , π). If the phase change in φ unwrap is indeed in this interval, the discrete gradient will be exact. Consequently, applying the discrete divergence gives the discrete Laplacian, i.e. Δ h φ apprx = div h (∇ h φ apprx ), which would also be exact as well as the resulting φ apprx , as solving the discrete Poisson problem does not introduce errors. However, in turn, any deviation of ∇ h φ apprx according to Equation [22] from ∇ h φ unwrap is an integer multiple of 2π/h x or 2π/h y and causes, as before, a global deviation of φ apprx from φ unwrap , again even in the case where this occurs only in a single pixel.
As for the Fourier-type methods, the spatial variants consist of only few computational steps. Lipschitz continuity can again be ensured for Equation [23], but, in contrast, fails for Equation [24] due to the fact the complex argument is not continuous. Nevertheless, as argued above, employing a suitably adapted version of Equation [23] leads to exact results if ∇ h φ apprx is exact, and this might still be the case in the presence of noise.