pylops.utils.signalprocessing.slope_estimate

pylops.utils.signalprocessing.slope_estimate(d, dz=1.0, dx=1.0, dy=None, smooth=5.0, eps=0.0, dips=False, anisotropies=None, batch_size=1000000)

Local slope estimation

Local slopes are estimated using the Structure Tensor algorithm [1]. Slopes are returned as \(\tan\theta\), defined in a RHS coordinate system with \(z\)-axis pointing upward.

Note

For stability purposes, it is important to ensure that the orders of magnitude of the samplings are similar.

Parameters:

dnumpy.ndarray

Input dataset of size \(n_z \times n_x\) for 2d or of size \(n_y \times n_x \times n_z\) for 3d.

dzfloat, optional

Sampling in \(z\)-axis, \(\Delta z\)

Warning

Since version 1.17.0, defaults to 1.0.

dxfloat, optional

Sampling in \(x\)-axis, \(\Delta x\)

Warning

Since version 1.17.0, defaults to 1.0.

dyfloat, optional

Added in version 2.9.0.

Sampling in \(y\)-axis, \(\Delta y\). Defaults to 1.0 when d is 3d; ignored when d is 2d.

smoothfloat or 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.

epsfloat, optional

Added in version 1.17.0.

Regularization term. All slopes where \(|g_{zx}| < \epsilon \max_{(x, z)} \{|g_{zx}|, |g_{zz}|, |g_{xx}|\}\) are set to zero. All anisotropies where \(\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.

dipsbool, optional

Added in version 2.0.0.

Return dips (True) instead of slopes (False).

anisotropiesbool, optional

Added in version 2.9.0.

Return anisotropies (True) or not (False). Ignored when d is 2d as anisotropies are always returned.

batch_sizeint, optional

Added in version 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:

slopesnumpy.ndarray or tuple

Estimated local slopes (in 2d) or set of local slopes along \(y\)-axis and \(x\)-axis (in 3d). The unit is that of \(\Delta z/\Delta x\) (and \(\Delta z/\Delta y\)).

Warning

Prior to version 1.17.0, always returned dips.

anisotropiesnumpy.ndarray

Estimated local linearities (\(1-\lambda_2/\lambda_1\)) (in 2d) or set of local linearities and planarities (\((\lambda_2-\lambda_3)/\lambda_1\)) in 3d, where \(\lambda_1 \ge \lambda_2 \ge \lambda_3\).

Note

Since 1.17.0, changed name from linearity to anisotropies. Definition remains the same.

Notes

In 2d, for each pixel of the input dataset \(\mathbf{d}\), the local gradients \(g_z = \frac{\partial \mathbf{d}}{\partial z}\) and \(g_x = \frac{\partial \mathbf{d}}{\partial x}\) are computed and used to define the following three quantities:

\[\begin{split}\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{aligned}\end{split}\]

They are then spatially smoothed and at each pixel their smoothed versions are arranged in a \(2 \times 2\) matrix called the smoothed gradient-square tensor:

\[\begin{split}\mathbf{G} = \begin{bmatrix} g_{zz} & g_{zx} \\ g_{zx} & g_{xx} \end{bmatrix}\end{split}\]

Local slopes can be expressed as \(p = \frac{\lambda_\text{max} - g_{zz}}{g_{zx}}\), where \(\lambda_\text{max}\) is the largest eigenvalue of \(\mathbf{G}\).

Similarly, local dips can be expressed as \(\tan(2\theta) = 2g_{zx} / (g_{zz} - g_{xx})\).

Moreover, a measure of local anisotropy can be defined as

\[a = 1-\lambda_\text{min}/\lambda_\text{max}\]

where \(\lambda_\text{min}\) is the smallest eigenvalue of \(\mathbf{G}\). A value of \(a = 0\) indicates perfect isotropy whereas \(a = 1\) indicates perfect anisotropy.

In 3d, the same procedure is applied to the local gradients \(g_y = \frac{\partial \mathbf{d}}{\partial y}\) and \(g_x = \frac{\partial \mathbf{d}}{\partial x}\) and \(g_z = \frac{\partial \mathbf{d}}{\partial z}\), which form a \(3 \times 3\) smoothed gradient-square tensor.

Local dips are computed as \(\tan(2\theta_x) = 2g_{zx} / (g_{zz} - g_{xx})\) and \(\tan(2\theta_y) = 2g_{zy} / (g_{zz} - g_{yy})\), whilst local slopes are defined \(p_x = -\frac{v_x}{v_z}\) and \(p_y = -\frac{v_y}{v_z}\), where \(v_y\), \(v_x\), and \(v_z\) are the components of the eigenvector of mathbf{G} associated with the largest eigenvalue.

Finally a measure of local linearity (same as anisotropy) is computed as

\[l = 1-\lambda_\text{min}/\lambda_\text{max}\]

whilst a measure of local planarity is computed as

\[l = (\lambda_2-\lambda_3)/\lambda_1\]

[1]

Van Vliet, L. J., Verbeek, P. W., “Estimators for orientation and anisotropy in digitized images”, Journal ASCI Imaging Workshop. 1995.

Examples using pylops.utils.signalprocessing.slope_estimate

PWD-based slope estimation and structural smoothing

PWD-based slope estimation and structural smoothing

Seislet transform

Seislet transform

Slope estimation via Structure Tensor algorithm

Slope estimation via Structure Tensor algorithm