{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Multi-Dimensional Convolution\nThis example shows how to use the :py:class:`pylops.waveeqprocessing.MDC` operator\nto convolve a 3D kernel with an input seismic data. The resulting data is\na blurred version of the input data and the problem of removing such blurring\nis reffered to as *Multi-dimensional Deconvolution (MDD)* and its implementation\nis discussed in more details in the **MDD** tutorial.\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\nfrom pylops.utils.tapers import taper3d\nfrom pylops.utils.wavelets import ricker\nfrom pylops.utils.seismicevents import makeaxis, hyperbolic2d\n\nplt.close('all')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating a set of hyperbolic events to be used as our MDC kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Input parameters\npar = {'ox':-300, 'dx':10, 'nx':61,\n       'oy':-500, 'dy':10, 'ny':101,\n       'ot':0, 'dt':0.004, 'nt':400,\n       'f0': 20, 'nfmax': 200}\n\nt0_m = 0.2\nvrms_m = 1100.0\namp_m = 1.0\n\nt0_G = (0.2, 0.5, 0.7)\nvrms_G = (1200., 1500., 2000.)\namp_G = (1., 0.6, 0.5)\n\n# Taper\ntap = taper3d(par['nt'], (par['ny'], par['nx']),\n              (5, 5), tapertype='hanning')\n\n# Create axis\nt, t2, x, y = makeaxis(par)\n\n# Create wavelet\nwav = ricker(t[:41], f0=par['f0'])[0]\n\n# Generate model\nm, mwav = hyperbolic2d(x, t, t0_m, vrms_m, amp_m, wav)\n\n# Generate operator\nG, Gwav = np.zeros((par['ny'], par['nx'], par['nt'])), \\\n          np.zeros((par['ny'], par['nx'], par['nt']))\nfor iy, y0 in enumerate(y):\n    G[iy], Gwav[iy] = hyperbolic2d(x-y0, t, t0_G, vrms_G, amp_G, wav)\nG, Gwav = G*tap, Gwav*tap\n\n# Add negative part to data and model\nm = np.concatenate((np.zeros((par['nx'], par['nt']-1)), m), axis=-1)\nmwav = np.concatenate((np.zeros((par['nx'], par['nt']-1)), mwav), axis=-1)\nGwav2 = np.concatenate((np.zeros((par['ny'], par['nx'], par['nt']-1)), Gwav), axis=-1)\n\n# Define MDC linear operator\nGwav_fft = np.fft.rfft(Gwav2, 2*par['nt']-1, axis=-1)\nGwav_fft = Gwav_fft[..., :par['nfmax']]\n\n# Move frequency/time to first axis\nm, mwav = m.T, mwav.T\nGwav_fft = Gwav_fft.transpose(2,0,1)\n\n# Create operator\nMDCop = pylops.waveeqprocessing.MDC(Gwav_fft, nt=2 * par['nt'] - 1, nv=1,\n                                    dt=0.004, dr=1., transpose=False,\n                                    dtype='float32')\n\n# Create data\nd = MDCop*m.flatten()\nd = d.reshape(2*par['nt']-1, par['ny'])\n\n# Apply adjoint operator to data\nmadj = MDCop.H*d.flatten()\nmadj = madj.reshape(2*par['nt']-1, par['nx'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally let's display the operator, input model, data and adjoint model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(1, 2, figsize=(9, 6))\naxs[0].imshow(Gwav2[int(par['ny']/2)].T, aspect='auto',\n              interpolation='nearest', cmap='gray',\n              vmin=-Gwav2.max(), vmax=Gwav2.max(),\n              extent=(x.min(), x.max(), t2.max(), t2.min()))\naxs[0].set_title('G - inline view', fontsize=15)\naxs[0].set_xlabel('r')\naxs[1].set_ylabel('t')\naxs[1].imshow(Gwav2[:, int(par['nx']/2)].T, aspect='auto',\n              interpolation='nearest', cmap='gray',\n              vmin=-Gwav2.max(), vmax=Gwav2.max(),\n              extent=(y.min(), y.max(), t2.max(), t2.min()))\naxs[1].set_title('G - inline view', fontsize=15)\naxs[1].set_xlabel('s')\naxs[1].set_ylabel('t')\nfig.tight_layout()\n\nfig, axs = plt.subplots(1, 3, figsize=(9, 6))\naxs[0].imshow(mwav, aspect='auto',\n              interpolation='nearest', cmap='gray',\n              vmin=-mwav.max(), vmax=mwav.max(),\n              extent=(x.min(), x.max(), t2.max(), t2.min()))\naxs[0].set_title(r'$m$', fontsize=15)\naxs[0].set_xlabel('r')\naxs[0].set_ylabel('t')\naxs[1].imshow(d, aspect='auto', interpolation='nearest', cmap='gray',\n              vmin=-d.max(), vmax=d.max(),\n              extent=(x.min(), x.max(), t2.max(), t2.min()))\naxs[1].set_title(r'$d$', fontsize=15)\naxs[1].set_xlabel('s')\naxs[1].set_ylabel('t')\naxs[2].imshow(madj, aspect='auto', interpolation='nearest', cmap='gray',\n              vmin=-madj.max(), vmax=madj.max(),\n              extent=(x.min(), x.max(), t2.max(), t2.min()))\naxs[2].set_title(r'$m_{adj}$', fontsize=15)\naxs[2].set_xlabel('s')\naxs[2].set_ylabel('t')\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
}