Expression voxel features

Background. Clinicians often read DCE images as ratios (how much a voxel enhances, how fast it washes out) rather than raw signal. An expression extractor computes such ratios per voxel from a formula string, without writing a class.

Purpose. You get a voxel-feature table with four derived columns and a map of arterial relative enhancement inside the ROI.

When to use. When each new column is simple arithmetic on your modalities; for anything more, see Custom features.

Key terms.

• voxel feature / extract – see Extracting voxel intensities.

• derived map – a voxel-wise image computed from other images, such as relative enhancement (LAP - pre_contrast) / pre_contrast.

ExpressionVoxelFeatures builds one column per formula. Names in the formula are modality keys on the subject (LAP, pre_contrast, …). eps keeps a zero denominator from breaking a ratio.

Clustering those columns is a later page (Clustering habitats from a derived map).

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.viz import plot_voxel_texture_slice
from habit.voxel_features import ExpressionVoxelFeatures

DATA = fetch_demo()
# All four DCE phases: unenhanced, arterial, portal-venous, delayed.
MODALITIES = ("pre_contrast", "LAP", "PVP", "delay_3min")
ROI = "LAP"
subject = cohort_from_directory(DATA, modalities=MODALITIES, roi=ROI)[0]
print(subject.subject_id, sorted(subject.images))

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).
subj001 ['LAP', 'PVP', 'delay_3min', 'pre_contrast']

Each key becomes one column of the voxel field. Only ROI voxels are evaluated; the formulas use the modality keys above.

extractor = ExpressionVoxelFeatures(
    features={
        "rel_enh_lap": "(LAP - pre_contrast) / (pre_contrast + eps)",
        "rel_enh_pvp": "(PVP - pre_contrast) / (pre_contrast + eps)",
        "rel_enh_delay": "(delay_3min - pre_contrast) / (pre_contrast + eps)",
        "washout": "(LAP - delay_3min) / (LAP + eps)",
    },
    roi=ROI,
)
field = extractor(subject)
print(f"{field.values.shape[0]} voxels, columns: {list(field.feature_names)}")
print(field.feature_frame().head())

34694 voxels, columns: ['rel_enh_lap', 'rel_enh_pvp', 'rel_enh_delay', 'washout']
   rel_enh_lap  rel_enh_pvp  rel_enh_delay   washout
0     1.388759     1.412178       1.330211  0.024510
1     1.324201     1.438356       1.269406  0.023576
2     1.453303     1.478360       1.243736  0.085422
3     1.353070     1.350877       1.234649  0.050326
4     1.188559     1.398305       1.175847  0.005808

Arterial relative enhancement inside the ROI.

Path("out").mkdir(exist_ok=True)
fig = plot_voxel_texture_slice(
    field,
    feature=0,
    anatomy=subject.image(ROI),
    roi_mask=subject.mask(ROI),
    title="relative enhancement (LAP)",
    crop_to="roi",
)
fig.savefig("out/expression_relative_enhancement.png", dpi=150, bbox_inches="tight")
plt.show()

F:\work\habit_project\habit\viz\voxel_texture.py:1019: 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.
  direction, spacing = resolve_display_geometry(