Introduction

Phase unwrapping is a problem that occurs in several fields as diverse as Synthetic Aperture Radar and MR Angiography. In all cases the problem is that the measured phase signal can only take on values in a $2\pi$ range, whilst the original phase signal can take on any value. That is, $\theta_{measured} = \theta_{original} \mathrm{mod} 2\pi$. This ``wrapping'' leads to artificial phase jumps being introduced near boundaries. To recover the original, smooth phase (up to some global $2\pi$ offset) requires an unwrapping algorithm.

Medical MR images require phase unwrapping for many applications such as field mapping (often used to remove geometric distortion in EPI [14,10,20]), chemical shift mapping [6] and velocity measurement [15] (used in MR angiography). Also, in medical imaging applications the phase maps are typically 3D, rather than the 2D maps acquired in applications such as optical and radar imaging. It is therefore important to have a method which takes into account the 3D nature and constraints of the problem. A 2D method may be applied on a slice-by-slice basis, but can result in phase wraps existing in the through-slice direction after unwrapping, which must be fixed by some other post-processing option. Furthermore, 3D images usually contain many more voxels than typical 2D images, and so speed is an issue, especially in applications such as automatic-shimming where rapid, automatic and robust methods for phase unwrapping are essential.

Furthermore, in 3D MRI, a large proportion of the voxels do not contain any interesting signal (e.g. air, bone, etc.) and so the signal to noise ratio (SNR) in these voxels is extremely low. Consequently the true phase information is swamped by noise. By disregarding such voxels when unwrapping, both robustness and speed can be greatly improved.

Existing unwrapping techniques can be categorised according to their:

• Dimensionality (1D, 2D, 3D, etc.);
• Application (Synthetic Aperture Radar, general optical interferometry, MR angiography, MR chemical shift mapping or MR field mapping);
• Approach (fitting functions, cost function optimisation, filtering, region growing/merging).

Methods based on fitting functions (usually truncated Taylor series [13,5] or polynomials [10,6]) can be easily generalised to work with data of any dimension. However, these methods impose considerable smoothness on the phase since the functions can not vary too rapidly, as otherwise no unwrapping will be performed. This is undesirable for mapping applications where there are small, rapidly changing regions that are not smooth and hence not estimated correctly.

Generalising other methods can be difficult because they use application specific knowledge, such as the coherence function in Synthetic Aperture Radar (SAR) imaging [4,11,12] or multiple MR images for chemical shift mapping [19] or prior expectations about expected velocities [9,21,3]. In addition difficulties can arise because the underlying method does not generalise naturally, such as for the filtering methods (e.g. 1D Markov [7] or 2D Kalman [12]).

Many of the existing unwrapping techniques proposed for MR applications were only written to deal with 2D slices. These methods are based on various approaches including region growing [16,18,1], cost function optimisation [17], neighbourhood voting/filtering [2] and phase path integration [8]. Although these approaches all generalise straightforwardly to 3D images they are generally inefficient and consequently unacceptably slow for many applications, especially rapid, automated shimming.

This report presents a fast, robust method for unwrapping phase maps that is general and can be applied to a single ``image'' of any dimension. The method uses a region-merging approach to find the minimum of a cost function that penalises phase differences across boundaries. It has been implemented specifically for 2D and 3D medical images (and is available at www.fmrib.ox.ac.uk/fsl), but the principles apply equally well to phase unwrapping problems of any dimension. In fact, the existing implementation would only require minimal modification to tackle any particular phase unwrapping problem.