{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Radon Transform\nThis example shows how to use the :py:class:`pylops.signalprocessing.Radon2D`\nand :py:class:`pylops.signalprocessing.Radon3D` operators to apply the Radon\nTransform to 2-dimensional or 3-dimensional signals, respectively.\nIn our implementation both linear, parabolic and hyperbolic parametrization\ncan be chosen.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating an empty 2d matrix of size $n_{p_x} \\times n_t$\nand add a single spike in it. We will see that applying the forward\nRadon operator will result in a single event (linear, parabolic or\nhyperbolic) in the resulting data vector.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt, nh = 41, 51\nnpx, pxmax = 41, 1e-2\n\ndt, dh = 0.005, 1\nt = np.arange(nt) * dt\nh = np.arange(nh) * dh\npx = np.linspace(0, pxmax, npx)\n\nx = np.zeros((npx, nt))\nx[4, nt // 2] = 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now define our operators for different parametric curves and apply\nthem to the input model vector. We also apply the adjoint to the resulting\ndata vector.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "RLop = pylops.signalprocessing.Radon2D(\n    t, h, px, centeredh=True, kind=\"linear\", interp=False, engine=\"numpy\"\n)\nRPop = pylops.signalprocessing.Radon2D(\n    t, h, px, centeredh=True, kind=\"parabolic\", interp=False, engine=\"numpy\"\n)\nRHop = pylops.signalprocessing.Radon2D(\n    t, h, px, centeredh=True, kind=\"hyperbolic\", interp=False, engine=\"numpy\"\n)\n\n# forward\nyL = RLop * x\nyP = RPop * x\nyH = RHop * x\n\n# adjoint\nxadjL = RLop.H * yL\nxadjP = RPop.H * yP\nxadjH = RHop.H * yH"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's now visualize the input model in the Radon domain, the data, and\nthe adjoint model the different parametric curves.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(2, 4, figsize=(10, 6), sharey=True)\naxs[0, 0].axis(\"off\")\naxs[1, 0].imshow(\n    x.T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 0].set(xlabel=r\"$p$ [s/km]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[1, 0].axis(\"tight\")\naxs[0, 1].imshow(\n    yL.T, vmin=-1, vmax=1, cmap=\"seismic_r\", extent=(h[0], h[-1], t[-1], t[0])\n)\naxs[0, 1].tick_params(labelleft=True)\naxs[0, 1].set(xlabel=r\"$x$ [m]\", ylabel=r\"$t$ [s]\", title=\"Linear data\")\naxs[0, 1].axis(\"tight\")\naxs[0, 2].imshow(\n    yP.T, vmin=-1, vmax=1, cmap=\"seismic_r\", extent=(h[0], h[-1], t[-1], t[0])\n)\naxs[0, 2].set(xlabel=r\"$x$ [m]\", title=\"Parabolic data\")\naxs[0, 2].axis(\"tight\")\naxs[0, 3].imshow(\n    yH.T, vmin=-1, vmax=1, cmap=\"seismic_r\", extent=(h[0], h[-1], t[-1], t[0])\n)\naxs[0, 3].set(xlabel=r\"$x$ [m]\", title=\"Hyperbolic data\")\naxs[0, 3].axis(\"tight\")\naxs[1, 1].imshow(\n    xadjL.T,\n    vmin=-20,\n    vmax=20,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 1].set(xlabel=r\"$p$ [s/km]\", title=\"Linear adjoint\")\naxs[1, 1].axis(\"tight\")\naxs[1, 2].imshow(\n    xadjP.T,\n    vmin=-20,\n    vmax=20,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 2].set(xlabel=r\"$p$ [s/km]\", title=\"Parabolic adjoint\")\naxs[1, 2].axis(\"tight\")\naxs[1, 3].imshow(\n    xadjH.T,\n    vmin=-20,\n    vmax=20,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 3].set(xlabel=r\"$p$ [s/km]\", title=\"Hyperbolic adjoint\")\naxs[1, 3].axis(\"tight\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As we can see in the bottom figures, the adjoint Radon transform is far\nfrom being close to the inverse Radon transform, i.e.\n$\\mathbf{R^H}\\mathbf{R} \\neq \\mathbf{I}$ (compared to the case of FFT\nwhere the adjoint and inverse are equivalent, i.e.\n$\\mathbf{F^H}\\mathbf{F} = \\mathbf{I}$). In fact when we apply the\nadjoint Radon Transform we obtain a *model* that\nis a smoothed version of the original model polluted by smearing and\nartifacts. In tutorial `sphx_glr_tutorials_radonfiltering.py` we will\nexploit a sparsity-promiting Radon transform to perform filtering of unwanted\nsignals from an input data.\n\nFinally we repeat the same exercise with 3d data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt, ny, nx = 21, 21, 11\nnpy, pymax = 13, 5e-3\nnpx, pxmax = 11, 5e-3\n\ndt, dy, dx = 0.005, 1, 1\nt = np.arange(nt) * dt\nhy = np.arange(ny) * dy\nhx = np.arange(nx) * dx\n\npy = np.linspace(0, pymax, npy)\npx = np.linspace(0, pxmax, npx)\n\nx = np.zeros((npy, npx, nt))\nx[npy // 2, npx // 2 - 2, nt // 2] = 1\n\nRLop = pylops.signalprocessing.Radon3D(\n    t, hy, hx, py, px, centeredh=True, kind=\"linear\", interp=False, engine=\"numpy\"\n)\nRPop = pylops.signalprocessing.Radon3D(\n    t, hy, hx, py, px, centeredh=True, kind=\"parabolic\", interp=False, engine=\"numpy\"\n)\nRHop = pylops.signalprocessing.Radon3D(\n    t, hy, hx, py, px, centeredh=True, kind=\"hyperbolic\", interp=False, engine=\"numpy\"\n)\n\n# forward\nyL = RLop * x.reshape(npy * npx, nt)\nyP = RPop * x.reshape(npy * npx, nt)\nyH = RHop * x.reshape(npy * npx, nt)\n\n# adjoint\nxadjL = RLop.H * yL\nxadjP = RPop.H * yP\nxadjH = RHop.H * yH\n\n# reshape\nyL = yL.reshape(ny, nx, nt)\nyP = yP.reshape(ny, nx, nt)\nyH = yH.reshape(ny, nx, nt)\nxadjL = xadjL.reshape(npy, npx, nt)\nxadjP = xadjP.reshape(npy, npx, nt)\nxadjH = xadjH.reshape(npy, npx, nt)\n\n# plotting\nfig, axs = plt.subplots(2, 4, figsize=(10, 6), sharey=True)\naxs[1, 0].imshow(\n    x[npy // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 0].set(xlabel=r\"$p_x$ [s/km]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[1, 0].axis(\"tight\")\naxs[0, 1].imshow(\n    yL[ny // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hx[0], hx[-1], t[-1], t[0]),\n)\naxs[0, 1].tick_params(labelleft=True)\naxs[0, 1].set(xlabel=r\"$x$ [m]\", ylabel=r\"$t$ [s]\", title=\"Linear data\")\naxs[0, 1].axis(\"tight\")\naxs[0, 2].imshow(\n    yP[ny // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hx[0], hx[-1], t[-1], t[0]),\n)\naxs[0, 2].set(xlabel=r\"$x$ [m]\", title=\"Parabolic data\")\naxs[0, 2].axis(\"tight\")\naxs[0, 3].imshow(\n    yH[ny // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hx[0], hx[-1], t[-1], t[0]),\n)\naxs[0, 3].set(xlabel=r\"$x$ [m]\", title=\"Hyperbolic data\")\naxs[0, 3].axis(\"tight\")\naxs[1, 1].imshow(\n    xadjL[npy // 2].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[0, 0].axis(\"off\")\naxs[1, 1].set(xlabel=r\"$p_x$ [s/km]\", title=\"Linear adjoint\")\naxs[1, 1].axis(\"tight\")\naxs[1, 2].imshow(\n    xadjP[npy // 2].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 2].set(xlabel=r\"$p_x$ [s/km]\", title=\"Parabolic adjoint\")\naxs[1, 2].axis(\"tight\")\naxs[1, 3].imshow(\n    xadjH[npy // 2].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),\n)\naxs[1, 3].set(xlabel=r\"$p_x$ [s/km]\", title=\"Hyperbolic adjoint\")\naxs[1, 3].axis(\"tight\")\nfig.tight_layout()\n\nfig, axs = plt.subplots(2, 4, figsize=(10, 6), sharey=True)\naxs[1, 0].imshow(\n    x[:, npx // 2 - 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),\n)\naxs[1, 0].set(xlabel=r\"$p_y$ [s/km]\", ylabel=r\"$t$ [s]\", title=\"Input model\")\naxs[1, 0].axis(\"tight\")\naxs[0, 1].imshow(\n    yL[:, nx // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hy[0], hy[-1], t[-1], t[0]),\n)\naxs[0, 1].tick_params(labelleft=True)\naxs[0, 1].set(xlabel=r\"$y$ [m]\", ylabel=r\"$t$ [s]\", title=\"Linear data\")\naxs[0, 1].axis(\"tight\")\naxs[0, 2].imshow(\n    yP[:, nx // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hy[0], hy[-1], t[-1], t[0]),\n)\naxs[0, 2].set(xlabel=r\"$y$ [m]\", title=\"Parabolic data\")\naxs[0, 2].axis(\"tight\")\naxs[0, 3].imshow(\n    yH[:, nx // 2].T,\n    vmin=-1,\n    vmax=1,\n    cmap=\"seismic_r\",\n    extent=(hy[0], hy[-1], t[-1], t[0]),\n)\naxs[0, 3].set(xlabel=r\"$y$ [m]\", title=\"Hyperbolic data\")\naxs[0, 3].axis(\"tight\")\naxs[1, 1].imshow(\n    xadjL[:, npx // 2 - 5].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),\n)\naxs[0, 0].axis(\"off\")\naxs[1, 1].set(xlabel=r\"$p_y$ [s/km]\", title=\"Linear adjoint\")\naxs[1, 1].axis(\"tight\")\naxs[1, 2].imshow(\n    xadjP[:, npx // 2 - 2].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),\n)\naxs[1, 2].set(xlabel=r\"$p_y$ [s/km]\", title=\"Parabolic adjoint\")\naxs[1, 2].axis(\"tight\")\naxs[1, 3].imshow(\n    xadjH[:, npx // 2 - 2].T,\n    vmin=-100,\n    vmax=100,\n    cmap=\"seismic_r\",\n    extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),\n)\naxs[1, 3].set(xlabel=r\"$p_y$ [s/km]\", title=\"Hyperbolic adjoint\")\naxs[1, 3].axis(\"tight\")\nfig.tight_layout()"
      ]
    }
  ],
  "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
}