How to use the pylops.signalprocessing function in pylops

# Blurring guassian operator
nh = [15, 25]
hz = np.exp(-0.1*np.linspace(-(nh[0]//2), nh[0]//2, nh[0])**2)
hx = np.exp(-0.03*np.linspace(-(nh[1]//2), nh[1]//2, nh[1])**2)
hz /= np.trapz(hz) # normalize the integral to 1
hx /= np.trapz(hx) # normalize the integral to 1
h = hz[:, np.newaxis] * hx[np.newaxis, :]

fig, ax = plt.subplots(1, 1, figsize=(5, 3))
him = ax.imshow(h)
ax.set_title('Blurring operator')
fig.colorbar(him, ax=ax)
ax.axis('tight')

Cop = pylops.signalprocessing.Convolve2D(Nz * Nx, h=h,
                                         offset=(nh[0] // 2,
                                                 nh[1] // 2),
                                         dims=(Nz, Nx), dtype='float32')

###############################################################################
# We first apply the blurring operator to the sharp image. We then
# try to recover the sharp input image by inverting the convolution operator
# from the blurred image. Note that when we perform inversion without any
# regularization, the deblurred image will show some ringing due to the
# instabilities of the inverse process. Using a L1 solver with a DWT
# preconditioner or TV regularization allows to recover sharper contrasts.
imblur = Cop * im.flatten()

imdeblur = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Cop, None,
                                                              imblur,

levels_size = np.flip(np.array([2 ** i for i in range(nlevels_max)]))
levels_cum = np.cumsum(levels_size)

plt.figure(figsize=(14, 5))
plt.imshow(seis.T, cmap='gray', vmin=-clip, vmax=clip)
for level in levels_cum:
    plt.axvline(level-0.5, color='w')
plt.title('Seislet transform')
plt.colorbar()
plt.axis('tight')

############################################
# 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.T),
                                      sampling=(dx, dt))
dwt = Wop * d.ravel()
dwt = dwt.reshape(nx, nt)

plt.figure(figsize=(14, 5))
plt.imshow(dwt.T, cmap='gray', vmin=-clip, vmax=clip)
for level in levels_cum:
    plt.axvline(level-0.5, color='w')
plt.title('Wavelet transform')
plt.colorbar()
plt.axis('tight')

############################################
# Finally we evaluate the compression capabilities of the Seislet transform.
# We zero-out all coefficients at the first two fine resolutions and keep those
# at coarser resolutions. We perform the inverse Seislet transform and asses

plt.title('TV inversion')

###############################################################################
# Finally, we repeat the same exercise on a 2-dimensional image. In this case
# we mock a medical imaging problem: the data is created by appling a 2D
# Fourier Transform to the input model and by randomly sampling 60% of
# its values.
x = np.load('../testdata/optimization/shepp_logan_phantom.npy')
x = x/x.max()
ny, nx = x.shape

perc_subsampling = 0.6
nxsub = int(np.round(ny*nx*perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(ny*nx))[:nxsub])
Rop = pylops.Restriction(ny*nx, iava, dtype=np.complex)
Fop = pylops.signalprocessing.FFT2D(dims=(ny, nx))

n = np.random.normal(0, 0., (ny, nx))
y = Rop*Fop*(x.flatten() + n.flatten())
yfft = Fop*(x.flatten() + n.flatten())
yfft = np.fft.fftshift(yfft.reshape(ny, nx))

ymask = Rop.mask(Fop*(x.flatten()) + n.flatten())
ymask = ymask.reshape(ny, nx)
ymask.data[:] = np.fft.fftshift(ymask.data)
ymask.mask[:] = np.fft.fftshift(ymask.mask)

fig, axs = plt.subplots(1, 3, figsize=(14, 5))
axs[0].imshow(x, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(np.abs(yfft), vmin=0, vmax=1, cmap='rainbow')

###############################################################################
# We repeat a similar exercise but using two dimensional signals and
# filters taking advantage of the
# :py:class:`pylops.signalprocessing.Convolve2D` operator.
nt = 51
nx = 81
dt = 0.004
t = np.arange(nt)*dt
x = np.zeros((nt, nx))
x[int(nt/2), int(nx/2)] = 1

nh = [11, 5]
h = np.ones((nh[0], nh[1]))

Cop = pylops.signalprocessing.Convolve2D(nt * nx, h=h,
                                         offset=(int(nh[0])/2, int(nh[1])/2),
                                         dims=(nt, nx), dtype='float32')
y = Cop*x.flatten()
xinv = Cop / y

y = y.reshape(nt, nx)
xinv = xinv.reshape(nt, nx)

fig, axs = plt.subplots(1, 3, figsize=(10, 3))
fig.suptitle('Convolve 2d data', fontsize=14,
             fontweight='bold', y=0.95)
axs[0].imshow(x, cmap='gray', vmin=-1, vmax=1)
axs[1].imshow(y, cmap='gray', vmin=-1, vmax=1)
axs[2].imshow(xinv, cmap='gray', vmin=-1, vmax=1)
axs[0].set_title('x')
axs[0].axis('tight')

pylops.signalprocessing.Radon2D(t, np.linspace(-par['dx']*winsize//2,
                                                   par['dx']*winsize//2,
                                                   winsize),
                                    px, centeredh=True, kind='linear',
                                    engine='numba')
Slid = pylops.signalprocessing.Sliding2D(Op, dims, dimsd,
                                         winsize, overlap,
                                         tapertype=None)

radon = Slid.H * data.flatten()
radon = radon.reshape(dims)

###############################################################################
# We now create  a similar operator but we also add a taper to the overlapping
# parts of the patches.
Slid = pylops.signalprocessing.Sliding2D(Op, dims, dimsd,
                                         winsize, overlap,
                                         tapertype='cosine')

reconstructed_data = Slid * radon.flatten()
reconstructed_data = reconstructed_data.reshape(dimsd)

###############################################################################
# We will see that our reconstructed signal presents some small artifacts.
# This is because we have not inverted our operator but simply applied
# the adjoint to estimate the representation of the input data in the Radon
# domain. We can do better if we use the inverse instead.
radoninv = pylops.LinearOperator(Slid, explicit=False).div(data.flatten(),
                                                           niter=10)
reconstructed_datainv = Slid * radoninv.flatten()

radoninv = radoninv.reshape(dims)

hy = np.arange(ny) * dy
hx = np.arange(nx) * dx

py = np.linspace(0, pymax, npy)
px = np.linspace(0, pxmax, npx)

x = np.zeros((npy, npx, nt))
x[npy//2, npx//2-2, nt//2] = 1

RLop = pylops.signalprocessing.Radon3D(t, hy, hx, py, px, centeredh=True,
                                       kind='linear', interp=False,
                                       engine='numpy')
RPop = pylops.signalprocessing.Radon3D(t, hy, hx, py, px, centeredh=True,
                                       kind='parabolic', interp=False,
                                       engine='numpy')
RHop = pylops.signalprocessing.Radon3D(t, hy, hx, py, px, centeredh=True,
                                       kind='hyperbolic', interp=False,
                                       engine='numpy')

# forward
yL = RLop * x.flatten()
yP = RPop * x.flatten()
yH = RHop * x.flatten()
yL = yL.reshape(ny, nx, nt)
yP = yP.reshape(ny, nx, nt)
yH = yH.reshape(ny, nx, nt)

# adjoint
xadjL = RLop.H * yL.flatten()
xadjP = RPop.H * yP.flatten()
xadjH = RHop.H * yH.flatten()
xadjL = xadjL.reshape(npy, npx, nt)