feat: added FourierRadon3D

docs/source/api/index.rst

@@ -110,6 +110,7 @@ Signal processing
     Radon2D
     Radon3D
     FourierRadon2D
+    FourierRadon3D
     ChirpRadon2D
     ChirpRadon3D
     Sliding1D

pylops/signalprocessing/__init__.py

@@ -56,6 +56,7 @@
 from .radon2d import *
 from .radon3d import *
 from .fourierradon2d import *
+from .fourierradon3d import *
 from .chirpradon2d import *
 from .chirpradon3d import *
 from .sliding1d import *
@@ -89,6 +90,8 @@
     "Bilinear",
     "Radon2D",
     "Radon3D",
+    "FourierRadon2D",
+    "FourierRadon3D",
     "ChirpRadon2D",
     "ChirpRadon3D",
     "Sliding1D",

pylops/signalprocessing/_fourierradon3d_numba.py

@@ -0,0 +1,42 @@
+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_3d(X, Y, f, py, px, hy, hx, flim0, flim1, npy, npx, nhy, nhx):
+    for ihy in prange(nhy):
+        for ihx in prange(nhx):
+            for ifr in range(flim0, flim1):
+                for ipy in range(npy):
+                    for ipx in range(npx):
+                        Y[ihy, ihx, ifr] += X[ipy, ipx, ifr] * np.exp(
+                            -1j
+                            * 2
+                            * np.pi
+                            * f[ifr]
+                            * (py[ipy] * hy[ihy] + px[ipx] * hx[ihx])
+                        )
+    return Y
+
+
+@jit(nopython=True, parallel=parallel, nogil=True, cache=True, fastmath=True)
+def _aradon_inner_3d(X, Y, f, py, px, hy, hx, flim0, flim1, npy, npx, nhy, nhx):
+    for ipy in prange(npy):
+        for ipx in range(npx):
+            for ifr in range(flim0, flim1):
+                for ihy in range(nhy):
+                    for ihx in range(nhx):
+                        X[ipy, ipx, ifr] += Y[ihy, ihx, ifr] * np.exp(
+                            1j
+                            * 2
+                            * np.pi
+                            * f[ifr]
+                            * (py[ipy] * hy[ihy] + px[ipx] * hx[ihx])
+                        )
+    return X

pylops/signalprocessing/_nonstatconvolve2d_cuda.py

@@ -100,9 +100,9 @@ def _matvec_rmatvec_call(
 ):
     """Caller for NonStationaryConvolve2D operator
 
-    Caller for cuda implementation of matvec and rmatvec for NonStationaryConvolve2D operato, with same signature
+    Caller for cuda implementation of matvec and rmatvec for NonStationaryConvolve2D operator, with same signature
     as numpy/numba counterparts. See :class:`pylops.signalprocessing.NonStationaryConvolve2D` for details about
-     input parameters.
+    input parameters.
 
     """
     _matvec_rmatvec[num_blocks, num_threads_per_blocks](

pylops/signalprocessing/fourierradon2d.py

@@ -26,8 +26,8 @@ class FourierRadon2D(LinearOperator):
     r"""2D Fourier Radon transform
 
     Apply Radon forward (and adjoint) transform using Fast
-    Fourier Transform to a 2-dimensional array of size
-    :math:`[n_x \times n_t]` (both in forward and adjoint mode).
+    Fourier Transform to a 2-dimensional array of size :math:`[n_{p_x} \times n_t]`
+    (and :math:`[n_x \times n_t]`).
 
     Note that forward and adjoint follow the same convention of the time-space
     implementation in :class:`pylops.signalprocessing.Radon2D`.
@@ -73,23 +73,23 @@ class FourierRadon2D(LinearOperator):
 
     .. math::
         \begin{bmatrix}
-            \mathbf{m}(p_{x,1}, \omega_i)  \\
-            \mathbf{m}(p_{x,2}, \omega_i)  \\
+            m(p_{x,1}, \omega_i)  \\
+            m(p_{x,2}, \omega_i)  \\
             \vdots          \\
-            \mathbf{m}(p_{x,N_p}, \omega_i)
+            m(p_{x,N_p}, \omega_i)
         \end{bmatrix}
         =
         \begin{bmatrix}
-            e^{-j \omega_i p_{x,1} x^l_1}  & e^{-j \omega_i p_{x,1} x^l_2} &  \ldots & e^{-j \omega_i p_{x,2} x^l_{N_x}}  \\
+            e^{-j \omega_i p_{x,1} x^l_1}  & e^{-j \omega_i p_{x,1} x^l_2} &  \ldots & e^{-j \omega_i p_{x,1} x^l_{N_x}}  \\
             e^{-j \omega_i p_{x,2} x^l_1}  & e^{-j \omega_i p_{x,2} x^l_2} &  \ldots & e^{-j \omega_i p_{x,2} x^l_{N_x}}  \\
             \vdots            & \vdots           &  \ddots & \vdots            \\
             e^{-j \omega_i p_{x,N_p} x^l_1}  & e^{-j \omega_i p_{x,N_p} x^l_2} &  \ldots & e^{-j \omega_i p_{x,N_p} x^l_{N_x}}  \\
         \end{bmatrix}
         \begin{bmatrix}
-            \mathbf{d}(x_1, \omega_i)  \\
-            \mathbf{d}(x_2, \omega_i)  \\
+            d(x_1, \omega_i)  \\
+            d(x_2, \omega_i)  \\
             \vdots          \\
-            \mathbf{d}(x_{N_x}, \omega_i)
+            d(x_{N_x}, \omega_i)
         \end{bmatrix}
 
     where :math:`l=1,2`. Similarly the forward mode is implemented by applying the