.. _api_ph_computer:

===========
PH Computer
===========

.. currentmodule:: medtda.ph_computer

The ``ph_computer`` module provides functions for computing persistent homology from images using cubical complexes.

.. contents:: Contents
   :local:
   :depth: 2

Overview
========

This module contains the core functions for computing persistent homology (PH) from image data. It uses the GUDHI library's cubical complex implementation.

**Key Features:**

* Sublevel and superlevel filtrations
* T-construction and V-construction cubical complexes
* Multi-dimensional homology (H0, H1, H2, ...)
* Automatic dimension detection
* Infinite persistence handling

Basic Usage
===========

Quick Start
-----------

::

    from medtda.ph_computer import compute_ph
    import numpy as np

    # Create a simple 2D image
    image = np.random.rand(50, 50)

    # Compute persistence barcodes
    barcodes = compute_ph(image, filtration_type='sublevel')

    # barcodes is a dict: {'H0': H0 pairs, 'H1': H1 pairs, ...}
    print(f"H0 features: {len(barcodes['H0'])}")
    print(f"H1 features: {len(barcodes['H1'])}")

With Custom Parameters
----------------------

::

    barcodes = compute_ph(
        image,
        filtration_type='superlevel',
        construction='V',
        max_dimension=2
    )

Functions
=========

compute_ph
----------

.. autofunction:: medtda.ph_computer.compute_ph

Compute persistent homology from an image.

**Signature:**

.. code-block:: python

   def compute_ph(
       image: np.ndarray,
       filtration_type: str = 'sublevel',
       construction: str = 'T',
       max_dimension: int = -1
   ) -> Dict[int, np.ndarray]

**Parameters:**

* **image** (*np.ndarray*) - Input image (2D or 3D)
* **filtration_type** (*str*) - ``'sublevel'`` or ``'superlevel'`` (default: ``'sublevel'``)
* **construction** (*str*) - ``'T'`` or ``'V'`` for cubical complex construction (default: ``'T'``)
* **max_dimension** (*int*) - Maximum homology dimension to compute, ``-1`` for auto (default: ``-1``)

**Returns:**

* **barcodes** (*dict*) - Dictionary mapping dimension to barcode array
  
  * Keys: homology dimensions (0, 1, 2, ...)
  * Values: NumPy arrays of shape ``(n_features, 2)`` with ``[birth, death]`` columns

**Example:**

::

    import numpy as np
    from medtda.ph_computer import compute_ph

    # 3D medical image
    image = np.random.rand(64, 64, 64)

    # Compute H0, H1, H2
    barcodes = compute_ph(
        image,
        filtration_type='sublevel',
        max_dimension=2
    )

    # Access barcodes by dimension
    h0_barcode = barcodes['H0']  # Connected components
    h1_barcode = barcodes['H1']  # Loops/tunnels
    h2_barcode = barcodes['H2']  # Voids/cavities

    # Each barcode is an array of [birth, death] pairs
    for birth, death in h1_barcode:
        persistence = death - birth
        print(f"H1 feature: birth={birth:.3f}, death={death:.3f}, persistence={persistence:.3f}")

clean_and_cap_ph
----------------

.. autofunction:: medtda.ph_computer.clean_and_cap_ph

Clean and cap persistence values to handle infinite persistence.

**Signature:**

.. code-block:: python

   def clean_and_cap_ph(
       ph_array: np.ndarray,
       cap_value: float
   ) -> np.ndarray

**Parameters:**

* **ph_array** (*np.ndarray*) - Raw persistence barcode array
* **cap_value** (*float*) - Value to cap infinite persistence

**Returns:**

* **cleaned_barcode** (*np.ndarray*) - Cleaned barcode with finite values

**Example:**

::

    from medtda.ph_computer import compute_ph, clean_and_cap_ph

    barcodes = compute_ph(image)
    
    # H0 often has infinite persistence features
    h0_raw = barcodes['H0']
    h0_clean = clean_and_cap_ph(h0_raw, cap_value=1.0)

Parameters
==========

Filtration Type
---------------

**sublevel**

Build filtration from low to high intensity values.

* Features appear at low intensities (birth)
* Features disappear at high intensities (death)
* Use for: bright features on dark background

**superlevel**

Build filtration from high to low intensity values.

* Features appear at high intensities
* Features disappear at low intensities  
* Use for: dark features on bright background

Example::

    # Dark lesions on bright CT scan
    barcodes_sublevel = compute_ph(image, filtration_type='sublevel')

    # Bright lesions on dark background
    barcodes_superlevel = compute_ph(image, filtration_type='superlevel')

Construction Method
-------------------

**T-construction (default)**

* Each voxel becomes a top-dimensional cell
* More features, finer detail
* Slightly slower computation

**V-construction**

* Each voxel becomes a vertex
* Fewer features, coarser detail
* Slightly faster computation

Example::

    # T-construction: finer detail
    barcodes_t = compute_ph(image, construction='T')

    # V-construction: faster computation
    barcodes_v = compute_ph(image, construction='V')

Max Dimension
-------------

**-1 (auto)**

Automatically determines max dimension from image:

* 2D images → max_dimension = 1 (H0, H1)
* 3D images → max_dimension = 2 (H0, H1, H2)

**0**

Compute only H0 (connected components)::

    barcodes = compute_ph(image, max_dimension=0)
    # Returns {0: h0_barcode}

**1**

Compute H0 and H1::

    barcodes = compute_ph(image, max_dimension=1)
    # Returns {0: h0_barcode, 1: h1_barcode}

**2**

Compute H0, H1, and H2::

    barcodes = compute_ph(image, max_dimension=2)
    # Returns {0: h0_barcode, 1: h1_barcode, 2: h2_barcode}

Homology Dimensions
===================

H0: Connected Components
-------------------------

Represents connected regions in the image::

    h0_barcode = barcodes['H0']
    
    # Number of significant connected components
    # (persistence > threshold)
    threshold = 0.1
    significant_components = np.sum(
        (h0_barcode[:, 1] - h0_barcode[:, 0]) > threshold
    )

H1: Loops and Tunnels
----------------------

Represents 1-dimensional holes::

    h1_barcode = barcodes['H1']
    
    # Features with high persistence indicate prominent loops
    persistent_loops = h1_barcode[
        (h1_barcode[:, 1] - h1_barcode[:, 0]) > 0.2
    ]

H2: Voids and Cavities
-----------------------

Represents 2-dimensional voids (3D images only)::

    h2_barcode = barcodes['H2']
    
    # Count significant cavities
    significant_voids = np.sum(
        (h2_barcode[:, 1] - h2_barcode[:, 0]) > 0.15
    )

Complete Examples
=================

Example 1: 2D Image Analysis
-----------------------------

::

    import numpy as np
    from medtda.ph_computer import compute_ph
    from medtda.loaders import load_image

    # Load 2D image
    image, metadata = load_image('xray.png')

    # Compute persistence
    barcodes = compute_ph(
        image,
        filtration_type='sublevel',
        construction='T',
        max_dimension=1  # H0 and H1 only for 2D
    )

    # Analyze connected components
    h0 = barcodes['H0']
    persistence_h0 = h0[:, 1] - h0[:, 0]
    
    print(f"Connected components: {len(h0)}")
    print(f"Mean persistence: {persistence_h0.mean():.3f}")
    print(f"Max persistence: {persistence_h0.max():.3f}")

    # Analyze loops
    h1 = barcodes['H1']
    persistence_h1 = h1[:, 1] - h1[:, 0]
    
    print(f"Loops detected: {len(h1)}")
    print(f"Significant loops (p>0.1): {np.sum(persistence_h1 > 0.1)}")

Example 2: 3D Medical Image
----------------------------

::

    import numpy as np
    from medtda.ph_computer import compute_ph
    from medtda import Preprocessor
    from medtda.loaders import load_image, load_mask

    # Load and preprocess
    image, _ = load_image('ct_scan.nii.gz')
    mask, _ = load_mask('organ_mask.nii.gz')

    preprocessor = Preprocessor(
        normalize=True,
        crop_to_roi=True
    )
    processed = preprocessor.preprocess(image, mask)

    # Compute all homology dimensions
    barcodes = compute_ph(
        processed,
        filtration_type='sublevel',
        construction='T',
        max_dimension=2
    )

    # Analyze 3D topology
    for dim, barcode in barcodes.items():
        persistence = barcode[:, 1] - barcode[:, 0]
        
        dim_names = {0: 'Connected Components', 1: 'Tunnels', 2: 'Voids'}
        print(f"\n{dim_names[dim]} (H{dim}):")
        print(f"  Total features: {len(barcode)}")
        print(f"  Mean persistence: {persistence.mean():.4f}")
        print(f"  Std persistence: {persistence.std():.4f}")
        print(f"  Max persistence: {persistence.max():.4f}")

Example 3: Comparing Filtrations
---------------------------------

::

    from medtda.ph_computer import compute_ph

    # Dark features (e.g., vessels on bright background)
    barcodes_sub = compute_ph(ct_image, filtration_type='sublevel')

    # Bright features (e.g., contrast-enhanced regions)
    barcodes_sup = compute_ph(ct_image, filtration_type='superlevel')

    # Compare
    print("Sublevel filtration:")
    print(f"  H1 features: {len(barcodes_sub[1])}")
    
    print("Superlevel filtration:")
    print(f"  H1 features: {len(barcodes_sup[1])}")

Performance Tips
================

1. **Choose appropriate max_dimension:**
   
   * For 2D images, max_dimension=1 is usually sufficient
   * Computing H2 on 3D images can be slow for large volumes

2. **Construction method:**
   
   * Use T-construction for detailed analysis
   * Use V-construction for faster computation on large images

3. **Image size:**
   
   * Downsampling or cropping to ROI can significantly speed up computation
   * Consider preprocessing with ``crop_to_roi=True``

4. **Infinite persistence:**
   
   * H0 always has one infinite persistence feature
   * Use ``clean_and_cap_ph()`` if needed for downstream analysis

See Also
========

* :doc:`barcodeextractor` - Higher-level barcode extraction class
* :doc:`featureextractor` - Complete feature extraction pipeline
* :doc:`../user_guide/persistent_homology` - PH concepts and usage guide
* :doc:`../theory/persistent_homology` - Mathematical background