Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

incorrect check on dips argument in in `slope_estimation #572

Closed
Beramos opened this issue Mar 4, 2024 · 2 comments · Fixed by #573
Closed

incorrect check on dips argument in in `slope_estimation #572

Beramos opened this issue Mar 4, 2024 · 2 comments · Fixed by #573

Comments

@Beramos
Copy link
Contributor

Beramos commented Mar 4, 2024

The documentation states the following:

image

However, the logic flow in the code is the opposite:

if not dips:
slopes = 0.5 * np.arctan2(2 * gzx, gzz - gxx)
else:
regdata = np.abs(gzx) > eps
slopes[regdata] = (l1 - gzz)[regdata] / gzx[regdata]

I double checked with the original implementation of dips_estimate

def dip_estimate(d, dz=1.0, dx=1.0, smooth=5, eps=0.0):
r"""Local dip estimation
Local dips are estimated using the *Structure Tensor* algorithm [1]_.
.. note:: For stability purposes, it is important to ensure that the orders
of magnitude of the samplings are similar.
Parameters
----------
d : :obj:`np.ndarray`
Input dataset of size :math:`n_z \times n_x`
dz : :obj:`float`, optional
Sampling in :math:`z`-axis, :math:`\Delta z`
dx : :obj:`float`, optional
Sampling in :math:`x`-axis, :math:`\Delta x`
smooth : :obj:`float` or :obj:`np.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. 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.
Returns
-------
slopes : :obj:`np.ndarray`
Estimated local dips. The unit is radians,
in the range of :math:`-\frac{\pi}{2}` to :math:`\frac{\pi}{2}`.
anisotropies : :obj:`np.ndarray`
Estimated local anisotropies: :math:`1-\lambda_\text{min}/\lambda_\text{max}`
Notes
-----
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}
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}
They are then spatially smoothed and at each pixel their smoothed versions are
arranged in a :math:`2 \times 2` matrix called the *smoothed
gradient-square tensor*:
.. math::
\mathbf{G} =
\begin{bmatrix}
g_{zz} & g_{zx} \\
g_{zx} & g_{xx}
\end{bmatrix}
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
.. math::
a = 1-\lambda_\text{min}/\lambda_\text{max}
where :math:`\lambda_\text{min}` is the smallest eigenvalue of :math:`\mathbf{G}`.
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.
"""
anisos = np.zeros_like(d)
gz, gx = np.gradient(d, dz, dx)
gzz, gzx, gxx = gz * gz, gz * gx, gx * gx
# smoothing
gzz = gaussian_filter(gzz, sigma=smooth)
gzx = gaussian_filter(gzx, sigma=smooth)
gxx = gaussian_filter(gxx, sigma=smooth)
gmax = max(gzz.max(), gxx.max(), np.abs(gzx).max())
if gmax <= eps:
return np.zeros_like(d), anisos
gzz /= gmax
gzx /= gmax
gxx /= gmax
lcommon1 = 0.5 * (gzz + gxx)
lcommon2 = 0.5 * np.sqrt((gzz - gxx) ** 2 + 4 * gzx**2)
l1 = lcommon1 + lcommon2
l2 = lcommon1 - lcommon2
regdata = l1 > eps
anisos[regdata] = 1 - l2[regdata] / l1[regdata]
dips = 0.5 * np.arctan2(2 * gzx, gzz - gxx)
return dips, anisos

And it is inline with the current documentation. It appears that the bug was introduced when dips_estimate and slope_estimate` merged.

I have a PR ready for this.

@mrava87
Copy link
Collaborator

mrava87 commented Mar 8, 2024

Thanks @Beramos for noticing this.

@cako, since you did quite some changes to this part of the code, do you agree with this? If so, I can go ahead and merge the PR

@mrava87
Copy link
Collaborator

mrava87 commented Mar 16, 2024

@cako?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants