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..fefd9ca4
--- /dev/null
+++ b/pylops/utils/_structuretensor.py
@@ -0,0 +1,263 @@
+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)
+
+    slopes = ncp.zeros_like(d)
+    anisos = ncp.zeros_like(d)
+
+    gz, gx = ncp.gradient(d, dz, dx)
+    gzz, gzx, gxx = gz * gz, gz * gx, gx * gx
+
+    # 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)
+
+    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 = l1 > eps
+    anisos[regdata] = 1 - l2[regdata] / l1[regdata]
+
+    if dips:
+        slopes = 0.5 * ncp.arctan2(2 * gzx, gzz - gxx)
+    else:
+        regdata = ncp.abs(gzx) > eps
+        slopes[regdata] = (l1 - gzz)[regdata] / gzx[regdata]
+
+    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, gzz = gx * gx, gz * gz
+    gyx, gyz, gxz = gy * gx, gy * gz, gx * gz
+
+    # smoothing
+    gxx = get_gaussian_filter(d)(gxx, sigma=smooth)
+    gzz = get_gaussian_filter(d)(gzz, sigma=smooth)
+    gyx = get_gaussian_filter(d)(gyx, sigma=smooth)
+    gyz = get_gaussian_filter(d)(gyz, sigma=smooth)
+    gxz = get_gaussian_filter(d)(gxz, sigma=smooth)
+
+    if not dips or anisotropies:
+        # additional paramets (not needed for dips)
+        gyy = gy * gy
+        gyy = get_gaussian_filter(d)(gyy, sigma=smooth)
+
+        # define batches
+        batch_size = int(batch_size)
+        ngrid_points = d.size
+        if batch_size is None or batch_size > ngrid_points:
+            batch_size = ngrid_points
+
+        batch_in = np.arange(0, ngrid_points, batch_size, dtype=np.int64)
+        batch_end = batch_in + batch_size
+        batch_end[-1] = min(batch_end[-1], ngrid_points)
+
+        # instantiated objects
+        vy = ncp.empty(d.size, dtype=d.dtype)
+        vx = ncp.empty(d.size, dtype=d.dtype)
+        vz = ncp.empty(d.size, dtype=d.dtype)
+
+        regdata = ncp.zeros(d.size, dtype=bool)
+        if anisotropies:
+            l1 = ncp.empty(d.size, dtype=d.dtype)
+            l2 = ncp.empty(d.size, dtype=d.dtype)
+            l3 = ncp.empty(d.size, dtype=d.dtype)
+
+        # compute eigenvalues/eigenvectors
+        for b_in, b_end in zip(batch_in, batch_end, strict=True):
+            # create matrices of second-order derivatives
+            G = ncp.empty(((b_end - b_in), 3, 3), dtype=d.dtype)
+
+            G[:, 0, 0] = gyy.ravel()[b_in:b_end]
+            G[:, 0, 1] = gyx.ravel()[b_in:b_end]
+            G[:, 0, 2] = gyz.ravel()[b_in:b_end]
+
+            G[:, 1, 0] = gyx.ravel()[b_in:b_end]
+            G[:, 1, 1] = gxx.ravel()[b_in:b_end]
+            G[:, 1, 2] = gxz.ravel()[b_in:b_end]
+
+            G[:, 2, 0] = gyz.ravel()[b_in:b_end]
+            G[:, 2, 1] = gxz.ravel()[b_in:b_end]
+            G[:, 2, 2] = gzz.ravel()[b_in:b_end]
+
+            evalues, evectors = ncp.linalg.eigh(G)
+
+            # extract the eigenvectors corresponding to the largest eigenvalue
+            idx = ncp.argmax(evalues, axis=1)
+            largest_evectors = evectors[np.arange(G.shape[0]), :, idx]
+
+            vy[b_in:b_end] = largest_evectors[:, 0]
+            vx[b_in:b_end] = largest_evectors[:, 1]
+            vz[b_in:b_end] = -largest_evectors[:, 2]
+
+            if anisotropies or eps > 0:
+                # re-order eigenvalues
+                evalues = ncp.sort(evalues, axis=1)
+                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 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]
+
+        linearity = linearity.reshape(d.shape)
+        planarity = planarity.reshape(d.shape)
+
+    if dips:
+        slopes_x = 0.5 * ncp.arctan2(2 * gxz, gzz - gxx)
+        slopes_y = 0.5 * ncp.arctan2(2 * gyz, gzz - gyy)
+        slopes_x = slopes_x.reshape(d.shape)
+        slopes_y = slopes_y.reshape(d.shape)
+    else:
+        slopes_x = -(vx / vz).reshape(d.shape)
+        slopes_y = -(vy / vz).reshape(d.shape)
+
+    if anisotropies:
+        return slopes_x, slopes_y, linearity, planarity
+    else:
+        return slopes_x, slopes_y
diff --git a/pylops/utils/signalprocessing.py b/pylops/utils/signalprocessing.py
index 2a0044e6..37a8221a 100644
--- a/pylops/utils/signalprocessing.py
+++ b/pylops/utils/signalprocessing.py
@@ -15,11 +15,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,
 )
@@ -154,10 +154,13 @@ 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,
-) -> tuple[NDArray, NDArray]:
+    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 +173,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 +187,11 @@ def slope_estimate(
         .. warning::
             Since version 1.17.0, defaults to 1.0.
 
+    dy : :obj:`float`, optional
+        .. versionadded:: 2.8.0
+
+        Sampling in :math:`y`-axis, :math:`\Delta y`. 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 +209,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.8.0
+
+        Return anisotropies (``True``) or not (``False``). Ignored when ``d``
+        is 2d as anisotropies are always returned.
+    batch_size : :obj:`int`, optional
+        .. versionadded:: 2.8.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:`y`-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,8 +249,8 @@ 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:
 
@@ -251,9 +276,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,45 +288,43 @@ 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)
+    """
+    if d.ndim == 2:
+        slopes, anisos = _structure_tensor_2d(d, dz, dx, smooth, eps, dips)
     else:
-        regdata = ncp.abs(gzx) > eps
-        slopes[regdata] = (l1 - gzz)[regdata] / gzx[regdata]
-
+        outs = _structure_tensor_3d(
+            d, dy, dx, dx, smooth, eps, dips, anisotropies, batch_size
+        )
+        slopes = (outs[0], outs[1])
+        if anisotropies:
+            anisos = (outs[2], outs[3])
+        else:
+            anisos = None
     return slopes, anisos
 
 
@@ -308,8 +332,10 @@ def dip_estimate(
     d: NDArray,
     dz: float = 1.0,
     dx: float = 1.0,
+    dy: float | None = None,
     smooth: int = 5,
     eps: float = 0.0,
+    batch_size: int | None = 1_000_000,
 ) -> tuple[NDArray, NDArray]:
     r"""Local dip estimation
 
@@ -326,6 +352,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.8.0
+
+        Sampling in :math:`y`-axis, :math:`\Delta y`. 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 +366,14 @@ 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.
+    batch_size : :obj:`int`, optional
+        .. versionadded:: 2.8.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
     -------
@@ -353,7 +392,9 @@ 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, 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."""