{
  "skill": "alterlab-pathml",
  "evals": [
    {
      "id": "he-preprocessing-pipeline",
      "prompt": "I have an Aperio SVS whole-slide image of H&E tissue. Build me a PathML preprocessing pipeline that does tissue detection and Macenko stain normalization, then iterate the processed tiles.",
      "expected_output": "Invokes alterlab-pathml. Loads the slide with HESlide(...) (or SlideData(filepath=..., slide_type=types.HE); NOT a from_slide factory), builds a Pipeline([TissueDetectionHE(), StainNormalizationHE(stain_estimation_method='macenko')]), runs it via wsi.run(pipeline), and iterates wsi.tiles accessing tile.image and tile.masks. Composable WSI preprocessing with stain normalization is a core PathML capability.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "StainNormalizationHE" }
      ]
    },
    {
      "id": "nucleus-segmentation-hovernet",
      "prompt": "I want to run simultaneous nucleus segmentation and classification on my pathology tiles using a pretrained model. Which PathML model do I use and how do I set up the workflow?",
      "expected_output": "Invokes alterlab-pathml. Recommends HoVer-Net for simultaneous nucleus segmentation and classification, using PathML's PyTorch integration and DataLoaders, with optional ONNX export for inference. Deep-learning nucleus detection/segmentation is an in-scope PathML ML workflow.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "HoVer-Net" }
      ]
    },
    {
      "id": "codex-multiplex-quantification",
      "prompt": "I have a multiplexed CODEX run. I need to collapse the channel data, segment cells, quantify marker expression, and export everything to an AnnData object for single-cell analysis.",
      "expected_output": "Invokes alterlab-pathml. Loads the CODEXSlide, collapses multi-run channel data, segments cells with the Mesmer model, quantifies per-cell marker expression, and exports to AnnData. Multiparametric/spatial-proteomics imaging (CODEX, Vectra, MERFISH) is a documented PathML capability.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "CODEX" }
      ]
    },
    {
      "id": "tissue-graph-construction",
      "prompt": "After segmenting cells on my H&E slides I want to build a spatial cell graph with extracted node features so I can train a graph neural network like HACTNet.",
      "expected_output": "Invokes alterlab-pathml. Constructs cell/tissue graphs from segmented objects, extracts node features, and prepares graph representations suitable for graph neural networks (e.g. HACTNet). Spatial graph construction for GNNs is an in-scope PathML capability.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "graph" }
      ]
    },
    {
      "id": "near-miss-histolab",
      "prompt": "I just need to extract a grid of H&E tiles from a single SVS slide at 20x and save them as PNGs — no segmentation, no normalization, no ML.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-histolab. The ask is simple H&E tile extraction with no deep-learning pipeline, multiplexed imaging, or graph construction. PathML is the heavier framework; for plain tiling histolab is the lighter tool, as the description states.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-histolab" }
      ]
    },
    {
      "id": "near-miss-scanpy",
      "prompt": "I already have a cells-by-genes AnnData matrix from my CODEX quantification. Cluster the cells with Leiden, run UMAP, and find marker genes per cluster.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-scanpy. The ask is downstream single-cell clustering/embedding/marker analysis on an existing AnnData matrix, not slide loading, segmentation, or image preprocessing. PathML produces the AnnData; scanpy analyzes it.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-scanpy" }
      ]
    }
  ]
}

# Data Management & Storage

## Overview

PathML provides efficient data management solutions for handling large-scale pathology datasets through HDF5 storage, tile management strategies, and optimized batch processing workflows. The framework enables seamless storage and retrieval of images, masks, features, and metadata in formats optimized for machine learning pipelines and downstream analysis.

## h5path Storage

PathML persists processed slides in its own **h5path** format (an HDF5 layout). The relevant calls are `SlideData.write(path)` and `SlideDataset.write(dir)` — there is no `to_hdf5`/`from_hdf5` method. h5path gives you compression, chunked storage, fast random access, and a hierarchical layout for images/masks/tiles/counts.

### Writing h5path

**Single slide:**
```python
from pathml.core import HESlide

# Load and process slide (run() handles tiling)
wsi = HESlide("slide.svs")
wsi.run(pipeline, tile_size=256, level=1)

# Write to disk (h5path); use the .h5path extension by convention
wsi.write("processed_slide.h5path")
```

**Multiple slides (SlideDataset):**
```python
from pathml.core import HESlide, SlideDataset
from dask.distributed import Client
import glob

slide_paths = glob.glob("data/*.svs")
dataset = SlideDataset([HESlide(p) for p in slide_paths])

client = Client(n_workers=8, threads_per_worker=2)
dataset.run(pipeline, client=client, tile_size=256, level=1)

# Writes one h5path file per slide into the directory
dataset.write("processed/")
```

### h5path File Structure

h5path files are organized hierarchically. The layout below is illustrative — the exact group names are internal and version-dependent, so read/write through the PathML API rather than hard-coding paths:

```
processed_dataset.h5path
├── slide_0/
│   ├── metadata/
│   │   ├── name
│   │   ├── level
│   │   ├── dimensions
│   │   └── ...
│   ├── tiles/
│   │   ├── tile_0/
│   │   │   ├── image  (H, W, C) array
│   │   │   ├── coords  (x, y)
│   │   │   └── masks/
│   │   │       ├── tissue
│   │   │       ├── nucleus
│   │   │       └── ...
│   │   ├── tile_1/
│   │   └── ...
│   └── features/
│       ├── tile_features  (n_tiles, n_features)
│       └── feature_names
├── slide_1/
└── ...
```

### Reading h5path back

**Reload a written slide** by constructing `SlideData` from the h5path file. To serve its tiles to a model, wrap it in `pathml.datasets.TileDataset`:

```python
from pathml.core import SlideData
from pathml.datasets import TileDataset

# Reload the processed slide
wsi = SlideData("processed_slide.h5path")

# Access tiles
for tile in wsi.tiles:
    image = tile.image
    masks = tile.masks

# Or stream tiles to a PyTorch DataLoader
tile_dataset = TileDataset("processed_slide.h5path")
```

(If your PathML version exposes a different reload entry point, confirm it in the core API docs.)

**Low-level access with h5py** is possible since h5path is HDF5 underneath, but the internal group layout is an implementation detail and can change between versions — prefer the PathML API for forward compatibility:

```python
import h5py

with h5py.File("processed_slide.h5path", "r") as f:
    print(list(f.keys()))  # inspect the layout for your version
```

## Tile Management

### Tile Generation Strategies

Tiling is driven by `generate_tiles` (a generator) or by the tiling arguments passed to `SlideData.run`/`SlideDataset.run`. Confirm the tile-size keyword (`shape` vs `tile_size`) for your installed version.

**Fixed-size tiles with no overlap** (`stride == shape`):
```python
for tile in wsi.generate_tiles(level=1, shape=256, stride=256, pad=False):
    ...
```
- Standard tile-based processing/classification; simple, no redundancy; edge effects at boundaries.

**Overlapping tiles** (`stride < shape`):
```python
for tile in wsi.generate_tiles(level=1, shape=256, stride=128):  # 50% overlap
    ...
```
- Better for segmentation/detection (reduces boundary artifacts); more tiles, redundant compute.

**Tissue-only tiling:** run a pipeline that includes `TissueDetectionHE`, then filter the resulting `wsi.tiles` by tissue coverage:
```python
wsi.run(pipeline, tile_size=256, level=1)  # pipeline includes TissueDetectionHE()

tile_area = 256 * 256
tissue_tiles = [
    tile for tile in wsi.tiles
    if tile.masks.get("tissue") is not None
    and tile.masks["tissue"].sum() / tile_area > 0.5
]
```

### Tile Stitching

PathML does not ship a general `stitch_tiles` utility; reassemble predictions yourself from each tile's `coords`. A simple non-overlapping reassembly:

```python
import numpy as np

# Allocate at the working level's dimensions (W, H from the backend)
W, H = wsi.shape
canvas = np.zeros((H, W), dtype=np.float32)

for tile in wsi.tiles:
    pred = model.predict(tile.image)  # (h, w)
    j, i = tile.coords                # top-left (x, y)
    h, w = pred.shape
    canvas[i:i + h, j:j + w] = pred   # average/blend instead for overlaps
```

For overlapping tiles, accumulate predictions and a per-pixel count, then divide to average.

## Dataset Organization

### Directory Structure for Large Projects

Organize pathology projects with consistent structure:

```
project/
├── raw_slides/
│   ├── cohort1/
│   │   ├── slide001.svs
│   │   ├── slide002.svs
│   │   └── ...
│   └── cohort2/
│       └── ...
├── processed/
│   ├── cohort1/
│   │   ├── slide001.h5
│   │   ├── slide002.h5
│   │   └── ...
│   └── cohort2/
│       └── ...
├── features/
│   ├── cohort1_features.h5
│   └── cohort2_features.h5
├── models/
│   ├── hovernet_checkpoint.pth
│   └── classifier.onnx
├── results/
│   ├── predictions/
│   ├── visualizations/
│   └── metrics.csv
└── metadata/
    ├── clinical_data.csv
    └── slide_manifest.csv
```

### Metadata Management

Store slide-level and cohort-level metadata:

```python
import pandas as pd

# Slide manifest
manifest = pd.DataFrame({
    'slide_id': ['slide001', 'slide002', 'slide003'],
    'path': ['raw_slides/cohort1/slide001.svs', ...],
    'cohort': ['cohort1', 'cohort1', 'cohort2'],
    'tissue_type': ['breast', 'breast', 'lung'],
    'scanner': ['Aperio', 'Hamamatsu', 'Aperio'],
    'magnification': [40, 40, 20],
    'staining': ['H&E', 'H&E', 'H&E']
})

manifest.to_csv('metadata/slide_manifest.csv', index=False)

# Clinical data
clinical = pd.DataFrame({
    'slide_id': ['slide001', 'slide002', 'slide003'],
    'patient_id': ['P001', 'P002', 'P003'],
    'age': [55, 62, 48],
    'diagnosis': ['invasive', 'in_situ', 'invasive'],
    'stage': ['II', 'I', 'III'],
    'outcome': ['favorable', 'favorable', 'poor']
})

clinical.to_csv('metadata/clinical_data.csv', index=False)

# Load and merge
manifest = pd.read_csv('metadata/slide_manifest.csv')
clinical = pd.read_csv('metadata/clinical_data.csv')
data = manifest.merge(clinical, on='slide_id')
```

## Batch Processing Strategies

### Sequential Processing

Process slides one at a time (memory-efficient):

```python
import glob
from pathml.core import HESlide
from pathml.preprocessing import Pipeline

slide_paths = glob.glob("raw_slides/**/*.svs", recursive=True)

for slide_path in slide_paths:
    wsi = HESlide(slide_path)
    wsi.run(pipeline, tile_size=256, level=1)

    output_path = slide_path.replace("raw_slides", "processed").replace(".svs", ".h5path")
    wsi.write(output_path)
    print(f"Processed: {slide_path}")
```

### Parallel Processing with Dask

Process multiple slides in parallel:

```python
from pathml.core import HESlide, SlideDataset
from dask.distributed import Client, LocalCluster
from pathml.preprocessing import Pipeline

# Start Dask cluster
cluster = LocalCluster(
    n_workers=8,
    threads_per_worker=2,
    memory_limit="8GB",
    dashboard_address=":8787",  # progress at localhost:8787
)
client = Client(cluster)

# Create dataset from slide objects
slide_paths = glob.glob("raw_slides/**/*.svs", recursive=True)
dataset = SlideDataset([HESlide(p) for p in slide_paths])

# Distribute processing (pass the client; tiling options go to run())
dataset.run(pipeline, client=client, tile_size=256, level=1)

# Write one h5path per slide into a directory
dataset.write("processed/")

client.close()
cluster.close()
```

### Batch Processing with Job Arrays

For HPC clusters (SLURM, PBS):

```python
# submit_jobs.py
import os
import glob

slide_paths = glob.glob('raw_slides/**/*.svs', recursive=True)

# Write slide list
with open('slide_list.txt', 'w') as f:
    for path in slide_paths:
        f.write(path + '\n')

# Create SLURM job script
slurm_script = """#!/bin/bash
#SBATCH --array=1-{n_slides}
#SBATCH --cpus-per-task=4
#SBATCH --mem=16G
#SBATCH --time=4:00:00
#SBATCH --output=logs/slide_%A_%a.out

# Get slide path for this array task
SLIDE_PATH=$(sed -n "${{SLURM_ARRAY_TASK_ID}}p" slide_list.txt)

# Run processing
python process_slide.py --slide_path $SLIDE_PATH
""".format(n_slides=len(slide_paths))

with open('submit_jobs.sh', 'w') as f:
    f.write(slurm_script)

# Submit: sbatch submit_jobs.sh
```

```python
# process_slide.py
import argparse
from pathml.core import HESlide
from pathml.preprocessing import Pipeline

parser = argparse.ArgumentParser()
parser.add_argument("--slide_path", type=str, required=True)
args = parser.parse_args()

# Load and process
wsi = HESlide(args.slide_path)

pipeline = Pipeline([...])
wsi.run(pipeline, tile_size=256, level=1)

# Save
output_path = args.slide_path.replace("raw_slides", "processed").replace(".svs", ".h5path")
wsi.write(output_path)

print(f"Processed: {args.slide_path}")
```

## Feature Extraction and Storage

### Extracting Features

```python
from pathml.core import SlideData
import torch
import numpy as np

# Load pre-trained model for feature extraction
model = torch.load('models/feature_extractor.pth')
model.eval()
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
model = model.to(device)

# Load processed slide from h5path
wsi = SlideData("processed/slide001.h5path")

# Extract features for each tile
features = []
coords = []

for tile in wsi.tiles:
    # Preprocess tile
    tile_tensor = torch.from_numpy(tile.image).permute(2, 0, 1).unsqueeze(0).float()
    tile_tensor = tile_tensor.to(device)

    # Extract features
    with torch.no_grad():
        feature_vec = model(tile_tensor).cpu().numpy().flatten()

    features.append(feature_vec)
    coords.append(tile.coords)

features = np.array(features)  # Shape: (n_tiles, feature_dim)
coords = np.array(coords)  # Shape: (n_tiles, 2)
```

### Storing Features in HDF5

```python
import h5py

# Save features
with h5py.File('features/slide001_features.h5', 'w') as f:
    f.create_dataset('features', data=features, compression='gzip')
    f.create_dataset('coords', data=coords)
    f.attrs['feature_dim'] = features.shape[1]
    f.attrs['num_tiles'] = features.shape[0]
    f.attrs['model'] = 'resnet50'

# Load features
with h5py.File('features/slide001_features.h5', 'r') as f:
    features = f['features'][:]
    coords = f['coords'][:]
    feature_dim = f.attrs['feature_dim']
```

### Feature Database for Multiple Slides

```python
# Create consolidated feature database
import h5py
import glob

feature_files = glob.glob('features/*_features.h5')

with h5py.File('features/all_features.h5', 'w') as out_f:
    for i, feature_file in enumerate(feature_files):
        slide_name = feature_file.split('/')[-1].replace('_features.h5', '')

        with h5py.File(feature_file, 'r') as in_f:
            features = in_f['features'][:]
            coords = in_f['coords'][:]

            # Store in consolidated file
            grp = out_f.create_group(f'slide_{i}')
            grp.create_dataset('features', data=features, compression='gzip')
            grp.create_dataset('coords', data=coords)
            grp.attrs['slide_name'] = slide_name

# Query features from all slides
with h5py.File('features/all_features.h5', 'r') as f:
    for slide_key in f.keys():
        slide_name = f[slide_key].attrs['slide_name']
        features = f[f'{slide_key}/features'][:]
        # Process...
```

## Data Versioning

### Version Control with DVC

Use Data Version Control (DVC) for large dataset management:

```bash
# Initialize DVC
dvc init

# Add data directory
dvc add raw_slides/
dvc add processed/

# Commit to git
git add raw_slides.dvc processed.dvc .gitignore
git commit -m "Add raw and processed slides"

# Push data to remote storage (S3, GCS, etc.)
dvc remote add -d storage s3://my-bucket/pathml-data
dvc push

# Pull data on another machine
git pull
dvc pull
```

### Checksums and Validation

Validate data integrity:

```python
import hashlib
import pandas as pd

def compute_checksum(file_path):
    """Compute MD5 checksum of file."""
    hash_md5 = hashlib.md5()
    with open(file_path, 'rb') as f:
        for chunk in iter(lambda: f.read(4096), b""):
            hash_md5.update(chunk)
    return hash_md5.hexdigest()

# Create checksum manifest
slide_paths = glob.glob('raw_slides/**/*.svs', recursive=True)
checksums = []

for slide_path in slide_paths:
    checksum = compute_checksum(slide_path)
    checksums.append({
        'path': slide_path,
        'checksum': checksum,
        'size_mb': os.path.getsize(slide_path) / 1e6
    })

checksum_df = pd.DataFrame(checksums)
checksum_df.to_csv('metadata/checksums.csv', index=False)

# Validate files
def validate_files(manifest_path):
    manifest = pd.read_csv(manifest_path)
    for _, row in manifest.iterrows():
        current_checksum = compute_checksum(row['path'])
        if current_checksum != row['checksum']:
            print(f"ERROR: Checksum mismatch for {row['path']}")
        else:
            print(f"OK: {row['path']}")

validate_files('metadata/checksums.csv')
```

## Performance Optimization

### Compression Settings

Optimize HDF5 compression for speed vs. size:

```python
import h5py

# Fast compression (less CPU, larger files)
with h5py.File('output.h5', 'w') as f:
    f.create_dataset(
        'images',
        data=images,
        compression='gzip',
        compression_opts=1  # Level 1-9, lower = faster
    )

# Maximum compression (more CPU, smaller files)
with h5py.File('output.h5', 'w') as f:
    f.create_dataset(
        'images',
        data=images,
        compression='gzip',
        compression_opts=9
    )

# Balanced (recommended)
with h5py.File('output.h5', 'w') as f:
    f.create_dataset(
        'images',
        data=images,
        compression='gzip',
        compression_opts=4,
        chunks=True  # Enable chunking for better I/O
    )
```

### Chunking Strategy

Optimize chunked storage for access patterns:

```python
# For tile-based access (access one tile at a time)
with h5py.File('tiles.h5', 'w') as f:
    f.create_dataset(
        'tiles',
        shape=(n_tiles, 256, 256, 3),
        dtype='uint8',
        chunks=(1, 256, 256, 3),  # One tile per chunk
        compression='gzip'
    )

# For channel-based access (access all tiles for one channel)
with h5py.File('tiles.h5', 'w') as f:
    f.create_dataset(
        'tiles',
        shape=(n_tiles, 256, 256, 3),
        dtype='uint8',
        chunks=(n_tiles, 256, 256, 1),  # All tiles for one channel
        compression='gzip'
    )
```

### Memory-Mapped Arrays

Use memory mapping for large arrays:

```python
import numpy as np

# Save as memory-mapped file
features_mmap = np.memmap(
    'features/features.mmap',
    dtype='float32',
    mode='w+',
    shape=(n_tiles, feature_dim)
)

# Populate
for i, tile in enumerate(wsi.tiles):
    features_mmap[i] = extract_features(tile)

# Flush to disk
features_mmap.flush()

# Load without reading into memory
features_mmap = np.memmap(
    'features/features.mmap',
    dtype='float32',
    mode='r',
    shape=(n_tiles, feature_dim)
)

# Access subset efficiently
subset = features_mmap[1000:2000]  # Only loads requested rows
```

## Best Practices

1. **Use HDF5 for processed data:** Save preprocessed tiles and features to HDF5 for fast access

2. **Separate raw and processed data:** Keep original slides separate from processed outputs

3. **Maintain metadata:** Track slide provenance, processing parameters, and clinical annotations

4. **Implement checksums:** Validate data integrity, especially after transfers

5. **Version datasets:** Use DVC or similar tools to version large datasets

6. **Optimize storage:** Balance compression level with I/O performance

7. **Organize by cohort:** Structure directories by study cohort for clarity

8. **Regular backups:** Back up both data and metadata to remote storage

9. **Document processing:** Keep logs of processing steps, parameters, and versions

10. **Monitor disk usage:** Track storage consumption as datasets grow

## Common Issues and Solutions

**Issue: HDF5 files very large**
- Increase compression level: `compression_opts=9`
- Store only necessary data (avoid redundant copies)
- Use appropriate data types (uint8 for images vs. float64)

**Issue: Slow HDF5 read/write**
- Optimize chunk size for access pattern
- Reduce compression level for faster I/O
- Use SSD storage instead of HDD
- Enable parallel HDF5 with MPI

**Issue: Running out of disk space**
- Delete intermediate files after processing
- Compress inactive datasets
- Move old data to archival storage
- Use cloud storage for less-accessed data

**Issue: Data corruption or loss**
- Implement regular backups
- Use RAID for redundancy
- Validate checksums after transfers
- Use version control (DVC)

## Additional Resources

- **HDF5 Documentation:** https://www.hdfgroup.org/solutions/hdf5/
- **h5py:** https://docs.h5py.org/
- **DVC (Data Version Control):** https://dvc.org/
- **Dask:** https://docs.dask.org/
- **PathML Data Management API:** https://pathml.readthedocs.io/en/latest/api_datasets_reference.html

# Graph Construction & Spatial Analysis

## Overview

PathML provides tools for constructing spatial graphs from tissue images to represent cellular and tissue-level relationships. Graph-based representations enable sophisticated spatial analysis, including neighborhood analysis, cell-cell interaction studies, and graph neural network applications. These graphs capture both morphological features and spatial topology for downstream computational analysis.

## Graph Types

PathML supports construction of multiple graph types:

### Cell Graphs
- Nodes represent individual cells
- Edges represent spatial proximity or biological interactions
- Node features include morphology, marker expression, cell type
- Suitable for single-cell spatial analysis

### Tissue Graphs
- Nodes represent tissue regions or superpixels
- Edges represent spatial adjacency
- Node features include tissue composition, texture features
- Suitable for tissue-level spatial patterns

### Spatial Transcriptomics Graphs
- Nodes represent spatial spots or cells
- Edges encode spatial relationships
- Node features include gene expression profiles
- Suitable for spatial omics analysis

## Graph Construction Workflow

PathML's graph module exposes builder classes (subclasses of `BaseGraphBuilder`) rather than a single `CellGraph` factory: `KNNGraphBuilder`, `RAGGraphBuilder` (region-adjacency, for superpixel/tissue graphs), and `MSTGraphBuilder` (minimum spanning tree). Tissue regions for tissue graphs come from the superpixel extractors (`SLICSuperpixelExtractor`, `ColorMergedSuperpixelExtractor`, etc.), and per-node features come from `GraphFeatureExtractor`. Verify class names and call signatures against the API docs for your installed version.

### From Segmentation to Graphs

Convert a cell instance-segmentation mask into a spatial cell graph with a builder. The builder takes the per-cell centroids and a feature matrix and returns a `pathml.graph.utils.Graph` (PyG-compatible) object.

```python
import numpy as np
from skimage.measure import regionprops
from pathml.graph import KNNGraphBuilder
from pathml.preprocessing import Pipeline, SegmentMIF

# 1. Segment cells
pipeline = Pipeline([
    SegmentMIF(model="mesmer", nuclear_channel=0, cytoplasm_channel=29),
])
slide.run(pipeline)

# 2. Get the instance mask and per-cell centroids + features
inst_map = slide.masks["cell_segmentation"]
props = regionprops(inst_map)
centroids = np.array([p.centroid[::-1] for p in props])  # (x, y) per cell
features = np.array([[p.area, p.eccentricity, p.solidity] for p in props])

# 3. Build a k-NN cell graph
builder = KNNGraphBuilder(k=5, thresh=50)  # k neighbors within a distance threshold
graph = builder.process(inst_map, features=features, centroids=centroids)

# 4. Access components (Graph is a torch_geometric-style object)
x = graph.node_features      # node feature matrix
edge_index = graph.edge_index  # (2, n_edges) connectivity
```

### Builder Choices

- `KNNGraphBuilder(k=..., thresh=...)` - connect each node to its k nearest neighbors within an optional distance threshold; good for cell graphs.
- `RAGGraphBuilder` - region-adjacency graph; connects touching regions, suited to superpixel/tissue graphs.
- `MSTGraphBuilder` - minimum spanning tree over node positions; a sparse, fully connected backbone.

For tissue graphs, first extract superpixels (e.g. `SLICSuperpixelExtractor` or `ColorMergedSuperpixelExtractor`), then pass the label map to `RAGGraphBuilder`.

## Node Features

PathML's `GraphFeatureExtractor` computes per-region features from an instance/label map; you can also compute features yourself with scikit-image and pass them to the builder as the `features` matrix.

### Morphological and Intensity Features (scikit-image)

`skimage.measure.regionprops` / `regionprops_table` covers the common morphology and per-channel intensity statistics:

```python
import numpy as np
from skimage.measure import regionprops_table

# Morphology
morph = regionprops_table(
    inst_map,
    properties=["area", "perimeter", "eccentricity", "solidity",
                "axis_major_length", "axis_minor_length", "orientation"],
)

# Per-channel mean intensity (image shape (H, W, C))
intensity = regionprops_table(
    inst_map,
    intensity_image=multichannel_image,
    properties=["intensity_mean"],
)

# Stack into a node feature matrix aligned with the builder's node order
node_features = np.column_stack([np.asarray(v) for v in morph.values()])
```

### Cell Type Annotations as Node Features

Append cell-type labels (e.g. from HoVer-Net) as an extra node-feature column or one-hot block before building the graph:

```python
cell_types = hovernet_type_predictions  # array of per-cell type ids, aligned to props order
onehot = np.eye(n_classes)[cell_types]
node_features = np.column_stack([node_features, onehot])
```

## Spatial Analysis

PathML builds the graph; for downstream spatial statistics (neighborhood enrichment, co-occurrence, interaction tests, spatial autocorrelation) move the per-cell coordinates + labels into AnnData and use **squidpy** — the dedicated, maintained tool for this. (`pathml.graph` does not provide `analyze_neighborhoods`, `spatial_clustering`, `cell_interaction_analysis`, or `spatial_statistics` functions; don't assume they exist.) See `multiparametric.md` for the AnnData-based squidpy workflow.

```python
import anndata as ad
import numpy as np
import squidpy as sq

# Wrap centroids + cell types in AnnData
adata = ad.AnnData(node_features)
adata.obsm["spatial"] = centroids
adata.obs["cell_type"] = cell_type_labels.astype(str)

# Spatial neighbors + neighborhood enrichment
sq.gr.spatial_neighbors(adata, coord_type="generic")
sq.gr.nhood_enrichment(adata, cluster_key="cell_type")
sq.pl.nhood_enrichment(adata, cluster_key="cell_type")
```

For plain spatial clustering of positions, scikit-learn's `DBSCAN`/`KMeans` on `centroids` works directly.

## Integration with Graph Neural Networks

### Convert to PyTorch Geometric Format

The `Graph` returned by a builder is already a `torch_geometric.data.Data`-style object, so it plugs straight into PyTorch Geometric (no conversion helper needed):

```python
import torch
from torch_geometric.nn import GCNConv
from pathml.graph import KNNGraphBuilder

graph = KNNGraphBuilder(k=5, thresh=50).process(inst_map, features=features, centroids=centroids)

x = graph.node_features        # node features
edge_index = graph.edge_index  # (2, n_edges)

class GNN(torch.nn.Module):
    def __init__(self, in_channels, hidden_channels, out_channels):
        super().__init__()
        self.conv1 = GCNConv(in_channels, hidden_channels)
        self.conv2 = GCNConv(hidden_channels, out_channels)

    def forward(self, data):
        x = self.conv1(data.node_features, data.edge_index).relu()
        return self.conv2(x, data.edge_index)

model = GNN(in_channels=x.shape[1], hidden_channels=64, out_channels=5)
output = model(graph)
```

### Dataset of Graphs for Multiple Slides

Use PathML's `EntityDataset` (from `pathml.datasets`) to serve multiple graphs, or collect them into a list and batch with PyG's `DataLoader`:

```python
from torch_geometric.loader import DataLoader
from pathml.graph import KNNGraphBuilder

builder = KNNGraphBuilder(k=5, thresh=50)
graphs = [
    builder.process(s_inst_map, features=s_features, centroids=s_centroids)
    for (s_inst_map, s_features, s_centroids) in per_slide_inputs
]

loader = DataLoader(graphs, batch_size=32, shuffle=True)
for batch in loader:
    output = model(batch)
    loss = criterion(output, batch.y)
    loss.backward()
    optimizer.step()
```

## Visualization

### Graph Visualization

Draw the graph directly from its `edge_index` and the cell `centroids` (no NetworkX conversion required):

```python
import matplotlib.pyplot as plt

edges = graph.edge_index.cpu().numpy()  # (2, n_edges)

fig, ax = plt.subplots(figsize=(12, 12))
# Edges
for s, d in edges.T:
    p1, p2 = centroids[s], centroids[d]
    ax.plot([p1[0], p2[0]], [p1[1], p2[1]], "b-", alpha=0.3, linewidth=0.5)
# Nodes colored by type
ax.scatter(centroids[:, 0], centroids[:, 1], c=cell_type_labels, cmap="tab10", s=20)
ax.set_aspect("equal")
ax.axis("off")
plt.title("Cell Graph")
plt.show()
```

### Overlay on Tissue Image

Same as above but `imshow` the tissue image first, then draw edges/nodes on the same axes using `centroids` in image (x, y) coordinates.

## Complete Workflow Example

```python
import numpy as np
from skimage.measure import regionprops, regionprops_table
from pathml.core import CODEXSlide
from pathml.preprocessing import Pipeline, CollapseRunsCODEX, SegmentMIF
from pathml.graph import KNNGraphBuilder

# 1. Load and preprocess slide
slide = CODEXSlide("path/to/codex", stain="IF")
pipeline = Pipeline([
    CollapseRunsCODEX(z=0),
    SegmentMIF(model="mesmer", nuclear_channel=0, cytoplasm_channel=29),
])
slide.run(pipeline)

# 2. Per-cell centroids + features from the instance mask
inst_map = slide.masks["cell_segmentation"]
props = regionprops(inst_map)
centroids = np.array([p.centroid[::-1] for p in props])  # (x, y)
morph = regionprops_table(
    inst_map, properties=["area", "perimeter", "eccentricity", "solidity"]
)
features = np.column_stack([np.asarray(v) for v in morph.values()])

# 3. Build the cell graph
graph = KNNGraphBuilder(k=6, thresh=50).process(
    inst_map, features=features, centroids=centroids
)

# 4. graph is PyG-ready (graph.node_features, graph.edge_index) -> feed to a GNN
#    For neighborhood/interaction statistics, hand centroids + labels to squidpy
#    (see "Spatial Analysis" above).
```

## Performance Considerations

**Large tissue sections:**
- Build graphs tile-by-tile, then merge
- Use sparse adjacency matrices
- Leverage GPU for feature extraction

**Memory efficiency:**
- Store only necessary edge features
- Use int32/float32 instead of int64/float64
- Batch process multiple slides

**Computational efficiency:**
- Parallelize feature extraction across cells
- Use KNN for faster neighbor queries
- Cache computed features

## Best Practices

1. **Choose an appropriate builder:** `KNNGraphBuilder` for cell graphs (tune `k`/`thresh`), `RAGGraphBuilder` for superpixel/tissue adjacency, `MSTGraphBuilder` for a sparse backbone

2. **Normalize features:** Scale morphological and intensity features for GNN compatibility

3. **Handle edge effects:** Exclude boundary cells or use tissue masks to define valid regions

4. **Validate graph construction:** Visualize graphs on small regions before large-scale processing

5. **Combine multiple feature types:** Morphology + intensity + texture provides rich representations

6. **Consider tissue context:** Tissue type affects appropriate graph parameters (connectivity, radius)

## Common Issues and Solutions

**Issue: Too many/few edges**
- Adjust `k` and the `thresh` distance cutoff on `KNNGraphBuilder`
- Verify pixel-to-micron conversion for biological relevance

**Issue: Memory errors with large graphs**
- Process tiles separately and merge graphs
- Use sparse representations and float32 features

**Issue: Missing cells at tissue boundaries**
- Use tissue masks to exclude invalid regions before building the graph

**Issue: Inconsistent feature scales**
- Normalize features: `(x - mean) / std`
- Use robust scaling for outliers

## Additional Resources

- **PathML Graph API:** https://pathml.readthedocs.io/en/latest/api_graph_reference.html
- **PyTorch Geometric:** https://pytorch-geometric.readthedocs.io/
- **NetworkX:** https://networkx.org/
- **Spatial Statistics:** Baddeley et al., "Spatial Point Patterns: Methodology and Applications with R"

# Image Loading & Formats

## Overview

PathML provides comprehensive support for loading whole-slide images (WSI) from 160+ proprietary medical imaging formats. The framework abstracts vendor-specific complexities through unified slide classes and interfaces, enabling seamless access to image pyramids, metadata, and regions of interest across different file formats.

## Supported Formats

PathML supports the following slide formats:

### Brightfield Microscopy Formats
- **Aperio SVS** (`.svs`) - Leica Biosystems
- **Hamamatsu NDPI** (`.ndpi`) - Hamamatsu Photonics
- **Leica SCN** (`.scn`) - Leica Biosystems
- **Zeiss ZVI** (`.zvi`) - Carl Zeiss
- **3DHISTECH** (`.mrxs`) - 3DHISTECH Ltd.
- **Ventana BIF** (`.bif`) - Roche Ventana
- **Generic tiled TIFF** (`.tif`, `.tiff`)

### Medical Imaging Standards
- **DICOM** (`.dcm`) - Digital Imaging and Communications in Medicine
- **OME-TIFF** (`.ome.tif`, `.ome.tiff`) - Open Microscopy Environment

### Multiparametric Imaging
- **CODEX** - Spatial proteomics imaging
- **Vectra** (`.qptiff`) - Multiplex immunofluorescence
- **MERFISH** - Multiplexed error-robust FISH

PathML leverages OpenSlide and other specialized libraries to handle format-specific nuances automatically.

## Core Classes for Loading Images

### SlideData

`SlideData` is the fundamental class for representing whole-slide images in PathML. You construct it directly (there is no `SlideData.from_slide` factory) or via a convenience subclass (`HESlide`, `VectraSlide`, `CODEXSlide`, `MultiparametricSlide`, `IHCSlide`).

**Loading from file:**
```python
from pathml.core import SlideData, HESlide, types

# Convenience subclass for H&E (recommended for brightfield H&E)
wsi = HESlide("path/to/slide.svs", name="example")

# Generic constructor. backend is a STRING ("openslide" | "bioformats" | "dicom");
# slide_type carries the stain/dimensionality (use a pathml.core.types instance).
wsi = SlideData("path/to/slide.svs", backend="openslide", slide_type=types.HE)

# OME-TIFF via the Bio-Formats backend
wsi = SlideData("path/to/slide.ome.tiff", backend="bioformats")
```

**Key attributes:**
- `wsi.slide` - Backend slide object (OpenSlide, BioFormats, etc.)
- `wsi.tiles` - Collection of image tiles (populated after a tiling/Pipeline run)
- `wsi.masks` - Slide-level masks
- `wsi.shape` - Image dimensions
- `wsi.name` - Slide name

**Methods:**
- `wsi.run(pipeline, ...)` - Run a preprocessing Pipeline (handles tiling + transforms)
- `wsi.generate_tiles()` - Generator over `Tile` objects
- `wsi.extract_region(location, size, ...)` - Read a specific region (delegates to the backend)
- `wsi.write(path)` - Write contents to disk in h5path format

### SlideType and types

`SlideType` describes a slide's stain/dimensionality (NOT the backend). The `pathml.core.types` module exposes pre-made instances you pass as `slide_type`:

```python
from pathml.core import types

types.HE       # H&E brightfield
types.IHC      # immunohistochemistry
types.IF       # immunofluorescence
types.CODEX    # CODEX multiplex IF
types.Vectra   # Vectra multiplex IF
```

The decoding backend is selected separately with the string `backend` argument (`"openslide"`, `"bioformats"`, or `"dicom"`).

### Specialized Slide Classes

PathML provides specialized slide classes for specific imaging modalities:

**CODEXSlide:**
```python
from pathml.core import CODEXSlide

# Load CODEX spatial proteomics data
codex_slide = CODEXSlide(
    path="path/to/codex_dir",
    stain="IF",  # Immunofluorescence
    backend="bioformats"
)
```

**VectraSlide:**
```python
from pathml.core import VectraSlide

# Load Vectra multiplex IF data (.qptiff); uses the Bio-Formats backend
vectra_slide = VectraSlide("path/to/vectra.qptiff")
```

**MultiparametricSlide:**
```python
from pathml.core import MultiparametricSlide

# Generic multiparametric imaging
mp_slide = MultiparametricSlide("path/to/multiparametric_data", backend="bioformats")
```

## Loading Strategies

### Tile-Based Loading

For large WSI files, tile-based loading enables memory-efficient processing. `generate_tiles` returns a generator over `Tile` objects; iterate it directly. To persist tiles on the slide, run a `Pipeline` (see `preprocessing.md`), after which they are available on `wsi.tiles`.

```python
from pathml.core import HESlide

# Load slide
wsi = HESlide("path/to/slide.svs")

# Iterate tiles lazily at a chosen level
for tile in wsi.generate_tiles(
    level=0,        # Pyramid level (0 = highest resolution)
    shape=256,      # Tile dimensions in pixels
    stride=256,     # Spacing between tiles (256 = no overlap)
    pad=False,      # Whether to pad edge tiles
):
    image = tile.image   # numpy array
    coords = tile.coords  # (i, j) coordinates
    # Process tile...
```

**Overlapping tiles:**
```python
# 50% overlap
for tile in wsi.generate_tiles(level=0, shape=256, stride=128):
    ...
```

Confirm the exact tiling keyword (`shape` vs `tile_size`) against your installed PathML version; it has differed across releases.

### Region-Based Loading

Extract specific regions of interest directly:

```python
# Read region at a specific location. extract_region delegates to the
# active backend (OpenSlide/Bio-Formats/DICOM); see your backend for the
# exact level/coordinate semantics.
region = wsi.extract_region(
    location=(10000, 15000),  # (x, y) coordinates
    size=(512, 512),          # width, height in pixels
)

# Returns a numpy array
```

### Pyramid Level Selection

Whole-slide images are stored in multi-resolution pyramids. Select the appropriate level based on desired magnification. Pyramid metadata lives on the backend object (`wsi.slide`); for OpenSlide it exposes `level_dimensions` and `level_downsamples`:

```python
# Inspect available levels via the OpenSlide backend
osh = wsi.slide.slide  # underlying openslide.OpenSlide handle
print(osh.level_dimensions)   # [(width0, height0), (width1, height1), ...]
print(osh.level_downsamples)  # [1.0, 4.0, 16.0, ...]

# Tile at a lower resolution for faster processing
for tile in wsi.generate_tiles(level=2, shape=256):  # level 2 (e.g. 16x downsampled)
    ...
```

**Common pyramid levels:**
- Level 0: Full resolution (e.g., 40x magnification)
- Level 1: 4x downsampled (e.g., 10x magnification)
- Level 2: 16x downsampled (e.g., 2.5x magnification)
- Level 3: 64x downsampled (thumbnail)

### Thumbnail Loading

Generate low-resolution thumbnails for visualization and quality control. PathML exposes `wsi.plot()` for a quick thumbnail view; for an array, pull one from the backend handle:

```python
import matplotlib.pyplot as plt

# Quick built-in thumbnail plot
wsi.plot()

# Or get an array via the OpenSlide backend handle
thumbnail = wsi.slide.slide.get_thumbnail((1024, 1024))
plt.imshow(thumbnail)
plt.axis("off")
plt.show()
```

## Batch Loading with SlideDataset

Process multiple slides efficiently using `SlideDataset`, which wraps a list of `SlideData` objects (not raw paths):

```python
from pathml.core import HESlide, SlideDataset
import glob

# Build slide objects, then a dataset
slide_paths = glob.glob("data/*.svs")
slides = [HESlide(p) for p in slide_paths]
dataset = SlideDataset(slides)

# Access individual slides
for slide in dataset.slides:
    print(slide.name)
```

**With preprocessing pipeline:**
```python
from pathml.preprocessing import Pipeline, StainNormalizationHE

pipeline = Pipeline([
    StainNormalizationHE(target="normalize"),
])

# Run across the dataset. PathML parallelizes with Dask via a client;
# pass a dask.distributed Client (see the distributed example below).
dataset.run(pipeline)
```

## Metadata Access

Extract slide metadata including acquisition parameters, magnification, and vendor-specific information. For OpenSlide-backed slides the vendor properties live on the backend handle:

```python
# OpenSlide properties dictionary
props = wsi.slide.slide.properties

print(props.get("openslide.objective-power"))  # Magnification
print(props.get("openslide.mpp-x"))            # Microns per pixel X
print(props.get("openslide.mpp-y"))            # Microns per pixel Y
print(props.get("openslide.vendor"))           # Scanner vendor

# Slide dimensions
print(wsi.shape)  # full-resolution dimensions
```

## Working with DICOM Slides

PathML supports DICOM WSI through specialized handling:

```python
from pathml.core import SlideData

# Load DICOM WSI (backend is the string "dicom")
dicom_slide = SlideData("path/to/slide.dcm", backend="dicom")
```

## Working with OME-TIFF

OME-TIFF provides an open standard for multi-dimensional imaging:

```python
from pathml.core import SlideData

# Load OME-TIFF via the Bio-Formats backend (requires a JDK / JVM)
ome_slide = SlideData("path/to/slide.ome.tiff", backend="bioformats")

# Inspect dimensions / channels
print(ome_slide.shape)
```

## Performance Considerations

### Memory Management

For large WSI files (often >1GB), use tile-based loading to avoid memory exhaustion:

```python
# Efficient: stream tiles one at a time
for tile in wsi.generate_tiles(level=1, shape=256):
    process_tile(tile)

# Inefficient: pulling the whole level-0 image into memory may crash
full_image = wsi.extract_region(location=(0, 0), size=wsi.shape)
```

### Distributed Processing

Use Dask for parallel processing across multiple workers:

```python
from pathml.core import HESlide, SlideDataset
from dask.distributed import Client

# Start a Dask client and pass it to run()
client = Client(n_workers=8, threads_per_worker=2)

dataset = SlideDataset([HESlide(p) for p in slide_paths])
dataset.run(pipeline, client=client)
```

### Level Selection

Balance resolution and performance by selecting appropriate pyramid levels:

- **Level 0:** Use for final analysis requiring maximum detail
- **Level 1-2:** Use for most preprocessing and model training
- **Level 3+:** Use for thumbnails, quality control, and rapid exploration

## Common Issues and Solutions

**Issue: Slide fails to load**
- Verify file format is supported
- Check file permissions and path
- Try different backend: `backend="bioformats"` or `backend="openslide"`

**Issue: Out of memory errors**
- Use tile-based loading instead of full-slide loading
- Process at lower pyramid level (e.g., level=1 or level=2)
- Reduce tile_size parameter
- Enable distributed processing with Dask

**Issue: Color inconsistencies across slides**
- Apply stain normalization preprocessing (see `preprocessing.md`)
- Check scanner metadata for calibration information
- Use `StainNormalizationHE` transform in preprocessing pipeline

**Issue: Metadata missing or incorrect**
- Different vendors store metadata in different locations
- Use `wsi.metadata` to inspect available fields
- Some formats may have limited metadata support

## Best Practices

1. **Always inspect pyramid structure** before processing: Check `level_dimensions` and `level_downsamples` to understand available resolutions

2. **Use appropriate pyramid levels**: Process at level 1-2 for most tasks; reserve level 0 for final high-resolution analysis

3. **Tile with overlap** for segmentation tasks: Use stride < tile_size to avoid edge artifacts

4. **Verify magnification consistency**: Check `openslide.objective-power` metadata when combining slides from different sources

5. **Handle vendor-specific formats**: Use specialized slide classes (CODEXSlide, VectraSlide) for multiparametric data

6. **Implement quality control**: Generate thumbnails and inspect for artifacts before processing

7. **Use distributed processing** for large datasets: Leverage Dask for parallel processing across multiple workers

## Example Workflows

### Loading and Inspecting a New Slide

```python
from pathml.core import HESlide

# Load slide
wsi = HESlide("path/to/slide.svs", name="example")

# Inspect properties
print(f"Shape: {wsi.shape}")
props = wsi.slide.slide.properties  # OpenSlide properties
print(f"Magnification: {props.get('openslide.objective-power')}")

# Quick thumbnail for QC
wsi.plot()
```

### Processing Multiple Slides

```python
from pathml.core import HESlide, SlideDataset
from pathml.preprocessing import Pipeline, TissueDetectionHE
from dask.distributed import Client
import glob

# Find all slides
slide_paths = glob.glob("data/slides/*.svs")

# Create pipeline
pipeline = Pipeline([TissueDetectionHE()])

# Build dataset from slide objects and run with a Dask client
dataset = SlideDataset([HESlide(p) for p in slide_paths])
client = Client(n_workers=8, threads_per_worker=2)
dataset.run(pipeline, client=client)

# Persist each processed slide to h5path
dataset.write("processed/")
```

### Loading CODEX Multiparametric Data

```python
from pathml.core import CODEXSlide
from pathml.preprocessing import Pipeline, CollapseRunsCODEX, SegmentMIF

# Load CODEX slide
codex = CODEXSlide("path/to/codex_dir", stain="IF")

# Create CODEX-specific pipeline. SegmentMIF channels are integer indices
# into the collapsed channel stack (see multiparametric.md).
pipeline = Pipeline([
    CollapseRunsCODEX(z=0),  # select z-plane
    SegmentMIF(
        model="mesmer",
        nuclear_channel=0,
        cytoplasm_channel=29,
    ),
])

codex.run(pipeline)
```

## Additional Resources

- **PathML Documentation:** https://pathml.readthedocs.io/
- **OpenSlide:** https://openslide.org/ (underlying library for WSI formats)
- **Bio-Formats:** https://www.openmicroscopy.org/bio-formats/ (alternative backend)
- **DICOM Standard:** https://www.dicomstandard.org/

# Machine Learning

## Overview

PathML provides comprehensive machine learning capabilities for computational pathology, including pre-built models for nucleus detection and segmentation, PyTorch-integrated training workflows, public dataset access, and ONNX-based inference deployment. The framework seamlessly bridges image preprocessing with deep learning to enable end-to-end pathology ML pipelines.

## Pre-Built Models

PathML includes state-of-the-art pre-trained models for nucleus analysis:

### HoVer-Net

**HoVer-Net** (Horizontal and Vertical Network) performs simultaneous nucleus instance segmentation and classification.

**Architecture:**
- Encoder-decoder structure with three prediction branches:
  - **Nuclear Pixel (NP)** - Binary segmentation of nuclear regions
  - **Horizontal-Vertical (HV)** - Distance maps to nucleus centroids
  - **Classification (NC)** - Nucleus type classification

**Nucleus types:**
1. Epithelial
2. Inflammatory
3. Connective/Soft tissue
4. Dead/Necrotic
5. Background

**Usage:**
```python
from pathml.ml import HoVerNet
import torch

# Construct the model (n_classes = number of nucleus types).
model = HoVerNet(n_classes=5)

# Move to GPU if available
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = model.to(device)

# Inference on a tile. HoVerNet returns a list/tuple of branch outputs
# (NP, HV, and the classification branch when n_classes is set).
tile_image = torch.from_numpy(tile.image).permute(2, 0, 1).unsqueeze(0).float()
tile_image = tile_image.to(device)

with torch.no_grad():
    outputs = model(tile_image)
```

For pretrained weights, download a published HoVer-Net checkpoint and load it with `model.load_state_dict(...)`; the constructor does not fetch weights itself. Confirm the exact constructor argument (`n_classes`) and the branch output layout against the API docs for your installed version.

**Post-processing** (the real function is `post_process_batch_hovernet`, not `hovernet_postprocess`):
```python
from pathml.ml import post_process_batch_hovernet

# Convert a batch of model outputs to instance + type maps
instance_map, type_map = post_process_batch_hovernet(outputs, n_classes=5)

# instance_map: each nucleus has a unique ID
# type_map: each nucleus assigned a type
```

### HACTNet

**HACTNet** is a graph-neural-network model that operates on hierarchical cell-graph + tissue-graph (HACT) representations of a tissue region for region/graph-level classification. It is not a per-pixel nucleus classifier; pair it with the graph builders in `pathml.graph` (see `graphs.md`).

```python
from pathml.ml import HACTNet

# HACTNet consumes batched cell-graph and tissue-graph inputs constructed
# from segmented tissue. Check the API docs for the exact constructor
# arguments and the expected graph batch format for your version.
model = HACTNet(...)
```

## Training Workflows

### Dataset Preparation

PathML provides PyTorch-compatible dataset classes in `pathml.datasets` (note: NOT `pathml.ml`):

**TileDataset** wraps an on-disk h5path slide so its tiles can be served to a PyTorch `DataLoader`:
```python
from pathml.datasets import TileDataset
from torch.utils.data import DataLoader

# Point at an h5path file written by SlideData.write(...)
tile_dataset = TileDataset("processed/slide001.h5path")

loader = DataLoader(tile_dataset, batch_size=32, shuffle=True, num_workers=4)
for images, masks, labels in loader:
    ...
```

`pathml.datasets` also provides `EntityDataset` (for graph/entity inputs) and ready-made downloadable DataModules, `PanNukeDataModule` and `DeepFocusDataModule`. There is no `PathMLDataModule`; build train/val/test `DataLoader`s yourself or use one of the provided DataModules. Confirm the `TileDataset.__getitem__` return tuple against your version.

### Training HoVer-Net

Complete workflow for training HoVer-Net on custom data:

PathML ships the composite HoVer-Net loss as `loss_hovernet` (combining the NP, HV, and classification branch terms) — use it rather than reimplementing the loss by hand.

```python
import torch
from torch.utils.data import DataLoader
from pathml.ml import HoVerNet, loss_hovernet, post_process_batch_hovernet
from pathml.datasets import PanNukeDataModule

# 1. Prepare data (PanNukeDataModule lives in pathml.datasets)
data_module = PanNukeDataModule(
    data_dir="path/to/pannuke",
    batch_size=8,
    nucleus_type_labels=True,
)
train_loader = data_module.train_dataloader

# 2. Initialize model
model = HoVerNet(n_classes=5)

# 3. Use PathML's HoVer-Net loss
def criterion(outputs, ground_truth, n_classes=5):
    return loss_hovernet(outputs, ground_truth, n_classes=n_classes)

# 4. Configure optimizer
optimizer = torch.optim.Adam(model.parameters(), lr=1e-4, weight_decay=1e-5)

scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(
    optimizer,
    mode='min',
    factor=0.5,
    patience=10
)

# 5. Training loop.
# PanNukeDataModule yields (images, masks, tissue_types) batches; the
# ground-truth tensors required by loss_hovernet (NP / HV / NC targets)
# are derived from the masks. Consult the PanNuke example in the PathML
# docs for the exact unpacking for your version.
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
model = model.to(device)

num_epochs = 100
for epoch in range(num_epochs):
    model.train()
    train_loss = 0.0
    n_batches = 0

    for images, masks, tissue_type in train_loader:
        images = images.float().to(device)
        ground_truth = masks.to(device)  # see docs for target construction

        optimizer.zero_grad()
        outputs = model(images)
        loss = loss_hovernet(outputs, ground_truth, n_classes=5)
        loss.backward()
        optimizer.step()

        train_loss += loss.item()
        n_batches += 1

    scheduler.step(train_loss)
    print(f"Epoch {epoch+1}/{num_epochs}  Train Loss: {train_loss / max(n_batches, 1):.4f}")

    if (epoch + 1) % 10 == 0:
        torch.save(model.state_dict(), f"hovernet_epoch_{epoch+1}.pth")
```

## Public Datasets

PathML bundles two downloadable DataModules in `pathml.datasets`: `PanNukeDataModule` and `DeepFocusDataModule`. (There is no built-in TCGA DataModule — download TCGA WSIs via the GDC portal and wrap them with `SlideDataset` / `TileDataset`.)

### PanNuke Dataset

**PanNuke** is a multi-tissue histology dataset with nucleus instance + type annotations across several cell categories, commonly used to train HoVer-Net.

```python
from pathml.datasets import PanNukeDataModule

# PanNukeDataModule will download the data to data_dir on first use.
pannuke = PanNukeDataModule(
    data_dir="path/to/pannuke",
    batch_size=16,
    nucleus_type_labels=True,  # include per-nucleus type labels
)

# DataLoaders are exposed as properties (not methods)
train_loader = pannuke.train_dataloader
valid_loader = pannuke.valid_dataloader
test_loader = pannuke.test_dataloader

# Each batch is a tuple: (images, masks, tissue_type)
for images, masks, tissue_type in train_loader:
    ...
```

Confirm the exact constructor arguments and batch tuple layout against the API docs for your installed version.

### Custom Dataset Integration

Create custom datasets for PathML workflows:

```python
from torch.utils.data import Dataset
import numpy as np
from pathlib import Path

class CustomPathologyDataset(Dataset):
    def __init__(self, data_dir, transform=None):
        self.data_dir = Path(data_dir)
        self.image_paths = list(self.data_dir.glob('images/*.png'))
        self.transform = transform

    def __len__(self):
        return len(self.image_paths)

    def __getitem__(self, idx):
        # Load image
        image_path = self.image_paths[idx]
        image = np.array(Image.open(image_path))

        # Load corresponding annotation
        annot_path = self.data_dir / 'annotations' / f'{image_path.stem}.npy'
        annotation = np.load(annot_path)

        # Apply transforms
        if self.transform:
            image = self.transform(image)

        return {
            'image': torch.from_numpy(image).permute(2, 0, 1).float(),
            'annotation': torch.from_numpy(annotation).long(),
            'path': str(image_path)
        }

# Use in PathML workflow
dataset = CustomPathologyDataset('path/to/data')
dataloader = DataLoader(dataset, batch_size=16, shuffle=True, num_workers=4)
```

## Data Augmentation

Apply augmentations to improve model generalization:

```python
import albumentations as A
from albumentations.pytorch import ToTensorV2

# Define augmentation pipeline
train_transform = A.Compose([
    A.RandomRotate90(p=0.5),
    A.Flip(p=0.5),
    A.ColorJitter(brightness=0.2, contrast=0.2, saturation=0.2, hue=0.1, p=0.5),
    A.GaussianBlur(blur_limit=(3, 7), p=0.3),
    A.ElasticTransform(alpha=1, sigma=50, alpha_affine=50, p=0.3),
    A.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),
    ToTensorV2()
])

val_transform = A.Compose([
    A.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),
    ToTensorV2()
])

# Apply augmentations inside your Dataset's __getitem__ (e.g. the
# CustomPathologyDataset above), or wrap a TileDataset and transform each
# tile before returning it. Pass the transform to whichever Dataset you use.
```

## Model Evaluation

### Metrics

Standard nucleus-segmentation metrics — Dice, Aggregated Jaccard Index (AJI), and Panoptic Quality (PQ = SQ x RQ) — are not bundled as a `pathml.ml.metrics` module. Use a published implementation (e.g. the HoVer-Net authors' `compute_stats` code, or a maintained instance-segmentation metrics package) and feed it the post-processed instance + type maps.

### Evaluation Loop

```python
import torch
from pathml.ml import post_process_batch_hovernet

model.eval()
all_pred_inst, all_pred_type = [], []

with torch.no_grad():
    for images, masks, tissue_type in test_loader:
        images = images.float().to(device)
        outputs = model(images)

        # Post-process the whole batch at once
        inst_map, type_map = post_process_batch_hovernet(outputs, n_classes=5)
        all_pred_inst.append(inst_map)
        all_pred_type.append(type_map)

# Pass all_pred_inst/all_pred_type plus the ground-truth instance/type maps
# to your chosen AJI / PQ implementation.
```

## ONNX Inference

Deploy models using ONNX for production inference:

### Export to ONNX

```python
import torch
from pathml.ml import HoVerNet

# Load trained model and weights
model = HoVerNet(n_classes=5)
model.load_state_dict(torch.load("hovernet_epoch_100.pth"))
model.eval()

# Create dummy input
dummy_input = torch.randn(1, 3, 256, 256)

# Export to ONNX
torch.onnx.export(
    model,
    dummy_input,
    'hovernet_model.onnx',
    export_params=True,
    opset_version=11,
    input_names=['input'],
    output_names=['np_output', 'hv_output', 'nc_output'],
    dynamic_axes={
        'input': {0: 'batch_size'},
        'np_output': {0: 'batch_size'},
        'hv_output': {0: 'batch_size'},
        'nc_output': {0: 'batch_size'}
    }
)
```

### ONNX Runtime Inference

```python
import onnxruntime as ort
import numpy as np

# Load ONNX model
session = ort.InferenceSession('hovernet_model.onnx')

# Prepare input
input_name = session.get_inputs()[0].name
tile_image = preprocess_tile(tile)  # Normalize, transpose to (1, 3, H, W)

# Run inference
outputs = session.run(None, {input_name: tile_image})

# Post-process with PathML's batch post-processor (rebuild the tensor
# structure it expects from the ONNX outputs first)
from pathml.ml import post_process_batch_hovernet
inst_map, type_map = post_process_batch_hovernet(outputs, n_classes=5)
```

### Batch Inference Pipeline

```python
from pathml.core import HESlide
from pathml.ml import post_process_batch_hovernet
import onnxruntime as ort

def run_onnx_inference_pipeline(slide_path, onnx_model_path):
    wsi = HESlide(slide_path)
    session = ort.InferenceSession(onnx_model_path)
    input_name = session.get_inputs()[0].name

    results = []
    for tile in wsi.generate_tiles(level=1, shape=256, stride=256):
        tile_array = preprocess_tile(tile.image)  # normalize -> (1, 3, H, W)
        outputs = session.run(None, {input_name: tile_array})
        inst_map, type_map = post_process_batch_hovernet(outputs, n_classes=5)
        results.append({
            "coords": tile.coords,
            "instance_map": inst_map,
            "type_map": type_map,
        })
    return results

results = run_onnx_inference_pipeline("slide.svs", "hovernet_model.onnx")
```

## Transfer Learning

Fine-tune pre-trained models on custom datasets:

```python
from pathml.ml import HoVerNet

# Load a model and a pretrained checkpoint
model = HoVerNet(n_classes=5)
model.load_state_dict(torch.load("pretrained_hovernet.pth"))

# Freeze encoder layers for initial training
for name, param in model.named_parameters():
    if 'encoder' in name:
        param.requires_grad = False

# Fine-tune only decoder and classification heads
optimizer = torch.optim.Adam(
    filter(lambda p: p.requires_grad, model.parameters()),
    lr=1e-4
)

# Train for a few epochs
train_for_n_epochs(model, train_loader, optimizer, num_epochs=10)

# Unfreeze all layers for full fine-tuning
for param in model.parameters():
    param.requires_grad = True

# Continue training with lower learning rate
optimizer = torch.optim.Adam(model.parameters(), lr=1e-5)
train_for_n_epochs(model, train_loader, optimizer, num_epochs=50)
```

## Best Practices

1. **Use pre-trained weights when available:**
   - Load a published checkpoint with `model.load_state_dict(...)` for better initialization
   - Fine-tune on domain-specific data

2. **Apply appropriate data augmentation:**
   - Rotate, flip for orientation invariance
   - Color jitter to handle staining variations
   - Elastic deformation for biological variability

3. **Monitor multiple metrics:**
   - Track detection, segmentation, and classification separately
   - Use domain-specific metrics (AJI, PQ) beyond standard accuracy

4. **Handle class imbalance:**
   - Weighted loss functions for rare cell types
   - Oversampling minority classes
   - Focal loss for hard examples

5. **Validate on diverse tissue types:**
   - Ensure generalization across different tissues
   - Test on held-out anatomical sites

6. **Optimize for inference:**
   - Export to ONNX for faster deployment
   - Batch tiles for efficient GPU utilization
   - Use mixed precision (FP16) when possible

7. **Save checkpoints regularly:**
   - Keep best model based on validation metrics
   - Save optimizer state for training resumption

## Common Issues and Solutions

**Issue: Poor segmentation at nucleus boundaries**
- Use HV maps (horizontal-vertical) to separate touching nuclei
- Increase weight of HV loss term
- Apply morphological post-processing

**Issue: Misclassification of similar cell types**
- Increase classification loss weight
- Add hierarchical classification (HACTNet)
- Augment training data for confused classes

**Issue: Training unstable or not converging**
- Reduce learning rate
- Use gradient clipping: `torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)`
- Check for data preprocessing issues

**Issue: Out of memory during training**
- Reduce batch size
- Use gradient accumulation
- Enable mixed precision training: `torch.cuda.amp`

**Issue: Model overfits to training data**
- Increase data augmentation
- Add dropout layers
- Reduce model capacity
- Use early stopping based on validation loss

## Additional Resources

- **PathML ML API:** https://pathml.readthedocs.io/en/latest/api_ml_reference.html
- **HoVer-Net Paper:** Graham et al., "HoVer-Net: Simultaneous Segmentation and Classification of Nuclei in Multi-Tissue Histology Images," Medical Image Analysis, 2019
- **PanNuke Dataset:** https://warwick.ac.uk/fac/cross_fac/tia/data/pannuke
- **PyTorch Lightning:** https://www.pytorchlightning.ai/
- **ONNX Runtime:** https://onnxruntime.ai/

# Multiparametric Imaging

## Overview

PathML provides specialized support for multiparametric imaging technologies that simultaneously measure multiple markers at single-cell resolution. These techniques include CODEX, Vectra multiplex immunofluorescence, MERFISH, and other spatial proteomics and transcriptomics platforms. PathML handles the unique data structures, processing requirements, and quantification workflows specific to each technology.

## Supported Technologies

### CODEX (CO-Detection by indEXing)
- Cyclic immunofluorescence imaging
- 40+ protein markers simultaneously
- Single-cell spatial proteomics
- Multi-cycle acquisition with antibody barcoding

### Vectra Polaris
- Multispectral multiplex immunofluorescence
- 6-8 markers per slide
- Spectral unmixing
- Whole-slide scanning

### MERFISH (Multiplexed Error-Robust FISH)
- Spatial transcriptomics
- 100s-1000s of genes
- Single-molecule resolution
- Error-correcting barcodes

### Other Platforms
- CycIF (Cyclic Immunofluorescence)
- IMC (Imaging Mass Cytometry)
- MIBI (Multiplexed Ion Beam Imaging)

## CODEX Workflows

### Loading CODEX Data

CODEX data is typically organized in multi-channel image stacks from multiple acquisition cycles:

```python
from pathml.core import CODEXSlide

# Load CODEX dataset (uses the Bio-Formats backend internally)
codex_slide = CODEXSlide("path/to/codex_directory", stain="IF")

# Inspect dimensions
print(f"Image shape: {codex_slide.shape}")
```

**CODEX directory structure:**
```
codex_directory/
├── cyc001_reg001/
│   ├── 1_00001_Z001_CH1.tif
│   ├── 1_00001_Z001_CH2.tif
│   └── ...
├── cyc002_reg001/
│   └── ...
└── channelnames.txt
```

### CODEX Preprocessing Pipeline

Complete pipeline for CODEX data processing:

```python
from pathml.preprocessing import Pipeline, CollapseRunsCODEX, SegmentMIF, QuantifyMIF

# Create CODEX-specific pipeline
codex_pipeline = Pipeline([
    # 1. Collapse the z-stack to a single multi-channel image
    CollapseRunsCODEX(z=0),

    # 2. Cell segmentation using Mesmer (integer channel indices)
    SegmentMIF(
        model="mesmer",
        nuclear_channel=0,        # e.g. DAPI index
        cytoplasm_channel=29,     # e.g. membrane/cytoplasm marker index
        image_resolution=0.377,   # microns per pixel
    ),

    # 3. Quantify per-cell marker expression across all channels
    QuantifyMIF(segmentation_mask="cell_segmentation"),
])

# Run pipeline on the slide
codex_slide.run(codex_pipeline)

# Access results
segmentation_mask = codex_slide.masks["cell_segmentation"]
cell_data = codex_slide.counts  # AnnData object (cells x markers)
```

### CollapseRunsCODEX

Consolidates multi-cycle CODEX acquisitions into a single multi-channel image:

```python
from pathml.preprocessing import CollapseRunsCODEX

transform = CollapseRunsCODEX(z=0)  # select the focal plane (0-indexed)
```

**Parameters:**
- `z`: which z-plane to extract from the CODEX z-stack (0-indexed)

**Output:** Single multi-channel image with all markers stacked along the channel axis.

### Cell Segmentation with Mesmer

DeepCell Mesmer provides accurate cell segmentation for multiparametric imaging:

```python
from pathml.preprocessing import SegmentMIF

transform = SegmentMIF(
    model="mesmer",            # DeepCell Mesmer model
    nuclear_channel=0,         # integer index of the nuclear marker (e.g. DAPI)
    cytoplasm_channel=29,      # integer index of a membrane/cytoplasm marker
    image_resolution=0.377,    # microns per pixel (important for accuracy)
)
```

`SegmentMIF` produces both `nuclear_segmentation` and `cell_segmentation` (whole-cell) masks. Tune behavior via `preprocess_kwargs`, `postprocess_kwargs_nuclear`, and `postprocess_kwargs_whole_cell`; there is no `compartment`/`min_cell_size`/`max_cell_size` argument.

**Choosing the cytoplasm channel** (pick the index of a marker that outlines cells):
- A pan-leukocyte marker (e.g. CD45) for immune-rich tissues
- A pan-cytokeratin marker (e.g. panCK) for epithelial tissues
- A universal membrane marker, or an averaged membrane channel

### Remote Segmentation

Use the DeepCell service for segmentation without a local GPU:

```python
from pathml.preprocessing import SegmentMIFRemote

transform = SegmentMIFRemote(
    model="mesmer",
    nuclear_channel=0,
    cytoplasm_channel=29,
)
```

### Marker Quantification

Extract single-cell marker expression from segmented images:

```python
from pathml.preprocessing import QuantifyMIF

# Quantify per-cell expression across all channels of the segmented mask.
transform = QuantifyMIF(segmentation_mask="cell_segmentation")
```

**Output:** an AnnData object written to `slide.counts`:
- `adata.X`: marker expression matrix (cells x channels)
- `adata.obs`: per-cell metadata (coordinates, area, etc.)
- `adata.var`: channel/marker metadata
- `adata.obsm["spatial"]`: cell centroid coordinates

### Integration with AnnData

Process multiple CODEX slides into unified AnnData object:

```python
from pathml.core import CODEXSlide, SlideDataset
from dask.distributed import Client
import anndata as ad

# Process multiple slides
slide_paths = ["slide1", "slide2", "slide3"]
dataset = SlideDataset([CODEXSlide(p, stain="IF") for p in slide_paths])

client = Client(n_workers=8, threads_per_worker=2)
dataset.run(codex_pipeline, client=client)

# Concatenate the per-slide AnnData objects (each lives on slide.counts)
combined_adata = ad.concat(
    [s.counts for s in dataset.slides],
    join="outer",
    label="Region",
    keys=slide_paths,
    index_unique="_",
)

# Save for downstream analysis
combined_adata.write("codex_dataset.h5ad")
```

## Vectra Workflows

### Loading Vectra Data

Vectra stores data in proprietary `.qptiff` format:

```python
from pathml.core import VectraSlide

# Load Vectra slide (.qptiff); VectraSlide uses the Bio-Formats backend
vectra_slide = VectraSlide("path/to/slide.qptiff")

print(f"Shape: {vectra_slide.shape}")
```

### Vectra Preprocessing

```python
from pathml.preprocessing import Pipeline, CollapseRunsVectra, SegmentMIF, QuantifyMIF

vectra_pipeline = Pipeline([
    # 1. Collapse the Vectra multi-channel data into a single image
    CollapseRunsVectra(),

    # 2. Cell segmentation (integer channel indices)
    SegmentMIF(
        model="mesmer",
        nuclear_channel=0,
        cytoplasm_channel=4,
        image_resolution=0.5,
    ),

    # 3. Quantification -> AnnData on slide.counts
    QuantifyMIF(segmentation_mask="cell_segmentation"),
])

vectra_slide.run(vectra_pipeline)
```

## Downstream Analysis

### Cell Type Annotation

Annotate cells based on marker expression:

```python
import anndata as ad
import numpy as np

# Load quantified data
adata = ad.read_h5ad('codex_dataset.h5ad')

# Define cell types by marker thresholds
def annotate_cell_types(adata, thresholds):
    cell_types = np.full(adata.n_obs, 'Unknown', dtype=object)

    # T cells: CD3+
    cd3_pos = adata[:, 'CD3'].X.flatten() > thresholds['CD3']
    cell_types[cd3_pos] = 'T cell'

    # CD4 T cells: CD3+ CD4+ CD8-
    cd4_tcells = (
        (adata[:, 'CD3'].X.flatten() > thresholds['CD3']) &
        (adata[:, 'CD4'].X.flatten() > thresholds['CD4']) &
        (adata[:, 'CD8'].X.flatten() < thresholds['CD8'])
    )
    cell_types[cd4_tcells] = 'CD4 T cell'

    # CD8 T cells: CD3+ CD8+ CD4-
    cd8_tcells = (
        (adata[:, 'CD3'].X.flatten() > thresholds['CD3']) &
        (adata[:, 'CD8'].X.flatten() > thresholds['CD8']) &
        (adata[:, 'CD4'].X.flatten() < thresholds['CD4'])
    )
    cell_types[cd8_tcells] = 'CD8 T cell'

    # B cells: CD20+
    b_cells = adata[:, 'CD20'].X.flatten() > thresholds['CD20']
    cell_types[b_cells] = 'B cell'

    # Macrophages: CD68+
    macrophages = adata[:, 'CD68'].X.flatten() > thresholds['CD68']
    cell_types[macrophages] = 'Macrophage'

    # Tumor cells: panCK+
    tumor = adata[:, 'panCK'].X.flatten() > thresholds['panCK']
    cell_types[tumor] = 'Tumor'

    return cell_types

# Apply annotation
thresholds = {
    'CD3': 0.5,
    'CD4': 0.4,
    'CD8': 0.4,
    'CD20': 0.3,
    'CD68': 0.3,
    'panCK': 0.5
}

adata.obs['cell_type'] = annotate_cell_types(adata, thresholds)

# Visualize cell type composition
import matplotlib.pyplot as plt
cell_type_counts = adata.obs['cell_type'].value_counts()
plt.figure(figsize=(10, 6))
cell_type_counts.plot(kind='bar')
plt.xlabel('Cell Type')
plt.ylabel('Count')
plt.title('Cell Type Composition')
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()
```

### Clustering

Unsupervised clustering to identify cell populations:

```python
import scanpy as sc

# Preprocessing for clustering
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.scale(adata, max_value=10)

# PCA
sc.tl.pca(adata, n_comps=50)

# Neighborhood graph
sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30)

# UMAP embedding
sc.tl.umap(adata)

# Leiden clustering
sc.tl.leiden(adata, resolution=0.5)

# Visualize
sc.pl.umap(adata, color=['leiden', 'CD3', 'CD8', 'CD20', 'panCK'])
```

### Spatial Visualization

Visualize cells in spatial context:

```python
import matplotlib.pyplot as plt

# Spatial scatter plot
fig, ax = plt.subplots(figsize=(15, 15))

# Color by cell type
cell_types = adata.obs['cell_type'].unique()
colors = plt.cm.tab10(np.linspace(0, 1, len(cell_types)))

for i, cell_type in enumerate(cell_types):
    mask = adata.obs['cell_type'] == cell_type
    coords = adata.obsm['spatial'][mask]
    ax.scatter(
        coords[:, 0],
        coords[:, 1],
        c=[colors[i]],
        label=cell_type,
        s=5,
        alpha=0.7
    )

ax.legend(markerscale=2)
ax.set_xlabel('X (pixels)')
ax.set_ylabel('Y (pixels)')
ax.set_title('Spatial Cell Type Distribution')
ax.axis('equal')
plt.tight_layout()
plt.show()
```

### Spatial Neighborhood Analysis

Analyze cell neighborhoods and interactions:

```python
import squidpy as sq

# Calculate spatial neighborhood enrichment
sq.gr.spatial_neighbors(adata, coord_type='generic', spatial_key='spatial')

# Neighborhood enrichment test
sq.gr.nhood_enrichment(adata, cluster_key='cell_type')

# Visualize interaction matrix
sq.pl.nhood_enrichment(adata, cluster_key='cell_type')

# Co-occurrence score
sq.gr.co_occurrence(adata, cluster_key='cell_type')
sq.pl.co_occurrence(
    adata,
    cluster_key='cell_type',
    clusters=['CD8 T cell', 'Tumor'],
    figsize=(8, 8)
)
```

### Spatial Autocorrelation

Test for spatial clustering of markers:

```python
# Moran's I spatial autocorrelation
sq.gr.spatial_autocorr(
    adata,
    mode='moran',
    genes=['CD3', 'CD8', 'PD1', 'PDL1', 'panCK']
)

# Visualize
results = adata.uns['moranI']
print(results.head())
```

## MERFISH and Other Platforms

PathML's most fully supported multiparametric workflows are CODEX and Vectra. MERFISH and related spatial-transcriptomics platforms do not have dedicated decoding/transcript-assignment transforms in the core `pathml.preprocessing` API at the time of writing, so do NOT assume classes like `MERFISHSlide`, `DecodeMERFISH`, or `AssignTranscripts` exist — verify against your installed version's API docs first.

A practical pattern for MERFISH-style data:
- Load the multi-channel image as a `MultiparametricSlide` (Bio-Formats backend).
- Segment cells with `SegmentMIF` (Mesmer) on the nuclear + a boundary channel (e.g. poly(T)).
- For barcode decoding and transcript-to-cell assignment, use a dedicated spatial-transcriptomics toolkit (e.g. a vendor pipeline or `squidpy`/`starfish`) and bring the resulting cell-by-gene matrix back into AnnData for downstream analysis.

## Quality Control

### Segmentation Quality

Inspect the instance segmentation mask directly with scikit-image (no PathML-specific QC helper is needed):

```python
import numpy as np
import matplotlib.pyplot as plt
from skimage.measure import regionprops

# segmentation_mask: integer label image (each cell a unique id)
props = regionprops(segmentation_mask)
cell_sizes = np.array([p.area for p in props])

print(f"Total cells: {len(cell_sizes)}")
print(f"Mean cell size: {cell_sizes.mean():.1f} pixels")

plt.hist(cell_sizes, bins=50)
plt.xlabel("Cell Size (pixels)")
plt.ylabel("Frequency")
plt.title("Cell Size Distribution")
plt.show()
```

### Marker Expression QC

```python
import scanpy as sc

# Load AnnData
adata = ad.read_h5ad('codex_dataset.h5ad')

# Calculate QC metrics
adata.obs['total_intensity'] = adata.X.sum(axis=1)
adata.obs['n_markers_detected'] = (adata.X > 0).sum(axis=1)

# Filter low-quality cells
adata = adata[adata.obs['total_intensity'] > 100, :]
adata = adata[adata.obs['n_markers_detected'] >= 3, :]

# Visualize
sc.pl.violin(adata, ['total_intensity', 'n_markers_detected'], multi_panel=True)
```

## Batch Processing

Process large multiparametric datasets efficiently:

```python
from pathml.core import SlideDataset
from pathml.preprocessing import Pipeline
from dask.distributed import Client
import glob

# Start Dask cluster
client = Client(n_workers=16, threads_per_worker=2, memory_limit='8GB')

# Find all CODEX slides
slide_dirs = glob.glob('data/codex_slides/*/')

# Create dataset
codex_slides = [CODEXSlide(d, stain='IF') for d in slide_dirs]
dataset = SlideDataset(codex_slides)

# Run pipeline in parallel (pass the Dask client to run())
dataset.run(codex_pipeline, client=client)

# Save the per-slide AnnData (each on slide.counts)
for i, slide in enumerate(dataset.slides):
    slide.counts.write(f"processed/slide_{i}.h5ad")

client.close()
```

## Integration with Other Tools

### Export to Spatial Analysis Tools

```python
# Export to Giotto
def export_to_giotto(adata, output_dir):
    import os
    os.makedirs(output_dir, exist_ok=True)

    # Expression matrix
    pd.DataFrame(
        adata.X.T,
        index=adata.var_names,
        columns=adata.obs_names
    ).to_csv(f'{output_dir}/expression.csv')

    # Cell coordinates
    pd.DataFrame(
        adata.obsm['spatial'],
        columns=['x', 'y'],
        index=adata.obs_names
    ).to_csv(f'{output_dir}/spatial_locs.csv')

# Export to Seurat
def export_to_seurat(adata, output_file):
    adata.write_h5ad(output_file)
    # Read in R with: library(Seurat); ReadH5AD(output_file)
```

## Best Practices

1. **Channel selection for segmentation:**
   - Use brightest, most consistent nuclear marker (usually DAPI)
   - Choose membrane/cytoplasm marker based on tissue type
   - Test multiple options to optimize segmentation

2. **Background subtraction:**
   - Apply before quantification to reduce autofluorescence
   - Use blank/control images to model background

3. **Quality control:**
   - Visualize segmentation on sample regions
   - Check cell size distributions for outliers
   - Validate marker expression ranges

4. **Cell type annotation:**
   - Start with canonical markers (CD3, CD20, panCK)
   - Use multiple markers for robust classification
   - Consider unsupervised clustering to discover populations

5. **Spatial analysis:**
   - Account for tissue architecture (epithelium, stroma, etc.)
   - Consider local density when interpreting interactions
   - Use permutation tests for statistical significance

6. **Batch effects:**
   - Include batch information in AnnData.obs
   - Apply batch correction if combining multiple experiments
   - Visualize batch effects with UMAP colored by batch

## Common Issues and Solutions

**Issue: Poor segmentation quality**
- Verify nuclear and cytoplasm channels are correctly specified
- Adjust image_resolution parameter to match actual resolution
- Try different cytoplasm markers
- Manually tune min/max cell size parameters

**Issue: Low marker intensity**
- Check for background subtraction artifacts
- Verify channel names match actual channels
- Inspect raw images for technical issues (focus, exposure)

**Issue: Cell type annotations don't match expectations**
- Adjust marker thresholds (too high/low)
- Visualize marker distributions to set data-driven thresholds
- Check for antibody specificity issues

**Issue: Spatial analysis shows no significant interactions**
- Increase neighborhood radius
- Check for sufficient cell numbers per type
- Verify spatial coordinates are correctly scaled

## Additional Resources

- **PathML Multiparametric API:** https://pathml.readthedocs.io/en/latest/api_preprocessing_reference.html
- **CODEX:** https://www.akoyabio.com/codex/
- **Vectra:** https://www.akoyabio.com/phenoimager/instruments/vectra-3-0/
- **DeepCell Mesmer:** https://www.deepcell.org/
- **Scanpy:** https://scanpy.readthedocs.io/ (single-cell analysis)
- **Squidpy:** https://squidpy.readthedocs.io/ (spatial omics analysis)

# Preprocessing Pipelines & Transforms

## Overview

PathML provides a modular preprocessing architecture based on composable transforms organized into pipelines. Transforms are individual operations that modify images, create masks, or extract features. Pipelines chain transforms together to create reproducible, scalable preprocessing workflows for computational pathology.

## Pipeline Architecture

### Pipeline Class

The `Pipeline` class composes a sequence of transforms applied consecutively:

A `Pipeline` is a reusable, serializable list of transforms. You apply it by passing it to a slide's or dataset's `run()` method (the slide/dataset drives execution, not the pipeline):

```python
from pathml.preprocessing import Pipeline, Transform1, Transform2

# Create pipeline
pipeline = Pipeline([
    Transform1(param1=value1),
    Transform2(param2=value2),
    # ... more transforms
])

# Run on a single slide
slide_data.run(pipeline)

# Run on a dataset (parallelized via a Dask client)
dataset.run(pipeline, client=client)
```

**Key features:**
- Sequential execution of transforms
- Automatic handling of tiles and masks
- Distributed processing support with Dask
- Reproducible workflows with serializable configuration

### Transform Base Class

All transforms inherit from the `Transform` base class and implement:
- `apply()` - Core transformation logic
- `input_type` - Expected input (tile, mask, etc.)
- `output_type` - Produced output

## Transform Categories

PathML provides transforms in six major categories:

1. **Image Modification** - Blur, rescale, histogram equalization
2. **Mask Creation** - Tissue detection, nucleus detection, thresholding
3. **Mask Modification** - Morphological operations on masks
4. **Stain Processing** - H&E stain normalization and separation
5. **Quality Control** - Artifact detection, white space labeling
6. **Specialized** - Multiparametric imaging, cell segmentation

## Image Modification Transforms

### Blur Operations

Apply various blurring kernels for noise reduction:

**MedianBlur:**
```python
from pathml.preprocessing import MedianBlur

# Apply median filter
transform = MedianBlur(kernel_size=5)
```
- Effective for salt-and-pepper noise
- Preserves edges better than Gaussian blur

**GaussianBlur:**
```python
from pathml.preprocessing import GaussianBlur

# Apply Gaussian blur
transform = GaussianBlur(kernel_size=5, sigma=1.0)
```
- Smooth noise reduction
- Adjustable sigma controls blur strength

**BoxBlur:**
```python
from pathml.preprocessing import BoxBlur

# Apply box filter
transform = BoxBlur(kernel_size=5)
```
- Fastest blur operation
- Uniform averaging within kernel

### Intensity Adjustments

**RescaleIntensity:**
```python
from pathml.preprocessing import RescaleIntensity

# Rescale intensity to [0, 255]
transform = RescaleIntensity(
    in_range=(0, 1.0),
    out_range=(0, 255)
)
```

**HistogramEqualization:**
```python
from pathml.preprocessing import HistogramEqualization

# Global histogram equalization
transform = HistogramEqualization()
```
- Enhances global contrast
- Spreads out intensity distribution

**AdaptiveHistogramEqualization (CLAHE):**
```python
from pathml.preprocessing import AdaptiveHistogramEqualization

# Contrast Limited Adaptive Histogram Equalization
transform = AdaptiveHistogramEqualization(
    clip_limit=0.03,
    tile_grid_size=(8, 8)
)
```
- Enhances local contrast
- Prevents over-amplification with clip_limit
- Better for images with varying local contrast

### Superpixel Processing

**SuperpixelInterpolation:**
```python
from pathml.preprocessing import SuperpixelInterpolation

# Divide into superpixels using SLIC
transform = SuperpixelInterpolation(
    n_segments=100,
    compactness=10.0
)
```
- Segments image into perceptually meaningful regions
- Useful for feature extraction and segmentation

## Mask Creation Transforms

### H&E Tissue and Nucleus Detection

**TissueDetectionHE:**
```python
from pathml.preprocessing import TissueDetectionHE

# Detect tissue regions in H&E slides
transform = TissueDetectionHE(
    use_saturation=True,  # Use HSV saturation channel
    threshold=10,  # Intensity threshold
    min_region_size=500  # Minimum tissue region size in pixels
)
```
- Creates binary tissue mask
- Filters small regions and artifacts
- Stores mask in `tile.masks['tissue']`

**NucleusDetectionHE:**
```python
from pathml.preprocessing import NucleusDetectionHE

# Detect nuclei in H&E images
transform = NucleusDetectionHE(
    stain='hematoxylin',  # Use hematoxylin channel
    threshold=0.3,
    min_nucleus_size=10
)
```
- Separates hematoxylin stain
- Thresholds to create nucleus mask
- Stores mask in `tile.masks['nucleus']`

### Binary Thresholding

**BinaryThreshold:**
```python
from pathml.preprocessing import BinaryThreshold

# Threshold using Otsu's method
transform = BinaryThreshold(
    method='otsu',  # 'otsu' or manual threshold value
    invert=False
)

# Or specify manual threshold
transform = BinaryThreshold(threshold=128)
```

### Foreground Detection

**ForegroundDetection:**
```python
from pathml.preprocessing import ForegroundDetection

# Detect foreground regions
transform = ForegroundDetection(
    threshold=0.5,
    min_region_size=1000,  # Minimum size in pixels
    use_saturation=True
)
```

## Mask Modification Transforms

Apply morphological operations to clean up masks:

**MorphOpen:**
```python
from pathml.preprocessing import MorphOpen

# Remove small objects and noise
transform = MorphOpen(
    kernel_size=5,
    mask_name='tissue'  # Which mask to modify
)
```
- Erosion followed by dilation
- Removes small objects and noise

**MorphClose:**
```python
from pathml.preprocessing import MorphClose

# Fill small holes
transform = MorphClose(
    kernel_size=5,
    mask_name='tissue'
)
```
- Dilation followed by erosion
- Fills small holes in mask

## Stain Normalization

### StainNormalizationHE

Normalize H&E staining across slides to account for variations in staining procedure and scanners:

```python
from pathml.preprocessing import StainNormalizationHE

# Normalize to reference stain vectors
transform = StainNormalizationHE(
    target="normalize",  # "normalize", "hematoxylin", or "eosin"
    stain_estimation_method="macenko",  # "macenko" or "vahadane"
)
```

**Target modes:**
- `"normalize"` - Normalize both stains to reference
- `"hematoxylin"` - Extract hematoxylin channel only
- `"eosin"` - Extract eosin channel only

**Stain estimation methods:**
- `"macenko"` - Macenko et al. 2009 method (faster, more stable)
- `"vahadane"` - Vahadane et al. 2016 method (more accurate, slower)

**Other documented parameters** (confirm against the API docs for your version):
```python
transform = StainNormalizationHE(
    target="normalize",
    stain_estimation_method="macenko",
    optical_density_threshold=0.15,   # OD cutoff separating tissue from background
    angular_percentile=0.01,          # robust percentile for stain-vector estimation
    regularizer=0.01,                 # regularization for the vahadane method
    background_intensity=245,          # assumed background intensity
    stain_matrix_target_od=None,       # reference stain matrix in OD space (optional)
    max_c_target=None,                 # reference max stain concentrations (optional)
)
```

Note: `StainNormalizationHE` does not take a `tissue_mask_name` argument; it estimates stain vectors from the tile's optical-density distribution.

**Workflow:**
1. Convert RGB to optical density (OD)
2. Estimate stain matrix (H&E vectors)
3. Decompose into stain concentrations
4. Normalize to reference stain distribution
5. Reconstruct normalized RGB image

**Example: detect tissue, then normalize:**
```python
from pathml.preprocessing import Pipeline, TissueDetectionHE, StainNormalizationHE

pipeline = Pipeline([
    TissueDetectionHE(),  # creates the 'tissue' mask
    StainNormalizationHE(
        target="normalize",
        stain_estimation_method="macenko",
    ),
])
```

## Quality Control Transforms

### Artifact Detection

**LabelArtifactTileHE:**
```python
from pathml.preprocessing import LabelArtifactTileHE

# Label tiles containing artifacts
transform = LabelArtifactTileHE(
    pen_threshold=0.5,  # Threshold for pen marking detection
    bubble_threshold=0.5  # Threshold for bubble detection
)
```
- Detects pen markings, bubbles, and other artifacts
- Labels affected tiles for filtering

**LabelWhiteSpaceHE:**
```python
from pathml.preprocessing import LabelWhiteSpaceHE

# Label tiles with excessive white space
transform = LabelWhiteSpaceHE(
    threshold=0.9,  # Fraction of white pixels
    mask_name='white_space'
)
```
- Identifies tiles with mostly background
- Useful for filtering uninformative tiles

## Multiparametric Imaging Transforms

### Cell Segmentation

**SegmentMIF:**
```python
from pathml.preprocessing import SegmentMIF

# Segment cells using the DeepCell Mesmer model.
# nuclear_channel / cytoplasm_channel are INTEGER indices into the
# (collapsed) channel stack, not marker name strings.
transform = SegmentMIF(
    model="mesmer",
    nuclear_channel=0,        # e.g. index of DAPI
    cytoplasm_channel=29,     # e.g. index of a membrane/cytoplasm marker
    image_resolution=0.377,   # microns per pixel
)
```
- Uses the DeepCell Mesmer model for whole-cell + nuclear segmentation
- Channels are specified by integer index into the channel stack
- Writes `cell_segmentation` and `nuclear_segmentation` masks
- Fine-grained control is via `preprocess_kwargs`, `postprocess_kwargs_nuclear`, and `postprocess_kwargs_whole_cell` (there is no `compartment` argument)

**SegmentMIFRemote:**
```python
from pathml.preprocessing import SegmentMIFRemote

# Remote inference using the DeepCell Kiosk API (no local GPU)
transform = SegmentMIFRemote(
    model="mesmer",
    nuclear_channel=0,
    cytoplasm_channel=29,
)
```
- Same functionality as SegmentMIF but offloads inference to the DeepCell service
- No local GPU required; suitable for batch processing

### Marker Quantification

**QuantifyMIF:**
```python
from pathml.preprocessing import QuantifyMIF

# Quantify marker expression per cell. The only argument is the name of
# the segmentation mask produced by SegmentMIF.
transform = QuantifyMIF(segmentation_mask="cell_segmentation")
```
- Extracts mean marker intensity per segmented cell across all channels
- Computes per-cell spatial/morphology info
- Writes an AnnData object to `slide.counts` for downstream single-cell analysis (read it back via `slide.counts`, or concatenate across a dataset with `anndata.concat([s.counts for s in dataset.slides], ...)`)

### CODEX/Vectra Specific

**CollapseRunsCODEX:**
```python
from pathml.preprocessing import CollapseRunsCODEX

# Consolidate a multi-cycle CODEX z-stack by selecting a focal plane.
transform = CollapseRunsCODEX(z=0)  # 'z' selects the z-plane index
```
- Collapses the CODEX z-stack/cycles into a single multi-channel image
- `z` selects the focal plane (0-indexed)

**CollapseRunsVectra:**
```python
from pathml.preprocessing import CollapseRunsVectra

# Process Vectra multiplex IF data
transform = CollapseRunsVectra(
    wavelengths=[520, 570, 620, 670, 780]  # Emission wavelengths
)
```

## Building Comprehensive Pipelines

### Basic H&E Preprocessing Pipeline

```python
from pathml.preprocessing import (
    Pipeline,
    TissueDetectionHE,
    StainNormalizationHE,
    NucleusDetectionHE,
    MedianBlur,
    LabelWhiteSpaceHE
)

pipeline = Pipeline([
    # 1. Quality control
    LabelWhiteSpaceHE(threshold=0.9),

    # 2. Noise reduction
    MedianBlur(kernel_size=3),

    # 3. Tissue detection
    TissueDetectionHE(min_region_size=500),

    # 4. Stain normalization
    StainNormalizationHE(
        target='normalize',
        stain_estimation_method='macenko',
        tissue_mask_name='tissue'
    ),

    # 5. Nucleus detection
    NucleusDetectionHE(threshold=0.3)
])
```

### CODEX Multiparametric Pipeline

```python
from pathml.preprocessing import (
    Pipeline,
    CollapseRunsCODEX,
    SegmentMIF,
    QuantifyMIF
)

codex_pipeline = Pipeline([
    # 1. Collapse the z-stack
    CollapseRunsCODEX(z=0),

    # 2. Cell segmentation (integer channel indices)
    SegmentMIF(
        model="mesmer",
        nuclear_channel=0,
        cytoplasm_channel=29,
        image_resolution=0.377,
    ),

    # 3. Quantify markers -> writes AnnData to slide.counts
    QuantifyMIF(segmentation_mask="cell_segmentation"),
])
```

### Advanced Pipeline with Quality Control

```python
from pathml.preprocessing import (
    Pipeline,
    LabelWhiteSpaceHE,
    LabelArtifactTileHE,
    TissueDetectionHE,
    MorphOpen,
    MorphClose,
    StainNormalizationHE,
    AdaptiveHistogramEqualization
)

advanced_pipeline = Pipeline([
    # Stage 1: Quality control
    LabelWhiteSpaceHE(threshold=0.85),
    LabelArtifactTileHE(pen_threshold=0.5, bubble_threshold=0.5),

    # Stage 2: Tissue detection
    TissueDetectionHE(threshold=10, min_region_size=1000),
    MorphOpen(kernel_size=5, mask_name='tissue'),
    MorphClose(kernel_size=7, mask_name='tissue'),

    # Stage 3: Stain normalization
    StainNormalizationHE(
        target='normalize',
        stain_estimation_method='vahadane',
        tissue_mask_name='tissue'
    ),

    # Stage 4: Contrast enhancement
    AdaptiveHistogramEqualization(clip_limit=0.03, tile_grid_size=(8, 8))
])
```

## Running Pipelines

### Single Slide Processing

```python
from pathml.core import HESlide

# Load slide
wsi = HESlide("slide.svs")

# Run the pipeline. SlideData.run handles tiling internally; pass tiling
# options through (e.g. tile_size / level) as supported by your version.
wsi.run(pipeline, tile_size=256, tile_stride=256, level=1)

# Access processed data
for tile in wsi.tiles:
    normalized_image = tile.image
    tissue_mask = tile.masks.get("tissue")
    nucleus_mask = tile.masks.get("nucleus")
```

### Batch Processing with Distributed Execution

```python
from pathml.core import HESlide, SlideDataset
from dask.distributed import Client
import glob

# Start Dask client
client = Client(n_workers=8, threads_per_worker=2, memory_limit="4GB")

# Create dataset from slide objects
slide_paths = glob.glob("data/*.svs")
dataset = SlideDataset([HESlide(p) for p in slide_paths])

# Run pipeline in parallel (tiling options passed through to run())
dataset.run(pipeline, client=client, tile_size=512, tile_stride=512, level=1)

# Persist results to h5path (one file per slide in the directory)
dataset.write("processed/")

client.close()
```

### Filtering Tiles After Processing

Filter on masks produced by the pipeline (e.g. keep only tissue tiles):

```python
wsi.run(pipeline, tile_size=256, level=1)

tissue_tiles = [
    tile for tile in wsi.tiles
    if tile.masks.get("tissue") is not None and tile.masks["tissue"].any()
]
```

## Performance Optimization

### Memory Management

```python
from pathml.core import HESlide, SlideDataset

# Process large datasets in batches
batch_size = 100
for i in range(0, len(slide_paths), batch_size):
    batch_paths = slide_paths[i:i + batch_size]
    batch_dataset = SlideDataset([HESlide(p) for p in batch_paths])
    batch_dataset.run(pipeline, client=client)
    batch_dataset.write(f"processed/batch_{i}/")
```

### GPU Acceleration

Certain transforms leverage GPU acceleration when available:

```python
import torch

# Check GPU availability
print(f"CUDA available: {torch.cuda.is_available()}")

# Transforms that benefit from GPU:
# - SegmentMIF (Mesmer deep learning model)
# - StainNormalizationHE (matrix operations)
```

### Parallel Workers Configuration

```python
from dask.distributed import Client

# CPU-bound tasks (image processing)
client = Client(
    n_workers=8,
    threads_per_worker=1,  # Use processes, not threads
    memory_limit='8GB'
)

# GPU tasks (deep learning inference)
client = Client(
    n_workers=2,  # Fewer workers for GPU
    threads_per_worker=4,
    processes=True
)
```

## Custom Transforms

Create custom preprocessing operations by subclassing `Transform`:

```python
from pathml.preprocessing.transforms import Transform
import numpy as np

class CustomTransform(Transform):
    def __init__(self, param1, param2):
        self.param1 = param1
        self.param2 = param2

    def apply(self, tile):
        # Access tile image
        image = tile.image

        # Apply custom operation
        processed = self.custom_operation(image, self.param1, self.param2)

        # Update tile
        tile.image = processed

        return tile

    def custom_operation(self, image, param1, param2):
        # Implement custom logic
        return processed_image

# Use in pipeline
pipeline = Pipeline([
    CustomTransform(param1=10, param2=0.5),
    # ... other transforms
])
```

## Best Practices

1. **Order transforms appropriately:**
   - Quality control first (LabelWhiteSpace, LabelArtifact)
   - Noise reduction early (Blur)
   - Tissue detection before stain normalization
   - Stain normalization before color-dependent operations

2. **Use tissue masks for stain normalization:**
   - Improves accuracy by excluding background
   - `TissueDetectionHE()` then `StainNormalizationHE(tissue_mask_name='tissue')`

3. **Apply morphological operations to clean masks:**
   - `MorphOpen` to remove small false positives
   - `MorphClose` to fill small gaps

4. **Leverage distributed processing for large datasets:**
   - Use Dask for parallel execution
   - Configure workers based on available resources

5. **Save intermediate results:**
   - Store processed data to HDF5 for reuse
   - Avoid reprocessing computationally expensive transforms

6. **Validate preprocessing on sample images:**
   - Visualize intermediate steps
   - Tune parameters on representative samples before batch processing

7. **Handle edge cases:**
   - Check for empty masks before downstream operations
   - Validate tile quality before expensive computations

## Common Issues and Solutions

**Issue: Stain normalization produces artifacts**
- Use tissue mask to exclude background
- Try different stain estimation method (macenko vs. vahadane)
- Verify optical density parameters match your images

**Issue: Out of memory during pipeline execution**
- Reduce number of Dask workers
- Decrease tile size
- Process images at lower pyramid level
- Enable memory_limit parameter in Dask client

**Issue: Tissue detection misses tissue regions**
- Adjust threshold parameter
- Use saturation channel: `use_saturation=True`
- Reduce min_region_size to capture smaller tissue fragments

**Issue: Nucleus detection is inaccurate**
- Verify stain separation quality (visualize hematoxylin channel)
- Adjust threshold parameter
- Apply stain normalization before nucleus detection

## Additional Resources

- **PathML Preprocessing API:** https://pathml.readthedocs.io/en/latest/api_preprocessing_reference.html
- **Stain Normalization Methods:**
  - Macenko et al. 2009: "A method for normalizing histology slides for quantitative analysis"
  - Vahadane et al. 2016: "Structure-Preserving Color Normalization and Sparse Stain Separation"
- **DeepCell Mesmer:** https://www.deepcell.org/ (cell segmentation model)