Note
Go to the end to download the full example code
Seislet transform#
This example shows how to use the pylops.signalprocessing.Seislet
operator. This operator the forward, adjoint and inverse Seislet transform
that is a modification of the well-know Wavelet transform where local slopes
are used in the prediction and update steps to further improve the prediction
of a trace from its previous (or subsequent) one and reduce the amount of
information passed to the subsequent scale. While this transform was initially
developed in the context of processing and compression of seismic data, it is
also suitable to any other oscillatory dataset such as GPR or Acoustic
recordings.
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import MaxNLocator
from mpl_toolkits.axes_grid1 import make_axes_locatable
import pylops
plt.close("all")
In this example we use the same benchmark
dataset
that was used in the original paper describing the Seislet transform. First,
local slopes are estimated using
pylops.utils.signalprocessing.slope_estimate
.
inputfile = "../testdata/sigmoid.npz"
d = np.load(inputfile)
d = d["sigmoid"]
nx, nt = d.shape
dx, dt = 0.008, 0.004
x, t = np.arange(nx) * dx, np.arange(nt) * dt
# slope estimation
slope, _ = pylops.utils.signalprocessing.slope_estimate(d.T, dt, dx, smooth=2.5)
slope = -slope.T # t-axis points down, reshape
# clip slopes above 80°
pmax = np.arctan(80 * np.pi / 180)
slope[slope > pmax] = pmax
slope[slope < -pmax] = -pmax
clip = 0.5 * np.max(np.abs(d))
clip_s = min(pmax, np.max(np.abs(slope)))
opts = dict(aspect=2, extent=(x[0], x[-1], t[-1], t[0]))
fig, axs = plt.subplots(1, 2, figsize=(14, 7), sharey=True, sharex=True)
axs[0].imshow(d.T, cmap="gray", vmin=-clip, vmax=clip, **opts)
axs[0].set(xlabel="Position [km]", ylabel="Time [s]", title="Data")
im = axs[1].imshow(slope.T, cmap="RdBu_r", vmin=-clip_s, vmax=clip_s, **opts)
axs[1].set(xlabel="Position [km]", title="Slopes")
fig.tight_layout()
pos = axs[1].get_position()
cbpos = [
pos.x0 + 0.1 * pos.width,
pos.y0 + 0.9 * pos.height,
0.8 * pos.width,
0.05 * pos.height,
]
cax = fig.add_axes(cbpos)
cb = fig.colorbar(im, cax=cax, orientation="horizontal")
cb.set_label("[s/km]")

Next the Seislet transform is computed.
Sop = pylops.signalprocessing.Seislet(slope, sampling=(dx, dt))
seis = Sop * d
nlevels_max = int(np.log2(nx))
levels_size = np.flip(np.array([2**i for i in range(nlevels_max)]))
levels_cum = np.cumsum(levels_size)
fig, ax = plt.subplots(figsize=(14, 6))
im = ax.imshow(
seis.T,
cmap="gray",
vmin=-clip,
vmax=clip,
aspect="auto",
interpolation="none",
extent=(1, seis.shape[0], t[-1], t[0]),
)
ax.xaxis.set_major_locator(MaxNLocator(nbins=20, integer=True))
for level in levels_cum:
ax.axvline(level + 0.5, color="w")
ax.set(xlabel="Scale", ylabel="Time [s]", title="Seislet transform")
cax = make_axes_locatable(ax).append_axes("right", size="2%", pad=0.1)
cb = fig.colorbar(im, cax=cax, orientation="vertical")
cb.formatter.set_powerlimits((0, 0))
fig.tight_layout()

We may also stretch the finer scales to be the width of the image
fig, axs = plt.subplots(2, nlevels_max // 2, figsize=(14, 7), sharex=True, sharey=True)
for i, ax in enumerate(axs.ravel()[:-1]):
curdata = seis[levels_cum[i] : levels_cum[i + 1], :].T
vmax = np.max(np.abs(curdata))
ax.imshow(curdata, vmin=-vmax, vmax=vmax, cmap="gray", interpolation="none", **opts)
ax.set(title=f"Scale {i+1}")
if i + 1 > nlevels_max // 2:
ax.set(xlabel="Position [km]")
curdata = seis[levels_cum[-1] :, :].T
vmax = np.max(np.abs(curdata))
axs[-1, -1].imshow(
curdata, vmin=-vmax, vmax=vmax, cmap="gray", interpolation="none", **opts
)
axs[0, 0].set(ylabel="Time [s]")
axs[1, 0].set(ylabel="Time [s]")
axs[-1, -1].set(xlabel="Position [km]", title=f"Scale {nlevels_max}")
fig.tight_layout()

As a comparison we also compute the Seislet transform fixing slopes to zero. This way we turn the Seislet tranform into a basic 1D Wavelet transform performed over the spatial axis.
Wop = pylops.signalprocessing.Seislet(np.zeros_like(slope), sampling=(dx, dt))
dwt = Wop * d
fig, ax = plt.subplots(figsize=(14, 6))
im = ax.imshow(
dwt.T,
cmap="gray",
vmin=-clip,
vmax=clip,
aspect="auto",
interpolation="none",
extent=(1, dwt.shape[0], t[-1], t[0]),
)
ax.xaxis.set_major_locator(MaxNLocator(nbins=20, integer=True))
for level in levels_cum:
ax.axvline(level + 0.5, color="w")
ax.set(xlabel="Scale", ylabel="Time [s]", title="Wavelet transform")
cax = make_axes_locatable(ax).append_axes("right", size="2%", pad=0.1)
cb = fig.colorbar(im, cax=cax, orientation="vertical")
cb.formatter.set_powerlimits((0, 0))
fig.tight_layout()

Again, we may decompress the finer scales
fig, axs = plt.subplots(2, nlevels_max // 2, figsize=(14, 7), sharex=True, sharey=True)
for i, ax in enumerate(axs.ravel()[:-1]):
curdata = dwt[levels_cum[i] : levels_cum[i + 1], :].T
vmax = np.max(np.abs(curdata))
ax.imshow(curdata, vmin=-vmax, vmax=vmax, cmap="gray", interpolation="none", **opts)
ax.set(title=f"Scale {i+1}")
if i + 1 > nlevels_max // 2:
ax.set(xlabel="Position [km]")
curdata = dwt[levels_cum[-1] :, :].T
vmax = np.max(np.abs(curdata))
axs[-1, -1].imshow(
curdata, vmin=-vmax, vmax=vmax, cmap="gray", interpolation="none", **opts
)
axs[0, 0].set(ylabel="Time [s]")
axs[1, 0].set(ylabel="Time [s]")
axs[-1, -1].set(xlabel="Position [km]", title=f"Scale {nlevels_max}")
fig.tight_layout()

Finally we evaluate the compression capabilities of the Seislet transform compared to the 1D Wavelet transform. We zero-out all but the strongest 25% of the components. We perform the inverse transforms and assess the compression error.
perc = 0.25
seis_strong_idx = np.argsort(-np.abs(seis.ravel()))
dwt_strong_idx = np.argsort(-np.abs(dwt.ravel()))
seis_strong = np.abs(seis.ravel())[seis_strong_idx]
dwt_strong = np.abs(dwt.ravel())[dwt_strong_idx]
fig, ax = plt.subplots()
ax.plot(range(1, len(seis_strong) + 1), seis_strong / seis_strong[0], label="Seislet")
ax.plot(
range(1, len(dwt_strong) + 1), dwt_strong / dwt_strong[0], "--", label="Wavelet"
)
ax.set(xlabel="n", ylabel="Coefficient strength [%]", title="Transform Coefficients")
ax.axvline(np.rint(len(seis_strong) * perc), color="k", label=f"{100*perc:.0f}%")
ax.legend()
fig.tight_layout()

seis1 = np.zeros_like(seis.ravel())
seis_strong_idx = seis_strong_idx[: int(np.rint(len(seis_strong) * perc))]
seis1[seis_strong_idx] = seis.ravel()[seis_strong_idx]
d_seis = Sop.inverse(seis1).reshape(Sop.dims)
dwt1 = np.zeros_like(dwt.ravel())
dwt_strong_idx = dwt_strong_idx[: int(np.rint(len(dwt_strong) * perc))]
dwt1[dwt_strong_idx] = dwt.ravel()[dwt_strong_idx]
d_dwt = Wop.inverse(dwt1).reshape(Wop.dims)
opts.update(dict(cmap="gray", vmin=-clip, vmax=clip))
fig, axs = plt.subplots(2, 3, figsize=(14, 7), sharex=True, sharey=True)
axs[0, 0].imshow(d.T, **opts)
axs[0, 0].set(title="Data")
axs[0, 1].imshow(d_seis.T, **opts)
axs[0, 1].set(title=f"Rec. from Seislet ({100*perc:.0f}% of coeffs.)")
axs[0, 2].imshow((d - d_seis).T, **opts)
axs[0, 2].set(title="Error from Seislet Rec.")
axs[1, 0].imshow(d.T, **opts)
axs[1, 0].set(ylabel="Time [s]", title="Data [Repeat]")
axs[1, 1].imshow(d_dwt.T, **opts)
axs[1, 1].set(title=f"Rec. from Wavelet ({100*perc:.0f}% of coeffs.)")
axs[1, 2].imshow((d - d_dwt).T, **opts)
axs[1, 2].set(title="Error from Wavelet Rec.")
for i in range(3):
axs[1, i].set(xlabel="Position [km]")
plt.tight_layout()
![Data, Rec. from Seislet (25% of coeffs.), Error from Seislet Rec., Data [Repeat], Rec. from Wavelet (25% of coeffs.), Error from Wavelet Rec.](../_images/sphx_glr_plot_seislet_007.png)
To conclude it is worth noting that the Seislet transform, differently to the Wavelet transform, is not orthogonal: in other words, its adjoint and inverse are not equivalent. While we have used the forward and inverse transformations, when used as linear operator in composition with other operators, the Seislet transform requires the adjoint be defined and that it also passes the dot-test pair that is. As shown below, this is the case when using the implementation in the PyLops package.
pylops.utils.dottest(Sop, verb=True)
Dot test passed, v^H(Opu)=-98.52334183020218 - u^H(Op^Hv)=-98.52334183020227
True
Total running time of the script: ( 0 minutes 13.599 seconds)