{
  "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\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\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(t, h, px, centeredh=True,\n                                       kind='linear', interp=False,\n                                       engine='numpy')\nRPop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,\n                                       kind='parabolic', interp=False,\n                                       engine='numpy')\nRHop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,\n                                       kind='hyperbolic', interp=False,\n                                       engine='numpy')\n\n# forward\nyL = RLop * x.ravel()\nyP = RPop * x.ravel()\nyH = RHop * x.ravel()\nyL = yL.reshape(nh, nt)\nyP = yP.reshape(nh, nt)\nyH = yH.reshape(nh, nt)\n\n# adjoint\nxadjL = RLop.H * yL.ravel()\nxadjP = RPop.H * yP.ravel()\nxadjH = RHop.H * yH.ravel()\nxadjL = xadjL.reshape(npx, nt)\nxadjP = xadjP.reshape(npx, nt)\nxadjH = xadjH.reshape(npx, nt)"
      ]
    },
    {
      "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))\naxs[0][0].imshow(x.T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[0][0].set_title('Input model')\naxs[0][0].axis('tight')\naxs[0][1].imshow(yL.T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(h[0], h[-1], t[-1], t[0]))\naxs[0][1].set_title('Linear data')\naxs[0][1].axis('tight')\naxs[0][2].imshow(yP.T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(h[0], h[-1], t[-1], t[0]))\naxs[0][2].set_title('Parabolic data')\naxs[0][2].axis('tight')\naxs[0][3].imshow(yH.T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(h[0], h[-1], t[-1], t[0]))\naxs[0][3].set_title('Hyperbolic data')\naxs[0][3].axis('tight')\naxs[1][1].imshow(xadjL.T, vmin=-20, vmax=20, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][0].axis('off')\naxs[1][1].set_title('Linear adjoint')\naxs[1][1].axis('tight')\naxs[1][2].imshow(xadjP.T, vmin=-20, vmax=20, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][2].set_title('Parabolic adjoint')\naxs[1][2].axis('tight')\naxs[1][3].imshow(xadjH.T, vmin=-20, vmax=20, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][3].set_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(t, hy, hx, py, px, centeredh=True,\n                                       kind='linear', interp=False,\n                                       engine='numpy')\nRPop = pylops.signalprocessing.Radon3D(t, hy, hx, py, px, centeredh=True,\n                                       kind='parabolic', interp=False,\n                                       engine='numpy')\nRHop = pylops.signalprocessing.Radon3D(t, hy, hx, py, px, centeredh=True,\n                                       kind='hyperbolic', interp=False,\n                                       engine='numpy')\n\n# forward\nyL = RLop * x.ravel()\nyP = RPop * x.ravel()\nyH = RHop * x.ravel()\nyL = yL.reshape(ny, nx, nt)\nyP = yP.reshape(ny, nx, nt)\nyH = yH.reshape(ny, nx, nt)\n\n# adjoint\nxadjL = RLop.H * yL.ravel()\nxadjP = RPop.H * yP.ravel()\nxadjH = RHop.H * yH.ravel()\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))\naxs[0][0].imshow(x[npy//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[0][0].set_title('Input model')\naxs[0][0].axis('tight')\naxs[0][1].imshow(yL[ny//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hx[0], hx[-1], t[-1], t[0]))\naxs[0][1].set_title('Linear data')\naxs[0][1].axis('tight')\naxs[0][2].imshow(yP[ny//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hx[0], hx[-1], t[-1], t[0]))\naxs[0][2].set_title('Parabolic data')\naxs[0][2].axis('tight')\naxs[0][3].imshow(yH[ny//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hx[0], hx[-1], t[-1], t[0]))\naxs[0][3].set_title('Hyperbolic data')\naxs[0][3].axis('tight')\naxs[1][1].imshow(xadjL[npy//2].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][0].axis('off')\naxs[1][1].set_title('Linear adjoint')\naxs[1][1].axis('tight')\naxs[1][2].imshow(xadjP[npy//2].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][2].set_title('Parabolic adjoint')\naxs[1][2].axis('tight')\naxs[1][3].imshow(xadjH[npy//2].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(px[0], px[-1], t[-1], t[0]))\naxs[1][3].set_title('Hyperbolic adjoint')\naxs[1][3].axis('tight')\nfig.tight_layout()\n\nfig, axs = plt.subplots(2, 4, figsize=(10, 6))\naxs[0][0].imshow(x[:, npx//2-2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(py[0], py[-1], t[-1], t[0]))\naxs[0][0].set_title('Input model')\naxs[0][0].axis('tight')\naxs[0][1].imshow(yL[:, nx//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hy[0], hy[-1], t[-1], t[0]))\naxs[0][1].set_title('Linear data')\naxs[0][1].axis('tight')\naxs[0][2].imshow(yP[:, nx//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hy[0], hy[-1], t[-1], t[0]))\naxs[0][2].set_title('Parabolic data')\naxs[0][2].axis('tight')\naxs[0][3].imshow(yH[:, nx//2].T, vmin=-1, vmax=1, cmap='seismic_r',\n                 extent=(hy[0], hy[-1], t[-1], t[0]))\naxs[0][3].set_title('Hyperbolic data')\naxs[0][3].axis('tight')\naxs[1][1].imshow(xadjL[:, npx//2-5].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(py[0], py[-1], t[-1], t[0]))\naxs[1][0].axis('off')\naxs[1][1].set_title('Linear adjoint')\naxs[1][1].axis('tight')\naxs[1][2].imshow(xadjP[:, npx//2-2].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(py[0], py[-1], t[-1], t[0]))\naxs[1][2].set_title('Parabolic adjoint')\naxs[1][2].axis('tight')\naxs[1][3].imshow(xadjH[:, npx//2-2].T, vmin=-100, vmax=100, cmap='seismic_r',\n                 extent=(py[0], py[-1], t[-1], t[0]))\naxs[1][3].set_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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}