Weighted Least-Squares (algo_ls_weighted)

This is a simplified weighted least-squares unwrap. It masks large wrapped gradients, then solves a weighted Poisson equation with Jacobi iterations. The code here uses binary weights and periodic neighbour rolls, so it matches the implementation in this repository rather than a full robust solver from the literature. See Ghiglia and Romero (1994) and Ghiglia and Pritt (1998).

API: xpunwrap.algo_ls_weighted

Derivation

1. Input: wrapped phase \(\phi_w\) in \([-\pi,\pi)\).

2. Wrapped gradients:

   \[g_x = \mathrm{wrap}(\phi_w(x+1,y) - \phi_w(x,y)), \quad g_y = \mathrm{wrap}(\phi_w(x,y+1) - \phi_w(x,y))\]

   In the code, forward differences are used and the last row/column is repeated before wrapping.

   The weighted divergence and Jacobi updates both use periodic neighbour access.

3. Gradient magnitude:

   \[m = \sqrt{g_x^2 + g_y^2}\]

4. Weights: Set \(w = 0\) where \(m\) exceeds a threshold, otherwise \(w=1\).

5. Weighted least-squares objective:

   \[\phi^\star = \arg\min_{\phi}\sum_{x,y} w\,\lvert \nabla \phi - g \rvert^2\]

6. Weighted divergence:

   \[f = \nabla \cdot (w \, g)\]

7. Weighted Poisson equation (normal equation):

   \[\nabla \cdot (w \nabla \phi) = f\]

   (This is a PDE, not an ODE.)

8. Jacobi iterations: The solver updates each pixel from its neighbours and local weights:

   \[\phi^{k+1} = \frac{ w_{x+}\phi_{x+} + w_{x-}\phi_{x-} + w_{y+}\phi_{y+} + w_{y-}\phi_{y-} - f }{w_{x+} + w_{x-} + w_{y+} + w_{y-}}\]

   The \(-f\) term above reflects the implementation’s discrete sign convention.

   The update uses periodic neighbour rolls in code.

9. Output: After a fixed number of iterations, \(\phi\) is returned. If the input was 2D, the leading singleton dimension is removed.

Note

Boundary conditions and padding. The Jacobi stencil uses xp.roll for all neighbour lookups, which wraps at domain edges (periodic boundary conditions). No zero-padding is applied. The border-weight mask (step 4) mitigates edge artifacts by setting weights to zero wherever the gradient magnitude exceeds border_thresh, effectively decoupling those pixels from the solve. This makes the algorithm more robust to non-periodic data than the pure FFT solvers, at the cost of Jacobi convergence speed (controlled by n_iter).

References

• D. C. Ghiglia and L. A. Romero, “Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods,” J. Opt. Soc. Am. A, vol. 11, no. 1, pp. 107-117, 1994.

• D. C. Ghiglia and M. D. Pritt, “Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software,” Wiley, 1998.