pyoptics.unwrap

Phase unwrapping in one and two dimensions.

Measured or simulated phases (np.angle() of a field, the output of a phase shifting algorithm, a reconstructed hologram) are only known modulo 2π: they are wrapped into (-π, π]. Unwrapping adds the multiples of 2π back, under the assumption that the true phase changes by less than π between neighbouring pixels.

Two families of 2D algorithms are implemented:

• unwrap_ls(): least squares unwrapping (Ghiglia and Romero 1994). The unwrapped phase φ is the function whose discrete gradient is closest, in a (weighted) least squares sense, to the wrapped differences of the data. The unweighted problem is a Poisson equation with Neumann boundary conditions that a single DCT solves. Weights (e.g. a modulation or amplitude map) and masks are handled by preconditioned conjugate gradients, with the DCT solver as preconditioner. Least squares solutions are smooth and robust, but in general not congruent to the data: wrap(φ) differs from the input wherever the data contain residues. congruent=True snaps the result to the nearest congruent phase.

• unwrap_quality_guided(): path following in the order of a reliability (quality) map (Herráez et al. 2002). The best pixel pairs are joined first, so errors are confined to unreliable regions. The result is congruent by construction.

Independent information about the phase avoids path following altogether:

• unwrap_with_reference(): snaps the wrapped phase to the nearest congruent phase of a coarse but unwrapped estimate, e.g. a phase from the transport of intensity equation (pyoptics.TIE), a model or a measurement at a longer wavelength. Every pixel is unwrapped on its own, so residues and disconnected regions do not matter; the estimate has to be correct to better than ±π.

• unwrap_two_wavelength(): two-wavelength (synthetic wavelength) unwrapping of phases measured at two wavelengths, e.g. in dual-wavelength holography or interferometry.

residues() locates the phase singularities (vortices, noise points) that make unwrapping path dependent, reliability() computes Herráez’s second difference quality map, and unwrap_1d() wraps numpy.unwrap.

Conventions

Arrays are indexed [y, x] as from np.meshgrid() / utils.grid2d(); the x direction is the last axis. Pixels outside a mask (or NaN in the input) are returned as NaN. The piston of the unwrapped phase is arbitrary; all functions fix it such that the result agrees with the wrapped input modulo 2π on average (and exactly, for congruent results), with the mean value of the result in (-π, π].

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 11, 107 (1994).

J. Gass, A. Dakoff and M. K. Kim, “Phase imaging without 2π ambiguity by multiwavelength digital holography”, Opt. Lett. 28, 1141 (2003).

M. A. Herráez, D. R. Burton, M. J. Lalor and M. A. Gdeisat, “Fast two-dimensional phase-unwrapping algorithm based on sorting by reliability following a noncontinuous path”, Appl. Opt. 41, 7437 (2002).

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

Functions

reliability	Reliability (quality) map after Herráez et al. (2002).
residues	Residue map (phase singularities) of a wrapped phase.
synthetic_wavelength	Synthetic (beat) wavelength λ1·λ2/|λ1 − λ2| of two wavelengths.
unwrap	Unwrap a 1D or 2D phase.
unwrap_1d	Unwrap along one axis (numpy.unwrap()).
unwrap_ls	Least squares phase unwrapping (Ghiglia and Romero 1994).
unwrap_quality_guided	Quality guided unwrapping by reliability sorting (Herráez et al. 2002).
unwrap_two_wavelength	Optical path difference from wrapped phases at two wavelengths.
unwrap_with_reference	Unwrap a phase with an unwrapped estimate of it.
wrap	Wrap a phase into [-π, π).

pyoptics.unwrap.wrap(phi)

Wrap a phase into [-π, π).

Parameters:

phiarray or number

Phase in radians.

Returns:

phi_wrappedarray or number

phi minus the multiple of 2π that brings it into [-π, π). NaNs stay NaN.

pyoptics.unwrap.unwrap_1d(phase_wrapped, axis=-1, discont=3.141592653589793)

Unwrap along one axis (numpy.unwrap()).

Parameters:

phase_wrappedarray

Wrapped phase.

axisint, optional

Axis along which to unwrap.

discontfloat, optional

Maximum jump between neighbouring samples that is considered physical.

Returns:

phasearray

Unwrapped phase; the first sample along axis is left unchanged.

pyoptics.unwrap.residues(phase_wrapped)

Residue map (phase singularities) of a wrapped phase.

The residue of an elementary loop of 2×2 pixels is the sum of the wrapped phase differences around it, counterclockwise (in x, then y; i.e. with y pointing up), divided by 2π. It is an integer: 0 for a smooth phase, +1 at the core of a vortex exp(+iφ) (φ the azimuth), -1 for exp(-iφ). Wherever residues exist, unwrapping depends on the path, and a congruent unwrapped phase has to contain cuts between residues of opposite sign.

Parameters:

phase_wrappedarray (Ny, Nx)

Wrapped phase.

Returns:

resinteger array (Ny-1, Nx-1)

res[i, j] belongs to the loop with corners [i, j], [i, j+1], [i+1, j+1], [i+1, j], i.e. it sits at the point between those four pixels. Loops that touch a NaN are 0.

pyoptics.unwrap.reliability(phase_wrapped, mask=None)

Reliability (quality) map after Herráez et al. (2002).

The reliability of a pixel is 1/D with the second difference

D² = H² + V² + D₁² + D₂²,

H = wrap(ψ[i, j-1] − ψ[i, j]) − wrap(ψ[i, j] − ψ[i, j+1]) and likewise V (along y) and D₁, D₂ (along the diagonals). Second differences are small for smooth (low curvature) phases and large for noise and near residues.

Parameters:

phase_wrappedarray (Ny, Nx)

Wrapped phase.

maskboolean array, optional

Valid pixels.

Returns:

Rarray

Reliability; border pixels, pixels outside the mask and pixels with an invalid neighbour have R = 0.

pyoptics.unwrap.unwrap_ls(phase_wrapped, weights=None, mask=None, method='auto', congruent=False, tol=1e-10, max_iter=500, return_info=False)

Least squares phase unwrapping (Ghiglia and Romero 1994).

Minimizes

Σ w_x (φ[i, j+1] − φ[i, j] − Δ_x ψ)² + Σ w_y (φ[i+1, j] − φ[i, j] − Δ_y ψ)²

with the wrapped differences Δψ = wrap(ψ[neighbour] − ψ) of the data ψ. If the true phase changes by less than π between neighbours, the wrapped differences are the true ones and the solution is exact (up to a piston). Edge weights are the smaller of the two pixel weights squared, w_x = min(w[i, j]², w[i, j+1]²), as in Ghiglia and Pritt (1998).

Parameters:

phase_wrappedarray (Ny, Nx)

Wrapped phase. NaNs are treated as masked pixels.

weightsarray (Ny, Nx), optional

Pixel weights ≥ 0, e.g. a modulation, amplitude or quality map (normalized or not). Pixels with weight 0 are excluded.

maskboolean array (Ny, Nx), optional

True for valid pixels.

methodstr, optional

“dct”: unweighted least squares, solved directly by one DCT (only without weights and mask). “pcg”: weighted least squares by preconditioned conjugate gradients, preconditioned with the DCT solver. “auto” (default): “dct” if there are no weights, mask or NaNs, otherwise “pcg”.

congruentbool, optional

If True, return the congruent phase closest to the least squares solution, ψ + 2π·round((φ_LS − ψ)/2π). It is exact wherever φ_LS is within π of the truth, and keeps the full detail (and noise) of the data.

tolfloat, optional

Relative residual at which “pcg” stops.

max_iterint, optional

Maximum number of “pcg” iterations.

return_infobool, optional

If True, also return a dict with the number of iterations and the final relative residual (“pcg”; 0 and 0.0 for “dct”).

Returns:

phasearray (Ny, Nx)

Unwrapped phase, NaN outside the valid region. The piston is chosen such that phase agrees with phase_wrapped modulo 2π on average, and the mean of phase is in (-π, π].

infodict

Only if return_info is True.

pyoptics.unwrap.unwrap_quality_guided(phase_wrapped, quality=None, mask=None)

Quality guided unwrapping by reliability sorting (Herráez et al. 2002).

Every pair of horizontally or vertically neighbouring pixels (an edge) gets the reliability q_1 + q_2 of its pixels. Edges are processed from the most to the least reliable one; each joins the groups of its two pixels, adding to the smaller group the multiple of 2π that makes the step across the edge smaller than π. Unreliable regions (noise, residues) are thus unwrapped last and cannot spoil the good parts. The path is not continuous, which makes the algorithm fast (O(N log N) for the sort).

Parameters:

phase_wrappedarray (Ny, Nx)

Wrapped phase. NaNs are treated as masked pixels.

qualityarray (Ny, Nx), optional

Pixel quality, larger is better (e.g. modulation or amplitude). The default is reliability(), the inverse second difference.

maskboolean array (Ny, Nx), optional

True for valid pixels.

Returns:

phasearray (Ny, Nx)

Unwrapped phase, congruent with phase_wrapped (phase − phase_wrapped is a multiple of 2π everywhere), NaN outside the valid region. Regions that are not connected to each other get independent pistons. The mean of the result is in (-π, π].

pyoptics.unwrap.unwrap_with_reference(phase_wrapped, reference, mask=None)

Unwrap a phase with an unwrapped estimate of it.

Each pixel gets the multiple of 2π that brings it closest to the reference,

φ = ψ + 2π·round((reference − ψ)/2π).

The result is congruent with the data and exact wherever the reference is within ±π of the true phase, independently of all other pixels: no path is followed, so residues, noise and disconnected regions do not spread errors. Typical references are a TIE phase (pyoptics.TIE), the phase at a synthetic wavelength (unwrap_two_wavelength()), a model of the sample, or a smoothed least squares solution.

Parameters:

phase_wrappedarray

Wrapped phase ψ (any shape). NaNs are treated as masked pixels.

referencearray or number

Unwrapped estimate of the phase, with the same piston as the data (e.g. both zero on the background).

maskboolean array, optional

True for valid pixels.

Returns:

phasearray

Unwrapped phase, NaN outside the mask. Unlike the path following methods, the piston is not changed: it is the reference’s.

pyoptics.unwrap.synthetic_wavelength(wavelength1, wavelength2)

Synthetic (beat) wavelength λ1·λ2/|λ1 − λ2| of two wavelengths.

Parameters:

wavelength1, wavelength2numbers

The two wavelengths (different).

Returns:

wavelength_syntheticnumber

pyoptics.unwrap.unwrap_two_wavelength(phase1, phase2, wavelength1, wavelength2, mask=None, refine=True)

Optical path difference from wrapped phases at two wavelengths.

Both phases belong to the same optical path difference (OPD), ψ_j = wrap(2π·OPD/λ_j). Their difference is the phase at the synthetic wavelength Λ = λ1·λ2/|λ1 − λ2| (synthetic_wavelength()),

OPD_Λ = Λ_s·wrap(ψ1 − ψ2)/2π,   Λ_s = λ1·λ2/(λ2 − λ1),

which is unambiguous for |OPD| < Λ/2, without any path following. Its noise is that of the single phases times Λ/λ, so with refine=True the coarse OPD is only used as reference to unwrap each phase (unwrap_with_reference()), and the two results are averaged. That needs the coarse OPD to be within ±λ/2 of the truth (hierarchical unwrapping, Gass et al. 2003).

Parameters:

phase1, phase2arrays

Wrapped phases at wavelength1 and wavelength2, in radians, with a common zero (e.g. both background compensated).

wavelength1, wavelength2numbers

The (vacuum) wavelengths, which must differ.

maskboolean array, optional

True for valid pixels.

refinebool, optional

If False, return the coarse synthetic wavelength OPD.

Returns:

opdarray

Optical path difference in the unit of the wavelengths, NaN outside the mask. The phase at λ_j is 2π·opd/λ_j. Dispersion of the sample (a different OPD at the two wavelengths) shows up as an error of Λ/λ times the OPD difference in the coarse result.

pyoptics.unwrap.unwrap(phase_wrapped, method='quality', **kwargs)

Unwrap a 1D or 2D phase.

Parameters:

phase_wrappedarray

Wrapped phase.

methodstr, optional

“quality” (unwrap_quality_guided(), default), “ls” (unwrap_ls()), “reference” (unwrap_with_reference(), needs reference=…) or “1d” (unwrap_1d(); the default for 1D input without a reference).

kwargs

Passed on to the selected function.

Returns:

phasearray

Unwrapped phase.