.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\04_designs\plot_03_pool_voxels.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_04_designs_plot_03_pool_voxels.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_04_designs_plot_03_pool_voxels.py:


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
(:doc:`/auto_examples/04_designs/plot_01_two_step`) clusters far fewer rows.

**Key terms.**

* **direct pooling** -- ``pool`` then ``fit`` with no ``partition``; each
  voxel is its own clustering unit.
* **pool** -- see :doc:`/auto_examples/04_designs/plot_01_two_step`.

Input: a cohort of at least two subjects. Output: one shared
:class:`~habit.contracts.HabitatModel` fitted on voxels, and one
:class:`~habit.contracts.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.

.. GENERATED FROM PYTHON SOURCE LINES 32-36

Load the cohort
---------------
Change ``DATA`` / ``MODALITIES`` / ``ROI`` to your preprocessed layout.
sphinx_gallery_thumbnail_number = 1

.. GENERATED FROM PYTHON SOURCE LINES 36-54

.. code-block:: Python

    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)}")





.. 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/<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']




.. GENERATED FROM PYTHON SOURCE LINES 55-59

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.

.. GENERATED FROM PYTHON SOURCE LINES 59-80

.. code-block:: Python

    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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    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]


.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>subject</th>
          <th>habitat_1_voxel_count</th>
          <th>habitat_1_volume_fraction</th>
          <th>habitat_2_voxel_count</th>
          <th>habitat_2_volume_fraction</th>
          <th>habitat_3_voxel_count</th>
          <th>habitat_3_volume_fraction</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>subj001</td>
          <td>11618.0</td>
          <td>0.334871</td>
          <td>15469.0</td>
          <td>0.445870</td>
          <td>7607.0</td>
          <td>0.21926</td>
        </tr>
        <tr>
          <th>1</th>
          <td>subj002</td>
          <td>5626.0</td>
          <td>0.570704</td>
          <td>2760.0</td>
          <td>0.279976</td>
          <td>1472.0</td>
          <td>0.14932</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 81-84

Intensities inside each habitat
-------------------------------
The histogram uses the displayed ROI image, not a feature column.

.. GENERATED FROM PYTHON SOURCE LINES 84-102

.. code-block:: Python

    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()




.. image-sg:: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_001.png
   :alt: pooled voxel intensities by habitat
   :srcset: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 103-105

One overlay per subject
-----------------------

.. GENERATED FROM PYTHON SOURCE LINES 105-118

.. code-block:: Python

    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()



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_002.png
         :alt: habitats (subj001), Axis 0 (axial-like) @ 96, Axis 1 (coronal-like) @ 165, Axis 2 (sagittal-like) @ 71
         :srcset: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_003.png
         :alt: habitats (subj002), Axis 0 (axial-like) @ 86, Axis 1 (coronal-like) @ 233, Axis 2 (sagittal-like) @ 69
         :srcset: /auto_examples/04_designs/images/sphx_glr_plot_03_pool_voxels_003.png
         :class: sphx-glr-multi-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    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(





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 5.131 seconds)


.. _sphx_glr_download_auto_examples_04_designs_plot_03_pool_voxels.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_03_pool_voxels.ipynb <plot_03_pool_voxels.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_03_pool_voxels.py <plot_03_pool_voxels.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_03_pool_voxels.zip <plot_03_pool_voxels.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_