foscat.scat_cov#

This is a softer alternative to iso_mean(). Instead of collapsing each orientation axis to a single mean value, it keeps the first nharm harmonics of the discrete Fourier transform along each orientation axis L. This preserves the amplitude of the angular variation while reducing the number of descriptors from L to nout = 1 + 2*nharm (with imaginary=True) per orientation axis.

Why imaginary=True is the recommended mode

With imaginary=False the projection basis is {1, cos(2π·l/L), cos(4π·l/L), …}. A field whose dominant orientation sits at exactly 90° (i.e. at the zero-crossing of the cosine) would give a near-zero first-harmonic coefficient even though the field is strongly anisotropic.

With imaginary=True the basis is {1, cos(2π·l/L), sin(2π·l/L), cos(4π·l/L), sin(4π·l/L), …}. The amplitude of the first harmonic,

\[A_1 = \sqrt{c_1^2 + s_1^2}\]

is rotation-invariant: it is the same regardless of the absolute orientation of the field. Use imaginary=True whenever you want a description that does not depend on the image orientation.

Reduction per statistic

Statistic	Input shape	Output shape (nharm=1)
S1, S2	`[…, L]`	`[…, nout]`
S3, S3P	`[…, L, L]`	`[…, L, nout]`
S4	`[…, L, L, L]`	`[…, L, L, nout]`

For S1/S2 the projection applies the Fourier basis directly on the single orientation axis.

For S3/S3P and S4 the projection is not a tensor product. Instead, the statistics are first reindexed by the relative-orientation differences (exactly as in iso_mean()), then the absolute orientation axis \(l_1\) is projected onto the Fourier basis:

\[\text{S3\_out}[\Delta l,\, k] = \sum_{l_1} \phi_k(l_1)\, 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)\, S4[l_1,\,(l_1+\Delta l_{12})\bmod L,\,(l_1+\Delta l_{13})\bmod L]\]

where \(\phi_0(l)=1\), \(\phi_1(l)=\cos(2\pi l/L)\), \(\phi_2(l)=\sin(2\pi l/L)\).

This preserves the relative-orientation axes (same as iso_mean) while also capturing how strongly the statistics vary as the whole frame rotates. The \(k=0\) component is identical to the iso_mean() result.

Parameters:

nharm (int, optional) – Number of harmonics to keep beyond the DC term. nharm=1 (default) keeps the mean and the first angular harmonic.
imaginary (bool, optional) – If False (default), keep only the cosine components {1, cos, cos(2·), …} — nout = 1 + nharm. If True, keep both cosine and sine components {1, cos, sin, cos(2·), sin(2·), …} — nout = 1 + 2·nharm. Recommended: ``True`` to obtain rotation-invariant amplitudes.

Returns:

scat_cov – A new statistics object with orientation axes compressed from L to nout.

Examples

>>> stat = scat_op.eval(image)
>>> stat_fft = stat.fft_ang(nharm=1, imaginary=True)
>>> # S2 shape:  (..., L)       -> (..., 3)      [DC, cos, sin]
>>> # S3 shape:  (..., L, L)    -> (..., L, 3)   [Δl axis kept, l1 compressed]
>>> # S4 shape:  (..., L, L, L) -> (..., L, L, 3)
>>>
>>> # Rotation-invariant first-harmonic amplitude for S2:
>>> import numpy as np
>>> A1_S2 = np.sqrt(stat_fft.S2[..., 1]**2 + stat_fft.S2[..., 2]**2)
>>>
>>> # Angular modulation of S3 at a given relative angle Δl:
>>> A1_S3_delta0 = np.sqrt(stat_fft.S3[..., 0, 1]**2 + stat_fft.S3[..., 0, 2]**2)

fft_ang_sigma(nharm=1, imaginary=False)[source]#

Propagate per-coefficient standard deviations through fft_ang.

When the input object contains standard deviations sigma_j (e.g. the variance output of eval() with calc_var=True), applying fft_ang() directly would be wrong because linear projection of sigma gives near-zero values for harmonic components (cosine/sine sum to zero over uniform sigma values).

The correct formula for a linear map \(y_k = \sum_j A_{jk}\,x_j\) with independent inputs of variance \(\sigma_j^2\) is:

\[\sigma_{y_k} = \sqrt{\sum_j A_{jk}^2 \, \sigma_j^2}\]

This method applies that formula element-wise using the squared projection matrices, ensuring that the output sigma is always positive and physically meaningful for use as the sigma argument of reduce_distance().

Parameters:

nharm (int, optional) – Passed through to the underlying projection matrices (same as in fft_ang()).
imaginary (bool, optional) – Passed through (same as in fft_ang()).

Returns:

scat_cov – A new statistics object with the same shape as fft_ang output but containing propagated standard deviations.

iso_std(repeat=False)[source]#

get_nscale()[source]#

get_norient()[source]#

add_data_from_log_slope(y, n, ds=3)[source]#

add_data_from_slope(y, n, ds=3)[source]#

up_grade(nscale, ds=3)[source]#

class foscat.scat_cov.funct(NORIENT=4, LAMBDA=1.2, KERNELSZ=3, slope=1.0, all_type='float32', nstep_max=20, padding='SAME', gpupos=0, mask_thres=None, mask_norm=False, isMPI=False, TEMPLATE_PATH=None, BACKEND='torch', use_2D=False, use_1D=False, return_data=False, DODIV=False, use_median=False, InitWave=None, silent=True, mpi_size=1, mpi_rank=0)[source]#

Bases: foscat.FoCUS.FoCUS

fill(im, nullval=hp.UNSEEN)[source]#

moments(list_scat)[source]#

stat_cfft(im, image2=None, upscale=False, smooth_scale=0, spin=0)[source]#

stat_cfft_cell_ids(im, image2=None, upscale=False, smooth_scale=0, spin=0, cell_ids=None)[source]#: Full-sky stat_cfft + per-scale slicing by cell_ids. If cell_ids is None, identical to stat_cfft. If provided, cell_ids must be NEST indices at the input resolution.

div_norm(complex_value, float_value)[source]#

eval(image1, image2=None, mask=None, norm=None, calc_var=False, cmat=None, cmat2=None, Jmax=None, out_nside=None, edge=True, nside=None, cell_ids=None, spin=0)[source]#

Calculates the scattering correlations for a batch of images. Mean are done over pixels. mean of modulus:

S1 = <|I * Psi_j3|>

Normalization : take the log

power spectrum:: S2 = <|I * Psi_j3|^2>

Normalization : take the log
orig. x modulus:: S3 = < (I * Psi)_j3 x (|I * Psi_j2| * Psi_j3)^* >

Normalization : divide by (S2_j2 * S2_j3)^0.5
modulus x modulus:: S4 = <(|I * psi1| * psi3)(|I * psi2| * psi3)^*>

Normalization : divide by (S2_j1 * S2_j2)^0.5

Parameters:

image1 (tensor) – Image on which we compute the scattering coefficients [Nbatch, Npix, 1, 1]
image2 (tensor) – Second image. If not None, we compute cross-scattering covariance coefficients.
mask
norm (None or str) – If None no normalization is applied, if ‘auto’ normalize by the reference S2, if ‘self’ normalize by the current S2.
spin (Integer) – If different from 0 compute spinned data (U,V to Divergence/Rotational spin==1) or (Q,U to E,B spin=2). This implies that the input data is 2*12*nside^2.

Returns:

S1, S2, S3, S4 normalized

clean_norm()[source]#

computer_filter(M, N, J, L)[source]#: This function is strongly inspire by the package SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel.

cut_high_k_off(data_f, dx, dy)[source]#: This function is strongly inspire by the package SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel.

get_dxdy(j, M, N)[source]#: This function is strongly inspire by the package SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel.

get_edge_masks(M, N, J, d0=1, in_mask=None, edge_dx=None, edge_dy=None)[source]#: This function is strongly inspire by the package SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel.

scattering_cov(data, data2=None, Jmax=None, if_large_batch=False, S4_criteria=None, use_ref=False, normalization='S2', edge=False, in_mask=None, pseudo_coef=1, get_variance=False, ref_sigma=None, iso_ang=False, return_table=False, fft_ang=False, fft_nharm=1, fft_imaginary=True)[source]#

Calculates the scattering correlations for a batch of images, including:

This function is strongly inspire by the package SihaoCheng/scattering_transform Done by Sihao Cheng and Rudy Morel.

orig. x orig.:: P00 = <(I * psi)(I * psi)*> = L2(I * psi)^2
orig. x modulus:: C01 = <(I * psi2)(|I * psi1| * psi2)*> / factor

when normalization == ‘P00’, factor = L2(I * psi2) * L2(I * psi1) when normalization == ‘P11’, factor = L2(I * psi2) * L2(|I * psi1| * psi2)
modulus x modulus:: C11_pre_norm = <(|I * psi1| * psi3)(|I * psi2| * psi3)> C11 = C11_pre_norm / factor

when normalization == ‘P00’, factor = L2(I * psi1) * L2(I * psi2) when normalization == ‘P11’, factor = L2(|I * psi1| * psi3) * L2(|I * psi2| * psi3)
modulus x modulus (auto):: P11 = <(|I * psi1| * psi2)(|I * psi1| * psi2)*>

Parameters:

data (numpy array or torch tensor) – image set, with size [N_image, x-sidelength, y-sidelength]
if_large_batch (Bool (=False)) – It is recommended to use “False” unless one meets a memory issue
C11_criteria (str or None (=None)) – Only C11 coefficients that satisfy this criteria will be computed. Any expressions of j1, j2, and j3 that can be evaluated as a Bool is accepted.The default “None” corresponds to “j1 <= j2 <= j3”.
use_ref (Bool (=False)) – When normalizing, whether or not to use the normalization factor computed from a reference field. For just computing the statistics, the default is False. However, for synthesis, set it to “True” will stablize the optimization process.
normalization (str 'P00' or 'P11' (=``’P00’``)) – Whether ‘P00’ or ‘P11’ is used as the normalization factor for C01 and C11.
remove_edge (Bool (=False)) – If true, the edge region with a width of rougly the size of the largest wavelet involved is excluded when taking the global average to obtain the scattering coefficients.

Returns:

‘P00’ (torch tensor with size [N_image, J, L] (# image, j1, l1)) – the power in each wavelet bands (the orig. x orig. term)
’S1’ (torch tensor with size [N_image, J, L] (# image, j1, l1)) – the 1st-order scattering coefficients, i.e., the mean of wavelet modulus fields
’C01’ (torch tensor with size [N_image, J, J, L, L] (# image, j1, j2, l1, l2)) – the orig. x modulus terms. Elements with j1 < j2 are all set to np.nan and not computed.
’C11’ (torch tensor with size [N_image, J, J, J, L, L, L] (# image, j1, j2, j3, l1, l2, l3)) – the modulus x modulus terms. Elements not satisfying j1 <= j2 <= j3 and the conditions defined in ‘C11_criteria’ are all set to np.nan and not computed.
’C11_pre_norm’ and ‘C11_pre_norm_iso’ (pre-normalized modulus x modulus terms.)
’P11’ (torch tensor with size [N_image, J, J, L, L] (# image, j1, j2, l1, l2)) – the modulus x modulus terms with the two wavelets within modulus the same. Elements not following j1 <= j3 are set to np.nan and not computed.
’P11_iso’ (torch tensor with size [N_image, J, J, L] (# image, j1, j2, l2-l1)) – ‘P11’ averaged over l1 while keeping l2-l1 constant.

purge_edge_mask()[source]#

to_gaussian(x, in_mask=None)[source]#

from_gaussian(x)[source]#

square(x)[source]#

sqrt(x)[source]#

reduce_mean(x)[source]#

reduce_mean_batch(x)[source]#

reduce_sum_batch(x)[source]#

reduce_distance(x, y, sigma=None)[source]#

reduce_sum(x)[source]#

ldiff(sig, x)[source]#

log(x)[source]#

eval_comp_fast(image1, image2=None, mask=None, norm=None, cmat=None, cmat2=None)[source]#

eval_fast(image1, image2=None, mask=None, norm=None, cmat=None, cmat2=None)[source]#

calc_matrix_orientation(noise_map, image2=None)[source]#

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)[source]#

Synthesise a new field whose scattering-covariance statistics match those of image_target.

This is the main high-level entry point for texture synthesis, denoising, inpainting, and component separation with FOSCAT. Internally it runs L-BFGS-B (scipy.optimize.fmin_l_bfgs_b) on the scattering-covariance loss, with an optional coarse-to-fine multi-resolution schedule controlled by nstep.

The optimisation minimises:

\[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 \(\Phi\) is the scattering-covariance operator, \(d\) is image_target, and \(\sigma_k^2\) is the per-coefficient variance of \(\Phi(d)\) (used only when use_variance=True).

Parameters:

image_target (array-like) – Reference field whose statistics are to be reproduced.
- HEALPix: shape (npix,) or (B, npix) with npix = 12 * nside**2.
- 2D: shape (H, W) or (B, H, W).
- 1D: shape (N,) or (B, N).
Normalise the input before calling (zero mean, unit variance is standard):
```
xnorm = (image - image.mean()) / image.std()
```
reference (array-like, optional) – Second reference field for cross-covariance synthesis. When provided the loss matches the cross-statistics \(\Phi_\times(u, d_2)\) against the target \(\Phi_\times(d_1, d_2)\), where d1 = image_target and d2 = reference. Both fields must have the same shape. Useful for component separation (e.g. CMB × dust template).
nstep (int, optional) – Number of resolution levels in the coarse-to-fine schedule. Default is 4.

The algorithm downsamples image_target by successive factors of 2 to build a resolution pyramid, then optimises from coarsest to finest, using each result as the warm start for the next level. This dramatically improves convergence for large maps.

Capped automatically when nstep > jmax - 1 (map too small for that many downsampling steps).
seed (int, optional) – Random seed for the Gaussian white-noise initial condition used at the coarsest resolution level. Default is 1234. Has no effect when input_image is provided. Change to generate independent realisations:
```
results = [scat_op.synthesis(xnorm, seed=s) for s in range(8)]
```
Jmax (int or None, optional) – Maximum wavelet scale index included in the loss. None (default) uses all scales available for the map size. Decremented by 1 at each coarser resolution level during the multi-resolution schedule. Reduce to constrain only small-scale statistics:
```
result = scat_op.synthesis(xnorm, Jmax=4)  # ignore large scales
```
edge (bool, optional) – Enable edge-aware boundary handling for non-periodic maps. Default is False (periodic/full-sphere boundaries). Set to True for rectangular images or partial-sky patches. Activated automatically when in_mask is provided.
to_gaussian (bool, optional) – Gaussianise image_target before computing target statistics, then invert at the end. Default is False. Useful for highly non-Gaussian fields (log-normal distributions, etc.) to decouple the histogram constraint from the scattering statistics.
use_variance (bool, optional) – Weight each loss term by the inverse variance of the corresponding coefficient in image_target. Default is True.
- True: loss is scale-invariant — large-scale and small-scale coefficients contribute equally regardless of amplitude.
- False: uniform weighting, dominated by the largest-amplitude statistics (usually coarse scales).
synthesised_N (int, optional) – Number of independent synthetic maps to produce simultaneously. Default is 1. All synthesised_N maps share the same loss, which amortises the per-iteration GPU cost. The output gains a leading batch dimension:
```
result = scat_op.synthesis(xnorm, synthesised_N=4)
# result.shape == (4, H, W)  for 2D
```
input_image (array-like or None, optional) – Warm-start initial condition. When None (default), starts from Gaussian white noise. When provided, input_image is downsampled to each resolution level and used as the initial guess at the coarsest level. Useful for iterative refinement or for injecting a prior:
```
result_v1 = scat_op.synthesis(xnorm, NUM_EPOCHS=100)
result_v2 = scat_op.synthesis(xnorm, NUM_EPOCHS=500,
                               input_image=result_v1)
```
grd_mask (array-like or None, optional) – Binary mask selecting which pixels the optimiser is allowed to update. Pixels where grd_mask = 0 are frozen to their values in image_target; only pixels where grd_mask = 1 move. Same shape as image_target. Downsampled at each resolution level. Useful for inpainting (reconstruct missing regions while keeping observed pixels fixed):
```
grd_mask = np.zeros_like(xnorm)
grd_mask[missing_pixels] = 1.0
result = scat_op.synthesis(xnorm, grd_mask=grd_mask)
```
in_mask (array-like or None, optional) – Binary mask marking invalid pixels in the input data that should not contribute to the reference statistics. Same shape as image_target. Pixels with in_mask = 0 are excluded from the scattering-covariance computation at every scale level (the mask is downsampled along with the map). Setting in_mask also enables edge=True internally. Useful for partial-sky CMB maps or images with missing regions:
```
mask = np.ones_like(xnorm)
mask[bad_pixels] = 0.0
result = scat_op.synthesis(xnorm, in_mask=mask)
```
iso_ang (bool, optional) – Use isotropically averaged statistics in the loss. Default is False.
- False: the full oriented S1–S4 statistics are used. The synthesised field can reproduce anisotropic structures (filaments, oriented textures).
- True: statistics are collapsed to their rotationally invariant content via iso_mean() before computing the loss. Reduces the number of constraints and accelerates convergence for statistically isotropic fields (e.g. CMB).
Angular-reduction summary with iso_ang=True:

Statistic

Input shape

Output shape

S1, S2

[…, L]

[…] (mean over L)

S3, S3P

[…, L, L]

[…, L] (by Δl = l₂−l₁)

S4

[…, L, L, L]

[…, L, L]

Note

iso_ang and fft_ang should not be used together. iso_ang is the harder reduction (mean only); fft_ang is softer and keeps angular variation.
fft_ang (bool, optional) – Use Fourier-compressed angular statistics in the loss. Default is False. Softer alternative to iso_ang: instead of collapsing each orientation axis to a single mean, keeps the first fft_nharm Fourier harmonics along each axis.

With fft_nharm=1, fft_imaginary=True (defaults), each orientation axis L is projected to 3 coefficients:
- index 0 — DC (mean, same as iso_ang);
- index 1 — cosine of first harmonic;
- index 2 — sine of first harmonic.
The amplitude \(A_1 = \sqrt{c_1^2 + s_1^2}\) is rotation-invariant: it does not depend on the absolute orientation of the image.

Shape reduction:

Statistic

Input shape

Output shape (nharm=1)

S1, S2

[…, L]

[…, 3]

S3, S3P

[…, L, L]

[…, 3, 3]

S4

[…, L, L, L]

[…, 3, 3, 3]

For S3/S4 the projection is the tensor product of independent 1-D Fourier projections on each orientation axis:
```
result = scat_op.synthesis(xnorm, fft_ang=True, NUM_EPOCHS=300)
```
fft_nharm (int, optional) – Number of Fourier harmonics to keep beyond the DC term when fft_ang=True. Default is 1. The number of output coefficients per orientation axis is 1 + 2*fft_nharm (with fft_imaginary=True) or 1 + fft_nharm (with fft_imaginary=False).
fft_imaginary (bool, optional) – Whether to include both cosine and sine components when fft_ang=True. Default is True.
- True (recommended): output per axis = 1 + 2*fft_nharm coefficients [DC, cos₁, sin₁, cos₂, sin₂, …]. The harmonic amplitude \(\sqrt{c_k^2 + s_k^2}\) is independent of the image orientation.
- False: output per axis = 1 + fft_nharm coefficients [DC, cos₁, cos₂, …]. A field oriented at 90° (zero of cosine) gives a near-zero first-harmonic coefficient even if strongly anisotropic — use only if orientation is fixed.
EVAL_FREQUENCY (int, optional) – Print the current loss every N L-BFGS-B iterations. Default is 100.
NUM_EPOCHS (int, optional) – Maximum number of L-BFGS-B iterations per resolution level. Default is 300. The optimiser may stop earlier when its convergence criterion is met. Total wall-clock time scales as nstep × NUM_EPOCHS.
scat_cov_method (str, optional) – Internal method for computing scattering covariances. Default is 'eval' (recommended), which uses eval() with norm='auto' and caches the normalisation after the first call. Any other value falls back to the legacy scattering_cov() path (2D only, kept for backward compatibility).
n_up (int, optional) – Number of extra upsampling steps beyond the target size, keeping the same Jmax. Default is 0 (no extra steps).

With n_up=1, after completing synthesis at N × N, the algorithm continues at 2N × 2N using the same scattering statistics as the N × N target. The result is a field that locally matches the statistics of the original target, embedded in a larger canvas. Only available for 2D maps.

The Jmax used during n_up steps is pinned to the value effective for the N × N target so that wavelet filters and the norm cache remain consistent:
```
result = scat_op.synthesis(xnorm_256, nstep=3, n_up=1)
# result.shape == (512, 512)
```

Returns:

numpy.ndarray – Synthesised field.

Shape equals image_target.shape when synthesised_N=1 and n_up=0.
When synthesised_N > 1, a leading dimension of size synthesised_N is added.
When n_up > 0 (2D only), spatial dimensions are multiplied by 2**n_up.

Examples

Minimal 2D synthesis:

import foscat.scat_cov2D as sc
import numpy as np

scat_op = sc.funct(NORIENT=4)
xnorm = (image - image.mean()) / image.std()
result = scat_op.synthesis(xnorm, seed=10, nstep=3, NUM_EPOCHS=300)

Batch synthesis with mask:

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)

Inpainting — reconstruct missing pixels, freeze observed ones:

grd_mask = np.zeros_like(xnorm)
grd_mask[hole_pixels] = 1.0
result = scat_op.synthesis(
    xnorm, grd_mask=grd_mask, nstep=3, edge=True, NUM_EPOCHS=500
)

Upsampled synthesis (n_up):

result = scat_op.synthesis(xnorm_256, nstep=3, n_up=1, NUM_EPOCHS=300)
# result.shape == (512, 512)

Cross-covariance — component separation:

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,
)

Notes

The optimiser is scipy.optimize.fmin_l_bfgs_b, a quasi-Newton method that uses a low-memory approximation of the inverse Hessian. It converges in tens to a few hundred iterations for typical scattering-covariance losses, far fewer than first-order methods.

For large maps (nside ≥ 256 for HEALPix, H/W ≥ 256 for 2D) always use nstep ≥ 3. Starting directly at full resolution wastes compute and converges poorly.

to_numpy(x)[source]#

Statistic	Input shape	Output shape
S1, S2	`[…, L]`	`[…]` (mean over L)
S3, S3P	`[…, L, L]`	`[…, L]` (by Δl = l₂−l₁)
S4	`[…, L, L, L]`	`[…, L, L]`

Statistic	Input shape	Output shape (nharm=1)
S1, S2	`[…, L]`	`[…, 3]`
S3, S3P	`[…, L, L]`	`[…, 3, 3]`
S4	`[…, L, L, L]`	`[…, 3, 3, 3]`