yoyololicon/kamui

Kamui

Kamui is a python package for robust and accurate phase unwrapping on 2-D, 3-D, or sparse data.

Kamui unwrap the phases by viewing the data points as vertices $V$ connected with edges $E$ and solving the following integer linear programming (ILP) problem:

$$\min_{k} w^T |k|,$$ $$\text{s.t.} Ak = -A\frac{x}{2\pi},$$

where $k_{i \in [0, M)} \in \mathbb{Z}$ is the edge ambiguities to be computed, $w_{i \in [0, M)} \in \mathbb{R}^+$ is the weights, $x_{i \in [0, M)} = (V_v - V_u + \pi) \pmod {2\pi} - \pi | (u, v) = E_i$ is the pseudo phase derivatives, $M = |E|$. $A_{ij} \in \{-1, 0, 1\} | i \in [0, N) \cap j \in [0, M)$ and $N$ is the number of elementary cycles enclosed by $E$.

This formulation is based on the fact that the true phase differences, $2\pi k + x$, should fulfill the irrotationality constraint, which means the summation of phase derivatives of each elementary cycles is zero. This is the general form of the network programming approach proposed in the paper "A novel phase unwrapping method based on network programming".

Unwrapping phase with Kamui can be computationally heavy due to the fact that ILP is NP-hard. Acceleration techniques, such as dividing the graph into subgraphs, will be implemented in the future.

Installation

pip install kamui

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.

Usage

For regular 2-D or 3-D data such as interferograms, use kamui.unwrap_dimensional:

import numpy as np

def unwrap_dimensional(
    x: np.ndarray,
    start_pixel: Optional[Union[Tuple[int, int], Tuple[int, int, int]]] = None,
    use_edgelist: bool = False,
    **kwargs
) -> np.ndarray:
    """
    Unwrap the phase of a 2-D or 3-D array.

    Parameters
    ----------
    x : 2-D or 3-D np.ndarray
        The phase to be unwrapped.
    start_pixel : (2,) or (3,) tuple
        the reference pixel to start unwrapping.
        Default to (0, 0) for 2-D data and (0, 0, 0) for 3-D data.
    use_edgelist : bool
        Whether to use the edgelist method.
        Default to False.
    kwargs : dict
        Other arguments passed to `kamui.unwrap_arbitrary`.

    Returns
    -------
    np.ndarray
        The unwrapped phase of the same shape as x.
    """

For sparse data, use kamui.unwrap_arbitrary:

import numpy as np

def unwrap_arbitrary(
    psi: np.ndarray,
    edges: np.ndarray,
    simplices: Iterable[Iterable[int]] = None,
    method: str = "ilp",
    period: float = 2 * np.pi,
    start_i: int = 0,
    **kwargs,
) -> np.ndarray:
    """
    Unwrap the phase of arbitrary data.

    Parameters
    ----------
    psi : 1D np.ndarray of shape (P,)
        The phase (vertices) to be unwrapped.
    edges : 2-D np.ndarray of shape (M, 2)
        The edges of the graph.
    simplices : Iterable[Iterable[int]] of length (N,)
        Each element is a list of vertices that form a simplex (a.k.a elementary cycle).
        The connections should be consistent with the edges.
        This is also used to compute automatic weights for each edge.
        If not provided and method is "ilp", an edgelist-based ILP solver will be used without weighting.
    method : str
        The method to be used. Valid options are "ilp" and "gc", where "gc" correponds to PUMA.
        Default to "ilp".
    period : float
        The period of the phase.
        Default to 2 * np.pi.
    start_i : int
        The index of the reference vertex to start unwrapping.
        Default to 0.
    kwargs : dict
        Other arguments passed to the solver.

    Returns
    -------
    np.ndarray
        The unwrapped phase of the same shape as psi.
    """

Examples

• Unwrapping on 2-D synthetic data

TODO

• subgraph division
• edges-based custom weighting
• vertices-based custom weighting

References

• scikit-image/scikit-image/#4622
• my medium blogpost
• A novel phase unwrapping method based on network programming
• Phase Unwrapping via Graph Cuts
• Edgelist phase unwrapping algorithm for time series InSAR analysis
• Time Series Phase Unwrapping Based on Graph Theory and Compressed Sensing

kamui's Issues

Unwrapping of cyclical data

I would love to give Kamui a try in my particular application where I am tasked with unwrapping weighted gridded data with one cyclical axis. However, I do not understand the algorithm well enough to figure out how to implement cyclical axes. An example of how to unwrap 2D data with one cyclical axis would be highly appreciated. Ideally, this functionality would be part of the library itself with an optional cyclical_axis keyword.

Unwrapping for Non-Continuous Data

Hi @yoyololicon ! As you said yesterday, I opened an issue here.. might be better. Thanks for the fast answer, btw! I tried the masking of the NaNs before performing the unwrapping already, but then the issue is that the unwrapping works well in some places and not in others. I think it may have something to do with the starting point of the unwrapping... I think some weights are usually used in other methods to solve this, but I don't have this option here, right? Maybe I could use the edges or starting pixel parameter to try and fix this?
I am using wrapped cyclical data from -pi to pi and the unwrapping works well until there are too many NaNs dividing two areas.

I'm really trying to find a PUMA code that solves this, since right now I'm using SNAPHU but I think PUMA could be better.

Thanks!!

Bug in weights for ILP

Hey, I've noticed what seems to be a typo in line 104 in core.py.
The c which you use as argument to np.tile is not defined earlier.
I'd appreciate if you could let me know what the correct argument here is so that I can use weights in my computations.

Questions about the computing time and `unwrap_arbitrary` inputs

Thank you a lot for creating this open phase-unwrapping package! This is indeed what I am searching hard.

I have some questions about the implementation of

1. I think you implement the edgelist algorithm. As far as I know, this algorithm minimizes L1 as the mcf method. Why this algorithm is NP-hard while mcf is not? Does it mean the edgelist is much slower than mcf?
2. How should I generate the edges and simplices in unwrap_arbitrary?

Thanks!