morphology — Minkowski functionals and tensors

morphology computes topological and geometrical statistics on excursion sets — the subsets of pixels that exceed a given intensity threshold ν. These statistics complement the power spectrum by probing non-Gaussian morphological features such as the shapes and connectivity of tSZ clusters.

Minkowski functionals

The three scalar functionals characterise excursion-set topology as a function of threshold ν (paper §4.1, Figure 6):

  • M0: area fraction — fraction of pixels above ν

  • M1: perimeter — total boundary length normalised by map area

  • M2: Euler characteristic — connected components minus holes

Requires pip install quantimpy:

from foregrounds_diffusion.morphology import compute_mfs
from foregrounds_diffusion.preprocessing import apply_maxmin_normalization

thresholds = np.linspace(0.05, 0.95, 19)
M0, M1, M2 = compute_mfs(maps_nhw, apply_maxmin_normalization, thresholds, n_jobs=-1)
# Each: (N, len(thresholds))

Minkowski tensors

compute_minkowski_tensors() is a tensorial extension that additionally captures anisotropy and orientation per excursion set. It returns the anisotropy index β = λ_min / λ_max ∈ [0, 1] (β = 1 is perfectly isotropic) and major-axis orientation θ for three tensor types (W012, W200, W201). W012 (boundary normal tensor) is recommended as the default:

from foregrounds_diffusion.morphology import compute_minkowski_tensors

results = compute_minkowski_tensors(
    maps_nhw, apply_maxmin_normalization, thresholds,
    tensor_types=["W012"], centred=True,
)
beta = results["W012"]["beta"]   # (N, len(thresholds))
theta = results["W012"]["theta"]

See 12 — Minkowski Tensors for β(ν) curves, polar orientation histograms, and DDPM/AGORA residuals.

API

Morphological statistics: Minkowski functionals and Minkowski tensors.

This module computes two families of morphological statistics on excursion sets — the subsets of pixels above a given intensity threshold ν.

Minkowski Functionals

The three scalar Minkowski functionals characterise the topology of excursion sets (paper §4.1, Figure 6):

  • M0 (area fraction): fraction of pixels above threshold ν.

  • M1 (perimeter): total boundary length normalised by map area.

  • M2 (Euler characteristic): number of connected components minus number of holes, measuring the connectivity of the excursion set.

Computed by compute_mfs() using the quantimpy library (Boelens & Tchelepi 2021). Requires pip install quantimpy.

Minkowski Tensors

Rank-2 tensorial generalisations of the functionals that additionally encode the anisotropy and orientation of morphological features. compute_minkowski_tensors() returns the anisotropy index β = λ_min / λ_max ∈ [0, 1] (β = 1 is perfectly isotropic) and the major-axis orientation θ per map per threshold, for three tensor types:

  • W012: boundary normal tensor via Sobel-estimated outward normals n⊗n. Probes isotropy of cluster boundary shapes. Recommended default.

  • W200: area inertia tensor r⊗r over interior pixels. Measures elongation of filled excursion regions.

  • W201: boundary position tensor r⊗r over boundary pixels.

See tutorial docs/tutorials/12_minkowski_tensors.ipynb for β(ν) curves and polar orientation histograms comparing AGORA and DDPM maps.

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

foregrounds_diffusion.morphology.compute_mfs(maps_nhw, norm_fn, thresholds, n_jobs=1)[source]

Compute Minkowski functionals M0, M1, M2 across a stack of maps.

Requires the quantimpy package (Boelens & Tchelepi 2021).

Parameters:

  • maps_nhw (ndarray, shape (N, H, W))

  • norm_fn (callable) – Normalisation function mapping a 2D array to [0, 1] (e.g. apply_maxmin_normalization()).

  • thresholds (array_like) – Intensity thresholds ν.

  • n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). -1 = all cores.

Returns:

M0, M1, M2 (ndarray, each shape (N, len(thresholds))) – M0 = area fraction, M1 = perimeter, M2 = Euler characteristic.

foregrounds_diffusion.morphology.compute_minkowski_tensors(maps_nhw, norm_fn, thresholds, tensor_types=('W012',), centred=True, n_jobs=1)[source]

Compute Minkowski tensor anisotropy indices across intensity thresholds.

For each map and threshold ν, the map is binarised to the excursion set K = {x : T(x) > ν}. The requested rank-2 Minkowski tensor(s) are computed from K, then eigendecomposed to give the anisotropy index β = λ_min/λ_max and the major-axis orientation θ.

Parameters:

  • maps_nhw (ndarray, shape (N, H, W)) – Stack of flat-sky patches.

  • norm_fn (callable) – Normalisation function applied to each 2D patch before thresholding (e.g. apply_maxmin_normalization). Pass lambda x: x to use raw pixel values.

  • thresholds (array_like) – Intensity thresholds ν in the normalised range.

  • tensor_types (tuple of str) – Which tensor(s) to compute. Any subset of:

    'W012'

    W^{0,2}_1 — interface normal tensor (default). Boundary normals n⊗n via Sobel gradients. Most sensitive to boundary shape.

    'W200'

    W^{2,0}_0 — area inertia tensor. Interior pixel positions r⊗r. Measures elongation of filled regions.

    'W201'

    W^{2,0}_1 — boundary position tensor. Boundary pixel positions r⊗r. Hybrid between W012 and W200.

  • centred (bool) – If True (default), subtract the excursion-set centroid from position vectors before forming the tensor. Affects W200 and W201 only.

  • n_jobs (int) – Number of parallel workers (joblib). 1 = serial (default). −1 = use all cores.

Returns:

results (dict) – Keyed by tensor type string. Each value is a dict:

'beta'ndarray, shape (N, len(thresholds))

Anisotropy index β = λ_min/λ_max ∈ [0, 1]. 1 = isotropic.

'theta'ndarray, shape (N, len(thresholds))

Major-axis orientation in radians ∈ (-π/2, π/2].

Notes

β = 1 and θ = 0 are returned for empty or full excursion sets (fewer than 2 interior pixels or fewer than 2 exterior pixels). As a sanity check, a Gaussian random field should give β ≈ 1 at all thresholds.

Examples

>>> thresholds = np.linspace(0.05, 0.95, 50)
>>> results = compute_minkowski_tensors(
...     maps, apply_maxmin_normalization, thresholds,
...     tensor_types=('W012', 'W200'))
>>> beta_W012 = results['W012']['beta']   # shape (N, 50)
>>> theta_W012 = results['W012']['theta'] # shape (N, 50)