Kamui is a python package for robust and accurate phase unwrapping on 2-D, 3-D, or sparse data.
Kamui views the data as a graph
The unwrapped phase
Here:
-
$M = \lvert E \rvert$ is the number of edges, and$w \in \mathbb{R}_{\ge 0}^{M}$ are the edge weights. -
$A \in \{-1, 0, 1\}^{N \times M}$ is the incidence matrix of the$N$ elementary cycles.$A_{ce} = 1$ if cycle$c$ traverses edge$e$ from$u$ to$v$ ,$-1$ if it traverses it from$v$ to$u$ , and$0$ if$e$ is not on$c$ . - The right-hand side is integral:
$\frac{1}{2\pi} (A x)_c$ is the residue of cycle$c$ , the number of whole turns that the wrapped differences around it add up to.
Kamui then integrates period=T, replace
On 2-D grids and planar meshes, every edge lies on at most two cycles. Kamui orients the cycles so that they traverse each shared edge in opposite directions, reversing any that do not; reversing a cycle negates both sides of its constraint, so the program does not change.
Otherwise, as on 3-D grids or with fractional weights, Kamui first solves the linear programming (LP) relaxation, with
pip install kamui
This also installs pylmcf, Python bindings for LEMON's network simplex under the Boost Software License, which Kamui uses wherever it applies.
Kamui also provides PUMA, a fast and robust phase unwrapping algorithm based on graph cuts as an alternative. To install PUMA, run
pip install kamui[extra]
However, it uses the original maxflow implementation by Vladimir Kolmogorov with GPL license. Please follow the licensing instruction in PyMaxflow if you use this version of Kamui.
The examples below build on each other and run as written.
import numpy as np
import kamui
# a smooth synthetic phase surface, wrapped into [-pi, pi)
y, x = np.mgrid[0:128, 0:128]
true_phase = 0.002 * (x - 40) ** 2 + 0.001 * (y - 70) ** 2 + 0.05 * x
wrapped = kamui.wrap_difference(true_phase)
unwrapped = kamui.unwrap_dimensional(wrapped)
# the result equals the true phase up to a constant multiple of 2 pi
assert np.allclose(unwrapped - unwrapped[0, 0], true_phase - true_phase[0, 0])The reference pixel keeps its wrapped value. It defaults to the first pixel; choose another with start_pixel, for example start_pixel=(64, 64).
Per-pixel quality weights, such as InSAR coherence, tell the solver where corrections are cheap. Each edge gets the sum of its two pixels' weights; merging_method="min" or "max" change that. Edges that touch a NaN weight get weight 0. Kamui uses the weights as given: multiplying them all by one factor does not change the result.
LEMON needs integer weights, and Kamui does not round them for you. Fractional weights are solved by HiGHS, which is several times slower. To use LEMON, scale and round the weights yourself; a larger factor keeps more resolution:
coherence = np.random.default_rng(0).uniform(0.2, 1.0, wrapped.shape)
unwrapped = kamui.unwrap_dimensional(wrapped, weights=np.round(coherence * 100))With solver="lemon", fractional weights raise an error instead of falling back to HiGHS.
unwrap_arbitrary takes per-edge weights and uses them exactly as given. Up to version 0.2, unwrap_dimensional rescaled the per-pixel weights linearly to kamui.prepare_weights still does that, so pass its result to unwrap_arbitrary to keep that behaviour:
edges, cycles = kamui.get_2d_edges_and_simplices(wrapped.shape)
edge_weights = kamui.prepare_weights(coherence, edges) # rescaled, then averaged per edge
unwrapped = kamui.unwrap_arbitrary(wrapped.ravel(), edges, cycles, weights=edge_weights)
unwrapped = unwrapped.reshape(wrapped.shape)For data that wraps around along an axis, such as the angle of a polar grid, mark that axis as cyclical so the first and last slices are connected:
theta = np.linspace(0, 2 * np.pi, 64, endpoint=False)
radius = np.linspace(1, 3, 48)
true_ring = 6 * np.sin(theta)[:, None] + 2.5 * radius[None, :]
unwrapped = kamui.unwrap_dimensional(kamui.wrap_difference(true_ring), cyclical_axis=0)
assert np.allclose(unwrapped - unwrapped[0, 0], true_ring - true_ring[0, 0])z, y, x = np.mgrid[0:16, 0:16, 0:16]
true_volume = 0.4 * x + 0.3 * y - 0.35 * z + 0.01 * (x - 8) ** 2
unwrapped = kamui.unwrap_dimensional(kamui.wrap_difference(true_volume))
assert np.allclose(unwrapped - unwrapped[0, 0, 0], true_volume - true_volume[0, 0, 0])For scattered points, use unwrap_arbitrary with the phase at each point, the edges between points and, optionally, the elementary cycles ("simplices"), such as the triangles of a Delaunay triangulation:
from scipy.spatial import Delaunay
rng = np.random.default_rng(0)
points = rng.uniform(0, 10, (2000, 2))
point_phase = 1.5 * points[:, 0] - points[:, 1]
psi = kamui.wrap_difference(point_phase)
triangles = Delaunay(points).simplices
# Delaunay adds long, thin triangles along the convex hull. Drop them, so that
# no edge is too long for its phase change to be resolved.
sides = triangles[:, [[0, 1], [1, 2], [2, 0]]]
lengths = np.linalg.norm(points[sides[..., 0]] - points[sides[..., 1]], axis=-1)
keep = lengths.max(axis=1) < 1.0
triangles = triangles[keep]
edges = np.unique(np.sort(sides[keep].reshape(-1, 2), axis=1), axis=0)
unwrapped = kamui.unwrap_arbitrary(psi, edges, triangles)
assert np.allclose(unwrapped - unwrapped[0], point_phase - point_phase[0])List each edge once, in one direction, as the np.unique(np.sort(...)) line above does. The simplex ILP raises ValueError if an edge appears twice, in the same or the opposite direction.
Without simplices, unwrap_arbitrary(psi, edges) uses the edgelist formulation, which needs only edges. That suits graphs without a natural cycle structure, such as interferogram networks in time series analysis.
| Solver | How to select it | Needs | Notes |
|---|---|---|---|
| Simplex ILP (default) |
unwrap_dimensional(x), or unwrap_arbitrary(psi, edges, simplices)
|
edges and elementary cycles | By default each edge costs the number of its cycles that have no residue; weights overrides that. 2-D grids and planar meshes with integer weights, the defaults included, are solved by LEMON, and everything else by HiGHS. solver="highs" forces HiGHS, and solver="lemon" raises an error instead of falling back. |
| Edgelist ILP |
use_edgelist=True, or unwrap_arbitrary(psi, edges)
|
edges only | Uniform weights unless weights is given. Integer weights, the default included, are solved by LEMON, and fractional ones by HiGHS; solver works as for the simplex ILP. |
| PUMA | method="gc" |
edges and pip install kamui[extra]
|
Minimizes a p, max_jump, and optional non-negative weights). Uses graph cuts; slower than the single-pass ILP/LEMON solvers on grids (#26). GPL; see Installation. |
The edgelist ILP needs no cycles, because it optimizes vertex offsets
The term inside the absolute value is the edge ambiguity
Its LP dual is a min-cost circulation on the graph itself, with a flow
PUMA instead uses graph cuts to minimize the weighted
Pass return_info=True to also get the solver's report, a scipy.optimize.OptimizeResult. It is returned even when the result is None, so info.message says why the solve failed.
unwrapped, info = kamui.unwrap_dimensional(wrapped, return_info=True)
print(info.success, info.fun, info.ilp_fallback)funis the weighted L1 cost of the corrections; for PUMA, it is the final energy.success,statusandmessagecome from the solver.solvernames the solver that ran:"lemon"or"highs".ilp_fallbacktells whether the LP optimum was fractional, so HiGHS solved the integer program instead.
With the same weights, the simplex and edgelist solvers minimize the same cost, so their fun values can be compared directly:
noisy = kamui.wrap_difference(true_phase + np.random.default_rng(1).normal(0, 1.0, wrapped.shape))
edges, cycles = kamui.get_2d_edges_and_simplices(noisy.shape)
_, simplex = kamui.unwrap_arbitrary(
noisy.ravel(), edges, cycles, weights=np.ones(len(edges)), return_info=True
)
_, edgelist = kamui.unwrap_arbitrary(noisy.ravel(), edges, return_info=True)
assert np.isclose(simplex.fun, edgelist.fun)unwrap_dimensional treats NaN pixels as missing data. It drops every edge and cycle that touches one, and the result stays NaN there. To mask out pixels, such as those with low coherence, set them to NaN.
Each connected component of the graph is unwrapped on its own, whether NaNs split the grid or the graph passed to unwrap_arbitrary is disconnected. The component that contains start_pixel (or start_i) is referenced to it, and every other component to its lowest-numbered vertex, so the offsets between components are arbitrary multiples of the period. A pixel with no valid neighbour keeps its wrapped value:
nan_phase = wrapped.copy()
nan_phase[40:60, :] = np.nan
unwrapped_nan = kamui.unwrap_dimensional(nan_phase)
assert np.isnan(unwrapped_nan[50, 0])
assert not np.isnan(unwrapped_nan[10, 0])unwrap_dimensional with default settings on noisy 2-D grids, measured on an Apple M1 Pro with SciPy 1.18 (HiGHS 1.12) and pylmcf 1.3.0. Memory is the peak footprint, which includes compressed memory:
| Grid | HiGHS | LEMON |
|---|---|---|
| 300×300 | 1.4 s, 0.3 GB | 0.09 s, 0.1 GB |
| 600×600 | 9.5 s, 1.3 GB | 0.52 s, 0.2 GB |
| 1000×1000 | 53 s, 3.6 GB | 1.7 s, 0.6 GB |
| 2000×2000 | — | 11 s, 3.2 GB |
| 4628×2562 | — | 49 s, 9.6 GB |
Both solvers reach the same optimal cost. The last row is the size of the interferograms in #11 and #12. LEMON's own solve takes most of the time and about half of the memory there.
The edgelist path (use_edgelist=True) and grids with integer weights also gain from LEMON: at 600×600, from 78 s to 18 s and from 12 s to 2.3 s.
If unwrapping is unexpectedly slow, check your SciPy build. On macOS (Apple Silicon), conda-forge's SciPy 1.15.0–1.15.2 builds solve kamui's programs hundreds of times slower than SciPy's PyPI wheels and other conda-forge releases, for example 5 s instead of 10 ms for a 32×32 grid.