{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 11. Radon filtering\nIn this example we will be taking advantage of the\n:py:class:`pylops.signalprocessing.Radon2D` operator to perform filtering of\nunwanted events from a seismic data. For those of you not familiar with seismic\ndata, let's imagine that we have a data composed of a certain number of flat\nevents and a parabolic event , we are after a transform that allows us to\nseparate such an event from the others and filter it out.\nThose of you with a geophysics background may immediately realize this\nis the case of seismic angle (or offset) gathers after migration and those\nevents with parabolic moveout are generally residual multiples that we would\nlike to suppress prior to performing further analysis of our data.\n\nThe Radon transform is actually a very good transform to perform such a\nseparation. We can thus devise a simple workflow that takes our data as input,\napplies a Radon transform, filters some of the events out and goes back to the\noriginal domain.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close(\"all\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first create a data composed on 3 linear events and a parabolic event.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "par = {\"ox\": 0, \"dx\": 2, \"nx\": 121, \"ot\": 0, \"dt\": 0.004, \"nt\": 100, \"f0\": 30}\n\n# linear events\nv = 1500  # m/s\nt0 = [0.1, 0.2, 0.3]  # s\ntheta = [0, 0, 0]\namp = [1.0, -2, 0.5]\n\n# parabolic event\ntp0 = [0.13]  # s\npx = [0]  # s/m\npxx = [5e-7]  # s\u00b2/m\u00b2\nampp = [0.7]\n\n# create axis\ntaxis, taxis2, xaxis, yaxis = pylops.utils.seismicevents.makeaxis(par)\n\n# create wavelet\nwav = ricker(taxis[:41], f0=par[\"f0\"])[0]\n\n# generate model\ny = (\n    pylops.utils.seismicevents.linear2d(xaxis, taxis, v, t0, theta, amp, wav)[1]\n    + pylops.utils.seismicevents.parabolic2d(xaxis, taxis, tp0, px, pxx, ampp, wav)[1]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now create the :py:class:`pylops.signalprocessing.Radon2D` operator.\nWe also apply its adjoint to the data to obtain a representation of those\n3 linear events overlapping to a parabolic event in the Radon domain.\nSimilarly, we feed the operator to a sparse solver like\n:py:class:`pylops.optimization.sparsity.FISTA` to obtain a sparse\nrepresention of the data in the Radon domain. At this point we try to filter\nout the unwanted event. We can see how this is much easier for the sparse\ntransform as each event has a much more compact representation in the Radon\ndomain than for the adjoint transform.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# radon operator\nnpx = 61\npxmax = 5e-4  # s/m\npx = np.linspace(-pxmax, pxmax, npx)\n\nRop = pylops.signalprocessing.Radon2D(\n    taxis, xaxis, px, kind=\"linear\", interp=\"nearest\", centeredh=False, dtype=\"float64\"\n)\n\n# adjoint Radon transform\nxadj = Rop.H * y\n\n# sparse Radon transform\nxinv, niter, cost = pylops.optimization.sparsity.fista(\n    Rop, y.ravel(), niter=15, eps=1e1\n)\nxinv = xinv.reshape(Rop.dims)\n\n# filtering\nxfilt = np.zeros_like(xadj)\nxfilt[npx // 2 - 3 : npx // 2 + 4] = xadj[npx // 2 - 3 : npx // 2 + 4]\n\nyfilt = Rop * xfilt\n\n# filtering on sparse transform\nxinvfilt = np.zeros_like(xinv)\nxinvfilt[npx // 2 - 3 : npx // 2 + 4] = xinv[npx // 2 - 3 : npx // 2 + 4]\n\nyinvfilt = Rop * xinvfilt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we visualize our results.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pclip = 0.7\nfig, axs = plt.subplots(1, 5, sharey=True, figsize=(12, 5))\naxs[0].imshow(\n    y.T,\n    cmap=\"gray\",\n    vmin=-pclip * np.abs(y).max(),\n    vmax=pclip * np.abs(y).max(),\n    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),\n)\naxs[0].set(xlabel=\"$x$ [m]\", ylabel=\"$t$ [s]\", title=\"Data\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    xadj.T,\n    cmap=\"gray\",\n    vmin=-pclip * np.abs(xadj).max(),\n    vmax=pclip * np.abs(xadj).max(),\n    extent=(px[0], px[-1], taxis[-1], taxis[0]),\n)\naxs[1].axvline(px[npx // 2 - 3], color=\"r\", linestyle=\"--\")\naxs[1].axvline(px[npx // 2 + 3], color=\"r\", linestyle=\"--\")\naxs[1].set(xlabel=\"$p$ [s/m]\", title=\"Radon\")\naxs[1].axis(\"tight\")\naxs[2].imshow(\n    yfilt.T,\n    cmap=\"gray\",\n    vmin=-pclip * np.abs(yfilt).max(),\n    vmax=pclip * np.abs(yfilt).max(),\n    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),\n)\naxs[2].set(xlabel=\"$x$ [m]\", title=\"Filtered data\")\naxs[2].axis(\"tight\")\naxs[3].imshow(\n    xinv.T,\n    cmap=\"gray\",\n    vmin=-pclip * np.abs(xinv).max(),\n    vmax=pclip * np.abs(xinv).max(),\n    extent=(px[0], px[-1], taxis[-1], taxis[0]),\n)\naxs[3].axvline(px[npx // 2 - 3], color=\"r\", linestyle=\"--\")\naxs[3].axvline(px[npx // 2 + 3], color=\"r\", linestyle=\"--\")\naxs[3].set(xlabel=\"$p$ [s/m]\", title=\"Sparse Radon\")\naxs[3].axis(\"tight\")\naxs[4].imshow(\n    yinvfilt.T,\n    cmap=\"gray\",\n    vmin=-pclip * np.abs(y).max(),\n    vmax=pclip * np.abs(y).max(),\n    extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]),\n)\n\naxs[4].set(xlabel=\"$x$ [m]\", title=\"Sparse filtered data\")\naxs[4].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As expected, the Radon domain is a suitable domain for this type of filtering\nand the sparse transform improves the ability to filter out parabolic events\nwith small curvature.\n\nOn the other hand, it is important to note that we have not been able to\ncorrectly preserve the amplitudes of each event. This is because the sparse\nRadon transform can only identify a sparsest response that explain the data\nwithin a certain threshold. For this reason a more suitable approach for\npreserving amplitudes could be to apply a parabolic Raodn transform with the\naim of reconstructing only the unwanted event and apply an adaptive\nsubtraction between the input data and the reconstructed unwanted event.\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
}