{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 1D, 2D and 3D Sliding\nThis example shows how to use the\n:py:class:`pylops.signalprocessing.Sliding1D`,\n:py:class:`pylops.signalprocessing.Sliding2D`\nand  :py:class:`pylops.signalprocessing.Sliding3D` operators\nto perform repeated transforms over small strides of a 1-, 2- or 3-dimensional\narray.\n\nFor the 1-d case, the transform that we apply in this example is the\n:py:class:`pylops.signalprocessing.FFT`.\n\nFor the 2- and 3-d cases, the transform that we apply in this example is the\n:py:class:`pylops.signalprocessing.Radon2D`\n(and :py:class:`pylops.signalprocessing.Radon3D`) but this operator has been\ndesign to allow a variety of transforms as long as they operate with signals\nthat are 2 or 3-dimensional in nature, respectively.\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 a 1-dimensional array of size $n_t$ and create\na sliding operator to compute its transformed representation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nwins = 4\nnwin = 26\nnover = 3\nnop = 64\ndimd = nwin * nwins - 3 * nover\n\nt = np.arange(dimd) * 0.004\ndata = np.sin(2 * np.pi * 20 * t)\n\nOp = pylops.signalprocessing.FFT(nwin, nfft=nop, real=True)\n\nnwins, dim, mwin_inends, dwin_inends = pylops.signalprocessing.sliding1d_design(\n    dimd, nwin, nover, (nop + 2) // 2\n)\nSlid = pylops.signalprocessing.Sliding1D(\n    Op.H,\n    dim,\n    dimd,\n    nwin,\n    nover,\n    tapertype=None,\n)\n\nx = Slid.H * data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now create a similar operator but we also add a taper to the overlapping\nparts of the patches and use it to reconstruct the original signal.\nThis is done by simply using the adjoint of the\n:py:class:`pylops.signalprocessing.Sliding1D` operator. Note that for non-\northogonal operators, this must be replaced by an inverse.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Slid = pylops.signalprocessing.Sliding1D(\n    Op.H, dim, dimd, nwin, nover, tapertype=\"cosine\"\n)\n\nreconstructed_data = Slid * x\n\nfig, axs = plt.subplots(1, 2, figsize=(15, 3))\naxs[0].plot(t, data, \"k\", label=\"Data\")\naxs[0].plot(t, reconstructed_data.real, \"--r\", label=\"Rec Data\")\naxs[0].legend()\naxs[1].set(xlabel=r\"$t$ [s]\", title=\"Original domain\")\nfor i in range(nwins):\n    axs[1].plot(Op.f, np.abs(x[i, :]), label=f\"Window {i+1}/{nwins}\")\naxs[1].set(xlabel=r\"$f$ [Hz]\", title=\"Transformed domain\")\naxs[1].legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now create a 2-dimensional array of size $n_x \\times n_t$\ncomposed of 3 parabolic events\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "par = {\"ox\": -140, \"dx\": 2, \"nx\": 140, \"ot\": 0, \"dt\": 0.004, \"nt\": 200, \"f0\": 20}\n\nv = 1500\nt0 = [0.2, 0.4, 0.5]\npx = [0, 0, 0]\npxx = [1e-5, 5e-6, 1e-20]\namp = [1.0, -2, 0.5]\n\n# Create axis\nt, t2, x, y = pylops.utils.seismicevents.makeaxis(par)\n\n# Create wavelet\nwav = pylops.utils.wavelets.ricker(t[:41], f0=par[\"f0\"])[0]\n\n# Generate model\n_, data = pylops.utils.seismicevents.parabolic2d(x, t, t0, px, pxx, amp, wav)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We want to divide this 2-dimensional data into small overlapping\npatches in the spatial direction and apply the adjoint of the\n:py:class:`pylops.signalprocessing.Radon2D` operator to each patch. This is\ndone by simply using the adjoint of the\n:py:class:`pylops.signalprocessing.Sliding2D` operator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "winsize = 36\noverlap = 10\nnpx = 61\npx = np.linspace(-5e-3, 5e-3, npx)\ndimsd = data.shape\n\n# Sliding window transform without taper\nOp = pylops.signalprocessing.Radon2D(\n    t,\n    np.linspace(-par[\"dx\"] * winsize // 2, par[\"dx\"] * winsize // 2, winsize),\n    px,\n    centeredh=True,\n    kind=\"linear\",\n    engine=\"numba\",\n)\n\nnwins, dims, mwin_inends, dwin_inends = pylops.signalprocessing.sliding2d_design(\n    dimsd, winsize, overlap, (npx, par[\"nt\"])\n)\nSlid = pylops.signalprocessing.Sliding2D(\n    Op, dims, dimsd, winsize, overlap, tapertype=None\n)\n\nradon = Slid.H * data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now create a similar operator but we also add a taper to the overlapping\nparts of the patches.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Slid = pylops.signalprocessing.Sliding2D(\n    Op, dims, dimsd, winsize, overlap, tapertype=\"cosine\"\n)\nreconstructed_data = Slid * radon\n\n# Reshape for plotting\nradon = radon.reshape(dims)\nreconstructed_data = reconstructed_data.reshape(dimsd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will see that our reconstructed signal presents some small artifacts.\nThis is because we have not inverted our operator but simply applied\nthe adjoint to estimate the representation of the input data in the Radon\ndomain. We can do better if we use the inverse instead.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "radoninv = pylops.LinearOperator(Slid, explicit=False).div(data.ravel(), niter=10)\nreconstructed_datainv = Slid * radoninv.ravel()\n\nradoninv = radoninv.reshape(dims)\nreconstructed_datainv = reconstructed_datainv.reshape(dimsd)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's finally visualize all the intermediate results as well as our final\ndata reconstruction after inverting the\n:py:class:`pylops.signalprocessing.Sliding2D` operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 10))\nim = axs[0][0].imshow(data.T, cmap=\"gray\")\naxs[0][0].set_title(\"Original data\")\nplt.colorbar(im, ax=axs[0][0])\naxs[0][0].axis(\"tight\")\nim = axs[0][1].imshow(radon.T, cmap=\"gray\")\naxs[0][1].set_title(\"Adjoint Radon\")\nplt.colorbar(im, ax=axs[0][1])\naxs[0][1].axis(\"tight\")\nim = axs[0][2].imshow(reconstructed_data.T, cmap=\"gray\")\naxs[0][2].set_title(\"Reconstruction from adjoint\")\nplt.colorbar(im, ax=axs[0][2])\naxs[0][2].axis(\"tight\")\naxs[1][0].axis(\"off\")\nim = axs[1][1].imshow(radoninv.T, cmap=\"gray\")\naxs[1][1].set_title(\"Inverse Radon\")\nplt.colorbar(im, ax=axs[1][1])\naxs[1][1].axis(\"tight\")\nim = axs[1][2].imshow(reconstructed_datainv.T, cmap=\"gray\")\naxs[1][2].set_title(\"Reconstruction from inverse\")\nplt.colorbar(im, ax=axs[1][2])\naxs[1][2].axis(\"tight\")\n\nfor i in range(0, 114, 24):\n    axs[0][0].axvline(i, color=\"w\", lw=1, ls=\"--\")\n    axs[0][0].axvline(i + winsize, color=\"k\", lw=1, ls=\"--\")\n    axs[0][0].text(\n        i + winsize // 2,\n        par[\"nt\"] - 10,\n        \"w\" + str(i // 24),\n        ha=\"center\",\n        va=\"center\",\n        weight=\"bold\",\n        color=\"w\",\n    )\n\nfor i in range(0, 305, 61):\n    axs[0][1].axvline(i, color=\"w\", lw=1, ls=\"--\")\n    axs[0][1].text(\n        i + npx // 2,\n        par[\"nt\"] - 10,\n        \"w\" + str(i // 61),\n        ha=\"center\",\n        va=\"center\",\n        weight=\"bold\",\n        color=\"w\",\n    )\n    axs[1][1].axvline(i, color=\"w\", lw=1, ls=\"--\")\n    axs[1][1].text(\n        i + npx // 2,\n        par[\"nt\"] - 10,\n        \"w\" + str(i // 61),\n        ha=\"center\",\n        va=\"center\",\n        weight=\"bold\",\n        color=\"w\",\n    )"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We notice two things, i)provided small enough patches and a transform\nthat can explain data *locally*, we have been able reconstruct our\noriginal data almost to perfection. ii) inverse is betten than adjoint as\nexpected as the adjoin does not only introduce small artifacts but also does\nnot respect the original amplitudes of the data.\n\nAn appropriate transform alongside with a sliding window approach will\nresult a very good approach for interpolation (or *regularization*) or\nirregularly sampled seismic data.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we do the same for a 3-dimensional array of size\n$n_y \\times n_x \\times n_t$ composed of 3 hyperbolic events\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "par = {\n    \"oy\": -13,\n    \"dy\": 2,\n    \"ny\": 14,\n    \"ox\": -17,\n    \"dx\": 2,\n    \"nx\": 18,\n    \"ot\": 0,\n    \"dt\": 0.004,\n    \"nt\": 50,\n    \"f0\": 30,\n}\n\nvrms = [200, 200]\nt0 = [0.05, 0.1]\namp = [1.0, -2]\n\n# Create axis\nt, t2, x, y = pylops.utils.seismicevents.makeaxis(par)\n\n# Create wavelet\nwav = pylops.utils.wavelets.ricker(t[:41], f0=par[\"f0\"])[0]\n\n# Generate model\n_, data = pylops.utils.seismicevents.hyperbolic3d(x, y, t, t0, vrms, vrms, amp, wav)\n\n# Sliding window plan\nwinsize = (5, 6)\noverlap = (2, 3)\nnpx = 21\npx = np.linspace(-5e-3, 5e-3, npx)\ndimsd = data.shape\n\n# Sliding window transform without taper\nOp = pylops.signalprocessing.Radon3D(\n    t,\n    np.linspace(-par[\"dy\"] * winsize[0] // 2, par[\"dy\"] * winsize[0] // 2, winsize[0]),\n    np.linspace(-par[\"dx\"] * winsize[1] // 2, par[\"dx\"] * winsize[1] // 2, winsize[1]),\n    px,\n    px,\n    centeredh=True,\n    kind=\"linear\",\n    engine=\"numba\",\n)\n\nnwins, dims, mwin_inends, dwin_inends = pylops.signalprocessing.sliding3d_design(\n    dimsd, winsize, overlap, (npx, npx, par[\"nt\"])\n)\nSlid = pylops.signalprocessing.Sliding3D(\n    Op, dims, dimsd, winsize, overlap, (npx, npx), tapertype=None\n)\n\nradon = Slid.H * data\n\nSlid = pylops.signalprocessing.Sliding3D(\n    Op, dims, dimsd, winsize, overlap, (npx, npx), tapertype=\"cosine\"\n)\n\nreconstructed_data = Slid * radon\n\nradoninv = pylops.LinearOperator(Slid, explicit=False).div(data.ravel(), niter=10)\nradoninv = radoninv.reshape(Slid.dims)\nreconstructed_datainv = Slid * radoninv\n\nfig, axs = plt.subplots(2, 3, sharey=True, figsize=(12, 7))\nim = axs[0][0].imshow(data[par[\"ny\"] // 2].T, cmap=\"gray\", vmin=-2, vmax=2)\naxs[0][0].set_title(\"Original data\")\nplt.colorbar(im, ax=axs[0][0])\naxs[0][0].axis(\"tight\")\nim = axs[0][1].imshow(\n    radon[nwins[0] // 2, :, :, npx // 2].reshape(nwins[1] * npx, par[\"nt\"]).T,\n    cmap=\"gray\",\n    vmin=-25,\n    vmax=25,\n)\naxs[0][1].set_title(\"Adjoint Radon\")\nplt.colorbar(im, ax=axs[0][1])\naxs[0][1].axis(\"tight\")\nim = axs[0][2].imshow(\n    reconstructed_data[par[\"ny\"] // 2].T, cmap=\"gray\", vmin=-1000, vmax=1000\n)\naxs[0][2].set_title(\"Reconstruction from adjoint\")\nplt.colorbar(im, ax=axs[0][2])\naxs[0][2].axis(\"tight\")\naxs[1][0].axis(\"off\")\nim = axs[1][1].imshow(\n    radoninv[nwins[0] // 2, :, :, npx // 2].reshape(nwins[1] * npx, par[\"nt\"]).T,\n    cmap=\"gray\",\n    vmin=-0.025,\n    vmax=0.025,\n)\naxs[1][1].set_title(\"Inverse Radon\")\nplt.colorbar(im, ax=axs[1][1])\naxs[1][1].axis(\"tight\")\nim = axs[1][2].imshow(\n    reconstructed_datainv[par[\"ny\"] // 2].T, cmap=\"gray\", vmin=-2, vmax=2\n)\naxs[1][2].set_title(\"Reconstruction from inverse\")\nplt.colorbar(im, ax=axs[1][2])\naxs[1][2].axis(\"tight\")\n\n\nfig, axs = plt.subplots(2, 3, figsize=(12, 7))\nim = axs[0][0].imshow(data[:, :, 25], cmap=\"gray\", vmin=-2, vmax=2)\naxs[0][0].set_title(\"Original data\")\nplt.colorbar(im, ax=axs[0][0])\naxs[0][0].axis(\"tight\")\nim = axs[0][1].imshow(\n    radon[nwins[0] // 2, :, :, :, 25].reshape(nwins[1] * npx, npx).T,\n    cmap=\"gray\",\n    vmin=-25,\n    vmax=25,\n)\naxs[0][1].set_title(\"Adjoint Radon\")\nplt.colorbar(im, ax=axs[0][1])\naxs[0][1].axis(\"tight\")\nim = axs[0][2].imshow(reconstructed_data[:, :, 25], cmap=\"gray\", vmin=-1000, vmax=1000)\naxs[0][2].set_title(\"Reconstruction from adjoint\")\nplt.colorbar(im, ax=axs[0][2])\naxs[0][2].axis(\"tight\")\naxs[1][0].axis(\"off\")\nim = axs[1][1].imshow(\n    radoninv[nwins[0] // 2, :, :, :, 25].reshape(nwins[1] * npx, npx).T,\n    cmap=\"gray\",\n    vmin=-0.025,\n    vmax=0.025,\n)\naxs[1][1].set_title(\"Inverse Radon\")\nplt.colorbar(im, ax=axs[1][1])\naxs[1][1].axis(\"tight\")\nim = axs[1][2].imshow(reconstructed_datainv[:, :, 25], cmap=\"gray\", vmin=-2, vmax=2)\naxs[1][2].set_title(\"Reconstruction from inverse\")\nplt.colorbar(im, ax=axs[1][2])\naxs[1][2].axis(\"tight\")\nplt.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
}