synthesis — Complete parameter reference#

Module: foscat.scat_cov, foscat.scat_cov2D, foscat.scat_cov1D

synthesis is the high-level entry point for field synthesis, denoising, and texture generation. It wraps the full L-BFGS-B optimisation loop and handles multi-resolution scheduling automatically.

result = scat_op.synthesis(
    image_target,
    reference     = None,
    nstep         = 4,
    seed          = 1234,
    Jmax          = None,
    edge          = False,
    to_gaussian   = False,
    use_variance  = True,
    synthesised_N = 1,
    input_image   = None,
    grd_mask      = None,
    in_mask       = None,
    iso_ang       = False,
    fft_ang       = False,
    fft_nharm     = 1,
    fft_imaginary = True,
    EVAL_FREQUENCY= 100,
    NUM_EPOCHS    = 300,
    scat_cov_method = 'eval',
    n_up          = 0,
)

How synthesis works#

The optimiser searches for a field \(u\) whose scattering-covariance statistics \(\Phi(u)\) match those of image_target:

\[u^* = \arg\min_u \; \mathcal{L}(u), \qquad \mathcal{L}(u) = \sum_k \frac{(\Phi(u)_k - \Phi(d)_k)^2}{\sigma_k^2}\]

where \(\sigma_k\) is the per-coefficient variance of image_target (when use_variance=True). Gradients are computed via PyTorch autograd and the minimiser is scipy.optimize.fmin_l_bfgs_b.

The multi-resolution schedule (nstep) runs this minimisation from coarse to fine resolution, using each step’s result as the warm start for the next.

Parameters#

image_target (required)#

The reference field whose statistics are to be reproduced.

  • HEALPix: ndarray shape (npix,) or (B, npix) with npix = 12·nside²

  • 2D: ndarray shape (H, W) or (B, H, W)

  • 1D: ndarray shape (N,) or (B, N)

The input should be normalised (zero mean, unit variance is standard):

xnorm = (image - np.mean(image)) / np.std(image)
result = scat_op.synthesis(xnorm, ...)

nstep (default: 4)#

Number of resolution levels in the multi-resolution optimisation.

The algorithm builds a resolution pyramid by repeated 2× downsampling of image_target. It then synthesises from coarsest to finest, using each result as the initial condition for the next level.

nstep

Resolutions visited (2D example, 256×256 target)

1

256×256 only (direct synthesis, slow convergence)

2

128×128 → 256×256

3

64×64 → 128×128 → 256×256

4

32×32 → 64×64 → 128×128 → 256×256

Capped automatically if nstep > jmax - 1 (the map is too small for that many downsampling steps). For a 256×256 2D map, the maximum useful nstep is about 6.

# Multiscale synthesis — recommended for maps larger than ~64×64
result = scat_op.synthesis(xnorm, nstep=3, NUM_EPOCHS=300)

# Single-scale (faster, but worse convergence at large resolutions)
result = scat_op.synthesis(xnorm, nstep=1, NUM_EPOCHS=1000)

n_up (default: 0)#

Number of extra upsampling steps beyond the target size, keeping the same Jmax (same wavelet scales).

With n_up=1, after completing synthesis at N×N, the algorithm continues at 2N×2N. The scattering statistics used are still those of the N×N target — only the domain is larger. The result is a field that locally matches the statistics of the original target, embedded in a larger canvas.

n_up

Output size (2D, N×N target)

Comment

0

N×N

Standard synthesis

1

2N×2N

Same statistics, larger domain

2

4N×4N

Two extra levels

The Jmax used during n_up steps is pinned to the value effective for the N×N target (computed from its size and KERNELSZ), so wavelet filters and the norm cache (P1_dic) remain consistent.

# Synthesise a 512×512 map whose statistics match a 256×256 target
result = scat_op.synthesis(xnorm_256, nstep=3, n_up=1, NUM_EPOCHS=300)
# result.shape == (512, 512)

NUM_EPOCHS (default: 300)#

Maximum number of L-BFGS-B iterations per resolution level.

The optimiser may stop earlier if the convergence criterion (factr) is reached. A good rule of thumb: convergence is typically achieved in 100–500 iterations per level for moderate-size maps.

The total wall-clock time scales as nstep × NUM_EPOCHS.

EVAL_FREQUENCY (default: 100)#

Print the current loss value every N L-BFGS-B iterations.

seed (default: 1234)#

Random seed for the Gaussian white noise used as the initial condition at the coarsest resolution level (k=0).

Change seed to generate statistically independent realisations of the same target:

results = [scat_op.synthesis(xnorm, seed=s, nstep=3) for s in range(10)]

Has no effect when input_image is provided (the initial condition is then taken from input_image).

synthesised_N (default: 1)#

Number of independent synthetic maps to produce in a single run (batch synthesis).

All synthesised_N maps are optimised simultaneously with the same loss, which amortises the per-iteration cost. The output has an extra leading dimension:

result = scat_op.synthesis(xnorm, synthesised_N=4, nstep=3)
# result.shape == (4, H, W)  for 2D

When synthesised_N > 1 the reference statistics are computed as the batch-mean of image_target (if batched) or from the single target map replicated for all synthesised_N maps.

use_variance (default: True)#

Controls the loss weighting.

  • True (recommended): each coefficient is divided by its standard deviation estimated from image_target. This makes the loss scale-invariant across spatial scales — large-scale and small-scale coefficients contribute equally regardless of their dynamic range.

    \[\mathcal{L}(u) = \sum_k \frac{(\Phi(u)_k - \Phi(d)_k)^2}{\sigma_k^2}\]

  • False: all coefficients are weighted equally. Dominated by the largest-amplitude statistics (usually the coarsest scales). Useful if you want to enforce exact coefficient values rather than normalised deviations.

Jmax (default: None)#

Maximum wavelet scale index included in the loss. When None, all scales available for the map size are used.

Reduce Jmax to match only small-scale statistics (fast textures) while leaving large-scale structure free:

# Use only scales j < 4 — ignore large-scale correlations
result = scat_op.synthesis(xnorm, Jmax=4, nstep=3)

During the multi-resolution schedule, Jmax at level k is automatically decremented by 1 for each coarser level: Jmax_k = Jmax - (nstep - 1 - k).

edge (default: False)#

Whether the map has non-periodic boundaries (limited spatial domain).

  • False: assumes periodic / full-sphere boundary conditions. The convolutions wrap around the edges.

  • True: activates an internal edge mask that down-weights pixels near the boundary at each resolution level. Use for rectangular images or partial-sky HEALPix patches where the field does not wrap.

Automatically set to True when in_mask is provided.

in_mask (default: None)#

A binary (or soft) mask that marks invalid pixels in the input data — pixels that should not contribute to the reference statistics.

Shape must match image_target. Pixels with in_mask = 0 are excluded from the scattering-covariance computation at every scale level (the mask is downsampled at each coarser level). Setting in_mask also enables edge=True internally.

Use case: partial-sky CMB analysis, images with missing regions, survey footprints.

mask = np.ones_like(image_target)
mask[bad_pixels] = 0.0
result = scat_op.synthesis(xnorm, in_mask=mask, NUM_EPOCHS=300)

grd_mask (default: None)#

A binary mask controlling which pixels of the synthesised map are free to move.

Pixels where grd_mask = 0 are frozen to their values in image_target; only pixels where grd_mask = 1 are updated by the optimiser.

Use case: inpainting (reconstruct missing regions while preserving observed pixels), or constrained synthesis where part of the field is fixed.

grd_mask = np.zeros_like(image_target)
grd_mask[missing_pixels] = 1.0   # only synthesise these pixels
result = scat_op.synthesis(xnorm, grd_mask=grd_mask, NUM_EPOCHS=500)

At each multi-resolution level the gradient mask is downsampled to match the current resolution.

reference (default: None)#

A second reference field for cross-scattering-covariance synthesis.

When provided, the loss matches the cross-statistics between the synthesised map and a fixed second field (rather than the auto-statistics of image_target):

\[\mathcal{L}(u) = \sum_k (\Phi_\text{cross}(u,\, d_2)_k - \Phi_\text{cross}(d_1,\, d_2)_k)^2\]

where \(d_1\) = image_target, \(d_2\) = reference.

Use case: component separation, where the synthesised field must match a specific covariance structure with a known companion field.

result = scat_op.synthesis(
    cmb_estimate,
    reference=dust_template,
    NUM_EPOCHS=300,
)

iso_ang (default: False)#

Whether to use isotropically averaged statistics in the loss.

  • False: the loss is computed on the full oriented statistics (S1, S2, S3, S3P, S4 with all orientation indices). The synthesised map can reproduce anisotropic structures (filaments, oriented textures).

  • True: the statistics are collapsed to their rotationally invariant content via iso_mean() before computing the loss:

    Stat

    Shape before iso

    Shape after iso

    Reduction

    S1, S2

    (..., L)

    (...)

    simple mean over L

    S3, S3P

    (..., L, L)

    (..., L)

    mean over l1 at fixed l2-l1

    S4

    (..., L, L, L)

    (..., L, L)

    mean over l1 at fixed (l2-l1, l3-l1)

    Use for fields that are statistically isotropic (no preferred direction). Reduces the number of constraints and accelerates convergence.

# Isotropic synthesis — suitable for CMB-like fields
result = scat_op.synthesis(xnorm, iso_ang=True, nstep=3)

Note

iso_ang and fft_ang should not be used together. iso_ang is the harder reduction (mean only); fft_ang keeps angular variation and is the recommended soft alternative.

fft_ang (default: False) / fft_nharm (default: 1) / fft_imaginary (default: True)#

A softer alternative to iso_ang that compresses the orientation axes to their first Fourier harmonics instead of collapsing them to a single mean. The three parameters work together:

Parameter

Type

Default

Role

fft_ang

bool

False

Enable Fourier angular compression

fft_nharm

int

1

Number of harmonics to keep beyond DC

fft_imaginary

bool

True

Keep both cos and sin components

What fft_ang=True keeps (with fft_nharm=1, fft_imaginary=True)#

Three coefficients per projected axis (basis functions \(\phi_k\)):

Index k

\(\phi_k(l)\)

Physical meaning

0

\(1/L\)

DC — mean over the projected axis (= iso_mean value)

1

\(\cos(2\pi l/L)\)

In-phase first harmonic

2

\(\sin(2\pi l/L)\)

Quadrature first harmonic

The amplitude \(A_1 = \sqrt{c_1^2 + s_1^2}\) is rotation-invariant regardless of how the image is oriented. This is why fft_imaginary=True is strongly recommended: with cosine only (fft_imaginary=False) a field oriented at 90° gives \(c_1 \approx 0\) even when strongly anisotropic.

Shapes and projection rules after fft_ang(nharm=1, imaginary=True) → nout = 3#

Statistic

Before

After

Projection axis

S1, S2

(…, L)

(…, 3)

the single orientation axis

S3, S3P

(…, L, L)

(…, L, 3)

l1 at fixed Δl = l2−l1

S4

(…, L, L, L)

(…, L, L, 3)

l1 at fixed (Δl12, Δl13)

For S3/S4 the projection is not a tensor product. The relative-orientation axes (Δl) are preserved exactly as in iso_mean, and only the absolute orientation axis l1 is projected:

\[\text{S3\_out}[\Delta l, k] = \sum_{l_1} \phi_k(l_1) \cdot S3[l_1,\,(l_1+\Delta l)\bmod L]\]
\[\text{S4\_out}[\Delta l_{12}, \Delta l_{13}, k] = \sum_{l_1} \phi_k(l_1) \cdot S4[l_1,\,(l_1+\Delta l_{12})\bmod L,\,(l_1+\Delta l_{13})\bmod L]\]

The \(k=0\) coefficient equals iso_mean exactly (\(\phi_0 = 1/L\), so \(\sum_{l_1}(1/L)\cdot S3[\ldots] = \text{mean}\)). Components \(k=1,2\) capture the angular modulation around that mean.

fft_ang vs iso_ang — choosing the right reduction#

iso_ang=True

fft_ang=True

S1/S2 output shape

(…) scalar

(…, 3)

S3/S3P output shape

(…, L)

(…, L, 3)

S4 output shape

(…, L, L)

(…, L, L, 3)

Angular information kept

mean only

mean + angular modulation amplitude

Number of constraints

minimum

moderate (+factor 3 on last axis)

Suited for

strictly isotropic fields

fields with preferred but variable orientation

Cost per iteration

lowest

slightly higher

 
# fft_ang synthesis — preserves angular variation amplitude
result = scat_op.synthesis(xnorm, fft_ang=True, nstep=3, NUM_EPOCHS=300)

# Keep two harmonics (DC + first two angular modes):
result = scat_op.synthesis(xnorm, fft_ang=True, fft_nharm=2, nstep=3)

# With fft_imaginary=False only if the image orientation is known and fixed:
result = scat_op.synthesis(xnorm, fft_ang=True, fft_imaginary=False, nstep=3)

to_gaussian (default: False)#

Apply a Gaussianisation transform to image_target before synthesis, then invert it at the end.

  • True: maps the histogram of image_target to a standard Gaussian before computing target statistics. The final output is mapped back to the original histogram. Useful for highly non-Gaussian fields (log-normal dust emission, etc.) where the scattering statistics alone do not fully constrain the one-point distribution.

  • False: the histogram of the synthesised map is controlled only by the scattering statistics (which capture some non-Gaussian content but are not guaranteed to reproduce the full marginal distribution).

input_image (default: None)#

A warm-start initial condition for the coarsest resolution level.

When None, the optimiser starts from Gaussian white noise (seeded by seed). When provided, the field input_image is downsampled to each resolution level and used as the initial guess at k=0.

Use case: iterative refinement, re-synthesis starting from a previous result, or injecting a prior.

# Refine a previous result
result_v1 = scat_op.synthesis(xnorm, nstep=3, NUM_EPOCHS=100)
result_v2 = scat_op.synthesis(xnorm, nstep=3, NUM_EPOCHS=500,
                               input_image=result_v1)

scat_cov_method (default: 'eval')#

Selects the internal method used to compute scattering covariances at each optimisation step.

  • 'eval' (default): uses funct.eval() with norm='auto'. The normalisation factors are computed from the target image on the first call and cached; subsequent iterations reuse the same normalisation. Works for all geometry types (2D, HEALPix, 1D). Supports all orientation reductions: iso_ang, fft_ang.

  • any other value (e.g. 'scattering_cov'): uses the direct scattering_cov() path. 2D only. Statistics are computed without the norm='auto' caching layer; normalisation is applied via the ref_sigma mechanism inside scattering_cov itself. Supports iso_ang and fft_ang (both with and without use_variance).

    Use this path when you want to bypass the eval normalisation cache — for example, in cases where the cached norm from the target image is not appropriate for the synthesised image (non-stationary textures, masked domains with very different statistics).

'eval'

'scattering_cov'

Geometry

2D, HEALPix, 1D

2D only

iso_ang

✓

✓ (native)

fft_ang

✓

✗ — raises ValueError

n_up

✓

✓ (ref recomputed at each resolution)

Normalisation

cached from target (norm='auto')

via ref_sigma at each step

use_variance

✓ (calc_var=True)

✓ (get_variance=True)

Note

fft_ang=True is incompatible with scat_cov_method != 'eval'. The scattering_cov path returns a flat concatenated tensor (all statistics serialised into a 1-D vector) rather than a structured scat_cov object, so the orientation-Fourier post-processing that fft_ang requires cannot be applied. A ValueError is raised immediately if this combination is attempted.

Return value#

  • If synthesised_N == 1: ndarray with the same shape as image_target (the leading batch dimension is squeezed).

  • If synthesised_N > 1: ndarray with a leading dimension of size synthesised_N.

  • If n_up > 0 (2D only): the output has shape (2^n_up · H, 2^n_up · W).

Complete examples#

Minimal 2D synthesis#

import foscat.scat_cov2D as sc
import numpy as np

scat_op = sc.funct(NORIENT=4)
xnorm = (image - np.mean(image)) / np.std(image)

result = scat_op.synthesis(xnorm, seed=10, nstep=3, NUM_EPOCHS=300)

Batch synthesis with masking#

mask = np.ones_like(xnorm)
mask[invalid] = 0.0

results = scat_op.synthesis(
    xnorm,
    in_mask       = mask,
    synthesised_N = 4,
    nstep         = 3,
    iso_ang       = True,
    NUM_EPOCHS    = 500,
)
# results.shape == (4, H, W)

Upsampled synthesis (n_up)#

# Synthesise a 512×512 map matching the statistics of a 256×256 target
result = scat_op.synthesis(
    xnorm_256,
    nstep      = 3,
    n_up       = 1,      # one extra 2× upsampling after full resolution
    NUM_EPOCHS = 300,
)
# result.shape == (512, 512)

Inpainting (frozen observed pixels)#

grd_mask = np.zeros_like(xnorm)
grd_mask[hole_pixels] = 1.0   # only these pixels move

result = scat_op.synthesis(
    xnorm,
    grd_mask   = grd_mask,
    nstep      = 3,
    edge       = True,
    NUM_EPOCHS = 500,
)

Cross-covariance synthesis (component separation)#

# Synthesise CMB-like field matching cross-statistics with a dust template
result = scat_op.synthesis(
    cmb_estimate,
    reference  = dust_template,
    nstep      = 3,
    iso_ang    = True,
    NUM_EPOCHS = 300,
)

Isotropic Gaussianised synthesis#

result = scat_op.synthesis(
    image_target,
    to_gaussian = True,
    iso_ang     = True,
    nstep       = 4,
    NUM_EPOCHS  = 500,
)