Note
Go to the end to download the full example code.
Pooling voxels across the cohort
Background. Direct pooling skips supervoxels: every ROI voxel of every subject goes into one matrix and one cohort model is fitted on it, so habitat ids are shared across patients.
Purpose. You get one shared habitat model fitted on voxels, a habitat map and volume fractions per subject, and a histogram of the arterial-phase intensities inside each habitat.
When to use. You want shared ids without the supervoxel step, and the pooled voxel matrix is small enough to cluster (it grows with every voxel of every subject). Otherwise the two-step design (Defining habitats in two steps) clusters far fewer rows.
Key terms.
direct pooling –
poolthenfitwith nopartition; each voxel is its own clustering unit.pool – see Defining habitats in two steps.
Input: a cohort of at least two subjects. Output: one shared
HabitatModel fitted on voxels, and one
HabitatMap per subject. The stage list is
pool then fit, with no partition. direct_pooling_habitat(...)
is a shortcut that builds the same stage list.
Load the cohort
Change DATA / MODALITIES / ROI to your preprocessed layout.
sphinx_gallery_thumbnail_number = 1
from pathlib import Path
import matplotlib.pyplot as plt
from habit.contracts import cohort_from_directory
from habit.datasets import fetch_demo
from habit.recipes import Study
from habit.spec import HabitatSpec, Spec, Stage
from habit.viz import plot_habitat_overlay
import numpy as np
DATA = fetch_demo()
# Three DCE phases: unenhanced, arterial, and portal-venous.
MODALITIES = ("pre_contrast", "LAP", "PVP")
ROI = "LAP"
cohort = cohort_from_directory(DATA, modalities=MODALITIES, roi=ROI)[:2]
print(f"Cohort: {list(cohort.subject_ids)}")
HABIT demo data (cached)
DATA (preprocessed root): C:\Users\dongm\.habit_data\demo-data-v1\preprocessed
On-disk inventory of this folder:
subjects (5): subj001, subj002, subj003, subj004, subj005
image series: LAP, PVP, delay_3min, pre_contrast
mask keys: LAP, PVP, delay_3min, pre_contrast
example image: images/subj001/delay_3min/WATER__BH_Ax_LAVA_Flex_3min_Series0012.nrrd
example mask: masks/subj001/delay_3min/WATER__BH_Ax_LAVA_Flex_10min_Series0017_mask.nrrd
Your own data must use the same folder tree (change IDs / series names):
DATA/
images/<subject_id>/<modality>/<one image file>
masks/<subject_id>/<roi>/<one mask file>
Then load it with the same call the demos use:
cohort = cohort_from_directory(DATA, modalities=("LAP",), roi="LAP")
Swap DATA / modalities / roi to match your tree. Mask key is often the
same as one image series (here LAP).
Cohort: ['subj001', 'subj002']
Fit one model on every ROI voxel
No partition: each ROI voxel is its own clustering unit. pool
puts the voxels of every subject together and one model is fit on them.
spec = HabitatSpec(
name="direct_pooling",
stages=(
# extract: one intensity column per DCE phase, inside the ROI.
Stage("extract", Spec("raw", {"modalities": list(MODALITIES), "roi": ROI})),
# No partition. pool stacks every subject's ROI voxels.
Stage("pool", Spec("pool")),
# fit: one shared k-means on those voxels. Count fixed at 3.
Stage("fit", Spec("kmeans", {"n_habitats": 3, "n_init": 10})),
Stage("assign", Spec("nearest_centroid")),
Stage("volume", Spec("volume")),
),
random_seed=0,
)
# One shared habitat_model, as in two-step, but its centroids were learnt
# on voxels rather than supervoxels.
result = Study(spec).fit_predict(cohort)
print(result.habitat_model.summary())
print(result.features.frame)
result.features.frame
Cohort.map[_ComputeUnits]: 0%| | 0/2 [00:00<?, ?it/s]
Cohort.map[_ComputeUnits]: 50%|█████ | 1/2 [00:01<00:01, 1.11s/it]
Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:02<00:00, 1.07s/it]
Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:02<00:00, 1.07s/it]
Cohort.map[_AssignPrecomputedUnits]: 0%| | 0/2 [00:00<?, ?it/s]
Cohort.map[_AssignPrecomputedUnits]: 50%|█████ | 1/2 [00:00<00:00, 3.31it/s]
Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00, 3.26it/s]
Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00, 3.26it/s]
HabitatModel kmeans-3fb1a8736ca6964a
habitats : 3
features (3) : pre_contrast, LAP, PVP
defining cohort : n=2
modalities : pre_contrast, LAP, PVP
cohort digest : cc1b16a34cc7478d...
produced by : habitat_model_fitter.kmeans
habit version : 3.0.0
random seed : 0
preprocessing state: inertia, validation
subject ... habitat_3_volume_fraction
0 subj001 ... 0.21926
1 subj002 ... 0.14932
[2 rows x 7 columns]
Intensities inside each habitat
The histogram uses the displayed ROI image, not a feature column.
Path("out").mkdir(exist_ok=True)
fig_hist, ax = plt.subplots(figsize=(6.2, 3.2))
colors = ["#4C78A8", "#F58518", "#54A24B"]
for habitat_id in sorted(
int(v) for v in np.unique(result.habitat_maps[0].label_array) if int(v) != 0
):
# Histogram uses the displayed ROI image (LAP), not a feature column.
values = cohort[0].image(ROI).data[
result.habitat_maps[0].label_array == habitat_id
]
ax.hist(values, bins=30, alpha=0.6, label=f"habitat {habitat_id}", color=colors[habitat_id - 1])
ax.set_xlabel(ROI)
ax.set_ylabel("voxels")
ax.set_title("pooled voxel intensities by habitat")
ax.legend()
fig_hist.savefig("out/pooling_intensity_hist.png", dpi=150, bbox_inches="tight")
plt.show()
One overlay per subject
for subject, habitat_map in zip(cohort, result.habitat_maps):
fig = plot_habitat_overlay(
subject.image(ROI),
habitat_map,
title=f"habitats ({habitat_map.subject_id})",
crop_to="labels",
)
fig.savefig(
f"out/pooling_{habitat_map.subject_id}.png",
dpi=150,
bbox_inches="tight",
)
plt.show()
F:\work\habit_project\habit\viz\habitat_overlay.py:806: UserWarning: Display geometry conflict: image/anatomy direction does not match mask/label direction. Using the mask/label direction so coronal/sagittal superior-up follows the labelled anatomy. Pass direction= to override.
resolved_direction, resolved_spacing = resolve_display_geometry(
F:\work\habit_project\habit\viz\habitat_overlay.py:806: UserWarning: Display geometry conflict: image/anatomy direction does not match mask/label direction. Using the mask/label direction so coronal/sagittal superior-up follows the labelled anatomy. Pass direction= to override.
resolved_direction, resolved_spacing = resolve_display_geometry(
Total running time of the script: (0 minutes 5.131 seconds)