{
  "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 of a two-dimensional array\nusing :py:func:`pylops.utils.signalprocessing.slope_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\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\nfrom scipy.ndimage import gaussian_filter\n\nimport pylops\nfrom pylops.signalprocessing.Seislet import _predict_trace\n\nplt.close('all')\nnp.random.seed(10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To start we import a 2d image and estimate the local slopes 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.5\nim = gaussian_filter(im, sigma=2)\n\nslopes, linearity = \\\n    pylops.utils.signalprocessing.slope_estimate(im, 1., 1., smooth=7)\n\nfig, axs = plt.subplots(1, 3, figsize=(12, 4))\niax = axs[0].imshow(im, cmap='viridis', origin='lower')\naxs[0].set_title('Data')\naxs[0].axis('tight')\niax = axs[1].imshow(-np.rad2deg(slopes), cmap='seismic', origin='lower',\n                    vmin=-90, vmax=90)\naxs[1].set_title('Slopes (degrees)')\nplt.colorbar(iax, ax=axs[1])\naxs[1].axis('tight')\niax = axs[2].imshow(linearity, cmap='hot', origin='lower', vmin=0, vmax=1)\naxs[2].set_title('Linearity')\nplt.colorbar(iax, ax=axs[2])\naxs[2].axis('tight')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We 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 = 4, 0.004\nx, t = np.arange(nx)*dx, np.arange(nt)*dt\n\nrefl = np.zeros(nt)\nit = np.sort(np.random.permutation(np.arange(10, nt-20))[:nt//8])\nrefl[it] = np.random.normal(0., 1., nt//8)\n\n# Wavelet\nntwav = 41\nf0 = 30\nwav = pylops.utils.wavelets.ricker(np.arange(ntwav)*0.004, f0)[0]\nwavc = np.argmax(wav)\n\n# Input trace\ntrace = np.convolve(refl, wav, mode='same')\n\n# Slopes\ntheta = np.linspace(0, 30, nx)\nslope = np.outer(np.ones(nt), np.deg2rad(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 = -pylops.utils.signalprocessing.slope_estimate(d, dt, dx,\n                                                          smooth=10)[0]\n\nfig, axs = plt.subplots(1, 3, figsize=(12, 5))\naxs[0].imshow(d, cmap='gray', vmin=-1, vmax=1,\n              extent=(x[0], x[-1], t[-1], t[0]))\naxs[0].set_title('Data')\naxs[0].axis('tight')\naxs[1].imshow(np.rad2deg(slope)*dx/dt, cmap='seismic', vmin=0, vmax=40,\n              extent=(x[0], x[-1], t[-1], t[0]))\naxs[1].set_title('True Slopes')\naxs[1].axis('tight')\niax = axs[2].imshow(np.rad2deg(slope_est)*dx/dt, cmap='seismic',\n                    vmin=0, vmax=40,\n                    extent=(x[0], x[-1], t[-1], t[0]))\naxs[2].set_title('Estimated Slopes')\nplt.colorbar(iax, ax=axs[2])\naxs[2].axis('tight')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As 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 nature.\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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}