Ebrains Sibra Atlas Tool & DICOM Prepocessor

README.md

@@ -13,6 +13,7 @@ a LangGraph-based multi-agent pipeline for automated classification of neuroimag
 | `CNNClassifier` | VGG16 / DenseNet169 / ResNet101 | Task-specific classification |
 | `SAM3Tool` | SAM3 frozen backbone + linear probe | Lesion segmentation (Dice = 0.836) |
 | `BiomedCLIPTool` | microsoft/BiomedCLIP (ViT-B/16, layer 6) | Zero-shot re-ranking for ambiguous multiclass cases |
+| `SiibraAtlasTool` | EBRAINS Julich-Brain v2.9 | Anatomical region assignment via MNI152 coordinates |
 
 **Pipeline flow** (linear - every node runs for every image):
 
@@ -47,7 +48,8 @@ MultiAgentMedClassifier/
 │   ├── medgemma_agent.py   # MedGemma: triage, bbox diagnosis, report
 │   ├── cnn_tool.py         # CNN classifier (VGG16 / DenseNet / ResNet)
 │   ├── sam3_tool.py        # SAM3 segmentation + linear probe head
-│   └── biomedclip_tool.py  # BiomedCLIP zero-shot / linear probe
+│   ├── biomedclip_tool.py  # BiomedCLIP zero-shot / linear probe
+│   └── sibra_tool.py       # EBRAINS siibra: lesion centroid → Julich-Brain region
 ├── pipeline/
 │   ├── graph.py            # LangGraph StateGraph assembly
 │   ├── nodes.py            # Node factory functions
@@ -64,9 +66,12 @@ MultiAgentMedClassifier/
 ├── prompts/
 │   ├── system_prompt.txt       # MedGemma radiologist persona + JSON schema
 │   └── system_prompt_bbox.txt  # Same schema, bbox-overlay context
+├── tests/
+│   └── test_atlas_enrichment.py  # Standalone test for atlas node on a single image
 ├── checkpoints/            # PyTorch state dicts: {arch}_{dataset}_final.pt
 ├── outputs/
 │   ├── explainability/     # Saliency maps: gradcam_pp_*.png, ig_*.png
+│   ├── segmentation/       # SAM3 binary masks (mask_*.png) and bbox overlays (guided_*.png)
 │   ├── fhir/               # FHIR R4 bundles: fhir_<id>.json
 │   └── eval/               # comparison_summary.csv, <task>_tumor_eval.jsonl
 ├── app.py                  # Gradio web GUI (python app.py → http://localhost:7860)
@@ -210,12 +215,25 @@ zero-shot mode (BiomedCLIP) instead of attempting a download.
 
 ## Usage
 
-**Single image:**
+**Single image (PNG/JPEG):**
 
 ```bash
 python run_pipeline.py --image scan.png --task binary_tumor
 ```
 
+**Single DICOM slice:**
+
+```bash
+python run_pipeline.py --image scan.dcm --task binary_tumor
+```
+
+**DICOM series directory:**
+
+```bash
+python run_pipeline.py --image /path/to/dicom_series --task binary_tumor
+```
+Raw DICOM inputs are converted to `outputs/preprocessed/*.png` before routing. For series inputs, the middle slice is selected as the representative 2D image and the chosen slice path is preserved in `metadata["dicom_path"]` for atlas enrichment.
+
 **With explainability (Grad-CAM++ + Integrated Gradients):**
 
 ```bash
@@ -304,6 +322,37 @@ The report is returned in `state["final_report"]` (plain text, ≤150 words). Th
 | `explainability_result` | Paths to `gradcam_pp_*.png` and `ig_*.png` (if enabled) |
 | `verification_result` | MedGemma post-hoc agreement check against Grad-CAM++ saliency map (if explainability enabled) |
 | `fhir_report` | FHIR R4 DiagnosticReport dict; saved to `outputs/fhir/fhir_<id>.json` |
+| `atlas_enrichment` | EBRAINS atlas assignment: `assigned_region`, `hemisphere`, `mni_coords`, `assignment_scores` (only on `sam3_then_cnn` path) |
+
+## Atlas Enrichment (EBRAINS / siibra)
+
+On the `sam3_then_cnn` path, after SAM3 produces a binary mask, the pipeline runs an optional anatomical assignment step using [siibra-python](https://siibra-python.readthedocs.io) and the [EBRAINS Julich-Brain Cytoarchitectonic Atlas v2.9](https://search.kg.ebrains.eu).
+
+The SAM3 mask centroid is mapped to MNI152 coordinates, which are then assigned to the nearest cytoarchitectonic region via a statistical probability map.
+
+**Coordinate accuracy (best → worst):**
+
+| Input | How it's used |
+|---|---|
+| `metadata["nifti_path"]` | NIfTI affine → exact MNI coords |
+| `metadata["dicom_path"]` | `ImagePositionPatient` + `PixelSpacing` → scanner-space coords (≈ MNI for pre-registered datasets) |
+| PNG only (default) | Pixel centroid normalised to MNI152 range; z locked to axial midplane (0 mm) |
+
+Pass coordinate metadata at inference time:
+```python
+initial_state("scan.png", "binary_tumor", metadata={"nifti_path": "scan.nii.gz"})
+```
+
+If the main input itself is DICOM, `initial_state()` now injects `metadata["dicom_path"]` automatically after preprocessing the slice or series into a PNG.
+
+**Test the node standalone (without running the full pipeline):**
+```bash
+python tests/test_atlas_enrichment.py --image data/scan.png --mask outputs/segmentation/mask_abc123.png
+```
+
+The `atlas_enrichment` result flows into the MedGemma report prompt and is serialised as an `ebrains-atlas-assignment` extension in the FHIR bundle.
+
+**Note:** z=0.0 (axial midplane) is assumed when no NIfTI or DICOM is available. Region assignments are anatomically meaningful but spatially approximate. Scores are Julich-Brain probability map values; scores < 0.1 occur near region boundaries or at the midplane.
 
 ## Explainability Methods
 

agents/dicom_tool.py

@@ -0,0 +1,212 @@
+# agents/dicom_tool.py
+"""
+DICOM ingestion and preprocessing for the neuroimaging pipeline.
+
+Handles:
+  - Single .dcm file (one slice)
+  - DICOM series directory (multi-slice → selects representative slice)
+  - Passthrough for PNG/JPEG (no-op)
+
+Output:
+  - PNG saved to outputs/preprocessed/<stem>_<hash>.png (input to existing pipeline)
+  - Metadata dict injected into NeuroimagingState
+"""
+
+import pydicom
+import numpy as np
+from PIL import Image
+from pathlib import Path
+import hashlib
+
+from pydicom.misc import is_dicom
+
+
+class DICOMPreprocessor:
+
+    def __init__(self, output_dir: str = "outputs/preprocessed"):
+        self.output_dir = Path(output_dir)
+        self.output_dir.mkdir(parents=True, exist_ok=True)
+
+    def prepare(self, input_path: str) -> dict:
+        """
+        Main entry point. Accepts a .dcm file, a directory of .dcm files,
+        or a PNG/JPEG (passthrough).
+
+        Returns:
+            {
+              "image_path":   str,   # path to PNG ready for CNN/MedGemma
+              "nifti_path":   None,  # populated if you add NIfTI conversion (see below)
+              "dicom_path":   str | None,
+              "dicom_metadata": dict | None
+            }
+        """
+        p = Path(input_path)
+
+        if p.is_dir():
+            return self._process_series(p)
+        elif self._is_dicom_file(p):
+            return self._process_single(p)
+        else:
+            # Already a PNG/JPEG — passthrough, no metadata
+            return {
+                "image_path": str(p),
+                "nifti_path": None,
+                "dicom_path": None,
+                "dicom_metadata": None,
+            }
+
+    # ── Single slice ──────────────────────────────────────────────────────────
+
+    def _process_single(self, dcm_path: Path) -> dict:
+        ds = pydicom.dcmread(str(dcm_path))
+        metadata = self._extract_metadata(ds)
+        png_path = self._pixel_array_to_png(ds, self._output_stem(dcm_path))
+        return {
+            "image_path":    str(png_path),
+            "nifti_path":    None,
+            "dicom_path":    str(dcm_path),
+            "dicom_metadata": metadata,
+        }
+
+    # ── Series (folder of slices) ─────────────────────────────────────────────
+
+    def _process_series(self, series_dir: Path) -> dict:
+        """
+        Load a DICOM series, sort slices by InstanceNumber,
+        select the middle slice as the representative 2D image,
+        and extract metadata from that selected slice.
+        """
+        dcm_files = sorted(f for f in series_dir.iterdir() if self._is_dicom_file(f))
+        if not dcm_files:
+            raise ValueError(f"No DICOM files found in {series_dir}")
+
+        # Sort by InstanceNumber (slice position)
+        slices = []
+        for f in dcm_files:
+            ds = pydicom.dcmread(str(f), stop_before_pixels=False)
+            instance = int(getattr(ds, "InstanceNumber", 0))
+            slices.append((instance, f, ds))
+        slices.sort(key=lambda x: x[0])
+
+        # Representative slice: middle of the series
+        # For tumour tasks the middle axial slice is usually most informative
+        # You could also select by maximum lesion area if a rough threshold is known
+        mid_idx  = len(slices) // 2
+        _, mid_path, mid_ds = slices[mid_idx]
+
+        metadata = self._extract_metadata(mid_ds)
+        metadata["total_slices"]      = len(slices)
+        metadata["selected_slice_idx"] = mid_idx
+
+        png_path = self._pixel_array_to_png(mid_ds, self._output_stem(series_dir))
+
+        return {
+            "image_path":    str(png_path),
+            "nifti_path":    None,           # see _series_to_nifti() below if needed
+            "dicom_path":    str(mid_path),
+            "dicom_metadata": metadata,
+        }
+
+    # ── Pixel array → normalised PNG ─────────────────────────────────────────
+
+    def _pixel_array_to_png(self, ds: "pydicom.Dataset", stem: str) -> Path:
+        """
+        Convert DICOM pixel array to a normalised 8-bit grayscale PNG.
+
+        Handles:
+          - Modality LUT (rescale slope/intercept for CT Hounsfield units)
+          - Window/level for soft tissue vs bone CT windows
+          - Inversion for some MRI sequences
+        """
+        pixel_array = ds.pixel_array.astype(np.float32)
+
+        # Apply RescaleSlope / RescaleIntercept (standard for CT)
+        slope     = float(getattr(ds, "RescaleSlope",     1.0))
+        intercept = float(getattr(ds, "RescaleIntercept", 0.0))
+        pixel_array = pixel_array * slope + intercept
+
+        # Apply window/level if present (typical for CT brain window)
+        window_center = getattr(ds, "WindowCenter", None)
+        window_width  = getattr(ds, "WindowWidth",  None)
+        if window_center and window_width:
+            wc = float(window_center[0] if hasattr(window_center, "__iter__") else window_center)
+            ww = float(window_width[0]  if hasattr(window_width,  "__iter__") else window_width)
+            lo = wc - ww / 2
+            hi = wc + ww / 2
+            pixel_array = np.clip(pixel_array, lo, hi)
+
+        if getattr(ds, "PhotometricInterpretation", "") == "MONOCHROME1":
+            pixel_array = pixel_array.max() - pixel_array
+
+        # Min-max normalise to 0–255
+        lo, hi = pixel_array.min(), pixel_array.max()
+        if hi > lo:
+            pixel_array = (pixel_array - lo) / (hi - lo) * 255.0
+        pixel_array = pixel_array.astype(np.uint8)
+
+        # Convert to RGB (your CNNs and MedGemma expect 3-channel input)
+        img = Image.fromarray(pixel_array, mode="L").convert("RGB")
+
+        out_path = self.output_dir / f"{stem}.png"
+        img.save(str(out_path))
+        return out_path
+
+    @staticmethod
+    def _is_dicom_file(path: Path) -> bool:
+        if not path.is_file():
+            return False
+        if path.suffix.lower() in {".dcm", ".dicom", ".ima"}:
+            return True
+        try:
+            return bool(is_dicom(str(path)))
+        except Exception:
+            return False
+
+    @staticmethod
+    def _output_stem(path: Path) -> str:
+        digest = hashlib.sha1(str(path.resolve()).encode("utf-8")).hexdigest()[:8]
+        return f"{path.stem}_{digest}" if path.is_file() else f"{path.name}_{digest}"
+
+    # ── Metadata extraction ───────────────────────────────────────────────────
+
+    @staticmethod
+    def _extract_metadata(ds: "pydicom.Dataset") -> dict:
+        """
+        Pull the clinically and technically relevant tags from the DICOM header.
+        These are injected into NeuroimagingState and passed to MedGemma's prompt.
+        """
+        def safe(tag):
+            val = getattr(ds, tag, None)
+            if val is None:
+                return None
+            # pydicom sequences are not JSON-serialisable — convert to str
+            return str(val) if not isinstance(val, (int, float, str)) else val
+
+        return {
+            # Patient context (anonymised in real deployments)
+            "patient_age":          safe("PatientAge"),       # e.g. "068Y"
+            "patient_sex":          safe("PatientSex"),       # "M" | "F" | "O"
+
+            # Acquisition parameters — critical for modality routing
+            "modality":             safe("Modality"),         # "MR" | "CT"
+            "series_description":   safe("SeriesDescription"),# "T2 FLAIR AX"
+            "sequence_name":        safe("SequenceVariant"),
+            "scanning_sequence":    safe("ScanningSequence"), # "EP" | "GR" | "SE"
+            "field_strength":       safe("MagneticFieldStrength"), # 1.5 | 3.0
+
+            # Scanner provenance
+            "manufacturer":         safe("Manufacturer"),
+            "manufacturer_model":   safe("ManufacturerModelName"),
+            "institution":          safe("InstitutionName"),
+
+            # Geometric parameters (needed for siibra coordinate mapping)
+            "slice_thickness_mm":   safe("SliceThickness"),
+            "pixel_spacing":        safe("PixelSpacing"),     # [row_mm, col_mm]
+            "image_orientation":    safe("ImageOrientationPatient"),
+            "image_position":       safe("ImagePositionPatient"), # MNI registration anchor
+
+            # Study context
+            "study_description":    safe("StudyDescription"),
+            "body_part":            safe("BodyPartExamined"),
+            "accession_number":     safe("AccessionNumber"),  # links to RIS/PACS
+        }