API Reference

ScatteringTransforms.jl — Native Julia implementation of wavelet scattering transforms

Surfaces: 1D signals, 2D and 3D gridded fields, scattered planar points (NUFFT), and the sphere — both on a structured grid and at scattered points. Each has a monogenic variant; the gridded ones also have a localized (Mallat) field output and a multi-resolution second order. The spherical and scattered paths run dependency-free by default, with fast paths supplied by extensions.

Quick Start

Names live in submodules and are reached by their full path; alias the package for brevity.

using ScatteringTransforms: ScatteringTransforms as ST

# 1D scattering (N = signal length, J = number of octaves)
signal = randn(1024)
st = ST.Scattering1D.ScatteringTransform1D(1024, 8; Q=1, max_order=2)
coeffs = st(signal)

# 2D planar scattering (N = image size, J = scales, L = orientations)
image = randn(256, 256)
st2d = ST.Scattering2D.ScatteringTransform2D((256, 256), 4; L=8, max_order=2)
coeffs2d = st2d(image)

Implementation Notes

• FFT-based convolutions for O(N log N) performance
• Frequency-domain Morlet filter banks
• Breadth-first CSR path list, walked grouped by first-order wavelet
• Modular extensions for spherical (NUFSHT) and GPU support

References

• Mallat (2012): Group invariant scattering. Comm. Pure Appl. Math.
• Bruna & Mallat (2013): Invariant Scattering Convolution Networks. IEEE PAMI.
• Cheng & Ménard (2021): How to quantify fields or textures? A guide to the scattering transform.

ScatteringTransform1D{T,V,M,P,Tree,FB,G}

1D scattering transform: a filter bank, the admissible path tree, a spectral plan, and the workspace the cascade runs in. Every array field is a type parameter, so the same struct holds CPU, GPU or static storage.

Fields

• filter_bank: pre-computed 1D filter bank
• tree: admissible scattering paths
• groups: (j1, children) from tree, longest-first — the order the cascade walks
• max_order: maximum scattering order (1 or 2)
• plan: spectral transform plan (in-core direct sum by default; FFTW fast path if loaded)
• buffer_input: complex buffer for real→complex promotion, and multiply scratch
• buffer_signal_fft: the input spectrum, read-only for the whole cascade
• buffer_conv: inverse-transform output
• buffer_mod: real modulus buffer for the localized-field path
• buffer_u1, buffer_u1_fft: the current first-order modulus and its spectrum, reused across that wavelet's children — one pair, not one per wavelet

ScatteringTransform2D{T,M,R}

2D planar scattering transform with oriented wavelets and pre-allocated workspace.

Type Parameters

• T: Real element type (Float32, Float64, ...)
• M: Complex matrix type for buffers (Matrix{Complex{T}}, CuMatrix{Complex{T}}, ...)
• R: Real matrix type for modulus buffers (Matrix{T}, ...)

Fields

• filter_bank: Pre-computed 2D filter bank
• max_order::Int: Maximum scattering order (1 or 2)
• plan: spectral transform plan (in-core direct sum by default; FFTW fast path if loaded)
• buffer_input: Complex matrix for real→complex promotion (zero alloc)
• buffer_signal_fft: Preserved copy of signal FFT (buffer_conv gets overwritten)
• buffer_conv: Complex matrix for IFFT output
• buffer_mod: Real matrix for modulus output
• buffer_u1, buffer_u1_fft: the current first-order modulus and its spectrum, reused across that wavelet's children — one pair, not one per wavelet

ScatteringTransform3D([T=Float64,] N, J; n_orient=6, max_order=2, spectral=AutoSpectralBackend())

Build a 3D volumetric scattering transform for N = (Nz, Ny, Nx) volumes over J scales and n_orient sphere directions. The element type is positional, as for zeros(T, …); omit it for Float64.

(st::ScatteringTransform3D)(volume) -> ScatteringCoefficients2D

Apply the 3D scattering transform. (Coefficients use the scales×orientations container.)

scattering_transform!(coeffs, st, signal)

In-place scattering transform. Fills pre-allocated S1/S2, returns the coefficients with S0 updated. Allocation-free; only allocates a new wrapper struct when S0 is a scalar (immutable).

scattering_transform!(coeffs, backend, st, signal)

Transform one signal on an explicit execution backend. SerialBackend runs the cascade in this task; ThreadedBackend (OhMyThreads extension) spreads the first-order wavelet groups across tasks, which is the only parallel axis available when there is a single field rather than a batch. The input transform is done once up front, so only the group loop is parallel.

scattering_transform2d!(coeffs, st, image)

In-place 2D scattering transform. Zero allocations for S1/S2 (buffers reused).

scattering_transform2d!(coeffs, backend, st, image)

Transform one image on an explicit execution backend. SerialBackend runs the cascade in this task; ThreadedBackend (OhMyThreads extension) spreads the first-order wavelet groups across tasks, which is the only parallel axis available for a single image.

scattering_transform3d!(coeffs, st, volume) -> coeffs

In-place 3D volumetric scattering transform; fills coeffs (a scales×orientations container) and returns it with S0 updated.

scattering_transform3d!(coeffs, backend, st, volume)

Transform one volume on an explicit execution backend — see the 2D counterpart.

scattering(st, x) -> ScatteringCoefficients

Non-mutating, allocation-tolerant, element-type-generic scattering transform — the autodiff-friendly counterpart of the in-place callable st(x). It composes the non-mutating Plans.forward_transform/Plans.inverse_transform with broadcast modulus and mean (no preallocated workspace, no mul!), so gradients flow through it via DifferentiationInterface (Mooncake/Zygote/Enzyme) and it accepts Dual/Float32 inputs. It returns the same coefficient container as st(x) and matches it numerically. Methods are added for the 1D/2D/3D transforms in their respective submodules.

Use st(x) (mutating, zero-alloc) for production forward passes; use scattering(st, x) when you need to differentiate the forward map (e.g. gradient-descent synthesis).

scattered_planar_scattering(x, y, ms, J; L=8, max_order=2, T=Float64,
                            spectral=SpectralBackends.AutoSpectralBackend(), period=nothing,
                            solve=false, weights=nothing, eps=nothing, maxiter=100, rtol=1e-8)

Build a 2D planar scattering transform for a scalar field sampled at scattered points (x, y), using the same oriented Morlet wavelet bank as the gridded [ScatteringTransform2D] but computing the wavelet convolutions on a uniform Fourier mode grid of size ms = (m1, m2) via a nonuniform DFT: analysis maps the scattered points to the mode grid, the wavelet multiply happens there, and synthesis evaluates the filtered field back at the points. Apply it to a length-M vector of samples.

spectral selects the transform: SpectralBackends.DirectSumSpectralBackend is the in-core, dependency-free exact NUDFT (always available, O(M·prod(ms))); Plans.FINUFFTBackend and Plans.NonuniformFFTsBackend select a specific fast library (using FINUFFT / using NonuniformFFTs); SpectralBackends.NUFFTSpectralBackend takes whichever of those is loaded; and SpectralBackends.AutoSpectralBackend (the default) picks a fast library if one is loaded, else the direct sum. period is the physical domain size per axis (the Fourier period); it defaults so a uniform 0:m-1 grid reproduces the gridded FFT transform exactly. solve=false uses the fast adjoint (type-1) — exact for adequately-sampled band-limited fields, approximate on gappy/irregular data; solve=true uses a conjugate-gradient least-squares inversion for the true band-limited coefficients (slower, needed for irregular sampling). weights (length M, summing to 1) sets the quadrature for the spatial mean; the default is the uniform sample mean. eps is the FINUFFT tolerance (ignored by the exact direct sum).

spherical_scattering(pts_theta, pts_phi, lmax, J; max_order=2,
                     spectral=SpectralBackends.AutoSpectralBackend(), rtol=1e-8, maxiter=500)

Build a spherical scattering transform for a scalar field at scattered points (θ, φ) on S², using smooth difference-of-Gaussians band-pass wavelets. spectral selects the spherical-harmonic transform: SpectralBackends.DirectSumSpectralBackend is the in-core, dependency-free exact least-squares transform (always available, O(M·(lmax+1)²)); SpectralBackends.NUFSHTSpectralBackend uses the NUFSHT fast path (needs using NUFSHT); SpectralBackends.AutoSpectralBackend (the default) picks NUFSHT if its extension is loaded, else the direct transform. Accurate analysis needs the sampling to resolve the band limit, i.e. roughly M ≳ (lmax+1)² well-distributed points.

structured_spherical_scattering(lmax, J; max_order=2,
                                spectral=SpectralBackends.AutoSpectralBackend(), T=Float64,
                                rtol=1e-8, maxiter=500)

Build a spherical scattering transform for a scalar field sampled on the structured equiangular grid (Nθ = lmax+1 colatitudes, Nφ = 2lmax+1 longitudes), using the same smooth difference-of-Gaussians band-pass wavelets as spherical_scattering. Apply it to a (Nθ, Nφ) grid of samples; obtain the grid points with structured_sphere_points. spectral selects the transform: SpectralBackends.DirectSumSpectralBackend (in-core, dependency-free) by default, or SpectralBackends.FSHTSpectralBackend (the fast exact SHT, needs using FastSphericalHarmonics); SpectralBackends.AutoSpectralBackend picks the fast path if its extension is loaded, else the direct SHT.

structured_sphere_points(lmax) -> (Θ, Φ)

Colatitudes Θ (length lmax+1) and longitudes Φ (length 2lmax+1) of the equiangular structured grid used by structured_spherical_scattering; sample a field as [f(θ, φ) for θ in Θ, φ in Φ]. In-core (no dependency); matches FastSphericalHarmonics.sph_points.

ReconstructionWorkspace(st)

Scratch for the linear wavelet inverse and for phase retrieval: the complex wavelet coefficient fields (nw of them — that is the representation's own size), the low-pass field, a spectrum accumulator, and the conjugated filters.

The conjugates are precomputed because iwavelet! would otherwise rebuild conj.(ψ_λ) for every wavelet on every iteration, and reconstruct_phase runs iters of them.

wavelet_transform(st, x) -> (; wavelet, lowpass)
wavelet_transform!(ws, st, x) -> (; wavelet, lowpass)

The linear (pre-modulus) wavelet layer underlying the scattering transform: the complex wavelet coefficient fields wavelet[λ] = x ⋆ ψ_λ (the continuous wavelet transform) and the low-pass field lowpass = x ⋆ φ. Together these are an exact, invertible representation of x (see iwavelet!); taking |wavelet[λ]| and averaging is the first scattering layer.

iwavelet(st, wavelet, lowpass) -> x
iwavelet(st, wt) -> x
iwavelet!(ws, st, wavelet, lowpass) -> x

Exact inverse of wavelet_transform! via the tight-frame conjugate-filter sum x̂ = Σ_λ ψ̂_λ^* · (x̂·ψ̂_λ) + φ̂^* · (x̂·φ̂). Returns the reconstructed real field. With the tight-frame bank (Σ|ψ̂_λ|²+|φ̂|² ≡ 1) this satisfies iwavelet(st, wavelet_transform(st, x)...) ≈ x to machine precision.

The accumulation is in place: the spectrum sum is built in one buffer rather than rebuilt per wavelet, which is what makes iters rounds of phase retrieval affordable.

reconstruct_phase(st, moduli; iters=200, init=nothing, seed_lowpass=nothing) -> x

Phase retrieval from the first-order moduli moduli[λ] = |x ⋆ ψ_λ| (the real fields the scattering transform averages), via Gerchberg–Saxton alternating projections:

1. take the linear wavelet transform of the current estimate;
2. re-impose the target magnitudes on each band-pass channel (keeping the recovered phase);
3. reconstruct with the exact frame inverse iwavelet!; repeat.

The low-pass channel carries no magnitude target, so it is taken from the current estimate each iteration (or from seed_lowpass if supplied). The reconstruction is determined only up to a global sign/phase. init (defaults to a random field matched in energy to the moduli) seeds the estimate; pass one for reproducibility. Returns the real reconstructed field.

One workspace is built up front and reused across all iters, so the loop allocates nothing.

synthesize(st, target; backend, init=nothing, iters=500, lr=0.05, loss=scattering_loss) -> (; field, losses)

Reconstruct a field whose scattering coefficients match target by gradient descent (Bruna & Mallat microcanonical synthesis): starting from init (random by default), minimize loss(scattering(st, x̂), target) with Adam. target may be a precomputed coefficient container or a field (its coefficients are taken first). The gradient is obtained through DifferentiationInterface, so backend is any ADTypes backend (e.g. AutoMooncake()); the synthesized result is a sample with matching statistics, not the exact original (the modulus discards local phase). Requires using DifferentiationInterface and an AD backend package.

scattering_loss(c, target) -> Real

Default synthesize objective: the normalized squared error between the coefficient container c and the target, summed over the first- and (when present) second-order coefficients, ‖S₁(c)−S₁(t)‖² + ‖S₂(c)−S₂(t)‖² divided by the target energy. Differentiable in c, so it composes with scattering(st, ·) under autodiff.

MonogenicScattering{D,T,FB,Tree,P,V}

Monogenic scattering transform on a D-dimensional grid: an isotropic MonogenicFilterBank, the scattering path tree (strictly-increasing scale), a spectral plan, and reusable workspace buffers.

MonogenicFilterBank{D,T,A,W,R,MV}

Isotropic band-pass wavelets wavelets (one per scale/sub-octave), the D scale-free Riesz multipliers riesz, and the complementary low-pass averaging, forming a tight frame. Every container is a type parameter.

build_monogenic_bank([T=Float64,] dims::NTuple{D,Int}, J; Q=1) -> MonogenicFilterBank

Build a D-dimensional isotropic Morlet-style monogenic filter bank: J octaves × Q sub-octaves of radial band-pass wavelets (center frequency ξ_j = ξ₀·2^{-j/Q}, widths from the Lostanlen/Kymatio rule, reusing Filters.Morlet1D for the radial profile), the Riesz multipliers, and the tight-frame complementary low-pass.

riesz_multipliers(dims, ::Type{T}=Float64) -> NTuple{D, Array{Complex{T},D}}

The D Riesz-transform frequency multipliers R_d(k) = -i k_d/|k| (with R_d(0)=0) over a grid of size dims. Scale-free (a ratio of frequencies), so one set serves every wavelet scale. They satisfy Σ_d |R_d(k)|² = 1 off the DC bin.

monogenic_amplitude(m0, riesz_components) -> A

Monogenic amplitude A = √(m₀² + Σ_d m_d²) from the band-pass field m0 and the tuple/vector of Riesz component fields. Broadcast, so it is CPU/GPU/autodiff-generic.

monogenic_components(st, x, j) -> (; bandpass, riesz, amplitude, phase)

The monogenic decomposition of x band-passed by isotropic wavelet j (1-based): the band-pass field bandpass = x ⋆ ψ_j, the D Riesz component fields riesz, the monogenic amplitude √(bandpass² + Σ|riesz|²), and the local monogenic phase = atan(‖riesz‖, bandpass). The Riesz vector's direction gives the local orientation (e.g. atan(riesz[2], riesz[1]) in 2D).

spherical_monogenic_scattering(pts_theta, pts_phi, lmax, J; max_order=2)

Build a monogenic spherical scattering transform for a scalar field at scattered points (θ, φ) on S². The nonlinearity is the spherical monogenic amplitude A_j = √(U⁰_j² + |U^R_j|²), where U⁰_j is the difference-of-Gaussians band-pass and U^R_j is the spin-1 Riesz field R = ð∘(−Δ_S)^{-1/2}. The Riesz energy |U^R_j|² = |∇_S g_j|² (with g_j = (−Δ_S)^{-1/2} U⁰_j) is evaluated using only spin-0 spherical-harmonic transforms via the identity |∇_S g|² = ½ Δ_S(g²) − g Δ_S g, so no spin-weighted synthesis is required. spectral selects the spherical-harmonic transform as in spherical_scattering (dependency-free direct SH transform by default, NUFSHT fast path when loaded).

spherical_monogenic_components(st, field, j) -> (; bandpass, riesz, amplitude, phase, orientation)

Pointwise spherical monogenic decomposition of field band-passed at scale j — the S² analogue of the planar monogenic_components. Returns the band-pass field bandpass = U⁰_j, the spin-1 Riesz tangent vector riesz = (u_θ, u_φ) (U^R_j = ð∘(−Δ_S)^{-1/2} U⁰_j), the monogenic amplitude √(U⁰² + ‖U^R‖²), the phase = atan(‖U^R‖, U⁰), and the local orientation = atan(u_φ, u_θ) of the Riesz vector. st is a spherical_monogenic_scattering transform. With the in-core direct SH plan (dependency-free) the Riesz field is synthesized as the surface gradient of g = (−Δ_S)^{-1/2}U⁰; with a NUFSHT-backed plan it uses spin-weighted synthesis — the two agree to solver accuracy.

structured_spherical_monogenic_scattering(lmax, J; max_order=2,
                                          spectral=SpectralBackends.AutoSpectralBackend(),
                                          T=Float64, rtol=1e-8, maxiter=500)

Structured-grid counterpart of spherical_monogenic_scattering. spectral selects the SH transform as in structured_spherical_scattering (dependency-free direct SHT by default, fast exact SHT with using FastSphericalHarmonics).

scattering_field(st, x; subsample) -> ScatteringField{1,2}D

Localized (Mallat) scattering transform: returns the per-path low-passed, subsampled fields. Methods are added by the per-dimension transform modules.

scattering_field!(field, st, x) -> field

In-place localized scattering transform into a pre-allocated ScatteringField.

ScatteringField1D{T,A,Tree}

Localized 1D scattering field. data is (M, npaths); column p is the localized field of path p at the subsampled resolution M = N ÷ s. data/the integer fields stay parametric.

ScatteringField2D{T,A,Tree}

Localized 2D scattering field. data is (My, Mx, npaths); slice [:, :, p] is the localized field of path p at subsampled resolution (My, Mx) = (Ny, Nx) ÷ s.

path_field(sf, p) -> view

Non-allocating view of path p's localized field (a vector for 1D, a matrix for 2D).

subsample_factor(sf) -> Int

The decimation factor s applied to produce the field resolution.

ScatteringCoefficients1D{T,V,M,S0}

Immutable container for 1D scattering coefficients. S0 can be scalar T (return new struct) or mutable container (update in place). Uses multiple dispatch for optimal S0 handling.

Type Parameters

• T: Element type
• V: 1D array type
• M: 2D array type
• S0: S0 storage type (T for scalar, AbstractVector{T} for mutable)

ScatteringCoefficients2D{T,V,M,S0}

Immutable container for 2D planar scattering coefficients. S0 can be scalar T or mutable container - dispatch handles both optimally.

zeroth_order(c) -> Real

Extract the zeroth-order (average / S0) scattering coefficient.

first_order(c) -> AbstractVector

Extract the first-order (S1) scattering coefficients.

second_order(c) -> AbstractMatrix

Extract the second-order (S2) scattering coefficients.

flatten1d(coeffs::ScatteringCoefficients1D{T}) -> Vector{T}

Flatten to vector: [S0; S1; vec(S2 upper triangular)]. Only includes unique S2 elements where j2 > j1 (saves ~50% space).

flatten2d(coeffs::ScatteringCoefficients2D{T}) -> Vector{T}

Flatten to vector: [S0; S1; vec(S2 upper triangular)].

flatten1d!(out, c) -> out

Zero-allocation flatten into a pre-allocated vector of length flatten_length(c).

flatten2d!(out, c) -> out

Zero-allocation flatten into a pre-allocated vector of length flatten_length(c).

flatten_length(c) -> Int

Length of the flattened coefficient vector [S0; S1; vec(S2 upper triangle)].

flat_length(n) -> Int
flat_row_s0() -> Int
flat_row_s1(j, n) -> Int
flat_row_s2(j1, j2, n) -> Int

Row layout of the flattened coefficient vector [S0; S1; vec(S2 upper triangle)] for n wavelets, as produced by flatten1d!/flatten2d!.

The batched paths write straight into flattened columns rather than filling a coefficient container first, so they need the layout as arithmetic. Defining it here keeps the one definition that flatten*! also walks — a second copy elsewhere would silently drift.

update_S0(c, val)

Update the zeroth-order coefficient storage with val and return the coefficients.

compute_shape_sparsity(S1, S2, meta) -> (; sparsity, shape)

Reduced second-order descriptors (in the spirit of the reduced wavelet scattering transform, Allys et al. 2019; Cheng & Ménard 2021), as J × J matrices over scale pairs (j1, j2) with j2 > j1:

• sparsity (s₂₁): the orientation-averaged ratio ⟨S₂ / S₁⟩ — how much energy cascades from scale j1 to the coarser scale j2 (large for sparse/intermittent fields).
• shape (s₂₂): the anisotropy of the cascade — the normalized second angular harmonic ⟨S₂ · cos(2 Δθ)⟩ / ⟨S₂⟩ over orientation pairs, where Δθ = θ₂ − θ₁. It is ≈ 0 for statistically isotropic fields and departs from zero when the field has oriented structure.

normalized_coefficients(c) -> (; S0, s1, s2)

Amplitude-normalized coefficients: s1[j] = S1[j]/S0 and s2[j1,j2] = S2[j1,j2]/S1[j1] (zero where S1[j1] == 0).

log_coefficients(c; pad=eps) -> (; S0, logS1, logS2)

Log-coefficients log(S1 + pad), log(S2 + pad) (the small pad keeps structural zeros finite). Useful to gaussianize the heavy-tailed coefficients of intermittent fields.

scattering_batch(st::Scattering1D.ScatteringTransform1D, X) -> Matrix

Apply a 1D scattering transform to a batch of signals X of size (N, B) (signals as columns), returning a (flatten_length, B) matrix whose column b is flatten1d(st(X[:, b])). The plan and all workspace buffers are reused across the batch (only the small scalar-S0 wrapper is re-allocated per column).

scattering_batch(st::Scattering2D.ScatteringTransform2D, X) -> Matrix

Apply a 2D scattering transform to a batch of images X of size (Ny, Nx, B), returning a (flatten_length, B) matrix whose column b is flatten2d(st(X[:, :, b])). Plan and workspace are reused across the batch.

scattering_batch(st::Scattering3D.ScatteringTransform3D, X) -> Matrix

Apply a 3D scattering transform to a batch of volumes X of size (Nz, Ny, Nx, B), returning a (flatten_length, B) matrix. Plan and workspace are reused across the batch.

scattering_batch(st::SubsampledScattering.MultiResolutionScattering, X) -> Matrix
scattering_batch(st::ScatteredPlanar.ScatteredPlanarScattering, X) -> Matrix

Transform every slice of X against one transform, returning a (flatten_length, B) matrix. X is (dims…, B) for a multi-resolution transform and (M, B) for the scattered planar one, where M is the plan's point count.

scattering_batch(st::SphericalCore.SphericalScattering, X) -> Matrix
scattering_batch(st::SphericalCore.SphericalMonogenicScattering, X) -> Matrix

Transform a stack of B fields sampled by the same spherical plan — X of size (field_size…, B), so (M, B) for a scattered point set and (nθ, nφ, B) on a structured grid.

Rows follow the flat layout of Coefficients.flat_length and its row accessors. The spherical cascade pairs each scale with strictly coarser ones (j2 < j1), so a pair lands in the row that layout assigns to the unordered pair, flat_row_s2(j2, j1, J).

scattering_batch!(out, st, X; workspace = nothing) -> out

In-place counterpart of scattering_batch: write the flattened coefficients of each slice of X into the preallocated (flatten_length, B) matrix out. Backend-dispatched ! methods (scattering_batch!(out, backend, st, X)) are added by the corresponding extensions.

The default transforms one slice at a time against the transform's own single-slice plan. Passing a batch_workspace instead runs the whole stack through one batched plan (Batched.batch_cascade!).

Per-slice is the default because it is faster on a CPU at every batch size measured — 1.2–1.5× serial and 4.5× threaded. The cascade is memory-bandwidth bound and issues O(nw + paths) operations against the same data, so what governs is whether that data stays resident: a slice does, a B-slice stack does not. The batched plan's amortised per-call overhead does not recover the difference. It wins on a device, where the batch is what fills the machine — which is why the GPU extension builds one.

batch_coeffs(st, T) -> coefficient container

A coefficient container for st whose S0 is a 1-element array, so a loop over a batch writes every field in place instead of rebuilding the struct once per slice.

flat_rows(st) -> Int

Rows a flattened coefficient column of st occupies — the height of scattering_batch's output.

batch_workspace(st, B; spectral = ..., fft_nthreads = 1) -> Batched.BatchWorkspace

Build the batched workspace for transform st at batch size B: a spectral plan over the whole (spatial…, B) stack plus the scratch Batched.batch_cascade! runs in.

The plan is rebuilt because a batched transform is a different plan from a single-slice one — that is exactly what makes it one execution instead of B. Build it once and pass it back in through scattering_batch! when transforming many stacks of the same size.

fft_nthreads is the FFT library's own thread count, baked into the plan. It defaults to 1: the cascade is memory-bandwidth bound and the transform is only ~a quarter of a cascade step, so threading inside the FFT is capped near 1.3× by Amdahl no matter how many cores it gets, while threading over the batch or wavelet axis parallelises the whole step. Measured on 8 threads it returns 1.08–1.17× on large stacks and 0.33× on small ones, where the library's per-execution task spawn (which also allocates) exceeds the transform itself. Raise it only for a single stack large enough to be worth it, transformed with no parallelism above.

BatchWorkspace{P,RA,CA,WV,RV,G}

Preallocated state for repeated batched transforms at a fixed batch size B: the batched plan, the (spatial…, B) scratch, the filters in whatever storage the plan expects, and the precomputed order-2 groups. Once built, a batched transform does no data-proportional allocation.

red is the spatial-reduction target, shaped (1…, 1, B) so sum! contracts every axis but the batch.

batch_cascade!(out, ws, X) -> out

Write the flattened scattering coefficients of every slice of X into the columns of out.

out has Coefficients.flat_length(n) rows; row assignment goes through Coefficients.flat_row_*, the same layout flatten1d!/flatten2d! produce, so batched and per-slice results are interchangeable.

MultiResolutionScattering{T,D}

Scattering transform with a multi-resolution second order, on a D-dimensional grid. Build one with SubsampledScattering1D, SubsampledScattering2D or SubsampledScattering3D; they differ only in which filter bank they build, and share this type and one cascade.

SubsampledScattering1D([T=Float64,] N, J; Q=1, max_order=2, oversampling=1,
                       spectral=AutoSpectralBackend())

1D scattering with a multi-resolution second order: the first-order modulus envelope is decimated by 2^(scale - oversampling) before the second wavelet transform. Opt-in and approximate — a large oversampling reproduces the exact ScatteringTransform1D; aggressive values trade a little accuracy for speed.

SubsampledScattering2D([T=Float64,] (Ny, Nx), J; L=8, max_order=2, oversampling=1,
                       spectral=AutoSpectralBackend())

2D counterpart of SubsampledScattering1D. Decimating by ds shrinks the second order's grid by ds².

SubsampledScattering3D([T=Float64,] (Nz, Ny, Nx), J; n_orient=6, max_order=2, oversampling=1,
                       spectral=AutoSpectralBackend())

3D counterpart of SubsampledScattering1D. Decimating by ds shrinks the second order's grid by ds³, which is where this path pays off most.

subsampled_scattering!(coeffs, st, field) -> coeffs

In-place multi-resolution scattering into a pre-allocated coefficient container. Allocation-free.

Grouped by first-order wavelet, so each U₁ is decimated and transformed once and reused by every one of its children, rather than once per order-2 path.

FilterBank1D{T,V,W,MV}

Complete 1D filter bank for the scattering transform. Every container is a type parameter (no hardcoded Vector): V the per-filter array type (CPU/GPU/static/…), W the wavelet collection, MV the metadata collection.

Fields

• wavelets::W: wavelet filters in the Fourier domain (W<:AbstractVector{V})
• averaging::V: low-pass averaging (scaling) filter
• meta::MV: per-wavelet WaveletMeta
• J::Int: number of octaves (scales)
• Q::Int: wavelets per octave

FilterBank2D{T,M,W,MV}

Complete 2D filter bank with oriented wavelets. Containers are type parameters (no hardcoded Vector): M the per-filter matrix type, W the wavelet collection, MV the metadata collection.

Fields

• wavelets::W: oriented wavelet filters (W<:AbstractVector{M})
• averaging::M: low-pass averaging filter
• meta::MV: per-wavelet WaveletMeta
• J::Int: number of scales
• L::Int: number of orientations

FilterBank3D{T,A<:AbstractArray{Complex{T},3}}

Complete 3D oriented Morlet filter bank: J scales × n_orient sphere directions, plus a low-pass averaging filter.

WaveletMeta{T}

Concrete per-wavelet metadata: a struct rather than a NamedTuple, so the container stays concretely typed.

Fields

• scale::Int: octave index j
• q::Int: sub-octave index within the octave (1D, 0..Q-1); 0 for 2D
• orient::Int: orientation index l (2D, 0..L-1); 0 for 1D
• j_eff::T: effective log-scale used to order paths. j + q/Q in 1D, T(j) in 2D. The second-order admissibility constraint is j_eff(child) > j_eff(parent) (frequency strictly decreasing) — which for 2D means scale strictly increasing over all orientation pairs.
• center_freq::T: wavelet center frequency
• theta::T: orientation angle in radians (2D); 0 for 1D

build_filter_bank1d(N::Int, J::Int; Q::Int=1) -> FilterBank1D

Build a 1D Morlet filter bank with dyadic scales.

Arguments

• N::Int: Signal length (FFT size)
• J::Int: Maximum scale (number of octaves)
• Q::Int: Wavelets per octave (default 1 for dyadic, 8 for high Q)

Returns

• FilterBank1D: Complete filter bank with J scales

build_filter_bank2d(N::NTuple{2,Int}, J::Int; L::Int=8) -> FilterBank2D

Build a 2D oriented Morlet filter bank.

Arguments

• N::NTuple{2,Int}: Image dimensions (Ny, Nx)
• J::Int: Number of dyadic scales
• L::Int: Number of orientations (default 8, evenly spaced)

Returns

• FilterBank2D: Complete 2D filter bank

build_filter_bank3d(N::NTuple{3,Int}, J::Int; n_orient::Int=6, T=Float64) -> FilterBank3D

Build a 3D oriented Morlet filter bank with J dyadic scales and n_orient near-uniform orientations on the sphere (Fibonacci spiral).

Morlet1D{T<:Real}

1D Morlet wavelet in frequency domain.

The Morlet wavelet is a complex sinusoid modulated by a Gaussian: ψ(x) = (1/√|Σ|) exp(-x²/(2σ²)) (exp(i k₀ x) - β)

where β = exp(-σ²k₀²/2) ensures zero mean (admissibility condition).

In frequency domain: Ψ(ω) = exp(-(ω-k₀)²σ²/2) - β exp(-ω²σ²/2)

Type Parameters

• T: Element type (Float32, Float64, etc.)

Fields

• center_freq::T: Center frequency k₀
• bandwidth::T: Standard deviation σ of Gaussian envelope
• beta::T: Correction factor for zero mean
• N::Int: Filter length (FFT size)

Morlet2D{T<:Real}

2D oriented Morlet wavelet in frequency domain.

The 2D Morlet wavelet is created by taking a 1D Morlet and rotating it to angle θ, with elliptical Gaussian envelope controlled by elongation.

Type Parameters

• T: Element type (Float32, Float64, etc.)

Fields

• center_freq::T: Center wavenumber |k₀|
• bandwidth_x::T: Bandwidth along major axis
• bandwidth_y::T: Bandwidth along minor axis (controls elongation)
• theta::T: Orientation angle in radians
• beta::T: Correction factor
• N::NTuple{2,Int}: Filter dimensions (Ny, Nx)

Morlet3D{T<:Real}

3D oriented Morlet wavelet in the frequency domain, a bump centered at k₀ n̂ for a unit direction n̂ on the sphere, with an anisotropic Gaussian envelope (std σ∥ along n̂, σ⊥ = σ∥/elongation perpendicular) and analytic on the half-space k·n̂ ≥ 0.

Fields

• center_freq::T: |k₀|
• sigma_par::T, sigma_perp::T: real-space envelope widths along / perpendicular to n̂
• direction::NTuple{3,T}: unit orientation n̂
• beta::T: zero-mean correction
• N::NTuple{3,Int}: grid dimensions

frequency_response(m::Morlet2D{T}) -> Matrix{Complex{T}}

Compute the 2D frequency response Ψ(kx, ky) of an oriented Morlet wavelet. Element type matches the wavelet's precision.

frequency_response(m::Morlet3D{T}) -> Array{Complex{T},3}

3D frequency response Ψ(kx,ky,kz) of an oriented Morlet wavelet.

fibonacci_directions(n, ::Type{T}=Float64) -> Vector{NTuple{3,T}}

n near-uniform unit directions on the sphere (Fibonacci spiral), used as 3D wavelet orientations.

ScatteringTree{IV,RV}

Flat CSR description of the scattering tree. Pure integer topology (indices into the filter bank + integer order labels) — outside the differentiable data path, so the values are integers, but the containers stay parametric (default builders produce Vector{Int} / Vector{UnitRange{Int}}, yet MArray/CuVector/Int32 storage is permitted).

Fields

• path_data: concatenated wavelet-index lists for all paths
• path_ptr: CSR offsets (length npaths+1); path p is path_data[path_ptr[p]:path_ptr[p+1]-1]
• order: scattering order of each path (0, 1, 2, …)
• by_order: by_order[o+1] is the contiguous range of path ids of order o

build_tree(j_eff, max_order) -> ScatteringTree

Enumerate all admissible paths up to max_order from the per-wavelet effective log-scales j_eff (typically [m.j_eff for m in filter_bank.meta]). Admissibility: j_eff strictly increasing along the path. Paths are laid out grouped by order (order 0, then 1, then 2, …), so by_order ranges are contiguous.

order2_groups(tree, nw) -> Vector{Tuple{Int,Vector{Int},Vector{Int}}}

Every first-order wavelet j1 paired with its admissible order-2 children j2 and the tree path ids of those (j1, j2) pairs, ordered longest-first. The path ids are what the localized-field output is indexed by; the coefficient cascade only needs (j1, j2).

The cascade walks this instead of the flat order-2 path list so that a first-order modulus is computed and transformed once per j1 and reused across its children, rather than recomputed per path. Longest-first ordering is for the parallel backends: the child count varies by roughly 3× across scales, so equal-sized chunks of a natural-order list load-imbalance badly.

order_range(tree, o) -> UnitRange

Contiguous range of path ids with scattering order o.

path_indices(tree, p) -> view

Wavelet indices of path p (empty for the order-0 root), as a non-allocating view.

npaths(tree) -> Int

Total number of paths (including the order-0 root).

AbstractScatteringPlan

Supertype for spectral transform plans. A plan implements forward_transform! and inverse_transform!; the in-core default is DirectSumPlan, with fast paths (FFTW, AbstractFFTs/device, FINUFFT, NonuniformFFTs) provided by extensions.

DirectSumPlan{T,V,D,S}

Direct-summation DFT plan: evaluates X_k = Σ_n x_n e^{-2πi kn/N} (and its inverse) by direct summation, separably over the leading D dimensions, leaving a trailing batch axis of length nbatch untouched. Memory is O(N) — a per-axis table of roots of unity and its conjugate, never an N×N matrix. scratch is nothing for D == 1 and a full-size complex array otherwise. Containers stay parametric.

forward_transform!(out, plan, x) -> out

In-place forward (fft-convention) spectral transform. Methods provided by concrete plans.

inverse_transform!(out, plan, x) -> out

In-place inverse (ifft-convention, 1/N-scaled) spectral transform.

forward_transform(plan, x) -> X̂
inverse_transform(plan, x) -> x

Non-mutating, allocating, element-type-generic spectral transforms — the autodiff-friendly counterparts of the in-place forward_transform!/inverse_transform!. They never touch the plan's preallocated scratch and never mutate their inputs, so they accept ForwardDiff.Dual/Float32 inputs and are differentiable by reverse-mode backends. Used by the non-mutating scattering(st, x) path; the in-place ! versions remain the production hot path.

DirectSumPlan implements these as dense per-axis DFT matrix-multiplies (W*x) — differentiable by every AD backend with no special rules. The matrices are built once on first use and cached on the plan, so the O(N²) build is not repeated per call.

inverse_transform(plan, x) -> x

See forward_transform.

make_plan(spectral, T, dims; nbatch=1, kwargs...) -> AbstractScatteringPlan

Build the spectral plan selected by spectral for arrays whose leading dimensions are dims and whose element type is Complex{T}, with a trailing batch axis of length nbatch.

SpectralBackends.DirectSumSpectralBackend is the dependency-free in-core default; SpectralBackends.FFTSpectralBackend requires using FFTW; SpectralBackends.AutoSpectralBackend takes the FFTW fast path when its extension is loaded and otherwise the in-core direct sum.

make_scattered_plan(spectral, x, y, ms, T; period, solve, maxiter, rtol, eps) -> AbstractScatteringPlan

Build the scattered/nonuniform planar plan selected by spectral over points (x, y) and a uniform mode grid of size ms. SpectralBackends.DirectSumSpectralBackend is the dependency-free exact NUDFT; FINUFFTBackend and NonuniformFFTsBackend select a specific fast library; SpectralBackends.NUFFTSpectralBackend takes whichever fast library is loaded, and SpectralBackends.AutoSpectralBackend falls back to the exact direct sum when neither is.

spectral_backend(plan) -> SpectralBackends.AbstractSpectralBackend

The tag that would rebuild plan. Each plan type declares its own, so a transform can be reconstructed faithfully on a remote worker; plans that cannot be rebuilt from a tag alone (a device-resident FFT plan needs its device too) throw rather than report a host plan.

task_local_plan(plan) -> plan

A plan equivalent to plan that is safe to use concurrently with it. Stateless plans (FFTW, AbstractFFTs) return themselves; plans carrying mutable scratch return a copy that shares their read-only tables and owns fresh scratch. Called once per task by the parallel backends.

with_fft_nthreads(f, n) -> f()

Run f with the FFT library's global thread count set to n, restoring it afterwards.

FFTW's thread count is process-global and is raised as a side effect of loading unrelated packages, so a plan built without pinning it inherits whatever was last set — and a plan built for more threads than it needs spawns (and allocates) a task per thread on every execution. Plan builders wrap construction in this so a plan's threading is a property of the plan, not of load order.

The default is a no-op: only the FFTW extension has a global count to set. It takes args... so the extension's fixed-arity method is strictly more specific and adds a method rather than overwriting this one — overwriting is an error during precompilation.

wavelet_convolve(signal_fft, filter_fft, plan)

Perform wavelet convolution via frequency-domain multiplication then inverse transform. Allocates output. For zero-allocation hot paths, use wavelet_convolve!.

wavelet_convolve!(out, signal_fft, filter_fft, plan, buffer)

Truly zero-allocation wavelet convolution.

Multiplies signal_fft .* filter_fft into buffer in-place, then applies the inverse spectral transform via Plans.inverse_transform!(out, plan, buffer) — writing directly into out.

out and buffer must both be pre-allocated complex arrays of the same size.

apply_modulus(signal)

Apply complex modulus |·| to get envelope. Allocates output. For zero-allocation hot paths, use apply_modulus!.

apply_modulus!(out, signal)

In-place modulus. Stores |signal| in pre-allocated out. Zero allocation.

modulus_mean(signal) -> Real

⟨|signal|⟩ in a single reduction. A scattering coefficient is the mean of a modulus, so the modulus field itself is never needed unless a coarser scale consumes it — this is the leaf case, which writes nothing.

modulus_mean!(out, signal) -> Real

Write |signal| into out and return ⟨|signal|⟩. Used where the modulus field is consumed downstream; the generic method is two device-friendly passes, the CPU method fuses them into one.

spatial_average(signal::AbstractArray{T}) -> T

Compute spatial average (global mean) for translation invariance. Type-stable: returns element type T.

task_workspace(st) -> st′

A transform equivalent to st that shares its read-only parts — filter bank, path tree, work list — but owns fresh buffers and a task-local spectral plan, so the two can run concurrently.

This is what lets a parallel backend give each task private scratch without duplicating the filter bank, which dominates a transform's memory (for a 256×256 J=4 L=8 transform, 33 MiB of the 38 MiB). Methods are defined per transform type.

fftw_plan(T, dims; nbatch, planning, fft_nthreads) -> AbstractScatteringPlan

Build the FFTW-backed plan. The real method lives in the FFTW extension; the definition here is a throwing stub, so an explicit FFTSpectralBackend() request dispatches straight to the extension when it is loaded and to an actionable error when it is not — with no capability lookup on either path.

abstractffts_plan(dummy_device_array; region=1:ndims(dummy_device_array))

Vendor-neutral device FFT plan constructor. Declaration only — the sole method is provided by the KernelAbstractions extension, which builds forward/inverse plans via AbstractFFTs.plan_fft / plan_ifft on dummy_device_array. Because those dispatch on the array type, the same builder yields cuFFT (CuArray), rocFFT (ROCArray), or FFTW (plain Array) plans — so the GPU scattering path is device-agnostic. region selects the transformed dimensions (e.g. (1, 2) for a batched (Ny, Nx, B) stack, leaving the batch axis untouched).

finufft_scattered_plan(x, y, ms, T; period, solve, maxiter, rtol, eps)
nonuniformffts_scattered_plan(x, y, ms, T; period, solve, maxiter, rtol, eps)

Fast-path scattered-planar plan constructors. The real methods live in the FINUFFT and NonuniformFFTs extensions; the definitions here are throwing stubs, so naming one of those backends explicitly costs a plain dispatch rather than a capability lookup.

See finufft_scattered_plan.

FINUFFT fast path for scattered/nonuniform planar points; requires using FINUFFT.

NonuniformFFTs.jl fast path for scattered/nonuniform planar points; requires using NonuniformFFTs.

resolve_backend(backend) -> AbstractLocalBackend

Concrete backends pass through unchanged — a request is honoured exactly or refused, never silently downgraded. AutoBackend resolves on real capability: ThreadedBackend when there is more than one thread and the OhMyThreads extension is loaded, otherwise SerialBackend.

check_available(backend) -> backend

Throw unless backend can actually execute, naming the package that would enable it. Called by the batch entry points so an unloadable request fails immediately instead of at a MethodError.

true when the OhMyThreads extension is loaded, i.e. ThreadedBackend can actually run.

true when the KernelAbstractions extension is loaded, i.e. GPUBackend can actually run.

true when the Distributed extension is loaded.

true when the MPI extension is loaded.

cascade!(S1, S2, st, signal_fft) -> (S1, S2)

Both scattering orders in one pass over the tree, grouped by first-order wavelet:

for (j1, children):  U₁ = |x ⋆ ψ_j1| ;  S1[j1] = ⟨U₁⟩
                     Û₁ = fft(U₁)     ;  S2[j1,j2] = ⟨|U₁ ⋆ ψ_j2|⟩  for each child

The first-order convolution is therefore evaluated once, not once for S1 and again for S2, and only one U₁/Û₁ pair is live at a time rather than one per wavelet. Wavelets with no admissible child skip the modulus buffer entirely, reducing to a single fused ⟨|·|⟩.

signal_fft is read only, so the caller's preserved signal spectrum survives the call.

cascade!(S1, S2, st, image_fft) -> (S1, S2)

Both scattering orders in one grouped pass — see the 1D cascade! for the scheme. Admissibility is j_eff strictly increasing, i.e. scale strictly increasing across all orientation pairs, which is what the tree encodes; same-scale different-orientation pairs are not order-2 paths.

cascade!(S1, S2, st, vol_fft) -> (S1, S2)

Both scattering orders in one grouped pass — see the 1D cascade! for the scheme.

(st::ScatteredPlanarScattering)(x) -> ScatteringCoefficients2D

Apply the scattered planar scattering transform to a length-M vector of samples at the plan's points. Allocates a coefficient container per call; use scattered_planar_scattering! to reuse one.

scattered_planar_scattering!(coeffs, st, x) -> coeffs

In-place counterpart of the callable: write into a preallocated ScatteringCoefficients2D. Allocation-free when coeffs carries a 1-element S0 (so update_S0 writes rather than rewraps).

Level{WV,P,CA,RA}

One resolution of the multi-resolution cascade: the wavelet bank and spectral plan on the decimated grid, plus the buffers the second order runs in there. Buffers live on the level, so a path costs no allocation.

AbstractSphericalPlan

Supertype for spherical spectral plans (scattered NUFSHT, structured SHT, …). A concrete plan must implement sphere_coeffs, sphere_apply!, and sphere_mean.

sphere_coeffs(plan, field) -> C

Spherical-harmonic analysis: the coefficients C of field (up to the plan's band limit). For a structured grid this is the fast forward SHT; for scattered points it is the exact (least-squares / CG) inversion — not the adjoint, which mis-scales the coefficients. Returned opaquely and consumed only by sphere_apply!/sphere_mean on the same plan.

sphere_coeffs!(C, plan, field) -> C

In-place sphere_coeffs: analyse field into the pre-allocated coefficient container C, which must have come from sphere_coeffs_buffer(plan). This is what lets the cascade re-analyse a first-order field without allocating a coefficient vector per scale.

sphere_coeffs_buffer(plan) -> C

A coefficient container of the right type and size for plan, suitable for sphere_coeffs!. Backends whose coefficients are not a plain vector (the FastSphericalHarmonics triangular layout, the NUFSHT dense spin layout) return their own shape.

sphere_apply!(out, plan, C, h) -> out

Apply the per-degree multiplier h(ℓ) to (a copy of) the coefficients C and synthesise the result into out — i.e. out = Σ_{ℓm} h(ℓ) · C_{ℓm} · Y_{ℓm}. Does not mutate C.

sphere_mean(plan, field) -> scalar

Spherical average of field under the plan's sampling: an unweighted sample mean for (quasi-uniform) scattered points, the exact quadrature integral for a structured grid. Provided by each backend.

SphericalScattering{T,P,V}

Spherical scattering transform over a backend plan::P. sigma2[k+1] is the Gaussian-transfer variance for the dyadic low-pass cutoff ℓ_k (k=0..J); band-pass wavelet j is lowpass(ℓ_j) − lowpass(ℓ_{j-1}).

SphericalMonogenicScattering{T,P}

Spherical monogenic scattering: shares the dyadic difference-of-Gaussians bands of SphericalScattering but replaces the analytic modulus with the spherical monogenic amplitude A_j = √(U⁰_j² + |∇_S g_j|²) (spin-0 Bochner identity — see monogenic_amplitude!).

SphericalWorkspace{A,C}

Scratch for one spherical cascade: the band-pass output band, the current first-order field u1, and two coefficient containers — C for the input field, C1 for the first-order field being re-analysed. Two is all the cascade ever needs, because it finishes every child of a scale before starting the next.

Build one with SphericalWorkspace(st, field) and reuse it across calls; a task that transforms concurrently needs its own (and its own Plans.task_local_plan of the spherical plan).

SphericalMonogenicWorkspace{A,C,S}

Scratch for one spherical monogenic cascade: the current amplitude field u1, a second amp for the order-2 amplitudes, two coefficient containers, and scratch — the (g, lapg, g2, lapg2) fields the Bochner identity needs (see monogenic_amplitude!).

spherical_scattering!(S1, S2, st, ws, field) -> (; S0, S1, S2)

In-place spherical scattering into pre-allocated S1/S2 using the workspace ws (see SphericalWorkspace). Allocation-free once ws exists.

Grouped by first-order scale so only one first-order field and one coefficient vector are live at a time: scale j1 is band-passed, averaged into S1[j1], analysed once, and then consumed by every coarser j2 < j1 before the next j1 begins.

spherical_monogenic_scattering!(S1, S2, st, ws, field) -> (; S0, S1, S2)

In-place spherical monogenic scattering — the monogenic counterpart of spherical_scattering!, grouped by first-order scale so one amplitude field is live at a time rather than all J.

monogenic_amplitude!(amp, st, C, j, w) -> amp

Spherical monogenic amplitude at scale j of the field whose (already-computed) SH coefficients are C, written into amp. w is a NamedTuple of scratch fields (g, lapg, g2, lapg2).

Uses only spin-0 transforms via the Bochner/product identity: with g_j = (−Δ_S)^{-1/2} U⁰_j,

|U^R_j|² = |∇_S g_j|² = ½ Δ_S(g_j²) − g_j · Δ_S g_j,

so A_j = √(U⁰_j² + |∇_S g_j|²). The Riesz vector itself (needed for orientation/phase) requires spin-1 synthesis and is handled separately by spherical_monogenic_components.

task_local(st) -> st

A copy of the spherical transform safe to run concurrently with the original: the design matrix, Gram factor and point set are read-only and shared, while the plan's analysis/synthesis scratch is duplicated by Plans.task_local_plan. Analysis writes through that scratch, so tasks sharing one plan would overwrite each other's coefficients.

Per-degree multiplier for the difference-of-Gaussians band-pass wavelet j.

dog_sigma2(lmax, J, T) -> Vector{T}

Gaussian-transfer variances for the J+1 dyadic low-pass cutoffs ℓ_k = lmax / 2^(J-k), k=0..J (ℓ_J = lmax). Band-pass wavelet j is lowpass(ℓ_j) − lowpass(ℓ_{j-1}).

structured_grid(lmax, T) -> (Θ, Φ)

Colatitudes Θ (length lmax+1) and longitudes Φ (length 2lmax+1) of the equiangular structured grid. Matches FastSphericalHarmonics.sph_points(lmax+1).

make_spherical_plan(spectral, θ, φ, lmax, T; rtol, maxiter) -> AbstractSphericalPlan

Build the scattered-sphere plan selected by spectral over points (θ, φ) and band limit lmax. SpectralBackends.DirectSumSpectralBackend is the dependency-free default; SpectralBackends.NUFSHTSpectralBackend (or SpectralBackends.AutoSpectralBackend once the NUFSHT extension is loaded) uses the NUFSHT fast path.

make_structured_plan(spectral, lmax, T; rtol, maxiter) -> AbstractSphericalPlan

Build the structured-sphere plan selected by spectral. SpectralBackends.DirectSumSpectralBackend is the dependency-free default (direct SHT on the grid); SpectralBackends.FSHTSpectralBackend (or SpectralBackends.AutoSpectralBackend once the FastSphericalHarmonics extension is loaded) uses the fast exact SHT.

nusht_spherical_plan(θ, φ, lmax, T; rtol, maxiter)

Fast-path scattered-sphere plan constructor. The real method lives in the NUFSHT extension; this is its throwing stub.

fsh_structured_plan(lmax, T)

Fast-path structured-sphere plan constructor. The real method lives in the FastSphericalHarmonics extension; this is its throwing stub.

plot_coefficients(c; …) — plot scattering coefficients. Requires `using CairoMakie`.

plot_filter_bank(fb) — plot a filter bank. Requires `using CairoMakie`.