### Install Pre-commit Hooks

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Install pre-commit hooks to ensure code quality and consistency before committing changes.

```console
$ pre-commit install
```

--------------------------------

### Clone Repository and Install Dev Requirements

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Clone the main branch of the repository and install development requirements using make.

```console
$ git clone https://github.com/PyLops/curvelops
$ git remote add upstream https://github.com/PyLops/curvelops
$ make dev-install
```

--------------------------------

### Install FFTW

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/installation.rst

Steps to download, extract, configure, compile, and install the FFTW library. Ensure you use the same compiler for FFTW and CurveLab.

```console
$ wget https://www.fftw.org/fftw-2.1.5.tar.gz
$ tar xvzf fftw-2.1.5.tar.gz
$ mkdir -p /home/$USER/opt/
$ mv fftw-2.1.5/ /home/$USER/opt/
$ cd /home/$USER/opt/fftw-2.1.5/
$ ./configure --with-pic --prefix=/home/$USER/opt/fftw-2.1.5 --with-gcc=$(which gcc)
$ make
$ make install
```

```console
$ make check
```

--------------------------------

### Install Curvelops

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/installation.rst

Install Curvelops using pip after setting FFTW and FDCT environment variables. FFTW should point to the FFTW installation directory, and FDCT to the CurveLab root.

```console
$ export FFTW=/path/to/fftw-2.1.5
$ export FDCT=/path/to/CurveLab-2.1.3
$ python3 -m pip install git+https://github.com/PyLops/curvelops@0.23.4
```

--------------------------------

### Install CurveLab

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/installation.rst

Instructions to extract, move, configure, compile, and test the CurveLab library. Set FFTW_DIR, CC, and CXX variables in makefile.opt.

```console
$ tar xvzf CurveLab-2.1.3.tar.gz
$ mkdir -p /home/$USER/opt/
$ mv CurveLab-2.1.3/ /home/$USER/opt/
$ cd /home/$USER/opt/CurveLab-2.1.3/
$ cp makefile.opt makefile.opt.bak
```

```console
$ cd /home/$USER/opt/CurveLab-2.1.3/
$ make clean
$ make lib
```

```console
$ make test
```

--------------------------------

### Install Curvelops in Editable Mode

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Install Curvelops in editable mode using pip after setting up dependencies in a separate environment.

```console
$ python3 -m pip install -e .
```

--------------------------------

### Install Curvelops with Pip

Source: https://github.com/pylops/curvelops/blob/main/README.md

Install curvelops using pip after setting environment variables for FFTW and FDCT. Ensure pip version is 10.0 or higher.

```bash
export FFTW=/path/to/fftw-2.1.5
export FDCT=/path/to/CurveLab-2.1.3
python3 -m pip install git+https://github.com/PyLops/curvelops@0.23.4
```

--------------------------------

### Compute and verify k-space vector

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Example usage of find_kmax to extract the direction vector and compare it against the ground truth.

```python
kvec = find_kmax(img_k, kx, kz)
kvec / np.linalg.norm(kvec)
```

```python
vec
```

--------------------------------

### Import necessary libraries for pylops

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Imports essential libraries including numpy, matplotlib, and specific modules from pylops for signal processing and utility functions. This setup is required for the subsequent examples.

```python
from math import ceil, floor

import matplotlib.pyplot as plt
import numpy as np
from pylops.signalprocessing import FFT2D
from pylops.utils.tapers import taper2d
from scipy.signal import filtfilt

plt.rcParams.update({"image.interpolation": "blackman"})
```

--------------------------------

### 2D Curvelet Transform Example

Source: https://github.com/pylops/curvelops/blob/main/README.md

Demonstrates a 2D curvelet transform and its inverse using curvelops. Requires numpy and curvelops libraries. Asserts that the reconstructed signal is close to the original.

```python
import numpy as np
import curvelops as cl

x = np.random.randn(100, 50)
FDCT = cl.FDCT2D(dims=x.shape)
c = FDCT @ x
xinv = FDCT.H @ c
np.testing.assert_allclose(x, xinv)
```

--------------------------------

### Sparse Reconstruction with Curvelets

Source: https://context7.com/pylops/curvelops/llms.txt

Example demonstrating sparse signal representation and reconstruction using curvelets, with a comparison to wavelets. Requires FDCT2D for curvelet transforms.

```python
import numpy as np
import matplotlib.pyplot as plt
from curvelops import FDCT2D
```

--------------------------------

### Load and prepare input data

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Loads a sigmoid dataset and calculates normalization parameters.

```python
inputfile = "../testdata/sigmoid.npz"

d = np.load(inputfile)
d = d["sigmoid"]
nx, nz = d.shape
dx, dz = 0.006, 0.005
x, z = np.arange(nx) * dx, np.arange(nz) * dz

vmax = 0.5 * np.max(np.abs(d))
```

--------------------------------

### Set up interactive sliders for curvelet visualization

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Single_Curvelet_Interactive.ipynb

Creates sliders for scale, wedge, x-index, and y-index to interactively control the displayed curvelet. Includes observer functions to update slider maximums based on selected scale and wedge.

```python
max_scale = DCT.nbscales
max_wedge = len(y_struct[0])
max_iy, max_ix = y_struct[0][0].shape
curr_scale = 1
curr_wedge = 1

slider_scale = IntSlider(
    min=1, max=max_scale, value=curr_scale, step=1, description="Scales"
)
slider_wedge = IntSlider(
    min=1, max=max_wedge, value=curr_wedge, step=1, description="Wedge"
)
slider_ix = IntSlider(
    min=1, max=max_ix, value=max_ix // 2 + 1, step=1, description="X Index"
)
slider_iy = IntSlider(
    min=1, max=max_iy, value=max_iy // 2 + 1, step=1, description="Y Index"
)


def handle_scale_change(change):
    global curr_scale
    curr_scale = change.new
    slider_wedge.max = len(y_struct[curr_scale - 1])
    global curr_wedge
    curr_wedge = slider_wedge.value
    A, B = y_struct[curr_scale - 1][curr_wedge - 1].shape
    slider_ix.max = B
    slider_iy.max = A


def handle_wedge_change(change):
    global curr_wedge
    curr_wedge = change.new
    A, B = y_struct[curr_scale - 1][curr_wedge - 1].shape
    slider_ix.max = B
    slider_iy.max = A


slider_scale.observe(handle_scale_change, names="value")
slider_wedge.observe(handle_wedge_change, names="value")

out = interactive_output(
    display_curvelet,
    {
        "scale": slider_scale,
        "wedge": slider_wedge,
        "ix": slider_ix,
        "iy": slider_iy,
    },
)
vbox1 = VBox([slider_scale, slider_wedge])
vbox2 = VBox([slider_ix, slider_iy])
ui = HBox([vbox1, vbox2])
display(ui, out)
```

--------------------------------

### Initialize FDCT2D with Custom Parameters

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Initializes the 2D Fast Discrete Curvelet Transform (FDCT2D) with specified parameters for number of scales, coarse angles, and whether to include all curvelets on the finest scale. Use `allcurvelets=False` unless a specific reason requires it.

```python
C2D = FDCT2D(
    logo.shape[:-1], nbscales=4, nbangles_coarse=8, allcurvelets=False
)
```

--------------------------------

### Run Development Checks

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Execute various development checks including tests, linting, type annotation checks, coverage, and documentation builds.

```console
$ make tests
$ make lint
$ make typeannot
$ make coverage
$ make doc
$ make servedoc
```

--------------------------------

### Import FDCT2D and curveshow from curvelops

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Import necessary components for using the Fast Discrete Curvelet Transform and its plotting utility.

```python
from curvelops import FDCT2D
from curvelops.plot import curveshow
```

--------------------------------

### Visualize directional variability in k-space

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Generates a grid of monochromatic signals and plots their k-space representation with annotated direction vectors.

```python
nthetas = 6
frequency = 5
thetas = np.linspace(-90, 60, nthetas)

cols = ceil(np.sqrt(nthetas))
rows = ceil(nthetas / cols)
fig, axes = plt.subplots(rows, cols, figsize=(2 * cols, 2 * rows))
for theta, ax in zip(thetas, axes.ravel()):
    img, x, z, _ = generate_monochromatic_2d(theta, freq=frequency, taper=True)
    img_k = F2D @ img
    kvec_max = find_kmax(img_k, kx, kz)

    ax.imshow(
        np.abs(img_k).T,
        vmin=0,
        cmap="turbo",
        extent=[kx[0], kx[-1], kz[-1], kz[0]],
    )
    ax.set(
        xlim=[-2 * frequency, 2 * frequency],
        ylim=[2 * frequency, -2 * frequency],
    )
    ax.axhline(0, ls=":", color="w", lw=1)
    ax.axvline(0, ls=":", color="w", lw=1)
    ax.annotate(
        "",
        xy=(kvec_max[0], kvec_max[1]),
        xytext=(0.0, 0.0),
        arrowprops=dict(edgecolor="k", facecolor="y"),
    )
for ax in axes.ravel():
    ax.axis("off")
fig.suptitle("Figure 4. 2D FFT of the monochromatic signals in Figure 2")
fig.tight_layout()
```

--------------------------------

### Import plotting utilities

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Imports necessary modules for plotting curvelet overlays.

```python
import matplotlib as mpl
from mpl_toolkits.axes_grid1 import make_axes_locatable

from curvelops.plot import (
    create_inset_axes_grid,
    create_axes_grid,
    overlay_disks,
)
```

--------------------------------

### Generate and visualize bichromatic signal

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Creates a signal with two distinct frequency components and displays them in spatial and Fourier domains.

```python
img_lo, *_ = generate_monochromatic_2d(
    -15,
    freq=0.15 * np.abs(kz).max(),
    taper=True,  # taper prevents Gibb's effect
)
img_hi, *_ = generate_monochromatic_2d(
    15, freq=0.75 * np.abs(kz).max(), taper=True
)
img = 0.4 * img_lo + 0.6 * img_hi
img_k = F2D @ img

vmax = 0.75 * np.abs(img).max()
vmax_k = 0.5 * np.abs(img_k).max()
```

```python
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
axes[0].imshow(
    img.T,
    vmin=-vmax,
    vmax=vmax,
    cmap="gray",
    extent=[x[0], x[-1], z[-1], z[0]],
)
axes[1].imshow(
    np.abs(img_k).T,
    vmin=0,
    vmax=vmax_k,
    cmap="turbo",
    extent=[kx[0], kx[-1], kz[-1], kz[0]],
)
axes[1].axhline(0, ls=":", color="w", lw=1)
axes[1].axvline(0, ls=":", color="w", lw=1)
for ax in axes.ravel():
    ax.axis("off")
fig.suptitle("Figure 5. Bichromatic signal in spatial (left) and Fourier (right) domains")
fig.tight_layout()
```

--------------------------------

### Initialize FDCT2D Transform

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Initializes the Fast Discrete Curvelet Transform (FDCT2D) for a 2D image. Configurable parameters include the number of scales, number of coarse angles, and whether to include all curvelets.

```python
C2D = FDCT2D(logo.shape[:-1], nbscales=3, nbangles_coarse=8, allcurvelets=True)
```

--------------------------------

### Build Low-pass and High-pass Filters

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Creates a frequency mask to separate signals based on Nyquist frequency thresholds with optional smoothing.

```python
def build_lowpass(kx, kz, nsmooth=None):
    if nsmooth is None:
        nsmooth = 0.01 * min(len(kx), len(kz))
        nsmooth = 2 * max(1, floor(nsmooth / 2)) + 1
    mask = np.zeros((len(kx), len(kz)))
    KX, KZ = np.meshgrid(kx, kz, indexing="ij")
    kxnyq = np.abs(kx).max()
    kznyq = np.abs(kz).max()
    inside_scale = (np.abs(KX) < 0.5 * kxnyq) & (np.abs(KZ) < 0.5 * kznyq)
    mask[inside_scale] = 1.0
    s = np.ones((nsmooth,)) / nsmooth
    mask = filtfilt(s, 1, mask, axis=0)
    mask = filtfilt(s, 1, mask, axis=1)
    return mask


lowpass = build_lowpass(kx, kz)
highpass = 1.0 - lowpass
```

--------------------------------

### Visualize Signal Separation by Wedges

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

This code visualizes the original signal, the k-space representation, and the separated signal using wedge masks. It plots the transformed image, its absolute k-space magnitude, and the corresponding wedge mask. Note that only four wedges are shown due to symmetry for real signals.

```python
fig, axes = plt.subplots(3, 5, figsize=(16, 10))
for i, _mask in enumerate(
    [lowpass + highpass, *wedges],
):
    _img_k = img_k * _mask
    _img = (F2D.H @ _img_k).real
    axes[0, i].imshow(
        _img.T,
        vmin=-vmax,
        vmax=vmax,
        cmap="gray",
        extent=[x[0], x[-1], z[-1], z[0]],
    )
    axes[0, i].text(
        0.5,
        0.95,
        f"Shape: {_img.shape}",
        horizontalalignment="center",
        verticalalignment="center",
        transform=axes[0, i].transAxes,
    )
    axes[1, i].imshow(
        np.abs(_img_k).T,
        vmin=0,
        vmax=vmax_k,
        cmap="turbo",
        extent=[kx[0], kx[-1], kz[-1], kz[0]],
    )
    kxmax = np.abs(kx).max()
    kzmax = np.abs(kz).max()
    axes[1, i].plot(
        0.5 * np.array([-kxmax, kxmax, kxmax, -kxmax, -kxmax]),
        0.5 * np.array([-kzmax, -kzmax, kzmax, kzmax, -kzmax]),
        color="w",
    )
    for s in [1, -1]:
        axes[1, i].plot(
            [s * kxmax, s * 0.5 * kxmax],
            [s * kzmax, s * 0.5 * kzmax],
            color="w",
        )
        axes[1, i].plot(
            [kxmax, 0.5 * kxmax],
            [s * kzmax, s * 0.5 * kzmax],
            color="w",
        )
        axes[1, i].plot(
            [-s * kxmax, -s * 0.5 * kxmax],
            [kzmax, 0.5 * kzmax],
            color="w",
        )
        axes[1, i].plot(
            [0, 0],
            [s * kzmax, s * 0.5 * kzmax],
            color="w",
        )
        axes[1, i].plot(
            [s * kxmax, s * 0.5 * kxmax],
            [0, 0],
            color="w",
        )
    axes[2, i].imshow(
        _mask.T,
        vmin=0,
        cmap="turbo",
        extent=[kx[0], kx[-1], kz[-1], kz[0]],
    )
for ax in axes.ravel():
    ax.axis("off")
fig.suptitle("Figure 9. Original tetrachromatic signal (leftmost column). Signal separated by \"wedges\"
in the Fourier domain (center-left to rightmost columns). Rows like Figures 6, 7 and 8.")
fig.tight_layout()
```

--------------------------------

### Load and Preprocess Image for Curvelet Transform

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Loads a PNG image and preprocesses it for curvelet transform. It swaps axes, normalizes pixel values to [0, 1], and converts to float. Ensure `matplotlib.pyplot` is imported.

```python
logo = plt.imread("../testdata/python.png", format="PNG")
logo = logo.swapaxes(0, 1)

fig, ax = plt.subplots()
ax.imshow(logo.swapaxes(0, 1))
ax.axis("off")

logo = logo.astype(float)
logo /= 255.0
```

--------------------------------

### Create a New Feature Branch

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Create a new branch for your feature development, following the convention of 'patch-some-cool-feature'.

```console
$ git checkout -b patch-some-cool-feature
```

--------------------------------

### Curvelops Module Overview

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/modules.rst

Overview of the primary modules available in the Curvelops library.

```APIDOC
## Module: curvelops.curvelops

### Description
Core functionality for the curvelops library.

## Module: curvelops.plot

### Description
Visualization utilities for plotting curvelops data.

## Module: curvelops.typing

### Description
Type definitions and annotations used throughout the library.

## Module: curvelops.utils

### Description
General utility functions and helper methods.
```

--------------------------------

### Create Inset Axes Grid

Source: https://context7.com/pylops/curvelops/llms.txt

Creates a grid of inset axes overlaid on an existing axis. Ideal for displaying detailed overlays or secondary plots within a main visualization. Allows specifying inset dimensions and projection.

```python
import numpy as np
import matplotlib.pyplot as plt
from curvelops.plot import create_inset_axes_grid

# Create main figure with image
fig, ax = plt.subplots(figsize=(8, 8))
ax.imshow(
    np.random.randn(100, 100),
    extent=[-2, 2, 2, -2],
    cmap="gray"
)

# Create 3x3 grid of polar inset axes
rows, cols = 3, 3
inset_axes = create_inset_axes_grid(
    ax, rows, cols,
    width=0.6,              # Width as fraction of column
    height=0.6,             # Height as fraction of row
    kwargs_inset_axes=dict(projection="polar")
)

# Plot in each inset
for irow in range(rows):
    for icol in range(cols):
        theta = np.linspace(0, 2*np.pi, 100)
        r = 0.5 + 0.3 * np.sin(3*theta + irow + icol)
        inset_axes[irow, icol].plot(theta, r)
        inset_axes[irow, icol].axis("off")

plt.show()
```

--------------------------------

### Import necessary libraries

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Single_Curvelet_Interactive.ipynb

Imports libraries for plotting, numerical operations, IPython display, and curvelops.

```python
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
from ipywidgets import HBox, IntSlider, VBox, interactive_output

from curvelops import FDCT2D
```

--------------------------------

### Overlay directional vectors on spatial images

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Visualizes the comparison between ground truth variability and k-space estimated variability.

```python
nthetas = 6
frequency = 5
thetas = np.linspace(-90, 60, nthetas)

cols = ceil(np.sqrt(nthetas))
rows = ceil(nthetas / cols)
fig, axes = plt.subplots(rows, cols, figsize=(2 * cols, 2 * rows))
for theta, ax in zip(thetas, axes.ravel()):
    img, x, z, vec = generate_monochromatic_2d(theta, freq=frequency)

    img_k = F2D @ img
    kvec_max = find_kmax(img_k, kx, kz, unit=True)

    im = ax.imshow(
        img.T, vmin=-1, vmax=1, cmap="gray", extent=[x[0], x[-1], z[-1], z[0]]
    )
    norm = 0.5
    offset = 0.1 * np.array([vec[1], -vec[0]])
    ax.annotate(
        "",
        xy=(-offset[0] + vec[0] * norm, -offset[1] + vec[1] * norm),
        xytext=-offset,
        arrowprops=dict(edgecolor="r", facecolor="r"),
    )
    ax.annotate(
        "",
        xy=(offset[0] + kvec_max[0] * norm, offset[1] + kvec_max[1] * norm),
        xytext=offset,
        arrowprops=dict(edgecolor="k", facecolor="y"),
    )
for ax in axes.ravel():
    ax.axis("off")
fig.suptitle("Figure 5. Monochromatic signals in Figure 2.\nRed arrow: direction of maximum variability. Yellow arrows: \ndirection of maximum variability estimated from 2D FFT")
fig.tight_layout()
```

--------------------------------

### Initialize FDCT2D operator and structure

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Single_Curvelet_Interactive.ipynb

Creates an FDCT2D operator for a given image size and initializes an empty structure to hold curvelet coefficients.

```python
nx = 300
nz = 350

# Create operator
DCT = FDCT2D((nx, nz), nbangles_coarse=8)

# Create empty structure for curvelet
y_struct = DCT.struct(np.zeros(DCT.shape[0]))
```

--------------------------------

### Build 8-Wedge Symmetric Mask in k-space

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

This function constructs a symmetric mask for 8 wedges in the k-space domain. It's used to split signals based on their angular direction. The coarsest scale (scale=0) is typically excluded from this splitting.

```python
def build_8wedges_sym(kx, kz, wedge=0, nsmooth=None):
    lopass = build_lowpass(kx, kz, nsmooth=nsmooth)
    if nsmooth is None:
        nsmooth = 0.01 * min(len(kx), len(kz))
        nsmooth = 2 * max(1, floor(nsmooth / 2)) + 1

    mask = np.zeros_like(lopass)
    KX, KZ = np.meshgrid(kx, kz, indexing="ij")
    kxnyq = np.abs(kx).max()
    kznyq = np.abs(kz).max()

    s = 1 if wedge == 0 or wedge == 1 else -1
    K1 = KX if wedge == 0 or wedge == 2 else KZ
    K2 = KZ if wedge == 0 or wedge == 2 else KX
    inside_wedge = (s * (np.abs(K1) - np.abs(K2)) > 0) & (s * KX * KZ > 0)
    mask[inside_wedge] = 1
    mask[inside_wedge] -= lopass[inside_wedge]
    s = np.ones((nsmooth,)) / nsmooth
    mask = filtfilt(s, 1, mask, axis=0)
    mask = filtfilt(s, 1, mask, axis=1)
    return mask


wedges = [build_8wedges_sym(kx, kz, wedge=wedge) for wedge in range(4)]
```

--------------------------------

### Initialize and apply FDCT2D

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Initialize the FDCT2D operator with image dimensions and curvelet parameters. The operator is then applied to an image, and the resulting coefficients are structured for easier access. Note that FDCT2D transposes image dimensions.

```python
C2D = FDCT2D(img.shape, nbscales=2, nbangles_coarse=8, allcurvelets=True)

# `FDCT2D` returns an unstructured unidimensional vector, we will reorganize the
# coefficients by structuring them as a list of lists. These in turn can be
# indexed as C[scale][wedge][zpos, xpos].
# Note that the `FDCT2D` transposes the dimensions of the original image.
# This is a remnant of the original implementation.
img_c = C2D.struct(C2D @ img)
```

--------------------------------

### Update Main Branch

Source: https://github.com/pylops/curvelops/blob/main/docssrc/source/contributing.rst

Ensure you are on the latest main branch by checking it out and pulling upstream changes.

```console
$ git checkout main
$ git pull upstream main
```

--------------------------------

### Apply 2D FFT to a monochromatic image

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Applies a 2D Fast Fourier Transform (FFT) to a generated monochromatic image using pylops.FFT2D. This transforms the image into the frequency domain (k-space) to analyze its spatial frequencies.

```python
img, x, z, vec = generate_monochromatic_2d(45)
dx, dz = x[1] - x[0], z[1] - z[0]

F2D = FFT2D(
    img.shape,
    sampling=(dx, dz),
    ifftshift_before=True,
    fftshift_after=True,
    norm="ortho",
    real=False,
)
kx = F2D.f1  # Unit is 1 / (unit of dx)
kz = F2D.f2  # Unit is 1 / (unit of dz)

img_k = F2D @ img
```

--------------------------------

### Visualize 2D monochromatic images

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Generates and displays a grid of 2D monochromatic images with varying angles. Each image shows directional variations indicated by a red arrow. This helps in understanding how image orientation relates to spatial frequencies.

```python
nthetas = 6
frequency = 5
thetas = np.linspace(-90, 60, nthetas)

cols = ceil(np.sqrt(nthetas))
rows = ceil(nthetas / cols)
fig, axes = plt.subplots(rows, cols, figsize=(2 * cols, 2 * rows))
for iax, (theta, ax) in enumerate(zip(thetas, axes.ravel())):
    img, x, z, vec = generate_monochromatic_2d(theta, freq=frequency)
    ax.imshow(
        img.T, vmin=-1, vmax=1, cmap="gray", extent=[x[0], x[-1], z[-1], z[0]]
    )
    ax.set_title(f"{theta:.0f}°")
    ax.annotate(
        f"",
        xy=(vec[0] * 0.5, vec[1] * 0.5),
        xytext=(0, 0),
        arrowprops=dict(edgecolor="r", facecolor="r"),
    )
    if iax == 0 or iax == cols:
        ax.set_ylabel("z\n" + r"$\downarrow$", rotation=0)
    if iax >= cols:
        ax.set_xlabel(r"x $\rightarrow$")
for ax in axes.ravel():
    ax.xaxis.set_ticks([])
    ax.yaxis.set_ticks([])
fig.suptitle("Figure 2. 2D monochromatic images, with direction\nof change in each image denoted by the red arrow.")
fig.tight_layout()
```

--------------------------------

### Apply FDCT2D and Structure Coefficients

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Applies the FDCT2D to each color channel (R, G, B) of the input image and structures the resulting coefficients. This is typically done after initializing the `FDCT2D` object.

```python
logo_r = C2D.struct(C2D @ logo[..., 0])
logo_g = C2D.struct(C2D @ logo[..., 1])
logo_b = C2D.struct(C2D @ logo[..., 2])
```

--------------------------------

### Display Curvelet Coefficients as Images

Source: https://context7.com/pylops/curvelops/llms.txt

Visualizes curvelet coefficients as images, with one figure per scale showing all angular wedges. Coefficients are normalized for better visualization. Supports spatial or frequency domain display.

```python
import numpy as np
import matplotlib.pyplot as plt
from curvelops import FDCT2D
from curvelops.plot import curveshow
from curvelops.utils import apply_along_wedges

# Create data and compute curvelet transform
data = np.random.randn(101, 101).astype(complex)
C = FDCT2D(data.shape, nbscales=3, nbangles_coarse=8)
y_struct = C.struct(C @ data)

# Normalize coefficients for visualization
y_norm = apply_along_wedges(
    y_struct,
    lambda w, *_: w / (np.abs(w).max() + 1e-10)
)

# Display curvelet coefficients (one figure per scale)
figs_axes = curveshow(
    y_norm,
    k_space=False,          # Show spatial domain (True for frequency)
    basesize=3,             # Figure size multiplier
    showaxis=False,         # Hide axis labels
    real=True,              # Show real part (False for imaginary)
    kwargs_imshow=dict(
        aspect="auto",
        vmin=-0.5,
        vmax=0.5,
        cmap="RdBu"
    )
)
plt.show()
```

--------------------------------

### Resample Coarsest Scale

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Demonstrates resampling the low-frequency scale to simulate curvelet coefficient behavior.

```python
fig, axes = plt.subplots(3, 3, figsize=(10, 10))
for i, _mask in enumerate(
    [lowpass + highpass, lowpass, highpass],
):
    _img_k = _mask * img_k
    if i == 1:  # Coarsest scale
        kxmax = 0.5 * np.abs(kx).max()
        kzmax = 0.5 * np.abs(kz).max()
        _img_k_lp = _img_k[np.abs(kx) < kxmax, :][:, np.abs(kx) < kzmax]
        F2D_lp = FFT2D(
            _img_k_lp.shape,
            sampling=(2 * dx, 2 * dz),
            ifftshift_before=True,
            fftshift_after=True,
            norm="ortho",
            real=False,
        )
        _img = (F2D_lp.H @ _img_k_lp).real
```

--------------------------------

### Perform Curvelet Transform and Sparse Reconstruction

Source: https://context7.com/pylops/curvelops/llms.txt

Initializes a 2D curvelet transform, performs a forward transform, and reconstructs the signal using a subset of the strongest coefficients.

```python
# Create test signal (e.g., image with curved features)
nx, nt = 128, 128
x = np.linspace(-1, 1, nx)
t = np.linspace(-1, 1, nt)
xm, tm = np.meshgrid(x, t, indexing="ij")
data = np.sin(5 * np.pi * (xm + 0.5 * tm**2)).astype(complex)

# Create curvelet transform
C = FDCT2D(data.shape, nbscales=4, nbangles_coarse=16, allcurvelets=True)
print(f"Transform: {data.shape} -> {C.shape[0]} coefficients")

# Forward transform
coeffs = C @ data

# Sparse reconstruction: keep only top N% of coefficients
def sparse_reconstruct(data, C, keep_fraction=0.1):
    """Reconstruct using only strongest coefficients."""
    y = C @ data
    threshold_idx = int(len(y) * keep_fraction)

    # Sort by magnitude and keep strongest
    sorted_idx = np.argsort(-np.abs(y))
    y_sparse = np.zeros_like(y)
    y_sparse[sorted_idx[:threshold_idx]] = y[sorted_idx[:threshold_idx]]

    # Inverse transform
    return C.H @ y_sparse

# Reconstruct with different sparsity levels
for keep_pct in [0.05, 0.10, 0.20]:
    data_rec = sparse_reconstruct(data, C, keep_pct)
    error = np.linalg.norm(data - data_rec.reshape(data.shape))
    print(f"Keep {keep_pct*100:.0f}%: L2 error = {error:.4f}")

# Visualize
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
axes[0].imshow(data.real, cmap="RdBu")
axes[0].set_title("Original")

data_10pct = sparse_reconstruct(data, C, 0.10).reshape(data.shape)
axes[1].imshow(data_10pct.real, cmap="RdBu")
axes[1].set_title("10% coefficients")

axes[2].imshow((data - data_10pct).real, cmap="RdBu")
axes[2].set_title("Reconstruction error")

plt.tight_layout()
plt.show()
```

--------------------------------

### Visualize Frequency Scales

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Plots the original, low-passed, and high-passed signals in both spatial and Fourier domains.

```python
fig, axes = plt.subplots(3, 3, figsize=(10, 10))
for i, _mask in enumerate(
    [lowpass + highpass, lowpass, highpass],
):
    _img_k = _mask * img_k
    _img = (F2D @ _img_k).real
    axes[0, i].imshow(
        _img.T,
        vmin=-vmax,
        vmax=vmax,
        cmap="gray",
        extent=[x[0], x[-1], z[-1], z[0]],
    )
    axes[1, i].imshow(
        np.abs(_img_k).T,
        vmin=0,
        vmax=vmax_k,
        cmap="turbo",
        extent=[kx[0], kx[-1], kz[-1], kz[0]],
    )
    kxnyq = 0.5 * np.abs(kx).max()
    kznyq = 0.5 * np.abs(kz).max()
    axes[1, i].plot(
        [-kxnyq, kxnyq, kxnyq, -kxnyq, -kxnyq],
        [-kznyq, -kznyq, kznyq, kznyq, -kznyq],
        color="w",
    )
    axes[2, i].imshow(
        _mask.T,
        vmin=0,
        cmap="turbo",
        extent=[kx[0], kx[-1], kz[-1], kz[0]],
    )
for ax in axes.ravel():
    ax.axis("off")
fig.suptitle("Figure 6. Original bichromatic signal (left column). Low-passed signal (center column).\nHigh-passed signal (right column). Spatial domain (top row). Fourier domain (middle row).\nMasks in Fourier domain which generate the signals in the top row (bottom row).\nDeep purple denotes zero value, crimson corresponds to unity.")
fig.tight_layout()
```

--------------------------------

### Compute Spatially Localized Energy within a Wedge

Source: https://context7.com/pylops/curvelops/llms.txt

Splits a wedge into a grid of subregions and computes the energy of each. This enables spatially-localized analysis of curvelet coefficients.

```python
import numpy as np
from curvelops.utils import energy_split, array_split_nd

# Create a wedge with varying energy
wedge = np.zeros((20, 30), dtype=complex)
wedge[:10, :15] = 1 + 1j      # High energy in top-left
wedge[10:, 15:] = 0.5 + 0.5j  # Lower energy in bottom-right

# Split into 2x3 grid and compute energy of each region
energy_grid = energy_split(wedge, rows=2, cols=3)
print(f"Energy grid shape: {energy_grid.shape}")
print(f"Energy values:\n{energy_grid}")
```

--------------------------------

### Plot curvelet coefficient disks

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Visualizes the strength of curvelet coefficients as disks overlaid on the original signal.

```python
rows, cols = 4, 5
fig, ax = plt.subplots(figsize=(8, 10))
ax.imshow(
    d.T, vmin=-vmax, vmax=vmax, extent=(x[0], x[-1], z[-1], z[0]), cmap="gray"
)
axesin = create_inset_axes_grid(
    ax, rows, cols, width=0.4, kwargs_inset_axes=dict(projection="polar")
)
overlay_disks(d_c, axesin, linewidth=0.0, cmap="turbo")
ax.set(xlabel="Position [km]", ylabel="Depth [km]")

divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
mpl.colorbar.ColorbarBase(
    cax, cmap="turbo", norm=mpl.colors.Normalize(vmin=0, vmax=1)
)
```

--------------------------------

### Overlay Curvelet Energy Disks

Source: https://context7.com/pylops/curvelops/llms.txt

Visualizes the energy distribution of curvelet coefficients as polar disks. Maps energy to colormap and overlays on a polar axis. Useful for analyzing frequency and scale content.

```python
import matplotlib.pyplot as plt
import numpy as np
from curvelops import FDCT2D
from curvelops.plot import create_axes_grid, overlay_disks
from curvelops.utils import apply_along_wedges

# Create signal and compute curvelets
x = np.random.randn(100, 100).astype(complex)
C = FDCT2D(x.shape, nbscales=4, nbangles_coarse=8)
y_struct = C.struct(C @ x)

# Create figure with polar axis
fig, axes = create_axes_grid(
    1, 1,
    kwargs_subplots=dict(projection="polar"),
    kwargs_figure=dict(figsize=(8, 8)),
)

# Overlay curvelet energy disks
overlay_disks(
    y_struct,
    axes,
    linewidth=0.5,          # Line separating scales
    linecolor="r",          # Color of scale separator
    map_cmap=True,          # Map energy to colormap
    cmap="gray_r",          # Colormap for energy
    alpha=1.0,              # Transparency
    pclip=0.5,              # Clip maximum amplitude
    normalize="all",        # Normalize across all wedges ("scale" for per-scale)
    annotate=True           # Show scale,wedge indices
)

plt.title("Curvelet Energy Distribution")
plt.show()
```

--------------------------------

### Visualize Curvelet Coefficients

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Visualizes curvelet coefficients using the `curveshow` function. It iterates through scales and wedges, displaying their shapes. Ensure `curveshow` is imported and `img_c` is properly formatted.

```python
figs_axes = curveshow(
    img_c,
    basesize=3,
    kwargs_imshow=dict(vmin=-1, vmax=1, extent=[x[0], x[-1], z[-1], z[0]]),
)
for c_scale, (fig, axes) in zip(img_c, figs_axes):
    for iwedge, (c_wedge, ax) in enumerate(
        zip(c_scale, np.atleast_1d(axes).ravel())
    ):
        ax.text(
            0.5,
            0.95,
            f"Shape: {c_wedge.shape}",
            horizontalalignment="center",
            verticalalignment="center",
            transform=ax.transAxes,
        )
```

--------------------------------

### Create a 2D Grid of Matplotlib Axes

Source: https://context7.com/pylops/curvelops/llms.txt

Creates a 2D grid of matplotlib axes with customizable figure and gridspec parameters. Useful for arranging multiple plots in a structured layout.

```python
import numpy as np
import matplotlib.pyplot as plt
from curvelops.plot import create_axes_grid
```

--------------------------------

### Generate 2D monochromatic images

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Creates a 2D image with sinusoidal variations along a specified direction. Optionally applies a 2D taper to the image. Useful for visualizing directional frequencies.

```python
def generate_monochromatic_2d(
    theta_deg: float = 0,
    nx: int = 101,
    nz: int = 101,
    freq: float = 10,
    taper: bool = False,
):
    x = np.linspace(-1, 1, nx)
    z = np.linspace(-1, 1, nz)
    dx, dz = x[1] - x[0], z[1] - z[0]
    xm, zm = np.meshgrid(x, z, indexing="ij")
    theta = np.deg2rad(theta_deg)
    vec = np.array([np.sin(theta), np.cos(theta)])
    img = np.cos(2 * np.pi * freq * (vec[0] * xm + vec[1] * zm))

    if taper:
        img *= taper2d(*img.shape, [nx // 5 + 1, nz // 5 + 1])
    return img, x, z, vec
```

--------------------------------

### Apply Function to Each Wedge in Curvelet Transform

Source: https://context7.com/pylops/curvelops/llms.txt

Applies a custom function to each curvelet wedge. Useful for normalization or analysis. Ensure the provided function signature matches `fun(wedge, wedge_idx, scale_idx, num_wedges, num_scales)`.

```python
import numpy as np
from curvelops import FDCT2D
from curvelops.utils import apply_along_wedges

# Create transform
data = np.random.randn(64, 64).astype(complex)
C = FDCT2D(data.shape, nbscales=3, nbangles_coarse=8)
y_struct = C.struct(C @ data)

# Get shape of each wedge
# Function signature: fun(wedge, wedge_idx, scale_idx, num_wedges, num_scales)
shapes = apply_along_wedges(y_struct, lambda w, *_: w.shape)
print("Wedge shapes by scale:")
for iscale, scale_shapes in enumerate(shapes):
    print(f"  Scale {iscale}: {scale_shapes}")

# Normalize each wedge by its maximum absolute value
y_normalized = apply_along_wedges(
    y_struct,
    lambda w, iw, iscale, nw, ns: w / (np.abs(w).max() + 1e-10)
)

# Compute energy per wedge
from curvelops.utils import energy
energies = apply_along_wedges(y_struct, lambda w, *_: energy(w))
print(f"Energies per scale: {[sum(e) for e in energies]}")
```

--------------------------------

### Create Grid of Axes

Source: https://context7.com/pylops/curvelops/llms.txt

Generates a grid of axes for plotting. Useful for organizing multiple plots in a structured layout. Supports customization of figure and grid parameters.

```python
import numpy as np
import matplotlib.pyplot as plt
from curvelops.plot import create_axes_grid

# Create 2x3 grid of axes
rows, cols = 2, 3
fig, axs = create_axes_grid(
    rows, cols,
    kwargs_figure=dict(figsize=(12, 8)),
    kwargs_gridspec=dict(wspace=0.3, hspace=0.3),
    kwargs_subplots=dict()  # e.g., projection="polar"
)

# Plot different data in each axis
for irow in range(rows):
    for icol in range(cols):
        freq = 2 + irow + icol**2
        x = np.linspace(0, 2*np.pi, 100)
        axs[irow, icol].plot(x, np.cos(freq * x))
        axs[irow, icol].set_title(f"Row {irow}, Col {icol}")

plt.tight_layout()
plt.show()
```

--------------------------------

### Compute Wedge Energy using Root-Mean-Square

Source: https://context7.com/pylops/curvelops/llms.txt

Computes the root-mean-square energy of a curvelet wedge. This is useful for analyzing coefficient distribution and importance.

```python
import numpy as np
from curvelops.utils import energy

# Energy is sqrt(mean(|x|^2))
wedge = np.array([[1+1j, 2+2j], [3+3j, 4+4j]])
e = energy(wedge)
print(f"Wedge energy: {e:.4f}")

# Manual calculation for verification
manual_energy = np.sqrt((wedge.real**2 + wedge.imag**2).sum() / wedge.size)
np.testing.assert_allclose(e, manual_energy)
```

--------------------------------

### Visualize 2D FFT spectrum

Source: https://github.com/pylops/curvelops/blob/main/notebooks/Desmystifying_Curvelets.ipynb

Displays the magnitude of the 2D FFT spectrum of a monochromatic image. The plot shows peaks corresponding to the dominant spatial frequencies and their directions in k-space. Helps in identifying the signal's frequency content.

```python
fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(
    np.abs(img_k).T,
    vmin=0,
    cmap="turbo",
    extent=[kx[0], kx[-1], kz[-1], kz[0]],
)
ax.set(
    xlim=[-2 * frequency, 2 * frequency], ylim=[2 * frequency, -2 * frequency],
    xlabel=r"$k_x \rightarrow$"
)
ax.set_ylabel(r"$k_z$" + "\n" + r"$\downarrow$", rotation=0)
ax.axhline(0, ls=":", color="w", lw=1)
ax.axvline(0, ls=":", color="w", lw=1)
fig.suptitle("Figure 3. 2D FFT of the monochromatic\nsignal which varies most at the 45° direction.")
```