Cleaned up

ihc_quantification_simplified.py → ihc_quantification.py

@@ -190,7 +190,6 @@ def segment_spheroids(
     C = 15
     block_size = round(min(gray.shape) * 0.14)
     block_size = block_size - 1 if block_size % 2 == 0 else block_size
-    gray = cv2.cvtColor(preview, cv2.COLOR_RGB2GRAY)
     tissue_thresh = cv2.adaptiveThreshold(
         gray, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY_INV, block_size, C
     )

ihc_quantification_gpt.py

@@ -1,853 +0,0 @@
-"""Quantify DAB staining per spheroid while excluding debris.
-
-The source image is never inpainted.  Spheroids and debris are detected on an
-8-bit preview, while H-DAB stain separation is performed on the original image
-values.  Results from multiple runs are merged into one wide CSV: each
-spheroid is a column and each measurement is a row.
-"""
-
-from __future__ import annotations
-
-import argparse
-import csv
-import hashlib
-import os
-import re
-import tempfile
-from dataclasses import dataclass
-from pathlib import Path
-
-import cv2
-import numpy as np
-import tifffile
-from skimage.color import hdx_from_rgb, separate_stains
-from skimage.morphology import remove_small_holes
-from skimage.util import img_as_float32
-
-from cleaning import _as_rgb_uint8, create_debris_mask
-from ihc_qc import DEFAULT_QC_MAX_DIMENSION, save_qc_image
-
-
-IMAGE_EXTENSIONS = {".bmp", ".jpeg", ".jpg", ".png", ".tif", ".tiff"}
-DEFAULT_OUTPUT_DIR = Path("out/quantification")
-DEFAULT_ADAPTIVE_BLOCK_FRACTION = 0.14
-DEFAULT_ADAPTIVE_OFFSET = 15.0
-DEFAULT_OPENING_RADIUS = 2
-DEFAULT_CLOSING_RADIUS = 25
-DEFAULT_MIN_SPHEROID_AREA_FRACTION = 0.0025
-DEFAULT_MAX_HOLE_AREA_FRACTION = 0.00025
-CONCENTRATION_PATTERN = re.compile(
-    r"^(?P<value>\d+(?:[.,]\d+)?)\s*(?P<unit>pM|nM|uM|µM|μM|mM)$",
-    re.IGNORECASE,
-)
-MAGNIFICATION_PATTERN = re.compile(r"^\d+(?:[.,]\d+)?x$", re.IGNORECASE)
-CSV_METRICS = (
-    "image_name",
-    "source_path",
-    "cell_line",
-    "culture_model",
-    "treatment",
-    "concentration_value",
-    "concentration_unit",
-    "antibody",
-    "magnification",
-    "image_number",
-    "spheroid_id",
-    "touches_image_border",
-    "centroid_x_px",
-    "centroid_y_px",
-    "total_spheroid_area_px",
-    "segmented_tissue_area_px",
-    "valid_measured_area_px",
-    "debris_excluded_area_px",
-    "debris_excluded_fraction",
-    "mean_dab_intensity",
-    "median_dab_intensity",
-    "integrated_dab_intensity",
-    "positive_dab_fraction",
-    "dab_positive_threshold",
-    "stain_separation",
-    "background_correction",
-)
-
-
-@dataclass(frozen=True)
-class FilenameMetadata:
-    cell_line: str = ""
-    culture_model: str = ""
-    treatment: str = ""
-    concentration_value: float | None = None
-    concentration_unit: str = ""
-    antibody: str = ""
-    magnification: str = ""
-    image_number: int | None = None
-
-
-@dataclass(frozen=True)
-class SpheroidMeasurement:
-    spheroid_id: int
-    touches_border: bool
-    centroid_x: float
-    centroid_y: float
-    envelope_area: int
-    tissue_area: int
-    valid_area: int
-    debris_area: int
-    debris_fraction: float
-    mean_dab: float
-    median_dab: float
-    integrated_dab: float
-    positive_fraction: float | None
-
-
-@dataclass(frozen=True)
-class ImageResult:
-    image_path: Path
-    image_id: str
-    metadata: FilenameMetadata
-    labels: np.ndarray
-    debris_mask: np.ndarray
-    dab: np.ndarray
-    measurements: tuple[SpheroidMeasurement, ...]
-    dab_positive_threshold: float | None
-    background_corrected: bool
-
-
-def parse_filename(image_path: str | Path) -> FilenameMetadata:
-    """Extract experimental metadata from a conventionally named image."""
-    tokens = Path(image_path).stem.split()
-    magnification_index = next(
-        (
-            index
-            for index, token in enumerate(tokens)
-            if MAGNIFICATION_PATTERN.fullmatch(token)
-        ),
-        None,
-    )
-
-    if magnification_index is None or magnification_index < 3:
-        return FilenameMetadata()
-
-    number_index = magnification_index + 1
-    if number_index != len(tokens) - 1 or not tokens[number_index].isdigit():
-        return FilenameMetadata()
-
-    experiment_tokens = tokens[2 : magnification_index - 1]
-    concentration_value = None
-    concentration_unit = ""
-    treatment_tokens: list[str] = []
-
-    for token in experiment_tokens:
-        concentration_match = CONCENTRATION_PATTERN.fullmatch(token)
-        if concentration_match and concentration_value is None:
-            concentration_value = float(
-                concentration_match.group("value").replace(",", ".")
-            )
-            unit = concentration_match.group("unit").replace("µ", "u").replace(
-                "μ", "u"
-            )
-            concentration_unit = {
-                "pm": "pM",
-                "nm": "nM",
-                "um": "uM",
-                "mm": "mM",
-            }[unit.lower()]
-        else:
-            treatment_tokens.append(token)
-
-    return FilenameMetadata(
-        cell_line=tokens[0],
-        culture_model=tokens[1],
-        treatment=" ".join(treatment_tokens),
-        concentration_value=concentration_value,
-        concentration_unit=concentration_unit,
-        antibody=tokens[magnification_index - 1],
-        magnification=tokens[magnification_index],
-        image_number=int(tokens[number_index]),
-    )
-
-
-def _select_rgb_plane(image: np.ndarray) -> np.ndarray:
-    """Select a single RGB plane without changing the source bit depth."""
-    image = np.asarray(image)
-
-    while image.ndim > 2 and 1 in image.shape:
-        image = np.squeeze(image)
-
-    if image.ndim != 3:
-        raise ValueError("Quantification requires a 2D RGB image.")
-
-    if image.shape[-1] in {3, 4}:
-        image = image[..., :3]
-    elif image.shape[0] in {3, 4}:
-        image = np.moveaxis(image[:3], 0, -1)
-    else:
-        raise ValueError("Quantification requires an RGB image with 3 channels.")
-
-    return image
-
-
-def read_original_rgb(image_path: str | Path) -> np.ndarray:
-    """Read an RGB image while preserving its original numeric values."""
-    image_path = Path(image_path)
-
-    if image_path.suffix.lower() in {".tif", ".tiff"}:
-        return _select_rgb_plane(tifffile.imread(image_path))
-
-    image = cv2.imread(str(image_path), cv2.IMREAD_UNCHANGED)
-    if image is None:
-        raise FileNotFoundError(f"Could not read image: {image_path}")
-
-    image = _select_rgb_plane(image)
-    return cv2.cvtColor(image, cv2.COLOR_BGR2RGB)
-
-
-def _odd_block_size(shape: tuple[int, int], requested: int | None) -> int:
-    shortest_side = min(shape)
-    if shortest_side < 3:
-        raise ValueError("Image is too small for adaptive segmentation.")
-
-    if requested is None:
-        block_size = round(shortest_side * DEFAULT_ADAPTIVE_BLOCK_FRACTION)
-    else:
-        block_size = requested
-
-    block_size = max(3, min(block_size, shortest_side))
-    if block_size % 2 == 0:
-        block_size -= 1
-    return block_size
-
-
-def segment_spheroids(
-    preview_rgb: np.ndarray,
-    *,
-    adaptive_block_size: int | None = None,
-    adaptive_offset: float = DEFAULT_ADAPTIVE_OFFSET,
-    opening_radius: int = DEFAULT_OPENING_RADIUS,
-    closing_radius: int = DEFAULT_CLOSING_RADIUS,
-    min_area_fraction: float = DEFAULT_MIN_SPHEROID_AREA_FRACTION,
-    max_hole_area_fraction: float = DEFAULT_MAX_HOLE_AREA_FRACTION,
-) -> np.ndarray:
-    """Return row-major spheroid labels using local contrast segmentation."""
-    if not 0 < min_area_fraction < 1:
-        raise ValueError("min_area_fraction must be between 0 and 1.")
-    if opening_radius < 0 or closing_radius < 0:
-        raise ValueError("Morphology radii must be non-negative.")
-    if not 0 <= max_hole_area_fraction < 1:
-        raise ValueError("max_hole_area_fraction must be between 0 and 1.")
-
-    gray = cv2.cvtColor(preview_rgb, cv2.COLOR_RGB2GRAY)
-    block_size = _odd_block_size(gray.shape, adaptive_block_size)
-    tissue_seed = cv2.adaptiveThreshold(
-        gray,
-        255,
-        cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
-        cv2.THRESH_BINARY_INV,
-        block_size,
-        adaptive_offset,
-    )
-
-    # Remove isolated dark specks before closing so they cannot form bridges
-    # between a spheroid and unrelated debris or vignetted image edges.
-    if opening_radius:
-        kernel_size = 2 * opening_radius + 1
-        opening_kernel = cv2.getStructuringElement(
-            cv2.MORPH_ELLIPSE,
-            (kernel_size, kernel_size),
-        )
-        tissue_seed = cv2.morphologyEx(
-            tissue_seed,
-            cv2.MORPH_OPEN,
-            opening_kernel,
-        )
-
-    if closing_radius:
-        kernel_size = 2 * closing_radius + 1
-        closing_kernel = cv2.getStructuringElement(
-            cv2.MORPH_ELLIPSE,
-            (kernel_size, kernel_size),
-        )
-        tissue_seed = cv2.morphologyEx(
-            tissue_seed,
-            cv2.MORPH_CLOSE,
-            closing_kernel,
-        )
-
-    max_hole_area = round(
-        preview_rgb.shape[0]
-        * preview_rgb.shape[1]
-        * max_hole_area_fraction
-    )
-    if max_hole_area:
-        tissue_seed = remove_small_holes(
-            tissue_seed.astype(bool),
-            max_size=max_hole_area,
-            connectivity=2,
-        ).astype(np.uint8)
-
-    count, initial_labels, stats, centroids = cv2.connectedComponentsWithStats(
-        tissue_seed,
-        connectivity=8,
-    )
-    min_area = round(preview_rgb.shape[0] * preview_rgb.shape[1] * min_area_fraction)
-    components = [
-        component
-        for component in range(1, count)
-        if stats[component, cv2.CC_STAT_AREA] >= min_area
-    ]
-    components.sort(
-        key=lambda component: (
-            centroids[component, 1],
-            centroids[component, 0],
-        )
-    )
-
-    labels = np.zeros(initial_labels.shape, dtype=np.int32)
-    for spheroid_id, component in enumerate(components, start=1):
-        labels[initial_labels == component] = spheroid_id
-
-    return labels
-
-
-def _quadratic_design(x: np.ndarray, y: np.ndarray) -> np.ndarray:
-    return np.column_stack((np.ones_like(x), x, y, x * x, x * y, y * y))
-
-
-def spheroid_envelope(labels: np.ndarray) -> np.ndarray:
-    """Fill each spheroid's external contour for size and background fitting."""
-    envelope = np.zeros(labels.shape, dtype=np.uint8)
-    for spheroid_id in range(1, int(labels.max()) + 1):
-        component = (labels == spheroid_id).astype(np.uint8)
-        contours, _ = cv2.findContours(
-            component,
-            cv2.RETR_EXTERNAL,
-            cv2.CHAIN_APPROX_SIMPLE,
-        )
-        cv2.drawContours(envelope, contours, -1, spheroid_id, cv2.FILLED)
-    return envelope
-
-
-def _fit_background_coefficients(
-    rgb: np.ndarray,
-    background_mask: np.ndarray,
-    *,
-    target_samples: int = 100_000,
-) -> np.ndarray | None:
-    """Fit a robust quadratic illumination surface to empty background."""
-    height, width = background_mask.shape
-    stride = max(1, int(np.sqrt(background_mask.size / target_samples)))
-    sampled_mask = background_mask[::stride, ::stride]
-    sampled_y, sampled_x = np.nonzero(sampled_mask)
-    if sampled_x.size < 1_000:
-        return None
-
-    sampled_y = sampled_y * stride
-    sampled_x = sampled_x * stride
-    x = sampled_x.astype(np.float64) / max(width - 1, 1) * 2 - 1
-    y = sampled_y.astype(np.float64) / max(height - 1, 1) * 2 - 1
-    design = _quadratic_design(x, y)
-
-    coefficients = np.zeros((3, design.shape[1]), dtype=np.float64)
-    for channel in range(3):
-        values = rgb[sampled_y, sampled_x, channel].astype(np.float64)
-        keep = np.ones(values.shape, dtype=bool)
-
-        for _ in range(4):
-            coefficients[channel], *_ = np.linalg.lstsq(
-                design[keep],
-                values[keep],
-                rcond=None,
-            )
-            residuals = values - design @ coefficients[channel]
-            center = np.median(residuals[keep])
-            mad = np.median(np.abs(residuals[keep] - center))
-            if mad <= np.finfo(np.float64).eps:
-                break
-            robust_sigma = 1.4826 * mad
-            new_keep = (
-                (residuals >= center - 2.5 * robust_sigma)
-                & (residuals <= center + 3.5 * robust_sigma)
-            )
-            if new_keep.sum() < 1_000 or np.array_equal(new_keep, keep):
-                break
-            keep = new_keep
-
-    return coefficients
-
-
-def correct_background_illumination(
-    rgb: np.ndarray,
-    tissue_mask: np.ndarray,
-    debris_mask: np.ndarray,
-) -> tuple[np.ndarray, bool]:
-    """Normalize smooth per-channel illumination using empty background."""
-    background_mask = ~tissue_mask & ~debris_mask
-    coefficients = _fit_background_coefficients(rgb, background_mask)
-    if coefficients is None:
-        return rgb, False
-
-    height, width = tissue_mask.shape
-    x = np.linspace(-1, 1, width, dtype=np.float32)[None, :]
-    y = np.linspace(-1, 1, height, dtype=np.float32)[:, None]
-    corrected = rgb.copy()
-
-    for channel in range(3):
-        c0, cx, cy, cxx, cxy, cyy = coefficients[channel]
-        surface = (
-            c0
-            + cx * x
-            + cy * y
-            + cxx * x * x
-            + cxy * x * y
-            + cyy * y * y
-        ).astype(np.float32)
-        surface = np.clip(surface, 0.05, 1.5)
-        corrected[..., channel] = np.clip(
-            corrected[..., channel] / surface,
-            0,
-            1,
-        )
-
-    return corrected, True
-
-
-def extract_dab(rgb: np.ndarray) -> np.ndarray:
-    """Extract the DAB component with scikit-image's fixed H-DAB matrix."""
-    stains = separate_stains(rgb, hdx_from_rgb)
-    dab = stains[..., 1].astype(np.float32)
-    return dab
-
-
-def _touches_border(component_mask: np.ndarray) -> bool:
-    return bool(
-        component_mask[0].any()
-        or component_mask[-1].any()
-        or component_mask[:, 0].any()
-        or component_mask[:, -1].any()
-    )
-
-
-def measure_spheroids(
-    labels: np.ndarray,
-    debris_mask: np.ndarray,
-    dab: np.ndarray,
-    dab_positive_threshold: float | None,
-) -> tuple[SpheroidMeasurement, ...]:
-    measurements: list[SpheroidMeasurement] = []
-    debris = debris_mask.astype(bool)
-    envelope_labels = spheroid_envelope(labels)
-
-    for spheroid_id in range(1, int(labels.max()) + 1):
-        spheroid = labels == spheroid_id
-        envelope = envelope_labels == spheroid_id
-        valid = spheroid & ~debris
-        tissue_area = int(spheroid.sum())
-        envelope_area = int(envelope.sum())
-        valid_area = int(valid.sum())
-        debris_area = tissue_area - valid_area
-        if not valid_area:
-            continue
-
-        y, x = np.nonzero(spheroid)
-        dab_values = dab[valid]
-        positive_fraction = None
-        if dab_positive_threshold is not None:
-            positive_fraction = float(
-                np.count_nonzero(dab_values >= dab_positive_threshold) / valid_area
-            )
-
-        measurements.append(
-            SpheroidMeasurement(
-                spheroid_id=spheroid_id,
-                touches_border=_touches_border(spheroid),
-                centroid_x=float(x.mean()),
-                centroid_y=float(y.mean()),
-                envelope_area=envelope_area,
-                tissue_area=tissue_area,
-                valid_area=valid_area,
-                debris_area=debris_area,
-                debris_fraction=debris_area / tissue_area,
-                mean_dab=float(dab_values.mean()),
-                median_dab=float(np.median(dab_values)),
-                integrated_dab=float(dab_values.sum(dtype=np.float64)),
-                positive_fraction=positive_fraction,
-            )
-        )
-
-    return tuple(measurements)
-
-
-def _safe_identifier(value: str) -> str:
-    identifier = re.sub(r"[^A-Za-z0-9._ -]+", "_", value).strip(" ._")
-    return identifier or "image"
-
-
-def quantify_image(
-    image_path: str | Path,
-    *,
-    dab_positive_threshold: float | None = None,
-    adaptive_block_size: int | None = None,
-    adaptive_offset: float = DEFAULT_ADAPTIVE_OFFSET,
-    opening_radius: int = DEFAULT_OPENING_RADIUS,
-    closing_radius: int = DEFAULT_CLOSING_RADIUS,
-    min_area_fraction: float = DEFAULT_MIN_SPHEROID_AREA_FRACTION,
-    max_hole_area_fraction: float = DEFAULT_MAX_HOLE_AREA_FRACTION,
-    background_correction: bool = True,
-) -> tuple[ImageResult, np.ndarray]:
-    image_path = Path(image_path)
-    original = read_original_rgb(image_path)
-    preview = _as_rgb_uint8(original)
-    labels = segment_spheroids(
-        preview,
-        adaptive_block_size=adaptive_block_size,
-        adaptive_offset=adaptive_offset,
-        opening_radius=opening_radius,
-        closing_radius=closing_radius,
-        min_area_fraction=min_area_fraction,
-        max_hole_area_fraction=max_hole_area_fraction,
-    )
-    if not labels.max():
-        raise ValueError(
-            "No spheroids were detected. Review the QC segmentation settings."
-        )
-
-    debris_mask = create_debris_mask(preview).astype(bool)
-    rgb = img_as_float32(original)
-    was_corrected = False
-    if background_correction:
-        rgb, was_corrected = correct_background_illumination(
-            rgb,
-            spheroid_envelope(labels) > 0,
-            debris_mask,
-        )
-
-    dab = extract_dab(rgb)
-    measurements = measure_spheroids(
-        labels,
-        debris_mask,
-        dab,
-        dab_positive_threshold,
-    )
-    if not measurements:
-        raise ValueError("No spheroid contained valid pixels after debris exclusion.")
-
-    result = ImageResult(
-        image_path=image_path.resolve(),
-        image_id=_safe_identifier(image_path.name),
-        metadata=parse_filename(image_path),
-        labels=labels,
-        debris_mask=debris_mask,
-        dab=dab,
-        measurements=measurements,
-        dab_positive_threshold=dab_positive_threshold,
-        background_corrected=was_corrected,
-    )
-    return result, preview
-
-
-def _format_float(value: float | None) -> str:
-    if value is None:
-        return ""
-    return f"{value:.8g}"
-
-
-def _measurement_values(
-    result: ImageResult,
-    measurement: SpheroidMeasurement,
-) -> dict[str, str]:
-    metadata = result.metadata
-    return {
-        "image_name": result.image_path.name,
-        "source_path": str(result.image_path),
-        "cell_line": metadata.cell_line,
-        "culture_model": metadata.culture_model,
-        "treatment": metadata.treatment,
-        "concentration_value": _format_float(metadata.concentration_value),
-        "concentration_unit": metadata.concentration_unit,
-        "antibody": metadata.antibody,
-        "magnification": metadata.magnification,
-        "image_number": (
-            str(metadata.image_number) if metadata.image_number is not None else ""
-        ),
-        "spheroid_id": str(measurement.spheroid_id),
-        "touches_image_border": str(measurement.touches_border).lower(),
-        "centroid_x_px": _format_float(measurement.centroid_x),
-        "centroid_y_px": _format_float(measurement.centroid_y),
-        "total_spheroid_area_px": str(measurement.envelope_area),
-        "segmented_tissue_area_px": str(measurement.tissue_area),
-        "valid_measured_area_px": str(measurement.valid_area),
-        "debris_excluded_area_px": str(measurement.debris_area),
-        "debris_excluded_fraction": _format_float(measurement.debris_fraction),
-        "mean_dab_intensity": _format_float(measurement.mean_dab),
-        "median_dab_intensity": _format_float(measurement.median_dab),
-        "integrated_dab_intensity": _format_float(measurement.integrated_dab),
-        "positive_dab_fraction": _format_float(measurement.positive_fraction),
-        "dab_positive_threshold": _format_float(result.dab_positive_threshold),
-        "stain_separation": "scikit-image H-DAB (hdx_from_rgb)",
-        "background_correction": (
-            "quadratic empty-background fit"
-            if result.background_corrected
-            else "none"
-        ),
-    }
-
-
-def _read_existing_results(
-    csv_path: Path,
-) -> tuple[list[str], dict[str, dict[str, str]], list[str]]:
-    if not csv_path.exists():
-        return [], {}, []
-
-    with csv_path.open(newline="", encoding="utf-8") as handle:
-        reader = csv.reader(handle)
-        rows = list(reader)
-
-    if not rows:
-        return [], {}, []
-    if not rows[0] or rows[0][0] != "metric":
-        raise ValueError(
-            f"Existing results file does not use the expected wide format: {csv_path}"
-        )
-
-    columns = rows[0][1:]
-    values: dict[str, dict[str, str]] = {}
-    metric_order: list[str] = []
-    for row in rows[1:]:
-        if not row:
-            continue
-        metric = row[0]
-        metric_order.append(metric)
-        padded = row[1:] + [""] * max(0, len(columns) - len(row[1:]))
-        values[metric] = dict(zip(columns, padded[: len(columns)], strict=True))
-
-    return columns, values, metric_order
-
-
-def update_results_csv(csv_path: str | Path, results: list[ImageResult]) -> None:
-    """Replace reprocessed images and append new spheroid columns atomically."""
-    csv_path = Path(csv_path)
-    csv_path.parent.mkdir(parents=True, exist_ok=True)
-    columns, values, previous_metric_order = _read_existing_results(csv_path)
-    source_paths = {str(result.image_path) for result in results}
-    source_row = values.get("source_path", {})
-    removed_columns = {
-        column for column in columns if source_row.get(column) in source_paths
-    }
-    columns = [column for column in columns if column not in removed_columns]
-    for metric_values in values.values():
-        for column in removed_columns:
-            metric_values.pop(column, None)
-
-    for result in results:
-        base_prefix = result.image_id
-        proposed_columns = [
-            f"{base_prefix}::spheroid_{measurement.spheroid_id:03d}"
-            for measurement in result.measurements
-        ]
-        if any(column in columns for column in proposed_columns):
-            digest = hashlib.sha1(
-                str(result.image_path).encode("utf-8")
-            ).hexdigest()[:8]
-            base_prefix = f"{base_prefix} [{digest}]"
-
-        for measurement in result.measurements:
-            column = f"{base_prefix}::spheroid_{measurement.spheroid_id:03d}"
-            columns.append(column)
-            for metric, value in _measurement_values(result, measurement).items():
-                values.setdefault(metric, {})[column] = value
-
-    metric_order = list(CSV_METRICS)
-    metric_order.extend(
-        metric
-        for metric in previous_metric_order
-        if metric not in metric_order
-    )
-    metric_order.extend(
-        metric for metric in values if metric not in metric_order
-    )
-
-    with tempfile.NamedTemporaryFile(
-        mode="w",
-        newline="",
-        encoding="utf-8",
-        dir=csv_path.parent,
-        prefix=f".{csv_path.name}.",
-        suffix=".tmp",
-        delete=False,
-    ) as handle:
-        temporary_path = Path(handle.name)
-        writer = csv.writer(handle)
-        writer.writerow(["metric", *columns])
-        for metric in metric_order:
-            metric_values = values.get(metric, {})
-            writer.writerow(
-                [metric, *(metric_values.get(column, "") for column in columns)]
-            )
-
-    os.replace(temporary_path, csv_path)
-
-
-def discover_images(inputs: list[str]) -> list[Path]:
-    images: list[Path] = []
-    seen: set[Path] = set()
-
-    for input_value in inputs:
-        input_path = Path(input_value)
-        if input_path.is_dir():
-            candidates = sorted(
-                path
-                for path in input_path.rglob("*")
-                if path.is_file() and path.suffix.lower() in IMAGE_EXTENSIONS
-            )
-        elif input_path.is_file():
-            if input_path.suffix.lower() not in IMAGE_EXTENSIONS:
-                raise ValueError(f"Unsupported image extension: {input_path}")
-            candidates = [input_path]
-        else:
-            raise FileNotFoundError(f"Input does not exist: {input_path}")
-
-        for candidate in candidates:
-            resolved = candidate.resolve()
-            if resolved not in seen:
-                images.append(candidate)
-                seen.add(resolved)
-
-    if not images:
-        raise ValueError("No supported images were found.")
-    return images
-
-
-def parse_args() -> argparse.Namespace:
-    parser = argparse.ArgumentParser(
-        description=(
-            "Quantify H-DAB signal per spheroid and update a wide results CSV."
-        )
-    )
-    parser.add_argument(
-        "images",
-        nargs="+",
-        help="Image files or directories. Directories are searched recursively.",
-    )
-    parser.add_argument(
-        "--output-dir",
-        type=Path,
-        default=DEFAULT_OUTPUT_DIR,
-        help=f"Results directory (default: {DEFAULT_OUTPUT_DIR}).",
-    )
-    parser.add_argument(
-        "--dab-positive-threshold",
-        type=float,
-        default=None,
-        help=(
-            "Fixed H-DAB threshold derived from negative controls. If omitted, "
-            "positive fractions are left blank."
-        ),
-    )
-    parser.add_argument(
-        "--adaptive-block-size",
-        type=int,
-        default=None,
-        help="Odd local-threshold window (default: 14%% of shortest image side).",
-    )
-    parser.add_argument(
-        "--adaptive-offset",
-        type=float,
-        default=DEFAULT_ADAPTIVE_OFFSET,
-        help=f"Local-threshold offset (default: {DEFAULT_ADAPTIVE_OFFSET:g}).",
-    )
-    parser.add_argument(
-        "--opening-radius",
-        type=int,
-        default=DEFAULT_OPENING_RADIUS,
-        help=f"Speck-removal radius in pixels (default: {DEFAULT_OPENING_RADIUS}).",
-    )
-    parser.add_argument(
-        "--closing-radius",
-        type=int,
-        default=DEFAULT_CLOSING_RADIUS,
-        help=(
-            "Spheroid-mask closing radius in pixels "
-            f"(default: {DEFAULT_CLOSING_RADIUS})."
-        ),
-    )
-    parser.add_argument(
-        "--min-spheroid-area-fraction",
-        type=float,
-        default=DEFAULT_MIN_SPHEROID_AREA_FRACTION,
-        help=(
-            "Minimum spheroid area as a fraction of image area "
-            f"(default: {DEFAULT_MIN_SPHEROID_AREA_FRACTION:g})."
-        ),
-    )
-    parser.add_argument(
-        "--max-hole-area-fraction",
-        type=float,
-        default=DEFAULT_MAX_HOLE_AREA_FRACTION,
-        help=(
-            "Fill internal mask holes up to this fraction of image area "
-            f"(default: {DEFAULT_MAX_HOLE_AREA_FRACTION:g})."
-        ),
-    )
-    parser.add_argument(
-        "--no-background-correction",
-        action="store_true",
-        help="Disable quadratic empty-background illumination correction.",
-    )
-    parser.add_argument(
-        "--qc-max-dimension",
-        type=int,
-        default=DEFAULT_QC_MAX_DIMENSION,
-        help=(
-            "Maximum width or height of each QC panel "
-            f"(default: {DEFAULT_QC_MAX_DIMENSION})."
-        ),
-    )
-    return parser.parse_args()
-
-
-def main() -> None:
-    args = parse_args()
-    image_paths = discover_images(args.images)
-    results: list[ImageResult] = []
-    qc_dir = args.output_dir / "qc"
-
-    if args.dab_positive_threshold is None:
-        print(
-            "Note: positive_dab_fraction will be blank until a fixed "
-            "--dab-positive-threshold is supplied."
-        )
-
-    for image_path in image_paths:
-        result, preview = quantify_image(
-            image_path,
-            dab_positive_threshold=args.dab_positive_threshold,
-            adaptive_block_size=args.adaptive_block_size,
-            adaptive_offset=args.adaptive_offset,
-            opening_radius=args.opening_radius,
-            closing_radius=args.closing_radius,
-            min_area_fraction=args.min_spheroid_area_fraction,
-            max_hole_area_fraction=args.max_hole_area_fraction,
-            background_correction=not args.no_background_correction,
-        )
-        qc_path = qc_dir / f"{_safe_identifier(result.image_path.stem)}_qc.png"
-        save_qc_image(
-            qc_path,
-            preview,
-            result.labels,
-            result.debris_mask,
-            result.dab,
-            max_dimension=args.qc_max_dimension,
-        )
-        results.append(result)
-        print(
-            f"{image_path}: {len(result.measurements)} spheroids; "
-            f"QC: {qc_path}"
-        )
-
-    csv_path = args.output_dir / "ihc_quantification.csv"
-    update_results_csv(csv_path, results)
-    print(f"Updated results: {csv_path}")
-
-
-if __name__ == "__main__":
-    main()