feat: added FourierRadon2d operator

CHANGELOG.md

@@ -2,7 +2,7 @@ Changelog
 =========
 
 # 2.3.1
-* Fixed bug in :py:mod:`pylops.utils.backend` (see [Issue #606](https://github.com/PyLops/pylops/issues/606))
+* Fixed bug in `pylops.utils.backend` (see [Issue #606](https://github.com/PyLops/pylops/issues/606))
 
 # 2.3.0
 

docs/source/api/index.rst

@@ -109,6 +109,7 @@ Signal processing
     Seislet
     Radon2D
     Radon3D
+    FourierRadon2D
     ChirpRadon2D
     ChirpRadon3D
     Sliding1D

examples/plot_fourierradon.py

@@ -0,0 +1,130 @@
+r"""
+Fourier Radon Transform
+=======================
+This example shows how to use the :py:class:`pylops.signalprocessing.FourierRadon2D`
+operator to apply the linear and parabolic Radon Transform to 2-dimensional signals.
+
+This operator provides a transformation equivalent to that of
+:py:class:`pylops.signalprocessing.Radon2D`, however since the shift-and-sum step
+is performed in the frequency domain, this is analytically correct (compared to
+performing to shifting the data via nearest or linear interpolation).
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pylops
+
+plt.close("all")
+
+###############################################################################
+# Let's start by creating a empty 2d matrix of size :math:`n_x \times n_t`
+# with a single linear event.
+
+par = {
+    "ot": 0,
+    "dt": 0.004,
+    "nt": 51,
+    "ox": -250,
+    "dx": 10,
+    "nx": 51,
+    "oy": -250,
+    "dy": 10,
+    "ny": 51,
+    "f0": 40,
+}
+theta = [10.0]
+t0 = [0.1]
+amp = [1.0]
+
+# Create axes
+t, t2, x, y = pylops.utils.seismicevents.makeaxis(par)
+dt, dx, dy = par["dt"], par["dx"], par["dy"]
+
+# Create wavelet
+wav, _, wav_c = pylops.utils.wavelets.ricker(t[:41], f0=par["f0"])
+
+# Generate data
+_, d = pylops.utils.seismicevents.linear2d(x, t, 1500.0, t0, theta, amp, wav)
+
+
+###############################################################################
+# We can now define our operators and apply the forward, adjoint and inverse
+# steps.
+nfft = int(2 ** np.ceil(np.log2(par["nt"])))
+npx, pxmax = 2 * par["nx"], 5e-4
+px = np.linspace(-pxmax, pxmax, npx)
+
+R2Op = pylops.signalprocessing.FourierRadon2D(
+    t, x, px, nfft, kind="linear", engine="numpy", dtype="float64"
+)
+dL_chirp = R2Op.H * d
+dadj_chirp = R2Op * dL_chirp
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+axs[0].imshow(d.T, vmin=-1, vmax=1, cmap="bwr_r", extent=(x[0], x[-1], t[-1], t[0]))
+axs[0].set(xlabel=r"$x$ [m]", ylabel=r"$t$ [s]", title="Input linear")
+axs[0].axis("tight")
+axs[1].imshow(
+    dL_chirp.T,
+    cmap="bwr_r",
+    vmin=-dL_chirp.max(),
+    vmax=dL_chirp.max(),
+    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),
+)
+axs[1].scatter(1e3 * np.sin(np.deg2rad(theta[0])) / 1500.0, t0[0], s=50, color="r")
+axs[1].set(xlabel=r"$p$ [s/km]", title="Radon")
+axs[1].axis("tight")
+axs[2].imshow(
+    dadj_chirp.T,
+    cmap="bwr_r",
+    vmin=-dadj_chirp.max(),
+    vmax=dadj_chirp.max(),
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon")
+axs[2].axis("tight")
+plt.tight_layout()
+
+
+###############################################################################
+# We repeat now the same with a parabolic event
+
+# Generate data
+pxx = [1e-6]
+_, d = pylops.utils.seismicevents.parabolic2d(x, t, t0, 0, np.array(pxx), amp, wav)
+
+# Radon transform
+npx, pxmax = 2 * par["nx"], 5e-6
+px = np.linspace(-pxmax, pxmax, npx)
+
+R2Op = pylops.signalprocessing.FourierRadon2D(
+    t, x, px, nfft, kind="parabolic", engine="numpy", dtype="float64"
+)
+dL_chirp = R2Op.H * d
+dadj_chirp = R2Op * dL_chirp
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+axs[0].imshow(d.T, vmin=-1, vmax=1, cmap="bwr_r", extent=(x[0], x[-1], t[-1], t[0]))
+axs[0].set(xlabel=r"$x$ [m]", ylabel=r"$t$ [s]", title="Input parabolic")
+axs[0].axis("tight")
+axs[1].imshow(
+    dL_chirp.T,
+    cmap="bwr_r",
+    vmin=-dL_chirp.max(),
+    vmax=dL_chirp.max(),
+    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),
+)
+axs[1].scatter(1e3 * pxx[0], t0[0], s=50, color="r")
+axs[1].set(xlabel=r"$p$ [s/km]", title="Radon")
+axs[1].axis("tight")
+axs[2].imshow(
+    dadj_chirp.T,
+    cmap="bwr_r",
+    vmin=-dadj_chirp.max(),
+    vmax=dadj_chirp.max(),
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon")
+axs[2].axis("tight")
+plt.tight_layout()

pylops/signalprocessing/__init__.py

@@ -28,6 +28,9 @@
     DTCWT                           Dual-Tree Complex Wavelet Transform.
     Radon2D	                        Two dimensional Radon transform.
     Radon3D	                        Three dimensional Radon transform.
+    FourierRadon2D	                Two dimensional Fourier Radon transform.
+    ChirpRadon2D	                Two dimensional Chirp Radon transform.
+    ChirpRadon3D	                Three dimensional Chirp Radon transform.
     Seislet                         Two dimensional Seislet operator.
     Sliding1D                       1D Sliding transform operator.
     Sliding2D                       2D Sliding transform operator.
@@ -52,6 +55,7 @@
 from .bilinear import *
 from .radon2d import *
 from .radon3d import *
+from .fourierradon2d import *
 from .chirpradon2d import *
 from .chirpradon3d import *
 from .sliding1d import *

pylops/signalprocessing/_fourierradon2d_cuda.py

@@ -0,0 +1,86 @@
+from cmath import exp
+from math import pi
+
+from numba import cuda
+
+
+@cuda.jit()
+def _radon_inner_2d_kernel(x, y, f, px, h, flim0, flim1, npx, nh):
+    """Cuda kernels for FourierRadon2D operator
+
+    Cuda implementation of the on-the-fly kernel creation and application for the
+    FourierRadon2D operator. See :class:`pylops.signalprocessing.FourierRadon2D`
+    for details about input parameters.
+
+    """
+    ih, ifr = cuda.grid(2)
+    if ih < nh and ifr >= flim0 and ifr <= flim1:
+        for ipx in range(npx):
+            y[ih, ifr] += x[ipx, ifr] * exp(-1j * 2 * pi * f[ifr] * px[ipx] * h[ih])
+
+
+@cuda.jit()
+def _aradon_inner_2d_kernel(x, y, f, px, h, flim0, flim1, npx, nh):
+    """Cuda kernels for FourierRadon2D operator
+
+    Cuda implementation of the on-the-fly kernel creation and application for the
+    FourierRadon2D operator. See :class:`pylops.signalprocessing.FourierRadon2D`
+    for details about input parameters.
+
+    """
+    ipx, ifr = cuda.grid(2)
+    if ipx < npx and ifr >= flim0 and ifr <= flim1:
+        for ih in range(nh):
+            x[ipx, ifr] += y[ih, ifr] * exp(1j * 2 * pi * f[ifr] * px[ipx] * h[ih])
+
+
+def _radon_inner_2d_cuda(
+    x,
+    y,
+    f,
+    px,
+    h,
+    flim0,
+    flim1,
+    npx,
+    nh,
+    num_blocks=(32, 32),
+    num_threads_per_blocks=(32, 32),
+):
+    """Caller for FourierRadon2D operator
+
+    Caller for cuda implementation of matvec kernel for FourierRadon2D operator.
+    See :class:`pylops.signalprocessing.FourierRadon2D` for details about
+    input parameters.
+
+    """
+    _radon_inner_2d_kernel[num_blocks, num_threads_per_blocks](
+        x, y, f, px, h, flim0, flim1, npx, nh
+    )
+    return y
+
+
+def _aradon_inner_2d_cuda(
+    x,
+    y,
+    f,
+    px,
+    h,
+    flim0,
+    flim1,
+    npx,
+    nh,
+    num_blocks=(32, 32),
+    num_threads_per_blocks=(32, 32),
+):
+    """Caller for FourierRadon2D operator
+
+    Caller for cuda implementation of rmatvec kernel for FourierRadon2D operator.
+    See :class:`pylops.signalprocessing.FourierRadon2D` for details about
+    input parameters.
+
+    """
+    _aradon_inner_2d_kernel[num_blocks, num_threads_per_blocks](
+        x, y, f, px, h, flim0, flim1, npx, nh
+    )
+    return x

pylops/signalprocessing/_fourierradon2d_numba.py

@@ -0,0 +1,30 @@
+import os
+
+import numpy as np
+from numba import jit, prange
+
+# detect whether to use parallel or not
+numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1"))
+parallel = True if numba_threads != 1 else False
+
+
+@jit(nopython=True, parallel=parallel, nogil=True, cache=True, fastmath=True)
+def _radon_inner_2d(X, Y, f, px, h, flim0, flim1, npx, nh):
+    for ih in prange(nh):
+        for ifr in range(flim0, flim1):
+            for ipx in range(npx):
+                Y[ih, ifr] += X[ipx, ifr] * np.exp(
+                    -1j * 2 * np.pi * f[ifr] * px[ipx] * h[ih]
+                )
+    return Y
+
+
+@jit(nopython=True, parallel=parallel, nogil=True, cache=True, fastmath=True)
+def _aradon_inner_2d(X, Y, f, px, h, flim0, flim1, npx, nh):
+    for ipx in prange(npx):
+        for ifr in range(flim0, flim1):
+            for ih in range(nh):
+                X[ipx, ifr] += Y[ih, ifr] * np.exp(
+                    1j * 2 * np.pi * f[ifr] * px[ipx] * h[ih]
+                )
+    return X