"""ROI warp validation checks.
NOTE: Runs inside a Nipype Function node; the thesis logger is not available at
runtime. ``warnings.warn()`` is used for soft failures and ``PipelineError``
is raised for hard failures.
"""
from typing import List, Optional
[docs]
def validate_warped_rois_task(
roi_paths: List[str],
reference_image: str,
original_roi_paths: Optional[List[str]] = None,
min_voxels: int = 10,
singularity_threshold: float = 1e-6,
volume_ratio_min: float = 0.5,
volume_ratio_max: float = 1.5,
) -> List[str]:
"""Validate warped ROI masks for voxel count, centroid, and volume consistency.
Performs four checks per non-empty ROI path:
* **(a) Voxel count** — raises ``PipelineError`` if the mask contains fewer
than ``min_voxels`` non-zero voxels (hard failure, halts pipeline).
* **(a2) Affine singularity** — raises ``PipelineError`` if the absolute
determinant of the affine's 3×3 rotation/scale block is below
``singularity_threshold``. A near-zero determinant causes
``probtrackx2`` to crash with ``inv(): matrix is singular`` (hard
failure, halts pipeline).
* **(b) Centroid plausibility** — emits a ``UserWarning`` if the mask
centroid falls outside the subject brain mask (soft failure, continues).
* **(c) Volume consistency** — emits a ``UserWarning`` if the warped mask
volume ratio (warped / original) falls outside
``[volume_ratio_min, volume_ratio_max]`` (soft failure, continues).
Only performed when ``original_roi_paths`` is provided.
Empty-string entries in ``roi_paths`` are silently skipped (absent masks).
Args:
roi_paths: Paths to warped ROI masks in subject space. Empty strings
are treated as absent and skipped.
reference_image: Path to the subject-space brain mask used for centroid
plausibility check (check b).
original_roi_paths: Corresponding pre-warp ROI paths used for volume
comparison (check c). Must have the same length as ``roi_paths``
when provided. ``None`` or empty-string entries skip the check for
that mask.
min_voxels: Minimum number of non-zero voxels required in each warped
mask. Fewer voxels triggers the hard failure (check a).
singularity_threshold: Minimum absolute determinant of the affine
3×3 rotation/scale block. Below this the mask is rejected as
degenerate (check a2).
volume_ratio_min: Lower bound for warped-vs-original voxel count ratio.
Below this a warning is emitted (check c).
volume_ratio_max: Upper bound for warped-vs-original voxel count ratio.
Above this a warning is emitted (check c).
Returns:
``roi_paths`` unchanged (passthrough for Nipype graph chaining).
Raises:
PipelineError: If any warped mask is empty (fewer than ``min_voxels``
non-zero voxels) or has a degenerate (near-singular) affine matrix.
"""
import warnings
import nibabel as nib
import numpy as np
from thesis.core.exceptions import PipelineError
# Load brain mask once (shared across all ROI checks)
brain_mask_data = None
if reference_image:
try:
brain_mask_data = nib.load(reference_image).get_fdata() # type: ignore[attr-defined]
except Exception: # pragma: no cover — missing mask is handled downstream
pass
for idx, roi_path in enumerate(roi_paths):
if not roi_path:
continue
# --- Load warped mask ---
img = nib.load(roi_path)
data = img.get_fdata() # type: ignore[attr-defined]
# (a) Voxel count check — hard failure
n_voxels = int(np.sum(data > 0))
if n_voxels < min_voxels:
raise PipelineError(
f"Warped ROI mask is empty (or near-empty): '{roi_path}' has "
f"{n_voxels} non-zero voxel(s), minimum required is {min_voxels}. "
"Check that the warp field is correctly specified and points to the "
"right subject."
)
# (a2) Affine singularity check — hard failure
# probtrackx2 calls inv() on the affine internally; a near-zero determinant
# causes 'inv(): matrix is singular' and an immediate crash.
affine_det = abs(float(np.linalg.det(img.affine[:3, :3]))) # type: ignore[attr-defined]
if affine_det < singularity_threshold:
raise PipelineError(
f"Warped ROI mask has a degenerate (near-singular) affine matrix: "
f"'{roi_path}' (det={affine_det:.2e}). "
"probtrackx2 will crash with 'inv(): matrix is singular'. "
"Check that the warp field produces a valid output geometry for "
"this subject."
)
# (b) Centroid plausibility — soft failure
if brain_mask_data is not None:
coords = np.argwhere(data > 0)
centroid = tuple(int(round(c)) for c in coords.mean(axis=0))
# Clamp to valid array indices
clamped = tuple(max(0, min(c, s - 1)) for c, s in zip(centroid, brain_mask_data.shape))
if brain_mask_data[clamped] == 0:
warnings.warn(
f"centroid of warped ROI '{roi_path}' at voxel {centroid} is "
"outside the subject brain mask. The warp may be off-target.",
UserWarning,
stacklevel=2,
)
# (c) Volume ratio — soft failure (only when original provided)
if original_roi_paths is not None and idx < len(original_roi_paths):
orig_path = original_roi_paths[idx]
if orig_path:
orig_data = nib.load(orig_path).get_fdata() # type: ignore[attr-defined]
orig_voxels = int(np.sum(orig_data > 0))
if orig_voxels > 0:
ratio = n_voxels / orig_voxels
if ratio < volume_ratio_min or ratio > volume_ratio_max:
warnings.warn(
f"volume ratio of warped ROI '{roi_path}' vs original "
f"'{orig_path}' is {ratio:.2f} "
f"(expected {volume_ratio_min}–{volume_ratio_max}). "
"The warp may have significantly changed the ROI size.",
UserWarning,
stacklevel=2,
)
return roi_paths
[docs]
def validate_warped_rois_passthrough_task(
seed: str,
waypoints_file: str,
stop_mask: str,
avoid_mask: str,
target_mask: str,
reference_image: str,
min_voxels: int = 10,
singularity_threshold: float = 1e-6,
volume_ratio_min: float = 0.5,
volume_ratio_max: float = 1.5,
) -> tuple:
"""Validate all warped ROI masks and pass them through unchanged.
Nipype-compatible wrapper around :func:`validate_warped_rois_task`. Collects
individual NIfTI mask paths and expands the waypoints text file into its
constituent paths, then delegates validation. Returns all five inputs unchanged
so this node can be inserted transparently between the transformer and any
downstream node.
Runs inside a Nipype Function node — all imports are local; ``warnings.warn``
is used for soft failures (no module-level logger).
Args:
seed: Path to the transformed seed mask, or empty string if absent.
waypoints_file: Path to text file listing transformed waypoint NIfTI paths,
or empty string if absent.
stop_mask: Path to transformed stop mask, or empty string.
avoid_mask: Path to transformed avoid mask, or empty string.
target_mask: Path to transformed target mask, or empty string.
reference_image: Subject-space brain mask for centroid plausibility check.
min_voxels: Minimum non-zero voxel count per warped mask.
singularity_threshold: Minimum absolute determinant of affine 3x3 block.
volume_ratio_min: Lower bound for warped-vs-original voxel ratio.
volume_ratio_max: Upper bound for warped-vs-original voxel ratio.
Returns:
5-tuple ``(seed, waypoints_file, stop_mask, avoid_mask, target_mask)``
unchanged.
Raises:
PipelineError: If any transformed mask has fewer than ``min_voxels``
non-zero voxels or has a degenerate affine matrix.
"""
from pathlib import Path
from thesis.workflows.hcp.operations.validation import validate_warped_rois_task
roi_paths: list = []
for p in (seed, stop_mask, avoid_mask, target_mask):
if p:
roi_paths.append(p)
if waypoints_file and Path(waypoints_file).exists():
with open(waypoints_file, "r", encoding="utf-8") as fh:
for line in fh:
stripped = line.strip()
if stripped:
roi_paths.append(stripped)
validate_warped_rois_task(
roi_paths,
reference_image,
min_voxels=min_voxels,
singularity_threshold=singularity_threshold,
volume_ratio_min=volume_ratio_min,
volume_ratio_max=volume_ratio_max,
)
return (seed, waypoints_file, stop_mask, avoid_mask, target_mask)
[docs]
def verify_waypoint_avoid_overlap_task(
seed: str,
waypoints_file: str,
stop_mask: str,
avoid_mask: str,
target_mask: str,
max_overlap_fraction: float = 0.9,
) -> tuple:
"""Fail-fast guard against waypoint masks landing inside the avoid mask.
ProbTrackX2 must satisfy two mutually exclusive constraints when a waypoint
overlaps the avoid mask: streamlines must pass through the waypoint AND
must not enter avoid. If a waypoint is wholly contained in avoid, every
streamline is killed at step 0 and the run silently produces
``waytotal=0``.
This check runs after DWI-grid resampling — the last point before the
masks are handed to ``probtrackx2`` — and raises ``PipelineError`` when
any waypoint has more than ``max_overlap_fraction`` of its voxels inside
avoid. Pure passthrough on success.
Args:
seed: Path to seed mask (unused, for signature compatibility).
waypoints_file: Text file listing waypoint mask paths (one per line),
or empty string if absent.
stop_mask: Path to stop mask (unused, for signature compatibility).
avoid_mask: Path to avoid mask, or empty string. When absent, the
check is a no-op.
target_mask: Path to target mask (unused, for signature compatibility).
max_overlap_fraction: Maximum allowed fraction of waypoint voxels
inside the avoid mask before failing. Default 0.9 (90%). The
threshold is deliberately high — anatomically valid waypoints
often clip the avoid mask near boundaries (e.g. the cerebral
peduncle near the brainstem can show ~30-40% overlap on healthy
subjects). The bug this guard targets is the catastrophic case
(e.g. an IC waypoint warped into the wrong hemisphere where the
entire mask lands inside the contralateral-hemisphere avoid),
which produces ~100% overlap.
Returns:
5-tuple of the five inputs unchanged.
Raises:
PipelineError: If any waypoint has more than ``max_overlap_fraction``
of its voxels inside the avoid mask.
"""
from pathlib import Path
import nibabel as nib
import numpy as np
from thesis.core.exceptions import PipelineError
if not avoid_mask or not Path(avoid_mask).exists():
return (seed, waypoints_file, stop_mask, avoid_mask, target_mask)
if not waypoints_file or not Path(waypoints_file).exists():
return (seed, waypoints_file, stop_mask, avoid_mask, target_mask)
avoid_data = nib.load(avoid_mask).get_fdata() > 0 # type: ignore[attr-defined]
with open(waypoints_file, "r", encoding="utf-8") as handle:
waypoint_paths = [line.strip() for line in handle if line.strip()]
blocked: list[str] = []
for wp_path in waypoint_paths:
if not wp_path or not Path(wp_path).exists():
continue
wp_data = nib.load(wp_path).get_fdata() > 0 # type: ignore[attr-defined]
wp_voxels = int(np.sum(wp_data))
if wp_voxels == 0:
continue
# Both masks share the DWI grid at this point; element-wise AND.
if wp_data.shape != avoid_data.shape:
continue
overlap = int(np.sum(wp_data & avoid_data))
fraction = overlap / wp_voxels
if fraction > max_overlap_fraction:
blocked.append(
f" - {wp_path}: {overlap}/{wp_voxels} voxels "
f"({fraction:.1%}) inside avoid mask"
)
if blocked:
raise PipelineError(
"Waypoint mask(s) overlap the avoid mask beyond the allowed "
f"{max_overlap_fraction:.0%} threshold. ProbTrackX2 cannot produce "
"any streamline that simultaneously passes through such a waypoint "
"and avoids the avoid mask — the run would silently emit "
"waytotal=0. Investigate the upstream warp / hemisphere routing "
"for this subject.\n" + "\n".join(blocked)
)
return (seed, waypoints_file, stop_mask, avoid_mask, target_mask)
