"""Module to read and return anatomical information, such as:
- surfaces, with class Surf
- brains, with class BrainSurf (both hemispheres)
"""
from collections import Counter
from logging import getLogger
from os import environ
from pathlib import Path
from re import compile
from struct import unpack
from numpy import (array, empty, vstack, around, dot, append, reshape,
meshgrid, asarray)
from ..utils import MissingDependency
try:
from nibabel.freesurfer import load, read_annot
from nibabel import load as nload
except ImportError as err:
load = read_annot = nload = MissingDependency(err)
lg = getLogger(__name__)
FS_AFFINE = array([[-1, 0, 0, 128],
[0, 0, -1, 128],
[0, 1, 0, 128],
[0, 0, 0, 1]])
HEMISPHERES = 'lh', 'rh'
def _read_geometry(surf_file):
"""Read a triangular format Freesurfer surface mesh.
Parameters
----------
surf_file : str
path to surface file
Returns
-------
coords : numpy.ndarray
nvtx x 3 array of vertex (x, y, z) coordinates
faces : numpy.ndarray
nfaces x 3 array of defining mesh triangles
Notes
-----
This function comes from nibabel, but it doesn't use numpy because numpy
doesn't return the correct values in Python 3.
"""
with open(surf_file, 'rb') as f:
filebytes = f.read()
assert filebytes[:3] == b'\xff\xff\xfe'
i0 = filebytes.index(b'\x0A\x0A') + 2
i1 = i0 + 4
vnum = unpack('>i', filebytes[i0:i1])[0]
i0 = i1
i1 += 4
fnum = unpack('>i', filebytes[i0:i1])[0]
i0 = i1
i1 += 4 * vnum * 3
verts = unpack('>' + 'f' * vnum * 3, filebytes[i0:i1])
i0 = i1
i1 += 4 * fnum * 3
faces = unpack('>' + 'i' * fnum * 3, filebytes[i0:i1])
verts = asarray(verts).reshape(vnum, 3)
faces = asarray(faces).reshape(fnum, 3)
return verts, faces
def _find_neighboring_regions(pos, mri_dat, region, approx, exclude_regions):
spot_size = approx * 2 + 1
x, y, z = meshgrid(range(spot_size), range(spot_size), range(spot_size))
neighb = vstack((reshape(x, (1, spot_size ** 3)),
reshape(y, (1, spot_size ** 3)),
reshape(z, (1, spot_size ** 3)))).T - approx
regions = []
for p in range(neighb.shape[0]):
d_type = mri_dat[pos[0] + neighb[p, 0], pos[1] + neighb[p, 1],
pos[2] + neighb[p, 2]]
label_index = region['index'].index(d_type)
regions.append(region['label'][label_index])
if exclude_regions:
excluded = compile('|'.join(exclude_regions))
regions = [x for x in regions if not excluded.search(x)]
return regions
[docs]def import_freesurfer_LUT(fs_lut=None):
"""Import Look-up Table with colors and labels for anatomical regions.
It's necessary that Freesurfer is installed and that the environmental
variable 'FREESURFER_HOME' is present.
Parameters
----------
fs_lut : str or Path
path to file called FreeSurferColorLUT.txt
Returns
-------
idx : list of int
indices of regions
label : list of str
names of the brain regions
rgba : numpy.ndarray
one row is a brain region and the columns are the RGB + alpha colors
"""
if fs_lut is not None:
lg.info('Reading user-specified lookuptable {}'.format(fs_lut))
fs_lut = Path(fs_lut)
else:
try:
fs_home = environ['FREESURFER_HOME']
except KeyError:
raise OSError('Freesurfer is not installed or FREESURFER_HOME is '
'not defined as environmental variable')
else:
fs_lut = Path(fs_home) / 'FreeSurferColorLUT.txt'
lg.info('Reading lookuptable in FREESURFER_HOME {}'.format(fs_lut))
idx = []
label = []
rgba = empty((0, 4))
with fs_lut.open('r') as f:
for l in f:
if len(l) <= 1 or l[0] == '#' or l[0] == '\r':
continue
(t0, t1, t2, t3, t4, t5) = [t(s) for t, s in
zip((int, str, int, int, int, int),
l.split())]
idx.append(t0)
label.append(t1)
rgba = vstack((rgba, array([t2, t3, t4, t5])))
return idx, label, rgba
[docs]class Surf:
"""Provide class Surf, with the positions and connections of vertices.
Parameters
----------
surf_file : str or Path
freesurfer file containing the surface
Attributes
----------
surf_file : path to file
freesurfer file containing the surface
vert : numpy.ndarray
vertices of the mesh
tri : numpy.ndarray
triangulation of the mesh
"""
def __init__(self, surf_file):
self.surf_file = surf_file
surf_vert, surf_tri = _read_geometry(self.surf_file)
self.vert = surf_vert
self.tri = surf_tri
self.n_vert = surf_vert.shape[0]
[docs]class Brain:
"""Class that contains the left and right hemispheres.
Parameters
----------
freesurfer_dir : str
subject-specific directory created by freesurfer
surf_type : str
'pial', 'smoothwm', 'inflated', 'white', or 'sphere'
"""
def __init__(self, freesurfer_dir, surf_type='pial'):
freesurfer_dir = Path(freesurfer_dir)
for hemi in HEMISPHERES:
surf_file = freesurfer_dir / 'surf' / (hemi + '.' + surf_type)
setattr(self, hemi, Surf(surf_file))
[docs]class Freesurfer:
"""Provide class Freesurfer, with the information from freesurfer.
Parameters
----------
freesurfer_dir : str or Path
subject-specific directory created by freesurfer
fs_lut : str or Path
path to file called FreeSurferColorLUT.txt
Notes
-----
It's necessary that Freesurfer is installed and that the environmental
variable 'FREESURFER_HOME' is present.
"""
def __init__(self, freesurfer_dir, fs_lut=None):
freesurfer_dir = Path(freesurfer_dir)
if not freesurfer_dir.exists():
raise OSError(str(freesurfer_dir) + ' does not exist')
self.dir = freesurfer_dir
try:
lut = import_freesurfer_LUT(fs_lut)
self.lookuptable = {'index': lut[0], 'label': lut[1],
'RGBA': lut[2]}
except OSError as err:
self.lookuptable = None
lg.warning('Could not find lookup table (see below for explanation)'
'. Some functions that rely on it might complain or '
'crash.')
lg.warning(err)
@property
def surface_ras_shift(self):
"""Freesurfer uses two coordinate systems: one for volumes ("RAS") and
one for surfaces ("tkReg", "tkRAS", and "Surface RAS").
To get from surface to volume coordinates, add this numbers.
To get from volume to surface coordinates, substract this numbers.
"""
T1_path = self.dir / 'mri' / 'T1.mgz'
assert T1_path.exists()
T1 = nload(str(T1_path))
return T1.header['Pxyz_c']
[docs] def find_brain_region(self, abs_pos, parc_type='aparc', max_approx=None,
exclude_regions=None):
"""Find the name of the brain region in which an electrode is located.
Parameters
----------
abs_pos : numpy.ndarray
3x0 vector with the position of interest.
parc_type : str
'aparc', 'aparc.a2009s', 'BA', 'BA.thresh', or 'aparc.DKTatlas40'
'aparc.DKTatlas40' is only for recent freesurfer versions
max_approx : int, optional
max approximation to define position of the electrode.
exclude_regions : list of str or empty list
do not report regions if they contain these substrings. None means
that it does not exclude any region.
Notes
-----
It determines the possible brain region in which one electrode is
present, based on Freesurfer segmentation. You can imagine that
sometimes one electrode is not perfectly located within one region,
but it's a few mm away. The parameter "approx" specifies this tolerance
where each value is one mm. It keeps on searching in larger and larger
spots until it finds at least one region which is not white matter. If
there are multiple regions, it returns the region with the most
detection.
Minimal value is 0, which means only if the electrode is in the
precise location.
If you want to exclude white matter regions with 'aparc', use
exclude_regions = ('White', 'WM', 'Unknown')
and with 'aparc.a2009s', use:
exclude_regions = ('White-Matter')
"""
# convert to freesurfer coordinates of the MRI
pos = around(dot(FS_AFFINE, append(abs_pos, 1)))[:3].astype(int)
lg.debug('Position in the MRI matrix: {}'.format(pos))
mri_dat, _ = self.read_seg(parc_type)
if max_approx is None:
max_approx = 3
for approx in range(max_approx + 1):
lg.debug('Trying approx {} out of {}'.format(approx, max_approx))
regions = _find_neighboring_regions(pos, mri_dat,
self.lookuptable, approx,
exclude_regions)
if regions:
break
if regions:
c_regions = Counter(regions)
return c_regions.most_common(1)[0][0], approx
else:
return '--not found--', approx
[docs] def read_label(self, hemi, parc_type='aparc'):
"""Read the labels (annotations) for each hemisphere.
Parameters
----------
hemi : str
'lh' or 'rh'
parc_type : str
'aparc', 'aparc.a2009s', 'BA', 'BA.thresh', or 'aparc.DKTatlas40'
'aparc.DKTatlas40' is only for recent freesurfer versions
Returns
-------
numpy.ndarray
value at each vertex, indicating the label
numpy.ndarray
RGB + alpha colors for each label
list of str
names of the labels
"""
parc_file = self.dir / 'label' / (hemi + '.' + parc_type + '.annot')
vert_val, region_color, region_name = read_annot(parc_file)
region_name = [x.decode('utf-8') for x in region_name]
return vert_val, region_color, region_name
[docs] def read_seg(self, parc_type='aparc'):
"""Read the MRI segmentation.
Parameters
----------
parc_type : str
'aparc' or 'aparc.a2009s'
Returns
-------
numpy.ndarray
3d matrix with values
numpy.ndarray
4x4 affine matrix
"""
seg_file = self.dir / 'mri' / (parc_type + '+aseg.mgz')
seg_mri = load(seg_file)
seg_aff = seg_mri.affine
seg_dat = seg_mri.get_data()
return seg_dat, seg_aff
[docs] def read_brain(self, surf_type='pial'):
"""Read the surface of both hemispheres.
Parameters
----------
surf_type : str
'pial', 'smoothwm', 'inflated', 'white', or 'sphere'
Returns
-------
instance of Brain
the surfaces of both brain hemispheres
"""
return Brain(self.dir, surf_type=surf_type)