Skip to content

Commit be916b1

Browse files
mrava87mrava87
authored
Merge pull request #773 from mrava87/feat-udct
Feat: UDCT operator
2 parents 7204a71 + 468d784 commit be916b1

19 files changed

Lines changed: 511 additions & 106 deletions

‎.github/workflows/buildft.yaml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ jobs:
2929
- name: Install dependencies and pylops
3030
run: uv sync --locked --extra deep --all-groups
3131
- name: Add more dependencies
32-
run: uv add PyWavelets spgl1 devito
32+
run: uv add PyWavelets spgl1 curvelets devito
3333
- name: Verify free threading
3434
run: |
3535
uv run python - <<'PY'

‎docs/source/api/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,7 @@ Signal processing
112112
DCT
113113
DTCWT
114114
Seislet
115+
UDCT
115116
Radon2D
116117
Radon3D
117118
FourierRadon2D

‎docs/source/gpu.rst‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -457,6 +457,10 @@ Signal processing:
457457
- |:white_check_mark:|
458458
- |:red_circle:|
459459
- |:red_circle:|
460+
* - :class:`pylops.signalprocessing.UDCT`
461+
- |:white_check_mark:|
462+
- |:red_circle:|
463+
- |:red_circle:|
460464
* - :class:`pylops.signalprocessing.Radon2D`
461465
- |:white_check_mark:|
462466
- |:red_circle:|

‎docs/source/installation.rst‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -433,14 +433,18 @@ them as part of the installation process of PyLops provided in the table below.
433433
- :py:class:`pylops.medical.CT2D`
434434
- |:white_check_mark:|
435435
- `pip install pylops[advanced]` / `uv add "pylops[advanced]"`
436-
* - dtcwt
437-
- :py:class:`pylops.signalprocessing.DTCWT`
436+
* - Curvelets
437+
- :py:class:`pylops.signalprocessing.UDCT`
438438
- |:white_check_mark:|
439439
- `pip install pylops[advanced]` / `uv add "pylops[advanced]"`
440440
* - Devito
441441
- :py:class:`pylops.waveeqprocessing.AcousticWave2D`
442442
- |:white_check_mark:|
443443
- `pip install pylops[advanced]` / `uv add "pylops[advanced]"`
444+
* - dtcwt
445+
- :py:class:`pylops.signalprocessing.DTCWT`
446+
- |:white_check_mark:|
447+
- `pip install pylops[advanced]` / `uv add "pylops[advanced]"`
444448
* - FFTW
445449
- :py:class:`pylops.signalprocessing.FFT`, :py:class:`pylops.signalprocessing.FFT2D`,
446450
:py:class:`pylops.signalprocessing.FFTND`

‎environment-dev-arm.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ dependencies:
2424
- autopep8
2525
- pip:
2626
- torch
27+
- curvelets>=1.2
2728
- devito
2829
# - dtcwt (until numpy>=2.0.0 is supported)
2930
- scikit-fmm

‎environment-dev-intel-mkl.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ dependencies:
2727
- mkl_fft
2828
- pip:
2929
- torch
30+
- curvelets>=1.2
3031
- devito
3132
# - dtcwt (until numpy>=2.0.0 is supported)
3233
- scikit-fmm

‎environment-dev.yml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,7 @@ dependencies:
2525
- autopep8
2626
- pip:
2727
- torch
28+
- curvelets>=1.2
2829
- devito
2930
# - dtcwt (until numpy>=2.0.0 is supported)
3031
- scikit-fmm

‎examples/plot_udct.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
"""
2+
Uniform Discrete Curvelet Transform
3+
===================================
4+
5+
This example shows how to use the :py:class:`pylops.signalprocessing.UDCT` operator.
6+
This operator performs the Uniform Discrete Curvelet Transform along all or a portion
7+
of the axes of multi-dimensional input array.
8+
9+
"""
10+
11+
from math import prod
12+
13+
import matplotlib.pyplot as plt
14+
import numpy as np
15+
from curvelets.utils import make_zone_plate
16+
17+
import pylops
18+
19+
plt.close("all")
20+
21+
###############################################################################
22+
# To begin with consider the same input of the first tutorial of the
23+
# `curvelets <https://curvelets.readthedocs.io/en/stable/auto_examples/plot_01_zone_plate.html>`_
24+
# library. However, we will use here the PyLops operator to apply the 2D UDCT.
25+
# A round-trip (i.e., forward-adjoing) is performed and the output is shown to be
26+
# equal to the input, meaning that the adjoint of this transform is also the inverse.
27+
28+
shape = (256, 256)
29+
x = make_zone_plate(shape)
30+
31+
Cop = pylops.signalprocessing.UDCT(shape)
32+
33+
y = Cop @ x
34+
xinv = Cop.H @ y
35+
36+
fig, axs = plt.subplots(1, 3, figsize=(10, 4))
37+
axs[0].imshow(x, cmap="gray", vmin=-1, vmax=1)
38+
axs[0].set_title("Input")
39+
axs[1].imshow(xinv, cmap="gray", vmin=-1, vmax=1)
40+
axs[1].set_title("Round-trip")
41+
axs[2].imshow(x - xinv, cmap="gray", vmin=-1, vmax=1)
42+
axs[2].set_title("Difference")
43+
for ax in axs:
44+
ax.axis("tight")
45+
ax.axis("off")
46+
fig.tight_layout()
47+
48+
49+
###############################################################################
50+
# Next, the 2D version of the UDCT can also be applied along the last two axes
51+
# of a N-dimensional input. This can be easily done with the aid of the
52+
# :py:class:`pylops.Kronecker` operator.
53+
54+
shape = (2, 3, 64, 32)
55+
x = make_zone_plate(shape[2:])
56+
57+
# replicate over the first two axes.
58+
x = np.tile(x[None, None], (shape[0], shape[1], 1, 1))
59+
60+
Cop = pylops.signalprocessing.UDCT(dims=shape[2:])
61+
Iop = pylops.Identity(prod(shape[:2]))
62+
Cop = pylops.Kronecker(Iop, Cop)
63+
64+
y = Cop @ x.ravel()
65+
xinv = (Cop.H @ y).reshape(shape)
66+
67+
fig, axs = plt.subplots(1, 3, figsize=(10, 4))
68+
axs[0].imshow(np.mean(x, axis=(0, 1)), cmap="gray", vmin=-1, vmax=1)
69+
axs[0].set_title("Input")
70+
axs[1].imshow(np.mean(xinv, axis=(0, 1)), cmap="gray", vmin=-1, vmax=1)
71+
axs[1].set_title("Round-trip")
72+
axs[2].imshow(
73+
np.mean(x, axis=(0, 1)) - np.mean(xinv, axis=(0, 1)), cmap="gray", vmin=-1, vmax=1
74+
)
75+
axs[2].set_title("Difference")
76+
for ax in axs:
77+
ax.axis("tight")
78+
ax.axis("off")
79+
fig.tight_layout()
80+
81+
82+
###############################################################################
83+
# And similarly, if we want to apply the UDCT over any two axes, the
84+
# :py:class:`pylops.Transpose` can be used to re-organize the
85+
# axes of the input prior to applying the transform
86+
shape = (2, 64, 3, 32)
87+
x = make_zone_plate((shape[1], shape[3]))
88+
89+
# replicate over the first two axes.
90+
x = np.tile(x[None, :, None], (shape[0], 1, shape[2], 1))
91+
92+
Cop = pylops.signalprocessing.UDCT(dims=(shape[1], shape[3]))
93+
Top = pylops.Transpose(shape, axes=(0, 2, 1, 3))
94+
Iop = pylops.Identity(prod((shape[0], shape[2])))
95+
Cop = pylops.Kronecker(Iop, Cop) @ Top
96+
97+
y = Cop @ x.ravel()
98+
xinv = (Cop.H @ y).reshape(shape)
99+
100+
fig, axs = plt.subplots(1, 3, figsize=(10, 4))
101+
axs[0].imshow(np.mean(x, axis=(0, 2)), cmap="gray", vmin=-1, vmax=1)
102+
axs[0].set_title("Input")
103+
axs[1].imshow(np.mean(xinv, axis=(0, 2)), cmap="gray", vmin=-1, vmax=1)
104+
axs[1].set_title("Round-trip")
105+
axs[2].imshow(
106+
np.mean(x, axis=(0, 2)) - np.mean(xinv, axis=(0, 2)), cmap="gray", vmin=-1, vmax=1
107+
)
108+
axs[2].set_title("Difference")
109+
for ax in axs:
110+
ax.axis("tight")
111+
ax.axis("off")
112+
fig.tight_layout()
113+
114+
115+
###############################################################################
116+
# Finally, if the :py:class:`pylops.signalprocessing.UDCT` is applied directly
117+
# to a 3-dimensional (or N-dimensional) input, the 3D/ND version of the
118+
# transform is used.
119+
120+
shape = (16, 32, 16)
121+
x = make_zone_plate(shape[1:])
122+
123+
# replicate over the first axis.
124+
x = np.tile(x[None], (shape[0], 1, 1))
125+
126+
Cop = pylops.signalprocessing.UDCT(dims=shape)
127+
128+
y = Cop @ x
129+
xinv = Cop.H @ y
130+
131+
fig, axs = plt.subplots(1, 3, figsize=(10, 4))
132+
axs[0].imshow(np.mean(x, axis=0), cmap="gray", vmin=-1, vmax=1)
133+
axs[0].set_title("Input")
134+
axs[1].imshow(np.mean(xinv, axis=0), cmap="gray", vmin=-1, vmax=1)
135+
axs[1].set_title("Round-trip")
136+
axs[2].imshow(np.mean(x, axis=0) - np.mean(xinv, axis=0), cmap="gray", vmin=-1, vmax=1)
137+
axs[2].set_title("Difference")
138+
for ax in axs:
139+
ax.axis("tight")
140+
ax.axis("off")
141+
fig.tight_layout()

‎pylops/linearoperator.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1439,7 +1439,7 @@ def __init__(self, A: LinearOperator, B: LinearOperator):
14391439
)
14401440
raise ValueError(msg)
14411441
if A.shape[1] != B.shape[0]:
1442-
msg = f"Cannot add {A} and {B}: shape mismatch"
1442+
msg = f"Cannot multiply {A} and {B}: shape mismatch"
14431443
raise ValueError(msg)
14441444
super().__init__(dtype=_get_dtype([A, B]), shape=(A.shape[0], B.shape[1]))
14451445
self.args = (A, B)

‎pylops/signalprocessing/__init__.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
DWTND N-dimensional Wavelet operator.
2828
DCT Discrete Cosine Transform.
2929
DTCWT Dual-Tree Complex Wavelet Transform.
30+
UDCT Uniform Discrete Curvelet Transform.
3031
Radon2D Two dimensional Radon transform.
3132
Radon3D Three dimensional Radon transform.
3233
FourierRadon2D Two dimensional Fourier Radon transform.
@@ -77,6 +78,7 @@
7778
from .dct import *
7879
from .dtcwt import *
7980
from .pwd2d import *
81+
from .udct import *
8082

8183
__all__ = [
8284
"FFT",
@@ -114,4 +116,5 @@
114116
"DTCWT",
115117
"PWSprayer2D",
116118
"PWSmoother2D",
119+
"UDCT",
117120
]

0 commit comments

Comments
 (0)