.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples\04_designs\plot_01_two_step.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_01_two_step.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_04_designs_plot_01_two_step.py:


Defining habitats in two steps
==============================

**Background.** The two-step design is the complete analysis used across
this Guide: each tumour is
first cut into supervoxels, then the supervoxels of all subjects are
clustered together once, so habitat ids mean the same thing in every
patient.

**Purpose.** You get one shared habitat model, a habitat map and a
volume-fraction table per subject, a supervoxel / habitat triptych, and an
elbow plot for choosing the habitat count.

**When to use.** You want habitats that are comparable across patients
(e.g. volume fractions as cohort features) and whole-cohort voxel
clustering would be slow or noisy.

**Key terms.**

* **supervoxel** -- a small patch of neighbouring voxels with similar
  features, clustered inside one subject first (``partition``); see
  :doc:`/auto_examples/02_stages/plot_12_supervoxels`.
* **pool** -- stacks every training subject's rows into one matrix so one
  model is fitted to the whole cohort.
* **SLIC** -- a supervoxel method that grows compact regions from a regular
  grid of seeds, trading intensity similarity against spatial distance.
* **elbow** -- see :doc:`/auto_quickstart/plot_quickstart_python`.

Input: a cohort of at least two subjects. Output: one shared
:class:`~habit.contracts.HabitatModel` and one
:class:`~habit.contracts.HabitatMap` per subject. The design is the
stage list: ``partition`` then ``pool`` then ``fit``.

SLIC is the partitioner here. ``Spec("kmeans", ...)`` in the
``partition`` stage is the same study with a different partitioner, not
a second habitat definition to compare on this page. A fixed count of 3
draws the map; a 2-10 search draws the elbow on the same cohort.
``two_step_habitat(...)`` is a shortcut that builds the same stage list.

.. GENERATED FROM PYTHON SOURCE LINES 43-47

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

.. GENERATED FROM PYTHON SOURCE LINES 47-64

.. 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_cluster_validation_from_report, plot_partition_triptych

    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 65-71

Fit three shared habitats
-------------------------
``partition`` splits each ROI into supervoxels. Those supervoxels are
the clustering units: one mean feature vector per supervoxel. ``pool``
puts the units of all subjects together and ``fit`` learns one shared
habitat model on them.

.. GENERATED FROM PYTHON SOURCE LINES 71-99

.. code-block:: Python

    first_stages = (
        # extract: one intensity column per DCE phase, voxels inside the ROI.
        Stage("extract", Spec("raw", {"modalities": list(MODALITIES), "roi": ROI})),
        # partition: SLIC supervoxels. These rows, not voxels, are clustered.
        Stage("partition", Spec("slic", {"n_supervoxels": 100})),
        # pool: stack every subject's supervoxels before a single fit.
        Stage("pool", Spec("pool")),
    )
    spec = HabitatSpec(
        name="two_step_slic",
        stages=first_stages + (
            # fit: one shared k-means. n_init=10 restarts; count is fixed at 3.
            Stage("fit", Spec("kmeans", {"n_habitats": 3, "n_init": 10})),
            # assign: nearest centroid paints habitat ids back onto voxels.
            Stage("assign", Spec("nearest_centroid")),
            # quantify: per-subject volume fractions of those ids.
            Stage("volume", Spec("volume")),
        ),
        # Seeds the stochastic stages (here the k-means fit) so a rerun paints the same map.
        random_seed=0,
    )
    # fit_predict: learn the shared centroids on the cohort, then label every
    # subject with them in the same call.
    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:23<00:23, 23.91s/it]
    Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:46<00:00, 23.40s/it]
    Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:46<00:00, 23.40s/it]

    Cohort.map[_AssignPrecomputedUnits]:   0%|          | 0/2 [00:00<?, ?it/s]
    Cohort.map[_AssignPrecomputedUnits]:  50%|█████     | 1/2 [00:00<00:00,  3.24it/s]
    Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00,  3.32it/s]
    Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00,  3.32it/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.271632
    1  subj002  ...                   0.167174

    [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>7925.0</td>
          <td>0.228426</td>
          <td>17345.0</td>
          <td>0.499942</td>
          <td>9424.0</td>
          <td>0.271632</td>
        </tr>
        <tr>
          <th>1</th>
          <td>subj002</td>
          <td>5041.0</td>
          <td>0.511361</td>
          <td>3169.0</td>
          <td>0.321465</td>
          <td>1648.0</td>
          <td>0.167174</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 100-102

Supervoxels, clustering units, and habitats
-------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 102-112

.. code-block:: Python

    fig = plot_partition_triptych(
        cohort[0].image(ROI),
        result.units[0],
        result.habitat_maps[0],
        axis=0,
    )
    Path("out").mkdir(exist_ok=True)
    fig.savefig("out/two_step_triptych.png", dpi=150, bbox_inches="tight")
    plt.show()




.. image-sg:: /auto_examples/04_designs/images/sphx_glr_plot_01_two_step_001.png
   :alt: Two-step partitions, Anatomy, Supervoxels, Habitats
   :srcset: /auto_examples/04_designs/images/sphx_glr_plot_01_two_step_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:1110: 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 113-116

Elbow on the same stages, with ``fit`` searching 2-10 habitats. The
chosen count is the knee, not a second overlay. A saved model is
:doc:`/auto_examples/05_apply/plot_04_apply_saved_model`.

.. GENERATED FROM PYTHON SOURCE LINES 116-132

.. code-block:: Python

    elbow_spec = HabitatSpec(
        name="two_step_slic_elbow",
        stages=first_stages + (
            Stage("fit", Spec("kmeans", {"min_habitats": 2, "max_habitats": 10, "validation": "elbow", "n_init": 10})),
            Stage("assign", Spec("nearest_centroid")),
        ),
        random_seed=0,
    )
    elbow = Study(elbow_spec).fit_predict(cohort)
    # The fitter stores its per-count scores inside the model, so the chosen
    # count can be audited later.
    report = (elbow.habitat_model.preprocessing_state or {}).get("selection_report")
    print(elbow.habitat_model.summary())
    fig_k = plot_cluster_validation_from_report(report)
    fig_k.savefig("out/two_step_elbow.png", dpi=150, bbox_inches="tight")
    plt.show()



.. image-sg:: /auto_examples/04_designs/images/sphx_glr_plot_01_two_step_002.png
   :alt: Cluster validation curves, elbow (selected k=4)
   :srcset: /auto_examples/04_designs/images/sphx_glr_plot_01_two_step_002.png
   :class: sphx-glr-single-img


.. 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:23<00:23, 23.36s/it]
    Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:46<00:00, 23.32s/it]
    Cohort.map[_ComputeUnits]: 100%|██████████| 2/2 [00:46<00:00, 23.32s/it]

    Cohort.map[_AssignPrecomputedUnits]:   0%|          | 0/2 [00:00<?, ?it/s]
    Cohort.map[_AssignPrecomputedUnits]:  50%|█████     | 1/2 [00:00<00:00,  5.00it/s]
    Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00,  5.00it/s]
    Cohort.map[_AssignPrecomputedUnits]: 100%|██████████| 2/2 [00:00<00:00,  4.99it/s]
    HabitatModel kmeans-3432f3b65bfec11a
      habitats           : 4
      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, selection_report, validation





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

   **Total running time of the script:** (1 minutes 36.125 seconds)


.. _sphx_glr_download_auto_examples_04_designs_plot_01_two_step.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_two_step.ipynb <plot_01_two_step.ipynb>`

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

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

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

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


.. only:: html

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

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