Batch Workflow

This example demonstrates how to process multiple medical images efficiently, including batch feature extraction, parallel processing, and handling large datasets.

Basic Batch Processing

Process Directory of Images

from medtda import FeatureExtractor
from medtda.loaders import load_image
from pathlib import Path
import numpy as np
import pandas as pd

# Setup
image_dir = Path('data/images')
output_file = 'features.csv'

# Create extractor
extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

# Process all images
results = []

for img_path in image_dir.glob('*.nii.gz'):
    print(f"Processing {img_path.name}...")

    # Extract features
    features = extractor.execute(str(img_path))

    # Store result
    results.append({
        'filename': img_path.name,
        **features  # flat dict of all feature values
    })

# Save to DataFrame
df = pd.DataFrame(results)
df.to_csv(output_file, index=False)

print(f"Processed {len(results)} images")
print(f"Saved to {output_file}")

With Masks

from medtda import FeatureExtractor
from pathlib import Path
import pandas as pd

image_dir = Path('data/images')
mask_dir = Path('data/masks')

extractor = FeatureExtractor(
    normalize=True,
    crop_to_roi=True,
    vectorization_method='PersStats'
)

results = []

for img_path in image_dir.glob('*.nii.gz'):
    # Find corresponding mask
    mask_path = mask_dir / img_path.name

    if not mask_path.exists():
        print(f"Warning: No mask for {img_path.name}, skipping")
        continue

    print(f"Processing {img_path.name}...")

    # Extract with mask
    features = extractor.execute(
        str(img_path),
        str(mask_path)
    )

    results.append({
        'filename': img_path.name,
        **features  # flat dict of all feature values
    })

df = pd.DataFrame(results)
print(f"Processed {len(results)} images with masks")

Error Handling

Robust Batch Processing

from medtda import FeatureExtractor
from pathlib import Path
import traceback

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

image_dir = Path('data/images')

successes = []
failures = []

for img_path in image_dir.glob('*.nii.gz'):
    try:
        print(f"Processing {img_path.name}...")
        features = extractor.execute(str(img_path))

        successes.append({
            'filename': img_path.name,
            **features  # flat dict of all feature values
        })

    except Exception as e:
        print(f"  ERROR: {str(e)}")
        failures.append({
            'filename': img_path.name,
            'error': str(e),
            'traceback': traceback.format_exc()
        })

print(f"\nSuccess: {len(successes)}")
print(f"Failures: {len(failures)}")

# Log failures
if failures:
    with open('failures.log', 'w') as f:
        for fail in failures:
            f.write(f"\n{fail['filename']}:\n")
            f.write(f"{fail['traceback']}\n")
            f.write("-" * 80 + "\n")

Validation and Filtering

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

def is_valid_pair(image_path, mask_path):
    """Check if image-mask pair is valid."""
    try:
        # Load
        image, img_meta = load_image(str(image_path))
        mask, mask_meta = load_mask(str(mask_path))

        # Check shapes match
        if image.shape != mask.shape:
            print(f"  Shape mismatch: {image.shape} vs {mask.shape}")
            return False

        # Check mask not empty
        if np.sum(mask > 0) == 0:
            print(f"  Empty mask")
            return False

        # Check ROI size
        roi_voxels = np.sum(mask > 0)
        if roi_voxels < 100:  # Minimum voxels
            print(f"  ROI too small: {roi_voxels} voxels")
            return False

        return True

    except Exception as e:
        print(f"  Validation error: {str(e)}")
        return False

# Process with validation
image_dir = Path('data/images')
mask_dir = Path('data/masks')

extractor = FeatureExtractor(
    normalize=True,
    crop_to_roi=True,
    vectorization_method='PersStats'
)

results = []
skipped = 0

for img_path in image_dir.glob('*.nii.gz'):
    mask_path = mask_dir / img_path.name

    if not mask_path.exists():
        skipped += 1
        continue

    print(f"Validating {img_path.name}...")

    if not is_valid_pair(img_path, mask_path):
        skipped += 1
        continue

    print(f"  Processing...")
    features = extractor.execute(str(img_path), str(mask_path))

    results.append({
        'filename': img_path.name,
        **features  # flat dict of all feature values
    })

print(f"\nProcessed: {len(results)}")
print(f"Skipped: {skipped}")

Parallel Processing

Using Multiprocessing

from medtda import FeatureExtractor
from pathlib import Path
from multiprocessing import Pool, cpu_count
import numpy as np

def process_image(img_path):
    """Process single image (for parallel execution)."""
    extractor = FeatureExtractor(
        normalize=True,
        vectorization_method='PersStats'
    )

    try:
        features = extractor.execute(str(img_path))
        return {
            'filename': img_path.name,
            **features,  # flat dict of all feature values
            'success': True
        }
    except Exception as e:
        return {
            'filename': img_path.name,
            'error': str(e),
            'success': False
        }

# Get all image paths
image_dir = Path('data/images')
image_paths = list(image_dir.glob('*.nii.gz'))

print(f"Processing {len(image_paths)} images...")
print(f"Using {cpu_count()} CPUs")

# Process in parallel
with Pool(processes=cpu_count()) as pool:
    results = pool.map(process_image, image_paths)

# Separate successes and failures
successes = [r for r in results if r['success']]
failures = [r for r in results if not r['success']]

print(f"\nSuccess: {len(successes)}")
print(f"Failures: {len(failures)}")

With Images and Masks

from medtda import FeatureExtractor
from pathlib import Path
from multiprocessing import Pool

def process_pair(paths):
    """Process image-mask pair."""
    img_path, mask_path = paths

    extractor = FeatureExtractor(
        normalize=True,
        crop_to_roi=True,
        vectorization_method='PersStats'
    )

    try:
        features = extractor.execute(str(img_path), str(mask_path))
        return {
            'filename': img_path.name,
            **features,  # flat dict of all feature values
            'success': True
        }
    except Exception as e:
        return {
            'filename': img_path.name,
            'error': str(e),
            'success': False
        }

# Create pairs
image_dir = Path('data/images')
mask_dir = Path('data/masks')

pairs = []
for img_path in image_dir.glob('*.nii.gz'):
    mask_path = mask_dir / img_path.name
    if mask_path.exists():
        pairs.append((img_path, mask_path))

print(f"Processing {len(pairs)} pairs...")

# Parallel processing
with Pool() as pool:
    results = pool.map(process_pair, pairs)

successes = [r for r in results if r['success']]
print(f"Processed {len(successes)} pairs")

Progress Tracking

Using tqdm

from medtda import FeatureExtractor
from pathlib import Path
from tqdm import tqdm

image_dir = Path('data/images')
image_paths = list(image_dir.glob('*.nii.gz'))

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

results = []

# Process with progress bar
for img_path in tqdm(image_paths, desc="Processing images"):
    features = extractor.execute(str(img_path))
    results.append({
        'filename': img_path.name,
        **features  # flat dict of all feature values
    })

print(f"Done! Processed {len(results)} images")

Custom Progress Reporting

from medtda import FeatureExtractor
from pathlib import Path
import time

image_dir = Path('data/images')
image_paths = list(image_dir.glob('*.nii.gz'))

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

results = []
total = len(image_paths)
start_time = time.time()

for i, img_path in enumerate(image_paths, 1):
    # Process
    features = extractor.execute(str(img_path))
    results.append({
        'filename': img_path.name,
        **features  # flat dict of all feature values
    })

    # Report progress
    if i % 10 == 0 or i == total:
        elapsed = time.time() - start_time
        rate = i / elapsed
        eta = (total - i) / rate if rate > 0 else 0

        print(f"[{i}/{total}] {i/total*100:5.1f}% | "
              f"Rate: {rate:.1f} img/s | "
              f"ETA: {eta/60:.1f} min")

print(f"\nComplete! Total time: {elapsed/60:.1f} min")

Different Vectorization Methods

Multiple Methods per Image

from medtda import FeatureExtractor
from pathlib import Path
import numpy as np
import pandas as pd

image_dir = Path('data/images')

# Configure extractor with multiple methods
extractor = FeatureExtractor(normalize=True)
extractor.enable_vectorization_methods([
    'PersStats',
    'BettiCurve',
    'PersImage'
])

all_results = []

for img_path in image_dir.glob('*.nii.gz'):
    print(f"Processing {img_path.name}...")

    features = extractor.execute(str(img_path))

    # features is a flat dict with all method features combined
    # Keys: 'PersStats_H0_mean', 'BettiCurve_H0_f0', 'PersImage_H0_f0', ...
    all_results.append({
        'filename': img_path.name,
        **features
    })

# Save all features to a single CSV
df = pd.DataFrame(all_results)
df.to_csv('features_all_methods.csv', index=False)
print(f"Saved {len(all_results)} results with combined features")

Compare Methods on Same Data

from medtda import FeatureExtractor
from pathlib import Path
import numpy as np

image_dir = Path('data/images')
image_paths = list(image_dir.glob('*.nii.gz'))[:10]  # First 10

methods = [
    'PersStats',
    'BettiCurve',
    'PersImage'
]

# Extract with each method
method_results = {method: [] for method in methods}

for method in methods:
    print(f"\nExtracting with {method}...")

    extractor = FeatureExtractor(
        normalize=True,
        vectorization_method=method
    )

    for img_path in image_paths:
        features = extractor.execute(str(img_path))
        # features is a flat dict; collect all values as a vector
        method_results[method].append(np.array(list(features.values())))

    print(f"  Extracted {len(method_results[method])} feature vectors")

# Convert to arrays
for method in methods:
    method_results[method] = np.array(method_results[method])
    print(f"{method}: shape {method_results[method].shape}")

Organized Output Structure

Save to Structured Directory

from medtda import FeatureExtractor
from pathlib import Path
import numpy as np
import json

image_dir = Path('data/images')
output_dir = Path('output')

# Create output structure
features_dir = output_dir / 'features'
metadata_dir = output_dir / 'metadata'
features_dir.mkdir(parents=True, exist_ok=True)
metadata_dir.mkdir(parents=True, exist_ok=True)

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats',
    return_barcodes=True
)

for img_path in image_dir.glob('*.nii.gz'):
    stem = img_path.stem.replace('.nii', '')  # Remove .nii.gz

    print(f"Processing {img_path.name}...")

    # Extract
    features, barcodes = extractor.execute(str(img_path))

    # Save features as array of all values
    feature_values = np.array(list(features.values()))
    feature_file = features_dir / f"{stem}_features.npy"
    np.save(feature_file, feature_values)

    # Save barcodes
    barcode_file = features_dir / f"{stem}_barcodes.npz"
    np.savez(barcode_file, **{f'H{i}': bc for i, bc in barcodes.items()})

    # Save metadata
    metadata = {
        'filename': img_path.name,
        'feature_dim': len(feature_values),
        'barcode_dims': {f'H{i}': len(bc) for i, bc in barcodes.items()}
    }

    metadata_file = metadata_dir / f"{stem}_metadata.json"
    with open(metadata_file, 'w') as f:
        json.dump(metadata, f, indent=2)

print(f"\nOutput saved to {output_dir}/")

HDF5 Storage

from medtda import FeatureExtractor
from pathlib import Path
import h5py
import numpy as np

image_dir = Path('data/images')
output_file = 'features.h5'

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

# Create HDF5 file
with h5py.File(output_file, 'w') as hf:
    # Create group
    features_group = hf.create_group('features')

    for i, img_path in enumerate(image_dir.glob('*.nii.gz')):
        print(f"Processing {img_path.name}...")

        features = extractor.execute(str(img_path))

        # Store all feature values as an array
        dataset_name = img_path.stem.replace('.nii', '')
        features_group.create_dataset(
            dataset_name,
            data=np.array(list(features.values()))
        )

print(f"Saved to {output_file}")

# Load later
with h5py.File(output_file, 'r') as hf:
    print(f"\nStored features: {list(hf['features'].keys())[:5]}...")

Dataset Creation

Train/Test Split

from medtda import FeatureExtractor
from pathlib import Path
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd

# Load paths and labels
image_dir = Path('data/images')

# Assuming labels in CSV file
labels_df = pd.read_csv('data/labels.csv')

# Split
train_files, test_files, train_labels, test_labels = train_test_split(
    labels_df['filename'],
    labels_df['label'],
    test_size=0.2,
    stratify=labels_df['label'],
    random_state=42
)

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

def extract_features(file_list):
    """Extract features for file list."""
    features = []
    for filename in file_list:
        img_path = image_dir / filename
        feat = extractor.execute(str(img_path))
        features.append(np.array(list(feat.values())))
    return np.array(features)

print("Extracting train features...")
X_train = extract_features(train_files)

print("Extracting test features...")
X_test = extract_features(test_files)

# Save
np.savez('dataset.npz',
         X_train=X_train, y_train=train_labels,
         X_test=X_test, y_test=test_labels)

print(f"Train: {X_train.shape}, Test: {X_test.shape}")

Cross-Validation Folds

from medtda import FeatureExtractor
from pathlib import Path
from sklearn.model_selection import KFold
import numpy as np
import pandas as pd

# Load data
labels_df = pd.read_csv('data/labels.csv')
image_dir = Path('data/images')

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

# Extract all features first
print("Extracting all features...")
X = []
y = []

for _, row in labels_df.iterrows():
    img_path = image_dir / row['filename']
    features = extractor.execute(str(img_path))
    X.append(np.array(list(features.values())))
    y.append(row['label'])

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

# Create CV folds
kfold = KFold(n_splits=5, shuffle=True, random_state=42)

for fold, (train_idx, val_idx) in enumerate(kfold.split(X)):
    X_train, X_val = X[train_idx], X[val_idx]
    y_train, y_val = y[train_idx], y[val_idx]

    # Save fold
    np.savez(f'fold_{fold}.npz',
             X_train=X_train, y_train=y_train,
             X_val=X_val, y_val=y_val)

    print(f"Fold {fold}: train={len(X_train)}, val={len(X_val)}")

Quality Control

Feature Statistics

from medtda import FeatureExtractor
from pathlib import Path
import numpy as np
import pandas as pd

image_dir = Path('data/images')

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats'
)

all_features = []
filenames = []

for img_path in image_dir.glob('*.nii.gz'):
    features = extractor.execute(str(img_path))
    all_features.append(np.array(list(features.values())))
    filenames.append(img_path.name)

all_features = np.array(all_features)

# Compute statistics
print("Feature Statistics:")
print(f"  Shape: {all_features.shape}")
print(f"  Mean: {all_features.mean(axis=0)[:5]}...")  # First 5
print(f"  Std:  {all_features.std(axis=0)[:5]}...")
print(f"  Min:  {all_features.min(axis=0)[:5]}...")
print(f"  Max:  {all_features.max(axis=0)[:5]}...")

# Check for outliers (features outside 3 sigma)
mean = all_features.mean(axis=0)
std = all_features.std(axis=0)
outliers = np.abs(all_features - mean) > 3 * std

outlier_counts = outliers.sum(axis=0)
print(f"\nOutliers per feature: {outlier_counts}")

# Identify outlier samples
outlier_samples = outliers.any(axis=1)
outlier_files = [f for f, is_outlier in zip(filenames, outlier_samples) if is_outlier]

print(f"\nOutlier samples: {len(outlier_files)}")
for f in outlier_files[:10]:
    print(f"  {f}")

Missing Data Report

from medtda import FeatureExtractor
from pathlib import Path
import numpy as np

image_dir = Path('data/images')

extractor = FeatureExtractor(
    normalize=True,
    vectorization_method='PersStats',
    return_barcodes=True
)

report = []

for img_path in image_dir.glob('*.nii.gz'):
    features, barcodes = extractor.execute(str(img_path))

    # Check for issues
    issues = []

    # Check for NaN
    feature_values = np.array(list(features.values()))
    if np.isnan(feature_values).any():
        issues.append('NaN in features')

    # Check for empty barcodes
    for dim, barcode in barcodes.items():
        if len(barcode) == 0:
            issues.append(f'Empty H{dim} barcode')

    # Check for suspiciously low feature counts
    n_features = features.get('PersStats_H0_count', 0)
    if n_features < 5:
        issues.append(f'Low feature count: {n_features}')

    if issues:
        report.append({
            'filename': img_path.name,
            'issues': ', '.join(issues)
        })

# Print report
if report:
    print(f"Found {len(report)} images with issues:\n")
    for item in report:
        print(f"{item['filename']}: {item['issues']}")
else:
    print("All images processed successfully!")

Complete Workflow Example

End-to-End Pipeline

from medtda import FeatureExtractor
from pathlib import Path
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import classification_report
import numpy as np
import pandas as pd
import json

# Configuration
config = {
    'image_dir': 'data/images',
    'mask_dir': 'data/masks',
    'labels_file': 'data/labels.csv',
    'output_dir': 'output',
    'test_size': 0.2,
    'n_jobs': 4
}

# Setup
image_dir = Path(config['image_dir'])
mask_dir = Path(config['mask_dir'])
output_dir = Path(config['output_dir'])
output_dir.mkdir(exist_ok=True)

# Load labels
labels_df = pd.read_csv(config['labels_file'])
print(f"Loaded {len(labels_df)} labels")

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

# Extract features
print("\nExtracting features...")
X = []
y = []
filenames = []

for _, row in labels_df.iterrows():
    img_path = image_dir / row['filename']
    mask_path = mask_dir / row['filename']

    if not img_path.exists() or not mask_path.exists():
        print(f"Skipping {row['filename']} (missing files)")
        continue

    try:
        features = extractor.execute(str(img_path), str(mask_path))
        X.append(np.array(list(features.values())))
        y.append(row['label'])
        filenames.append(row['filename'])
    except Exception as e:
        print(f"Error with {row['filename']}: {e}")

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

print(f"Extracted {len(X)} feature vectors")
print(f"Feature dimension: {X.shape[1]}")

# Split data
X_train, X_test, y_train, y_test, files_train, files_test = train_test_split(
    X, y, filenames,
    test_size=config['test_size'],
    stratify=y,
    random_state=42
)

print(f"\nTrain: {len(X_train)}, Test: {len(X_test)}")

# Train model
print("\nTraining model...")
clf = RandomForestClassifier(n_estimators=100, n_jobs=config['n_jobs'], random_state=42)
clf.fit(X_train, y_train)

# Evaluate
train_score = clf.score(X_train, y_train)
test_score = clf.score(X_test, y_test)

print(f"Train accuracy: {train_score:.3f}")
print(f"Test accuracy: {test_score:.3f}")

# Cross-validation
cv_scores = cross_val_score(clf, X_train, y_train, cv=5, n_jobs=config['n_jobs'])
print(f"CV accuracy: {cv_scores.mean():.3f} ± {cv_scores.std():.3f}")

# Detailed metrics
y_pred = clf.predict(X_test)
print("\nClassification Report:")
print(classification_report(y_test, y_pred))

# Save results
results = {
    'config': config,
    'train_size': len(X_train),
    'test_size': len(X_test),
    'feature_dim': X.shape[1],
    'train_accuracy': float(train_score),
    'test_accuracy': float(test_score),
    'cv_accuracy': float(cv_scores.mean()),
    'cv_std': float(cv_scores.std())
}

with open(output_dir / 'results.json', 'w') as f:
    json.dump(results, f, indent=2)

# Save features
np.savez(output_dir / 'dataset.npz',
         X_train=X_train, y_train=y_train,
         X_test=X_test, y_test=y_test)

# Save file lists
pd.DataFrame({'filename': files_train, 'label': y_train}).to_csv(
    output_dir / 'train_files.csv', index=False)
pd.DataFrame({'filename': files_test, 'label': y_test}).to_csv(
    output_dir / 'test_files.csv', index=False)

print(f"\nAll results saved to {output_dir}/")

Next Steps

• Basic Usage - Basic feature extraction

• CLI Workflows - Using the CLI for batch processing

• Batch Processing - Complete batch processing guide

• FeatureExtractor - FeatureExtractor API

See Also

• Advanced Preprocessing - Custom preprocessing