Multiscale#

Overview#

The multiscale module in SplineOps models signals and images as a hierarchy of spline approximations at progressively coarser resolutions. On top of this model, it provides reduction (spline filtering + dyadic decimation) and expansion (upsampling + spline interpolation), which are the building blocks of pyramid and wavelet transforms [1], [2], [3], [4].

Spline Representation#

A 1D discrete signal \(\{f[k]\}\) can be modeled as the continuous function

\[f(x) \;=\; \sum_{k} c[k]\, \phi\bigl(x - k\bigr),\]

where \(\phi(x)\) is a polynomial spline basis function (e.g., a B-spline of degree 3), and \(c[k]\) are the spline coefficients determined from the samples \(f[k]\). This representation allows the application of downsampling and upsampling filters directly to the spline model.

Pyramid Decomposition#

Two key operators are proposed:

  • Reduce: it filters the signal (or image) and downsamples by dyadic factors, thus producing a coarse approximation.

  • Expand: it upsamples and interpolates the coarse approximation back to the original resolution.

When applied iteratively, these operations create a pyramid structure (approximation at multiple scales). In 2D, the same concept applies along rows and columns.

Subpixel Registration#

Method [5] does subpixel registration as least-squares matching of image intensities under a global transform (affine, optionally restricted to rigid/similarity) with an optional contrast change, and solves it with a modified Levenberg-Marquardt optimizer.

Its key ingredient is a spline pyramid: starting from a dyadic, least-squares fine-to-coarse decomposition built with cubic-spline filtering, the optimizer estimates the transform at the coarsest level and propagates the parameters down the pyramid so that finer levels only apply small corrections. The spline model (cubic) is used consistently for resampling and for computing exact spatial derivatives, which further stabilizes the coarse-to-fine.

A compact Python implementation following this approach is available in the GitHub repository glichtner/pystackreg.

Wavelet Decomposition#

We construct a spline-based multiscale basis (wavelets) by capturing the detail lost at each reduction step. For the explicitly supported shape and scale combinations, the wavelet (detail) coefficients together with the final coarse approximation allow perfect reconstruction (synthesis).

The next figure, from Wavelet Decomposition, shows a three-level 2D Haar decomposition: the coarse approximation in the top-left corner and the horizontal, vertical, and diagonal detail sub-bands at each scale.

../_images/sphx_glr_02_wavelet_decomposition_001.png

At each scale (analysis):

  • an approximation is obtained (the reduced signal or image);

  • a corresponding detail or wavelet sub-band is formed (the difference or “error” relative to the expanded approximation).

The application of this decomposition over multiple scales yields a so-called wavelet representation, where the stored approximation plus the detail coefficients can be used to perfectly reconstruct the original data (synthesis).

Implementation Details#

  • Reduce and expand perform the core downsampling and upsampling with spline filters. Whole-array axis operations replace row-by-row and column-by-column Python dispatch.

  • Haar split/merge operations write into preallocated scale regions. Small independent planes are transformed in vectorized cache-sized groups; large planes are dispatched directly so a batch does not pay for full-array transpose copies. Spline wavelets retain short, explicit loops over filter taps.

  • Various spline degrees (e.g., degree 3) control how data are dispatched between the approximation channel and sub-bands.

Supported shapes and boundaries#

The pyramid functions accept finite, non-empty real 1-D signals and 2-D arrays and use their documented mirror mappings. Floating inputs preserve their precision; integer inputs promote to float64, and booleans are rejected. Reducing an odd length returns floor(n / 2) samples; expanding that result therefore does not recover the dropped extent. A singleton is preserved exactly.

The multiscale Haar and spline-wavelet classes currently support non-empty 2D arrays whose two dimensions are divisible by 2**scales. Unsupported odd or too-small scale regions are rejected instead of silently losing samples. Perfect reconstruction is tested for even square and rectangular arrays for Haar and the cubic spline-wavelet implementation. The order-1 spline filter reconstructs the audited float64 cases within 2e-7. The inherited order-5 filter coefficients contain only roughly five to six significant digits and produce errors up to about 2e-3 in the current randomized rectangular audit. Order 5 is therefore an approximate research implementation, not a perfect-reconstruction transform. The test suite records that limitation so it cannot silently become a stronger claim. The upstream DeconvolutionLab2 source contains the same table; because no authoritative higher-precision taps were found, SplineOps retains those values and exposes exact_reconstruction=False with a "bounded-approximation" contract.

These APIs never infer batch or channel dimensions. reduce_2d, expand_2d, and multi-scale wavelet analysis/synthesis accept explicit spatial_axes and transform every remaining slice independently. Execution adapts to the selected plane working set, while the divisibility and reconstruction contracts apply only to selected axes. See Performance evidence for the reproducible row/column-oracle comparison and Reusable workflow recipes for a batched example.

Multiscale Examples#

References#