Cubic Spline Interpolator

.gitignore

@@ -50,3 +50,6 @@ ASV
 .asv/
 asv.conf.json
 benchmarks/
+
+# Virtual Environments
+.venv**

docs/source/api/index.rst

@@ -64,7 +64,7 @@ Basic operators
     Imag
     Conj
     ToCupy
-    
+
 
 Smoothing and derivatives
 ~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -100,6 +100,7 @@ Signal processing
     NonStationaryFilters1D
     NonStationaryFilters2D
     Interp
+    InterpCubicSpline
     Bilinear
     FFT
     FFT2D

examples/plot_restriction.py

@@ -47,20 +47,41 @@
 # We then create the restriction and interpolation operators and display
 # the original signal as well as the subsampled signal.
 
-Rop = pylops.Restriction(nt, iava, dtype="float64")
+Rop = pylops.Restriction(
+    nt,
+    iava,
+    dtype="float64",
+)
 NNop, iavann = pylops.signalprocessing.Interp(
-    nt, iava + 0.4, kind="nearest", dtype="float64"
+    nt,
+    iava + 0.4,
+    kind="nearest",
+    dtype="float64",
 )
 LIop, iavali = pylops.signalprocessing.Interp(
-    nt, iava + 0.4, kind="linear", dtype="float64"
+    nt,
+    iava + 0.4,
+    kind="linear",
+    dtype="float64",
 )
 SIop, iavasi = pylops.signalprocessing.Interp(
-    nt, iava + 0.4, kind="sinc", dtype="float64"
+    nt,
+    iava + 0.4,
+    kind="sinc",
+    dtype="float64",
+)
+CuSIop, iavacusi = pylops.signalprocessing.Interp(
+    nt,
+    iava + 0.4,
+    kind="cubic_spline",
+    dtype="float64",
 )
+
 y = Rop * x
 ynn = NNop * x
 yli = LIop * x
 ysi = SIop * x
+ycusi = CuSIop * x
 ymask = Rop.mask(x)
 
 # Visualize data
@@ -70,6 +91,7 @@
 plt.plot(iavann, ynn, ".r", ms=25, label="NN interp samples")
 plt.plot(iavali, yli, ".m", ms=20, label="Linear interp samples")
 plt.plot(iavasi, ysi, ".y", ms=15, label="Sinc interp samples")
+plt.plot(iavasi, ysi, "cyan", linestyle="none", marker="D", ms=4, label="Cubic Spline interp samples")
 plt.legend(loc="right")
 plt.title("Data restriction")
 
@@ -79,6 +101,7 @@
 subax.plot(iavann, ynn, ".r", ms=25)
 subax.plot(iavali, yli, ".m", ms=20)
 subax.plot(iavasi, ysi, ".y", ms=15)
+subax.plot(iavasi, ysi, "cyan", linestyle="none", marker="D", ms=4)
 subax.set_xlim([120, 127])
 subax.set_ylim([-0.5, 0.5])
 plt.tight_layout()

pylops/signalprocessing/__init__.py

@@ -16,6 +16,7 @@
     NonStationaryFilters1D          1D nonstationary filter estimation operator.
     NonStationaryFilters2D          2D nonstationary filter estimation operator.
     Interp                          Interpolation operator.
+    InterpCubicSpline               Cubic Spline Interpolation operator.
     Bilinear                        Bilinear interpolation operator.
     FFT                             One dimensional Fast-Fourier Transform.
     FFT2D                           Two dimensional Fast-Fourier Transform.
@@ -55,6 +56,7 @@
 from .nonstatconvolve3d import *
 from .shift import *
 from .interp import *
+from .interpspline import *
 from .bilinear import *
 from .radon2d import *
 from .radon3d import *
@@ -90,6 +92,7 @@
     "NonStationaryFilters1D",
     "NonStationaryFilters2D",
     "Interp",
+    "InterpCubicSpline",
     "Bilinear",
     "Radon2D",
     "Radon3D",

pylops/signalprocessing/_interp_utils.py

@@ -0,0 +1,46 @@
+import warnings
+
+import numpy as np
+
+from pylops.utils.typing import NumericNDArray
+
+
+def ensure_iava_is_unique(
+    iava: NumericNDArray,
+    axis: int | None = None,
+) -> None:
+    _, count = np.unique(
+        iava,
+        axis=axis,
+        return_counts=True,
+    )
+
+    if np.any(count > 1):
+        raise ValueError("Found repeated/non-unique values in iava.")
+
+    return
+
+
+def clip_iava_above_last_sample_index(
+    iava: NumericNDArray,
+    sample_size: int,
+) -> None:
+    # ensure that samples are not beyond the last sample, in that case set to
+    # penultimate sample and raise a warning
+    last_sample_index = sample_size - 1
+    outside = iava >= last_sample_index
+    if np.any(outside):
+        warnings.warn(
+            f"At least one value in iava is beyond the penultimate sample index "
+            f"{last_sample_index}. Out-of-bound-values are forced below penultimate "
+            f"sample."
+        )
+
+        # NOTE: ``numpy.nextafter(x, -np.inf)`` gives the closest float-value that is
+        #       less than ``x``, i.e., this logic clips ``iava`` to the highest possible
+        #       value that is still below the last sample
+        iava[np.where(outside)] = np.nextafter(last_sample_index, -np.inf)
+
+    ensure_iava_is_unique(iava=iava)
+
+    return

pylops/signalprocessing/bilinear.py

@@ -3,23 +3,16 @@
 import logging
 
 import numpy as np
-import numpy.typing as npt
 
 from pylops import LinearOperator
+from pylops.signalprocessing._interp_utils import ensure_iava_is_unique
 from pylops.utils.backend import get_add_at, get_array_module, to_numpy
 from pylops.utils.decorators import reshaped
 from pylops.utils.typing import DTypeLike, InputDimsLike, IntNDArray, NDArray
 
 logger = logging.getLogger(__name__)
 
 
-def _checkunique(iava: npt.ArrayLike) -> None:
-    """Check that vector as only unique values"""
-    _, count = np.unique(iava, axis=1, return_counts=True)
-    if np.any(count > 1):
-        raise ValueError("Repeated values in iava array")
-
-
 class Bilinear(LinearOperator):
     r"""Bilinear interpolation operator.
 
@@ -140,7 +133,10 @@ def __init__(
 
         ncp = get_array_module(iava)
         # check non-unique pairs (works only with numpy arrays)
-        _checkunique(to_numpy(iava))
+        ensure_iava_is_unique(
+            iava=to_numpy(iava),
+            axis=1,
+        )
 
         # find indices and weights
         self.iava_t = ncp.floor(iava[0]).astype(int)

pylops/signalprocessing/interp.py

@@ -1,24 +1,21 @@
 __all__ = ["Interp"]
 
-import warnings
-from typing import Tuple, Union
+from typing import Literal, Tuple, Union
 
 import numpy as np
-import numpy.typing as npt
 
 from pylops import LinearOperator, aslinearoperator
 from pylops.basicoperators import Diagonal, MatrixMult, Restriction, Transpose
+from pylops.signalprocessing._interp_utils import (
+    clip_iava_above_last_sample_index,
+    ensure_iava_is_unique,
+)
+from pylops.signalprocessing.interpspline import InterpCubicSpline
 from pylops.utils._internal import _value_or_sized_to_tuple
 from pylops.utils.backend import get_array_module
 from pylops.utils.typing import DTypeLike, InputDimsLike, IntNDArray
 
 
-def _checkunique(iava: npt.ArrayLike) -> None:
-    _, count = np.unique(iava, return_counts=True)
-    if np.any(count > 1):
-        raise ValueError("Repeated values in iava array")
-
-
 def _nearestinterp(
     dims: Union[int, InputDimsLike],
     iava: IntNDArray,
@@ -27,8 +24,8 @@ def _nearestinterp(
 ):
     """Nearest neighbour interpolation."""
     iava = np.round(iava).astype(int)
-    _checkunique(iava)
-    return Restriction(dims, iava, axis=axis, dtype=dtype), iava
+    ensure_iava_is_unique(iava=iava)
+    return (Restriction(dims, iava, axis=axis, dtype=dtype), iava)
 
 
 def _linearinterp(
@@ -43,21 +40,14 @@ def _linearinterp(
     if np.issubdtype(iava.dtype, np.integer):
         iava = iava.astype(np.float64)
 
-    lastsample = dims[axis]
+    sample_size = dims[axis]
     dimsd = list(dims)
     dimsd[axis] = len(iava)
     dimsd = tuple(dimsd)
 
     # ensure that samples are not beyond the last sample, in that case set to
     # penultimate sample and raise a warning
-    outside = iava >= lastsample - 1
-    if sum(outside) > 0:
-        warnings.warn(
-            "At least one value is beyond the penultimate sample, "
-            "forced to be at penultimate sample"
-        )
-    iava[outside] = lastsample - 1 - 1e-10
-    _checkunique(iava)
+    clip_iava_above_last_sample_index(iava=iava, sample_size=sample_size)
 
     # find indices and weights
     iva_l = ncp.floor(iava).astype(int)
@@ -81,7 +71,9 @@ def _sincinterp(
 ):
     """Sinc interpolation."""
     ncp = get_array_module(iava)
-    _checkunique(iava)
+
+    # TODO: is ``iava`` bound to an integer dtype
+    ensure_iava_is_unique(iava=iava)
 
     # create sinc interpolation matrix
     nreg = dims[axis]
@@ -111,7 +103,7 @@ def Interp(
     dims: Union[int, InputDimsLike],
     iava: IntNDArray,
     axis: int = -1,
-    kind: str = "linear",
+    kind: Literal["linear", "nearest", "sinc", "cubic_spline"] = "linear",
     dtype: DTypeLike = "float64",
     name: str = "I",
 ) -> Tuple[LinearOperator, IntNDArray]:
@@ -136,11 +128,19 @@ def Interp(
       cost as it involves multiplying the input data by a matrix of size
       :math:`N \times M`.
 
+    - *Cubic Spline interpolation* relies on a cubic spline, i.e., a 2-times
+      continuously differentiable piecewise third order polynomial with equally spaced
+      knots. It is interpolated at the locations ``iava`` by evaluating the respective
+      polynomial fitted between ``np.floor(iava)`` and ``np.floor(iava) + 1``.
+      It offers an excellent tradeoff between accuracy and computational complexity
+      and its results oscillate less than those obtained from sinc interpolation.
+      It can also be accessed directly via :class:`pylops.singalprocessing.InterpCubicSpline`.
+
     .. note:: The vector ``iava`` should contain unique values. If the same
       index is repeated twice an error will be raised. This also applies when
       values beyond the last element of the input array for
-      *linear interpolation* as those values are forced to be just before this
-      element.
+      *linear interpolation* and *cubic spline interpolation* as those values are forced
+      to be just before this element.
 
     Parameters
     ----------
@@ -153,8 +153,13 @@ def Interp(
 
         Axis along which interpolation is applied.
     kind : :obj:`str`, optional
-        Kind of interpolation (``nearest``, ``linear``, and ``sinc`` are
-        currently supported)
+        Kind of interpolation.
+        Currently, ``"nearest"``, ``"linear"``, ``"sinc"``, and ``"cubic_spline"`` are
+        available.
+
+        .. versionadded:: 2.0.0
+
+        The ``"cubic_spline"``-interpolation was added.
     dtype : :obj:`str`, optional
         Type of elements in input array.
     name : :obj:`str`, optional
@@ -165,7 +170,7 @@ def Interp(
     Returns
     -------
     op : :obj:`pylops.LinearOperator`
-        Linear intepolation operator
+        Linear interpolation operator
     iava : :obj:`list` or :obj:`numpy.ndarray`
         Corrected indices of locations of available samples
         (samples at ``M-1`` or beyond are forced to be at ``M-1-eps``)
@@ -216,6 +221,11 @@ def Interp(
     :math:`i,j` possible combinations.
 
     """
+
+    kind = kind.lower()  # type: ignore
+    if kind not in {"nearest", "linear", "sinc", "cubic_spline"}:
+        raise NotImplementedError(f"{kind} interpolation could not be found.")
+
     dims = _value_or_sized_to_tuple(dims)
 
     if kind == "nearest":
@@ -224,11 +234,19 @@ def Interp(
         interpop, iava, dims, dimsd = _linearinterp(dims, iava, axis=axis, dtype=dtype)
     elif kind == "sinc":
         interpop, dims, dimsd = _sincinterp(dims, iava, axis=axis, dtype=dtype)
-    else:
-        raise NotImplementedError("kind is not correct...")
+    elif kind == "cubic_spline":
+        interpop = InterpCubicSpline(
+            dims=dims,
+            iava=iava,
+            axis=axis,
+            dtype=dtype,
+            name=name,
+        )
+        iava = interpop.iava
+
     # add dims and dimsd to composite operators (not needed for neareast as
     # interpop is a Restriction operator already
-    if kind != "nearest":
+    if kind not in {"nearest", "cubic_spline"}:
         interpop = aslinearoperator(interpop)
         interpop.dims = dims
         interpop.dimsd = dimsd