feat(mri): add mask_brain (intensity threshold + morphological opening)
Browse files
src/pipelines/mri_pipeline.py
CHANGED
|
@@ -13,8 +13,10 @@ from __future__ import annotations
|
|
| 13 |
|
| 14 |
import os
|
| 15 |
|
|
|
|
| 16 |
import numpy as np
|
| 17 |
import pyarrow as pa
|
|
|
|
| 18 |
|
| 19 |
from src.core.logger import get_logger
|
| 20 |
|
|
@@ -48,3 +50,32 @@ def is_valid_volume(volume: np.ndarray | None) -> bool:
|
|
| 48 |
if not np.all(np.isfinite(volume)):
|
| 49 |
return False
|
| 50 |
return True
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 13 |
|
| 14 |
import os
|
| 15 |
|
| 16 |
+
import nibabel as nib
|
| 17 |
import numpy as np
|
| 18 |
import pyarrow as pa
|
| 19 |
+
from scipy import ndimage as scipy_ndimage
|
| 20 |
|
| 21 |
from src.core.logger import get_logger
|
| 22 |
|
|
|
|
| 50 |
if not np.all(np.isfinite(volume)):
|
| 51 |
return False
|
| 52 |
return True
|
| 53 |
+
|
| 54 |
+
|
| 55 |
+
def mask_brain(
|
| 56 |
+
volume: np.ndarray,
|
| 57 |
+
intensity_threshold: float | None = None,
|
| 58 |
+
) -> np.ndarray:
|
| 59 |
+
"""Build a brain mask from a 3-D MRI volume.
|
| 60 |
+
|
| 61 |
+
Two-step pipeline:
|
| 62 |
+
1. Intensity threshold: keep voxels above `intensity_threshold`. When
|
| 63 |
+
`None`, use the volume's mean as a robust auto-threshold (works on
|
| 64 |
+
the synthetic fixture where brain ≫ background; for real data the
|
| 65 |
+
caller should pass an Otsu or BET-derived threshold explicitly).
|
| 66 |
+
2. Morphological opening (`scipy.ndimage.binary_opening`) to remove
|
| 67 |
+
isolated noise voxels and disconnected fragments.
|
| 68 |
+
|
| 69 |
+
Args:
|
| 70 |
+
volume: 3-D numeric `np.ndarray` (must satisfy `is_valid_volume`).
|
| 71 |
+
intensity_threshold: Voxel-intensity floor. `None` → use `volume.mean()`.
|
| 72 |
+
|
| 73 |
+
Returns:
|
| 74 |
+
A boolean `np.ndarray` of the same shape as `volume`. True = brain.
|
| 75 |
+
"""
|
| 76 |
+
if intensity_threshold is None:
|
| 77 |
+
intensity_threshold = float(volume.mean())
|
| 78 |
+
|
| 79 |
+
raw = volume > intensity_threshold
|
| 80 |
+
cleaned = scipy_ndimage.binary_opening(raw, iterations=1)
|
| 81 |
+
return cleaned.astype(bool)
|
tests/pipelines/test_mri_pipeline.py
CHANGED
|
@@ -3,10 +3,14 @@ from __future__ import annotations
|
|
| 3 |
|
| 4 |
from pathlib import Path
|
| 5 |
|
|
|
|
| 6 |
import numpy as np
|
| 7 |
import pytest
|
| 8 |
|
| 9 |
-
from src.pipelines.mri_pipeline import
|
|
|
|
|
|
|
|
|
|
| 10 |
|
| 11 |
|
| 12 |
FIXTURE_DIR = Path(__file__).parent.parent / "fixtures" / "mri_sample"
|
|
@@ -45,3 +49,57 @@ class TestIsValidVolume:
|
|
| 45 |
def test_rejects_non_array(self) -> None:
|
| 46 |
assert is_valid_volume([[[1, 2]], [[3, 4]]]) is False
|
| 47 |
assert is_valid_volume(None) is False
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 3 |
|
| 4 |
from pathlib import Path
|
| 5 |
|
| 6 |
+
import nibabel as nib
|
| 7 |
import numpy as np
|
| 8 |
import pytest
|
| 9 |
|
| 10 |
+
from src.pipelines.mri_pipeline import (
|
| 11 |
+
is_valid_volume,
|
| 12 |
+
mask_brain,
|
| 13 |
+
)
|
| 14 |
|
| 15 |
|
| 16 |
FIXTURE_DIR = Path(__file__).parent.parent / "fixtures" / "mri_sample"
|
|
|
|
| 49 |
def test_rejects_non_array(self) -> None:
|
| 50 |
assert is_valid_volume([[[1, 2]], [[3, 4]]]) is False
|
| 51 |
assert is_valid_volume(None) is False
|
| 52 |
+
|
| 53 |
+
|
| 54 |
+
class TestMaskBrain:
|
| 55 |
+
def _load_subject(self, sid: str) -> np.ndarray:
|
| 56 |
+
return nib.load(FIXTURE_DIR / f"{sid}.nii.gz").get_fdata()
|
| 57 |
+
|
| 58 |
+
def test_returns_bool_mask_of_same_shape(self) -> None:
|
| 59 |
+
vol = self._load_subject("subject_0")
|
| 60 |
+
mask = mask_brain(vol)
|
| 61 |
+
assert isinstance(mask, np.ndarray)
|
| 62 |
+
assert mask.dtype == bool
|
| 63 |
+
assert mask.shape == vol.shape
|
| 64 |
+
|
| 65 |
+
def test_mask_separates_brain_from_background(self) -> None:
|
| 66 |
+
"""Default threshold should keep the spherical-brain center voxels in."""
|
| 67 |
+
vol = self._load_subject("subject_0")
|
| 68 |
+
mask = mask_brain(vol)
|
| 69 |
+
# The fixture's brain region (radius 3 around center) intensity is ~10;
|
| 70 |
+
# background is ~0.1. Some brain voxels MUST survive the mask.
|
| 71 |
+
assert mask.sum() > 0
|
| 72 |
+
# The center voxel (always brain) MUST be in the mask.
|
| 73 |
+
center = tuple(s // 2 for s in vol.shape)
|
| 74 |
+
assert mask[center]
|
| 75 |
+
|
| 76 |
+
def test_mask_drops_low_intensity_background(self) -> None:
|
| 77 |
+
"""Voxels with intensity well below the brain core must be excluded."""
|
| 78 |
+
vol = self._load_subject("subject_0")
|
| 79 |
+
mask = mask_brain(vol, intensity_threshold=5.0)
|
| 80 |
+
# Background voxels (intensity ~0.1) must NOT be in the mask.
|
| 81 |
+
bg_voxel = (0, 0, 0)
|
| 82 |
+
assert mask[bg_voxel] == False # noqa: E712
|
| 83 |
+
|
| 84 |
+
def test_explicit_threshold_overrides_default(self) -> None:
|
| 85 |
+
vol = self._load_subject("subject_0")
|
| 86 |
+
# A very high threshold should produce far fewer mask voxels.
|
| 87 |
+
mask_default = mask_brain(vol)
|
| 88 |
+
mask_strict = mask_brain(vol, intensity_threshold=100.0)
|
| 89 |
+
assert mask_strict.sum() < mask_default.sum()
|
| 90 |
+
|
| 91 |
+
def test_does_not_mutate_input(self) -> None:
|
| 92 |
+
vol = self._load_subject("subject_0")
|
| 93 |
+
original = vol.copy()
|
| 94 |
+
_ = mask_brain(vol)
|
| 95 |
+
np.testing.assert_array_equal(vol, original)
|
| 96 |
+
|
| 97 |
+
def test_morphological_cleanup_removes_isolated_voxels(self) -> None:
|
| 98 |
+
"""A single bright voxel surrounded by background must be removed by the
|
| 99 |
+
opening-style morphological cleanup."""
|
| 100 |
+
vol = np.zeros((8, 8, 8), dtype=np.float64)
|
| 101 |
+
vol[4, 4, 4] = 100.0
|
| 102 |
+
mask = mask_brain(vol, intensity_threshold=50.0)
|
| 103 |
+
# Without cleanup, the single voxel would survive. With morphological
|
| 104 |
+
# opening, it must be removed.
|
| 105 |
+
assert mask.sum() == 0
|