How to use the pylops.signalprocessing.Radon2D function in pylops

np.abs(taxis[1] - taxis[1]))
            Pop = FFT2D(dims=dims, nffts=nffts,
                        sampling=sampling)
            Pop = Pop.H
        SIop = Rop * Pop
    elif 'radon' in kind:
        prec = True
        dotcflag = 0
        kindradon = kind.split('-')[-1]
        if ndims == 3:
            Pop = Radon3D(taxis, spataxis, spat1axis, paxis, p1axis,
                          centeredh=centeredh, kind=kindradon, engine=engine)
            dimsp = (paxis.size, p1axis.size, taxis.size)

        else:
            Pop = Radon2D(taxis, spataxis, paxis,
                          centeredh=centeredh, kind=kindradon, engine=engine)
            dimsp = (paxis.size, taxis.size)
        SIop = Rop * Pop
    elif kind == 'sliding':
        prec = True
        dotcflag = 0
        if ndims == 3:
            dspat = np.abs(spataxis[1]-spataxis[0])
            dspat1 = np.abs(spat1axis[1]-spat1axis[0])
            nspat, nspat1 = spataxis.size, spat1axis.size
            spataxis_local = np.linspace(-dspat * nwin[0] // 2,
                                         dspat * nwin[0] // 2,
                                         nwin[0])
            spat1axis_local = np.linspace(-dspat1 * nwin[1] // 2,
                                          dspat1 * nwin[1] // 2,
                                          nwin[1])

def test_Radon2D(par):
    """Dot-test and sparse inverse for Radon2D operator
    """
    dt, dh = 0.005, 1
    t = np.arange(par['nt']) * dt
    h = np.arange(par['nhx']) * dh
    px = np.linspace(0, par['pxmax'], par['npx'])
    x = np.zeros((par['npx'], par['nt']))
    x[2, par['nt']//2] = 1

    Rop = Radon2D(t, h, px, centeredh=par['centeredh'],
                  interp=par['interp'], kind=par['kind'],
                  onthefly=par['onthefly'], engine=par['engine'],
                  dtype='float64')
    assert dottest(Rop, par['nhx']*par['nt'], par['npx']*par['nt'],
                   complexflag=0)
    y = Rop * x.flatten()
    y = y.reshape(par['nhx'], par['nt'])

    xinv, _, _ = FISTA(Rop, y.flatten(), 30, eps=1e0, returninfo=True)
    assert_array_almost_equal(x.flatten(), xinv, decimal=1)

# where :math:`\theta` is the angle between the x-axis (:math:`x`) and
# the perpendicular to the summation line and :math:`r` is the distance
# from the origin of the summation line.
@jit(nopython=True)
def radoncurve(x, r, theta):
    return (r - ny//2)/(np.sin(np.deg2rad(theta))+1e-15) + np.tan(np.deg2rad(90 - theta))*x  + ny//2

x = np.load('../testdata/optimization/shepp_logan_phantom.npy')
x = x / x.max()
ny, nx = x.shape

ntheta = 150
theta = np.linspace(0., 180., ntheta, endpoint=False)

RLop = \
    pylops.signalprocessing.Radon2D(np.arange(ny), np.arange(nx),
                                    theta, kind=radoncurve,
                                    centeredh=True, interp=False,
                                    engine='numba', dtype='float64')

y = RLop.H * x.T.ravel()
y = y.reshape(ntheta, ny).T

###############################################################################
# We can now first perform the adjoint, which in the medical imaging literature
# is also referred to as back-projection.
#
# This is the first step of a common reconstruction technique, named filtered
# back-projection, which simply applies a correction filter in the
# frequency domain to the adjoint model.
xrec = RLop*y.T.ravel()
xrec = xrec.reshape(nx, ny).T

nt, nh = 41, 51
npx, pxmax = 41, 1e-2

dt, dh = 0.005, 1
t = np.arange(nt) * dt
h = np.arange(nh) * dh
px = np.linspace(0, pxmax, npx)

x = np.zeros((npx, nt))
x[4, nt//2] = 1

###############################################################################
# We can now define our operators for different parametric curves and apply
# them to the input model vector. We also apply the adjoint to the resulting
# data vector.
RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='linear', interp=False,
                                       engine='numpy')
RPop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='parabolic', interp=False,
                                       engine='numpy')
RHop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='hyperbolic', interp=False,
                                       engine='numpy')

# forward
yL = RLop * x.flatten()
yP = RPop * x.flatten()
yH = RHop * x.flatten()
yL = yL.reshape(nh, nt)
yP = yP.reshape(nh, nt)
yH = yH.reshape(nh, nt)

# patches in the spatial direction and apply the adjoint of the
# :py:class:`pylops.signalprocessing.Radon2D` operator to each patch. This is
# done by simply using the adjoint of the
# :py:class:`pylops.signalprocessing.Sliding2D` operator
nwins = 5
winsize = 36
overlap = 10
npx = 61
px = np.linspace(-5e-3, 5e-3, npx)

dimsd = data.shape
dims = (nwins*npx, par['nt'])

# Sliding window transform without taper
Op = \
    pylops.signalprocessing.Radon2D(t, np.linspace(-par['dx']*winsize//2,
                                                   par['dx']*winsize//2,
                                                   winsize),
                                    px, centeredh=True, kind='linear',
                                    engine='numba')
Slid = pylops.signalprocessing.Sliding2D(Op, dims, dimsd,
                                         winsize, overlap,
                                         tapertype=None)

radon = Slid.H * data.flatten()
radon = radon.reshape(dims)

###############################################################################
# We now create  a similar operator but we also add a taper to the overlapping
# parts of the patches.
Slid = pylops.signalprocessing.Sliding2D(Op, dims, dimsd,
                                         winsize, overlap,

px = np.linspace(0, pxmax, npx)

x = np.zeros((npx, nt))
x[4, nt//2] = 1

###############################################################################
# We can now define our operators for different parametric curves and apply
# them to the input model vector. We also apply the adjoint to the resulting
# data vector.
RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='linear', interp=False,
                                       engine='numpy')
RPop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='parabolic', interp=False,
                                       engine='numpy')
RHop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='hyperbolic', interp=False,
                                       engine='numpy')

# forward
yL = RLop * x.flatten()
yP = RPop * x.flatten()
yH = RHop * x.flatten()
yL = yL.reshape(nh, nt)
yP = yP.reshape(nh, nt)
yH = yH.reshape(nh, nt)

# adjoint
xadjL = RLop.H * yL.flatten()
xadjP = RPop.H * yP.flatten()
xadjH = RHop.H * yH.flatten()
xadjL = xadjL.reshape(npx, nt)

# We can now create the :py:class:`pylops.signalprocessing.Radon2D` operator.
# We also apply its adjoint to the data to obtain a representation of those
# 3 linear events overlapping to a parabolic event in the Radon domain.
# Similarly, we feed the operator to a sparse solver like
# :py:class:`pylops.optimization.sparsity.FISTA` to obtain a sparse
# represention of the data in the Radon domain. At this point we try to filter
# out the unwanted event. We can see how this is much easier for the sparse
# transform as each event has a much more compact representation in the Radon
# domain than for the adjoint transform.

# radon operator
npx = 61
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)

Rop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, kind='linear',
                                      interp='nearest',  centeredh=False,
                                      dtype='float64')

# adjoint Radon transform
xadj = Rop.H * y.flatten()
xadj = xadj.reshape(npx, par['nt'])

# sparse Radon transform
xinv, niter, cost = \
    pylops.optimization.sparsity.FISTA(Rop, y.flatten(), 15,
                                       eps=1e1, returninfo=True)
xinv = xinv.reshape(npx, par['nt'])

# filtering
xfilt = np.zeros_like(xadj)
xfilt[npx//2-3:npx//2+4] = xadj[npx//2-3:npx//2+4]

dt, dh = 0.005, 1
t = np.arange(nt) * dt
h = np.arange(nh) * dh
px = np.linspace(0, pxmax, npx)

x = np.zeros((npx, nt))
x[4, nt//2] = 1

###############################################################################
# We can now define our operators for different parametric curves and apply
# them to the input model vector. We also apply the adjoint to the resulting
# data vector.
RLop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='linear', interp=False,
                                       engine='numpy')
RPop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='parabolic', interp=False,
                                       engine='numpy')
RHop = pylops.signalprocessing.Radon2D(t, h, px, centeredh=True,
                                       kind='hyperbolic', interp=False,
                                       engine='numpy')

# forward
yL = RLop * x.flatten()
yP = RPop * x.flatten()
yH = RHop * x.flatten()
yL = yL.reshape(nh, nt)
yP = yP.reshape(nh, nt)
yH = yH.reshape(nh, nt)

# adjoint
xadjL = RLop.H * yL.flatten()

axs[1].set_xlim(-0.1, 0.1)
axs[1].set_ylim(50, 0)
axs[2].plot(cost, 'k', lw=3)
axs[2].set_title('FISTA convergence')

###############################################################################
# We see how adding prior information to the inversion can help improving the
# estimate of the regularized seismic data. Nevertheless, in both cases the
# reconstructed data is not perfect. A better sparsyfing transform could in
# fact be chosen here to be the linear
# :py:class:`pylops.signalprocessing.Radon2D` transform in spite of the
# :py:class:`pylops.FFT2` transform.
npx = 40
pxmax = 1e-3
px = np.linspace(-pxmax, pxmax, npx)
Radop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, engine='numba')

RRop = Rop*Radop

# adjoint
Xadj_fromx = Radop.H*x.flatten()
Xadj_fromx = Xadj_fromx.reshape(npx, par['nt'])

Xadj = RRop.H*y.flatten()
Xadj = Xadj.reshape(npx, par['nt'])

# L1 inverse
xl1, Xl1, cost = \
    pylops.waveeqprocessing.SeismicInterpolation(y, par['nx'], iava,
                                                 kind='radon-linear',
                                                 spataxis=xaxis,
                                                 taxis=taxis, paxis=px,