.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples\06_matching\plot_01_label_switching.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_06_matching_plot_01_label_switching.py: Why habitat ids must be matched =============================== **Background.** Clustering algorithms such as k-means number their clusters in arbitrary order, so two fits of the same tumour can give the same region different habitat ids. Any comparison by id (Dice, volume per habitat, a cohort table) needs the ids matched first. **Purpose.** You get two k-means maps of one demo subject shown with raw and with matched ids, the voxel-overlap table that decides the pairing, and a per-habitat Dice table. **When to use.** Whenever two maps label the same voxels (a restart, a perturbed image, a second reader). For maps of different patients see :doc:`/auto_examples/06_matching/plot_03_prototype_steps`. **Key terms.** * **habitat** -- a sub-region inside the tumour (the ROI) whose voxels behave alike across the input images; each ROI voxel gets a habitat id (1, 2, 3, ...). * **label switching** -- independent clusterings number the same habitat differently (habitat 1 in one fit can be habitat 3 in another), so ids must be matched before comparing. * **overlap table** -- for every pair of habitats (one from each map), the number of voxels they share. * **Hungarian assignment** -- an exact algorithm (Kuhn-Munkres, SciPy's ``linear_sum_assignment``) that picks the one-to-one pairing of rows and columns with the best total; here, the largest total shared voxels. * **Dice** -- overlap between two label maps or regions (0 = none, 1 = identical). Cluster the same tumour twice with the same features and the same ``k``, changing only the k-means random seed. The two maps describe the same tissue, but k-means numbers its clusters in arbitrary order, so habitat 1 of one run is habitat 2 or 3 of the other. This is **label switching**. Comparing the maps by raw id reports disagreement that is not there. Because both maps label the same voxels, the ids are matched by voxel overlap: count the voxels shared by every pair of habitats, then pick the one-to-one pairing with the largest total (Hungarian assignment). .. GENERATED FROM PYTHON SOURCE LINES 46-51 Two k-means runs on one subject ------------------------------- Relative enhancement in the arterial (LAP) and portal-venous (PVP) phases, ``k = 3``, seeds 0 and 1. sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 51-89 .. code-block:: Python from pathlib import Path import matplotlib.pyplot as plt import numpy as np from habit.contracts import Cohort, cohort_from_directory from habit.datasets import fetch_demo from habit.habitat_model import KMeansHabitatModelFitter from habit.pipeline import voxel_units from habit.precision import align_habitat_map, habitat_stability from habit.viz import plot_habitat_label_compare, plot_label_overlap_matrix from habit.voxel_features import ExpressionVoxelFeatures # Change DATA / MODALITIES / ROI and the expressions to your own layout. DATA = fetch_demo() MODALITIES = ("pre_contrast", "LAP", "PVP") ROI = "LAP" subject = cohort_from_directory(DATA, modalities=MODALITIES, roi=ROI)[0] extractor = ExpressionVoxelFeatures( features={ "rel_enh_lap": "(LAP - pre_contrast) / (pre_contrast + eps)", "rel_enh_pvp": "(PVP - pre_contrast) / (pre_contrast + eps)", }, roi=ROI, ) # One row per ROI voxel; both runs cluster exactly these rows. units = voxel_units(extractor(subject)) runs = [] for seed in (0, 1): # n_init=1: a single k-means start per seed, so only the seed differs. fitter = KMeansHabitatModelFitter(n_habitats=3, n_init=1) fitter.set_random_state(seed) model = fitter.fit([units], cohort=Cohort([subject], name=subject.subject_id)) runs.append((model, model.assigner()(units))) print(f"seed {seed} centroids:\n{np.asarray(model.centroids).round(3)}") (model_a, map_a), (model_b, map_b) = runs .. rst-class:: sphx-glr-script-out .. code-block:: none 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/// masks/// 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). seed 0 centroids: [[1.484 1.733] [0.98 0.928] [1.616 1.071]] seed 1 centroids: [[0.983 0.926] [1.619 1.074] [1.48 1.735]] .. GENERATED FROM PYTHON SOURCE LINES 90-92 Same centroids, different order: row 1 of seed 0 is row 3 of seed 1. Raw ids therefore disagree on most voxels. .. GENERATED FROM PYTHON SOURCE LINES 92-97 .. code-block:: Python labels_a = np.asarray(map_a.label_array) labels_b = np.asarray(map_b.label_array) roi = labels_a > 0 print(f"voxels with the same raw id: {np.mean(labels_a[roi] == labels_b[roi]):.1%}") .. rst-class:: sphx-glr-script-out .. code-block:: none voxels with the same raw id: 0.0% .. GENERATED FROM PYTHON SOURCE LINES 98-102 Raw ids side by side -------------------- ``align_labels=False`` draws both maps as they came out of k-means: the colours differ although the regions are the same. .. GENERATED FROM PYTHON SOURCE LINES 102-115 .. code-block:: Python Path("out").mkdir(exist_ok=True) image = subject.image(ROI) fig = plot_habitat_label_compare( image, map_a, map_b, titles=("seed 0", "seed 1 (raw ids)"), align_labels=False, crop_to="labels", ) fig.savefig("out/switching_raw_ids.png", dpi=150, bbox_inches="tight") plt.show() .. image-sg:: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_001.png :alt: Habitat label compare, seed 0, seed 1 (raw ids), Disagreement :srcset: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none F:\work\habit_project\habit\viz\habitat_core.py:948: 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( .. GENERATED FROM PYTHON SOURCE LINES 116-122 The overlap table decides the pairing ------------------------------------- Each cell counts voxels in moving habitat ``i`` (seed 1) and reference habitat ``j`` (seed 0). The outlined cells are the Hungarian pairs: every habitat gets exactly one partner and the outlined total is the largest possible. .. GENERATED FROM PYTHON SOURCE LINES 122-126 .. code-block:: Python fig = plot_label_overlap_matrix(map_a, map_b, reference_name="seed 0", moving_name="seed 1") fig.savefig("out/switching_overlap_matrix.png", dpi=150, bbox_inches="tight") plt.show() .. image-sg:: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_002.png :alt: Voxel overlap (outlined: Hungarian pairs) :srcset: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 127-132 Rename seed 1 into the seed-0 ids --------------------------------- Both runs share a ``model_id`` (it digests the spec and subject, not the random seed), so :func:`~habit.precision.align_habitat_map` would treat them as one definition and do nothing. ``force=True`` aligns anyway. .. GENERATED FROM PYTHON SOURCE LINES 132-148 .. code-block:: Python print(f"same model_id: {model_a.model_id == model_b.model_id}") aligned_b = align_habitat_map(map_a, map_b, force=True) same = np.mean(labels_a[roi] == np.asarray(aligned_b.label_array)[roi]) print(f"voxels with the same id after matching: {same:.1%}") fig = plot_habitat_label_compare( image, map_a, aligned_b, titles=("seed 0", "seed 1 (matched ids)"), align_labels=False, crop_to="labels", ) fig.savefig("out/switching_matched_ids.png", dpi=150, bbox_inches="tight") plt.show() .. image-sg:: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_003.png :alt: Habitat label compare, seed 0, seed 1 (matched ids), Disagreement :srcset: /auto_examples/06_matching/images/sphx_glr_plot_01_label_switching_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none same model_id: True voxels with the same id after matching: 99.6% F:\work\habit_project\habit\viz\habitat_core.py:948: 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( .. GENERATED FROM PYTHON SOURCE LINES 149-154 Dice per habitat ---------------- :func:`~habit.precision.habitat_stability` pairs the ids by the same overlap rule and scores each pair. Pass the *original* moving map; it matches internally. .. GENERATED FROM PYTHON SOURCE LINES 154-157 .. code-block:: Python dice = habitat_stability(map_a, [map_b]) print(dice.to_string(index=False)) dice .. rst-class:: sphx-glr-script-out .. code-block:: none perturbation habitat_id matched_id dice n_reference n_matched 0 1 3 0.995383 7160 7136 0 2 1 0.996708 14096 14151 0 3 2 0.995493 13438 13407 .. raw:: html
perturbation habitat_id matched_id dice n_reference n_matched
0 0 1 3 0.995383 7160 7136
1 0 2 1 0.996708 14096 14151
2 0 3 2 0.995493 13438 13407


.. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.507 seconds) .. _sphx_glr_download_auto_examples_06_matching_plot_01_label_switching.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_label_switching.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_label_switching.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_label_switching.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_