.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/plot_sigmoid_disks.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_plot_sigmoid_disks.py: 6. Multiscale Local Directions ============================== This example shows how to use the Curvelet transform to visualize local, multiscale preferrential directions in an image. Inspired by `Kymatio's Scattering disks `__. .. GENERATED FROM PYTHON SOURCE LINES 8-10 .. code-block:: default # sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 11-27 .. code-block:: default import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import numpy.typing as npt from mpl_toolkits.axes_grid1 import make_axes_locatable from pylops.signalprocessing import FFT2D from curvelops import FDCT2D from curvelops.plot import ( create_axes_grid, create_inset_axes_grid, overlay_arrows, overlay_disks, ) from curvelops.utils import array_split_nd, ndargmax .. GENERATED FROM PYTHON SOURCE LINES 28-30 Input ===== .. GENERATED FROM PYTHON SOURCE LINES 32-41 .. code-block:: default inputfile = "../testdata/sigmoid.npz" data = np.load(inputfile) data = data["sigmoid"] nx, nz = data.shape dx, dz = 0.005, 0.004 x, z = np.arange(nx) * dx, np.arange(nz) * dz .. GENERATED FROM PYTHON SOURCE LINES 42-57 .. code-block:: default aspect = dz / dx figsize_aspect = aspect * nz / nx opts_space = dict( extent=(x[0], x[-1], z[-1], z[0]), cmap="gray", interpolation="lanczos", aspect=aspect, ) vmax = 0.5 * np.max(np.abs(data)) fig, ax = plt.subplots(figsize=(8, figsize_aspect * 8)) ax.imshow(data.T, vmin=-vmax, vmax=vmax, **opts_space) ax.set(xlabel="Position [km]", ylabel="Depth [km]", title="Data") fig.tight_layout() .. image-sg:: /gallery/images/sphx_glr_plot_sigmoid_disks_001.png :alt: Data :srcset: /gallery/images/sphx_glr_plot_sigmoid_disks_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-60 Understanding Curvelet Disks ============================ .. GENERATED FROM PYTHON SOURCE LINES 62-63 First we create and apply curvelet transform. .. GENERATED FROM PYTHON SOURCE LINES 63-66 .. code-block:: default Cop = FDCT2D(data.shape, nbscales=4, nbangles_coarse=8, allcurvelets=False) d_c = Cop.struct(Cop @ data) .. GENERATED FROM PYTHON SOURCE LINES 67-76 Each wedge is mapped to a region of the scattering disk. The first number refers to the scale, the second to the wedge index, zero-indexed. The disks have the most energy in the direction perpendicular to the directions of minimum change. The following disk is computed with the entire image. We observe that with energy mostly along the top-bottom direction, the directions in the image will be mostly along the left-right direction, which matches the input data. .. GENERATED FROM PYTHON SOURCE LINES 76-86 .. code-block:: default rows, cols = 1, 1 fig, axes = create_axes_grid( rows, cols, kwargs_subplots=dict(projection="polar"), kwargs_figure=dict(figsize=(4, 4)), ) overlay_disks(d_c, axes, annotate=True) .. image-sg:: /gallery/images/sphx_glr_plot_sigmoid_disks_002.png :alt: plot sigmoid disks :srcset: /gallery/images/sphx_glr_plot_sigmoid_disks_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 87-93 Multiscale Local Directions ============================ The power of the curvelet transform is to provide dip information varying with location and scale. Below we will compute preferrential local directions using an approach based on the 2D FFT that does not differentiate between scales. .. GENERATED FROM PYTHON SOURCE LINES 95-126 .. code-block:: default rows, cols = 5, 6 def local_single_scale_dips(data: npt.NDArray, rows: int, cols: int) -> npt.NDArray: kvecs = np.empty((rows, cols, 2)) d_split = array_split_nd(data.T, rows, cols) for irow in range(kvecs.shape[0]): for icol in range(kvecs.shape[1]): d_loc = d_split[irow][icol].T Fop_loc = FFT2D( d_loc.shape, sampling=[dx, dz], norm="ortho", real=False, ifftshift_before=True, fftshift_after=True, engine="scipy", ) d_k_loc = Fop_loc @ d_loc kx_loc = Fop_loc.f1 kz_loc = Fop_loc.f2 kx_locmax, kz_locmax = ndargmax(np.abs(d_k_loc[:, kz_loc > 0])) k = np.array([kx_loc[kx_locmax], kz_loc[kz_loc > 0][kz_locmax]]) kvecs[irow, icol, :] = k / np.linalg.norm(k) return kvecs .. GENERATED FROM PYTHON SOURCE LINES 127-158 .. code-block:: default diskcmap = "turbo" rows, cols = 5, 6 kvecs = local_single_scale_dips(data, rows, cols) kvecs *= 0.15 * min(x[-1] - x[0], z[-1] - z[0]) fig, ax = plt.subplots(figsize=(8, figsize_aspect * 8)) ax.imshow(data.T, vmin=-vmax, vmax=vmax, **opts_space) ax.set(xlabel="Position [km]", ylabel="Depth [km]") divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.1) mpl.colorbar.ColorbarBase( cax, cmap=plt.get_cmap(diskcmap), norm=mpl.colors.Normalize(vmin=0, vmax=1), alpha=0.8, ) # Local single-scale directions overlay_arrows(kvecs, ax) # Local multsicale directions axesin = create_inset_axes_grid( ax, rows, cols, height=0.6, width=0.6, kwargs_inset_axes=dict(projection="polar"), ) overlay_disks(d_c, axesin, linewidth=0.0, cmap=diskcmap) fig.tight_layout() .. image-sg:: /gallery/images/sphx_glr_plot_sigmoid_disks_003.png :alt: plot sigmoid disks :srcset: /gallery/images/sphx_glr_plot_sigmoid_disks_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 6.142 seconds) .. _sphx_glr_download_gallery_plot_sigmoid_disks.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_sigmoid_disks.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_sigmoid_disks.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_