add laterality index, functional connectivity, surface-level analyses support, and other bug fixes

funROI/__init__.py

@@ -2,51 +2,40 @@
 funROI: A package for functional region of interest analysis in fMRI data.
 """
 
-from .settings import Settings
-
-_settings = Settings()
-
-set_bids_data_folder = _settings.set_bids_data_folder
-"""Set the BIDS data folder path."""
-
-get_bids_data_folder = _settings.get_bids_data_folder
-"""Get the BIDS data folder path."""
-
-set_bids_deriv_folder = _settings.set_bids_deriv_folder
-"""Set the BIDS derivatives folder path."""
-
-get_bids_deriv_folder = _settings.get_bids_deriv_folder
-"""Get the BIDS derivatives folder path."""
-
-set_bids_preprocessed_folder = _settings.set_bids_preprocessed_folder
-"""Set the BIDS preprocessed folder path."""
-
-get_bids_preprocessed_folder = _settings.get_bids_preprocessed_folder
-"""Get the BIDS preprocessed folder path."""
-
-get_bids_preprocessed_folder_relative = (
-    _settings.get_bids_preprocessed_folder_relative
+from . import analysis, datasets, first_level
+from .analysis import (
+    EffectEstimator,
+    FunctionalConnectivityEstimator,
+    FROIGenerator,
+    LateralityIndexAnalyzer,
+    load_preprocessed_bold_for_fc,
+    OverlapEstimator,
+    ParcelsGenerator,
+    preprocess_bold_for_fc,
+    SurfaceParcelsGenerator,
+    SpatialCorrelationEstimator,
 )
-"""Get the BIDS preprocessed folder path relative to the data folder."""
-
-set_analysis_output_folder = _settings.set_analysis_output_folder
-"""Set the analysis output folder path."""
-
-get_analysis_output_folder = _settings.get_analysis_output_folder
-"""Get the analysis output folder path."""
-
-reset_settings = _settings.reset
-"""Reset all path settings to None."""
-
-from .first_level import *
-from .analysis import *
-from .parcels import ParcelsConfig
 from .froi import FROIConfig
+from .parcels import ParcelsConfig, SurfaceParcelsConfig
+from .settings import (
+    Settings,
+    get_analysis_output_folder,
+    get_bids_data_folder,
+    get_bids_deriv_folder,
+    get_bids_preprocessed_folder,
+    get_bids_preprocessed_folder_relative,
+    reset_settings,
+    set_analysis_output_folder,
+    set_bids_data_folder,
+    set_bids_deriv_folder,
+    set_bids_preprocessed_folder,
+)
 
 __all__ = [
-    "first_level",
+    "Settings",
     "analysis",
-    "set_bids_data_folder",
+    "datasets",
+    "first_level",
     "set_bids_data_folder",
     "get_bids_data_folder",
     "set_bids_deriv_folder",
@@ -57,7 +46,17 @@
     "set_analysis_output_folder",
     "get_analysis_output_folder",
     "reset_settings",
-    "FROIConfig",
     "ParcelsConfig",
-    "datasets",
+    "SurfaceParcelsConfig",
+    "FROIConfig",
+    "ParcelsGenerator",
+    "SurfaceParcelsGenerator",
+    "FROIGenerator",
+    "EffectEstimator",
+    "FunctionalConnectivityEstimator",
+    "preprocess_bold_for_fc",
+    "load_preprocessed_bold_for_fc",
+    "SpatialCorrelationEstimator",
+    "OverlapEstimator",
+    "LateralityIndexAnalyzer",
 ]

funROI/_registry.py

@@ -0,0 +1,88 @@
+from pathlib import Path
+from typing import Any, Mapping
+import json
+
+import pandas as pd
+
+
+def _normalize_record_value(value: Any) -> Any:
+    if isinstance(value, Path):
+        return str(value)
+    if isinstance(value, Mapping):
+        return json.dumps(
+            {
+                str(key): _normalize_record_value(val)
+                for key, val in value.items()
+            },
+            sort_keys=True,
+        )
+    if isinstance(value, (list, tuple)):
+        return json.dumps(
+            [_normalize_record_value(item) for item in value]
+        )
+    return value
+
+
+def _row_matches(row: pd.Series, criteria: Mapping[str, Any]) -> bool:
+    for column, expected in criteria.items():
+        expected = _normalize_record_value(expected)
+        actual = row[column]
+        if expected is None:
+            if not pd.isna(actual):
+                return False
+        elif actual != expected:
+            return False
+    return True
+
+
+def build_record_frame(criteria: Mapping[str, Any]) -> pd.DataFrame:
+    return pd.DataFrame(
+        [{key: _normalize_record_value(value) for key, value in criteria.items()}]
+    )
+
+
+def get_or_create_record_id(
+    info_path: Path, criteria: Mapping[str, Any], create: bool = False
+) -> int:
+    record = build_record_frame(criteria)
+    if not info_path.exists():
+        record_id = 0
+        if create:
+            info_path.parent.mkdir(parents=True, exist_ok=True)
+            record.assign(id=record_id).to_csv(info_path, index=False)
+        return record_id
+
+    info = pd.read_csv(info_path)
+    matches = info.apply(lambda row: _row_matches(row, criteria), axis=1)
+    if matches.any():
+        return int(info.loc[matches, "id"].iloc[0])
+
+    record_id = 0 if info.empty else int(info["id"].max()) + 1
+    if create:
+        info = pd.concat(
+            [info, record.assign(id=record_id)], ignore_index=True
+        )
+        info.to_csv(info_path, index=False)
+    return record_id
+
+
+def append_record(info_path: Path, new_row: pd.DataFrame) -> int:
+    info_path.parent.mkdir(parents=True, exist_ok=True)
+    new_row = new_row.copy()
+
+    if not info_path.exists():
+        record_id = 0
+        info = new_row.assign(id=record_id)
+    else:
+        info = pd.read_csv(info_path)
+        record_id = 0 if info.empty else int(info["id"].max()) + 1
+        info = pd.concat(
+            [info, new_row.assign(id=record_id)], ignore_index=True
+        )
+
+    info.to_csv(info_path, index=False)
+    return record_id
+
+
+def get_record_folder(root: Path, prefix: str, record_id: int) -> Path:
+    return root / f"{prefix}_{record_id:04d}"

funROI/_surface.py

@@ -0,0 +1,119 @@
+from pathlib import Path
+from typing import Dict, Mapping, Union
+
+import nibabel as nib
+import numpy as np
+from nilearn.surface import SurfaceImage, load_surf_data, load_surf_mesh
+
+
+SURFACE_HEMIS = ("L", "R")
+SURFACE_PARTS = {"L": "left", "R": "right"}
+
+
+def _as_path_dict(
+    paths: Mapping[str, Union[str, Path]]
+) -> Dict[str, Path]:
+    return {hemi: Path(paths[hemi]) for hemi in SURFACE_HEMIS}
+
+
+def is_surface_image(img) -> bool:
+    return isinstance(img, SurfaceImage)
+
+
+def write_gifti(path: Union[str, Path], data: np.ndarray) -> None:
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    data = np.asarray(data, dtype=np.float32)
+    if data.ndim == 1:
+        arrays = [data]
+    else:
+        arrays = [data[:, i] for i in range(data.shape[1])]
+    img = nib.gifti.GiftiImage(
+        darrays=[nib.gifti.GiftiDataArray(data=array) for array in arrays]
+    )
+    nib.save(img, path)
+
+
+def load_surface_numeric_data(path: Union[str, Path]) -> np.ndarray:
+    return np.asarray(load_surf_data(path), dtype=np.float32)
+
+
+def load_surface_image(
+    data_paths: Mapping[str, Union[str, Path]],
+    mesh_paths: Mapping[str, Union[str, Path]],
+) -> SurfaceImage:
+    data_paths = _as_path_dict(data_paths)
+    mesh_paths = _as_path_dict(mesh_paths)
+    return SurfaceImage(
+        mesh={
+            SURFACE_PARTS[hemi]: load_surf_mesh(mesh_paths[hemi])
+            for hemi in SURFACE_HEMIS
+        },
+        data={
+            SURFACE_PARTS[hemi]: load_surface_numeric_data(data_paths[hemi])
+            for hemi in SURFACE_HEMIS
+        },
+    )
+
+
+def get_surface_data_parts(img: SurfaceImage) -> Dict[str, np.ndarray]:
+    return {
+        hemi: np.asarray(img.data.parts[SURFACE_PARTS[hemi]])
+        for hemi in SURFACE_HEMIS
+    }
+
+
+def flatten_surface_parts(parts: Mapping[str, np.ndarray]) -> np.ndarray:
+    return np.concatenate(
+        [np.asarray(parts[hemi]).reshape(-1) for hemi in SURFACE_HEMIS]
+    )
+
+
+def flatten_image_data(img) -> np.ndarray:
+    if is_surface_image(img):
+        return flatten_surface_parts(get_surface_data_parts(img))
+    return np.asarray(img.get_fdata()).reshape(-1)
+
+
+def surface_hemi_sizes(img: SurfaceImage) -> Dict[str, int]:
+    parts = get_surface_data_parts(img)
+    return {hemi: int(parts[hemi].shape[0]) for hemi in SURFACE_HEMIS}
+
+
+def surface_flat_to_parts(
+    flat_data: np.ndarray, reference_img: SurfaceImage
+) -> Dict[str, np.ndarray]:
+    flat_data = np.asarray(flat_data).reshape(-1)
+    sizes = surface_hemi_sizes(reference_img)
+    parts = {}
+    start = 0
+    for hemi in SURFACE_HEMIS:
+        stop = start + sizes[hemi]
+        parts[hemi] = flat_data[start:stop].astype(np.float32, copy=False)
+        start = stop
+    if start != flat_data.size:
+        raise ValueError(
+            "Surface data length does not match the reference mesh size."
+        )
+    return parts
+
+
+def surface_image_from_flat(
+    flat_data: np.ndarray, reference_img: SurfaceImage
+) -> SurfaceImage:
+    parts = surface_flat_to_parts(flat_data, reference_img)
+    return SurfaceImage(
+        mesh=reference_img.mesh,
+        data={
+            SURFACE_PARTS[hemi]: parts[hemi] for hemi in SURFACE_HEMIS
+        },
+    )
+
+
+def save_surface_image(
+    img: SurfaceImage, data_paths: Mapping[str, Union[str, Path]]
+) -> None:
+    data_paths = _as_path_dict(data_paths)
+    parts = get_surface_data_parts(img)
+    for hemi in SURFACE_HEMIS:
+        write_gifti(data_paths[hemi], parts[hemi])

funROI/analysis/__init__.py

@@ -1,14 +1,24 @@
 from .parcels_gen import ParcelsGenerator
+from .surface_parcels_gen import SurfaceParcelsGenerator
 from .froi_gen import FROIGenerator
 from .effect import EffectEstimator
+from .fconn import (
+    FunctionalConnectivityEstimator,
+    load_preprocessed_bold_for_fc,
+    preprocess_bold_for_fc,
+)
 from .spcorr import SpatialCorrelationEstimator
 from .overlap import OverlapEstimator
 from .li import LateralityIndexAnalyzer
 
 __all__ = [
     "ParcelsGenerator",
+    "SurfaceParcelsGenerator",
     "FROIGenerator",
     "EffectEstimator",
+    "FunctionalConnectivityEstimator",
+    "preprocess_bold_for_fc",
+    "load_preprocessed_bold_for_fc",
     "SpatialCorrelationEstimator",
     "OverlapEstimator",
     "LateralityIndexAnalyzer",