12 — Minkowski Tensors

Purpose: Quantify morphological anisotropy in CIB and tSZ patches using rank-2 Minkowski tensors — the tensorial generalisation of the scalar Minkowski functionals computed in docs/05_plots.ipynb.

For each intensity threshold ν, the map is binarised to the excursion set K = {x : T(x) > ν} and a rank-2 tensor W is computed from the geometry of K. Eigendecomposing W yields:

  • β = λ_min / λ_max ∈ [0, 1]: anisotropy index (1 = isotropic, 0 = maximally elongated)

  • θ ∈ (-π/2, π/2]: orientation of the major axis

Three tensor types are computed to expose their different sensitivities:

  • W012 (W^{0,2}_1): boundary normal tensor — sums n⊗n over boundary pixels. Most sensitive to the shape of cluster boundaries. Recommended default.

  • W200 (W^{2,0}_0): area inertia tensor — sums r⊗r over interior pixels. Measures elongation of filled regions.

  • W201 (W^{2,0}_1): boundary position tensor — sums r⊗r over boundary pixels. Hybrid between W012 and W200.

Comparing β(ν) and θ(ν) between Agora, DDPM, and Gaussian baseline tests whether the diffusion model reproduces the morphological anisotropy of the true foreground fields — a property that scalar statistics (power spectra, moments, MFs) are blind to.

Sanity check: A statistically isotropic Gaussian random field should give β ≈ 1 on average and θ uniformly distributed — deviations indicate genuine non-Gaussian anisotropy.

Inputs:

  • Test maps: data/low_pass/2mJy/*.npy

  • DDPM samples: data/low_pass/2mJy/new_samples_*.npy

  • Norm params: data/low_pass/2mJy/norm_params_2mJy.npy

Outputs:

  • plots/minkowski_tensors_beta.pdf — β(ν) curves, all tensor types, CIB and tSZ

  • plots/minkowski_tensors_theta.pdf — θ polar histograms at fixed thresholds

  • plots/minkowski_tensors_residuals.pdf — (β_Agora − β_DDPM) / σ_Agora for W012

Key module functions:

  • foregrounds_diffusion.statistics.compute_minkowski_tensors

  • foregrounds_diffusion.statistics.MINKOWSKI_TENSOR_DESCRIPTIONS

Reference: Schroder-Turk et al. (2013), New J. Phys. 15 083028

Load modules and set configuration.

THRESHOLDS is in the normalised [0, 1] range (after apply_maxmin_normalization);

avoid ν = 0 and ν = 1 where the excursion set is trivially full or empty.

1 Configuration and tensor descriptions

Three tensor types are computed at 30 intensity thresholds spanning ν ∈ [0.05, 0.95] of the normalised range:

Tensor

Symbol

Physical interpretation

W012

W^{0,2}_1

Boundary normal tensor: sums n⊗n over boundary pixels (Sobel normals). Most sensitive to the shape of cluster boundaries.

W200

W^{2,0}_0

Area inertia tensor: sums r⊗r over interior pixels. Measures elongation of the filled excursion region.

W201

W^{2,0}_1

Boundary position tensor: sums r⊗r over boundary pixels. Hybrid between W012 and W200.

The eigendecomposition W = λ₁v₁v₁ᵀ + λ₂v₂v₂ᵀ (λ₁ ≤ λ₂) yields:

  • β = λ₁/λ₂ ∈ [0, 1]: anisotropy index. β = 1 means the excursion set is statistically isotropic; β → 0 means it is highly elongated.

  • θ ∈ (−π/2, π/2]: orientation of the major axis.

[ ]:

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

from foregrounds_diffusion.preprocessing import apply_maxmin_normalization
from foregrounds_diffusion.morphology import (
    compute_minkowski_tensors,
    MINKOWSKI_TENSOR_DESCRIPTIONS,
)

PTSRC            = 2
PIXEL_RES_ARCMIN = 1.40625
N_MAPS           = 100     # maps per source; reduce for quick tests

# 30 thresholds in the normalised [0, 1] range, avoiding the extreme tails
# where excursion sets are near-empty or near-full
THRESHOLDS = np.linspace(0.05, 0.95, 30)

TENSOR_TYPES = ('W012', 'W200', 'W201')

PATCHES_DIR  = Path(f'data/low_pass/{PTSRC}mJy')
PROJECT_ROOT = Path('/home/apb86/cmb_foregrounds_diffusion')

print('Tensor descriptions:')
for tt, desc in MINKOWSKI_TENSOR_DESCRIPTIONS.items():
    print(f'  {tt}: {desc}')

Load normalised patches. The Minkowski tensor code re-normalises internally

via norm_fn, so the absolute pixel values here are less critical than for

power spectra. Using the standard z-score / min-max loading ensures

consistency with the train/test split.

2 Load and denormalise maps

The same loading pattern as notebooks 06–11. For the Minkowski tensor computation the maps are re-normalised to [0, 1] by apply_maxmin_normalization inside compute_minkowski_tensors, regardless of the input units. This ensures that the threshold parameter ν is dimensionless and comparable across maps with very different amplitudes (e.g. tSZ clusters of different masses).

[ ]:

# ---------------------------------------------------------------------------
# Load and denormalise maps  (identical pattern to notebooks 10 and 11)
# ---------------------------------------------------------------------------
norm_params = np.load(PATCHES_DIR / f'norm_params_{PTSRC}mJy.npy')
cib_mean, cib_std, tsz_mean, tsz_std = norm_params
print(f'CIB — mean: {cib_mean:.4f} µK  std: {cib_std:.4f} µK')
print(f'tSZ — mean: {tsz_mean:.4f} µK  std: {tsz_std:.4f} µK')

cib_maps = np.load(PATCHES_DIR / f'CIB_map_150GHz_256_st6_zscore_{PTSRC}mJy_lp.npy')
tsz_maps = np.load(PATCHES_DIR / f'tSZ3_map_150GHz_256_st6_zscore_{PTSRC}mJy_lp.npy')

rng = np.random.default_rng(seed=42)
indices  = rng.permutation(len(cib_maps))
test_idx = indices[int(0.8 * len(cib_maps)):]

cib_test = cib_maps[test_idx, :, :, 0] * cib_std + cib_mean
tsz_test = tsz_maps[test_idx, :, :, 0] * tsz_std + tsz_mean
print(f'Test set: {len(cib_test)} patches')

DDPM samples shape: (N, 2, H, W) channels-first.

Index [:, 0] = CIB channel, [:, 1] = tSZ channel.

Denormalise: CIB = raw * (cib_max - cib_min) + cib_min;

tSZ = raw * tsz_std + tsz_mean.

3 Load DDPM samples and Gaussian baseline

[ ]:

ddpm_raw = np.load(
    PROJECT_ROOT / 'data' / 'low_pass' / f'{PTSRC}mJy' /
    f'new_samples_cib_tsz_{PTSRC}mJy_zero_norm_6x6_w_au_lp.npy'
)  # (N, 2, H, W)
ddpm_cib = ddpm_raw[:, 0] * cib_std + cib_mean
ddpm_tsz = ddpm_raw[:, 1] * tsz_std + tsz_mean
print(f'DDPM samples: {len(ddpm_cib)} patches')

gauss_maps = np.load(PATCHES_DIR / f'gaussian_cib_tsz_{PTSRC}mJy_lp.npy')
gauss_cib  = gauss_maps[:, 0] if gauss_maps.shape[1] == 2 else gauss_maps[:, :, :, 0]
gauss_tsz  = gauss_maps[:, 1] if gauss_maps.shape[1] == 2 else gauss_maps[:, :, :, 1]
print(f'Gaussian samples: {len(gauss_cib)} patches')

N = min(N_MAPS, len(cib_test), len(ddpm_cib), len(gauss_cib))
print(f'Using {N} maps per source')

compute_minkowski_tensors(maps_nhw, norm_fn, thresholds, tensor_types, centred)

norm_fn : normalisation applied per-map before thresholding

(use apply_maxmin_normalization so ν ∈ [0,1] is meaningful)

centred : if True, subtract excursion-set centroid from r⊗r tensors

(affects W200, W201 only; default True)

n_jobs : pass n_jobs=-1 to use all CPU cores via joblib

4 Compute Minkowski tensors

compute_minkowski_tensors loops over N maps and T thresholds. For each (map, threshold) pair:

  1. Apply norm_fn (here apply_maxmin_normalization) to map the patch to [0,1].

  2. Binarise: binary = normalised_map > ν.

  3. Skip trivially empty (< 2 interior pixels) or trivially full (< 2 exterior pixels) excursion sets — these are uninformative and return β = 1, θ = 0.

  4. Compute the requested tensor(s) from the binary mask.

  5. Eigendecompose to get β, θ.

Returns a dict:

{
  'W012': {'beta': ndarray(N, T), 'theta': ndarray(N, T)},
  'W200': {'beta': ...,           'theta': ...},
  'W201': {'beta': ...,           'theta': ...},
}

Runtime: ≈ 5–15 min for 100 maps × 30 thresholds × 3 tensor types on CPU. Use n_jobs=-1 to parallelise over maps with joblib.

[ ]:

# ---------------------------------------------------------------------------
# Compute Minkowski tensors
# NOTE: ~5-15 min for 100 maps × 30 thresholds × 3 tensor types on CPU.
#       Reduce N_MAPS or THRESHOLDS for a quick test.
# ---------------------------------------------------------------------------
norm_fn = apply_maxmin_normalization

print('Computing Agora CIB tensors...')
mt_agora_cib = compute_minkowski_tensors(
    cib_test[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Computing Agora tSZ tensors...')
mt_agora_tsz = compute_minkowski_tensors(
    tsz_test[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Computing DDPM CIB tensors...')
mt_ddpm_cib = compute_minkowski_tensors(
    ddpm_cib[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Computing DDPM tSZ tensors...')
mt_ddpm_tsz = compute_minkowski_tensors(
    ddpm_tsz[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Computing Gaussian CIB tensors...')
mt_gauss_cib = compute_minkowski_tensors(
    gauss_cib[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Computing Gaussian tSZ tensors...')
mt_gauss_tsz = compute_minkowski_tensors(
    gauss_tsz[:N], norm_fn, THRESHOLDS, tensor_types=TENSOR_TYPES)

print('Done.')

5 Anisotropy index β(ν): comparing CIB and tSZ morphology

Plot β(ν) — the mean anisotropy index as a function of threshold ν — in a 2-row (CIB, tSZ) × 3-column (W012, W200, W201) layout. Error bands show the standard deviation across the N test maps.

A statistically isotropic random field (e.g. the Gaussian baseline) has β ≈ 1 at all thresholds by symmetry. Non-isotropic structure — elongated tSZ filaments, asymmetric CIB source clusters — drives β < 1. The DDPM should reproduce the Agora β(ν) curve if it has learnt the correct cluster morphology; the Gaussian baseline provides a reference for how much anisotropy is generated purely by the patch geometry.

[ ]:

# ---------------------------------------------------------------------------
# Plot β(ν): anisotropy index vs threshold
# Layout: 2 rows (CIB, tSZ) × 3 cols (W012, W200, W201)
# ---------------------------------------------------------------------------
sources_cib = [
    ('Agora',    mt_agora_cib, 'C0', '-'),
    ('DDPM',     mt_ddpm_cib,  'C1', '-'),
    ('Gaussian', mt_gauss_cib, 'C2', '--'),
]
sources_tsz = [
    ('Agora',    mt_agora_tsz, 'C0', '-'),
    ('DDPM',     mt_ddpm_tsz,  'C1', '-'),
    ('Gaussian', mt_gauss_tsz, 'C2', '--'),
]

tensor_titles = {
    'W012': r'$W^{0,2}_1$ — boundary normals',
    'W200': r'$W^{2,0}_0$ — area inertia',
    'W201': r'$W^{2,0}_1$ — boundary positions',
}

fig, axes = plt.subplots(2, 3, figsize=(15, 8), sharex=True, sharey=True)

for row, (channel, sources) in enumerate([('CIB', sources_cib), ('tSZ', sources_tsz)]):
    for col, tt in enumerate(TENSOR_TYPES):
        ax = axes[row, col]
        for label, mt, color, ls in sources:
            beta = mt[tt]['beta']          # (N, T)
            mean = beta.mean(axis=0)
            std  = beta.std(axis=0)
            ax.plot(THRESHOLDS, mean, color=color, ls=ls, lw=1.5, label=label)
            ax.fill_between(THRESHOLDS, mean - std, mean + std,
                            alpha=0.2, color=color)
        ax.axhline(1, color='k', lw=0.6, ls=':', alpha=0.5)
        ax.set_ylim(0, 1.05)
        ax.set_xlim(THRESHOLDS[0], THRESHOLDS[-1])
        ax.grid(True, ls=':', lw=0.4, alpha=0.5)
        if row == 0:
            ax.set_title(tensor_titles[tt], fontsize=12)
        if col == 0:
            ax.set_ylabel(fr'{channel} — $\beta = \lambda_{{\min}}/\lambda_{{\max}}$',
                          fontsize=11)
        if row == 1:
            ax.set_xlabel(r'Threshold $\nu$', fontsize=11)
        if row == 0 and col == 2:
            ax.legend(fontsize=10)

plt.suptitle('Minkowski tensor anisotropy index β(ν)\n'
             r'(β = 1 = isotropic; β → 0 = maximally elongated)',
             fontsize=13)
plt.tight_layout()
Path('plots').mkdir(exist_ok=True)
plt.savefig('plots/minkowski_tensors_beta.pdf', dpi=200, bbox_inches='tight')
plt.show()

6 Orientation distribution θ at representative thresholds

Polar histograms of the major-axis orientation θ at three representative excursion levels (low ν = 0.2, mid ν = 0.5, high ν = 0.8). A uniform distribution around the circle indicates statistical isotropy. Any preferred direction (e.g. a peak near 0° or 90°) would indicate a systematic alignment of cluster features — which should not appear in either the AGORA simulations or a well-trained DDPM unless the simulation has an inherent preferred axis (e.g. from the scanning strategy or the flat-sky projection).

[ ]:

# ---------------------------------------------------------------------------
# Plot θ distributions: polar histograms at three representative thresholds
# A statistically isotropic field should give a uniform distribution.
# Deviations indicate preferred orientations in the morphology.
# ---------------------------------------------------------------------------
PROBE_THRESHOLDS = [0.2, 0.5, 0.8]   # low / mid / high excursion level
probe_idx = [int(np.argmin(np.abs(THRESHOLDS - t))) for t in PROBE_THRESHOLDS]

fig, axes = plt.subplots(
    2, len(PROBE_THRESHOLDS),
    figsize=(12, 7),
    subplot_kw={'projection': 'polar'},
)

theta_bins = np.linspace(-np.pi / 2, np.pi / 2, 19)  # 18 bins × 10°
bin_centres = 0.5 * (theta_bins[:-1] + theta_bins[1:])
bin_width   = theta_bins[1] - theta_bins[0]

for row, (channel, mt_agora, mt_ddpm, mt_gauss) in enumerate([
    ('CIB', mt_agora_cib, mt_ddpm_cib, mt_gauss_cib),
    ('tSZ', mt_agora_tsz, mt_ddpm_tsz, mt_gauss_tsz),
]):
    for col, (t_val, t_idx) in enumerate(zip(PROBE_THRESHOLDS, probe_idx)):
        ax = axes[row, col]
        for label, mt, color in [
            ('Agora',    mt_agora, 'C0'),
            ('DDPM',     mt_ddpm,  'C1'),
            ('Gaussian', mt_gauss, 'C2'),
        ]:
            theta_vals = mt['W012']['theta'][:, t_idx]  # (N,)
            counts, _ = np.histogram(theta_vals, bins=theta_bins)
            counts = counts / counts.sum()              # normalise to density
            ax.bar(bin_centres, counts, width=bin_width * 0.85,
                   alpha=0.5, color=color, label=label)
        # Uniform reference
        ax.axhline(1 / len(bin_centres), color='k', lw=0.8, ls='--', alpha=0.6)
        ax.set_theta_zero_location('N')
        ax.set_theta_direction(-1)
        ax.set_thetalim(-np.pi / 2, np.pi / 2)
        ax.set_title(f'{channel}  ν = {t_val:.1f}', fontsize=10, pad=8)
        if row == 0 and col == 0:
            ax.legend(fontsize=9, loc='upper right',
                      bbox_to_anchor=(1.5, 1.15))

plt.suptitle(r'Orientation θ of major axis (W012, $\nu$ = threshold)'
             '\nUniform = isotropic; peaked = preferred orientation',
             fontsize=12)
plt.tight_layout()
plt.savefig('plots/minkowski_tensors_theta.pdf', dpi=200, bbox_inches='tight')
plt.show()

7 Standardised residuals: (β_Agora − β_DDPM) / σ_Agora

Residuals for the W012 tensor (most physically interpretable) show whether any systematic anisotropy bias persists across the threshold range, in units of the Agora sample variance. Large residuals at high thresholds (rare excursion sets) may indicate that the DDPM does not accurately reproduce the statistics of the most extreme cluster configurations.

[ ]:

# ---------------------------------------------------------------------------
# Residuals: (β_Agora − β_DDPM) / σ_Agora  for W012
# Shows where the DDPM anisotropy deviates from Agora in units of Agora σ.
# ---------------------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)

for ax, channel, mt_agora, mt_ddpm in [
    (axes[0], 'CIB', mt_agora_cib, mt_ddpm_cib),
    (axes[1], 'tSZ', mt_agora_tsz, mt_ddpm_tsz),
]:
    beta_agora = mt_agora['W012']['beta']   # (N, T)
    beta_ddpm  = mt_ddpm['W012']['beta']
    sigma = beta_agora.std(axis=0) + 1e-10
    resid = (beta_agora.mean(axis=0) - beta_ddpm.mean(axis=0)) / sigma

    ax.plot(THRESHOLDS, resid, color='C0', lw=1.5)
    ax.fill_between(THRESHOLDS, resid, alpha=0.3, color='C0')
    ax.axhline(0,  color='k',    lw=0.8, ls='--')
    ax.axhline(1,  color='gray', lw=0.6, ls=':')
    ax.axhline(-1, color='gray', lw=0.6, ls=':')
    ax.set_xlabel(r'Threshold $\nu$', fontsize=11)
    ax.set_title(fr'{channel} — W012 $\beta$ residual', fontsize=12)
    ax.grid(True, ls=':', lw=0.4, alpha=0.5)

axes[0].set_ylabel(
    r'$(\bar{\beta}^\mathrm{Agora} - \bar{\beta}^\mathrm{DDPM})\,/\,\sigma^\mathrm{Agora}$',
    fontsize=11)

plt.suptitle('W012 anisotropy residuals (Agora − DDPM)', fontsize=13)
plt.tight_layout()
plt.savefig('plots/minkowski_tensors_residuals.pdf', dpi=200, bbox_inches='tight')
plt.show()

for channel, mt_agora, mt_ddpm in [
    ('CIB', mt_agora_cib, mt_ddpm_cib),
    ('tSZ', mt_agora_tsz, mt_ddpm_tsz),
]:
    sigma = mt_agora['W012']['beta'].std(axis=0) + 1e-10
    resid = (mt_agora['W012']['beta'].mean(0) - mt_ddpm['W012']['beta'].mean(0)) / sigma
    print(f'{channel}: {(np.abs(resid) < 1).mean():.0%} of thresholds within 1σ')