Merge pull request #29 from Deep-MI/remove-freesurfer

Remove freesurfer dependencies and add tests

Author: m-reuter

.github/workflows/pytest.yml

@@ -12,14 +12,21 @@ on:
       - '**.py'
   workflow_dispatch:
 
+env:
+  SUBJECTS_DIR: /home/runner/work/BrainPrint/BrainPrint/data
+  SUBJECT_ID: test
+  DESTINATION_DIR: /home/runner/work/BrainPrint/BrainPrint/data
+
 jobs:
   pytest:
     timeout-minutes: 30
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu, macos, windows]
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+#        os: [ubuntu, macos, windows]
+#        python-version: [3.8, 3.9, "3.10", "3.11"]
+        os: [ubuntu]
+        python-version: ["3.10"]
     name: ${{ matrix.os }} - py${{ matrix.python-version }}
     runs-on: ${{ matrix.os }}-latest
     defaults:
@@ -32,12 +39,25 @@ jobs:
         uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python-version }}
+          architecture: 'x64'
       - name: Install package
         run: |
           python -m pip install --progress-bar off --upgrade pip setuptools wheel
           python -m pip install --progress-bar off .[test]
+      - name: Create data folders
+        run: |
+          mkdir -p data/test/mri
+          mkdir -p data/test/surf
+          mkdir -p data/test/temp
       - name: Display system information
         run: brainprint-sys_info --developer
+      - name: Download files
+        run: |
+          wget https://surfer.nmr.mgh.harvard.edu/pub/data/tutorial_data/buckner_data/tutorial_subjs/good_output/mri/aseg.mgz -O data/test/mri/aseg.mgz
+          wget https://surfer.nmr.mgh.harvard.edu/pub/data/tutorial_data/buckner_data/tutorial_subjs/good_output/surf/lh.white -O data/test/surf/lh.white
+          wget https://surfer.nmr.mgh.harvard.edu/pub/data/tutorial_data/buckner_data/tutorial_subjs/good_output/surf/rh.white -O data/test/surf/rh.white
+          wget https://surfer.nmr.mgh.harvard.edu/pub/data/tutorial_data/buckner_data/tutorial_subjs/good_output/surf/lh.pial -O data/test/surf/lh.pial
+          wget https://surfer.nmr.mgh.harvard.edu/pub/data/tutorial_data/buckner_data/tutorial_subjs/good_output/surf/rh.pial -O data/test/surf/rh.pial
       - name: Run pytest
         run: pytest brainprint --cov=brainprint --cov-report=xml --cov-config=pyproject.toml
       - name: Upload to codecov

README.md

@@ -12,7 +12,7 @@ cortical parcellations or label files).
 
 ## Installation
 
-Use the following code to install the latest release of LaPy into your local
+Use the following code to install the latest release into your local
 Python package directory:
 
 `python3 -m pip install brainprint`
@@ -88,7 +88,14 @@ asymmetry calculation is performed and/or for the eigenvectors (CLI `--evecs` fl
 
 ## Changes
 
-There are some changes in functionality in comparison to the original [BrainPrint](https://github.com/Deep-MI/BrainPrint-legacy)
+Since version 0.5.0, some changes break compatibility with earlier versions (0.4.0 and lower) as well as the [original BrainPrint](https://github.com/Deep-MI/BrainPrint-legacy). These changes include:
+
+- for the creation of surfaces from voxel-based segmentations, we have replaced FreeSurfer's marching cube algorithm by scikit-image's marching cube algorithm. Similarly, other FreeSurfer binaries have been replaced by custom Python functions. As a result, a parallel FreeSurfer installation is no longer a requirement for running the brainprint software.
+- we have changed / removed the following composite structures from the brainprint shape descriptor: the left and right *striatum* (composite of caudate, putamen, and nucleus accumbens) and the left and right *ventricles* (composite of lateral, inferior lateral, 3rd ventricle, choroid plexus, and CSF) have been removed; the left and right *cerebellum-white-matter* and *cerebellum-cortex* have been merged into left and right *cerebellum*.
+
+As a result of these changes, numerical values for the brainprint shape descriptor that are obtained from version 0.5.0 and higher are expected to differ from earlier versions when applied to the same data, but should remain highly correlated with earlier results.
+
+There are some changes in version 0.4.0 (and lower) in functionality in comparison to the original [BrainPrint](https://github.com/Deep-MI/BrainPrint-legacy)
 scripts:
 
 - currently no support for tetrahedral meshes

brainprint/asymmetry.py

@@ -25,25 +25,10 @@ def compute_asymmetry(
     dict[str, float]
         {left_label}_{right_label}, distance.
     """
-    # Define structures
-
-    # combined and individual aseg labels:
-    # - Left  Striatum: left  Caudate + Putamen + Accumbens
-    # - Right Striatum: right Caudate + Putamen + Accumbens
-    # - CorpusCallosum: 5 subregions combined
-    # - Cerebellum: brainstem + (left+right) cerebellum WM and GM
-    # - Ventricles: (left+right) lat.vent + inf.lat.vent + choroidplexus + 3rdVent + CSF
-    # - Lateral-Ventricle: lat.vent + inf.lat.vent + choroidplexus
-    # - 3rd-Ventricle: 3rd-Ventricle + CSF
 
     structures_left_right = [
-        ("Left-Striatum", "Right-Striatum"),
         ("Left-Lateral-Ventricle", "Right-Lateral-Ventricle"),
-        (
-            "Left-Cerebellum-White-Matter",
-            "Right-Cerebellum-White-Matter",
-        ),
-        ("Left-Cerebellum-Cortex", "Right-Cerebellum-Cortex"),
+        ("Left-Cerebellum", "Right-Cerebellum"),
         ("Left-Thalamus-Proper", "Right-Thalamus-Proper"),
         ("Left-Caudate", "Right-Caudate"),
         ("Left-Putamen", "Right-Putamen"),

brainprint/brainprint.py

@@ -16,8 +16,6 @@
 from .utils.utils import (
     create_output_paths,
     export_brainprint_results,
-    test_freesurfer,
-    validate_environment,
     validate_subject_dir,
 )
 
@@ -76,7 +74,6 @@ def compute_surface_brainprint(
     Parameters
     ----------
     path : Path
-
         Path to the *.vtk* surface path.
     return_eigenvectors : bool, optional
         Whether to store eigenvectors in the result (default is True).
@@ -249,8 +246,6 @@ def run_brainprint(
         - Eigenvectors
         - Distances
     """  # noqa: E501
-    validate_environment()
-    test_freesurfer()
     subject_dir = validate_subject_dir(subjects_dir, subject_id)
     destination = create_output_paths(
         subject_dir=subject_dir,
@@ -301,8 +296,6 @@ def __init__(
         asymmetry: bool = False,
         asymmetry_distance: str = "euc",
         keep_temp: bool = False,
-        environment_validation: bool = True,
-        freesurfer_validation: bool = True,
         use_cholmod: bool = False,
     ) -> None:
         """
@@ -353,11 +346,6 @@ def __init__(
         self._eigenvectors = None
         self._distances = None
 
-        if environment_validation:
-            validate_environment()
-        if freesurfer_validation:
-            test_freesurfer()
-
     def run(self, subject_id: str, destination: Path = None) -> dict[str, Path]:
         """
         Run Brainprint analysis for a specified subject.

brainprint/cli/__init__.py

@@ -1,6 +1,7 @@
 """
 BrainPrint analysis CLI.
 """
+
 from ..brainprint import run_brainprint
 from .parser import parse_options
 

brainprint/cli/help_text.py

@@ -1,6 +1,7 @@
 """
 Help text strings for the :mod:`brainprint.cli` module.
 """
+
 CLI_DESCRIPTION: str = (
     "This program conducts a brainprint analysis of FreeSurfer output."
 )
@@ -20,7 +21,9 @@
 ASYM_DISTANCE: str = (
     "Distance measurement to use for asymmetry calculation (default: euc)"
 )
-CHOLMOD: str = "Use cholesky decomposition (faster) instead of LU decomposition (slower). May require manual install of scikit-sparse library. Default is LU decomposition."
+CHOLMOD: str = (
+    "Use cholesky decomposition (faster) instead of LU decomposition (slower). May require manual install of scikit-sparse library. Default is LU decomposition."
+)
 KEEP_TEMP: str = (
     "Whether to keep the temporary files directory or not, by default False"
 )
@@ -44,14 +47,11 @@
 
 CorpusCallosum                  [251, 252, 253, 254, 255]
 Cerebellum                      [7, 8, 16, 46, 47]
-Ventricles                      [4, 5, 14, 24, 31, 43, 44, 63]
 3rd-Ventricle                   [14, 24]
 4th-Ventricle                   15
 Brain-Stem                      16
-Left-Striatum                   [11, 12, 26]
 Left-Lateral-Ventricle          [4, 5, 31]
-Left-Cerebellum-White-Matter    7
-Left-Cerebellum-Cortex          8
+Left-Cerebellum                 [7, 8]
 Left-Thalamus-Proper            10
 Left-Caudate                    11
 Left-Putamen                    12
@@ -60,10 +60,8 @@
 Left-Amygdala                   18
 Left-Accumbens-area             26
 Left-VentralDC                  28
-Right-Striatum                  [50, 51, 58]
 Right-Lateral-Ventricle         [43, 44, 63]
-Right-Cerebellum-White-Matter   46
-Right-Cerebellum-Cortex         47
+Right-Cerebellum                [46, 47]
 Right-Thalamus-Proper           49
 Right-Caudate                   50
 Right-Putamen                   51

brainprint/cli/utils.py

@@ -1,6 +1,7 @@
 """
 Utility functions for the :mod:`brainprint.cli` module.
 """
+
 from . import help_text
 
 

brainprint/surfaces.py

@@ -1,12 +1,14 @@
 """
 Utility module holding surface generation related functions.
 """
-import uuid
+import os
 from pathlib import Path
 
+import nibabel as nb
+import numpy as np
 from lapy import TriaMesh
-
-from .utils.utils import run_shell_command
+from scipy import sparse as sp
+from skimage.measure import marching_cubes
 
 
 def create_aseg_surface(
@@ -30,51 +32,76 @@ def create_aseg_surface(
         Path to the generated surface in VTK format.
     """
     aseg_path = subject_dir / "mri/aseg.mgz"
-    norm_path = subject_dir / "mri/norm.mgz"
-    temp_name = f"temp/aseg.{uuid.uuid4()}"
+    temp_name = "temp/aseg.{indices}".format(indices="_".join(indices))
     indices_mask = destination / f"{temp_name}.mgz"
-    # binarize on selected labels (creates temp indices_mask)
-    # always binarize first, otherwise pretess may scale aseg if labels are
-    # larger than 255 (e.g. aseg+aparc, bug in mri_pretess?)
-    binarize_template = "mri_binarize --i {source} --match {match} --o {destination}"
-    binarize_command = binarize_template.format(
-        source=aseg_path, match=" ".join(indices), destination=indices_mask
-    )
-    run_shell_command(binarize_command)
 
-    label_value = "1"
-    # if norm exist, fix label (pretess)
-    if norm_path.is_file():
-        pretess_template = (
-            "mri_pretess {source} {label_value} {norm_path} {destination}"
-        )
-        pretess_command = pretess_template.format(
-            source=indices_mask,
-            label_value=label_value,
-            norm_path=norm_path,
-            destination=indices_mask,
-        )
-        run_shell_command(pretess_command)
+    # binarize on selected labels (creates temp indices_mask)
+    aseg = nb.load(aseg_path)
+    indices_num = [int(x) for x in indices]
+    aseg_data_bin = np.isin(aseg.get_fdata(), indices_num).astype(np.float32)
+    aseg_bin = nb.MGHImage(dataobj=aseg_data_bin, affine=aseg.affine)
+    nb.save(img=aseg_bin, filename=indices_mask)
+
+    # legacy code for applying mask smoothing
+    # from scipy import ndimage as sn
+    # k = 1.0 / np.sqrt(np.array([
+    #    [[3, 2, 3], [2, 1, 2], [3, 2, 3]],
+    #    [[2, 1, 2], [1, 1, 1], [2, 1, 1]],
+    #    [[3, 2, 3], [2, 1, 2], [3, 2, 3]],
+    # ]))
+    # aseg_data_bin = sn.convolve(aseg_data_bin, k)
+    # aseg_data_bin = np.round(aseg_data_bin / np.sum(k))
+    # nb.save(img=nb.MGHImage(dataobj=aseg_data_bin, affine=aseg.affine), \
+    #         filename=str(indices_mask).replace(".mgz", "-filter.mgz"))
+
+    # legacy code for running FreeSurfer's mri_pretess
+    # import subprocess
+    # subprocess.run(["cp", str(indices_mask), \
+    #                 str(indices_mask).replace(".mgz", "-no_pretess.mgz")])
+    ##subprocess.run(["mri_pretess", \
+    ##                   str(indices_mask).replace(".mgz", "-no_pretess.mgz"), \
+    ##                   "pretess" , \
+    ##                   str(indices_mask).replace(".mgz", "-no_pretess.mgz"), \
+    ##                   str(indices_mask)])
+    ##subprocess.run(["mri_pretess", \
+    ##                   str(indices_mask).replace(".mgz", "-no_pretess.mgz"), \
+    ##                   "pretess" , \
+    ##                   str(subject_dir / "mri/norm.mgz"), str(indices_mask)])
+    # aseg_data_bin = nb.load(indices_mask).get_fdata()
 
     # runs marching cube to extract surface
-    surface_name = f"{temp_name}.surf"
-    surface_path = destination / surface_name
-    extraction_template = "mri_mc {source} {label_value} {destination}"
-    extraction_command = extraction_template.format(
-        source=indices_mask, label_value=label_value, destination=surface_path
+    vertices, trias, _, _ = marching_cubes(
+        volume=aseg_data_bin, level=0.5, allow_degenerate=False, method="lorensen"
     )
-    run_shell_command(extraction_command)
+
+    # convert to surface RAS
+    vertices = np.matmul(
+        aseg.header.get_vox2ras_tkr(),
+        np.append(vertices, np.ones((vertices.shape[0], 1)), axis=1).transpose(),
+    ).transpose()[:, 0:3]
+
+    # create tria mesh
+    aseg_mesh = TriaMesh(v=vertices, t=trias)
+
+    # keep largest connected component
+    comps = sp.csgraph.connected_components(aseg_mesh.adj_sym, directed=False)
+    if comps[0] > 1:
+        comps_largest = np.argmax(np.unique(comps[1], return_counts=True)[1])
+        vtcs_remove = np.where(comps[1] != comps_largest)
+        tria_keep = np.sum(np.isin(aseg_mesh.t, vtcs_remove), axis=1) == 0
+        aseg_mesh.t = aseg_mesh.t[tria_keep, :]
+
+    # remove free vertices
+    aseg_mesh.rm_free_vertices_()
 
     # convert to vtk
     relative_path = "surfaces/aseg.final.{indices}.vtk".format(
         indices="_".join(indices)
     )
+
     conversion_destination = destination / relative_path
-    conversion_template = "mris_convert {source} {destination}"
-    conversion_command = conversion_template.format(
-        source=surface_path, destination=conversion_destination
-    )
-    run_shell_command(conversion_command)
+    os.makedirs(os.path.dirname(conversion_destination), exist_ok=True)
+    aseg_mesh.write_vtk(filename=conversion_destination)
 
     return conversion_destination
 
@@ -98,25 +125,21 @@ def create_aseg_surfaces(subject_dir: Path, destination: Path) -> dict[str, Path
     # Define aseg labels
 
     # combined and individual aseg labels:
-    # - Left  Striatum: left  Caudate + Putamen + Accumbens
-    # - Right Striatum: right Caudate + Putamen + Accumbens
     # - CorpusCallosum: 5 subregions combined
     # - Cerebellum: brainstem + (left+right) cerebellum WM and GM
-    # - Ventricles: (left+right) lat.vent + inf.lat.vent + choroidplexus + 3rdVent + CSF
+    # - Left-Cerebellum: left cerebellum WM and GM
+    # - Right-Cerebellum: right cerebellum WM and GM
     # - Lateral-Ventricle: lat.vent + inf.lat.vent + choroidplexus
     # - 3rd-Ventricle: 3rd-Ventricle + CSF
 
     aseg_labels = {
         "CorpusCallosum": ["251", "252", "253", "254", "255"],
         "Cerebellum": ["7", "8", "16", "46", "47"],
-        "Ventricles": ["4", "5", "14", "24", "31", "43", "44", "63"],
         "3rd-Ventricle": ["14", "24"],
         "4th-Ventricle": ["15"],
         "Brain-Stem": ["16"],
-        "Left-Striatum": ["11", "12", "26"],
         "Left-Lateral-Ventricle": ["4", "5", "31"],
-        "Left-Cerebellum-White-Matter": ["7"],
-        "Left-Cerebellum-Cortex": ["8"],
+        "Left-Cerebellum": ["7", "8"],
         "Left-Thalamus-Proper": ["10"],
         "Left-Caudate": ["11"],
         "Left-Putamen": ["12"],
@@ -125,10 +148,8 @@ def create_aseg_surfaces(subject_dir: Path, destination: Path) -> dict[str, Path
         "Left-Amygdala": ["18"],
         "Left-Accumbens-area": ["26"],
         "Left-VentralDC": ["28"],
-        "Right-Striatum": ["50", "51", "58"],
         "Right-Lateral-Ventricle": ["43", "44", "63"],
-        "Right-Cerebellum-White-Matter": ["46"],
-        "Right-Cerebellum-Cortex": ["47"],
+        "Right-Cerebellum": ["46", "47"],
         "Right-Thalamus-Proper": ["49"],
         "Right-Caudate": ["50"],
         "Right-Putamen": ["51"],
@@ -249,5 +270,5 @@ def surf_to_vtk(source: Path, destination: Path) -> Path:
     Path
         Resulting *.vtk* file.
     """
-    TriaMesh.read_fssurf(source).write_vtk(destination)
+    TriaMesh.read_fssurf(source).write_vtk(str(destination))
     return destination