Three Dimensional phase unwrapping algorithms: a comparison.
Hussein S. Abdul-Rahman, Munther A. Gdeisat, David R. Burton, Michael J. Lalor
General Engineering Research Institute (GERI), Liverpool John Moores University, UK
, , ,
Abstract. In this paper, three three-dimensional phase unwrapping algorithms are discussed and explained. The first algorithm is the one proposed by Cusack in 2002 to unwrap MRI images. The second algorithm proposed by Huntley in 2001. The final one is the one proposed by the authors in 2005. The major steps of each algorithm is discussed and explained in this paper. A computer generated objects were developed to compare between these algorithms. Finally the above three algorithms were compared with a well known two dimensional phase unwrapping algorithms.
1. Introduction
Phase unwrapping is one of the most important processing step in many advanced imaging technologies such as optical interferometry, satellite radar interferometry (SAR), and magnetic resonance imaging (MRI) where the required data is encoded in the form of a phase distribution. The measure phase is normally wrapped between – π and π. The wrapped phase is not usable until the 2π phase discontinuities are resolved using a phase unwrapping algorithm.
During the last three decades, the field of two-dimensional phase unwrapping has been intensively studied and plenty of journal papers have been published and the field has its own textbook1. Nowadays, three-dimensional phase unwrapping algorithms have increasing interests to unwrap a wrapped phase volumes. Three-dimensional phase unwrapping field is still new area of research and only a few algorithms have been proposed.
Huntley in 2001 proposed a three dimensional algorithm which extends the two-dimensional residue-balancing method into three dimensions. In this method; all residues in the phase volume are identified and connected together to form singularity loops, then these loops are set as prohibited areas during unwrapping paths2.
Cusack et al in 2002 proposed another three-dimensional phase unwrapping algorithm, which uses a quality measurement depending on calculating the residue field. This algorithm used many iteration to unwrap the whole volume, where in each iteration only those voxels (a single element in the 3D volume) whose quality exceed a certain threshold are unwrapped and other are left to the next iteration where the threshold is lower3.
The authors in 2005 proposed another quality-guided three-dimensional phase unwrapping algorithm called three-dimensional sorting reliabilities algorithm (3D-SRA). This algorithm depends on finding the quality of each edge, which connect two voxels together in x, y or z direction. The unwrapping path is then carried out from high quality edges to low quality edges4.
The above three algorithms follow different strategies to unwrap phase volumes, Cusack’s algorithm as well as 3D-SRA depends on a quality measure to guide the unwrapping. However, both of them follow different approaches to construct the unwrapping path as will be explained in subsections 2.1 and 2.2 later. Huntley’s algorithm on the other hand, unwrap the phase volumes without relying on a quality measure as will be explained in subsection 2.3 in this paper. The results, discussions and comparisons are shown in section 3.
2. Algorithms
2.1. Cusack’s Algorithm
Cusack et al present a new algorithm for robust three-dimensional phase unwrapping, in which unwrapping is guided, so that it initially works on less noisy regions3.
Fig. 1 below outlines the main steps in this algorithm. As shown in the figure, this algorithm calculates the amount of noise in each voxel and unwrap the less noisy parts first by assigning a noise threshold that allows certain voxels to be unwrapped and leave the other voxels to next iterations.
The following subsections describe each stage of this algorithm separately.
Figure 1 Flow chart for Cusack’s algorithm
2.1.1. Calculating noise estimation field.
This algorithm adopted two main methods to calculate the noise estimation field in order to guide the unwrapping procedure. The first method was based on the magnitude of the signal in the phase map. The second estimation was derived from the phase singularities (poles). Phase poles can be easily identified by unwrapping around 2x2 voxels in x,y and z directions.
To identify poles in z direction, we have to unwrap around 2x2 loop in xy plane with coordinates: (i,j,k), (i,j+1,k) , (i+1,j+1,k), (i+1,j,k) and (i,j,k) using Eq. 1 below. A pole is encountered if the result of the equation is either +1 or -1.If we find a pole in this square, then this pole is positioned at (i+0.5, j+0.5,k). Similar calculation is carried out in both y and x directions using Eqs. 2 and 3 respectively.
After identifying all poles in the phase data volume, the noise estimation field values is calculated. The noise estimator is assigned so that it is high in the voxels close to the pole and lower in the voxels further away.
Where the operator R[.] rounds its argument to the nearest intger.
2.1.2. Setting the noise threshold.
To control the threshold, the minimum and maximum noise estimator vales were found and the number of steps is assigned. First, the threshold is set to be small to unwrap only the good voxels. Then the threshold is increased by single step and the unwrapping is repeated again and again until the whole volume is unwrapped.
2.1.3. Flood-fill unwrapping.
In this algorithm a threshold-guided flood fill algorithm is used to unwrap the whole volume. The guided flood fill phase unwrapping can be summarized as follows:
1- Select a starting voxel C.
2- for each voxel of the six neighbors of the voxel C:
a. If this neighbor has a noise field estimator less than or equal the current threshold, unwrap it and repeat step 2.
b. If this neighbor has a noise field estimator greater than the current threshold leave it and repeat step 2 for the remaining neighbors.
3- When all voxels encountered below the current threshold are completed, increase noise field estimator and repeat step 2 until all voxels are unwrapped.
2.2. Three-dimensional sorting reliabilities algorithm (3D-SRA)
Three-dimensional sorting reliabilities algorithm unwraps the phase data volume beginning with high quality voxels first leaving the bad quality voxels to the last. Unlike Cusack’s algorithm, three-dimensional sorting reliabilities algorithm unwraps the phase volume following a discrete path. This discrete path assures to select the best path.
The algorithm can be outlined in the following steps:
· Determine the reliability of all voxels except the voxels at the border surfaces.
· Calculate the horizontal, vertical and normal edges’ reliability.
· Sort all of these edges according to their reliabilities.
· Unwrap voxels according to the edges reliabilities; so that the voxels which connect the highest reliable edges are unwrapped first.
· Unwrap the border surfaces.
2.2.1. Determining voxels’ reliabilities.
Many methods can be applied to determine the reliability (quality) of each voxel individually. Determining the quality of each voxel is carried out by calculating the quality map to the input data volume. Quality maps can define to be arrays of values that define the quality or goodness of each voxel of the given phase data. Different algorithms already exist to construct a quality map such as; phase derivative variance, maximum phase gradient and second phase difference and many other algorithms with can be applied in both two-dimensional and three-dimensional phase unwrapping algorithms1. Three-dimensional sorting reliabilities algorithm can use and of the mentioned three-dimensional quality maps to guide the unwrapping process.
2.2.2. Determining edges’ reliabilities.
An edge is an intersection of two voxels that are connected in the x, y or z directions. Edges can be classified as vertical, horizontal and normal edges as shown in Fig. 2. The reliability of an edge is defined as the summation of the reliability of the two voxels that the edge connects as shown in Fig. 2. The reliability value of an edge that connects a border voxel with another voxel in the phase volume is set to zero.
Figure 2 The Horizotal, vertical and Normal Edges and their reliabilities.
2.2.3. Unwrapping Path.
An unwrapping path cannot be defined relative to the reliability of the voxels. Instead it is defined by looking at the value of the reliability of the edges. The definition of the unwrapping path is relatively simple. All the edges are stored in an array and sorted by value of reliability. Those edges with higher reliability are resolved first as will be explained below.
Figs. 3(a) - 3(d) illustrate the principle of the algorithm. Suppose that we want to unwrap the voxels, whose edges reliability are already sorted and their orders are shown in Fig 3(a). The integer numbers represent the order of edges in the sorted array. The edge that has the order 1 is processed first then the edge with order 2 and so on.
Initially all voxels are considered as not belonging to any group, which is represented by the white circles in Fig. 3(a). Both voxels surrounding the edge with order 1 in Fig. 3(a) are unwrapped first with respect to each other, and both voxels are joined together in single group as illustrated in Fig 3(a).
The voxels forming the edge 6 in Fig. 3(b) are to be processed in this stage. As we can see one of these voxels has its own group but the other pixel does not belong to any group. The ungrouped voxel is unwrapped with respect to the grouped one, and set to be a member of the group as shown in Fig. 3(c).
Edge number 7 in Fig. 3(c) connects two grouped voxels, so in this case the group of fewer members is unwrapped with respect to the group with more members, and both groups merge into a single group as illustrated in Fig. 3(d). Note that unwrapping a voxel or a group of voxels with respect to another group may require the addition or subtraction of multiples of 2p.
From the previous explanation, the unwrapping path can be summarised as follows:
1- If both voxels do not belong to any group and have not been unwrapped before; the voxels are unwrapped with respect to each other and gathered into a single group.
2- If one of the voxels has been processed before and belongs to a group but the other has not, then the voxel that has not been processed before is unwrapped with respect to the other voxels in the group and now joins this group.
3- If both voxels have been processed before and both belong to different groups, then the two groups are unwrapped with respect to each other. The smallest group is unwrapped with respect to the largest group. Then the two groups are joined together to construct a single group.
2.3. Huntley’s algorithm
Huntley extends the classical problem of phase unwrapping in two dimensions, how to create a path-independent unwrapped map, to the case of three-dimensional phase distribution2. In two dimensions the path dependence problem arises from isolated phase singularity points, in three dimensions the phase singularities forms closed loops in space. In two dimensions, path independence is achieved when branch-cut lines are placed between singular points of opposite sing. In three dimensions, path independence is achieved by placing surface-cuts on every phase singularity loop to prevent unwrapping through these loops. Unlike the previous two algorithms, Huntley’s algorithms does not rely on a quality measure to guide the unwrapping path, Instead this algorithm identify singularity loops in the phase volume to avoid passing through them in the unwrapping.
Huntley’s algorithm can be outlined in following steps:
1- Identifying singularity points (residues) in phase volume in all dimensions.
2- Connecting residues with each other to form singularity loops and placing surface-cuts on every single loop.
3- Unwrapping the phase volume without cutting the surface-cuts.
2.3.1. Identifying singularity points.
Singularity points or residues in Huntley’s algorithm are equivalent to the poles in Cusack’s algorithm and they can be extracted using Eqs. 1, 2 and 3 to identify all residue types in x, y and z direction. Here, residue is a vector quantity that has a direction. Figure 4 shows the direction of positive sz, sy and sx. Note that all of shown residues entering the cubes shown in the figure.
Figure 3 Unwrapping path for the three-dimensional sorting reliabilities algorithm.
Figure 4 The Direction of positive residues (a) positive Sz (b) positive Sy and (c) positive Sx.
2.3.2. Identifying singularity loops.
Identifying singularity loops in the phase data volume depends on two facts outlined and discussed in [2]. The first fact is that all residues must form a close loop in volume space or if not must terminate at the volume border. The second fact is that the number of residues entering any cube in the volume space must equal to the number of residues leaving that cube. Depending on these facts, identifying the singularity loops is kind of searching algorithm that connect appropriate residues with each other. Fig. 5 shows a part of a loop where each residue enters any and of the cubes is connected to other residue that leaves the same cube.
After Identify singularity loops, the entire surfaces bounded with the loops itself are set as cut-surfaces that must be avoided during the unwrapping stage.
Figure 5 Forming a singularity loop.
2.3.3. Unwrapping around surface-cuts.
Unwrapping the whole volume is carried out using the three-dimensional flood-fill algorithm described earlier in section 2.1.3. The only difference is that, instead of depending on the noise threshold to unwrap or leave a given voxel, we leave in this case those voxel cases passing through the surface-cuts and unwrap all other voxels. This stage of the algorithm is explained in details in [2].