Cost Function Formulation

The approach taken here to solve this problem starts by dividing the whole (n-dimensional) volume into $N$ regions, each of which contains no phase wrappings. Each region is then treated as a single entity by the algorithm, with a single offset. As the number of regions to be dealt with is much less than the number of voxels this provides a considerable speed increase.

Now consider the interface, or boundary, between two regions, $A$ and $B$.

The unwrapped phase, $\theta$ in region $A$ will be related to the wrapped phase, $\phi$ by: $\theta_{Aj} = \phi_{Aj} + 2\pi M_A$ where $M_A$ is an integer offset that is the same for all voxels in region $A$.

A cost function suitable for unwrapping is one that penalises phase differences along this interface. One good candidate is the sum of squared difference in phase. That is:

$$C_{AB} = \sum_{j, k \in {\cal N}(j)} \left( \theta_{Aj} - \theta_{Bk} \right)^2$$ (2)

$$= \sum_{j,k \in {\cal N}(j)} \left( \phi_{Aj} - \phi_{Bk} + 2 \pi M_{AB} \right)^2$$ (3)

where $M_{AB} = M_A - M_B$ and $M_A$ is the integer offset applied to region $A$. The summation is taken over two indices, $j$ and $k$, where $j$ is the index of a voxel in region $A$, while $k$ is the index of a voxel in region $B$ that is also in the (simply connected) neighbourhood of the voxel with index $j$ (i.e. $j \in {\cal I}(A)$ and $k \in {\cal N}(j) \cap {\cal I}(B)$).

The total cost over the whole volume is the sum over all the interfaces:

$$C = \sum_{A,B} C_{AB}.$$ (4)

Differentiating with respect to the individual parameters, $M_A$, gives

$$\frac{\partial C_{AB}}{\partial M_A} = 8 \pi^2 N_{AB} M_{AB} + 4 \pi P_{AB} = 0$$ (5)

which yields the equations for the minimum cost solution:

$$M_{AB} = M_A - M_B = \frac{- P_{AB}}{2 \pi N_{AB}}.$$ (6)

where $N_{AB}$ is the number of interfacing voxel pairs and $P_{AB} = \sum_{j,k \in {\cal N}(j)} ( \phi_{Aj} - \phi_{Bk} )$.

However, not only does this tend to generate more equations (equal to the number of distinct interfaces) but, more importantly, only integer values can be taken by the parameters $M_A$. Therefore, equation 6 cannot be solved exactly as it is an integer programming problem. Finding the optimal solution (minimum cost) for such a problem is a difficult and often time-consuming task. For example, there are $2^N$ ``neighbouring'' solutions of the continuous solution (where $N$ is the number of regions).