Feat: pad_to next methods

docs/source/api/others.rst

@@ -102,7 +102,11 @@ Sliding and Patching
     sliding3d.sliding3d_design
     patch2d.patch2d_design
     patch3d.patch3d_design
-
+    sliding1d.sliding1d_pad_to_next
+    sliding2d.sliding2d_pad_to_next
+    sliding3d.sliding3d_pad_to_next
+    patch2d.patch2d_pad_to_next
+    patch3d.patch3d_pad_to_next
 
 
 Synthetics

examples/plot_patching.py

@@ -11,6 +11,7 @@
 that are 2- or 3-dimensional in nature, respectively.
 
 """
+
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -107,6 +108,47 @@
 
 ###############################################################################
 # We repeat now the same exercise in 3d
+
+
+def plot_3d(data, par, x, y, t, title):
+    fig, axs = plt.subplots(1, 3, figsize=(12, 5))
+    fig.suptitle(title, fontsize=12, fontweight="bold", y=0.95)
+    axs[0].imshow(
+        data[par["ny"] // 2].T,
+        aspect="auto",
+        interpolation="nearest",
+        vmin=-2,
+        vmax=2,
+        cmap="gray",
+        extent=(x.min(), x.max(), t.max(), t.min()),
+    )
+    axs[0].set_xlabel(r"$x(m)$")
+    axs[0].set_ylabel(r"$t(s)$")
+    axs[1].imshow(
+        data[:, par["nx"] // 2].T,
+        aspect="auto",
+        interpolation="nearest",
+        vmin=-2,
+        vmax=2,
+        cmap="gray",
+        extent=(y.min(), y.max(), t.max(), t.min()),
+    )
+    axs[1].set_xlabel(r"$y(m)$")
+    axs[1].set_ylabel(r"$t(s)$")
+    axs[2].imshow(
+        data[:, :, par["nt"] // 2],
+        aspect="auto",
+        interpolation="nearest",
+        vmin=-2,
+        vmax=2,
+        cmap="gray",
+        extent=(x.min(), x.max(), y.max(), x.min()),
+    )
+    axs[2].set_xlabel(r"$x(m)$")
+    axs[2].set_ylabel(r"$y(m)$")
+    plt.tight_layout()
+
+
 par = {
     "oy": -60,
     "dy": 2,
@@ -134,43 +176,7 @@
 # Generate model
 _, data = pylops.utils.seismicevents.hyperbolic3d(x, y, t, t0, vrms, vrms, amp, wav)
 
-
-fig, axs = plt.subplots(1, 3, figsize=(12, 5))
-fig.suptitle("Original data", fontsize=12, fontweight="bold", y=0.95)
-axs[0].imshow(
-    data[par["ny"] // 2].T,
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(x.min(), x.max(), t.max(), t.min()),
-)
-axs[0].set_xlabel(r"$x(m)$")
-axs[0].set_ylabel(r"$t(s)$")
-axs[1].imshow(
-    data[:, par["nx"] // 2].T,
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(y.min(), y.max(), t.max(), t.min()),
-)
-axs[1].set_xlabel(r"$y(m)$")
-axs[1].set_ylabel(r"$t(s)$")
-axs[2].imshow(
-    data[:, :, par["nt"] // 2],
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(x.min(), x.max(), y.max(), x.min()),
-)
-axs[2].set_xlabel(r"$x(m)$")
-axs[2].set_ylabel(r"$y(m)$")
-plt.tight_layout()
+plot_3d(data, par, x, y, t, "Original data")
 
 ###############################################################################
 # Let's create now the :py:class:`pylops.signalprocessing.Patch3D` operator
@@ -201,39 +207,76 @@
 )
 reconstructed_data = np.real(Patch * fftdata)
 
-fig, axs = plt.subplots(1, 3, figsize=(12, 5))
-fig.suptitle("Reconstructed data", fontsize=12, fontweight="bold", y=0.95)
-axs[0].imshow(
-    reconstructed_data[par["ny"] // 2].T,
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(x.min(), x.max(), t.max(), t.min()),
+plot_3d(reconstructed_data, par, x, y, t, "Reconstructed data")
+
+###############################################################################
+# However, in practice the shapes of data to be patched do not always conform
+# with the patch design, meaning that some of the edges do not fit within the
+# patches. In this case, we will use the
+# :func:`pylops.signalprocessing.patch3d_pad_to_next` method to pad the input
+# before feeding it to the patching operator.
+par = {
+    "oy": -66,
+    "dy": 2,
+    "ny": 66,
+    "ox": -54,
+    "dx": 2,
+    "nx": 54,
+    "ot": 0,
+    "dt": 0.004,
+    "nt": 100,
+    "f0": 20,
+}
+
+v = 1500
+t0 = [0.05, 0.2, 0.3]
+vrms = [500, 700, 1700]
+amp = [1.0, -2, 0.5]
+
+# Create axis
+t, t2, x, y = pylops.utils.seismicevents.makeaxis(par)
+
+# Create wavelet
+wav = pylops.utils.wavelets.ricker(t[:41], f0=par["f0"])[0]
+
+# Generate model
+_, data = pylops.utils.seismicevents.hyperbolic3d(x, y, t, t0, vrms, vrms, amp, wav)
+
+plot_3d(reconstructed_data, par, x, y, t, "Data")
+
+# Patching operator
+dimsd = data.shape
+nwins, dims, mwin_inends, dwin_inends = pylops.signalprocessing.patch3d_design(
+    dimsd, nwin, nover, (128, 128, 65)
 )
-axs[0].set_xlabel(r"$x(m)$")
-axs[0].set_ylabel(r"$t(s)$")
-axs[1].imshow(
-    reconstructed_data[:, par["nx"] // 2].T,
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(y.min(), y.max(), t.max(), t.min()),
+Patch = pylops.signalprocessing.Patch3D(
+    Op.H, dims, dimsd, nwin, nover, nop, tapertype=None
 )
-axs[1].set_xlabel(r"$y(m)$")
-axs[1].set_ylabel(r"$t(s)$")
-axs[2].imshow(
-    reconstructed_data[:, :, par["nt"] // 2],
-    aspect="auto",
-    interpolation="nearest",
-    vmin=-2,
-    vmax=2,
-    cmap="gray",
-    extent=(x.min(), x.max(), y.max(), x.min()),
+fftdata = Patch.H * data
+
+Patch = pylops.signalprocessing.Patch3D(
+    Op.H, dims, dimsd, nwin, nover, nop, tapertype="hanning"
 )
-axs[2].set_xlabel(r"$x(m)$")
-axs[2].set_ylabel(r"$y(m)$")
-plt.tight_layout()
+reconstructed_data = np.real(Patch * fftdata)
+
+plot_3d(reconstructed_data, par, x, y, t, "Reconstructed data without padding")
+
+# Pad to next
+data_pad, nwins, dims, mwin_inends, dwin_inends = (
+    pylops.signalprocessing.patch3d.patch3d_pad_to_next(
+        data, nwin, nover, (128, 128, 65)
+    )
+)
+
+dimsd_pad = data_pad.shape
+Patch = pylops.signalprocessing.Patch3D(
+    Op.H, dims, dimsd_pad, nwin, nover, nop, tapertype=None
+)
+fftdata = Patch.H * data_pad
+
+Patch = pylops.signalprocessing.Patch3D(
+    Op.H, dims, dimsd_pad, nwin, nover, nop, tapertype="hanning"
+)
+reconstructed_data = np.real(Patch * fftdata)[: dimsd[0], : dimsd[1], : dimsd[2]]
+
+plot_3d(reconstructed_data, par, x, y, t, "Reconstructed data with padding")

pylops/signalprocessing/patch2d.py

@@ -1,5 +1,6 @@
 __all__ = [
     "patch2d_design",
+    "patch2d_pad_to_next",
     "Patch2D",
 ]
 
@@ -45,7 +46,7 @@ def patch2d_design(
     Parameters
     ----------
     dimsd : :obj:`tuple`
-        Shape of 2-dimensional data.
+        Shape of the 2-dimensional data.
     nwin : :obj:`tuple`
         Number of samples of window.
     nover : :obj:`tuple`
@@ -61,7 +62,7 @@ def patch2d_design(
     nwins : :obj:`tuple`
         Number of windows.
     dims : :obj:`tuple`
-        Shape of 2-dimensional model.
+        Shape of the 2-dimensional model.
     mwins_inends : :obj:`tuple`
         Start and end indices for model patches (stored as tuple of tuples).
     dwins_inends : :obj:`tuple`
@@ -101,6 +102,67 @@ def patch2d_design(
     return nwins, dims, mwins_inends, dwins_inends
 
 
+def patch2d_pad_to_next(
+    inpt: NDArray,
+    nwin: tuple[int, int],
+    nover: tuple[int, int],
+    nop: tuple[int, int],
+) -> tuple[
+    NDArray,
+    tuple[int, int],
+    tuple[int, int],
+    tuple[tuple[NDArray, NDArray], tuple[NDArray, NDArray]],
+    tuple[tuple[NDArray, NDArray], tuple[NDArray, NDArray]],
+]:
+    """Pad input to next set of patches
+
+    Pad ``inpt`` to the next set of patches such that the padded input is completely
+    filled by overlapping patches.
+
+    Parameters
+    ----------
+     inpt : :obj:`numpy.ndarray`
+        3-dimensional input data.
+    nwin : :obj:`tuple`
+        Number of samples of window.
+    nover : :obj:`tuple`
+        Number of samples of overlapping part of window.
+    nop : :obj:`tuple`
+        Size of model in the transformed domain.
+
+    Returns
+    -------
+    inpt_pad : :obj:`numpy.ndarray`
+        3-dimensional input data after padding along the first and second dimensions
+    nwins : :obj:`tuple`
+        Number of windows.
+    dims : :obj:`tuple`
+        Shape of the 2-dimensional model.
+    mwins_inends : :obj:`tuple`
+        Start and end indices for model patches (stored as tuple of tuples).
+    dwins_inends : :obj:`tuple`
+        Start and end indices for data patches (stored as tuple of tuples).
+
+    """
+    # Identify current sliding design
+    dimsd = inpt.shape
+    nwins, dims, mwins_inends, dwins_inends = patch2d_design(dimsd, nwin, nover, nop)
+
+    # Pad to next slice
+    pad = [0, 0]
+    if dwins_inends[0][1][-1] != dimsd[0]:
+        pad[0] = dwins_inends[0][1][-1] - nover[0] + nwin[0] - dimsd[0]
+    if dwins_inends[1][1][-1] != dimsd[1]:
+        pad[1] = dwins_inends[1][1][-1] - nover[1] + nwin[1] - dimsd[1]
+    inpt_pad = np.pad(inpt, ((0, pad[0]), (0, pad[1])))
+    if pad[0] > 0 or pad[1] > 0:
+        dimsd_pad = inpt_pad.shape
+        nwins, dims, mwins_inends, dwins_inends = patch2d_design(
+            dimsd_pad, nwin, nover, nop
+        )
+    return inpt_pad, nwins, dims, mwins_inends, dwins_inends
+
+
 class Patch2D(LinearOperator):
     """2D Patch transform operator.
 
@@ -132,11 +194,11 @@ class Patch2D(LinearOperator):
     Op : :obj:`pylops.LinearOperator`
         Transform operator
     dims : :obj:`tuple`
-        Shape of 2-dimensional model. Note that ``dims[0]`` and ``dims[1]``
+        Shape of the 2-dimensional model. Note that ``dims[0]`` and ``dims[1]``
         should be multiple of the model size of the transform in their
         respective dimensions
     dimsd : :obj:`tuple`
-        Shape of 2-dimensional data
+        Shape of the 2-dimensional data
     nwin : :obj:`tuple`
         Number of samples of window
     nover : :obj:`tuple`