PH Computer

The ph_computer module provides functions for computing persistent homology from images using cubical complexes.

Overview 

This module contains the core functions for computing persistent homology (PH) from image data. It uses the GUDHI library’s cubical complex implementation.

Key Features:

Sublevel and superlevel filtrations
T-construction and V-construction cubical complexes
Multi-dimensional homology (H0, H1, H2, …)
Automatic dimension detection
Infinite persistence handling

Basic Usage 

Quick Start 

from medtda.ph_computer import compute_ph
import numpy as np

# Create a simple 2D image
image = np.random.rand(50, 50)

# Compute persistence barcodes
barcodes = compute_ph(image, filtration_type='sublevel')

# barcodes is a dict: {'H0': H0 pairs, 'H1': H1 pairs, ...}
print(f"H0 features: {len(barcodes['H0'])}")
print(f"H1 features: {len(barcodes['H1'])}")

With Custom Parameters 

barcodes = compute_ph(
    image,
    filtration_type='superlevel',
    construction='V',
    max_dimension=2
)

Functions 

compute_ph 

medtda.ph_computer.compute_ph(image_array, maxdim=None, construction='T', filtration='sublevel', max_value=None)[source]

Compute persistent homology of an image array using cubical complexes.

Parameters:

image_array (numpy.ndarray) – 2D, 3D, or 4D image array.
maxdim (int, optional) – Maximum homology dimension to compute. If None (default), automatically determined from image dimensions: - 2D images: maxdim=1 (compute H0, H1) - 3D images: maxdim=2 (compute H0, H1, H2) - 4D images: maxdim=3 (compute H0, H1, H2, H3)
construction ({'T', 'V'}, default='T') – Type of cubical complex construction: - ‘T’: T-construction (default, pixels/voxels as top-cells, 8-neighborhood in 2D) - ‘V’: V-construction (pixels/voxels as 0-cells, 4-neighborhood in 2D)
filtration ({'sublevel', 'superlevel'}, default='sublevel') – Filtration type: - ‘sublevel’: Standard sublevel set filtration - ‘superlevel’: Superlevel set filtration (image is negated)
max_value (float, optional) – Maximum value for capping infinite persistence. If None, uses the maximum value in the image array.

Returns:

ph_capped (numpy.ndarray) – Cleaned and capped persistence array with columns [dimension, birth, death].
dgms (list of numpy.ndarray) – List of persistence diagrams in GUDHI format (one per dimension).

Raises:

ValueError – If construction is not ‘T’ or ‘V’, or filtration is not ‘sublevel’ or ‘superlevel’.

Notes

For superlevel filtration, the image is negated internally, and the max_value is also negated for proper capping.

Examples

>>> import numpy as np
>>> from medtda.ph_computer import compute_ph
>>> image = np.random.rand(50, 50)
>>> ph_array, diagrams = compute_ph(image, maxdim=1, construction='T')

Compute persistent homology from an image.

Signature:

def compute_ph(
    image: np.ndarray,
    filtration_type: str = 'sublevel',
    construction: str = 'T',
    max_dimension: int = -1
) -> Dict[int, np.ndarray]

Parameters:

image (np.ndarray) - Input image (2D or 3D)
filtration_type (str) - 'sublevel' or 'superlevel' (default: 'sublevel')
construction (str) - 'T' or 'V' for cubical complex construction (default: 'T')
max_dimension (int) - Maximum homology dimension to compute, -1 for auto (default: -1)

Returns:

barcodes (dict) - Dictionary mapping dimension to barcode array
- Keys: homology dimensions (0, 1, 2, …)
- Values: NumPy arrays of shape (n_features, 2) with [birth, death] columns

Example:

import numpy as np
from medtda.ph_computer import compute_ph

# 3D medical image
image = np.random.rand(64, 64, 64)

# Compute H0, H1, H2
barcodes = compute_ph(
    image,
    filtration_type='sublevel',
    max_dimension=2
)

# Access barcodes by dimension
h0_barcode = barcodes['H0']  # Connected components
h1_barcode = barcodes['H1']  # Loops/tunnels
h2_barcode = barcodes['H2']  # Voids/cavities

# Each barcode is an array of [birth, death] pairs
for birth, death in h1_barcode:
    persistence = death - birth
    print(f"H1 feature: birth={birth:.3f}, death={death:.3f}, persistence={persistence:.3f}")

clean_and_cap_ph 

medtda.ph_computer.clean_and_cap_ph(ph_array, cap_value)[source]

Clean and cap a persistence homology array.

Removes zero-persistence bars and caps infinite values.

Parameters:

ph_array (array-like) – Persistence homology array with columns [dimension, birth, death].
cap_value (float) – Value to use for capping infinite death times.

Returns:

Cleaned and capped persistence array.

Return type:

numpy.ndarray

Clean and cap persistence values to handle infinite persistence.

Signature:

def clean_and_cap_ph(
    ph_array: np.ndarray,
    cap_value: float
) -> np.ndarray

Parameters:

ph_array (np.ndarray) - Raw persistence barcode array
cap_value (float) - Value to cap infinite persistence

Returns:

cleaned_barcode (np.ndarray) - Cleaned barcode with finite values

Example:

from medtda.ph_computer import compute_ph, clean_and_cap_ph

barcodes = compute_ph(image)

# H0 often has infinite persistence features
h0_raw = barcodes['H0']
h0_clean = clean_and_cap_ph(h0_raw, cap_value=1.0)

Parameters 

Filtration Type 

sublevel

Build filtration from low to high intensity values.

Features appear at low intensities (birth)
Features disappear at high intensities (death)
Use for: bright features on dark background

superlevel

Build filtration from high to low intensity values.

Features appear at high intensities
Features disappear at low intensities
Use for: dark features on bright background

Example:

# Dark lesions on bright CT scan
barcodes_sublevel = compute_ph(image, filtration_type='sublevel')

# Bright lesions on dark background
barcodes_superlevel = compute_ph(image, filtration_type='superlevel')

Construction Method 

T-construction (default)

Each voxel becomes a top-dimensional cell
More features, finer detail
Slightly slower computation

V-construction

Each voxel becomes a vertex
Fewer features, coarser detail
Slightly faster computation

Example:

# T-construction: finer detail
barcodes_t = compute_ph(image, construction='T')

# V-construction: faster computation
barcodes_v = compute_ph(image, construction='V')

Max Dimension 

-1 (auto)

Automatically determines max dimension from image:

2D images → max_dimension = 1 (H0, H1)
3D images → max_dimension = 2 (H0, H1, H2)

0

Compute only H0 (connected components):

barcodes = compute_ph(image, max_dimension=0)
# Returns {0: h0_barcode}

1

Compute H0 and H1:

barcodes = compute_ph(image, max_dimension=1)
# Returns {0: h0_barcode, 1: h1_barcode}

2

Compute H0, H1, and H2:

barcodes = compute_ph(image, max_dimension=2)
# Returns {0: h0_barcode, 1: h1_barcode, 2: h2_barcode}

Homology Dimensions 

H0: Connected Components 

Represents connected regions in the image:

h0_barcode = barcodes['H0']

# Number of significant connected components
# (persistence > threshold)
threshold = 0.1
significant_components = np.sum(
    (h0_barcode[:, 1] - h0_barcode[:, 0]) > threshold
)

H1: Loops and Tunnels 

Represents 1-dimensional holes:

h1_barcode = barcodes['H1']

# Features with high persistence indicate prominent loops
persistent_loops = h1_barcode[
    (h1_barcode[:, 1] - h1_barcode[:, 0]) > 0.2
]

H2: Voids and Cavities 

Represents 2-dimensional voids (3D images only):

h2_barcode = barcodes['H2']

# Count significant cavities
significant_voids = np.sum(
    (h2_barcode[:, 1] - h2_barcode[:, 0]) > 0.15
)

Complete Examples 

Example 1: 2D Image Analysis 

import numpy as np
from medtda.ph_computer import compute_ph
from medtda.loaders import load_image

# Load 2D image
image, metadata = load_image('xray.png')

# Compute persistence
barcodes = compute_ph(
    image,
    filtration_type='sublevel',
    construction='T',
    max_dimension=1  # H0 and H1 only for 2D
)

# Analyze connected components
h0 = barcodes['H0']
persistence_h0 = h0[:, 1] - h0[:, 0]

print(f"Connected components: {len(h0)}")
print(f"Mean persistence: {persistence_h0.mean():.3f}")
print(f"Max persistence: {persistence_h0.max():.3f}")

# Analyze loops
h1 = barcodes['H1']
persistence_h1 = h1[:, 1] - h1[:, 0]

print(f"Loops detected: {len(h1)}")
print(f"Significant loops (p>0.1): {np.sum(persistence_h1 > 0.1)}")

Example 2: 3D Medical Image 

import numpy as np
from medtda.ph_computer import compute_ph
from medtda import Preprocessor
from medtda.loaders import load_image, load_mask

# Load and preprocess
image, _ = load_image('ct_scan.nii.gz')
mask, _ = load_mask('organ_mask.nii.gz')

preprocessor = Preprocessor(
    normalize=True,
    crop_to_roi=True
)
processed = preprocessor.preprocess(image, mask)

# Compute all homology dimensions
barcodes = compute_ph(
    processed,
    filtration_type='sublevel',
    construction='T',
    max_dimension=2
)

# Analyze 3D topology
for dim, barcode in barcodes.items():
    persistence = barcode[:, 1] - barcode[:, 0]

    dim_names = {0: 'Connected Components', 1: 'Tunnels', 2: 'Voids'}
    print(f"\n{dim_names[dim]} (H{dim}):")
    print(f"  Total features: {len(barcode)}")
    print(f"  Mean persistence: {persistence.mean():.4f}")
    print(f"  Std persistence: {persistence.std():.4f}")
    print(f"  Max persistence: {persistence.max():.4f}")

Example 3: Comparing Filtrations 

from medtda.ph_computer import compute_ph

# Dark features (e.g., vessels on bright background)
barcodes_sub = compute_ph(ct_image, filtration_type='sublevel')

# Bright features (e.g., contrast-enhanced regions)
barcodes_sup = compute_ph(ct_image, filtration_type='superlevel')

# Compare
print("Sublevel filtration:")
print(f"  H1 features: {len(barcodes_sub[1])}")

print("Superlevel filtration:")
print(f"  H1 features: {len(barcodes_sup[1])}")

Performance Tips 

Choose appropriate max_dimension:
- For 2D images, max_dimension=1 is usually sufficient
- Computing H2 on 3D images can be slow for large volumes
Construction method:
- Use T-construction for detailed analysis
- Use V-construction for faster computation on large images
Image size:
- Downsampling or cropping to ROI can significantly speed up computation
- Consider preprocessing with crop_to_roi=True
Infinite persistence:
- H0 always has one infinite persistence feature
- Use clean_and_cap_ph() if needed for downstream analysis