{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Slope estimation via Structure Tensor algorithm\n\nThis example shows how to estimate local slopes or local dips of a two-dimensional\narray using :py:func:`pylops.utils.signalprocessing.slope_estimate` and\n:py:func:`pylops.utils.signalprocessing.dip_estimate`.\n\nKnowing the local slopes of an image (or a seismic data) can be useful for\na variety of tasks in image (or geophysical) processing such as denoising,\nsmoothing, or interpolation. When slopes are used with the\n:py:class:`pylops.signalprocessing.Seislet` operator, the input dataset can be\ncompressed and the sparse nature of the Seislet transform can also be used to\nprecondition sparsity-promoting inverse problems.\n\nWe will show examples of a variety of different settings, including a comparison\nwith the original implementation in [1].\n\n.. [1] van Vliet, L. J., Verbeek, P. W., \"Estimators for orientation and\n anisotropy in digitized images\", Journal ASCI Imaging Workshop. 1995.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib.image import imread\nfrom matplotlib.ticker import FuncFormatter, MultipleLocator\nfrom mpl_toolkits.axes_grid1 import make_axes_locatable\n\nimport pylops\nfrom pylops.signalprocessing.seislet import _predict_trace\nfrom pylops.utils.signalprocessing import dip_estimate, slope_estimate\n\nplt.close(\"all\")\nnp.random.seed(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python logo\nTo start we import a 2d image and estimate the local dips of the image.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "im = np.load(\"../testdata/python.npy\")[..., 0]\nim = im / 255.0 - 0.5\n\nangles, anisotropy = dip_estimate(im, smooth=7)\nangles = -np.rad2deg(angles)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharex=True, sharey=True)\niax = axs[0].imshow(im, cmap=\"viridis\", origin=\"lower\")\naxs[0].set_title(\"Data\")\ncax = make_axes_locatable(axs[0]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncax.axis(\"off\")\n\niax = axs[1].imshow(angles, cmap=\"twilight_shifted\", origin=\"lower\", vmin=-90, vmax=90)\naxs[1].set_title(\"Angle of incline\")\ncax = make_axes_locatable(axs[1]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(\n iax,\n ticks=MultipleLocator(30),\n format=FuncFormatter(lambda x, pos: \"{:.0f}\u00b0\".format(x)),\n cax=cax,\n orientation=\"vertical\",\n)\n\niax = axs[2].imshow(anisotropy, cmap=\"Reds\", origin=\"lower\", vmin=0, vmax=1)\naxs[2].set_title(\"Anisotropy\")\ncax = make_axes_locatable(axs[2]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(iax, cax=cax, orientation=\"vertical\")\nfig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Seismic data\nWe can now repeat the same using some seismic data. We will first define\na single trace and a slope field, apply such slope field to the trace\nrecursively to create the other traces of the data and finally try to recover\nthe underlying slope field from the data alone.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Reflectivity model\nnx, nt = 2**7, 121\ndx, dt = 0.01, 0.004\nx, t = np.arange(nx) * dx, np.arange(nt) * dt\n\nnspike = nt // 8\nrefl = np.zeros(nt)\nit = np.sort(np.random.permutation(range(10, nt - 20))[:nspike])\nrefl[it] = np.random.normal(0.0, 1.0, nspike)\n\n# Wavelet\nntwav = 41\nf0 = 30\ntwav = np.arange(ntwav) * dt\nwav, *_ = pylops.utils.wavelets.ricker(twav, f0)\n\n# Input trace\ntrace = np.convolve(refl, wav, mode=\"same\")\n\n# Slopes\ntheta = np.deg2rad(np.linspace(0, 30, nx))\nslope = np.outer(np.ones(nt), np.tan(theta) * dt / dx)\n\n# Model data\nd = np.zeros((nt, nx))\ntr = trace.copy()\nfor ix in range(nx):\n tr = _predict_trace(tr, t, dt, dx, slope[:, ix])\n d[:, ix] = tr\n\n# Estimate slopes\nslope_est, _ = slope_estimate(d, dt, dx, smooth=10)\nslope_est *= -1" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(2, 2, figsize=(6, 6), sharex=True, sharey=True)\n\nopts = dict(aspect=\"auto\", extent=(x[0], x[-1], t[-1], t[0]))\niax = axs[0, 0].imshow(d, cmap=\"gray\", vmin=-1, vmax=1, **opts)\naxs[0, 0].set(title=\"Data\", ylabel=\"Time [s]\")\ncax = make_axes_locatable(axs[0, 0]).append_axes(\"right\", size=\"5%\", pad=0.05)\nfig.colorbar(iax, cax=cax, orientation=\"vertical\")\n\nopts.update(dict(cmap=\"cividis\", vmin=np.min(slope), vmax=np.max(slope)))\niax = axs[0, 1].imshow(slope, **opts)\naxs[0, 1].set(title=\"True Slope\")\ncax = make_axes_locatable(axs[0, 1]).append_axes(\"right\", size=\"5%\", pad=0.05)\nfig.colorbar(iax, cax=cax, orientation=\"vertical\")\ncax.set_ylabel(\"[s/km]\")\n\niax = axs[1, 0].imshow(np.abs(slope - slope_est), **opts)\naxs[1, 0].set(\n title=\"Estimate absolute error\", ylabel=\"Time [s]\", xlabel=\"Position [km]\"\n)\ncax = make_axes_locatable(axs[1, 0]).append_axes(\"right\", size=\"5%\", pad=0.05)\nfig.colorbar(iax, cax=cax, orientation=\"vertical\")\ncax.set_ylabel(\"[s/km]\")\n\niax = axs[1, 1].imshow(slope_est, **opts)\naxs[1, 1].set(title=\"Estimated Slope\", xlabel=\"Position [km]\")\ncax = make_axes_locatable(axs[1, 1]).append_axes(\"right\", size=\"5%\", pad=0.05)\nfig.colorbar(iax, cax=cax, orientation=\"vertical\")\ncax.set_ylabel(\"[s/km]\")\n\nfig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Concentric circles\nThe original paper by van Vliet and Verbeek [1] has an example with concentric\ncircles. We recover their original images and compare our implementation with\ntheirs.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def rgb2gray(rgb):\n return np.dot(rgb[..., :3], [0.2989, 0.5870, 0.1140])\n\n\ncircles_input = rgb2gray(imread(\"../testdata/slope_estimate/concentric.png\"))\ncircles_angles = rgb2gray(imread(\"../testdata/slope_estimate/concentric_angles.png\"))\n\nangles, anisos_sm0 = dip_estimate(circles_input, smooth=0)\nangles_sm0 = np.rad2deg(angles)\n\nangles, anisos_sm4 = dip_estimate(circles_input, smooth=4)\nangles_sm4 = np.rad2deg(angles)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(2, 3, figsize=(6, 4), sharex=True, sharey=True)\naxs[0, 0].imshow(circles_input, cmap=\"gray\", aspect=\"equal\")\naxs[0, 0].set(title=\"Original Image\")\ncax = make_axes_locatable(axs[0, 0]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncax.axis(\"off\")\n\naxs[1, 0].imshow(-circles_angles, cmap=\"twilight_shifted\")\naxs[1, 0].set(title=\"Original Angles\")\ncax = make_axes_locatable(axs[1, 0]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncax.axis(\"off\")\n\nim = axs[0, 1].imshow(angles_sm0, cmap=\"twilight_shifted\", vmin=-90, vmax=90)\ncax = make_axes_locatable(axs[0, 1]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(\n im,\n ticks=MultipleLocator(30),\n format=FuncFormatter(lambda x, pos: \"{:.0f}\u00b0\".format(x)),\n cax=cax,\n orientation=\"vertical\",\n)\naxs[0, 1].set(title=\"Angles (smooth=0)\")\n\nim = axs[1, 1].imshow(angles_sm4, cmap=\"twilight_shifted\", vmin=-90, vmax=90)\ncax = make_axes_locatable(axs[1, 1]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(\n im,\n ticks=MultipleLocator(30),\n format=FuncFormatter(lambda x, pos: \"{:.0f}\u00b0\".format(x)),\n cax=cax,\n orientation=\"vertical\",\n)\naxs[1, 1].set(title=\"Angles (smooth=4)\")\n\nim = axs[0, 2].imshow(anisos_sm0, cmap=\"Reds\", vmin=0, vmax=1)\ncax = make_axes_locatable(axs[0, 2]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\naxs[0, 2].set(title=\"Anisotropy (smooth=0)\")\n\nim = axs[1, 2].imshow(anisos_sm4, cmap=\"Reds\", vmin=0, vmax=1)\ncax = make_axes_locatable(axs[1, 2]).append_axes(\"right\", size=\"5%\", pad=0.05)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\naxs[1, 2].set(title=\"Anisotropy (smooth=4)\")\n\nfor ax in axs.ravel():\n ax.axis(\"off\")\nfig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Core samples\nThe original paper by van Vliet and Verbeek [1] also has an example with images\nof core samples. Since the original paper does not have a scale with which to\nplot the angles, we have chosen ours it to match their image as closely as\npossible.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "core_input = rgb2gray(imread(\"../testdata/slope_estimate/core_sample.png\"))\ncore_angles = rgb2gray(imread(\"../testdata/slope_estimate/core_sample_orientation.png\"))\ncore_aniso = rgb2gray(imread(\"../testdata/slope_estimate/core_sample_anisotropy.png\"))\n\n\nangles, anisos_sm4 = dip_estimate(core_input, smooth=4)\nangles_sm4 = np.rad2deg(angles)\n\nangles, anisos_sm8 = dip_estimate(core_input, smooth=8)\nangles_sm8 = np.rad2deg(angles)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1, 6, figsize=(10, 6))\n\naxs[0].imshow(core_input, cmap=\"gray_r\", aspect=\"equal\")\naxs[0].set(title=\"Original\\nImage\")\ncax = make_axes_locatable(axs[0]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncax.axis(\"off\")\n\naxs[1].imshow(-core_angles, cmap=\"YlGnBu_r\")\naxs[1].set(title=\"Original\\nAngles\")\ncax = make_axes_locatable(axs[1]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncax.axis(\"off\")\n\nim = axs[2].imshow(angles_sm8, cmap=\"YlGnBu_r\", vmin=-49, vmax=-11)\ncax = make_axes_locatable(axs[2]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncb = fig.colorbar(\n im,\n ticks=MultipleLocator(30),\n format=FuncFormatter(lambda x, pos: \"{:.0f}\u00b0\".format(x)),\n cax=cax,\n orientation=\"vertical\",\n)\naxs[2].set(title=\"Angles\\n(smooth=8)\")\n\nim = axs[3].imshow(angles_sm4, cmap=\"YlGnBu_r\", vmin=-49, vmax=-11)\ncax = make_axes_locatable(axs[3]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncb = fig.colorbar(\n im,\n ticks=MultipleLocator(30),\n format=FuncFormatter(lambda x, pos: \"{:.0f}\u00b0\".format(x)),\n cax=cax,\n orientation=\"vertical\",\n)\naxs[3].set(title=\"Angles\\n(smooth=4)\")\n\nim = axs[4].imshow(anisos_sm8, cmap=\"Reds\", vmin=0, vmax=1)\ncax = make_axes_locatable(axs[4]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\naxs[4].set(title=\"Anisotropy\\n(smooth=8)\")\n\nim = axs[5].imshow(anisos_sm4, cmap=\"Reds\", vmin=0, vmax=1)\ncax = make_axes_locatable(axs[5]).append_axes(\"right\", size=\"20%\", pad=0.05)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\naxs[5].set(title=\"Anisotropy\\n(smooth=4)\")\n\nfor ax in axs.ravel():\n ax.axis(\"off\")\nfig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Final considerations\nAs you can see the Structure Tensor algorithm is a very fast, general purpose\nalgorithm that can be used to estimate local slopes to input datasets of\nvery different natures.\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.15" } }, "nbformat": 4, "nbformat_minor": 0 }