Merge pull request #99 from rmarkello/map_civet

[ENH] Adds function to map CIVET to FreeSurfer

Author: rmarkello

docs/api.rst

@@ -149,6 +149,7 @@ Functions to download atlases and templates
    :toctree: generated/
 
     fetch_cammoun2012
+    fetch_civet
     fetch_conte69
     fetch_fsaverage
     fetch_pauli2018
@@ -196,6 +197,24 @@ Functions to generate (pseudo-random) datasets
    spin_data
    spin_parcels
 
+.. _ref_civet:
+
+:mod:`netneurotools.civet` - CIVET compatibility functions
+----------------------------------------------------------
+
+.. automodule:: netneurotools.civet
+   :no-members:
+   :no-inherited-members:
+
+.. currentmodule:: netneurotools.civet
+
+.. autosummary::
+   :template: function.rst
+   :toctree: generated/
+
+   read_civet
+   civet_to_freesurfer
+
 .. _ref_utils:
 
 :mod:`netneurotools.utils` - Miscellaneous, grab bag utilities

netneurotools/civet.py

@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+"""
+Functions for working with CIVET data (ugh)
+"""
+
+import nibabel as nib
+import numpy as np
+from scipy.interpolate import griddata
+
+from .datasets import fetch_civet, fetch_fsaverage
+
+_MNI305to152 = np.array([[0.9975, -0.0073, 0.0176, -0.0429],
+                         [0.0146, 1.0009, -0.0024, 1.5496],
+                         [-0.0130, -0.0093, 0.9971, 1.1840],
+                         [0.0000, 0.0000, 0.0000, 1.0000]])
+
+
+def read_civet(fname):
+    """
+    Reads a CIVET-style .obj geometry file
+
+    Parameters
+    ----------
+    fname : str or os.PathLike
+        Filepath to .obj file
+
+    Returns
+    -------
+    vertices : (N, 3)
+    triangles : (T, 3)
+    """
+
+    k, polygons = 0, []
+    with open(fname, 'r') as src:
+        n_vert = int(src.readline().split()[6])
+        vertices = np.zeros((n_vert, 3))
+        for i, line in enumerate(src):
+            if i < n_vert:
+                vertices[i] = [float(i) for i in line.split()]
+            elif i >= (2 * n_vert) + 5:
+                if not line.strip():
+                    k = 1
+                elif k == 1:
+                    polygons.extend([int(i) for i in line.split()])
+
+    triangles = np.reshape(np.asarray(polygons), (-1, 3))
+
+    return vertices, triangles
+
+
+def civet_to_freesurfer(brainmap, surface='mid', version='v1',
+                        freesurfer='fsaverage6', method='nearest',
+                        data_dir=None):
+    """
+    Projects `brainmap` in CIVET space to `freesurfer` fsaverage space
+
+    Uses a nearest-neighbor projection based on the geometry of the vertices
+
+    Parameters
+    ----------
+    brainmap : array_like
+        CIVET brainmap to be converted to freesurfer space
+    surface : {'white', 'mid'}, optional
+        Which CIVET surface to use for geometry of `brainmap`. Default: 'mid'
+    version : {'v1', 'v2'}, optional
+        Which CIVET version to use for geometry of `brainmap`. Default: 'v1'
+    freesurfer : str, optional
+        Which version of FreeSurfer space to project data to. Must be one of
+        {'fsaverage', 'fsaverage3', 'fsaverage4', 'fsaverage5', 'fsaverage6'}.
+        Default: 'fsaverage6'
+    method : {'nearest', 'linear'}, optional
+        What method of interpolation to use when projecting the data between
+        surfaces. Default: 'nearest'
+    data_dir : str, optional
+        Path to use as data directory. If not specified, will check for
+        environmental variable 'NNT_DATA'; if that is not set, will use
+        `~/nnt-data` instead. Default: None
+
+    Returns
+    -------
+    data : np.ndarray
+        Provided `brainmap` mapped to FreeSurfer
+    """
+
+    brainmap = np.asarray(brainmap)
+    densities = (81924, 327684)
+    n_vert = brainmap.shape[0]
+    if n_vert not in densities:
+        raise ValueError('Unable to interpret `brainmap` space; provided '
+                         'array must have length in {}. Received: {}'
+                         .format(densities, n_vert))
+
+    n_vert = n_vert // 2
+    icbm = fetch_civet(density='41k' if n_vert == 40962 else '164k',
+                       version=version, data_dir=data_dir, verbose=0)[surface]
+    fsavg = fetch_fsaverage(version=freesurfer, data_dir=data_dir, verbose=0)
+    fsavg = fsavg['pial' if surface == 'mid' else surface]
+
+    data = []
+    for n, hemi in enumerate(('lh', 'rh')):
+        sl = slice(n_vert * n, n_vert * (n + 1))
+        vert_cv, _ = read_civet(getattr(icbm, hemi))
+        vert_fs = nib.affines.apply_affine(
+            _MNI305to152, nib.freesurfer.read_geometry(getattr(fsavg, hemi))[0]
+        )
+        data.append(griddata(vert_cv, brainmap[sl], vert_fs, method=method))
+
+    return np.hstack(data)

netneurotools/data/osf.json

@@ -113,6 +113,33 @@
             "md5": "8f75b95c0e47ae935d10745baefa2c49"
         }
     },
+    "tpl-civet": {
+        "v1": {
+            "civet41k": {
+                "url": [
+                    "mb37e",
+                    "601daffd84ecf800fe031868"
+                ],
+                "md5": "b27219c876464992e1b61da1c60d8d6e"
+            }
+        },
+        "v2": {
+            "civet41k": {
+                "url": [
+                    "mb37e",
+                    "601dafe77ad0a80119d9483c"
+                ],
+                "md5": "a47b015e471c6a800d236f107fda5b4a"
+            },
+            "civet164k": {
+                "url": [
+                    "mb37e",
+                    "601dafe87ad0a8011ad94938"
+                ],
+                "md5": "02537ea65d5366acd8de729022a34bab"
+            }
+        }
+    },
     "ds-connectomes": {
         "celegans": {
             "url": [

netneurotools/datasets/__init__.py

@@ -6,12 +6,13 @@
     'fetch_cammoun2012', 'fetch_pauli2018', 'fetch_fsaverage', 'fetch_conte69',
     'fetch_connectome', 'available_connectomes', 'fetch_vazquez_rodriguez2019',
     'fetch_mirchi2018', 'make_correlated_xy', 'fetch_schaefer2018',
-    'fetch_hcp_standards', 'fetch_voneconomo', 'fetch_mmpall'
+    'fetch_hcp_standards', 'fetch_voneconomo', 'fetch_mmpall', 'fetch_civet'
 ]
 
 from .fetchers import (fetch_cammoun2012, fetch_pauli2018, fetch_fsaverage,
                        fetch_conte69, fetch_connectome, available_connectomes,
                        fetch_vazquez_rodriguez2019, fetch_schaefer2018,
-                       fetch_hcp_standards, fetch_voneconomo, fetch_mmpall)
+                       fetch_hcp_standards, fetch_voneconomo, fetch_mmpall,
+                       fetch_civet)
 from .generators import (make_correlated_xy)
 from .mirchi import (fetch_mirchi2018)

netneurotools/datasets/fetchers.py

@@ -16,7 +16,7 @@
 from .utils import _get_data_dir, _get_dataset_info
 from ..utils import check_fs_subjid
 
-ANNOT = namedtuple('Surface', ('lh', 'rh'))
+SURFACE = namedtuple('Surface', ('lh', 'rh'))
 
 
 def fetch_cammoun2012(version='MNI152NLin2009aSym', data_dir=None, url=None,
@@ -137,7 +137,7 @@ def fetch_cammoun2012(version='MNI152NLin2009aSym', data_dir=None, url=None,
     if version == 'MNI152NLin2009aSym':
         keys += ['info']
     elif version in ('fslr32k', 'fsaverage', 'fsaverage5', 'fsaverage6'):
-        data = [ANNOT(*data[i:i + 2]) for i in range(0, len(data), 2)]
+        data = [SURFACE(*data[i:i + 2]) for i in range(0, len(data), 2)]
     else:
         data = [data[::2][i:i + 2] for i in range(0, len(data) // 2, 2)]
         # deal with the fact that last scale is split into three files :sigh:
@@ -212,7 +212,7 @@ def fetch_conte69(data_dir=None, url=None, resume=True, verbose=1):
         data[-1] = json.load(src)
 
     # bundle hemispheres together
-    data = [ANNOT(*data[:-1][i:i + 2]) for i in range(0, 6, 2)] + [data[-1]]
+    data = [SURFACE(*data[:-1][i:i + 2]) for i in range(0, 6, 2)] + [data[-1]]
 
     return Bunch(**dict(zip(keys + ['info'], data)))
 
@@ -339,7 +339,7 @@ def fetch_fsaverage(version='fsaverage', data_dir=None, url=None, resume=True,
                             files=[(op.join(dataset_name, f), url, opts)
                                    for f in filenames])
 
-    data = [ANNOT(*data[i:i + 2]) for i in range(0, len(keys) * 2, 2)]
+    data = [SURFACE(*data[i:i + 2]) for i in range(0, len(keys) * 2, 2)]
 
     return Bunch(**dict(zip(keys, data)))
 
@@ -564,7 +564,7 @@ def fetch_schaefer2018(version='fsaverage', data_dir=None, url=None,
     data = _fetch_files(data_dir, files=files, resume=resume, verbose=verbose)
 
     if suffix == 'annot':
-        data = [ANNOT(*data[i:i + 2]) for i in range(0, len(keys) * 2, 2)]
+        data = [SURFACE(*data[i:i + 2]) for i in range(0, len(keys) * 2, 2)]
 
     return Bunch(**dict(zip(keys, data)))
 
@@ -678,7 +678,7 @@ def fetch_mmpall(version='fslr32k', data_dir=None, url=None, resume=True,
     files = [(op.join(dataset_name, version, f), url, opts) for f in filenames]
     data = _fetch_files(data_dir, files=files, resume=resume, verbose=verbose)
 
-    return ANNOT(*data)
+    return SURFACE(*data)
 
 
 def fetch_voneconomo(data_dir=None, url=None, resume=True, verbose=1):
@@ -734,6 +734,92 @@ def fetch_voneconomo(data_dir=None, url=None, resume=True, verbose=1):
     ] + ['atl-vonEconomoKoskinas_info.csv']
     files = [(op.join(dataset_name, f), url, opts) for f in filenames]
     data = _fetch_files(data_dir, files=files, resume=resume, verbose=verbose)
-    data = [ANNOT(*data[:-1:2])] + [ANNOT(*data[1:-1:2])] + [data[-1]]
+    data = [SURFACE(*data[:-1:2])] + [SURFACE(*data[1:-1:2])] + [data[-1]]
+
+    return Bunch(**dict(zip(keys, data)))
+
+
+def fetch_civet(density='41k', version='v1', data_dir=None, url=None,
+                resume=True, verbose=1):
+    """
+    Fetches CIVET surface files
+
+    Parameters
+    ----------
+    density : {'41k', '164k'}, optional
+        Which density of the CIVET-space geometry files to fetch. The
+        high-resolution '164k' surface only exists for version 'v2'
+    version : {'v1, 'v2'}, optional
+        Which version of the CIVET surfaces to use. Default: 'v2'
+    data_dir : str, optional
+        Path to use as data directory. If not specified, will check for
+        environmental variable 'NNT_DATA'; if that is not set, will use
+        `~/nnt-data` instead. Default: None
+    url : str, optional
+        URL from which to download data. Default: None
+    resume : bool, optional
+        Whether to attempt to resume partial download, if possible. Default:
+        True
+    verbose : int, optional
+        Modifies verbosity of download, where higher numbers mean more updates.
+        Default: 1
+
+    Returns
+    -------
+    filenames : :class:`sklearn.utils.Bunch`
+        Dictionary-like object with keys ['mid', 'white'] containing geometry
+        files for CIVET surface. Note for version 'v1' the 'mid' and 'white'
+        files are identical.
+
+    References
+    ----------
+    Y. Ad-Dab’bagh, O. Lyttelton, J.-S. Muehlboeck, C. Lepage, D. Einarson, K.
+    Mok, O. Ivanov, R. Vincent, J. Lerch, E. Fombonne, A. C. Evans, The CIVET
+    image-processing environment: A fully automated comprehensive pipeline for
+    anatomical neuroimaging research. Proceedings of the 12th Annual Meeting of
+    the Organization for Human Brain Mapping (2006).
+
+    Notes
+    -----
+    License: https://github.com/aces/CIVET_Full_Project/blob/master/LICENSE
+    """
+
+    densities = ['41k', '164k']
+    if density not in densities:
+        raise ValueError('The density of CIVET requested "{}" does not exist. '
+                         'Must be one of {}'.format(density, densities))
+    versions = ['v1', 'v2']
+    if version not in versions:
+        raise ValueError('The version of CIVET requested "{}" does not exist. '
+                         'Must be one of {}'.format(version, versions))
+
+    if version == 'v1' and density == '164k':
+        raise ValueError('The "164k" density CIVET surface only exists for '
+                         'version "v2"')
+
+    dataset_name = 'tpl-civet'
+    keys = ['mid', 'white']
+
+    data_dir = _get_data_dir(data_dir=data_dir)
+    info = _get_dataset_info(dataset_name)[version]['civet{}'.format(density)]
+    if url is None:
+        url = info['url']
+
+    opts = {
+        'uncompress': True,
+        'md5sum': info['md5'],
+        'move': '{}.tar.gz'.format(dataset_name)
+    }
+    filenames = [
+        op.join(dataset_name, version, 'civet{}'.format(density),
+                'tpl-civet_space-ICBM152_hemi-{}_den-{}_{}.obj'
+                .format(hemi, density, surf))
+        for surf in keys for hemi in ['L', 'R']
+    ]
+
+    data = _fetch_files(data_dir, resume=resume, verbose=verbose,
+                        files=[(f, url, opts) for f in filenames])
+
+    data = [SURFACE(*data[i:i + 2]) for i in range(0, len(keys) * 2, 2)]
 
     return Bunch(**dict(zip(keys, data)))

netneurotools/tests/test_civet.py

@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+For testing netneurotools.civet functionality
+"""
+
+import numpy as np
+import pytest
+
+from netneurotools import civet, datasets
+
+
+@pytest.fixture(scope='module')
+def civet_surf(tmp_path_factory):
+    tmpdir = str(tmp_path_factory.getbasetemp())
+    return datasets.fetch_civet(data_dir=tmpdir, verbose=0)['mid']
+
+
+def test_read_civet(civet_surf):
+    vertices, triangles = civet.read_civet(civet_surf.lh)
+    assert len(vertices) == 40962
+    assert len(triangles) == 81920
+    assert np.all(triangles.max(axis=0) < vertices.shape[0])
+
+
+def test_civet_to_freesurfer():
+    brainmap = np.random.rand(81924)
+    out = civet.civet_to_freesurfer(brainmap)
+    out2 = civet.civet_to_freesurfer(brainmap, method='linear')
+    assert out.shape[0] == out2.shape[0] == 81924
+
+    with pytest.raises(ValueError):
+        civet.civet_to_freesurfer(np.random.rand(10))

netneurotools/tests/test_datasets.py

@@ -183,6 +183,18 @@ def test_get_dataset_info(dset, expected):
         utils._get_dataset_info('notvalid')
 
 
+@pytest.mark.parametrize('version', [
+    'v1', 'v2'
+])
+def test_fetch_civet(tmpdir, version):
+    civet = datasets.fetch_civet(version=version, data_dir=tmpdir, verbose=0)
+    for key in ('mid', 'white'):
+        assert key in civet
+        for hemi in ('lh', 'rh'):
+            assert hasattr(civet[key], hemi)
+            assert os.path.isfile(getattr(civet[key], hemi))
+
+
 def test_get_data_dir(tmpdir):
     data_dir = utils._get_data_dir(tmpdir)
     assert os.path.isdir(data_dir)