Persistent Homology

This guide explains how persistent homology is computed in MedTDA and how to configure the computation parameters.

Overview

Persistent homology captures the topological features (connected components, holes, voids) in an image across multiple scales. MedTDA computes persistent homology using cubical complexes, which are well-suited for gridded image data.

Key Concepts:

Filtration - A sequence of nested topological spaces
Homology Dimensions - H₀ (components), H₁ (loops/holes), H₂ (voids/cavities)
Persistence Barcodes - Birth and death times of topological features
Cubical Complexes - Topological representation of image pixels/voxels

For mathematical background, see Persistent Homology.

Persistent Homology Parameters

MedTDA provides three main parameters for controlling PH computation:

from medtda import FeatureExtractor

extractor = FeatureExtractor(
    filtration_type='sublevel',  # or 'superlevel'
    construction='T',             # or 'V'
    max_dimension=2              # Compute H0, H1, H2
)

Filtration Types

The filtration determines how the topological space grows.

Sublevel Filtration

Filtration parameter: intensity threshold growing from low to high

Features appear when intensity exceeds threshold
Captures dark-to-bright structures
Most common choice for medical images

extractor = FeatureExtractor(
    filtration_type='sublevel'  # Default
)

Use for:

General medical image analysis
Detecting bright structures (e.g., enhanced tumors on MRI)
Standard persistent homology applications

Example: In a brain MRI:

Low threshold → only brightest pixels (CSF, vessels)
High threshold → most of the image included
Features: bright connected regions, their boundaries

Superlevel Filtration

Filtration parameter: intensity threshold growing from high to low

Features appear when intensity falls below threshold
Captures bright-to-dark structures
Less common but useful for specific applications

extractor = FeatureExtractor(
    filtration_type='superlevel'
)

Use for:

Analyzing dark structures (e.g., air in CT, hypo-intense lesions)
Studying intensity decrease patterns
Specific research applications

Example: In CT scan:

High threshold → only darkest pixels (air, lung)
Low threshold → entire image
Features: dark connected regions

Choosing Filtration Type

# For bright structures (tumors, vessels, enhancement)
extractor = FeatureExtractor(filtration_type='sublevel')

# For dark structures (air, cysts, hypo-intense regions)
extractor = FeatureExtractor(filtration_type='superlevel')

Rule of thumb: Use sublevel unless you have a specific reason for superlevel.

Cubical Complex Construction

The construction method determines how pixels/voxels are connected.

T-Construction

Pixels/voxels represent top-cells (cubes of highest dimension):

Uses 8-neighborhood connectivity in 2D (26-neighborhood in 3D)
Captures more topological features
Standard choice for most applications

extractor = FeatureExtractor(
    construction='T'  # Default
)

Properties:

Larger homology groups
More detected features
Richer topological representation
Matches standard image analysis conventions

V-Construction

Pixels/voxels represent 0-cells (vertices):

Uses 4-neighborhood connectivity in 2D (6-neighborhood in 3D)
More conservative connectivity, fewer detected features
Alternative construction method

extractor = FeatureExtractor(
    construction='V'
)

Properties:

Smaller homology groups
Fewer detected features
More conservative topological representation
Experimental/research use

Choosing Construction Method

# Standard/recommended approach
extractor = FeatureExtractor(construction='T')

# Experimental/alternative approach
extractor = FeatureExtractor(construction='V')

Recommendation: Use 'T' unless conducting methodological research.

Homology Dimensions

Persistent homology computes features in multiple dimensions.

Dimension Meanings

H₀ (0-dimensional homology):

Connected components
Separate regions/structures in the image
Birth: region appears
Death: merges with another region

H₁ (1-dimensional homology):

Loops, holes, tunnels
Circular structures, boundaries
Birth: loop forms
Death: loop fills in

H₂ (2-dimensional homology):

Voids, cavities, enclosed volumes
3D hollow structures
Birth: void appears
Death: void fills

H₃+ (higher dimensions):

4D and higher features
Rare in medical imaging
Computational cost increases significantly

Setting Maximum Dimension

Control which dimensions to compute:

# Compute only H0 (components)
extractor = FeatureExtractor(max_dimension=0)

# Compute H0 and H1 (components and loops)
extractor = FeatureExtractor(max_dimension=1)

# Compute H0, H1, H2 (components, loops, voids)
extractor = FeatureExtractor(max_dimension=2)

# Auto-detect based on image dimensionality
extractor = FeatureExtractor(max_dimension=-1)  # Default

Auto-Detection Logic

When max_dimension=-1 (default):

2D images → compute H₀ and H₁
3D images → compute H₀, H₁, and H₂
4D images → compute H₀, H₁, H₂

This provides reasonable defaults for different image types.

Choosing Maximum Dimension

Recommendations:

# 2D images: H0 and H1 only
extractor = FeatureExtractor(max_dimension=1)

# 3D images, fast computation: H0 and H1
extractor = FeatureExtractor(max_dimension=1)

# 3D images, complete analysis: H0, H1, H2
extractor = FeatureExtractor(max_dimension=2)

# Let MedTDA decide
extractor = FeatureExtractor(max_dimension=-1)

Performance vs Information:

Higher dimensions → more features → more information
Higher dimensions → longer computation time → more memory
H₂ can be expensive for large 3D volumes

Persistence Barcodes

The output of persistent homology is a set of persistence barcodes.

Barcode Format

Each barcode is a 2D array of (birth, death) pairs:

extractor = FeatureExtractor(
    return_barcodes=True  # Return barcodes
)

features, barcodes = extractor.execute('image.nii.gz', 'mask.nii.gz')

# barcodes is a dict: {'H0': array, 'H1': array, 'H2': array}
print(barcodes['H0'].shape)  # (n_features_H0, 2)
print(barcodes['H1'].shape)  # (n_features_H1, 2)

Each row is [birth_value, death_value]:

Birth - Filtration value when feature appears
Death - Filtration value when feature disappears
Persistence - Death - Birth (lifespan of feature)

Interpreting Barcodes

import numpy as np

# Get H1 barcode (loops/holes)
barcode_h1 = barcodes['H1']

# Compute persistence (lifespan)
persistence = barcode_h1[:, 1] - barcode_h1[:, 0]

# Find most persistent features
top_indices = np.argsort(persistence)[-5:]  # Top 5
top_features = barcode_h1[top_indices]

print("Most persistent H1 features:")
for birth, death in top_features:
    print(f"  Birth: {birth:.3f}, Death: {death:.3f}, "
          f"Persistence: {death-birth:.3f}")

Long-persistence features → Important topological structures Short-persistence features → Noise, small-scale variations

Complete PH Configuration Examples

Basic Configuration

Default settings for general use:

extractor = FeatureExtractor(
    filtration_type='sublevel',
    construction='T',
    max_dimension=-1  # Auto-detect
)

Fast Configuration

Optimized for speed:

extractor = FeatureExtractor(
    filtration_type='sublevel',
    construction='T',
    max_dimension=1,  # Skip H2 computation
    crop_to_roi=True,  # Reduce image size
    spacing=(2.0, 2.0, 2.0)  # Downsample
)

Comprehensive Configuration

Complete analysis with all features:

extractor = FeatureExtractor(
    # Preprocessing
    normalize=True,
    spacing=(1.0, 1.0, 1.0),
    crop_to_roi=True,

    # Persistent Homology
    filtration_type='sublevel',
    construction='T',
    max_dimension=2,  # Compute H0, H1, H2

    # Vectorization
    vectorization_method='PersStats',
    return_barcodes=True  # Also return barcodes
)

features, barcodes = extractor.execute('image.nii.gz', 'mask.nii.gz')

Dark Structure Analysis

Analyzing dark structures (superlevel):

extractor = FeatureExtractor(
    normalize=True,
    filtration_type='superlevel',  # Dark structures
    construction='T',
    max_dimension=1,
    vectorization_method='PersStats'
)

Multiple Filtrations

Compare sublevel and superlevel:

# Bright structures
extractor_sub = FeatureExtractor(
    filtration_type='sublevel',
    vectorization_method='PersStats'
)
features_bright = extractor_sub.execute(image, mask)

# Dark structures
extractor_super = FeatureExtractor(
    filtration_type='superlevel',
    vectorization_method='PersStats'
)
features_dark = extractor_super.execute(image, mask)

# Combine both
combined_features = {**features_bright, **features_dark}

Performance Considerations

Computation Time

Factors affecting computation time:

Image size - Linear to cubic scaling
Max dimension - H₂ much slower than H₁
ROI size - Cropping significantly speeds up computation
Image complexity - More features → longer computation

Optimization Strategies

# Strategy 1: Limit dimension
extractor = FeatureExtractor(max_dimension=1)  # Skip H2

# Strategy 2: Downsample
extractor = FeatureExtractor(spacing=(2.0, 2.0, 2.0))

# Strategy 3: Crop to ROI
extractor = FeatureExtractor(crop_to_roi=True, roi_padding=1)

# Strategy 4: All optimizations
extractor = FeatureExtractor(
    max_dimension=1,
    spacing=(2.0, 2.0, 2.0),
    crop_to_roi=True
)

Typical Computation Times

Approximate times on a standard workstation:

2D image (512×512): < 1 second
3D image (128×128×64), H₀+H₁: 1-5 seconds
3D image (128×128×64), H₀+H₁+H₂: 5-30 seconds
3D image (256×256×128), H₀+H₁: 10-60 seconds
3D image (256×256×128), H₀+H₁+H₂: 1-10 minutes

Times vary significantly based on image content and complexity.

Memory Usage

Memory scales with:

Image size (storing cubical complex)
Number of features detected
Maximum dimension

Typical usage:

2D images: < 1 GB
Small 3D (128³): 1-2 GB
Large 3D (256³): 4-16 GB

Use ROI cropping to reduce memory requirements.

Troubleshooting

Computation is very slow

Reduce max_dimension to 1
Enable crop_to_roi
Increase spacing to downsample
Check image size

Out of memory error

Enable crop_to_roi with small roi_padding
Increase spacing to reduce image size
Set max_dimension=1 instead of 2
Process on a machine with more RAM

Too few features detected

Check normalization is enabled
Try different filtration_type
Verify mask is correct
Check image preprocessing

Too many features detected

Increase noise filtering in preprocessing
Use higher spacing (smoother image)
Features with low persistence are noise - vectorization handles this

Features seem incorrect

Verify filtration_type matches your structures (bright vs dark)
Check preprocessing (normalization, windowing)
Visualize barcodes to inspect raw PH output
Ensure mask is properly aligned with image

Best Practices

Use sublevel filtration for most applications
Use T-construction as standard choice
Set max_dimension explicitly instead of auto-detect for reproducibility
Enable return_barcodes during development to inspect raw output
Document PH parameters for reproducibility
Start with max_dimension=1 for exploration, add H₂ if needed
Optimize with ROI cropping for large images

Next Steps

Vectorization - Converting barcodes to feature vectors
Persistent Homology - Mathematical background
Barcodes and Diagrams - Understanding barcodes
PH Computer - PHComputer API reference
Basic Usage - Practical examples