diff --git a/.github/workflows/flake8.yaml b/.github/workflows/flake8.yaml
deleted file mode 100644
index 6d9916d0..00000000
--- a/.github/workflows/flake8.yaml
+++ /dev/null
@@ -1,24 +0,0 @@
-# This workflow runs Flake8 on the PR
-# For more information see: https://github.com/marketplace/actions/python-flake8-lint
-name: PyLops-flake8
-
-on: [push, pull_request]
-
-jobs:
-  flake8-lint:
-    runs-on: ubuntu-latest
-    name: Lint
-    steps:
-      - name: Check out source repository
-        uses: actions/checkout@v4
-      - name: Set up Python environment
-        uses: actions/setup-python@v5
-        with:
-          python-version: "3.11"
-      - name: flake8 Lint
-        uses: py-actions/flake8@v2
-        with:
-          ignore: "E203,E501,W503,E402"
-          max-line-length: "88"
-          path: "pylops"
-          args: "--per-file-ignores=__init__.py:F401,F403,F405"
diff --git a/examples/plot_slopeest.py b/examples/plot_slopeest.py
index 3d9ba4d2..7fa9e226 100755
--- a/examples/plot_slopeest.py
+++ b/examples/plot_slopeest.py
@@ -2,8 +2,9 @@
 Slope estimation via Structure Tensor algorithm
 ===============================================
 
-This example shows how to estimate local slopes or local dips of a two-dimensional
-array using :py:func:`pylops.utils.signalprocessing.slope_estimate` and
+This example shows how to estimate local slopes or local dips of two-
+or three-dimensional arrays using
+:py:func:`pylops.utils.signalprocessing.slope_estimate` and
 :py:func:`pylops.utils.signalprocessing.dip_estimate`.
 
 Knowing the local slopes of an image (or a seismic data) can be useful for
@@ -34,6 +35,14 @@
 plt.close("all")
 np.random.seed(10)
 
+
+def create_colorbar(fig, ax, im):
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size="5%", pad=0.1)
+    cb = fig.colorbar(im, cax=cax, orientation="vertical")
+    return cax, cb
+
+
 ###############################################################################
 # Python logo
 # -----------
@@ -280,6 +289,136 @@ def rgb2gray(rgb):
     ax.axis("off")
 fig.tight_layout()
 
+
+###############################################################################
+# 3D Seismic data
+# ---------------
+# Finally, we consider once again a seismic dataset. However, in this case
+# we create a 3D version of it by linearly shifting traces over the
+# y-direction. We then apply the 3D version of the Structure Tensor algorithm
+# to such a data and display the slopes along the y- and x-directions.
+
+# Input
+inputfile = "../testdata/sigmoid.npz"
+
+sigmoid = np.load(inputfile)["sigmoid"]
+ny, nx, nt = 20, *sigmoid.shape
+dy, dx, dt = 0.008, 0.008, 0.004
+
+# Axes
+y, x, t = np.arange(ny) * dy, np.arange(nx) * dx, np.arange(nt) * dt
+
+# 3D Input
+sigmoid3d = np.zeros((ny, nx, nt), dtype=sigmoid.dtype)
+for i in range(ny):
+    sigmoid3d[i, ...] = np.roll(sigmoid, i // 2, axis=1)
+
+clip = 0.5 * np.max(np.abs(sigmoid3d))
+fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True, width_ratios=(3, 1))
+fig.suptitle("Sigmoid model")
+ax[0].imshow(
+    sigmoid3d[ny // 2].T,
+    aspect="auto",
+    vmin=-clip,
+    vmax=clip,
+    cmap="gray",
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+ax[0].set_ylabel("Time [s]")
+ax[0].set_xlabel("X [km]")
+im1 = ax[1].imshow(
+    sigmoid3d[:, nx // 2].T,
+    aspect="auto",
+    vmin=-clip,
+    vmax=clip,
+    cmap="gray",
+    extent=(y[0], y[-1], t[-1], t[0]),
+)
+ax[1].set_xlabel("Y [km]")
+fig.tight_layout()
+
+# Slope estimation
+slopes, anisotropies = slope_estimate(
+    sigmoid3d / sigmoid3d.max(),
+    dy=dy,
+    dx=dx,
+    dz=dt,
+    smooth=3,
+    eps=0.0,
+    dips=True,
+    anisotropies=True,
+    batch_size=500_000,
+)
+
+st_slope3d_y, st_slope3d_x = slopes
+st_linearity, st_planarity = anisotropies
+
+# Revert t-axis downward
+st_slope3d_y *= -1
+st_slope3d_x *= -1
+
+# Clip slopes above 80°
+pmax = np.arctan(80 * np.pi / 180)
+st_slope3d_y[st_slope3d_y > pmax] = pmax
+st_slope3d_y[st_slope3d_y < -pmax] = -pmax
+
+st_slope3d_x[st_slope3d_x > pmax] = pmax
+st_slope3d_x[st_slope3d_x < -pmax] = -pmax
+
+# Visualize
+v = 1  # np.max(np.abs(dips))
+clip_sy = min(pmax, np.max(np.abs(st_slope3d_y)))
+clip_sx = min(pmax, np.max(np.abs(st_slope3d_x)))
+
+fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True, width_ratios=(3, 1))
+fig.suptitle("Slopes along x")
+ax[0].imshow(
+    st_slope3d_y[ny // 2].T,
+    aspect="auto",
+    cmap="RdBu_r",
+    vmin=-clip_sy,
+    vmax=clip_sy,
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+ax[0].set_ylabel("Time [s]")
+ax[0].set_xlabel("X [km]")
+im1 = ax[1].imshow(
+    st_slope3d_y[:, nx // 2].T,
+    aspect="auto",
+    cmap="RdBu_r",
+    vmin=-clip_sy,
+    vmax=clip_sy,
+    extent=(y[0], y[-1], t[-1], t[0]),
+)
+ax[1].set_xlabel("Y [km]")
+create_colorbar(fig, ax[1], im1)
+fig.tight_layout()
+
+fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True, width_ratios=(3, 1))
+fig.suptitle("Slopes along y")
+ax[0].imshow(
+    st_slope3d_x[ny // 2].T,
+    aspect="auto",
+    cmap="RdBu_r",
+    vmin=-clip_sx,
+    vmax=clip_sx,
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+ax[0].set_ylabel("Time [s]")
+ax[0].set_xlabel("X [km]")
+im1 = ax[1].imshow(
+    st_slope3d_x[:, nx // 2].T,
+    aspect="auto",
+    cmap="RdBu_r",
+    vmin=-clip_sx,
+    vmax=clip_sx,
+    extent=(y[0], y[-1], t[-1], t[0]),
+)
+ax[1].set_xlabel("Y [km]")
+create_colorbar(fig, ax[1], im1)
+fig.tight_layout()
+
+
 ###############################################################################
 # Final considerations
 # --------------------
diff --git a/pylops/utils/_structuretensor.py b/pylops/utils/_structuretensor.py
new file mode 100644
index 00000000..37f3a4ec
--- /dev/null
+++ b/pylops/utils/_structuretensor.py
@@ -0,0 +1,266 @@
+import numpy as np
+
+from pylops.utils.backend import (
+    get_array_module,
+    get_gaussian_filter,
+)
+from pylops.utils.typing import NDArray
+
+
+def _structure_tensor_2d(
+    d: NDArray,
+    dz: float = 1.0,
+    dx: float = 1.0,
+    smooth: float = 5.0,
+    eps: float = 0.0,
+    dips: bool = False,
+) -> tuple[NDArray, NDArray]:
+    r"""2D Structure Tensor local slope estimation
+
+    Parameters
+    ----------
+    d : :obj:`numpy.ndarray`
+        Input dataset of size :math:`n_z \times n_x`
+    dz : :obj:`float`, optional
+        Sampling in :math:`z`-axis, :math:`\Delta z`
+
+        .. warning::
+            Since version 1.17.0, defaults to 1.0.
+
+    dx : :obj:`float`, optional
+        Sampling in :math:`x`-axis, :math:`\Delta x`
+
+        .. warning::
+            Since version 1.17.0, defaults to 1.0.
+
+    smooth : :obj:`float` or :obj:`numpy.ndarray`, optional
+        Standard deviation for Gaussian kernel. The standard deviations of the
+        Gaussian filter are given for each axis as a sequence, or as a single number,
+        in which case it is equal for all axes.
+
+        .. warning::
+            Default changed in version 1.17.0 to 5 from previous value of 20.
+
+    eps : :obj:`float`, optional
+        .. versionadded:: 1.17.0
+
+        Regularization term. All slopes where
+        :math:`|g_{zx}| < \epsilon \max_{(x, z)} \{|g_{zx}|, |g_{zz}|, |g_{xx}|\}`
+        are set to zero. All anisotropies where :math:`\lambda_\text{max} < \epsilon`
+        are also set to zero. See Notes. When using with small values of ``smooth``,
+        start from a very small number (e.g. 1e-10) and start increasing by a power
+        of 10 until results are satisfactory.
+
+    dips : :obj:`bool`, optional
+        .. versionadded:: 2.0.0
+
+        Return dips (``True``) instead of slopes (``False``).
+
+    Returns
+    -------
+    slopes : :obj:`numpy.ndarray`
+        Estimated local slopes. The unit is that of
+        :math:`\Delta z/\Delta x`.
+
+        .. warning::
+            Prior to version 1.17.0, always returned dips.
+
+    anisotropies : :obj:`numpy.ndarray`
+        Estimated local anisotropies: :math:`1-\lambda_\text{min}/\lambda_\text{max}`
+
+        .. note::
+            Since 1.17.0, changed name from ``linearity`` to ``anisotropies``.
+            Definition remains the same.
+
+    """
+    ncp = get_array_module(d)
+
+    gz, gx = ncp.gradient(d, dz, dx)
+    gzz, gzx, gxx = gz * gz, gz * gx, gx * gx
+
+    # smoothing
+    gaussian_filter = get_gaussian_filter(d)
+    gzz = gaussian_filter(gzz, sigma=smooth)
+    gzx = gaussian_filter(gzx, sigma=smooth)
+    gxx = gaussian_filter(gxx, sigma=smooth)
+
+    anisos = ncp.zeros_like(d)
+    gmax = max(gzz.max(), gxx.max(), ncp.abs(gzx).max())
+    if gmax <= eps:
+        return ncp.zeros_like(d), anisos
+
+    gzz /= gmax
+    gzx /= gmax
+    gxx /= gmax
+
+    lcommon1 = 0.5 * (gzz + gxx)
+    lcommon2 = 0.5 * ncp.sqrt((gzz - gxx) ** 2 + 4 * gzx**2)
+    l1 = lcommon1 + lcommon2
+    l2 = lcommon1 - lcommon2
+
+    regdata_aniso = l1 > eps
+    anisos[regdata_aniso] = 1 - l2[regdata_aniso] / l1[regdata_aniso]
+
+    if dips:
+        slopes = 0.5 * ncp.arctan2(2 * gzx, gzz - gxx)
+    else:
+        slopes = ncp.zeros_like(d)
+        regdata_slope = ncp.abs(gzx) > eps
+        slopes[regdata_slope] = (l1 - gzz)[regdata_slope] / gzx[regdata_slope]
+
+    return slopes, anisos
+
+
+def _structure_tensor_3d(
+    d: NDArray,
+    dy: float = 1.0,
+    dx: float = 1.0,
+    dz: float = 1.0,
+    smooth: float = 5.0,
+    eps: float = 0.0,
+    dips: bool = False,
+    anisotropies: bool = False,
+    batch_size: int | None = 1_000_000,
+) -> tuple[NDArray, NDArray] | tuple[NDArray, NDArray, NDArray, NDArray]:
+    r"""3D Structure Tensor local slope estimation
+
+    Parameters
+    ----------
+    d : :obj:`numpy.ndarray`
+        Input dataset of size :math:`n_y \times n_x \times n_z`
+    dy : :obj:`float`, optional
+        Sampling in :math:`y`-axis, :math:`\Delta y`
+    dx : :obj:`float`, optional
+        Sampling in :math:`x`-axis, :math:`\Delta x`
+    dz : :obj:`float`, optional
+        Sampling in :math:`z`-axis, :math:`\Delta z`
+    smooth : :obj:`float` or :obj:`numpy.ndarray`, optional
+        Standard deviation for Gaussian kernel. The standard deviations of the
+        Gaussian filter are given for each axis as a sequence, or as a single number,
+        in which case it is equal for all axes.
+    eps : :obj:`float`, optional
+        Regularization term.
+    dips : :obj:`bool`, optional
+        Return dips (``True``) instead of slopes (``False``).
+    anisotropies : :obj:`bool`, optional
+        Return local linearity and planarity measures (``True``) or not (``False``).
+    batch_size : :obj:`int`, optional
+        Number of grid points being processed together if ``dips==False``
+        and/or ``anisotropies=True``; this is done to avoid forming
+        the smoothed gradient-square tensor for all grid points at once
+        and computing the corresponding eigenvalues and eigenvectors.
+        If ``None``, operates on all points at once.
+
+    Returns
+    -------
+    slopes_x : :obj:`numpy.ndarray`
+        Estimated local slopes along the :math:`x`-axis
+    slopes_y : :obj:`numpy.ndarray`
+        Estimated local slopes along the :math:`y`-axis
+    linearities : :obj:`numpy.ndarray`
+        Local linearity measure
+    planarities : :obj:`numpy.ndarray`
+        Local planarity measure
+
+    """
+    ncp = get_array_module(d)
+
+    gy, gx, gz = ncp.gradient(d, dy, dx, dz)
+
+    gxx, gyy, gzz = gx * gx, gy * gy, gz * gz
+    gyx, gyz, gxz = gy * gx, gy * gz, gx * gz
+
+    # smoothing
+    gaussian_filter = get_gaussian_filter(d)
+    gxx = gaussian_filter(gxx, sigma=smooth)
+    gyy = gaussian_filter(gyy, sigma=smooth)
+    gzz = gaussian_filter(gzz, sigma=smooth)
+    gyx = gaussian_filter(gyx, sigma=smooth)
+    gyz = gaussian_filter(gyz, sigma=smooth)
+    gxz = gaussian_filter(gxz, sigma=smooth)
+
+    if dips:
+        slopes_x = (0.5 * ncp.arctan2(2 * gxz, gzz - gxx)).reshape(d.shape)
+        slopes_y = (0.5 * ncp.arctan2(2 * gyz, gzz - gyy)).reshape(d.shape)
+        if not anisotropies:
+            return slopes_x, slopes_y
+    else:
+        slopes_x = slopes_y = ncp.empty(0, dtype=d.dtype)  # needed for typing only
+
+    # batch calculation for structure tensor (needed when dips=False or anisotropies=True)
+    bsize = d.size if batch_size is None else min(int(batch_size), d.size)
+    batch_in = np.arange(0, d.size, bsize, dtype=np.int64)
+    batch_end = np.minimum(batch_in + bsize, d.size)
+
+    if not dips:
+        vy = ncp.empty(d.size, dtype=d.dtype)
+        vx = ncp.empty(d.size, dtype=d.dtype)
+        vz = ncp.empty(d.size, dtype=d.dtype)
+    else:
+        vy = vx = vz = ncp.empty(0, dtype=d.dtype)  # needed for typing only
+
+    if anisotropies or eps > 0:
+        regdata = ncp.zeros(d.size, dtype=bool)
+        l1 = ncp.empty(d.size, dtype=d.dtype)
+        l2 = ncp.empty(d.size, dtype=d.dtype)
+        l3 = ncp.empty(d.size, dtype=d.dtype)
+    else:
+        regdata = ncp.empty(0, dtype=bool)
+        l1 = l2 = l3 = ncp.empty(0, dtype=d.dtype)  # needed for typing only
+
+    # flatten smoothed gradient tensors for batch slicing
+    gyy_r, gxx_r, gzz_r = gyy.ravel(), gxx.ravel(), gzz.ravel()
+    gyx_r, gyz_r, gxz_r = gyx.ravel(), gyz.ravel(), gxz.ravel()
+
+    # compute eigenvalues/eigenvectors in batches
+    for b_in, b_end in zip(batch_in, batch_end, strict=True):
+        G = ncp.empty((b_end - b_in, 3, 3), dtype=d.dtype)
+        G[:, 0, 0], G[:, 0, 1], G[:, 0, 2] = (
+            gyy_r[b_in:b_end],
+            gyx_r[b_in:b_end],
+            gyz_r[b_in:b_end],
+        )
+        G[:, 1, 0], G[:, 1, 1], G[:, 1, 2] = (
+            gyx_r[b_in:b_end],
+            gxx_r[b_in:b_end],
+            gxz_r[b_in:b_end],
+        )
+        G[:, 2, 0], G[:, 2, 1], G[:, 2, 2] = (
+            gyz_r[b_in:b_end],
+            gxz_r[b_in:b_end],
+            gzz_r[b_in:b_end],
+        )
+
+        evalues, evectors = ncp.linalg.eigh(G)
+
+        if not dips:
+            # largest eigenvalue eigenvector is evectors[:, :, 2] (eigh sorts in ascending order)
+            vy[b_in:b_end] = evectors[:, 0, 2]
+            vx[b_in:b_end] = evectors[:, 1, 2]
+            vz[b_in:b_end] = -evectors[:, 2, 2]
+
+        if anisotropies or eps > 0:
+            l1[b_in:b_end] = evalues[:, 2]
+            l2[b_in:b_end] = evalues[:, 1]
+            l3[b_in:b_end] = evalues[:, 0]
+            regdata[b_in:b_end] = l1[b_in:b_end] > eps
+
+    if not dips:
+        slopes_x = -(vx / vz).reshape(d.shape)
+        slopes_y = -(vy / vz).reshape(d.shape)
+
+    if anisotropies:
+        linearity = ncp.zeros(d.size, dtype=d.dtype)
+        planarity = ncp.zeros(d.size, dtype=d.dtype)
+
+        linearity[regdata] = 1 - l2[regdata] / l1[regdata]
+        planarity[regdata] = (l2[regdata] - l3[regdata]) / l1[regdata]
+
+        return (
+            slopes_x,
+            slopes_y,
+            linearity.reshape(d.shape),
+            planarity.reshape(d.shape),
+        )
+
+    return slopes_x, slopes_y
diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py
index 2a0044e6..11b1e678 100644
--- a/pylops/utils/signalprocessing.py
+++ b/pylops/utils/signalprocessing.py
@@ -8,6 +8,7 @@
 
 import warnings
 from collections.abc import Sequence
+from typing import Literal, overload
 
 import numpy as np
 
@@ -15,11 +16,11 @@
 from pylops.optimization.leastsquares import preconditioned_inversion
 from pylops.utils._internal import _value_or_sized_to_tuple
 from pylops.utils._pwd2d import _conv_allpass, _triangular_smoothing_from_boxcars
+from pylops.utils._structuretensor import _structure_tensor_2d, _structure_tensor_3d
 from pylops.utils.backend import (
     get_array_module,
     get_csr_matrix,
     get_dia_matrix,
-    get_gaussian_filter,
     get_normalize_axis_index,
     get_toeplitz,
 )
@@ -95,7 +96,7 @@ def nonstationary_convmtx(
     H: NDArray,
     n: int,
     hc: int = 0,
-    pad: tuple[int] = (0, 0),
+    pad: tuple[int, ...] = (0, 0),
     sparse: bool = False,
 ) -> NDArray:
     r"""Convolution matrix from a bank of filters
@@ -150,14 +151,54 @@ def nonstationary_convmtx(
     return C
 
 
+@overload
 def slope_estimate(
     d: NDArray,
     dz: float = 1.0,
     dx: float = 1.0,
+    dy: None = None,
     smooth: float = 5.0,
     eps: float = 0.0,
     dips: bool = False,
-) -> tuple[NDArray, NDArray]:
+    anisotropies: Literal[False] | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[NDArray, NDArray]: ...
+@overload
+def slope_estimate(
+    d: NDArray,
+    dz: float,
+    dx: float,
+    dy: float | None = None,
+    smooth: float = 5.0,
+    eps: float = 0.0,
+    dips: bool = False,
+    anisotropies: Literal[False] | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[tuple[NDArray, NDArray], None]: ...
+@overload
+def slope_estimate(
+    d: NDArray,
+    dz: float = 1.0,
+    dx: float = 1.0,
+    dy: float | None = None,
+    smooth: float = 5.0,
+    eps: float = 0.0,
+    dips: bool = False,
+    *,
+    anisotropies: Literal[True],
+    batch_size: int | None = 1_000_000,
+) -> tuple[tuple[NDArray, NDArray], tuple[NDArray, NDArray]]: ...
+def slope_estimate(
+    d: NDArray,
+    dz: float = 1.0,
+    dx: float = 1.0,
+    dy: float | None = None,
+    smooth: float = 5.0,
+    eps: float = 0.0,
+    dips: bool = False,
+    anisotropies: bool | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[NDArray | tuple[NDArray, NDArray], NDArray | tuple[NDArray, NDArray] | None]:
     r"""Local slope estimation
 
     Local slopes are estimated using the *Structure Tensor* algorithm [1]_.
@@ -170,7 +211,8 @@ def slope_estimate(
     Parameters
     ----------
     d : :obj:`numpy.ndarray`
-        Input dataset of size :math:`n_z \times n_x`
+        Input dataset of size :math:`n_z \times n_x` for 2d or
+        of size :math:`n_y \times n_x \times n_z` for 3d.
     dz : :obj:`float`, optional
         Sampling in :math:`z`-axis, :math:`\Delta z`
 
@@ -183,6 +225,11 @@ def slope_estimate(
         .. warning::
             Since version 1.17.0, defaults to 1.0.
 
+    dy : :obj:`float`, optional
+        .. versionadded:: 2.9.0
+
+        Sampling in :math:`y`-axis, :math:`\Delta y`. Defaults to 1.0 when ``d``
+        is 3d; ignored when ``d`` is 2d.
     smooth : :obj:`float` or :obj:`numpy.ndarray`, optional
         Standard deviation for Gaussian kernel. The standard deviations of the
         Gaussian filter are given for each axis as a sequence, or as a single number,
@@ -200,23 +247,39 @@ def slope_estimate(
         are also set to zero. See Notes. When using with small values of ``smooth``,
         start from a very small number (e.g. 1e-10) and start increasing by a power
         of 10 until results are satisfactory.
-
     dips : :obj:`bool`, optional
         .. versionadded:: 2.0.0
 
         Return dips (``True``) instead of slopes (``False``).
+    anisotropies : :obj:`bool`, optional
+        .. versionadded:: 2.9.0
+
+        Return anisotropies (``True``) or not (``False``). Ignored when ``d``
+        is 2d as anisotropies are always returned.
+    batch_size : :obj:`int`, optional
+        .. versionadded:: 2.9.0
+
+        Number of grid points being processed together if ``dips==False``
+        and/or ``anisotropies=True``; this is done to avoid forming
+        the smoothed gradient-square tensor for all grid points at once
+        and computing the corresponding eigenvalues and eigenvectors.
+        If ``None``, operates on all points at once.
 
     Returns
     -------
-    slopes : :obj:`numpy.ndarray`
-        Estimated local slopes. The unit is that of
-        :math:`\Delta z/\Delta x`.
+    slopes : :obj:`numpy.ndarray` or :obj:`tuple`
+        Estimated local slopes (in 2d) or set of local slopes
+        along :math:`y`-axis and :math:`x`-axis (in 3d). The unit
+        is that of :math:`\Delta z/\Delta x` (and :math:`\Delta z/\Delta y`).
 
         .. warning::
             Prior to version 1.17.0, always returned dips.
 
     anisotropies : :obj:`numpy.ndarray`
-        Estimated local anisotropies: :math:`1-\lambda_\text{min}/\lambda_\text{max}`
+        Estimated local linearities (:math:`1-\lambda_2/\lambda_1`)
+        (in 2d) or set of local linearities and planarities
+        (:math:`(\lambda_2-\lambda_3)/\lambda_1`) in 3d, where
+        :math:`\lambda_1 \ge \lambda_2 \ge \lambda_3`.
 
         .. note::
             Since 1.17.0, changed name from ``linearity`` to ``anisotropies``.
@@ -224,17 +287,17 @@ def slope_estimate(
 
     Notes
     -----
-    For each pixel of the input dataset :math:`\mathbf{d}` the local gradients
-    :math:`g_z = \frac{\partial \mathbf{d}}{\partial z}` and
+    In 2d, for each pixel of the input dataset :math:`\mathbf{d}`, the
+    local gradients :math:`g_z = \frac{\partial \mathbf{d}}{\partial z}` and
     :math:`g_x = \frac{\partial \mathbf{d}}{\partial x}` are computed
     and used to define the following three quantities:
 
     .. math::
-        \begin{align}
+        \begin{aligned}
         g_{zz} &= \left(\frac{\partial \mathbf{d}}{\partial z}\right)^2\\
         g_{xx} &= \left(\frac{\partial \mathbf{d}}{\partial x}\right)^2\\
         g_{zx} &= \frac{\partial \mathbf{d}}{\partial z}\cdot\frac{\partial \mathbf{d}}{\partial x}
-        \end{align}
+        \end{aligned}
 
     They are then spatially smoothed and at each pixel their smoothed versions are
     arranged in a :math:`2 \times 2` matrix called the *smoothed
@@ -251,9 +314,10 @@ def slope_estimate(
     :math:`p = \frac{\lambda_\text{max} - g_{zz}}{g_{zx}}`,
     where :math:`\lambda_\text{max}` is the largest eigenvalue of :math:`\mathbf{G}`.
 
-    Similarly, local dips can be expressed as :math:`\tan(2\theta) = 2g_{zx} / (g_{zz} - g_{xx})`.
+    Similarly, local dips can be expressed as
+    :math:`\tan(2\theta) = 2g_{zx} / (g_{zz} - g_{xx})`.
 
-    Moreover, we can obtain a measure of local anisotropy, defined as
+    Moreover, a measure of local anisotropy can be defined as
 
     .. math::
         a = 1-\lambda_\text{min}/\lambda_\text{max}
@@ -262,55 +326,89 @@ def slope_estimate(
     A value of :math:`a = 0`  indicates perfect isotropy whereas :math:`a = 1`
     indicates perfect anisotropy.
 
-    .. [1] Van Vliet, L. J.,  Verbeek, P. W., "Estimators for orientation and
-        anisotropy in digitized images", Journal ASCI Imaging Workshop. 1995.
-
-    """
-    ncp = get_array_module(d)
-
-    slopes = ncp.zeros_like(d)
-    anisos = ncp.zeros_like(d)
+    In 3d, the same procedure is applied to the
+    local gradients :math:`g_y = \frac{\partial \mathbf{d}}{\partial y}` and
+    :math:`g_x = \frac{\partial \mathbf{d}}{\partial x}` and
+    :math:`g_z = \frac{\partial \mathbf{d}}{\partial z}`, which form a
+    :math:`3 \times 3` *smoothed gradient-square tensor*.
 
-    gz, gx = ncp.gradient(d, dz, dx)
-    gzz, gzx, gxx = gz * gz, gz * gx, gx * gx
+    Local dips are computed as :math:`\tan(2\theta_x) = 2g_{zx} / (g_{zz} - g_{xx})`
+    and :math:`\tan(2\theta_y) = 2g_{zy} / (g_{zz} - g_{yy})`, whilst local
+    slopes are defined :math:`p_x = -\frac{v_x}{v_z}` and :math:`p_y = -\frac{v_y}{v_z}`,
+    where :math:`v_y`, :math:`v_x`, and :math:`v_z` are the components of the eigenvector
+    of `\mathbf{G}` associated with the largest eigenvalue.
 
-    # smoothing
-    gzz = get_gaussian_filter(d)(gzz, sigma=smooth)
-    gzx = get_gaussian_filter(d)(gzx, sigma=smooth)
-    gxx = get_gaussian_filter(d)(gxx, sigma=smooth)
+    Finally a measure of local linearity (same as anisotropy) is computed as
 
-    gmax = max(gzz.max(), gxx.max(), ncp.abs(gzx).max())
-    if gmax <= eps:
-        return ncp.zeros_like(d), anisos
+    .. math::
+        l = 1-\lambda_\text{min}/\lambda_\text{max}
 
-    gzz /= gmax
-    gzx /= gmax
-    gxx /= gmax
+    whilst a measure of local planarity is computed as
 
-    lcommon1 = 0.5 * (gzz + gxx)
-    lcommon2 = 0.5 * ncp.sqrt((gzz - gxx) ** 2 + 4 * gzx**2)
-    l1 = lcommon1 + lcommon2
-    l2 = lcommon1 - lcommon2
+    .. math::
+        l = (\lambda_2-\lambda_3)/\lambda_1
 
-    regdata = l1 > eps
-    anisos[regdata] = 1 - l2[regdata] / l1[regdata]
+    .. [1] Van Vliet, L. J.,  Verbeek, P. W., "Estimators for orientation and
+        anisotropy in digitized images", Journal ASCI Imaging Workshop. 1995.
 
-    if dips:
-        slopes = 0.5 * ncp.arctan2(2 * gzx, gzz - gxx)
-    else:
-        regdata = ncp.abs(gzx) > eps
-        slopes[regdata] = (l1 - gzz)[regdata] / gzx[regdata]
+    """
+    if d.ndim == 2:
+        return _structure_tensor_2d(d, dz, dx, smooth, eps, dips)
 
-    return slopes, anisos
+    dy_3d = 1.0 if dy is None else dy
+    anisotropies_3d = bool(anisotropies)
+    outs = _structure_tensor_3d(
+        d, dy_3d, dx, dz, smooth, eps, dips, anisotropies_3d, batch_size
+    )
+    slopes_3d = (outs[0], outs[1])
+    anisos_3d = (outs[2], outs[3]) if anisotropies_3d and len(outs) == 4 else None
+    return slopes_3d, anisos_3d
 
 
+@overload
 def dip_estimate(
     d: NDArray,
     dz: float = 1.0,
     dx: float = 1.0,
+    dy: None = None,
     smooth: int = 5,
     eps: float = 0.0,
-) -> tuple[NDArray, NDArray]:
+    anisotropies: Literal[False] | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[NDArray, NDArray]: ...
+@overload
+def dip_estimate(
+    d: NDArray,
+    dz: float,
+    dx: float,
+    dy: float | None = None,
+    smooth: int = 5,
+    eps: float = 0.0,
+    anisotropies: Literal[False] | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[tuple[NDArray, NDArray], None]: ...
+@overload
+def dip_estimate(
+    d: NDArray,
+    dz: float,
+    dx: float,
+    dy: float | None = None,
+    smooth: int = 5,
+    eps: float = 0.0,
+    *,
+    anisotropies: Literal[True],
+    batch_size: int | None = 1_000_000,
+) -> tuple[tuple[NDArray, NDArray], tuple[NDArray, NDArray]]: ...
+def dip_estimate(
+    d: NDArray,
+    dz: float = 1.0,
+    dx: float = 1.0,
+    dy: float | None = None,
+    smooth: int = 5,
+    eps: float = 0.0,
+    anisotropies: bool | None = None,
+    batch_size: int | None = 1_000_000,
+) -> tuple[NDArray | tuple[NDArray, NDArray], NDArray | tuple[NDArray, NDArray] | None]:
     r"""Local dip estimation
 
     Local dips are estimated using the *Structure Tensor* algorithm [1]_.
@@ -326,6 +424,11 @@ def dip_estimate(
         Sampling in :math:`z`-axis, :math:`\Delta z`
     dx : :obj:`float`, optional
         Sampling in :math:`x`-axis, :math:`\Delta x`
+    dy : :obj:`float`, optional
+        .. versionadded:: 2.9.0
+
+        Sampling in :math:`y`-axis, :math:`\Delta y`. Defaults to 1.0 when ``d``
+        is 3d; ignored when ``d`` is 2d.
     smooth : :obj:`float` or :obj:`numpy.ndarray`, optional
         Standard deviation for Gaussian kernel. The standard deviations of the
         Gaussian filter are given for each axis as a sequence, or as a single number,
@@ -335,6 +438,19 @@ def dip_estimate(
         are also set to zero. See Notes. When using with small values of ``smooth``,
         start from a very small number (e.g. 1e-10) and start increasing by a power
         of 10 until results are satisfactory.
+    anisotropies : :obj:`bool`, optional
+        .. versionadded:: 2.9.0
+
+        Return anisotropies (``True``) or not (``False``). Ignored when ``d``
+        is 2d as anisotropies are always returned.
+    batch_size : :obj:`int`, optional
+        .. versionadded:: 2.9.0
+
+        Number of grid points being processed together if ``dips==False``
+        and/or ``anisotropies=True``; this is done to avoid forming
+        the smoothed gradient-square tensor for all grid points at once
+        and computing the corresponding eigenvalues and eigenvectors.
+        If ``None``, operates on all points at once.
 
     Returns
     -------
@@ -342,7 +458,10 @@ def dip_estimate(
         Estimated local dips. The unit is radians,
         in the range of :math:`-\frac{\pi}{2}` to :math:`\frac{\pi}{2}`.
     anisotropies : :obj:`numpy.ndarray`
-        Estimated local anisotropies: :math:`1-\lambda_\text{min}/\lambda_\text{max}`
+        Estimated local linearities (:math:`1-\lambda_2/\lambda_1`)
+        (in 2d) or set of local linearities and planarities
+        (:math:`(\lambda_2-\lambda_3)/\lambda_1`) in 3d, where
+        :math:`\lambda_1 \ge \lambda_2 \ge \lambda_3`.
 
     Notes
     -----
@@ -353,7 +472,17 @@ def dip_estimate(
         anisotropy in digitized images", Journal ASCI Imaging Workshop. 1995.
 
     """
-    dips, anisos = slope_estimate(d, dz=dz, dx=dx, smooth=smooth, eps=eps, dips=True)
+    dips, anisos = slope_estimate(
+        d,
+        dz=dz,
+        dx=dx,
+        dy=dy,
+        smooth=smooth,
+        eps=eps,
+        dips=True,
+        anisotropies=anisotropies,
+        batch_size=batch_size,
+    )
     return dips, anisos
 
 
diff --git a/pytests/test_signalutils.py b/pytests/test_signalutils.py
index 5feffbbf..6324ef9e 100644
--- a/pytests/test_signalutils.py
+++ b/pytests/test_signalutils.py
@@ -47,6 +47,40 @@
 np.random.seed(10)
 
 
+def _plane_wave_2d(x, y, f, c, theta):
+    """2D Plane wave modelling"""
+    # Define x-y grid
+    Y, X = np.meshgrid(y, x, indexing="ij")
+
+    # Slowness vector
+    p = (np.cos(np.deg2rad(theta)) / c, np.sin(np.deg2rad(theta)) / c)
+
+    # Construct plane wave
+    pw = np.exp(-1j * (2 * np.pi * f * (-(p[0] * Y + p[1] * X))))
+    pw = np.real(pw)
+
+    return pw
+
+
+def _plane_wave_3d(y, x, z, f, c, theta, phi):
+    """2D Plane wave modelling"""
+    # Define y-x-z grid
+    Y, X, Z = np.meshgrid(y, x, z, indexing="ij")
+
+    # Slowness vector
+    p = (
+        np.sin(np.deg2rad(theta)) * np.cos(np.deg2rad(phi)) / c,
+        np.sin(np.deg2rad(theta)) * np.sin(np.deg2rad(phi)) / c,
+        np.cos(np.deg2rad(theta)) / c,
+    )  # slowness vector
+
+    # Construct plane wave
+    pw = np.exp(-1j * (2 * np.pi * f * (-(p[0] * Y + p[1] * X + p[2] * Z))))
+    pw = np.real(pw)
+
+    return pw
+
+
 @pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)])
 @pytest.mark.parametrize("sparse", [False, True])
 def test_convmtx(par, sparse):
@@ -112,7 +146,65 @@ def test_nonstationary_convmtx(par, sparse):
     assert_array_almost_equal(y, y1, decimal=4)
 
 
-def test_slope_estimation_dips():
+@pytest.mark.parametrize("angle", [-45, -20, 0, 20, 45])
+def test_slope_estimation_analytical_2d(angle):
+    """Slope estimation using the Structure tensor algorithm for
+    2D plane wave - test against analytical solution."""
+
+    # Define x and y axes
+    ox, dx, nx = 0, 5, 101
+    oy, dy, ny = 0, 5, 101
+    x, y = np.arange(nx) * dx + ox, np.arange(ny) * dy + oy
+
+    # Compute plane wave
+    f = 10  # frequency
+    c = 1500  # Velocity
+    pw = _plane_wave_2d(x, y, f, c, angle)
+
+    # Slopes
+    slopes, _ = slope_estimate(
+        pw,
+        smooth=11,
+        eps=0.0,
+        dips=False,
+        anisotropies=False,
+    )
+
+    assert_array_almost_equal(np.median(slopes), np.tan(np.deg2rad(angle)), decimal=2)
+
+
+@pytest.mark.parametrize("angle", [-45, -20, 0, 20, 45])
+def test_slope_estimation_analytical_3d(angle):
+    """Slope estimation using the Structure tensor algorithm for
+    3D plane wave - test against analytical solution."""
+
+    # Define x and y axes
+    oy, dy, ny = 0, 5, 21
+    ox, dx, nx = 0, 5, 51
+    oz, dz, nz = 0, 5, 51
+    y, x, z = np.arange(ny) * dy + oy, np.arange(nx) * dx + ox, np.arange(nz) * dz + oz
+
+    # Compute plane wave
+    f = 10  # frequency
+    c = 1500  # Velocity
+    pw = _plane_wave_3d(y, x, z, f, c, angle, phi=0.0)
+
+    # Slopes
+    slopes, _ = slope_estimate(
+        pw,
+        dy=1.0,
+        smooth=11,
+        eps=0.0,
+        dips=False,
+        anisotropies=False,
+    )
+
+    assert_array_almost_equal(
+        np.median(slopes[1]), np.tan(np.deg2rad(angle)), decimal=2
+    )
+
+
+def test_slope_estimation_reg():
     """Slope estimation using the Structure tensor algorithm should
     apply regularisation (some slopes are set to zero)
     while dips should not use regularisation."""