import openslide                                          # Library purpose-built for multi-resolution medical image formats (SVS, TIFF, DICOM)

slide = openslide.OpenSlide("patient_123.svs")            # Opens the file handle without loading pixel data — the pyramid index is read lazily

print(f"Levels: {slide.level_count}")                    # Number of resolution levels in the pyramid (typically 5–6 for a 40x scan)
print(f"Level 0 dimensions: {slide.level_dimensions[0]}")  # Full-res dimensions at 40x — e.g. (100000, 80000); just a tuple, no memory allocated

thumbnail = slide.get_thumbnail((1024, 1024))             # Reads a downsampled preview by selecting the appropriate pyramid level internally — useful for quick UI previews

region = slide.read_region((20000, 30000), 0, (1024, 1024))  # Reads a 1024x1024 pixel region from level 0 at (x=20000, y=30000) — only this patch is decoded into RAM

import openslide                                          # WSI reader that navigates the multi-resolution pyramid without loading the full slide
import numpy as np                                        # Used for pixel-level thresholding and coordinate arithmetic on the tissue mask
from PIL import Image                                     # Converts OpenSlide's RGBA region into a standard RGB PIL image for downstream use

def extract_tissue_tiles(svs_path, tile_size=256, threshold=0.1):
                                                          # tile_size=256: standard patch size for WSI models (balances context vs batch size)
                                                          # threshold=0.1: at least 10% of the patch must contain tissue to keep it
    slide = openslide.OpenSlide(svs_path)                 # Opens the slide handle — only metadata is loaded, not pixel data

    level = slide.level_count - 1                         # Uses the lowest-resolution pyramid level for the tissue mask (fastest to read)
    thumbnail = slide.read_region((0, 0), level, slide.level_dimensions[level])
                                                          # Reads the full downsampled thumbnail — e.g. 6250x5000 instead of 100000x80000
    thumbnail = thumbnail.convert("RGB")                  # OpenSlide returns RGBA; strip alpha for NumPy processing
    thumbnail_np = np.array(thumbnail)                    # Convert PIL image to a (H, W, 3) array for thresholding

    gray = np.mean(thumbnail_np, axis=2)                  # Average across RGB channels to get a single-channel intensity map
    tissue_mask = gray < 220                              # Binary mask: pixels darker than 220 (on 0–255 scale) are considered tissue, not background glass

    scale = slide.level_dimensions[0][0] / slide.level_dimensions[level][0]
                                                          # Downscale factor from level 0 to the thumbnail level (e.g. 16x for 5 levels down)

    tiles = []                                            # Accumulator for valid (tissue-containing) PIL Image patches
    step = int(tile_size * scale)                         # Step size at the thumbnail level corresponding to one tile at level 0

    for y in range(0, tissue_mask.shape[0] * int(scale) - tile_size, tile_size):
                                                          # Iterate y in level-0 pixel coordinates; shape[0] * scale recovers the full-res height
        for x in range(0, tissue_mask.shape[1] * int(scale) - tile_size, tile_size):
                                                          # Iterate x in level-0 coordinates; ensures tiles cover the entire slide area
            tx, ty = int(x / scale), int(y / scale)       # Map level-0 coordinates back to thumbnail coordinates for the mask lookup
            if tx >= tissue_mask.shape[1] or ty >= tissue_mask.shape[0]:
                continue                                  # Skip tiles whose mapped position falls outside the thumbnail bounds

            patch = tissue_mask[ty:ty + 1, tx:tx + 1]     # Sample a single-pixel row from the mask at the tile's location
            if np.mean(patch) > threshold:                # If the fraction of tissue pixels in this mask sample exceeds 10%, read the full tile
                tile = slide.read_region((x, y), 0, (tile_size, tile_size))
                                                          # Reads only the needed 256x256 region from level 0 — no surrounding pixels decoded
                tiles.append(tile)                        # Append the PIL Image tile to the results list

    return tiles                                          # Returns only tissue-bearing patches; empty-background tiles are discarded

import numpy as np                                        # All stain normalization is linear algebra on optical density vectors

def stain_normalize(source, target, alpha=1, beta=0.15):
                                                          # source: tile from a hospital with unknown stain appearance
                                                          # target: reference tile whose stain appearance we want to match
                                                          # alpha, beta: percentile trimming params (unused in this simplified version)
    source_od = -np.log(np.maximum(source / 255.0, 1e-6))
                                                          # Convert source RGB [0-255] to optical density: OD = -log(I / I_0)
                                                          # OD space makes stain separation linear because Beer-Lambert law is additive
    target_od = -np.log(np.maximum(target / 255.0, 1e-6))
                                                          # Same conversion for the target reference slide

    source_reshape = source_od.reshape((-1, 3))           # Flatten spatial dims (H*W, 3) so each row is a pixel's OD vector
    target_reshape = target_od.reshape((-1, 3))           # Same reshape for target to compute stain vectors per channel

    source_mask = np.mean(source_reshape, axis=1) > 0.05  # Keep only pixels whose mean OD > 0.05 — filters out glass background (low OD)
    target_mask = np.mean(target_reshape, axis=1) > 0.05  # Same foreground mask for the target; stain vectors should come from tissue only

    source_foreground = source_reshape[source_mask]       # Foreground-only pixel OD matrix for SVD — glass pixels would distort the stain basis
    target_foreground = target_reshape[target_mask]       # Target foreground; both matrices are (N_fg, 3) where N_fg is foreground pixel count

    _, _, source_V = np.linalg.svd(source_foreground, full_matrices=False)
                                                          # SVD on source foreground: V (3x3) rows are the principal stain directions
                                                          # First two rows correspond to Hematoxylin (blue/purple) and Eosin (pink/red)
    _, _, target_V = np.linalg.svd(target_foreground, full_matrices=False)
                                                          # SVD on target foreground: target's own stain basis to map into

    source_stains = source_od @ source_V.T                # Project every pixel's OD vector into the source stain basis (H&E space)
    target_stains = target_od @ target_V.T                # Project every pixel's OD vector into the target stain basis

    for i in range(3):                                    # Match the mean and std of each stain channel (H, E, residual) to the target
        source_std = np.std(source_stains[:, i])           # Standard deviation of the i-th stain component in the source
        target_std = np.std(target_stains[:, i])           # Target standard deviation — the spread we want to match
        source_mean = np.mean(source_stains[:, i])         # Mean of the i-th stain component in the source
        target_mean = np.mean(target_stains[:, i])         # Target mean — the intensity centroid we shift toward
        source_stains[:, i] = (source_stains[:, i] - source_mean) * \
            (target_std / source_std) + target_mean
                                                          # Z-score the source stain, then scale by target std and shift by target mean
                                                          # This aligns the statistical distribution of H&E intensities between labs

    normalized_od = source_stains @ target_V              # Project the adjusted stain coefficients back into RGB OD space via target basis
    normalized = np.exp(-normalized_od).reshape(source.shape)
                                                          # Invert Beer-Lambert: OD -> transmission -> RGB; reshape back to original (H, W, 3)

    return np.clip(normalized * 255, 0, 255).astype(np.uint8)
                                                          # Scale OD back to [0-255], clip to valid range, and return uint8 for model ingestion

class WSIDataset(Dataset):                                # Extends PyTorch Dataset — yields tiles on demand without holding pixel data in memory
    def __init__(self, slide_paths, tile_size=256, transform=None):
        self.slide_paths = slide_paths                    # List of SVS/TIFF file paths; each can be 1–2 GB on disk
        self.tile_size = tile_size                        # Standard 256x256 or 512x512 at 20x magnification
        self.transform = transform                        # Augmentation pipeline (flip, rotate, color jitter, stain normalize)
        self.tiles = []                                   # Flat list of (slide_path, x, y) tuples — pixel data is never stored here
        self._precompute_tiles()                          # Runs tissue detection on all slides upfront; slow but avoids repeated mask computation

    def _precompute_tiles(self):                          # Scans every slide once to build the coordinate index
        for path in self.slide_paths:                     # Iterate over all SVS files in the cohort (hundreds to thousands)
            slide = openslide.OpenSlide(path)             # Open handle — reads pyramid metadata but no pixel data yet
            for x, y in self._tissue_coords(slide):       # Generator yielding (x, y) top-left corners of tissue-bearing regions
                self.tiles.append((path, x, y))           # Store only the reference and coordinate — ~24 bytes per tile vs ~1 MB for the pixels

    def __getitem__(self, idx):                           # Called by the DataLoader for each batch item
        path, x, y = self.tiles[idx]                      # Look up the precomputed coordinate for this sample index
        slide = openslide.OpenSlide(path)                 # Open the slide handle again (or fetch from LRU cache in production)
        tile = slide.read_region((x, y), 0, (self.tile_size, self.tile_size))
                                                          # Lazy decode only the 256x256 region at level 0 — no adjacent pixels loaded
        tile = tile.convert("RGB")                        # OpenSlide returns RGBA; strip alpha channel for 3-channel model input

        if self.transform:                                # Apply stain normalization + augmentation if provided
            tile = self.transform(tile)                   # Usually a torchvision Compose: ToTensor -> MacenkoNormalize -> RandomFlip

        return tile                                       # Returns a single PIL Image (or Tensor if transform converts it) for the batch collator

def detect_artifacts(tile):                               # Filters out tiles with blur (out-of-focus scan) or excessive edges (tissue fold/ink)
    laplacian = cv2.Laplacian(np.array(tile), cv2.CV_64F) # Laplacian edge detection — high-frequency content indicates sharp focus
    focus_score = laplacian.var()                          # Variance of the Laplacian: low variance means the tile is blurry (no fine detail)

    gray = cv2.cvtColor(np.array(tile), cv2.COLOR_RGB2GRAY) # Single-channel intensity for edge detection with Canny
    edges = cv2.Canny(gray, 50, 150)                       # Canny with low=50, high=150 — tuned for H&E stain edge contrast
    edge_density = np.mean(edges > 0)                      # Fraction of pixels that are edge pixels in the binary Canny output

    if focus_score < 100 or edge_density > 0.3:            # Blur threshold: < 100 means out-of-focus (common in WSI scanning artifacts)
                                                           # Edge-threshold: > 0.3 means tissue fold, ink mark, or bubble (too many sharp boundaries)
        return False                                        # Reject this tile — it would confuse the model with non-tissue features
    return True                                             # Accept the tile as valid tissue for training or inference

from functools import lru_cache                            # lru_cache imported but not used in this snippet — could decorate get_slide for automatic eviction

class SlideCache:                                          # Manual LRU-like cache for OpenSlide handles to cap OS file descriptor usage
    def __init__(self, max_slides=10):                     # max_slides=10: typical physical memory limit — each handle reserves significant RAM for the pyramid index
        self.max_slides = max_slides                       # Upper bound on concurrently open slide handles
        self.cache = {}                                    # Dict mapping file path -> OpenSlide object; insertion order tracks recency

    def get_slide(self, path):                             # Request a slide by path; returns cached handle if available
        if path not in self.cache:                         # Cache miss — need to open a new handle
            if len(self.cache) >= self.max_slides:         # At capacity: evict the oldest (first inserted) slide handle
                oldest = next(iter(self.cache))            # dict insertion order is preserved in Python 3.7+; first key is least recently opened
                self.cache[oldest].close()                 # Close the OpenSlide file handle to release OS file descriptors and allocated memory
                del self.cache[oldest]                     # Remove the evicted entry from the dict
            self.cache[path] = openslide.OpenSlide(path)   # Open a new handle — reads the pyramid index but doesn't decode pixels
        return self.cache[path]                            # Return the (possibly just-opened) slide handle for region reads