Python implementation for StatisticsSurface (#669)

* Initial work * Run make format * Add brainstat to runtime dependencies * Fix typo in docstring Co-authored-by: omar-rifai <[email protected]> * Fix type mismatch warning * Fix type mismatch in function signature * Fix unresolved reference error * Fix extra hash in comment (PEP8) * Fix typo in docstring * Remove unused import * Use set literal * Use lowercase for inner function variable names * Fix typo in docstring * Add missing apostrophe * Use lowercase for inner function variable names * Fix type mismatch between plot functions * Silence broad try except warning * Ensure nilearn is installed with plotting capabilities * verbose + print -> cprint * improve docstring of clinica_surfstat * Rewrite... * fix missing Path * Start adding unit tests * Add more unit tests * Add more unit tests again * Linting... * fix * Linting... * Get rid of Templates * Typo in file name * Fix broken tests * Get rid of private call * some fixes * Fix broken test * Initial work * Run make format * Add brainstat to runtime dependencies * Fix typo in docstring Co-authored-by: omar-rifai <[email protected]> * Fix type mismatch warning * Fix type mismatch in function signature * Fix unresolved reference error * Fix extra hash in comment (PEP8) * Fix typo in docstring * Remove unused import * Use set literal * Use lowercase for inner function variable names * Fix typo in docstring * Add missing apostrophe * Use lowercase for inner function variable names * Fix type mismatch between plot functions * Silence broad try except warning * Ensure nilearn is installed with plotting capabilities * verbose + print -> cprint * improve docstring of clinica_surfstat * Rewrite... * fix missing Path * Start adding unit tests * Add more unit tests * Add more unit tests again * Linting... * fix * Linting... * Get rid of Templates * Typo in file name * Fix broken tests * Get rid of private call * some fixes * Fix broken test * Run formatter * Update lock file * Enable testing with Python 3.10 Add explicit install of a prerelease version of VTK with wheels provided for Python 3.10. * Use VTK 9.2 compatible BrainSpace * improve non regression test * update lock file * Drop version restrictions * Remove unused import * Remove unused imports * Remove self-resolving imports * Fix inconsistency between argname and docstring * Add static typing to write dispatcher * Add static typing to plot dispatcher * Fix typos and wrap docstring * Replace function call by set literal * Remove unused import * Enable passing a surface file to clinica_surfstat * use snake case style for attributes * Update clinica/pipelines/statistics_surface/_model.py Co-authored-by: Ghislain Vaillant <[email protected]> * remove comment * Turn GLM results classes into dataclasses * Refactor Results classes with better serializing API * Update clinica/pipelines/statistics_surface/clinica_surfstat.py Co-authored-by: Ghislain Vaillant <[email protected]> * Update clinica/pipelines/statistics_surface/_model.py Co-authored-by: Ghislain Vaillant <[email protected]> * Update clinica/pipelines/statistics_surface/_model.py Co-authored-by: Ghislain Vaillant <[email protected]> * Update clinica/pipelines/statistics_surface/_model.py Co-authored-by: Ghislain Vaillant <[email protected]> * Refactor _is_categorical for single responsability * Some fixes... * Replace GLMFactory class with create_glm_model function * Remove use of kwargs and pass parameters explicitely * Update minimum version of brainspace * Fix tests - reuse lambda sum aggregation * Add types to docstrings to comply with Numpydoc specs * Simplify _read_and_check_tsv_file Co-authored-by: Ghislain Vaillant <[email protected]> Co-authored-by: Ghislain Vaillant <[email protected]> Co-authored-by: omar-rifai <[email protected]>
aramis-lab · Sep 20, 2022 · 07d045f · 07d045f
1 parent 32bfdb0
commit 07d045f
Show file tree

Hide file tree

Showing 12 changed files with 2,970 additions and 107 deletions.
diff --git a/clinica/pipelines/statistics_surface/_inputs.py b/clinica/pipelines/statistics_surface/_inputs.py
@@ -0,0 +1,171 @@
+"""This file contains functions for loading data from disk
+and performing some checks on them.
+"""
+from os import PathLike
+from pathlib import Path
+from typing import Dict, Tuple
+
+import numpy as np
+import pandas as pd
+from nilearn.surface import Mesh
+
+TSV_FIRST_COLUMN = "participant_id"
+TSV_SECOND_COLUMN = "session_id"
+
+
+def _read_and_check_tsv_file(tsv_file: PathLike) -> pd.DataFrame:
+    """This function reads the TSV file provided and performs some basic checks.
+
+    Parameters
+    ----------
+    tsv_file : PathLike
+        Path to the TSV file to open.
+
+    Returns
+    -------
+    tsv_data : pd.DataFrame
+        DataFrame obtained from the file.
+    """
+    try:
+        return pd.read_csv(tsv_file, sep="\t").set_index(
+            [TSV_FIRST_COLUMN, TSV_SECOND_COLUMN]
+        )
+    except FileNotFoundError:
+        raise FileNotFoundError(f"File {tsv_file} does not exist.")
+    except KeyError:
+        raise ValueError(
+            f"The TSV data should have at least two columns: {TSV_FIRST_COLUMN} and {TSV_SECOND_COLUMN}"
+        )
+
+
+def _get_t1_freesurfer_custom_file_template(base_dir: PathLike) -> str:
+    """Returns a Template for the path to the desired surface file.
+
+    Parameters
+    ----------
+    base_dir : PathLike
+        Base directory to seach for the template.
+
+    Returns
+    -------
+    template_path : str
+        Path to the t1 freesurfer template.
+    """
+    return str(base_dir) + (
+        "/%(subject)s/%(session)s/t1/freesurfer_cross_sectional/%(subject)s_%(session)s"
+        "/surf/%(hemi)s.thickness.fwhm%(fwhm)s.fsaverage.mgh"
+    )
+
+
+def _build_thickness_array(
+    input_dir: PathLike,
+    surface_file: str,
+    df: pd.DataFrame,
+    fwhm: float,
+) -> np.ndarray:
+    """This function builds the cortical thickness array.
+
+    Parameters
+    ----------
+    input_dir : PathLike
+        Input directory.
+
+    surface_file : str
+        Template for the path to the surface file of interest.
+
+    df : pd.DataFrame
+        Subjects DataFrame.
+
+    fwhm : float
+        Smoothing parameter only used to retrieve the right surface file.
+
+    Returns
+    -------
+    thickness : np.ndarray
+        Cortical thickness. Hemispheres and subjects are stacked.
+    """
+    from nibabel.freesurfer.mghformat import load
+
+    thickness = []
+    for idx, row in df.iterrows():
+        subject = row[TSV_FIRST_COLUMN]
+        session = row[TSV_SECOND_COLUMN]
+        parts = []
+        for hemi in ["lh", "rh"]:
+            query = {"subject": subject, "session": session, "fwhm": fwhm, "hemi": hemi}
+            parts.append(
+                load(str(Path(input_dir) / Path(surface_file % query))).get_fdata()
+            )
+        combined = np.vstack(parts)
+        thickness.append(combined.flatten())
+    thickness = np.vstack(thickness)
+    if thickness.shape[0] != len(df):
+        raise ValueError(
+            f"Unexpected shape for thickness array : {thickness.shape}. "
+            f"Expected {len(df)} rows."
+        )
+    return thickness
+
+
+def _get_average_surface(fsaverage_path: PathLike) -> Tuple[Dict, Mesh]:
+    """This function extracts the average surface and the average mesh
+    from the path to the fsaverage templates.
+
+    .. note::
+
+        Note that the average surface is returned as a dictionary
+        with 'coord' and 'tri' as keys, while the average mesh is
+        returned as a Nilearn Mesh object (basically a NamedTuple
+        with 'coordinates' and 'faces' attributes). The surface
+        isn't returned as a Nilearn Surface object for compatibility
+        with BrainStats.
+
+    .. warning::
+
+        There is an issue with faces having  a value of 0 as index.
+        This is most likely a bug in BrainStat as MATLAB indexing
+        starts at 1 while Python starts at zero.
+
+    Parameters
+    ----------
+    fsaverage_path : PathLike
+        Path to the fsaverage templates.
+
+    Returns
+    -------
+    average_surface : dict
+        Average surface as a dictionary for BrainStat compatibility.
+
+    average_mesh : nilearn.surface.Mesh
+        Average mesh as a Nilearn Mesh object.
+    """
+    import copy
+
+    from nilearn.surface import Mesh, load_surf_mesh
+
+    meshes = [
+        load_surf_mesh(str(fsaverage_path / Path(f"{hemi}.pial")))
+        for hemi in ["lh", "rh"]
+    ]
+    coordinates = np.vstack([mesh.coordinates for mesh in meshes])
+    faces = np.vstack(
+        [meshes[0].faces, meshes[1].faces + meshes[0].coordinates.shape[0]]
+    )
+    average_mesh = Mesh(
+        coordinates=coordinates,
+        faces=copy.deepcopy(faces),
+    )
+    ##################
+    # UGLY HACK !!! Need investigation
+    ##################
+    # Uncomment the following line if getting an error
+    # with negative values in bincount in Brainstat.
+    # Not sure, but might be a bug in BrainStat...
+    #
+    faces += 1
+    #################
+    average_surface = {
+        "coord": coordinates,
+        "tri": faces,
+    }
+    return average_surface, average_mesh