Skip to content

Commit

Permalink
test: added tests for FourierRadon operators
Browse files Browse the repository at this point in the history
  • Loading branch information
@mrava87
mrava87 committed Sep 6, 2024
1 parent 69ca796 commit 66f4b0f
Show file tree
Hide file tree
Showing 6 changed files with 292 additions and 103 deletions.
8 changes: 8 additions & 0 deletions docs/source/gpu.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -281,6 +281,14 @@ Signal processing:
- |:white_check_mark:|
- |:red_circle:|
- |:red_circle:|
* - :class:`pylops.signalprocessing.FourierRadon2D`
- |:white_check_mark:|
- |:white_check_mark:|
- |:red_circle:|
* - :class:`pylops.signalprocessing.FourierRadon3D`
- |:white_check_mark:|
- |:white_check_mark:|
- |:red_circle:|
* - :class:`pylops.signalprocessing.ChirpRadon2D`
- |:white_check_mark:|
- |:white_check_mark:|
Expand Down
238 changes: 216 additions & 22 deletions examples/plot_fourierradon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,14 @@
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.
and :py:class:`pylops.signalprocessing.FourierRadon3D` operators to apply the linear
and parabolic Radon Transform to 2-dimensional or 3-dimensional signals, respectively.
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).
These operators provides transformations equivalent to those of
:py:class:`pylops.signalprocessing.Radon2D` and :py:class:`pylops.signalprocessing.Radon3D`,
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
Expand Down Expand Up @@ -49,7 +51,7 @@


###############################################################################
# We can now define our operators and apply the forward, adjoint and inverse
# We can now define our operators and apply the forward and adjoint
# steps.
nfft = int(2 ** np.ceil(np.log2(par["nt"])))
npx, pxmax = 2 * par["nx"], 5e-4
Expand All @@ -58,28 +60,28 @@
R2Op = pylops.signalprocessing.FourierRadon2D(
t, x, px, nfft, kind="linear", engine="numpy", dtype="float64"
)
dL_chirp = R2Op.H * d
dadj_chirp = R2Op * dL_chirp
dL = R2Op.H * d
dadj = R2Op * dL

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,
dL.T,
cmap="bwr_r",
vmin=-dL_chirp.max(),
vmax=dL_chirp.max(),
vmin=-dL.max(),
vmax=dL.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,
dadj.T,
cmap="bwr_r",
vmin=-dadj_chirp.max(),
vmax=dadj_chirp.max(),
vmin=-dadj.max(),
vmax=dadj.max(),
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon")
Expand All @@ -101,30 +103,222 @@
R2Op = pylops.signalprocessing.FourierRadon2D(
t, x, px, nfft, kind="parabolic", engine="numpy", dtype="float64"
)
dL_chirp = R2Op.H * d
dadj_chirp = R2Op * dL_chirp
dL = R2Op.H * d
dadj = R2Op * dL

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,
dL.T,
cmap="bwr_r",
vmin=-dL_chirp.max(),
vmax=dL_chirp.max(),
vmin=-dL.max(),
vmax=dL.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,
dadj.T,
cmap="bwr_r",
vmin=-dadj_chirp.max(),
vmax=dadj_chirp.max(),
vmin=-dadj.max(),
vmax=dadj.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()

###############################################################################
# Finally we repeat the same exercise with 3d data.

par = {
"ot": 0,
"dt": 0.004,
"nt": 51,
"ox": -100,
"dx": 10,
"nx": 21,
"oy": -200,
"dy": 10,
"ny": 41,
"f0": 20,
}
theta = [30]
phi = [10]
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"]

# Generate linear data
pxx = np.sin(np.deg2rad(theta[0])) * np.cos(np.deg2rad(phi[0])) / 1500.0
pyy = np.sin(np.deg2rad(theta[0])) * np.sin(np.deg2rad(phi[0])) / 1500.0
_, d = pylops.utils.seismicevents.linear3d(x, y, t, 1500.0, t0, theta, phi, amp, wav)

# Linear Radon
npy, pymax = par["ny"], 5e-4
npx, pxmax = par["nx"], 5e-4
py = np.linspace(-pymax, pymax, npy)
px = np.linspace(-pxmax, pxmax, npx)

R3Op = pylops.signalprocessing.FourierRadon3D(
t, y, x, py, px, nfft, kind=("linear", "linear"), engine="numpy", dtype="float64"
)
dL = R3Op.H * d
dadj = R3Op * dL

fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
axs[0].imshow(
d[par["ny"] // 2].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 3d - y")
axs[0].axis("tight")
axs[1].imshow(
dL[np.argmin(np.abs(pyy - py))].T,
cmap="bwr_r",
vmin=-dL.max(),
vmax=dL.max(),
extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),
)
axs[1].scatter(1e3 * pxx, t0[0], s=50, color="r")
axs[1].set(xlabel=r"$p_x$ [s/km]", title="Radon 3d - y")
axs[1].axis("tight")
axs[2].imshow(
dadj[par["ny"] // 2].T,
cmap="bwr_r",
vmin=-dadj.max(),
vmax=dadj.max(),
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon 3d - y")
axs[2].axis("tight")
plt.tight_layout()

fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
axs[0].imshow(
d[:, par["nx"] // 2].T,
vmin=-1,
vmax=1,
cmap="bwr_r",
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[0].set(xlabel=r"$y$ [m]", ylabel=r"$t$ [s]", title="Input linear 3d - x")
axs[0].axis("tight")
axs[1].imshow(
dL[:, np.argmin(np.abs(pxx - px))].T,
cmap="bwr_r",
vmin=-dL.max(),
vmax=dL.max(),
extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),
)
axs[1].scatter(1e3 * pyy, t0[0], s=50, color="r")
axs[1].set(xlabel=r"$p_y$ [s/km]", title="Radon 3d - x")
axs[1].axis("tight")
axs[2].imshow(
dadj[:, par["nx"] // 2].T,
cmap="bwr_r",
vmin=-dadj.max(),
vmax=dadj.max(),
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[2].set(xlabel=r"$y$ [m]", title="Adj Radon 3d - x")
axs[2].axis("tight")
plt.tight_layout()

# Generate parabolic data
pxx = [1e-6]
pyy = [2e-6]
_, d = pylops.utils.seismicevents.parabolic3d(
x, y, t, t0, 0, 0, np.array(pxx), np.array(pyy), amp, wav
)

# Parabolic Radon
npy, pymax = par["ny"], 5e-6
npx, pxmax = par["nx"], 5e-6
py = np.linspace(-pymax, pymax, npy)
px = np.linspace(-pxmax, pxmax, npx)

R3Op = pylops.signalprocessing.FourierRadon3D(
t,
y,
x,
py,
px,
nfft,
kind=("parabolic", "parabolic"),
engine="numpy",
dtype="float64",
)
dL = R3Op.H * d
dadj = R3Op * dL

fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
axs[0].imshow(
d[par["ny"] // 2].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 3d - y")
axs[0].axis("tight")
axs[1].imshow(
dL[np.argmin(np.abs(pyy - py))].T,
cmap="bwr_r",
vmin=-dL.max(),
vmax=dL.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_x$ [s/km]", title="Radon 3d - y")
axs[1].axis("tight")
axs[2].imshow(
dadj[par["ny"] // 2].T,
cmap="bwr_r",
vmin=-dadj.max(),
vmax=dadj.max(),
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon 3d - y")
axs[2].axis("tight")
plt.tight_layout()

fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
axs[0].imshow(
d[:, par["nx"] // 2].T,
vmin=-1,
vmax=1,
cmap="bwr_r",
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[0].set(xlabel=r"$y$ [m]", ylabel=r"$t$ [s]", title="Input parabolic 3d - x")
axs[0].axis("tight")
axs[1].imshow(
dL[:, np.argmin(np.abs(pxx - px))].T,
cmap="bwr_r",
vmin=-dL.max(),
vmax=dL.max(),
extent=(1e3 * py[0], 1e3 * py[-1], t[-1], t[0]),
)
axs[1].scatter(1e3 * pyy[0], t0[0], s=50, color="r")
axs[1].set(xlabel=r"$p_y$ [s/km]", title="Radon 3d - x")
axs[1].axis("tight")
axs[2].imshow(
dadj[:, par["nx"] // 2].T,
cmap="bwr_r",
vmin=-dadj.max(),
vmax=dadj.max(),
extent=(x[0], x[-1], t[-1], t[0]),
)
axs[2].set(xlabel=r"$y$ [m]", title="Adj Radon 3d - x")
axs[2].axis("tight")
plt.tight_layout()
2 changes: 2 additions & 0 deletions pylops/signalprocessing/fourierradon2d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,7 @@ def _register_multiplications(self, engine: str) -> None:
@reshaped
def _matvec_numpy(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
self.f = ncp.asarray(self.f)
x = ncp.fft.rfft(x, n=self.nfft, axis=-1)

H, PX, F = ncp.meshgrid(
Expand All @@ -201,6 +202,7 @@ def _matvec_numpy(self, x: NDArray) -> NDArray:
@reshaped
def _rmatvec_numpy(self, y: NDArray) -> NDArray:
ncp = get_array_module(y)
self.f = ncp.asarray(self.f)
y = ncp.fft.rfft(y, n=self.nfft, axis=-1)

PX, H, F = ncp.meshgrid(
Expand Down
Loading

0 comments on commit 66f4b0f

Please sign in to comment.