rows, cols = 6, 4
dx, dy = make_differentiation_matrices(rows, cols, boundary_conditions="periodic")
fig, axes = plt.subplots(1, 2)
axr = axes.ravel()
axr[0].imshow((dx.todense()), cmap="coolwarm")
axr[0].grid(True)
axr[0].set_title("Dx")
axr[1].imshow((dy.todense()), cmap="coolwarm")
axr[1].grid(True)
axr[1].set_title("Dy")
These matrices are used to differentiate the phase and obtain the gradient vector field $\vec{\phi} = [\phi_x, \phi_y]^T$.
The only adjustment fo this problem is to rewrap any derivatives outside the range $[-\pi, \pi]$:
Now the shrinkage operator:
The 0th axis is taken to be (x, y), and the magnitude is shrunk together (hence the sum along axis=0).
When $p = 0$, this is equal to
$$ S_0(x) = \text{max}(0, |x| - 1/|x|)\text{sign}(x) $$p = 0
xmin, xmax = -1.5, 1.5
xx = np.linspace(xmin, xmax).reshape((1, -1))
fig, ax = plt.subplots()
ax.plot(xx.ravel(), p_shrink(xx, p=p, epsilon=0).ravel())
ax.set_title(f"p = {p} shrinkage operator")
ax.set_ylim((xmin, xmax))
For the ADMM step of solving the linear system
$$||D \Phi - \phi||_2^2 = 0$$
the $D$ matrices can be diagonalized using the discrete cosine transform (DCT). Thus, instead of a linear solver, we compute the DCT of the right hand side $$ D_x^T \phi_x + D_x^T \phi_y $$ divide by the eigenvalues of the operator, then take the IDCT.
The make_laplace_kernel gives the inverse of the eigenvalues, so this can be multiplied by DCT(rhs).
def plot_diff(a, b, zero_mean=True, cmap="plasma", titles=None, vm=None, figsize=None):
if zero_mean:
print("subtracting mean")
aa = a - a.mean()
bb = b - b.mean()
else:
aa = a
bb = b
if titles is None:
titles = ["a", "b", "diff: a - b"]
if vm is None:
vmax = max(np.max(aa), np.max(bb))
vmin = min(np.min(aa), np.min(bb))
else:
vmin, vmax = -vm, vm
fig, axes = plt.subplots(1, 3, sharex=True, sharey=True, figsize=figsize)
axim = axes[0].imshow(aa, cmap=cmap, vmin=vmin, vmax=vmax)
fig.colorbar(axim, ax=axes[0])
axes[0].set_title(titles[0])
axim2 = axes[1].imshow(bb, cmap=cmap, vmin=vmin, vmax=vmax)
fig.colorbar(axim2, ax=axes[1])
axes[1].set_title(titles[1])
axim3 = axes[2].imshow(aa - bb, cmap=cmap)
fig.colorbar(axim3, ax=axes[2])
if len(titles) >= 3:
t = titles[2]
else:
t = f"Diff: {titles[0]} - {titles[1]}"
axes[2].set_title(t)
def make_ramp_test(r=20, c=10, sigma=0.0):
rramp = 2 * np.arange(-r // 2, r // 2) / 100.0 * 10
cramp = 2 * np.arange(-c // 2, c // 2) / 100.0 * 10
# steady ramp from top to bot in both X, Y directions
X, Y = np.meshgrid(cramp - cramp[0], rramp - rramp[0])
f0 = np.abs(X) + np.abs(Y)
f0 = np.max(f0) - f0
# noise level
# Unwrapped interferogram phase is f
f = f0 + sigma * np.random.randn(r, c)
# Wrapped version of unwrapped truth phase:
f_wrapped = ((f + np.pi) % (2 * np.pi)) - np.pi
return f, f_wrapped
f, f_wrapped = make_ramp_test()
plot_diff(f, f_wrapped)
Dx, Dy = make_differentiation_matrices(*f_wrapped.shape, boundary_conditions="neumann")
phi_x, phi_y = est_wrapped_gradient(f_wrapped, Dx=Dx, Dy=Dy)
fig, axes = plt.subplots(1, 2)
axim = axes[0].imshow(phi_x, cmap="plasma")
fig.colorbar(axim, ax=axes[0])
axim = axes[1].imshow(phi_y, cmap="plasma")
fig.colorbar(axim, ax=axes[1])
assert np.allclose(phi_x[:, :-1], -0.2)
assert np.allclose(phi_x[:, -1], 0)
assert np.allclose(phi_y[:-1, :], -0.2)
assert np.allclose(phi_y[-1, :], 0)
%%time
F = unwrap(f_wrapped, c=1.6, p=0, debug=True)
Note that the unwrapping should produce the same as the truth unwrapped, up to a constant. Therefore, we subtract the means to compare them:
plot_diff(F, f, zero_mean=True)
assert np.allclose(F - F.mean(), f - f.mean(), atol=1e-6)
igram = load_interferogram("../test/data/20150328_20150409.int")
igram_phase = np.angle(igram)
print(f"igrams phase dtype, shape: {igram_phase.dtype}, {igram.shape}")
rows, cols = igram.shape
%%time
unw_result = unwrap(
igram_phase, max_iters=100, c=1.3, p=0, lmbda=1, tol=np.pi / 10, debug=True
)
magphase_data = np.fromfile("../test/data/20150328_20150409_snaphu.unw", dtype="float32")
unw_snaphu = magphase_data.reshape((rows, 2 * cols))[:, cols:]
plt.figure()
plt.imshow(igram_phase, cmap="plasma")
plt.colorbar()
plt.title("Wrapped Phase")
plot_diff(
unw_result,
unw_snaphu,
zero_mean=True,
titles=["Algo. Result", "SNAPHU"],
figsize=(15, 4),
)
differences = (unw_result - unw_snaphu).reshape(-1)
differences -= differences.mean()
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].hist(differences, bins=51)
axes[0].set_xlabel("ours - snaphu")
axes[0].set_ylabel("pixel count")
axes[1].hist(differences, bins=51)
axes[1].set_xlabel("ours - snaphu")
axes[1].set_ylabel("pixel count (log scale)")
axes[1].set_yscale("log", nonpositive="clip")
