Basic Usage

This example demonstrates the fundamental usage patterns of MedTDA for feature extraction.

Simple Feature Extraction 

Minimal Example 

Extract features from a single image with default settings:

from medtda import FeatureExtractor

# Create extractor with defaults
extractor = FeatureExtractor()

# Extract features
features = extractor.execute('path/to/image.nii.gz')

print(f"Number of features extracted: {len(features)}")

Output:

Number of features extracted: 114

With Custom Parameters 

Configure preprocessing and vectorization:

from medtda import FeatureExtractor

extractor = FeatureExtractor(
    # Preprocessing
    normalize=True,
    normalize_method='robust',

    # Persistent homology
    filtration_type='sublevel',
    max_dimension=2,

    # Vectorization
    vectorization_method='PersImage'
)

features = extractor.execute('ct_scan.nii.gz')

# features is a flat dict: {'PersImage_H0_f0': ..., 'PersImage_H0_f1': ..., ...}
print(f"Number of persistence image features: {len(features)}")

Output:

Number of persistence image features: 1200

Using Masks 

Binary Mask 

Process only the region of interest:

from medtda import FeatureExtractor

extractor = FeatureExtractor(
    normalize=True,
    crop_to_roi=True,
    roi_padding=5
)

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

What this does:

Loads image and mask
Crops to liver bounding box (with 5-pixel padding)
Applies mask (sets background to 0)
Normalizes intensities
Computes persistent homology
Returns feature vector

Multi-Label Mask 

Extract specific label from multi-label segmentation:

from medtda import FeatureExtractor

# Multi-label mask: 0=background, 1=liver, 2=tumor
extractor = FeatureExtractor(
    label=2,  # Extract tumor only
    crop_to_roi=True,
    normalize=True
)

features = extractor.execute(
    image='ct_abdomen.nii.gz',
    mask='organs_mask.nii.gz'  # Multi-label
)

print("Extracted features from tumor region only")

Different Image Formats 

3D NIFTI Images 

Most common medical imaging format:

from medtda import FeatureExtractor
from medtda.loaders import load_image

# Load and inspect
image, metadata = load_image('brain_mri.nii.gz')
print(f"Image shape: {image.shape}")
print(f"Spacing: {metadata.get('spacing')}")

# Extract features
extractor = FeatureExtractor()
features = extractor.execute(image)  # Can pass array directly

Output:

Image shape: (256, 256, 180)
Spacing: (1.0, 1.0, 1.0)

2D PNG/TIFF Images 

For microscopy or 2D radiographs:

from medtda import FeatureExtractor

# 2D images → only H0 and H1 dimensions available
extractor = FeatureExtractor(
    max_dimension=1,  # H0 and H1 only
    normalize=True
)

features = extractor.execute('microscopy_image.png')

print("2D image processed - H0 and H1 features extracted")

NumPy Arrays 

Work directly with arrays:

import numpy as np
from medtda import FeatureExtractor

# Create or load image as numpy array
image_array = np.random.rand(100, 100, 100)

# Process array directly
extractor = FeatureExtractor(normalize=False)
features = extractor.execute(image_array)

Accessing Barcodes 

Get Raw Persistence Barcodes 

from medtda import FeatureExtractor

extractor = FeatureExtractor(
    vectorization_method='PersStats',
    return_barcodes=True  # Return barcodes too
)

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

# Access barcodes by dimension
h0_barcode = barcodes['H0']  # Connected components
h1_barcode = barcodes['H1']  # Loops

print(f"H0 features: {len(h0_barcode)}")
print(f"H1 features: {len(h1_barcode)}")

# Analyze persistence
h1_persistence = h1_barcode[:, 1] - h1_barcode[:, 0]
print(f"Mean H1 persistence: {h1_persistence.mean():.3f}")

Output:

H0 features: 47
H1 features: 12
Mean H1 persistence: 0.234

Filter Barcodes by Persistence 

import numpy as np

# Get barcodes
_, barcodes = extractor.execute('image.nii.gz')

# Filter H1 features
h1 = barcodes['H1']
persistence = h1[:, 1] - h1[:, 0]

# Keep only significant features
threshold = 0.1
significant_features = h1[persistence > threshold]

print(f"Total H1 features: {len(h1)}")
print(f"Significant (p > {threshold}): {len(significant_features)}")

Multiple Vectorization Methods 

Single Method 

from medtda import FeatureExtractor

extractor = FeatureExtractor(
    vectorization_method='PersImage'
)

features = extractor.execute('image.nii.gz')

# features is a dict with keys like 'PersImage_H0_f0', 'PersImage_H0_f1', etc.
print(f\"Extracted {len(features)} persistence image features\")

Multiple Methods 

from medtda import FeatureExtractor

extractor = FeatureExtractor()

# Enable multiple methods
extractor.enable_vectorization_methods([
    'PersStats',
    'BettiCurve',
    'PersImage'
])

features = extractor.execute('image.nii.gz')

# features is a flat dict with all method features combined
# Keys: 'PersStats_H0_mean', 'BettiCurve_H0_f0', 'PersImage_H0_f0', ...
stats_features = {k: v for k, v in features.items() if k.startswith('PersStats')}
betti_features = {k: v for k, v in features.items() if k.startswith('BettiCurve')}
pi_features = {k: v for k, v in features.items() if k.startswith('PersImage')}
print(f"PersStats features: {len(stats_features)}")
print(f"BettiCurve features: {len(betti_features)}")
print(f"PersImage features: {len(pi_features)}")

Output:

PersStats features: 26
BettiCurve features: 200
PersImage features: 800

Configure Method Parameters 

from medtda import FeatureExtractor

extractor = FeatureExtractor()
extractor.enable_vectorization_methods([
    'BettiCurve',
    'PersImage'
])

# Configure individual methods
extractor.set_vectorization_params(
    'BettiCurve',
    resolution=200  # Higher resolution
)

extractor.set_vectorization_params(
    'PersImage',
    resolution=30,    # Finer grid
    bandwidth=0.1     # Smaller bandwidth
)

features = extractor.execute('image.nii.gz')

Complete Workflow Examples 

Example 1: CT Liver Analysis 

from medtda import FeatureExtractor
from medtda.loaders import load_image, load_mask
import numpy as np

# Load data
image, img_meta = load_image('liver_ct.nii.gz')
mask, mask_meta = load_mask('liver_segmentation.nii.gz')

print(f"Image shape: {image.shape}")
print(f"Image spacing: {img_meta.get('spacing')}")

# Configure for CT
extractor = FeatureExtractor(
    # CT-specific preprocessing
    spacing=(1.0, 1.0, 1.0),      # Resample to isotropic
    window=(40, 400),              # Soft tissue window
    normalize=True,
    normalize_method='robust',

    # ROI processing
    crop_to_roi=True,
    roi_padding=5,

    # PH parameters
    filtration_type='sublevel',
    max_dimension=2,

    # Vectorization
    vectorization_method='PersImage'
)

# Extract features
features = extractor.execute(image, mask)

# features is a flat dict: {'PersImage_H0_f0': ..., 'PersImage_H1_f0': ..., ...}
feature_vector = np.array(list(features.values()))
print(f"Feature vector size: {len(feature_vector)}")

# Could now use for ML:
# model.predict(feature_vector.reshape(1, -1))

Example 2: Brain MRI Tumor Classification 

from medtda import FeatureExtractor
import numpy as np

# Configure for brain MRI
extractor = FeatureExtractor(
    # Preprocessing
    normalize=True,
    normalize_method='zscore',  # Good for MRI
    crop_to_roi=True,

    # Use multiple methods for classification
    vectorization_method=[
        'PersStats',
        'BettiCurve',
        'PersImage'
    ]
)

# Process
features = extractor.execute(
    'brain_t1.nii.gz',
    'tumor_mask.nii.gz'
)

# features is a flat dict with all method features combined
# Keys: 'PersStats_H0_mean', 'BettiCurve_H0_f0', 'PersImage_H0_f0', ...
all_features = np.array(list(features.values()))

print(f"Combined feature vector: {len(all_features)} dimensions")

Example 3: Microscopy Image Analysis 

from medtda import FeatureExtractor
from medtda.loaders import load_image

# Load 2D microscopy image
image, metadata = load_image('histology.tiff')

# Configure for 2D
extractor = FeatureExtractor(
    normalize=True,
    max_dimension=1,  # Only H0 and H1 for 2D
    filtration_type='sublevel',
    vectorization_method='PersImage',
    return_barcodes=True
)

# Extract
features, barcodes = extractor.execute(image)

# Analyze
h0_persistence = barcodes['H0'][:, 1] - barcodes['H0'][:, 0]
h1_persistence = barcodes['H1'][:, 1] - barcodes['H1'][:, 0]

print(f"Cells/nuclei (H0): {np.sum(h0_persistence > 0.05)}")
print(f"Glands/structures (H1): {np.sum(h1_persistence > 0.1)}")

Example 4: Quick Feature Inspection 

from medtda import FeatureExtractor

# Fast inspection with minimal features
extractor = FeatureExtractor(
    normalize=True,
    max_dimension=1,  # Faster
    vectorization_method='PersStats'  # Fastest
)

features = extractor.execute('scan.nii.gz')

# Access named statistics directly
print("Quick topology summary:")
print(f"  H0 features detected: {features.get('PersStats_H0_count', 'N/A')}")
print(f"  H0 mean persistence: {features.get('PersStats_H0_mean', 'N/A')}")
print(f"  H0 max persistence: {features.get('PersStats_H0_max', 'N/A')}")

Using with Machine Learning 

Scikit-learn Integration 

from medtda import FeatureExtractor
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
import numpy as np

# Configure extractor
extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersImage'
)

# Process training data
X = []
y = []

for image_path, label in training_data:
    features = extractor.execute(image_path)
    feature_vector = np.array(list(features.values()))
    X.append(feature_vector)
    y.append(label)

X = np.array(X)
y = np.array(y)

# Train model
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

clf = RandomForestClassifier(n_estimators=100)
clf.fit(X_train, y_train)

# Evaluate
score = clf.score(X_test, y_test)
print(f"Accuracy: {score:.3f}")

Tips and Best Practices 

Start with Default Settings

extractor = FeatureExtractor()
features = extractor.execute(image)

Always Normalize

extractor = FeatureExtractor(normalize=True)

Use ROI Cropping for Speed

extractor = FeatureExtractor(crop_to_roi=True)

Check Image Properties First

from medtda.loaders import load_image

image, metadata = load_image('scan.nii.gz')
print(f"Shape: {image.shape}")
print(f"Value range: [{image.min()}, {image.max()}]")
print(f"Spacing: {metadata.get('spacing')}")

Start with Fast Methods

# For initial exploration
extractor = FeatureExtractor(
    max_dimension=1,  # Skip H2
   vectorization_method='PersStats'  # Fastest

Common Pitfalls 

Issue: Shape Mismatch

# Wrong: mask doesn't match image
features = extractor.execute('image.nii.gz', 'wrong_mask.nii.gz')
# Error: ValueError: Image and mask shapes don't match

# Fix: Verify compatibility
from medtda.loaders import load_image, load_mask, validate_compatibility

image, _ = load_image('image.nii.gz')
mask, _ = load_mask('mask.nii.gz')
validate_compatibility(image, mask)

Issue: Empty Barcodes

# If ROI is too small or homogeneous
features, barcodes = extractor.execute(image, mask)

# Check for empty barcodes
for dim, barcode in barcodes.items():
    if len(barcode) == 0:
        print(f"Warning: No H{dim} features detected")

Issue: Memory Error on Large Images

# Problem: 3D image too large

# Solution 1: Downsample
extractor = FeatureExtractor(spacing=(2.0, 2.0, 2.0))

# Solution 2: Crop to ROI
extractor = FeatureExtractor(crop_to_roi=True)

# Solution 3: Process 2D slices instead

Next Steps 

Advanced Preprocessing - More preprocessing control
Batch Workflow - Process multiple images
User Guide - Complete user guide