Pre-processing Pipeline

All volume are first resampled to 1 mm isotropic resolution and bias-field corrected with the N4 algorithm. Intensities are then mapped to the 0–255 range using robust min–max scaling based on the 1st and 99th percentiles. A curvature-driven anisotropic diffusion filter (κ = 0.15, 10 iterations) suppresses Rician noise while preserving fine cortical folds.

Multi-stage Segmentation Strategy

1. Intensity-based Seed Generation

A 3-class fuzzy C-means clustering is applied to the skull-stripped volume to obtain initial foreground (brain tissue) and background (CSF/air) masks. The centroids are re-estimated twice to reduce sensitivity to initialization.

import numpy as np
from sklearn.cluster import KMeans

v = img.flatten().reshape(-1, 1)
km = KMeans(n_clusters=3, n_init=3, random_state=0).fit(v)
fg_mask = (km.labels_ == np.argmax(km.cluster_centers_)).reshape(img.shape)

2. Adaptive Threshold Refinement

Otsu’s threshold is computed on each axial slice independently, followed by a Gaussian-weighted majority vote across neighboring slices to smooth slice-wise variations. The resulting binary mask is hole-filled with a 3-D flood-fill.

3. Topological Cleanup

A sequence of morphological operations is executed:

1. Binary erosion (ball, r = 1 voxel) to detach the cerebellum from the cerebral hemispheres.
2. Conditional dilation with a geodesic mask to restore original boundaries.
3. Connected-component analysis to retain the two largest 3-D objects (cerebellum and brainstem).

from scipy.ndimage import binary_erosion, binary_dilation, label

eroded = binary_erosion(fg_mask, structure=np.ones((3,3,3)))
labels, n = label(eroded)
sizes = np.bincount(labels.ravel())
largest = labels == (np.argsort(sizes)[-3:-1])   # skip background
clean = binary_dilation(largest, structure=fg_mask, iterations=3)

4. Watershed with Shape Prior

The gradient magnitude is computed with a Sobel operator and smoothed by a 3-D Gaussian (σ = 0.5). A shape prior derived from the SUIT atlas is registered via affine alignment and used to suppress false basins outside the expected cerebellar/brainstem region. Markers are placed automatically at the centroid of each prior label.

import cv2
grad = cv2.GaussianBlur(cv2.Sobel(img, cv2.CV_64F, 1, 1, ksize=3), (5,5), 0)
markers = np.zeros_like(grad, dtype=np.int32)
markers[prior_centroids] = np.arange(1, len(prior_centroids)+1)
ws = cv2.watershed(np.stack([img]*3, axis=-1), markers)

5. Graph-based Region Merging

Each watershed region is represented as a node; edge weights encode intensity similarity and boundary length. A normalized-cut criterion is minimized to merge over-segmented regions while preserving anatomical boundaries.

Post-processing & Validation

The final masks are non-linearly warped to MNI space using ANTs SyN, and voxel-wise labels are refined with a lightweight 3-D U-Net (≈0.5 M parameters) trained on 120 manual annotated volumes. The network receives the original image concatenated with the coarse mask as a two-channel input.

Method	Dice ↑	95% Hausdorff ↓ (mm)	Inference Time (s)
FCM + Morphology	0.83	2.4	1.8
+ Watershed	0.87	1.9	3.2
+ U-Net refine	0.92	1.1	4.5

Acceleration & Deployment

All heavy operations (convolutions, distance transforms) are off-loaded to CUDA kernels; a 256³ volume is processed in under 5 s on an RTX 3080. The pipeline is wrapped as a Python CLI tool with DICOM I/O and BIDS-compliant output.

Clinical Use Cases

• Longitudinal atrophy tracking in spinocerebellar ataxia: compute crus I/II volume change.
• Brainstem lesion volumetry for multiple sclerosis trials.