morphology — Minkowski functionals and tensors

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)