diff --git a/config.yaml b/config.yaml index ccd6f0e..8ed7274 100644 --- a/config.yaml +++ b/config.yaml @@ -89,4 +89,10 @@ organ_adjacency_map: # affine reference, default as liver affine_reference_file_name: liver.nii.gz -if_save_combined_label: True \ No newline at end of file +if_save_combined_label: True + +# optional read-only vertebral instance audit +vertebrae_instance_analysis: + enabled: false + output_dir_name: vertebrae_analysis + use_reference_as_ct: false diff --git a/docs/config.md b/docs/config.md index 7a6b0fb..256c7e2 100644 --- a/docs/config.md +++ b/docs/config.md @@ -65,4 +65,56 @@ ├── liver.nii.gz ... └── veins.nii.gz - ``` \ No newline at end of file + ``` + +7. `vertebrae_instance_analysis`: optional read-only vertebral-instance audit. + The block is disabled by default: + + ```yaml + vertebrae_instance_analysis: + enabled: false + output_dir_name: vertebrae_analysis + use_reference_as_ct: false + ``` + + When enabled, ShapeKit analyzes the anatomical-name vertebral masks before + any organ or vertebral postprocessing and writes one canonical JSON report + per case: + + ``` + OUTPUT + ├── case_001 + │ └── segmentations + └── vertebrae_analysis + └── case_001.json + ``` + + The report is diagnostic only. It never changes, deletes, fills, merges, or + relabels segmentation voxels, and it is not used by the existing + postprocessing pipeline. + + `use_reference_as_ct` must remain `false` unless + `affine_reference_file_name` identifies a voxel-aligned Computed Tomography + image with finite intensity values. ShapeKit does not infer CT provenance + from a filename or intensity heuristic. If CT use is explicitly requested + but the loaded reference fails the shape or finite-value checks, the audit + logs a warning and falls back to geometry-only evidence. + + The audit directory name must be one relative path component. Configuration + mistakes stop batch startup with a clear error. An existing symbolic link + is never accepted as the audit directory, and the resolved audit directory + must remain below the resolved output root. + + Each invocation owns only its uniquely named temporary file. A successful + report is atomically replaced and represents the most recent successfully + completed write. A failed invocation is logged with a traceback and does + not remove or truncate an existing successful report, including one written + concurrently by another invocation. Segmentation processing continues + unchanged after an audit failure. + + Audit JSON files currently retain the owner-only permissions created by + `mkstemp`. ShapeKit does not impose a different report permission policy. + + Audit files are not backfilled for cases excluded by + `--continue_prediction`; rerun without resume filtering to audit previously + completed cases. diff --git a/docs/vertebrae_instance_analysis.md b/docs/vertebrae_instance_analysis.md new file mode 100644 index 0000000..7b362a3 --- /dev/null +++ b/docs/vertebrae_instance_analysis.md @@ -0,0 +1,305 @@ +# Read-only vertebral instance analysis + +`utils.vertebrae_instance_analysis` provides a diagnostic API for examining +physical vertebral instances, protected thick-core candidates, separation +between EDT-derived core extents, and anatomical-name sequence consistency. + +The analyzer is intentionally separate from ShapeKit's vertebra +post-processing pipeline. It does not delete, fill, merge, relabel, or save +segmentation voxels, and importing it does not change ShapeKit runtime +behavior. + +## Public API + +```python +from utils.vertebrae_instance_analysis import ( + VertebralInstanceAnalysisConfig, + analyze_vertebral_instances, +) + +report = analyze_vertebral_instances( + segmentation_dict, + affine=affine, + ordered_anatomical_names=( + "vertebrae_L5", + "vertebrae_L4", + "vertebrae_L3", + "vertebrae_L2", + "vertebrae_L1", + "vertebrae_T12", + "vertebrae_T11", + "vertebrae_T10", + "vertebrae_T9", + "vertebrae_T8", + "vertebrae_T7", + "vertebrae_T6", + "vertebrae_T5", + "vertebrae_T4", + "vertebrae_T3", + "vertebrae_T2", + "vertebrae_T1", + "vertebrae_C7", + "vertebrae_C6", + "vertebrae_C5", + "vertebrae_C4", + "vertebrae_C3", + "vertebrae_C2", + "vertebrae_C1", + ), + ct=None, + config=VertebralInstanceAnalysisConfig(), +) +``` + +`ordered_anatomical_names` is always interpreted +**inferior-to-superior**. The API uses anatomical names rather than combined +label values, so it does not depend on a particular numeric class map. + +## Input contract + +### `segmentation_dict` + +- A mapping from anatomical names to three-dimensional arrays. +- Nonzero values are treated as foreground. +- All requested masks must have the same shape. +- Missing and empty masks are accepted. +- Keys not listed in `ordered_anatomical_names` are ignored. +- Input arrays are never modified. +- Overlapping requested masks are reported as ambiguous. Their overlap is not + resolved using key order or a numeric-label priority. + +### `affine` + +- A finite, invertible 4×4 voxel-to-world affine. +- Physical spacing is derived from the affine. +- Axis permutations, axis sign changes, and in-plane rotations that preserve + superior-inferior alignment are supported. +- Mild superior-inferior obliquity up to 20 degrees is supported. Material + shear or greater SI obliquity produces an `unsupported_affine` unresolved + report instead of silently using an inappropriate voxel axis. +- The final affine row must be approximately `[0, 0, 0, 1]`; a malformed + homogeneous row raises `ValueError`. +- The 20-degree limit is intentionally conservative because discrete EDT and + persistence measurements can change near decision thresholds on more + oblique grids. Arbitrary orthogonal rotations are not claimed as supported. + +### `ct` + +CT means Computed Tomography input and is optional. When supplied, it must +already be aligned voxel-for-voxel with the segmentation masks, have the same +shape, and contain only finite intensity values. Shape mismatch or non-finite +values raise `ValueError`. + +CT values contribute only bone-support confidence. CT never changes: + +- the vertebral union; +- thick-core boundaries; +- physical-instance boundaries; +- thick-core candidate boundaries or core-separation measurements. + +When CT is absent, the analyzer reports `ct_evidence: "unavailable"`, sets +per-instance confidence mode to `geometry_only`, and uses a stricter +geometry-only confidence threshold. When CT is present, low bone support is +reported explicitly as `low_ct_bone_support`; it can keep a candidate +unresolved, but it still cannot change candidate boundaries. + +### `config` + +All anatomical distances and volumes are configurable in physical units. +Dimensionless confidence and shape thresholds are also explicit in +`VertebralInstanceAnalysisConfig`. + +All numeric configuration values must be finite. In particular, +`bone_hu_threshold` accepts finite positive, zero, or negative CT intensity +thresholds, while non-finite values raise `ValueError`. + +The defaults are provisional, conservative diagnostic starting values. They +are not population-validated constants or claims of clinical validity or +population-level generalization. The effective configuration is included in +every report. + +## Analysis stages + +1. Validate masks, anatomical-name order, CT shape, and affine geometry. +2. Construct the unmodified union of requested vertebral masks. +3. Estimate a smooth physical spine trajectory from the largest component of + each available named mask, using a polynomial least-squares fit followed by + at most one distance-thresholded refit. This is a deterministic provisional + estimator, not a fully robust trajectory estimator. +4. Run a spacing-aware distance transform and identify thick interior + candidates. +5. Measure volume, persistence, compactness, trajectory distance, name + composition, and optional CT bone support. +6. Analyze a label-independent thick-area profile for candidate peaks and + valleys between EDT-derived thick cores. +7. Sort selected candidates by physical inferior-to-superior position. +8. Report duplicate, internal-missing, nonmonotonic, abnormal-spacing, overlap, + and ambiguous-identity findings. + +No stage produces a proposed correction. + +## Protected and unresolved cores + +A thick-core candidate is reported as `protected_high_confidence` only when +multiple +independent signals agree: + +- sufficient thick-core physical volume; +- sufficient but not implausibly long superior-inferior persistence; +- compact transverse geometry; +- proximity to the estimated physical spine trajectory; +- a decisive anatomical-name vote; +- adequate CT bone support when CT is supplied; +- adequate separation evidence from neighboring thick cores. + +Otherwise it remains unresolved. Unresolved statuses include: + +- `unresolved_low_confidence`; +- `unresolved_mixed_identity`; +- `unresolved_overlap`; +- `unresolved_boundary_truncated`. + +The words “protected” and “core” are diagnostic. They mean that a thick, +compact interior candidate has high-confidence identity evidence. They do not +guarantee that the region is a clinical vertebral-body core, and this module +still does not edit any mask. + +## Rejected thick-core candidates + +Every EDT-derived thick component is accounted for. A component farther than +`trajectory_tube_radius_mm` from the provisional trajectory is excluded from +accepted vertebral instances but retained in `rejected_candidates` with: + +- a deterministic candidate ID; +- centroid in world millimetres; +- thick-core voxel count and physical volume; +- trajectory distance; +- status `rejected_unresolved`; +- reason `core_outside_trajectory_tube`. + +At least one rejected candidate also produces an unresolved +`off_trajectory_core` anomaly, so the report cannot claim an unqualified +`continuous_sequence`. The anomaly links to the rejected records through +`affected_candidate_ids`; it does not represent them as accepted instances. + +## Core-separation terminology + +`inferior_core_separation_mm` and `superior_core_separation_mm` measure the +physical separation between extents of neighboring EDT-derived thick cores. +They are threshold-dependent diagnostic measurements. They are not anatomical +intervertebral foreground gaps or disc-space measurements. + +## Partial field of view and missing identities + +Missing endpoint names are not automatically treated as missing anatomy. + +The analyzer reports an internal missing identity only when two consecutive, +high-confidence physical cores skip one or more anatomical names between +them. It does not infer missing identities: + +- above the most superior confident core; +- below the most inferior confident core; +- through an unresolved physical core. + +Foreground touching the physical inferior or superior array face is reported +as possible boundary truncation. Absence of boundary contact does not prove +that the anatomical field of view is complete; endpoint omissions can +therefore produce `extent_uncertain`. + +Duplicate or missing patterns that could reflect transitional anatomy remain +diagnostic findings. The analyzer does not force a standard identity onto +transitional anatomy. + +## Output contract + +The return value contains JSON-compatible Python primitives only: + +- dictionaries with string keys; +- lists; +- strings; +- integers; +- finite floats; +- booleans; +- `None`. + +It contains no NumPy arrays, masks, voxel-index lists, NIfTI objects, or +corrected labels. + +For canonical deterministic serialization: + +```python +import json + +canonical_json = json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + allow_nan=False, +) +``` + +Instances are ordered inferior-to-superior and receive stable report-local IDs +such as `instance_001`. Anomalies and name-composition entries also use stable +sorting. + +## Relationship to automatic vertebra re-identification + +This module is an independent audit layer. It neither imports nor wraps an +automatic vertebra re-identification implementation. + +It conceptually overlaps with upstream PR #5 in constructing a joint vertebral +mask, applying a physical distance transform, and ordering thick cores. PR #5 +uses those signals for automatic offset voting, watershed rebuilding, cleanup, +and runtime relabeling. This module instead accepts configurable anatomical +names, never relabels voxels, preserves uncertainty, and returns deterministic +structured diagnostic evidence. It does not duplicate PR #5's automatic +correction path. + +Potential future callers may use the report to decide whether additional +human review or a separately validated correction method is appropriate. Such +integration is outside this module's contract. + +In particular, this analyzer does not: + +- run a fixed global label offset; +- perform nearest-core watershed relabeling; +- suppress secondary connected components; +- fill holes; +- rebuild individual or combined NIfTI outputs; +- approve an automatic correction. + +## Tests + +The synthetic and API-contract suite requires no medical data: + +```bash +PYTHONDONTWRITEBYTECODE=1 \ +python -B -m unittest discover -s tests \ + -p 'test_vertebrae_instance_analysis.py' -v +``` + +The suite covers thick cores with thin appendages, posterior-like negative +controls, separated and merged instances, rejected off-trajectory candidates, +sequence anomalies, complete and partial overlaps, partial field of view, +optional and invalid CT, anisotropic spacing, supported mild obliquity, +unsupported geometry, fragmented masks, input non-mutation, deterministic +serialization, numeric-label independence, and multiprocessing consistency. + +## Limitations + +- Geometry alone cannot establish clinical anatomical identity in every case. +- CT evidence is limited to local intensity support; it is not a learned shape + model. +- Strongly oblique or sheared grids are reported as unsupported rather than + resampled internally. +- Distance, occupancy, and connected-component arrays scale with the cropped + foreground volume. Component coordinates are extracted from deterministic + `ndimage.find_objects` slices rather than rescanning the complete crop once + per component. Clinical runtime and multiprocessing capacity still require + workload-specific benchmarking. +- Transitional anatomy, merged vertebrae, severe leakage, fractures, implants, + and incomplete scans may remain unresolved. +- The analyzer has not been demonstrated to generalize to unseen clinical + populations. +- Diagnostic sequence anomalies are not equivalent to known segmentation + errors or ground truth. diff --git a/main.py b/main.py index a09b67a..7aa87f4 100644 --- a/main.py +++ b/main.py @@ -3,6 +3,11 @@ from multiprocessing import cpu_count from utils.organs_postprocessing import * from utils.vertebrae_postprocessing import postprocessing_vertebrae +from utils.vertebrae_instance_audit import ( + ordered_vertebral_anatomical_names, + parse_vertebrae_instance_audit_config, + run_vertebrae_instance_audit, +) import logging import yaml import traceback @@ -32,6 +37,14 @@ reference_file_name = affine_reference_file_name # affine info data_type = np.int16 save_combined_label_bool = bool(config['if_save_combined_label']) +vertebrae_audit_config = parse_vertebrae_instance_audit_config( + config.get('vertebrae_instance_analysis') +) +ordered_vertebral_names = ( + ordered_vertebral_anatomical_names(class_map.values()) + if vertebrae_audit_config.enabled + else () +) ############################################################## @@ -224,9 +237,18 @@ def main(input_path, input_folder_name, output_path=None): target_axcodes=nib.aff2axcodes(img.affine) ) - segmentation = combine_segmentation_dict(segmentation_dict, class_map) patient_id = os.path.basename(input_path) + run_vertebrae_instance_audit( + segmentation_dict, + reference_img=img, + ordered_anatomical_names=ordered_vertebral_names, + output_root=output_path, + patient_id=patient_id, + config=vertebrae_audit_config, + logger=logging, + ) + segmentation = combine_segmentation_dict(segmentation_dict, class_map) postprocessed_segmentation_dict = process_organs( segmentation_dict, img, diff --git a/tests/test_vertebrae_instance_analysis.py b/tests/test_vertebrae_instance_analysis.py new file mode 100644 index 0000000..d75ce2f --- /dev/null +++ b/tests/test_vertebrae_instance_analysis.py @@ -0,0 +1,1079 @@ +"""Synthetic and API-contract tests for vertebrae_instance_analysis.""" + +from __future__ import annotations + +import hashlib +import itertools +import json +import math +import multiprocessing +import unittest +from pathlib import Path +from typing import Dict, Mapping, Sequence, Tuple + +import numpy as np + +from utils.vertebrae_instance_analysis import ( + VertebralInstanceAnalysisConfig, + analyze_vertebral_instances, +) + + +NAMES = ("vertebrae_L3", "vertebrae_L2", "vertebrae_L1") + + +def _canonical(report: Mapping[str, object]) -> str: + return json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + allow_nan=False, + ) + + +def _ellipsoid( + shape: Sequence[int], + center: Sequence[float], + radii: Sequence[float] = (11.0, 12.0, 8.0), +) -> np.ndarray: + grid = np.indices(shape, dtype=np.float64) + normalized = np.zeros(shape, dtype=np.float64) + for axis in range(3): + normalized += ((grid[axis] - center[axis]) / radii[axis]) ** 2 + return normalized <= 1.0 + + +def _physical_ellipsoid( + shape: Sequence[int], + affine: np.ndarray, + center_world: Sequence[float], + radii_mm: Sequence[float] = (11.0, 12.0, 8.0), +) -> np.ndarray: + grid = np.indices(shape, dtype=np.float64).reshape(3, -1).T + world = grid @ affine[:3, :3].T + affine[:3, 3] + normalized = np.zeros(len(world), dtype=np.float64) + for axis in range(3): + normalized += ( + (world[:, axis] - center_world[axis]) / radii_mm[axis] + ) ** 2 + return (normalized <= 1.0).reshape(tuple(shape)) + + +def _body_with_thin_posterior( + shape: Sequence[int], center: Sequence[float] +) -> np.ndarray: + body = _ellipsoid(shape, center) + x, y, z = np.indices(shape) + process = ( + (np.abs(x - center[0]) <= 2) + & (y >= int(round(center[1] + 9))) + & (y <= int(round(center[1] + 20))) + & (np.abs(z - center[2]) <= 2) + ) + return body | process + + +def _standard_masks( + *, + shape: Sequence[int] = (48, 48, 112), + centers_z: Sequence[float] = (24.0, 52.0, 80.0), + names: Sequence[str] = NAMES, + posterior: bool = False, +) -> Dict[str, np.ndarray]: + result: Dict[str, np.ndarray] = {} + for name, z_value in zip(names, centers_z): + center = (24.0, 22.0, z_value) + result[name] = ( + _body_with_thin_posterior(shape, center) + if posterior + else _ellipsoid(shape, center) + ) + return result + + +def _fragmented_mask( + *, + shape: Sequence[int] = (112, 112, 112), + positions: Sequence[float] = (12.0, 38.0, 64.0, 90.0), +) -> np.ndarray: + mask = np.zeros(tuple(shape), dtype=bool) + for center in itertools.product(positions, repeat=3): + mask |= _ellipsoid(shape, center, radii=(9.0, 9.0, 9.0)) + return mask + + +def _rotation_about_x_affine( + angle_degrees: float, + pivot_world: Sequence[float], +) -> np.ndarray: + angle = math.radians(angle_degrees) + rotation = np.asarray( + [ + [1.0, 0.0, 0.0], + [0.0, math.cos(angle), -math.sin(angle)], + [0.0, math.sin(angle), math.cos(angle)], + ], + dtype=np.float64, + ) + affine = np.eye(4, dtype=np.float64) + affine[:3, :3] = rotation + pivot = np.asarray(pivot_world, dtype=np.float64) + affine[:3, 3] = pivot - rotation @ pivot + return affine + + +def _codes(report: Mapping[str, object]) -> Sequence[str]: + return [item["anomaly_code"] for item in report["anomalies"]] # type: ignore[index] + + +def _semantic_signature(report: Mapping[str, object]) -> Tuple[object, ...]: + instances = report["instances"] # type: ignore[index] + return ( + report["observed_sequence_inferior_to_superior"], + [item["status"] for item in instances], + tuple(_codes(report)), + report["overall_status"], + ) + + +def _array_digest(array: np.ndarray) -> str: + hasher = hashlib.sha256() + hasher.update(str(array.shape).encode("ascii")) + hasher.update(array.dtype.str.encode("ascii")) + hasher.update(np.ascontiguousarray(array).tobytes()) + return hasher.hexdigest() + + +def _reorient( + array: np.ndarray, + affine: np.ndarray, + permutation: Sequence[int], + flips: Sequence[bool], +) -> Tuple[np.ndarray, np.ndarray]: + transformed = np.transpose(array, axes=permutation) + for axis, should_flip in enumerate(flips): + if should_flip: + transformed = np.flip(transformed, axis=axis) + mapping = np.eye(4, dtype=np.float64) + mapping[:3, :3] = 0.0 + mapping[:3, 3] = 0.0 + for new_axis, old_axis in enumerate(permutation): + if flips[new_axis]: + mapping[old_axis, new_axis] = -1.0 + mapping[old_axis, 3] = transformed.shape[new_axis] - 1 + else: + mapping[old_axis, new_axis] = 1.0 + return transformed, affine @ mapping + + +def _multiprocessing_worker(payload: Tuple[Dict[str, np.ndarray], np.ndarray]) -> str: + masks, affine = payload + return _canonical( + analyze_vertebral_instances( + masks, + affine=affine, + ordered_anatomical_names=NAMES, + ) + ) + + +class VertebralInstanceAnalysisTests(unittest.TestCase): + def setUp(self) -> None: + self.affine = np.eye(4, dtype=np.float64) + + def analyze( + self, + masks: Mapping[str, np.ndarray], + *, + names: Sequence[str] = NAMES, + affine: np.ndarray = None, + ct: np.ndarray = None, + config: VertebralInstanceAnalysisConfig = None, + ) -> Mapping[str, object]: + return analyze_vertebral_instances( + masks, + affine=self.affine if affine is None else affine, + ordered_anatomical_names=names, + ct=ct, + config=config, + ) + + def test_thick_core_survives_thin_posterior_appendage(self) -> None: + masks = _standard_masks(centers_z=(52.0,), names=(NAMES[1],), posterior=True) + report = self.analyze(masks) + self.assertEqual(len(report["instances"]), 1) + instance = report["instances"][0] + self.assertEqual(instance["protected_core_identity"], NAMES[1]) + self.assertEqual(instance["status"], "protected_high_confidence") + + def test_two_separated_vertebrae_and_core_separation(self) -> None: + masks = _standard_masks(centers_z=(32.0, 68.0), names=NAMES[:2]) + report = self.analyze(masks, names=NAMES[:2]) + self.assertEqual( + report["observed_sequence_inferior_to_superior"], list(NAMES[:2]) + ) + self.assertGreater( + report["instances"][0]["superior_core_separation_mm"], + 0.0, + ) + self.assertEqual(report["overall_status"], "continuous_sequence") + + def test_core_separation_is_not_foreground_gap(self) -> None: + masks = _standard_masks( + centers_z=(32.0, 68.0), + names=NAMES[:2], + ) + report = self.analyze(masks, names=NAMES[:2]) + lower_z = np.flatnonzero(np.any(masks[NAMES[0]], axis=(0, 1))) + upper_z = np.flatnonzero(np.any(masks[NAMES[1]], axis=(0, 1))) + foreground_center_gap_mm = float(upper_z.min() - lower_z.max()) + core_separation = report["instances"][0][ + "superior_core_separation_mm" + ] + self.assertIsNotNone(core_separation) + self.assertGreater(core_separation, foreground_center_gap_mm) + self.assertNotIn("superior_gap_mm", report["instances"][0]) + + def test_off_trajectory_thick_candidate_is_reported(self) -> None: + shape = (150, 64, 80) + central = _ellipsoid( + shape, + (30.0, 30.0, 40.0), + radii=(13.0, 13.0, 9.0), + ) + remote = _ellipsoid( + shape, + (105.0, 30.0, 40.0), + radii=(10.0, 10.0, 8.0), + ) + report = self.analyze( + {NAMES[0]: central | remote}, + names=NAMES[:1], + ) + self.assertEqual(len(report["instances"]), 1) + self.assertEqual(len(report["rejected_candidates"]), 1) + rejected = report["rejected_candidates"][0] + self.assertEqual( + rejected["reasons"], + ["core_outside_trajectory_tube"], + ) + self.assertGreater(rejected["core_voxel_count"], 0) + self.assertGreater(rejected["core_volume_mm3"], 0.0) + self.assertGreater( + rejected["trajectory_distance_mm"], + VertebralInstanceAnalysisConfig().trajectory_tube_radius_mm, + ) + self.assertEqual(len(rejected["centroid_world_mm"]), 3) + self.assertIn("off_trajectory_core", _codes(report)) + anomaly = next( + item + for item in report["anomalies"] + if item["anomaly_code"] == "off_trajectory_core" + ) + self.assertEqual(anomaly["affected_instance_ids"], []) + self.assertEqual( + anomaly["affected_candidate_ids"], + [rejected["candidate_id"]], + ) + self.assertEqual(report["overall_status"], "unresolved") + + def test_duplicate_identity(self) -> None: + shape = (48, 48, 112) + duplicate = _ellipsoid(shape, (24.0, 22.0, 30.0)) + duplicate |= _ellipsoid(shape, (24.0, 22.0, 72.0)) + report = self.analyze({NAMES[0]: duplicate}) + self.assertIn("duplicate_identity", _codes(report)) + self.assertEqual( + report["observed_sequence_inferior_to_superior"], + [NAMES[0], NAMES[0]], + ) + + def test_internal_missing_identity(self) -> None: + masks = _standard_masks( + centers_z=(28.0, 76.0), names=(NAMES[0], NAMES[2]) + ) + report = self.analyze(masks) + self.assertIn("missing_internal_identity", _codes(report)) + missing = [ + anomaly + for anomaly in report["anomalies"] + if anomaly["anomaly_code"] == "missing_internal_identity" + ] + self.assertEqual(missing[0]["affected_anatomical_names"], [NAMES[1]]) + + def test_endpoint_missing_labels_are_not_internal_missing(self) -> None: + names = ("vertebrae_L4",) + NAMES + ("vertebrae_T12",) + masks = _standard_masks() + report = self.analyze(masks, names=names) + self.assertNotIn("missing_internal_identity", _codes(report)) + self.assertEqual(report["field_of_view_status"], "extent_uncertain") + + def test_nonmonotonic_identity(self) -> None: + shape = (48, 48, 100) + masks = { + NAMES[1]: _ellipsoid(shape, (24.0, 22.0, 28.0)), + NAMES[0]: _ellipsoid(shape, (24.0, 22.0, 70.0)), + } + report = self.analyze(masks, names=NAMES[:2]) + self.assertIn("nonmonotonic_identity", _codes(report)) + + def test_transitional_duplicate_missing_pattern_is_unresolved(self) -> None: + shape = (48, 48, 120) + masks = { + NAMES[0]: ( + _ellipsoid(shape, (24.0, 22.0, 24.0)) + | _ellipsoid(shape, (24.0, 22.0, 58.0)) + ), + NAMES[2]: _ellipsoid(shape, (24.0, 22.0, 92.0)), + } + report = self.analyze(masks) + self.assertIn("duplicate_identity", _codes(report)) + self.assertIn("missing_internal_identity", _codes(report)) + self.assertEqual(report["overall_status"], "unresolved") + ambiguous = [ + anomaly + for anomaly in report["anomalies"] + if anomaly["anomaly_code"] == "ambiguous_identity" + and "transitional anatomy" in anomaly["explanation"] + ] + self.assertEqual(len(ambiguous), 1) + self.assertEqual(ambiguous[0]["status"], "unresolved") + + def test_mixed_adjacent_identity_is_unresolved(self) -> None: + shape = (48, 48, 96) + body = _ellipsoid(shape, (24.0, 22.0, 48.0)) + x = np.indices(shape)[0] + masks = { + NAMES[0]: body & (x <= 24), + NAMES[1]: body & (x > 24), + } + report = self.analyze(masks, names=NAMES[:2]) + self.assertEqual(len(report["instances"]), 1) + self.assertIsNone(report["instances"][0]["protected_core_identity"]) + self.assertEqual( + report["instances"][0]["status"], "unresolved_mixed_identity" + ) + self.assertIn("ambiguous_identity", _codes(report)) + + def test_overlapping_masks_are_unresolved(self) -> None: + shape = (48, 48, 96) + body = _ellipsoid(shape, (24.0, 22.0, 48.0)) + report = self.analyze( + {NAMES[0]: body, NAMES[1]: body.copy()}, names=NAMES[:2] + ) + self.assertGreater(report["input_overlap_voxel_count"], 0) + self.assertEqual(report["instances"][0]["status"], "unresolved_overlap") + self.assertIn("overlapping_input_masks", _codes(report)) + + def test_partially_overlapping_masks_are_unresolved(self) -> None: + shape = (48, 48, 96) + body = _ellipsoid(shape, (24.0, 22.0, 48.0)) + x = np.indices(shape)[0] + masks = { + NAMES[0]: body & (x <= 26), + NAMES[1]: body & (x >= 22), + } + report = self.analyze(masks, names=NAMES[:2]) + expected_overlap = int( + np.count_nonzero(masks[NAMES[0]] & masks[NAMES[1]]) + ) + self.assertGreater(expected_overlap, 0) + self.assertEqual( + report["input_overlap_voxel_count"], + expected_overlap, + ) + self.assertEqual( + report["instances"][0]["status"], + "unresolved_overlap", + ) + + def test_merged_or_weak_core_separation_is_unresolved(self) -> None: + shape = (56, 56, 112) + lower = _ellipsoid(shape, (28.0, 26.0, 30.0)) + upper = _ellipsoid(shape, (28.0, 26.0, 78.0)) + x, y, z = np.indices(shape) + bridge = ( + ((x - 28.0) ** 2 + (y - 26.0) ** 2 <= 7.0**2) + & (z >= 30) + & (z <= 78) + ) + masks = { + NAMES[0]: lower | (bridge & (z < 54)), + NAMES[1]: upper | (bridge & (z >= 54)), + } + report = self.analyze(masks, names=NAMES[:2]) + self.assertEqual(len(report["instances"]), 1) + self.assertNotEqual( + report["instances"][0]["status"], "protected_high_confidence" + ) + self.assertTrue( + { + "possible_merged_instance", + "multiple_body_profile_peaks_in_one_core", + } + & set(report["instances"][0]["reasons"]) + ) + + def test_thick_posterior_like_candidate_is_unresolved(self) -> None: + shape = (72, 88, 96) + elongated = _ellipsoid( + shape, + (36.0, 44.0, 48.0), + radii=(7.0, 20.0, 8.0), + ) + report = self.analyze( + {NAMES[1]: elongated}, + names=(NAMES[1],), + ) + self.assertEqual(len(report["instances"]), 1) + instance = report["instances"][0] + self.assertNotEqual( + instance["status"], + "protected_high_confidence", + ) + self.assertIsNone(instance["protected_core_identity"]) + self.assertIn( + "core_compactness_below_threshold", + instance["reasons"], + ) + + def test_inferior_partial_field_of_view(self) -> None: + shape = (48, 48, 96) + masks = { + NAMES[0]: _ellipsoid(shape, (24.0, 22.0, 2.0)), + NAMES[1]: _ellipsoid(shape, (24.0, 22.0, 38.0)), + } + report = self.analyze(masks, names=NAMES[:2]) + self.assertEqual( + report["field_of_view_status"], "inferior_boundary_truncated" + ) + self.assertEqual( + report["instances"][0]["status"], "unresolved_boundary_truncated" + ) + self.assertIsNone(report["instances"][0]["protected_core_identity"]) + self.assertNotIn("missing_internal_identity", _codes(report)) + + def test_superior_partial_field_of_view(self) -> None: + shape = (48, 48, 96) + masks = { + NAMES[0]: _ellipsoid(shape, (24.0, 22.0, 52.0)), + NAMES[1]: _ellipsoid(shape, (24.0, 22.0, 94.0)), + } + report = self.analyze(masks, names=NAMES[:2]) + self.assertEqual( + report["field_of_view_status"], "superior_boundary_truncated" + ) + self.assertEqual( + report["instances"][-1]["status"], "unresolved_boundary_truncated" + ) + self.assertIsNone(report["instances"][-1]["protected_core_identity"]) + self.assertNotIn("missing_internal_identity", _codes(report)) + + def test_empty_input(self) -> None: + report = self.analyze({}) + self.assertEqual(report["overall_status"], "empty_input") + self.assertEqual(report["instances"], []) + self.assertEqual(report["shape"], []) + + def test_ct_absent_and_present(self) -> None: + masks = _standard_masks(centers_z=(52.0,), names=(NAMES[1],)) + geometry = self.analyze(masks) + ct = np.full(next(iter(masks.values())).shape, -100.0, dtype=np.float32) + ct[next(iter(masks.values()))] = 300.0 + supported = self.analyze(masks, ct=ct) + self.assertEqual(geometry["ct_evidence"], "unavailable") + self.assertEqual(supported["ct_evidence"], "used") + self.assertEqual(len(geometry["instances"]), len(supported["instances"])) + for geometry_instance, supported_instance in zip( + geometry["instances"], supported["instances"] + ): + for field in ( + "centroid_world_mm", + "core_voxel_count", + "core_volume_mm3", + "maximum_internal_thickness_mm", + "persistence_mm", + "inferior_core_separation_mm", + "superior_core_separation_mm", + ): + self.assertEqual( + geometry_instance[field], supported_instance[field] + ) + self.assertIsNone(geometry["instances"][0]["ct_bone_support_fraction"]) + self.assertGreaterEqual( + supported["instances"][0]["ct_bone_support_fraction"], 0.99 + ) + + def test_ct_shape_mismatch(self) -> None: + masks = _standard_masks(centers_z=(52.0,), names=(NAMES[1],)) + with self.assertRaisesRegex(ValueError, "ct shape differs"): + self.analyze(masks, ct=np.zeros((3, 4, 5), dtype=np.float32)) + + def test_low_ct_bone_support_is_explicit_and_geometry_is_stable( + self, + ) -> None: + masks = _standard_masks( + centers_z=(52.0,), + names=(NAMES[1],), + ) + geometry = self.analyze(masks) + low_ct = np.full( + next(iter(masks.values())).shape, + -100.0, + dtype=np.float32, + ) + low_support = self.analyze(masks, ct=low_ct) + geometry_instance = geometry["instances"][0] + low_instance = low_support["instances"][0] + for field in ( + "centroid_world_mm", + "core_voxel_count", + "core_volume_mm3", + "maximum_internal_thickness_mm", + "persistence_mm", + "inferior_core_separation_mm", + "superior_core_separation_mm", + ): + self.assertEqual( + geometry_instance[field], + low_instance[field], + ) + self.assertEqual(low_support["ct_evidence"], "used") + self.assertEqual( + low_instance["ct_bone_support_fraction"], + 0.0, + ) + self.assertIn("low_ct_bone_support", low_instance["reasons"]) + self.assertNotEqual( + low_instance["status"], + "protected_high_confidence", + ) + + def test_nonfinite_ct_is_rejected(self) -> None: + masks = _standard_masks( + centers_z=(52.0,), + names=(NAMES[1],), + ) + ct = np.zeros( + next(iter(masks.values())).shape, + dtype=np.float32, + ) + ct[0, 0, 0] = np.nan + with self.assertRaisesRegex(ValueError, "non-finite"): + self.analyze(masks, ct=ct) + + def test_anisotropic_spacing_is_physically_stable(self) -> None: + center_worlds = ((24.0, 22.0, 28.0), (24.0, 22.0, 72.0)) + affine_a = np.eye(4, dtype=np.float64) + shape_a = (48, 48, 104) + masks_a = { + name: _physical_ellipsoid(shape_a, affine_a, center) + for name, center in zip(NAMES[:2], center_worlds) + } + affine_b = np.diag([2.0, 1.0, 2.0, 1.0]) + shape_b = (24, 48, 52) + masks_b = { + name: _physical_ellipsoid(shape_b, affine_b, center) + for name, center in zip(NAMES[:2], center_worlds) + } + report_a = self.analyze(masks_a, names=NAMES[:2], affine=affine_a) + report_b = self.analyze(masks_b, names=NAMES[:2], affine=affine_b) + self.assertEqual(_semantic_signature(report_a), _semantic_signature(report_b)) + for first, second in zip(report_a["instances"], report_b["instances"]): + error = np.linalg.norm( + np.asarray(first["centroid_world_mm"]) + - np.asarray(second["centroid_world_mm"]) + ) + self.assertLessEqual(error, np.linalg.norm([2.0, 1.0, 2.0])) + + def test_axis_permutations_and_sign_flips(self) -> None: + base_masks = _standard_masks(centers_z=(30.0, 76.0), names=NAMES[:2]) + base = self.analyze(base_masks, names=NAMES[:2]) + base_centroids = [ + np.asarray(instance["centroid_world_mm"]) for instance in base["instances"] + ] + for permutation in itertools.permutations(range(3)): + for flips in itertools.product((False, True), repeat=3): + transformed: Dict[str, np.ndarray] = {} + transformed_affine = None + for name, mask in base_masks.items(): + new_mask, new_affine = _reorient( + mask, self.affine, permutation, flips + ) + transformed[name] = new_mask + transformed_affine = new_affine + report = self.analyze( + transformed, + names=NAMES[:2], + affine=transformed_affine, + ) + self.assertEqual(_semantic_signature(base), _semantic_signature(report)) + for expected, instance in zip(base_centroids, report["instances"]): + np.testing.assert_allclose( + expected, + instance["centroid_world_mm"], + atol=1e-6, + rtol=0, + ) + + def test_unsupported_affine_shear(self) -> None: + affine = np.eye(4, dtype=np.float64) + affine[0, 1] = 0.25 + report = self.analyze( + _standard_masks(centers_z=(52.0,), names=(NAMES[1],)), + affine=affine, + ) + self.assertEqual(report["overall_status"], "unresolved") + self.assertEqual(_codes(report), ["unsupported_affine"]) + self.assertEqual(report["instances"], []) + + def test_equivalent_29_degree_phantom_is_unresolved(self) -> None: + shape = (64, 64, 80) + affine = _rotation_about_x_affine( + 29.0, + pivot_world=(32.0, 32.0, 40.0), + ) + centers = ((32.0, 32.0, 28.0), (32.0, 32.0, 52.0)) + masks = { + name: _physical_ellipsoid(shape, affine, center) + for name, center in zip(NAMES[:2], centers) + } + report = self.analyze( + masks, + names=NAMES[:2], + affine=affine, + ) + self.assertEqual(_codes(report), ["unsupported_affine"]) + self.assertEqual(report["overall_status"], "unresolved") + self.assertEqual(report["instances"], []) + + def test_invalid_affine_homogeneous_row_raises(self) -> None: + affine = np.eye(4, dtype=np.float64) + affine[3, 0] = 0.2 + with self.assertRaisesRegex(ValueError, "homogeneous row"): + self.analyze( + _standard_masks( + centers_z=(52.0,), + names=(NAMES[1],), + ), + affine=affine, + ) + + def test_orthogonal_in_plane_rotation_has_unique_axis_codes(self) -> None: + angle = math.radians(35.0) + affine = np.asarray( + [ + [math.cos(angle), -math.sin(angle), 0.0, 0.0], + [math.sin(angle), math.cos(angle), 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 1.0], + ], + dtype=np.float64, + ) + report = self.analyze( + _standard_masks(centers_z=(30.0, 76.0), names=NAMES[:2]), + names=NAMES[:2], + affine=affine, + ) + axis_families = [ + "x" if code in ("L", "R") else "y" if code in ("P", "A") else "z" + for code in report["orientation_axcodes"] + ] + self.assertEqual(set(axis_families), {"x", "y", "z"}) + self.assertEqual( + report["observed_sequence_inferior_to_superior"], list(NAMES[:2]) + ) + + def test_equivalent_axis_aligned_and_20_degree_phantoms_match( + self, + ) -> None: + shape = (64, 64, 80) + aligned_affine = np.eye(4, dtype=np.float64) + oblique_affine = _rotation_about_x_affine( + 20.0, + pivot_world=(32.0, 32.0, 40.0), + ) + centers = ((32.0, 32.0, 24.0), (32.0, 32.0, 56.0)) + aligned_masks = { + name: _physical_ellipsoid(shape, aligned_affine, center) + for name, center in zip(NAMES[:2], centers) + } + oblique_masks = { + name: _physical_ellipsoid(shape, oblique_affine, center) + for name, center in zip(NAMES[:2], centers) + } + aligned = self.analyze( + aligned_masks, + names=NAMES[:2], + affine=aligned_affine, + ) + oblique = self.analyze( + oblique_masks, + names=NAMES[:2], + affine=oblique_affine, + ) + self.assertEqual( + _semantic_signature(aligned), + _semantic_signature(oblique), + ) + self.assertNotIn("unsupported_affine", _codes(oblique)) + self.assertEqual( + oblique["observed_sequence_inferior_to_superior"], + list(NAMES[:2]), + ) + + def test_fragmented_volume_reports_every_thick_component(self) -> None: + fragmented = _fragmented_mask() + report = self.analyze( + {NAMES[0]: fragmented}, + names=NAMES[:1], + ) + total_candidates = ( + len(report["instances"]) + + len(report["rejected_candidates"]) + ) + self.assertEqual(total_candidates, 64) + for candidate in report["rejected_candidates"]: + self.assertEqual( + candidate["reasons"], + ["core_outside_trajectory_tube"], + ) + + def test_fragmented_candidate_order_is_deterministic(self) -> None: + fragmented = _fragmented_mask( + shape=(88, 88, 88), + positions=(12.0, 38.0, 64.0), + ) + masks = {NAMES[0]: fragmented} + first = self.analyze(masks, names=NAMES[:1]) + second = self.analyze(masks, names=NAMES[:1]) + self.assertEqual(_canonical(first), _canonical(second)) + accepted_centroids = [ + tuple(item["centroid_world_mm"]) + for item in first["instances"] + ] + rejected_centroids = [ + tuple(item["centroid_world_mm"]) + for item in first["rejected_candidates"] + ] + self.assertEqual( + accepted_centroids, + sorted( + accepted_centroids, + key=lambda value: (value[2], value[1], value[0]), + ), + ) + self.assertEqual( + rejected_centroids, + sorted( + rejected_centroids, + key=lambda value: (value[2], value[1], value[0]), + ), + ) + + def test_all_raised_exception_paths_preserve_inputs(self) -> None: + base_masks = { + name: mask.astype(np.uint8) * (index + 2) + for index, (name, mask) in enumerate( + _standard_masks().items() + ) + } + bad_shape_masks = dict(base_masks) + bad_shape_masks[NAMES[1]] = np.zeros( + (4, 5, 6), + dtype=np.uint8, + ) + bad_dimension_masks = dict(base_masks) + bad_dimension_masks[NAMES[1]] = np.zeros( + (4, 5), + dtype=np.uint8, + ) + nonfinite_ct = np.zeros( + next(iter(base_masks.values())).shape, + dtype=np.float32, + ) + nonfinite_ct[0, 0, 0] = np.inf + bad_row = np.eye(4, dtype=np.float64) + bad_row[3, 1] = 0.5 + singular = np.zeros((4, 4), dtype=np.float64) + singular[3, 3] = 1.0 + cases = ( + ( + "ct_shape", + base_masks, + np.eye(4), + np.zeros((3, 4, 5), dtype=np.float32), + NAMES, + None, + ), + ( + "ct_nonfinite", + base_masks, + np.eye(4), + nonfinite_ct, + NAMES, + None, + ), + ( + "homogeneous_row", + base_masks, + bad_row, + None, + NAMES, + None, + ), + ( + "singular_affine", + base_masks, + singular, + None, + NAMES, + None, + ), + ( + "mask_shape", + bad_shape_masks, + np.eye(4), + None, + NAMES, + None, + ), + ( + "mask_dimension", + bad_dimension_masks, + np.eye(4), + None, + NAMES, + None, + ), + ( + "empty_names", + base_masks, + np.eye(4), + None, + (), + None, + ), + ( + "invalid_config", + base_masks, + np.eye(4), + None, + NAMES, + VertebralInstanceAnalysisConfig(core_radius_mm=-1.0), + ), + ( + "unsupported_configured_obliquity", + base_masks, + np.eye(4), + None, + NAMES, + VertebralInstanceAnalysisConfig( + max_si_axis_obliquity_degrees=29.0 + ), + ), + ) + for ( + case_name, + masks, + affine, + ct, + names, + config, + ) in cases: + arrays = list(masks.values()) + if ct is not None: + arrays.append(ct) + arrays.append(affine) + before = [ + (id(array), array.shape, array.dtype, _array_digest(array)) + for array in arrays + ] + with self.subTest(case=case_name): + with self.assertRaises(ValueError): + self.analyze( + masks, + affine=affine, + ct=ct, + names=names, + config=config, + ) + after = [ + ( + id(array), + array.shape, + array.dtype, + _array_digest(array), + ) + for array in arrays + ] + self.assertEqual(before, after) + + def test_nonfinite_bone_thresholds_preserve_inputs(self) -> None: + masks = { + name: mask.astype(np.uint8) * (index + 2) + for index, (name, mask) in enumerate( + _standard_masks().items() + ) + } + ct = np.full( + next(iter(masks.values())).shape, + 300.0, + dtype=np.float32, + ) + arrays = [*masks.values(), ct] + before = [ + (id(array), array.shape, array.dtype, _array_digest(array)) + for array in arrays + ] + for threshold in ( + float("nan"), + float("inf"), + float("-inf"), + ): + with self.subTest(threshold=threshold): + with self.assertRaisesRegex( + ValueError, + "Configuration values must be finite: " + "bone_hu_threshold", + ): + self.analyze( + masks, + ct=ct, + config=VertebralInstanceAnalysisConfig( + bone_hu_threshold=threshold + ), + ) + after = [ + ( + id(array), + array.shape, + array.dtype, + _array_digest(array), + ) + for array in arrays + ] + self.assertEqual(before, after) + + def test_finite_negative_bone_threshold_is_accepted(self) -> None: + masks = _standard_masks( + centers_z=(52.0,), + names=(NAMES[1],), + ) + ct = np.full( + next(iter(masks.values())).shape, + -100.0, + dtype=np.float32, + ) + report = self.analyze( + masks, + ct=ct, + config=VertebralInstanceAnalysisConfig( + bone_hu_threshold=-200.0 + ), + ) + self.assertEqual( + report["effective_config"]["bone_hu_threshold"], + -200.0, + ) + self.assertEqual( + report["instances"][0]["ct_bone_support_fraction"], + 1.0, + ) + _canonical(report) + + def test_zero_input_mutation(self) -> None: + masks = { + name: mask.astype(np.uint8) * (index + 2) + for index, (name, mask) in enumerate(_standard_masks().items()) + } + before = { + name: (id(array), array.shape, array.dtype, _array_digest(array)) + for name, array in masks.items() + } + self.analyze(masks) + after = { + name: (id(array), array.shape, array.dtype, _array_digest(array)) + for name, array in masks.items() + } + self.assertEqual(before, after) + + def test_twenty_repeated_serializations_are_identical(self) -> None: + masks = _standard_masks() + serializations = {_canonical(self.analyze(masks)) for _ in range(20)} + self.assertEqual(len(serializations), 1) + + def test_numeric_label_independence(self) -> None: + binary = _standard_masks() + low_values = { + name: mask.astype(np.uint8) * (index + 1) + for index, (name, mask) in enumerate(binary.items()) + } + high_values = { + name: mask.astype(np.uint16) * (200 + index) + for index, (name, mask) in enumerate(binary.items()) + } + self.assertEqual( + _canonical(self.analyze(low_values)), + _canonical(self.analyze(high_values)), + ) + + def test_multiprocessing_consistency(self) -> None: + payload = (_standard_masks(), self.affine) + expected = _multiprocessing_worker(payload) + context = multiprocessing.get_context("spawn") + with context.Pool(processes=2) as pool: + observed = pool.map(_multiprocessing_worker, [payload, payload]) + self.assertEqual(observed, [expected, expected]) + + def test_report_contains_only_primitives_and_no_corrected_output( + self, + ) -> None: + report = self.analyze(_standard_masks()) + forbidden_key_terms = ("corrected", "segmentation", "voxel_indices") + + def inspect(value: object) -> None: + self.assertNotIsInstance( + value, + (np.ndarray, np.generic, tuple, set), + ) + if isinstance(value, dict): + for key, child in value.items(): + self.assertIsInstance(key, str) + self.assertFalse( + any(term in key.lower() for term in forbidden_key_terms) + ) + inspect(child) + elif isinstance(value, list): + for child in value: + inspect(child) + else: + if isinstance(value, float): + self.assertTrue(math.isfinite(value)) + self.assertIsInstance( + value, (str, int, float, bool, type(None)) + ) + + inspect(report) + self.assertIsInstance(_canonical(report), str) + + def test_module_has_no_image_writing_api(self) -> None: + module_path = ( + Path(__file__).resolve().parents[1] + / "utils" + / "vertebrae_instance_analysis.py" + ) + source = module_path.read_text(encoding="utf-8") + forbidden = ( + "nib.save", + "nibabel.save", + "to_filename(", + "Nifti1Image(", + "SimpleITK.WriteImage", + ) + for token in forbidden: + self.assertNotIn(token, source) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_vertebrae_instance_audit.py b/tests/test_vertebrae_instance_audit.py new file mode 100644 index 0000000..46de7ac --- /dev/null +++ b/tests/test_vertebrae_instance_audit.py @@ -0,0 +1,1000 @@ +"""Tests for the opt-in vertebral instance batch audit adapter.""" + +from __future__ import annotations + +import ast +import hashlib +import json +import logging +import multiprocessing +import os +import subprocess +import sys +import tempfile +import threading +import unittest +from pathlib import Path +from unittest import mock + +import nibabel as nib +import numpy as np +import yaml + +from utils.vertebrae_instance_audit import ( + VertebraeInstanceAuditConfig, + audit_output_path, + ordered_vertebral_anatomical_names, + parse_vertebrae_instance_audit_config, + run_vertebrae_instance_audit, +) + + +NAMES = ( + "vertebrae_L3", + "vertebrae_L2", + "vertebrae_L1", +) + + +class _ReferenceImage: + def __init__(self, data: np.ndarray) -> None: + self.affine = np.eye(4, dtype=np.float64) + self.dataobj = data + + +class _GeometryOnlyReference: + affine = np.eye(4, dtype=np.float64) + + @property + def dataobj(self): + raise AssertionError("disabled CT data was accessed") + + +def _ellipsoid( + shape=(40, 40, 80), + center=(20.0, 20.0, 40.0), +) -> np.ndarray: + grid = np.ogrid[tuple(slice(0, length) for length in shape)] + normalized = ( + ((grid[0] - center[0]) / 11.0) ** 2 + + ((grid[1] - center[1]) / 12.0) ** 2 + + ((grid[2] - center[2]) / 8.0) ** 2 + ) + return normalized <= 1.0 + + +def _masks() -> dict: + return { + NAMES[0]: _ellipsoid(center=(20.0, 20.0, 24.0)), + NAMES[1]: _ellipsoid(center=(20.0, 20.0, 54.0)), + } + + +def _array_digest(array: np.ndarray) -> str: + hasher = hashlib.sha256() + hasher.update(str(array.shape).encode("ascii")) + hasher.update(array.dtype.str.encode("ascii")) + hasher.update(np.ascontiguousarray(array).tobytes()) + return hasher.hexdigest() + + +def _file_digest(path: Path) -> str: + return hashlib.sha256(path.read_bytes()).hexdigest() + + +_ACTUAL_MAIN_SCRIPT = r""" +import json +import os +import sys +import types + +# scikit-image is optional in the review environment. These import-only +# stubs fail loudly if ShapeKit executes any stubbed operation. target_organs +# is empty in this integration fixture, so no such operation is expected. +skimage = types.ModuleType("skimage") +morphology = types.ModuleType("skimage.morphology") +measure = types.ModuleType("skimage.measure") + +def _unexpected_stub_call(*args, **kwargs): + raise AssertionError("an import-only scikit-image stub was executed") + +morphology.disk = _unexpected_stub_call +morphology.convex_hull_image = _unexpected_stub_call +measure.label = _unexpected_stub_call +measure.regionprops = _unexpected_stub_call +skimage.morphology = morphology +skimage.measure = measure +sys.modules["skimage"] = skimage +sys.modules["skimage.morphology"] = morphology +sys.modules["skimage.measure"] = measure + +sys.path.insert(0, os.environ["SHAPEKIT_REPOSITORY"]) +sys.argv = [ + "main.py", + "--log_folder", + os.environ["SHAPEKIT_LOG_FOLDER"], +] + +import main + +analyzer_before = "utils.vertebrae_instance_analysis" in sys.modules +if os.environ["SHAPEKIT_AUDIT_MODE"] == "failure": + import utils.vertebrae_instance_audit as audit_adapter + + def _synthetic_analyzer_failure(): + raise RuntimeError("synthetic actual-main analyzer failure") + + audit_adapter._load_analyzer = _synthetic_analyzer_failure + +main.main( + os.environ["SHAPEKIT_CASE_PATH"], + "case_001", + os.environ["SHAPEKIT_OUTPUT_ROOT"], +) +analyzer_after = "utils.vertebrae_instance_analysis" in sys.modules +print( + "__SHAPEKIT_TEST_RESULT__" + + json.dumps( + { + "analyzer_before": analyzer_before, + "analyzer_after": analyzer_after, + }, + sort_keys=True, + ) +) +""" + + +def _create_actual_main_case(root: Path) -> Path: + segmentation_folder = root / "input" / "case_001" / "segmentations" + segmentation_folder.mkdir(parents=True) + shape = (40, 40, 80) + affine = np.diag([1.0, 1.0, 1.5, 1.0]) + reference = np.zeros(shape, dtype=np.float32) + reference[1:3, 1:3, 1:3] = 300.0 + vertebra = _ellipsoid(shape=shape, center=(20.0, 20.0, 40.0)) + nib.save( + nib.Nifti1Image(reference, affine), + segmentation_folder / "liver.nii.gz", + ) + nib.save( + nib.Nifti1Image(vertebra.astype(np.uint8), affine), + segmentation_folder / "vertebrae_L1.nii.gz", + ) + return segmentation_folder.parent + + +def _run_actual_main( + *, + repository: Path, + root: Path, + case_path: Path, + mode: str, +) -> dict: + raw_config = yaml.safe_load( + (repository / "config.yaml").read_text(encoding="utf-8") + ) + raw_config["target_organs"] = [] + if mode == "absent": + raw_config.pop("vertebrae_instance_analysis", None) + else: + raw_config["vertebrae_instance_analysis"] = { + "enabled": mode != "disabled", + "output_dir_name": "vertebrae_analysis", + "use_reference_as_ct": mode == "ct", + } + + run_folder = root / f"run_{mode}" + output_root = root / f"output_{mode}" + log_folder = root / f"logs_{mode}" + cache_folder = root / f"cache_{mode}" + run_folder.mkdir() + output_root.mkdir() + (run_folder / "config.yaml").write_text( + yaml.safe_dump(raw_config, sort_keys=False), + encoding="utf-8", + ) + + environment = os.environ.copy() + environment.update( + { + "PYTHONDONTWRITEBYTECODE": "1", + "PYTHONPYCACHEPREFIX": str(cache_folder), + "SHAPEKIT_REPOSITORY": str(repository), + "SHAPEKIT_LOG_FOLDER": str(log_folder), + "SHAPEKIT_AUDIT_MODE": mode, + "SHAPEKIT_CASE_PATH": str(case_path), + "SHAPEKIT_OUTPUT_ROOT": str(output_root), + } + ) + completed = subprocess.run( + [sys.executable, "-B", "-c", _ACTUAL_MAIN_SCRIPT], + cwd=run_folder, + env=environment, + check=False, + capture_output=True, + text=True, + ) + if completed.returncode != 0: + raise AssertionError( + f"actual main.main() failed for {mode}:\n" + f"stdout:\n{completed.stdout}\n" + f"stderr:\n{completed.stderr}" + ) + marker = "__SHAPEKIT_TEST_RESULT__" + result_line = next( + line for line in completed.stdout.splitlines() + if line.startswith(marker) + ) + nifti_hashes = { + str(path.relative_to(output_root)): _file_digest(path) + for path in sorted(output_root.rglob("*.nii.gz")) + } + audit_path = ( + output_root / "vertebrae_analysis" / "case_001.json" + ) + return { + "output_root": output_root, + "nifti_hashes": nifti_hashes, + "audit_path": audit_path, + "process_result": json.loads(result_line[len(marker):]), + "log_folder": log_folder, + } + + +def _spawn_audit_worker(payload) -> str: + output_root, patient_id = payload + masks = _masks() + reference = _GeometryOnlyReference() + status = run_vertebrae_instance_audit( + masks, + reference_img=reference, + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id=patient_id, + config=VertebraeInstanceAuditConfig(enabled=True), + logger=logging.getLogger(f"audit-test-{patient_id}"), + ) + return status + + +class VertebraeInstanceAuditTests(unittest.TestCase): + def test_repository_config_is_disabled_by_default(self) -> None: + repository = Path(__file__).resolve().parents[1] + raw = yaml.safe_load( + (repository / "config.yaml").read_text(encoding="utf-8") + ) + self.assertEqual( + raw["vertebrae_instance_analysis"], + { + "enabled": False, + "output_dir_name": "vertebrae_analysis", + "use_reference_as_ct": False, + }, + ) + self.assertFalse( + parse_vertebrae_instance_audit_config( + raw["vertebrae_instance_analysis"] + ).enabled + ) + + def test_missing_and_disabled_config_do_nothing(self) -> None: + masks = _masks() + reference = _GeometryOnlyReference() + for config in ( + parse_vertebrae_instance_audit_config(None), + VertebraeInstanceAuditConfig(enabled=False), + ): + with self.subTest(config=config): + with tempfile.TemporaryDirectory() as temporary: + output_root = Path(temporary) / "output" + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + side_effect=AssertionError( + "analyzer imported while disabled" + ), + ): + status = run_vertebrae_instance_audit( + masks, + reference_img=reference, + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ) + self.assertEqual(status, "disabled") + self.assertFalse(output_root.exists()) + + def test_enabled_geometry_only_writes_canonical_json(self) -> None: + masks = _masks() + before = { + name: (id(mask), _array_digest(mask)) + for name, mask in masks.items() + } + with tempfile.TemporaryDirectory() as temporary: + output_root = Path(temporary) + config = VertebraeInstanceAuditConfig(enabled=True) + status = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ) + self.assertEqual(status, "written") + path = audit_output_path( + output_root, + config, + "case_001", + ) + payload = path.read_bytes() + report = json.loads(payload) + expected = ( + json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + ensure_ascii=True, + allow_nan=False, + ) + + "\n" + ).encode("utf-8") + self.assertEqual(payload, expected) + self.assertEqual(report["ct_evidence"], "unavailable") + self.assertNotIn("segmentations", path.parts) + after = { + name: (id(mask), _array_digest(mask)) + for name, mask in masks.items() + } + self.assertEqual(before, after) + + def test_ct_is_passed_only_when_explicitly_requested(self) -> None: + masks = _masks() + ct = np.full( + next(iter(masks.values())).shape, + 300.0, + dtype=np.float32, + ) + captured = [] + + def fake_analyzer( + segmentation_dict, + *, + affine, + ordered_anatomical_names, + ct=None, + ): + captured.append(ct) + return {"ct_supplied": ct is not None} + + with tempfile.TemporaryDirectory() as temporary: + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + return_value=fake_analyzer, + ): + for index, use_ct in enumerate((False, True)): + status = run_vertebrae_instance_audit( + masks, + reference_img=_ReferenceImage(ct), + ordered_anatomical_names=NAMES, + output_root=temporary, + patient_id=f"case_{index}", + config=VertebraeInstanceAuditConfig( + enabled=True, + use_reference_as_ct=use_ct, + ), + logger=mock.Mock(), + ) + self.assertEqual(status, "written") + self.assertIsNone(captured[0]) + self.assertIs(captured[1], ct) + + def test_invalid_requested_ct_falls_back_to_geometry(self) -> None: + masks = _masks() + shape = next(iter(masks.values())).shape + invalid_ct_values = ( + np.zeros((3, 4, 5), dtype=np.float32), + np.full(shape, np.nan, dtype=np.float32), + ) + for index, invalid_ct in enumerate(invalid_ct_values): + captured = [] + + def fake_analyzer( + segmentation_dict, + *, + affine, + ordered_anatomical_names, + ct=None, + ): + captured.append(ct) + return {"ct_supplied": ct is not None} + + with self.subTest(index=index): + with tempfile.TemporaryDirectory() as temporary: + logger = mock.Mock() + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + return_value=fake_analyzer, + ): + status = run_vertebrae_instance_audit( + masks, + reference_img=_ReferenceImage(invalid_ct), + ordered_anatomical_names=NAMES, + output_root=temporary, + patient_id=f"case_{index}", + config=VertebraeInstanceAuditConfig( + enabled=True, + use_reference_as_ct=True, + ), + logger=logger, + ) + self.assertEqual(status, "written") + self.assertEqual(captured, [None]) + logger.warning.assert_called() + + def test_failure_preserves_existing_report_and_handoff( + self, + ) -> None: + masks = _masks() + before = { + name: (id(mask), _array_digest(mask)) + for name, mask in masks.items() + } + with tempfile.TemporaryDirectory() as temporary: + config = VertebraeInstanceAuditConfig(enabled=True) + path = audit_output_path(temporary, config, "case_001") + path.parent.mkdir(parents=True) + existing_payload = b'{"previous_success":true}\n' + path.write_bytes(existing_payload) + logger = mock.Mock() + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + side_effect=RuntimeError("synthetic analysis failure"), + ): + status = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=temporary, + patient_id="case_001", + config=config, + logger=logger, + ) + self.assertEqual(status, "failed") + self.assertEqual(path.read_bytes(), existing_payload) + self.assertEqual(list(path.parent.glob("*.tmp")), []) + logger.exception.assert_called() + + handoff = {} + + def downstream(segmentation_dict): + handoff["object"] = segmentation_dict + handoff["hashes"] = { + name: _array_digest(mask) + for name, mask in segmentation_dict.items() + } + + downstream(masks) + self.assertIs(handoff["object"], masks) + self.assertEqual( + before, + { + name: (id(mask), handoff["hashes"][name]) + for name, mask in masks.items() + }, + ) + + def test_write_failure_preserves_final_and_removes_temporary_output( + self, + ) -> None: + masks = _masks() + with tempfile.TemporaryDirectory() as temporary: + config = VertebraeInstanceAuditConfig(enabled=True) + path = audit_output_path(temporary, config, "case_001") + path.parent.mkdir(parents=True) + existing_payload = b'{"previous_success":true}\n' + path.write_bytes(existing_payload) + with mock.patch( + "utils.vertebrae_instance_audit.os.replace", + side_effect=OSError("synthetic write failure"), + ): + status = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=temporary, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ) + self.assertEqual(status, "failed") + self.assertEqual(path.read_bytes(), existing_payload) + self.assertEqual(list(path.parent.glob("*.tmp")), []) + + def test_same_patient_concurrent_failure_cannot_delete_success( + self, + ) -> None: + masks = _masks() + config = VertebraeInstanceAuditConfig(enabled=True) + failure_started = threading.Event() + success_finished = threading.Event() + statuses = {} + logger = mock.Mock() + + def coordinated_analyzer(*args, **kwargs): + if threading.current_thread().name == "failing-audit": + failure_started.set() + if not success_finished.wait(timeout=10): + raise AssertionError("successful worker did not finish") + raise RuntimeError("synthetic concurrent failure") + if not failure_started.wait(timeout=10): + raise AssertionError("failing worker did not start") + return {"worker": "successful-audit"} + + with tempfile.TemporaryDirectory() as temporary: + output_root = Path(temporary) + path = audit_output_path( + output_root, + config, + "case_001", + ) + + def successful_worker(): + statuses["success"] = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=logger, + ) + success_finished.set() + + def failing_worker(): + statuses["failure"] = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=logger, + ) + + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + return_value=coordinated_analyzer, + ): + failure_thread = threading.Thread( + target=failing_worker, + name="failing-audit", + ) + success_thread = threading.Thread( + target=successful_worker, + name="successful-audit", + ) + failure_thread.start() + success_thread.start() + failure_thread.join(timeout=15) + success_thread.join(timeout=15) + + self.assertFalse(failure_thread.is_alive()) + self.assertFalse(success_thread.is_alive()) + self.assertEqual( + statuses, + {"success": "written", "failure": "failed"}, + ) + payload = path.read_bytes() + report = json.loads(payload) + self.assertEqual(report, {"worker": "successful-audit"}) + self.assertEqual( + payload, + b'{"worker":"successful-audit"}\n', + ) + self.assertEqual(list(path.parent.glob("*.tmp")), []) + + def test_later_failure_preserves_same_and_other_patient_reports( + self, + ) -> None: + masks = _masks() + config = VertebraeInstanceAuditConfig(enabled=True) + with tempfile.TemporaryDirectory() as temporary: + output_root = Path(temporary) + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + return_value=lambda *args, **kwargs: {"success": True}, + ): + for patient_id in ("case_001", "case_002"): + self.assertEqual( + run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id=patient_id, + config=config, + logger=mock.Mock(), + ), + "written", + ) + first_path = audit_output_path( + output_root, + config, + "case_001", + ) + second_path = audit_output_path( + output_root, + config, + "case_002", + ) + first_payload = first_path.read_bytes() + second_payload = second_path.read_bytes() + + with mock.patch( + "utils.vertebrae_instance_audit._load_analyzer", + side_effect=RuntimeError("later failure"), + ): + self.assertEqual( + run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ), + "failed", + ) + self.assertEqual(first_path.read_bytes(), first_payload) + self.assertEqual(second_path.read_bytes(), second_payload) + self.assertEqual(list(first_path.parent.glob("*.tmp")), []) + + def test_symlinked_audit_directories_are_rejected(self) -> None: + masks = _masks() + before = { + name: _array_digest(mask) + for name, mask in masks.items() + } + config = VertebraeInstanceAuditConfig(enabled=True) + for destination_kind in ("segmentation", "outside"): + with self.subTest(destination_kind=destination_kind): + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + output_root = root / "output" + output_root.mkdir() + if destination_kind == "segmentation": + target = ( + output_root + / "case_001" + / "segmentations" + ) + else: + target = root / "outside" + target.mkdir(parents=True) + audit_directory = ( + output_root / config.output_dir_name + ) + try: + audit_directory.symlink_to( + target, + target_is_directory=True, + ) + except (NotImplementedError, OSError) as error: + self.skipTest( + f"symbolic links are unavailable: {error}" + ) + + status = run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ) + self.assertEqual(status, "failed") + self.assertEqual(list(target.glob("*.json")), []) + self.assertEqual(list(target.glob("*.tmp")), []) + self.assertEqual( + before, + { + name: _array_digest(mask) + for name, mask in masks.items() + }, + ) + + def test_existing_real_audit_directory_is_accepted(self) -> None: + masks = _masks() + config = VertebraeInstanceAuditConfig(enabled=True) + with tempfile.TemporaryDirectory() as temporary: + output_root = Path(temporary) + audit_directory = ( + output_root / config.output_dir_name + ) + audit_directory.mkdir() + self.assertEqual( + run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=output_root, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ), + "written", + ) + self.assertTrue((audit_directory / "case_001.json").is_file()) + + def test_repeated_runs_write_byte_identical_json(self) -> None: + masks = _masks() + with tempfile.TemporaryDirectory() as temporary: + config = VertebraeInstanceAuditConfig(enabled=True) + path = audit_output_path(temporary, config, "case_001") + payloads = [] + for _ in range(2): + self.assertEqual( + run_vertebrae_instance_audit( + masks, + reference_img=_GeometryOnlyReference(), + ordered_anatomical_names=NAMES, + output_root=temporary, + patient_id="case_001", + config=config, + logger=mock.Mock(), + ), + "written", + ) + payloads.append(path.read_bytes()) + self.assertEqual(payloads[0], payloads[1]) + + def test_spawned_cases_use_distinct_deterministic_paths(self) -> None: + with tempfile.TemporaryDirectory() as temporary: + payloads = [ + (temporary, "case_001"), + (temporary, "case_002"), + ] + context = multiprocessing.get_context("spawn") + with context.Pool(processes=2) as pool: + statuses = pool.map(_spawn_audit_worker, payloads) + self.assertEqual(statuses, ["written", "written"]) + output_dir = Path(temporary) / "vertebrae_analysis" + paths = sorted(output_dir.glob("*.json")) + self.assertEqual( + [path.name for path in paths], + ["case_001.json", "case_002.json"], + ) + for path in paths: + payload = path.read_bytes() + report = json.loads(payload) + self.assertEqual( + payload, + ( + json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + ensure_ascii=True, + allow_nan=False, + ) + + "\n" + ).encode("utf-8"), + ) + + def test_ordering_uses_names_not_numeric_class_keys(self) -> None: + first = { + 900: "vertebrae_C1", + -4: "vertebrae_L1", + 12: "liver", + 3: "vertebrae_T12", + 1: "vertebrae_L5", + } + second = { + 1: "vertebrae_T12", + 800: "vertebrae_L5", + 2: "vertebrae_C1", + 44: "vertebrae_L1", + 9: "liver", + } + expected = ( + "vertebrae_L5", + "vertebrae_L1", + "vertebrae_T12", + "vertebrae_C1", + ) + self.assertEqual( + ordered_vertebral_anatomical_names(first.values()), + expected, + ) + self.assertEqual( + ordered_vertebral_anatomical_names(second.values()), + expected, + ) + + def test_malformed_configuration_is_rejected(self) -> None: + invalid = ( + [], + {"enabled": "false"}, + {"use_reference_as_ct": 1}, + {"output_dir_name": ""}, + {"output_dir_name": ".."}, + {"output_dir_name": "/absolute"}, + {"output_dir_name": "nested/path"}, + {"output_dir_name": " leading"}, + {"unknown": True}, + ) + for raw in invalid: + with self.subTest(raw=raw): + with self.assertRaises(ValueError): + parse_vertebrae_instance_audit_config(raw) + + def test_patient_identifiers_with_control_characters_are_rejected( + self, + ) -> None: + config = VertebraeInstanceAuditConfig(enabled=True) + for patient_id in ( + "case\ninjected", + "case\tinjected", + "case\x1finjected", + "case\u0085injected", + ): + with self.subTest(patient_id=repr(patient_id)): + with self.assertRaises(ValueError): + audit_output_path(".", config, patient_id) + + def test_main_calls_audit_before_combination_as_unused_result( + self, + ) -> None: + repository = Path(__file__).resolve().parents[1] + tree = ast.parse( + (repository / "main.py").read_text(encoding="utf-8") + ) + main_function = next( + node + for node in tree.body + if isinstance(node, ast.FunctionDef) and node.name == "main" + ) + + call_positions = {} + audit_statement = None + for statement_index, statement in enumerate(main_function.body): + for node in ast.walk(statement): + if not isinstance(node, ast.Call): + continue + function = node.func + if isinstance(function, ast.Name): + call_positions.setdefault( + function.id, + statement_index, + ) + if function.id == "run_vertebrae_instance_audit": + audit_statement = statement + + self.assertLess( + call_positions["read_all_segmentations"], + call_positions["run_vertebrae_instance_audit"], + ) + self.assertLess( + call_positions["run_vertebrae_instance_audit"], + call_positions["combine_segmentation_dict"], + ) + self.assertLess( + call_positions["combine_segmentation_dict"], + call_positions["process_organs"], + ) + self.assertIsInstance(audit_statement, ast.Expr) + + def test_actual_main_modes_preserve_segmentation_output_bytes( + self, + ) -> None: + repository = Path(__file__).resolve().parents[1] + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + case_path = _create_actual_main_case(root) + input_hashes = { + str(path.relative_to(case_path)): _file_digest(path) + for path in sorted(case_path.rglob("*.nii.gz")) + } + results = { + mode: _run_actual_main( + repository=repository, + root=root, + case_path=case_path, + mode=mode, + ) + for mode in ( + "absent", + "disabled", + "geometry", + "ct", + "failure", + ) + } + + expected_hashes = results["absent"]["nifti_hashes"] + self.assertTrue(expected_hashes) + for mode, result in results.items(): + with self.subTest(mode=mode): + self.assertEqual( + result["nifti_hashes"], + expected_hashes, + ) + + for mode in ("absent", "disabled"): + result = results[mode] + self.assertFalse( + ( + result["output_root"] + / "vertebrae_analysis" + ).exists() + ) + self.assertFalse( + result["process_result"]["analyzer_before"] + ) + self.assertFalse( + result["process_result"]["analyzer_after"] + ) + + for mode, expected_ct_evidence in ( + ("geometry", "unavailable"), + ("ct", "used"), + ): + audit_path = results[mode]["audit_path"] + payload = audit_path.read_bytes() + report = json.loads(payload) + self.assertEqual( + report["ct_evidence"], + expected_ct_evidence, + ) + self.assertEqual( + payload, + ( + json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + ensure_ascii=True, + allow_nan=False, + ) + + "\n" + ).encode("utf-8"), + ) + + self.assertFalse(results["failure"]["audit_path"].exists()) + self.assertEqual( + list( + results["failure"]["audit_path"].parent.glob( + "*.tmp" + ) + ), + [], + ) + failure_log = ( + results["failure"]["log_folder"] / "debug.log" + ).read_text(encoding="utf-8") + self.assertIn( + "synthetic actual-main analyzer failure", + failure_log, + ) + self.assertIn("patient=case_001", failure_log) + self.assertIn("report=", failure_log) + self.assertIn("requested_ct_mode=geometry-only", failure_log) + self.assertIn("effective_ct_mode=geometry-only", failure_log) + + self.assertEqual( + input_hashes, + { + str(path.relative_to(case_path)): _file_digest(path) + for path in sorted(case_path.rglob("*.nii.gz")) + }, + ) + + +if __name__ == "__main__": + unittest.main() diff --git a/utils/vertebrae_instance_analysis.py b/utils/vertebrae_instance_analysis.py new file mode 100644 index 0000000..0a08704 --- /dev/null +++ b/utils/vertebrae_instance_analysis.py @@ -0,0 +1,1238 @@ +"""Read-only physical vertebral-instance and sequence analysis. + +This module deliberately does not post-process segmentations. It accepts +binary masks keyed by anatomical name and returns a deterministic, +JSON-compatible diagnostic report. No function in this module writes image +files or returns a corrected label map. + +``ordered_anatomical_names`` is always interpreted inferior-to-superior. +Geometry is measured in physical millimetres from the supplied affine. CT is +optional and, when present, contributes confidence evidence only; it never +changes instance boundaries. +""" + +from __future__ import annotations + +import itertools +import math +from dataclasses import asdict, dataclass +from typing import Dict, List, Literal, Mapping, Optional, Sequence, Tuple, TypedDict + +import numpy as np +from scipy import ndimage +from scipy import signal as scipy_signal + + +class NameCount(TypedDict): + anatomical_name: str + voxel_count: int + fraction_of_core: float + + +class VertebralInstanceReport(TypedDict): + instance_id: str + rank_inferior_to_superior: int + centroid_world_mm: List[float] + core_voxel_count: int + core_volume_mm3: float + maximum_internal_thickness_mm: float + persistence_mm: float + compactness: float + trajectory_distance_mm: float + current_name_composition: List[NameCount] + protected_core_identity: Optional[str] + identity_confidence: float + confidence_mode: Literal["geometry_only", "geometry_and_ct"] + ct_bone_support_fraction: Optional[float] + inferior_core_separation_mm: Optional[float] + superior_core_separation_mm: Optional[float] + status: Literal[ + "protected_high_confidence", + "unresolved_low_confidence", + "unresolved_mixed_identity", + "unresolved_overlap", + "unresolved_boundary_truncated", + ] + reasons: List[str] + + +class SequenceAnomalyReport(TypedDict): + anomaly_code: Literal[ + "duplicate_identity", + "missing_internal_identity", + "nonmonotonic_identity", + "off_trajectory_core", + "abnormal_spacing", + "ambiguous_identity", + "overlapping_input_masks", + "unsupported_affine", + ] + affected_instance_ids: List[str] + affected_candidate_ids: List[str] + affected_anatomical_names: List[str] + status: Literal["detected", "unresolved"] + explanation: str + + +class RejectedCandidateReport(TypedDict): + candidate_id: str + centroid_world_mm: List[float] + core_voxel_count: int + core_volume_mm3: float + trajectory_distance_mm: float + status: Literal["rejected_unresolved"] + reasons: List[str] + + +class VertebralAnalysisReport(TypedDict): + schema_version: str + ordered_names_direction: Literal["inferior_to_superior"] + shape: List[int] + spacing_mm: List[float] + orientation_axcodes: List[str] + ct_evidence: Literal["used", "unavailable"] + field_of_view_status: Literal[ + "not_truncated_at_array_boundary", + "inferior_boundary_truncated", + "superior_boundary_truncated", + "both_boundaries_truncated", + "extent_uncertain", + ] + effective_config: Dict[str, float] + input_overlap_voxel_count: int + instances: List[VertebralInstanceReport] + rejected_candidates: List[RejectedCandidateReport] + observed_sequence_inferior_to_superior: List[Optional[str]] + anomalies: List[SequenceAnomalyReport] + overall_status: Literal[ + "continuous_sequence", + "anomaly_detected", + "unresolved", + "empty_input", + ] + + +@dataclass(frozen=True) +class VertebralInstanceAnalysisConfig: + """Configurable physical and confidence thresholds. + + Defaults are conservative starting values for diagnostics and synthetic + validation. They are not population-level performance claims. + """ + + trajectory_tube_radius_mm: float = 55.0 + trajectory_outlier_mm: float = 35.0 + core_radius_mm: float = 5.0 + min_core_volume_mm3: float = 120.0 + min_core_persistence_mm: float = 5.0 + max_core_persistence_mm: float = 45.0 + profile_smoothing_mm: float = 2.0 + min_instance_peak_distance_mm: float = 14.0 + min_instance_peak_prominence: float = 0.08 + min_compactness: float = 0.30 + max_trajectory_distance_mm: float = 30.0 + min_label_vote_fraction: float = 0.55 + min_label_vote_margin: float = 0.15 + geometry_confidence_threshold: float = 0.75 + ct_augmented_confidence_threshold: float = 0.70 + bone_hu_threshold: float = 150.0 + min_bone_support_fraction: float = 0.35 + spacing_outlier_mad: float = 3.5 + max_si_axis_obliquity_degrees: float = 20.0 + + +@dataclass +class _Trajectory: + x_coefficients: np.ndarray + y_coefficients: np.ndarray + + def point_at_world_z(self, world_z: float) -> np.ndarray: + return np.asarray( + [ + np.polyval(self.x_coefficients, world_z), + np.polyval(self.y_coefficients, world_z), + world_z, + ], + dtype=np.float64, + ) + + +@dataclass +class _Candidate: + centroid_world: np.ndarray + core_voxel_count: int + core_volume_mm3: float + maximum_internal_thickness_mm: float + persistence_mm: float + compactness: float + trajectory_distance_mm: float + composition: List[NameCount] + proposed_identity: Optional[str] + confidence: float + ct_bone_support_fraction: Optional[float] + status: str + reasons: List[str] + min_world_z: float + max_world_z: float + center_si_index: float + inferior_core_separation_mm: Optional[float] = None + superior_core_separation_mm: Optional[float] = None + + +@dataclass +class _RejectedCandidate: + centroid_world: np.ndarray + core_voxel_count: int + core_volume_mm3: float + trajectory_distance_mm: float + reasons: List[str] + + +_SCHEMA_VERSION = "1.1" +_ORTHOGONALITY_TOLERANCE = 1e-3 +_MAX_SUPPORTED_SI_OBLIQUITY_DEGREES = 20.0 +_ROUND_DECIMALS = 6 + + +def _rounded(value: float) -> float: + return round(float(value), _ROUND_DECIMALS) + + +def _rounded_vector(values: np.ndarray) -> List[float]: + return [_rounded(value) for value in np.asarray(values).tolist()] + + +def _validate_config(config: VertebralInstanceAnalysisConfig) -> None: + if not np.isfinite(config.bone_hu_threshold): + raise ValueError( + "Configuration values must be finite: bone_hu_threshold" + ) + positive = { + "trajectory_tube_radius_mm": config.trajectory_tube_radius_mm, + "trajectory_outlier_mm": config.trajectory_outlier_mm, + "core_radius_mm": config.core_radius_mm, + "min_core_volume_mm3": config.min_core_volume_mm3, + "min_core_persistence_mm": config.min_core_persistence_mm, + "max_core_persistence_mm": config.max_core_persistence_mm, + "profile_smoothing_mm": config.profile_smoothing_mm, + "min_instance_peak_distance_mm": config.min_instance_peak_distance_mm, + "max_trajectory_distance_mm": config.max_trajectory_distance_mm, + "spacing_outlier_mad": config.spacing_outlier_mad, + } + invalid = [name for name, value in positive.items() if not np.isfinite(value) or value <= 0] + if invalid: + raise ValueError("Configuration values must be positive: " + ", ".join(invalid)) + fractions = { + "min_instance_peak_prominence": config.min_instance_peak_prominence, + "min_compactness": config.min_compactness, + "min_label_vote_fraction": config.min_label_vote_fraction, + "min_label_vote_margin": config.min_label_vote_margin, + "geometry_confidence_threshold": config.geometry_confidence_threshold, + "ct_augmented_confidence_threshold": config.ct_augmented_confidence_threshold, + "min_bone_support_fraction": config.min_bone_support_fraction, + } + invalid = [ + name + for name, value in fractions.items() + if not np.isfinite(value) or value < 0 or value > 1 + ] + if invalid: + raise ValueError("Configuration values must lie in [0, 1]: " + ", ".join(invalid)) + if config.max_core_persistence_mm < config.min_core_persistence_mm: + raise ValueError("max_core_persistence_mm must be >= min_core_persistence_mm") + if not ( + 0 + < config.max_si_axis_obliquity_degrees + <= _MAX_SUPPORTED_SI_OBLIQUITY_DEGREES + ): + raise ValueError( + "max_si_axis_obliquity_degrees must lie in (0, 20]" + ) + + +def _validate_affine( + affine: np.ndarray, config: VertebralInstanceAnalysisConfig +) -> Tuple[np.ndarray, np.ndarray, List[str], int, int, Optional[str]]: + matrix = np.asarray(affine, dtype=np.float64) + if matrix.shape != (4, 4): + raise ValueError("affine must be a 4x4 matrix") + if not np.all(np.isfinite(matrix)): + raise ValueError("affine contains non-finite values") + expected_homogeneous_row = np.asarray([0.0, 0.0, 0.0, 1.0]) + if not np.allclose( + matrix[3, :], + expected_homogeneous_row, + rtol=0.0, + atol=1e-8, + ): + raise ValueError( + "affine homogeneous row must be approximately [0, 0, 0, 1]" + ) + linear = matrix[:3, :3] + determinant = float(np.linalg.det(linear)) + if abs(determinant) < 1e-8: + raise ValueError("affine is singular") + spacing = np.linalg.norm(linear, axis=0) + if np.any(spacing <= 0): + raise ValueError("affine contains a zero-length voxel axis") + directions = linear / spacing + gram = directions.T @ directions + off_diagonal = gram - np.eye(3) + shear = float(np.max(np.abs(off_diagonal))) + orientation = _orientation_axcodes(directions) + si_axis = int(np.argmax(np.abs(directions[2, :]))) + si_sign = 1 if directions[2, si_axis] >= 0 else -1 + alignment = float(np.clip(abs(directions[2, si_axis]), 0.0, 1.0)) + obliquity_degrees = math.degrees(math.acos(alignment)) + unsupported_reason = None + if shear > _ORTHOGONALITY_TOLERANCE: + unsupported_reason = ( + "Affine voxel axes are materially non-orthogonal " + f"(maximum normalized dot product {shear:.6f})." + ) + elif ( + obliquity_degrees + > config.max_si_axis_obliquity_degrees + 1e-6 + ): + unsupported_reason = ( + "No voxel axis is sufficiently aligned with physical superior-inferior " + f"direction (obliquity {obliquity_degrees:.3f} degrees)." + ) + return matrix, spacing, orientation, si_axis, si_sign, unsupported_reason + + +def _orientation_axcodes(directions: np.ndarray) -> List[str]: + negative = ("L", "P", "I") + positive = ("R", "A", "S") + # Assign each voxel axis to a unique world axis. Independent argmax calls + # can emit duplicate codes for valid in-plane rotations near 45 degrees. + assignment = max( + itertools.permutations(range(3)), + key=lambda candidate: sum( + abs(float(directions[candidate[voxel_axis], voxel_axis])) + for voxel_axis in range(3) + ), + ) + codes: List[str] = [] + for voxel_axis in range(3): + world_axis = assignment[voxel_axis] + sign_positive = directions[world_axis, voxel_axis] >= 0 + codes.append(positive[world_axis] if sign_positive else negative[world_axis]) + return codes + + +def _shape_and_masks( + segmentation_dict: Mapping[str, np.ndarray], + ordered_names: Sequence[str], +) -> Tuple[List[int], Dict[str, np.ndarray]]: + requested: Dict[str, np.ndarray] = {} + shape: Optional[Tuple[int, ...]] = None + for name in ordered_names: + if name not in segmentation_dict: + continue + array = np.asarray(segmentation_dict[name]) + if array.ndim != 3: + raise ValueError(f"Mask {name!r} is not three-dimensional") + if shape is None: + shape = array.shape + elif array.shape != shape: + raise ValueError(f"Mask {name!r} shape differs from other vertebral masks") + requested[name] = array + if shape is None: + for value in segmentation_dict.values(): + array = np.asarray(value) + if array.ndim == 3: + shape = array.shape + break + return list(shape) if shape is not None else [], requested + + +def _apply_affine(affine: np.ndarray, coordinates: np.ndarray) -> np.ndarray: + coords = np.asarray(coordinates, dtype=np.float64) + return coords @ affine[:3, :3].T + affine[:3, 3] + + +def _largest_component_centroid(mask: np.ndarray) -> Optional[np.ndarray]: + structure = ndimage.generate_binary_structure(3, 1) + components, count = ndimage.label(mask, structure=structure) + if count == 0: + return None + sizes = np.bincount(components.ravel()) + sizes[0] = 0 + component_id = int(np.argmax(sizes)) + centroid = ndimage.center_of_mass(mask, components, component_id) + return np.asarray(centroid, dtype=np.float64) + + +def _fit_trajectory( + masks: Mapping[str, np.ndarray], + ordered_names: Sequence[str], + affine: np.ndarray, + outlier_mm: float, +) -> Optional[_Trajectory]: + seeds: List[np.ndarray] = [] + for name in ordered_names: + mask = masks.get(name) + if mask is None or not np.any(mask): + continue + centroid = _largest_component_centroid(mask != 0) + if centroid is not None: + seeds.append(_apply_affine(affine, centroid[None, :])[0]) + if not seeds: + return None + points = np.asarray(seeds, dtype=np.float64) + unique_z = np.unique(np.round(points[:, 2], decimals=6)) + degree = min(2, len(unique_z) - 1) + if degree <= 0: + return _Trajectory( + x_coefficients=np.asarray([float(np.median(points[:, 0]))]), + y_coefficients=np.asarray([float(np.median(points[:, 1]))]), + ) + x_coefficients = np.polyfit(points[:, 2], points[:, 0], degree) + y_coefficients = np.polyfit(points[:, 2], points[:, 1], degree) + predicted = np.column_stack( + [ + np.polyval(x_coefficients, points[:, 2]), + np.polyval(y_coefficients, points[:, 2]), + ] + ) + residuals = np.linalg.norm(points[:, :2] - predicted, axis=1) + keep = residuals <= outlier_mm + if int(np.count_nonzero(keep)) >= degree + 1 and not np.all(keep): + kept = points[keep] + kept_unique_z = np.unique(np.round(kept[:, 2], decimals=6)) + kept_degree = min(degree, len(kept_unique_z) - 1) + if kept_degree > 0: + x_coefficients = np.polyfit(kept[:, 2], kept[:, 0], kept_degree) + y_coefficients = np.polyfit(kept[:, 2], kept[:, 1], kept_degree) + return _Trajectory(x_coefficients, y_coefficients) + + +def _bbox(mask: np.ndarray) -> Tuple[slice, slice, slice]: + coordinates = np.argwhere(mask) + low = np.min(coordinates, axis=0) + high = np.max(coordinates, axis=0) + 1 + return tuple( # type: ignore[return-value] + slice(int(low[axis]), int(high[axis])) for axis in range(3) + ) + + +def _compactness(world: np.ndarray) -> float: + if len(world) < 5: + return 0.0 + transverse = np.asarray(world[:, :2], dtype=np.float64) + covariance = np.cov(transverse, rowvar=False) + eigenvalues = np.linalg.eigvalsh(covariance) + if eigenvalues[-1] <= 1e-8: + return 0.0 + return float(math.sqrt(max(float(eigenvalues[0]), 0.0) / float(eigenvalues[-1]))) + + +def _clip_ratio(value: float, reference: float) -> float: + if reference <= 0: + return 0.0 + return float(np.clip(value / reference, 0.0, 1.0)) + + +def _confidence( + *, + max_radius_mm: float, + volume_mm3: float, + persistence_mm: float, + compactness: float, + trajectory_distance_mm: float, + vote_fraction: float, + vote_margin: float, + bone_support: Optional[float], + config: VertebralInstanceAnalysisConfig, +) -> float: + thickness_score = _clip_ratio( + max_radius_mm - config.core_radius_mm, + config.core_radius_mm, + ) + volume_score = _clip_ratio(volume_mm3, 2.0 * config.min_core_volume_mm3) + persistence_score = _clip_ratio( + persistence_mm, config.min_core_persistence_mm + ) + compactness_score = _clip_ratio(compactness, 2.0 * config.min_compactness) + trajectory_score = 1.0 - _clip_ratio( + trajectory_distance_mm, config.max_trajectory_distance_mm + ) + vote_score = 0.65 * vote_fraction + 0.35 * _clip_ratio( + vote_margin, config.min_label_vote_margin + ) + geometry = ( + 0.15 * thickness_score + + 0.20 * volume_score + + 0.15 * persistence_score + + 0.15 * compactness_score + + 0.15 * trajectory_score + + 0.20 * vote_score + ) + if bone_support is None: + return float(np.clip(geometry, 0.0, 1.0)) + bone_score = _clip_ratio(bone_support, config.min_bone_support_fraction) + return float(np.clip(0.85 * geometry + 0.15 * bone_score, 0.0, 1.0)) + + +def _profile( + thick: np.ndarray, + spacing: np.ndarray, + si_axis: int, + config: VertebralInstanceAnalysisConfig, +) -> Tuple[np.ndarray, np.ndarray]: + transverse_axes = tuple(axis for axis in range(3) if axis != si_axis) + area = np.count_nonzero(thick, axis=transverse_axes).astype(np.float64) + scale = max(float(np.percentile(area, 95)), 1.0) + normalized = np.clip(area / scale, 0.0, 1.5) + sigma = config.profile_smoothing_mm / max(float(spacing[si_axis]), 1e-6) + smoothed = ndimage.gaussian_filter1d( + normalized, sigma=max(sigma, 0.5), mode="nearest" + ) + minimum_distance = max( + 1, + int( + round( + config.min_instance_peak_distance_mm + / max(float(spacing[si_axis]), 1e-6) + ) + ), + ) + peaks, _ = scipy_signal.find_peaks( + smoothed, + distance=minimum_distance, + prominence=config.min_instance_peak_prominence, + ) + return smoothed, np.asarray(sorted(int(peak) for peak in peaks), dtype=int) + + +def _core_valley_quality( + profile: np.ndarray, left_index: float, right_index: float +) -> float: + low = int(round(min(left_index, right_index))) + high = int(round(max(left_index, right_index))) + low = max(low, 0) + high = min(high, len(profile) - 1) + if high <= low: + return 0.0 + valley = float(np.min(profile[low : high + 1])) + endpoint = min(float(profile[low]), float(profile[high])) + if endpoint <= 1e-8: + return 0.0 + return float(np.clip(1.0 - valley / endpoint, 0.0, 1.0)) + + +def _field_of_view_status( + union: np.ndarray, + si_axis: int, + si_sign: int, + protected_identities: Sequence[Optional[str]], + ordered_names: Sequence[str], +) -> str: + inferior_touch, superior_touch = _si_boundary_touches( + union, si_axis, si_sign + ) + if inferior_touch and superior_touch: + return "both_boundaries_truncated" + if inferior_touch: + return "inferior_boundary_truncated" + if superior_touch: + return "superior_boundary_truncated" + confident = [name for name in protected_identities if name is not None] + if confident and ordered_names: + order = {name: index for index, name in enumerate(ordered_names)} + observed_indices = [order[name] for name in confident] + if min(observed_indices) > 0 or max(observed_indices) < len(ordered_names) - 1: + return "extent_uncertain" + return "not_truncated_at_array_boundary" + + +def _si_boundary_touches( + union: np.ndarray, si_axis: int, si_sign: int +) -> Tuple[bool, bool]: + inferior_index = 0 if si_sign > 0 else union.shape[si_axis] - 1 + superior_index = union.shape[si_axis] - 1 if si_sign > 0 else 0 + inferior_touch = bool( + np.any(np.take(union, indices=inferior_index, axis=si_axis)) + ) + superior_touch = bool( + np.any(np.take(union, indices=superior_index, axis=si_axis)) + ) + return inferior_touch, superior_touch + + +def _anomaly( + code: str, + instance_ids: Sequence[str], + names: Sequence[str], + status: str, + explanation: str, + candidate_ids: Sequence[str] = (), +) -> SequenceAnomalyReport: + return { + "anomaly_code": code, # type: ignore[typeddict-item] + "affected_instance_ids": list(instance_ids), + "affected_candidate_ids": list(candidate_ids), + "affected_anatomical_names": list(names), + "status": status, # type: ignore[typeddict-item] + "explanation": explanation, + } + + +def _empty_report( + *, + shape: Sequence[int], + spacing: np.ndarray, + orientation: Sequence[str], + ct_used: bool, + config: VertebralInstanceAnalysisConfig, + anomaly: Optional[SequenceAnomalyReport] = None, +) -> VertebralAnalysisReport: + anomalies = [] if anomaly is None else [anomaly] + return { + "schema_version": _SCHEMA_VERSION, + "ordered_names_direction": "inferior_to_superior", + "shape": [int(value) for value in shape], + "spacing_mm": _rounded_vector(spacing), + "orientation_axcodes": list(orientation), + "ct_evidence": "used" if ct_used else "unavailable", + "field_of_view_status": "extent_uncertain", + "effective_config": { + key: _rounded(value) for key, value in asdict(config).items() + }, + "input_overlap_voxel_count": 0, + "instances": [], + "rejected_candidates": [], + "observed_sequence_inferior_to_superior": [], + "anomalies": anomalies, + "overall_status": "unresolved" if anomalies else "empty_input", + } + + +def analyze_vertebral_instances( + segmentation_dict: Mapping[str, np.ndarray], + *, + affine: np.ndarray, + ordered_anatomical_names: Sequence[str], + ct: Optional[np.ndarray] = None, + config: Optional[VertebralInstanceAnalysisConfig] = None, +) -> VertebralAnalysisReport: + """Return a deterministic read-only vertebral-instance audit. + + The returned object contains JSON-compatible primitives only. It never + contains segmentation arrays, voxel-index lists, or proposed corrections. + """ + + effective = config or VertebralInstanceAnalysisConfig() + _validate_config(effective) + names = tuple(str(name) for name in ordered_anatomical_names) + if not names: + raise ValueError("ordered_anatomical_names must not be empty") + if len(set(names)) != len(names): + raise ValueError("ordered_anatomical_names must contain unique names") + + ( + affine_array, + spacing, + orientation, + si_axis, + si_sign, + unsupported_reason, + ) = _validate_affine(affine, effective) + shape, arrays = _shape_and_masks(segmentation_dict, names) + if ct is not None: + ct_array = np.asarray(ct) + if ct_array.ndim != 3: + raise ValueError("ct must be a three-dimensional array") + if shape and list(ct_array.shape) != shape: + raise ValueError("ct shape differs from vertebral masks") + if not np.all(np.isfinite(ct_array)): + raise ValueError("ct contains non-finite values") + else: + ct_array = None + + if unsupported_reason is not None: + return _empty_report( + shape=shape, + spacing=spacing, + orientation=orientation, + ct_used=ct_array is not None, + config=effective, + anomaly=_anomaly( + "unsupported_affine", + [], + [], + "unresolved", + unsupported_reason, + ), + ) + if not shape: + return _empty_report( + shape=[], + spacing=spacing, + orientation=orientation, + ct_used=ct_array is not None, + config=effective, + ) + + occupancy = np.zeros(tuple(shape), dtype=np.uint16) + binary_masks: Dict[str, np.ndarray] = {} + for name in names: + array = arrays.get(name) + if array is None: + continue + mask = np.asarray(array != 0) + binary_masks[name] = mask + occupancy += mask + union = occupancy > 0 + overlap = occupancy > 1 + overlap_count = int(np.count_nonzero(overlap)) + if not np.any(union): + report = _empty_report( + shape=shape, + spacing=spacing, + orientation=orientation, + ct_used=ct_array is not None, + config=effective, + ) + report["input_overlap_voxel_count"] = overlap_count + return report + + trajectory = _fit_trajectory( + binary_masks, + names, + affine_array, + effective.trajectory_outlier_mm, + ) + if trajectory is None: + return _empty_report( + shape=shape, + spacing=spacing, + orientation=orientation, + ct_used=ct_array is not None, + config=effective, + anomaly=_anomaly( + "ambiguous_identity", + [], + [], + "unresolved", + "No trajectory seeds could be estimated from the named masks.", + ), + ) + + crop = _bbox(union) + starts = np.asarray([axis_slice.start or 0 for axis_slice in crop], dtype=int) + union_crop = np.asarray(union[crop]) + distance = ndimage.distance_transform_edt(union_crop, sampling=spacing) + thick = distance >= effective.core_radius_mm + profile, peaks = _profile(thick, spacing, si_axis, effective) + structure = ndimage.generate_binary_structure(3, 1) + component_map, component_count = ndimage.label(thick, structure=structure) + component_slices = ndimage.find_objects( + component_map, + max_label=component_count, + ) + voxel_volume = float(abs(np.linalg.det(affine_array[:3, :3]))) + name_order = {name: index for index, name in enumerate(names)} + candidates: List[_Candidate] = [] + rejected_candidates: List[_RejectedCandidate] = [] + + for component_id, component_slice in enumerate( + component_slices, + start=1, + ): + if component_slice is None: + continue + component_starts = np.asarray( + [axis_slice.start or 0 for axis_slice in component_slice], + dtype=int, + ) + component_region = component_map[component_slice] + component_local_coords = np.argwhere( + component_region == component_id + ) + local_coords = component_local_coords + component_starts + if not len(local_coords): + continue + global_coords = local_coords + starts + world = _apply_affine(affine_array, global_coords) + centroid_world = np.mean(world, axis=0) + trajectory_point = trajectory.point_at_world_z(float(centroid_world[2])) + trajectory_distance = float( + np.linalg.norm(centroid_world[:2] - trajectory_point[:2]) + ) + count = int(len(global_coords)) + volume = count * voxel_volume + if trajectory_distance > effective.trajectory_tube_radius_mm: + rejected_candidates.append( + _RejectedCandidate( + centroid_world=centroid_world, + core_voxel_count=count, + core_volume_mm3=volume, + trajectory_distance_mm=trajectory_distance, + reasons=["core_outside_trajectory_tube"], + ) + ) + continue + + component_distance = distance[tuple(local_coords.T)] + max_radius = float(np.max(component_distance)) + # Include one voxel's physical SI extent. A coordinate range measures + # center-to-center distance and otherwise underestimates persistence + # more severely on coarse/anisotropic grids. + persistence = float( + np.ptp(world[:, 2]) + + np.sum(np.abs(affine_array[2, :3])) + ) + compactness = _compactness(world) + core_overlap = bool(np.any(overlap[tuple(global_coords.T)])) + composition: List[NameCount] = [] + counts: List[Tuple[str, int]] = [] + for name in names: + mask = binary_masks.get(name) + if mask is None: + continue + name_count = int(np.count_nonzero(mask[tuple(global_coords.T)])) + if name_count: + counts.append((name, name_count)) + composition.append( + { + "anatomical_name": name, + "voxel_count": name_count, + "fraction_of_core": _rounded(name_count / count), + } + ) + counts.sort(key=lambda item: (-item[1], name_order[item[0]], item[0])) + composition.sort( + key=lambda item: ( + name_order[item["anatomical_name"]], + item["anatomical_name"], + ) + ) + winner = counts[0][0] if counts else None + winner_count = counts[0][1] if counts else 0 + runner_count = counts[1][1] if len(counts) > 1 else 0 + vote_fraction = winner_count / count if count else 0.0 + vote_margin = (winner_count - runner_count) / count if count else 0.0 + + bone_support: Optional[float] + if ct_array is None: + bone_support = None + else: + ct_values = ct_array[tuple(global_coords.T)] + bone_support = float( + np.count_nonzero(ct_values >= effective.bone_hu_threshold) / count + ) + confidence = _confidence( + max_radius_mm=max_radius, + volume_mm3=volume, + persistence_mm=persistence, + compactness=compactness, + trajectory_distance_mm=trajectory_distance, + vote_fraction=vote_fraction, + vote_margin=vote_margin, + bone_support=bone_support, + config=effective, + ) + reasons: List[str] = [] + if volume < effective.min_core_volume_mm3: + reasons.append("core_volume_below_threshold") + if persistence < effective.min_core_persistence_mm: + reasons.append("core_persistence_below_threshold") + if persistence > effective.max_core_persistence_mm: + reasons.append("possible_merged_instance") + if compactness < effective.min_compactness: + reasons.append("core_compactness_below_threshold") + if trajectory_distance > effective.max_trajectory_distance_mm: + reasons.append("core_far_from_trajectory") + if winner is None: + reasons.append("no_anatomical_name_vote") + if vote_fraction < effective.min_label_vote_fraction: + reasons.append("mixed_identity_vote_fraction") + if vote_margin < effective.min_label_vote_margin: + reasons.append("mixed_identity_vote_margin") + if bone_support is not None and bone_support < effective.min_bone_support_fraction: + reasons.append("low_ct_bone_support") + local_si_min = int(np.min(local_coords[:, si_axis])) + local_si_max = int(np.max(local_coords[:, si_axis])) + contained_peaks = peaks[ + (peaks >= local_si_min) & (peaks <= local_si_max) + ] + if len(contained_peaks) > 1: + reasons.append("multiple_body_profile_peaks_in_one_core") + + threshold = ( + effective.geometry_confidence_threshold + if ct_array is None + else effective.ct_augmented_confidence_threshold + ) + if core_overlap: + status = "unresolved_overlap" + reasons.append("overlapping_input_masks") + elif ( + "mixed_identity_vote_fraction" in reasons + or "mixed_identity_vote_margin" in reasons + or "no_anatomical_name_vote" in reasons + ): + status = "unresolved_mixed_identity" + elif reasons or confidence < threshold: + status = "unresolved_low_confidence" + if confidence < threshold: + reasons.append("confidence_below_threshold") + else: + status = "protected_high_confidence" + + candidates.append( + _Candidate( + centroid_world=centroid_world, + core_voxel_count=count, + core_volume_mm3=volume, + maximum_internal_thickness_mm=2.0 * max_radius, + persistence_mm=persistence, + compactness=compactness, + trajectory_distance_mm=trajectory_distance, + composition=composition, + proposed_identity=winner if status == "protected_high_confidence" else None, + confidence=confidence, + ct_bone_support_fraction=bone_support, + status=status, + reasons=sorted(set(reasons)), + min_world_z=float(np.min(world[:, 2])), + max_world_z=float(np.max(world[:, 2])), + center_si_index=float(np.mean(local_coords[:, si_axis])), + ) + ) + + candidates.sort( + key=lambda item: ( + float(item.centroid_world[2]), + float(item.centroid_world[1]), + float(item.centroid_world[0]), + ) + ) + inferior_touch, superior_touch = _si_boundary_touches( + union, si_axis, si_sign + ) + boundary_candidates: List[Tuple[_Candidate, str]] = [] + if candidates and inferior_touch: + boundary_candidates.append((candidates[0], "inferior")) + if candidates and superior_touch: + boundary_candidates.append((candidates[-1], "superior")) + for candidate, boundary_name in boundary_candidates: + candidate.status = "unresolved_boundary_truncated" + candidate.proposed_identity = None + candidate.reasons = sorted( + set( + candidate.reasons + + [f"physical_instance_touches_{boundary_name}_scan_boundary"] + ) + ) + + for index in range(len(candidates) - 1): + inferior = candidates[index] + superior = candidates[index + 1] + core_separation = max( + 0.0, + superior.min_world_z - inferior.max_world_z, + ) + inferior.superior_core_separation_mm = core_separation + superior.inferior_core_separation_mm = core_separation + quality = _core_valley_quality( + profile, inferior.center_si_index, superior.center_si_index + ) + if quality < effective.min_instance_peak_prominence: + for candidate in (inferior, superior): + if candidate.status == "protected_high_confidence": + candidate.status = "unresolved_low_confidence" + candidate.proposed_identity = None + candidate.reasons = sorted( + set( + candidate.reasons + + ["weak_core_separation_evidence"] + ) + ) + + instance_reports: List[VertebralInstanceReport] = [] + for rank, candidate in enumerate(candidates, start=1): + instance_reports.append( + { + "instance_id": f"instance_{rank:03d}", + "rank_inferior_to_superior": rank, + "centroid_world_mm": _rounded_vector(candidate.centroid_world), + "core_voxel_count": candidate.core_voxel_count, + "core_volume_mm3": _rounded(candidate.core_volume_mm3), + "maximum_internal_thickness_mm": _rounded( + candidate.maximum_internal_thickness_mm + ), + "persistence_mm": _rounded(candidate.persistence_mm), + "compactness": _rounded(candidate.compactness), + "trajectory_distance_mm": _rounded( + candidate.trajectory_distance_mm + ), + "current_name_composition": candidate.composition, + "protected_core_identity": candidate.proposed_identity, + "identity_confidence": _rounded(candidate.confidence), + "confidence_mode": ( + "geometry_only" if ct_array is None else "geometry_and_ct" + ), + "ct_bone_support_fraction": ( + None + if candidate.ct_bone_support_fraction is None + else _rounded(candidate.ct_bone_support_fraction) + ), + "inferior_core_separation_mm": ( + None + if candidate.inferior_core_separation_mm is None + else _rounded(candidate.inferior_core_separation_mm) + ), + "superior_core_separation_mm": ( + None + if candidate.superior_core_separation_mm is None + else _rounded(candidate.superior_core_separation_mm) + ), + "status": candidate.status, # type: ignore[typeddict-item] + "reasons": candidate.reasons, + } + ) + + rejected_candidates.sort( + key=lambda item: ( + float(item.centroid_world[2]), + float(item.centroid_world[1]), + float(item.centroid_world[0]), + item.core_voxel_count, + ) + ) + rejected_candidate_reports: List[RejectedCandidateReport] = [] + for rank, candidate in enumerate(rejected_candidates, start=1): + rejected_candidate_reports.append( + { + "candidate_id": f"rejected_candidate_{rank:03d}", + "centroid_world_mm": _rounded_vector( + candidate.centroid_world + ), + "core_voxel_count": candidate.core_voxel_count, + "core_volume_mm3": _rounded(candidate.core_volume_mm3), + "trajectory_distance_mm": _rounded( + candidate.trajectory_distance_mm + ), + "status": "rejected_unresolved", + "reasons": candidate.reasons, + } + ) + + identities = [item["protected_core_identity"] for item in instance_reports] + fov_status = _field_of_view_status( + union, si_axis, si_sign, identities, names + ) + anomalies: List[SequenceAnomalyReport] = [] + if rejected_candidate_reports: + anomalies.append( + _anomaly( + "off_trajectory_core", + [], + [], + "unresolved", + "One or more thick-core candidates lie outside the " + "configured trajectory tube and were not accepted as " + "vertebral instances.", + candidate_ids=[ + item["candidate_id"] + for item in rejected_candidate_reports + ], + ) + ) + if overlap_count: + affected = [ + item["instance_id"] + for item in instance_reports + if item["status"] == "unresolved_overlap" + ] + anomalies.append( + _anomaly( + "overlapping_input_masks", + affected, + [], + "unresolved", + f"{overlap_count} voxels belong to more than one anatomical-name mask.", + ) + ) + + identity_to_instances: Dict[str, List[str]] = {} + for item in instance_reports: + identity = item["protected_core_identity"] + if identity is not None: + identity_to_instances.setdefault(identity, []).append(item["instance_id"]) + for identity in names: + matching = identity_to_instances.get(identity, []) + if len(matching) > 1: + anomalies.append( + _anomaly( + "duplicate_identity", + matching, + [identity], + "detected", + f"Multiple high-confidence physical cores vote for {identity}.", + ) + ) + + for left, right in zip(instance_reports, instance_reports[1:]): + left_name = left["protected_core_identity"] + right_name = right["protected_core_identity"] + if left_name is None or right_name is None: + continue + left_index = name_order[left_name] + right_index = name_order[right_name] + if right_index < left_index: + anomalies.append( + _anomaly( + "nonmonotonic_identity", + [left["instance_id"], right["instance_id"]], + [left_name, right_name], + "detected", + "Confident identities decrease while physical instances move superiorly.", + ) + ) + elif right_index - left_index > 1: + missing = list(names[left_index + 1 : right_index]) + anomalies.append( + _anomaly( + "missing_internal_identity", + [left["instance_id"], right["instance_id"]], + missing, + "detected", + "Confident neighboring physical cores skip internal anatomical names.", + ) + ) + + sequence_codes = {item["anomaly_code"] for item in anomalies} + if ( + "duplicate_identity" in sequence_codes + and "missing_internal_identity" in sequence_codes + ): + anomalies.append( + _anomaly( + "ambiguous_identity", + [item["instance_id"] for item in instance_reports], + [ + name + for name in identities + if name is not None + ], + "unresolved", + "Combined duplicate and missing identities may reflect " + "transitional anatomy or another nonstandard sequence.", + ) + ) + + unresolved_instances = [ + item for item in instance_reports if item["status"] != "protected_high_confidence" + ] + for item in unresolved_instances: + anomalies.append( + _anomaly( + "ambiguous_identity", + [item["instance_id"]], + [ + entry["anatomical_name"] + for entry in item["current_name_composition"] + ], + "unresolved", + "Physical core identity remains unresolved: " + + ", ".join(item["reasons"]), + ) + ) + if union.any() and not instance_reports: + anomalies.append( + _anomaly( + "ambiguous_identity", + [], + [], + "unresolved", + "Vertebral foreground exists but no thick core passed candidate selection.", + ) + ) + + if len(instance_reports) >= 5: + centroids = np.asarray( + [item["centroid_world_mm"] for item in instance_reports], dtype=np.float64 + ) + spacings = np.linalg.norm(np.diff(centroids, axis=0), axis=1) + median = float(np.median(spacings)) + mad = float(np.median(np.abs(spacings - median))) + scale = max(1.4826 * mad, 1.0) + for index, value in enumerate(spacings): + if abs(float(value) - median) > effective.spacing_outlier_mad * scale: + left = instance_reports[index] + right = instance_reports[index + 1] + anomalies.append( + _anomaly( + "abnormal_spacing", + [left["instance_id"], right["instance_id"]], + [ + name + for name in ( + left["protected_core_identity"], + right["protected_core_identity"], + ) + if name is not None + ], + "unresolved", + f"Inter-core spacing {_rounded(value)} mm is a robust local outlier.", + ) + ) + + anomaly_order = { + "abnormal_spacing": 0, + "ambiguous_identity": 1, + "duplicate_identity": 2, + "missing_internal_identity": 3, + "nonmonotonic_identity": 4, + "off_trajectory_core": 5, + "overlapping_input_masks": 6, + "unsupported_affine": 7, + } + anomalies.sort( + key=lambda item: ( + anomaly_order[item["anomaly_code"]], + item["affected_instance_ids"], + item["affected_anatomical_names"], + ) + ) + if any(item["status"] == "unresolved" for item in anomalies): + overall_status = "unresolved" + elif any(item["status"] == "detected" for item in anomalies): + overall_status = "anomaly_detected" + elif fov_status != "not_truncated_at_array_boundary": + overall_status = "unresolved" + else: + overall_status = "continuous_sequence" + + return { + "schema_version": _SCHEMA_VERSION, + "ordered_names_direction": "inferior_to_superior", + "shape": shape, + "spacing_mm": _rounded_vector(spacing), + "orientation_axcodes": orientation, + "ct_evidence": "used" if ct_array is not None else "unavailable", + "field_of_view_status": fov_status, # type: ignore[typeddict-item] + "effective_config": { + key: _rounded(value) for key, value in asdict(effective).items() + }, + "input_overlap_voxel_count": overlap_count, + "instances": instance_reports, + "rejected_candidates": rejected_candidate_reports, + "observed_sequence_inferior_to_superior": identities, + "anomalies": anomalies, + "overall_status": overall_status, # type: ignore[typeddict-item] + } + + +__all__ = [ + "VertebralInstanceAnalysisConfig", + "VertebralAnalysisReport", + "VertebralInstanceReport", + "RejectedCandidateReport", + "SequenceAnomalyReport", + "analyze_vertebral_instances", +] diff --git a/utils/vertebrae_instance_audit.py b/utils/vertebrae_instance_audit.py new file mode 100644 index 0000000..11eb6e9 --- /dev/null +++ b/utils/vertebrae_instance_audit.py @@ -0,0 +1,431 @@ +"""Opt-in batch adapter for read-only vertebral instance audit reports.""" + +from __future__ import annotations + +import json +import os +import re +import tempfile +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable, Literal, Mapping, Sequence, Tuple + +import numpy as np + + +@dataclass(frozen=True) +class VertebraeInstanceAuditConfig: + """Validated batch-integration configuration.""" + + enabled: bool = False + output_dir_name: str = "vertebrae_analysis" + use_reference_as_ct: bool = False + + +AuditRunStatus = Literal["disabled", "written", "failed"] + +_CONFIG_KEYS = { + "enabled", + "output_dir_name", + "use_reference_as_ct", +} +_VERTEBRA_NAME = re.compile(r"^vertebrae_([LTC])([1-9][0-9]*)$") +_REGION_ORDER = {"L": 0, "T": 1, "C": 2} +_REGION_MAXIMUM = {"L": 5, "T": 12, "C": 7} +_REPORT_SUFFIX = ".json" + + +def parse_vertebrae_instance_audit_config( + raw_config: object, +) -> VertebraeInstanceAuditConfig: + """Validate the optional nested configuration block.""" + + if raw_config is None: + return VertebraeInstanceAuditConfig() + if not isinstance(raw_config, Mapping): + raise ValueError( + "vertebrae_instance_analysis must be a mapping" + ) + + unknown = sorted( + (key for key in raw_config if key not in _CONFIG_KEYS), + key=str, + ) + if unknown: + raise ValueError( + "Unknown vertebrae_instance_analysis configuration keys: " + + ", ".join(str(key) for key in unknown) + ) + + enabled = raw_config.get("enabled", False) + use_reference_as_ct = raw_config.get( + "use_reference_as_ct", + False, + ) + if type(enabled) is not bool: + raise ValueError( + "vertebrae_instance_analysis.enabled must be a boolean" + ) + if type(use_reference_as_ct) is not bool: + raise ValueError( + "vertebrae_instance_analysis.use_reference_as_ct " + "must be a boolean" + ) + + output_dir_name = raw_config.get( + "output_dir_name", + "vertebrae_analysis", + ) + if not isinstance(output_dir_name, str): + raise ValueError( + "vertebrae_instance_analysis.output_dir_name " + "must be a string" + ) + if ( + not output_dir_name + or output_dir_name != output_dir_name.strip() + or output_dir_name in {".", ".."} + or os.path.isabs(output_dir_name) + or "/" in output_dir_name + or "\\" in output_dir_name + ): + raise ValueError( + "vertebrae_instance_analysis.output_dir_name must be " + "one nonempty relative path component" + ) + + return VertebraeInstanceAuditConfig( + enabled=enabled, + output_dir_name=output_dir_name, + use_reference_as_ct=use_reference_as_ct, + ) + + +def ordered_vertebral_anatomical_names( + anatomical_names: Iterable[str], +) -> Tuple[str, ...]: + """Return canonical vertebral names in inferior-to-superior order.""" + + parsed = [] + seen = set() + for anatomical_name in anatomical_names: + if not isinstance(anatomical_name, str): + raise ValueError("class-map anatomical names must be strings") + if not anatomical_name.startswith("vertebrae_"): + continue + match = _VERTEBRA_NAME.fullmatch(anatomical_name) + if match is None: + raise ValueError( + "Unsupported vertebral anatomical name: " + f"{anatomical_name}" + ) + region, level_text = match.groups() + level = int(level_text) + if level > _REGION_MAXIMUM[region]: + raise ValueError( + "Unsupported vertebral anatomical name: " + f"{anatomical_name}" + ) + if anatomical_name in seen: + raise ValueError( + "Duplicate vertebral anatomical name: " + f"{anatomical_name}" + ) + seen.add(anatomical_name) + parsed.append((region, level, anatomical_name)) + + if not parsed: + raise ValueError( + "No canonical vertebral anatomical names were configured" + ) + + parsed.sort( + key=lambda item: ( + _REGION_ORDER[item[0]], + -item[1], + item[2], + ) + ) + return tuple(item[2] for item in parsed) + + +def audit_output_path( + output_root: os.PathLike, + config: VertebraeInstanceAuditConfig, + patient_id: str, +) -> Path: + """Build the deterministic report path outside case segmentations.""" + + if ( + not isinstance(patient_id, str) + or not patient_id + or patient_id in {".", ".."} + or "/" in patient_id + or "\\" in patient_id + or "\x00" in patient_id + or any(not character.isprintable() for character in patient_id) + ): + raise ValueError( + "patient_id must be one printable nonempty path component" + ) + return ( + Path(output_root) + / config.output_dir_name + / f"{patient_id}{_REPORT_SUFFIX}" + ) + + +def _prepare_audit_output_path( + output_root: os.PathLike, + config: VertebraeInstanceAuditConfig, + patient_id: str, +) -> Tuple[Path, Path]: + """Create and validate a real audit directory below output_root.""" + + lexical_path = audit_output_path(output_root, config, patient_id) + resolved_output_root = Path(output_root).resolve(strict=False) + audit_directory = ( + resolved_output_root / config.output_dir_name + ) + if audit_directory.is_symlink(): + raise ValueError("audit output directory must not be a symlink") + + audit_directory.mkdir(parents=True, exist_ok=True) + if audit_directory.is_symlink(): + raise ValueError("audit output directory must not be a symlink") + + resolved_audit_directory = audit_directory.resolve(strict=True) + try: + resolved_audit_directory.relative_to(resolved_output_root) + except ValueError as error: + raise ValueError( + "audit output directory escapes the resolved output root" + ) from error + + report_path = ( + resolved_audit_directory / lexical_path.name + ) + return report_path, resolved_output_root + + +def _load_analyzer(): + # Keep the analyzer and its SciPy imports out of disabled batch runs. + from .vertebrae_instance_analysis import analyze_vertebral_instances + + return analyze_vertebral_instances + + +def _reference_ct_or_none( + *, + reference_img, + segmentation_dict: Mapping[str, np.ndarray], + ordered_names: Sequence[str], + logger, + patient_id: str, +): + requested_masks = [ + np.asarray(segmentation_dict[name]) + for name in ordered_names + if name in segmentation_dict + ] + if not requested_masks: + logger.warning( + "[ShapeKit][vertebrae-instance-analysis] " + "CT requested for %s, but no vertebral mask is available " + "for alignment validation; using geometry-only mode.", + patient_id, + ) + return None + + try: + ct_array = np.asanyarray(reference_img.dataobj) + except Exception: + logger.warning( + "[ShapeKit][vertebrae-instance-analysis] " + "Failed to load the requested CT reference for %s; " + "using geometry-only mode.", + patient_id, + exc_info=True, + ) + return None + + expected_shape = requested_masks[0].shape + try: + valid_masks = all( + mask.ndim == 3 and mask.shape == expected_shape + for mask in requested_masks + ) + valid_ct = ( + valid_masks + and ct_array.ndim == 3 + and ct_array.shape == expected_shape + and bool(np.all(np.isfinite(ct_array))) + ) + except Exception: + valid_ct = False + if not valid_ct: + logger.warning( + "[ShapeKit][vertebrae-instance-analysis] " + "The explicitly requested CT reference for %s does not " + "satisfy the 3-D shape/finiteness contract; using " + "geometry-only mode.", + patient_id, + ) + return None + return ct_array + + +def _canonical_json_bytes(report: Mapping[str, object]) -> bytes: + return ( + json.dumps( + report, + sort_keys=True, + separators=(",", ":"), + ensure_ascii=True, + allow_nan=False, + ) + + "\n" + ).encode("utf-8") + + +def _atomic_write( + path: Path, + payload: bytes, + *, + resolved_output_root: Path, +) -> None: + if path.parent.is_symlink(): + raise ValueError("audit output directory must not be a symlink") + resolved_parent = path.parent.resolve(strict=True) + try: + resolved_parent.relative_to(resolved_output_root) + except ValueError as error: + raise ValueError( + "audit output directory escapes the resolved output root" + ) from error + if resolved_parent != path.parent: + raise ValueError("audit output directory must not be a symlink") + + descriptor, temporary_name = tempfile.mkstemp( + prefix=f".{path.stem}.", + suffix=".tmp", + dir=str(resolved_parent), + ) + temporary_path = Path(temporary_name) + try: + with os.fdopen(descriptor, "wb") as stream: + stream.write(payload) + stream.flush() + os.fsync(stream.fileno()) + os.replace(str(temporary_path), str(path)) + finally: + try: + temporary_path.unlink() + except FileNotFoundError: + pass + + +def run_vertebrae_instance_audit( + segmentation_dict: Mapping[str, np.ndarray], + *, + reference_img, + ordered_anatomical_names: Sequence[str], + output_root: os.PathLike, + patient_id: str, + config: VertebraeInstanceAuditConfig, + logger, +) -> AuditRunStatus: + """Run and save an optional audit without affecting segmentation.""" + + if not config.enabled: + return "disabled" + + requested_ct_mode = ( + "reference-as-ct" + if config.use_reference_as_ct + else "geometry-only" + ) + effective_ct_mode = "not-evaluated" + report_path = None + try: + report_path = audit_output_path( + output_root, + config, + patient_id, + ) + report_path, resolved_output_root = _prepare_audit_output_path( + output_root, + config, + patient_id, + ) + except Exception: + logger.exception( + "[ShapeKit][vertebrae-instance-analysis] " + "Invalid audit output path for patient=%s, report=%s, " + "requested_ct_mode=%s, effective_ct_mode=%s; segmentation " + "postprocessing will continue unchanged.", + patient_id, + report_path, + requested_ct_mode, + effective_ct_mode, + ) + return "failed" + + try: + ct_array = None + if config.use_reference_as_ct: + ct_array = _reference_ct_or_none( + reference_img=reference_img, + segmentation_dict=segmentation_dict, + ordered_names=ordered_anatomical_names, + logger=logger, + patient_id=patient_id, + ) + effective_ct_mode = ( + "CT-supported" + if ct_array is not None + else "geometry-only" + ) + analyze_vertebral_instances = _load_analyzer() + report = analyze_vertebral_instances( + segmentation_dict, + affine=reference_img.affine, + ordered_anatomical_names=ordered_anatomical_names, + ct=ct_array, + ) + payload = _canonical_json_bytes(report) + _atomic_write( + report_path, + payload, + resolved_output_root=resolved_output_root, + ) + logger.info( + "[ShapeKit][vertebrae-instance-analysis] " + "Wrote %s audit for %s to %s.", + effective_ct_mode, + patient_id, + report_path, + ) + return "written" + except Exception: + logger.exception( + "[ShapeKit][vertebrae-instance-analysis] " + "Audit failed for patient=%s, report=%s, " + "requested_ct_mode=%s, effective_ct_mode=%s; segmentation " + "postprocessing will continue unchanged.", + patient_id, + report_path, + requested_ct_mode, + effective_ct_mode, + ) + return "failed" + + +__all__ = [ + "AuditRunStatus", + "VertebraeInstanceAuditConfig", + "audit_output_path", + "ordered_vertebral_anatomical_names", + "parse_vertebrae_instance_audit_config", + "run_vertebrae_instance_audit", +]