|
| 1 | +# -*- coding: utf-8 -*- |
| 2 | +""" |
| 3 | +Functions for working with CIVET data (ugh) |
| 4 | +""" |
| 5 | + |
| 6 | +import nibabel as nib |
| 7 | +import numpy as np |
| 8 | +from scipy.interpolate import griddata |
| 9 | + |
| 10 | +from .datasets import fetch_civet, fetch_fsaverage |
| 11 | + |
| 12 | +_MNI305to152 = np.array([[0.9975, -0.0073, 0.0176, -0.0429], |
| 13 | + [0.0146, 1.0009, -0.0024, 1.5496], |
| 14 | + [-0.0130, -0.0093, 0.9971, 1.1840], |
| 15 | + [0.0000, 0.0000, 0.0000, 1.0000]]) |
| 16 | + |
| 17 | + |
| 18 | +def read_civet(fname): |
| 19 | + """ |
| 20 | + Reads a CIVET-style .obj geometry file |
| 21 | +
|
| 22 | + Parameters |
| 23 | + ---------- |
| 24 | + fname : str or os.PathLike |
| 25 | + Filepath to .obj file |
| 26 | +
|
| 27 | + Returns |
| 28 | + ------- |
| 29 | + vertices : (N, 3) |
| 30 | + triangles : (T, 3) |
| 31 | + """ |
| 32 | + |
| 33 | + k, polygons = 0, [] |
| 34 | + with open(fname, 'r') as src: |
| 35 | + n_vert = int(src.readline().split()[6]) |
| 36 | + vertices = np.zeros((n_vert, 3)) |
| 37 | + for i, line in enumerate(src): |
| 38 | + if i < n_vert: |
| 39 | + vertices[i] = [float(i) for i in line.split()] |
| 40 | + elif i >= (2 * n_vert) + 5: |
| 41 | + if not line.strip(): |
| 42 | + k = 1 |
| 43 | + elif k == 1: |
| 44 | + polygons.extend([int(i) for i in line.split()]) |
| 45 | + |
| 46 | + triangles = np.reshape(np.asarray(polygons), (-1, 3)) |
| 47 | + |
| 48 | + return vertices, triangles |
| 49 | + |
| 50 | + |
| 51 | +def civet_to_freesurfer(brainmap, surface='mid', version='v1', |
| 52 | + freesurfer='fsaverage6', method='nearest', |
| 53 | + data_dir=None): |
| 54 | + """ |
| 55 | + Projects `brainmap` in CIVET space to `freesurfer` fsaverage space |
| 56 | +
|
| 57 | + Uses a nearest-neighbor projection based on the geometry of the vertices |
| 58 | +
|
| 59 | + Parameters |
| 60 | + ---------- |
| 61 | + brainmap : array_like |
| 62 | + CIVET brainmap to be converted to freesurfer space |
| 63 | + surface : {'white', 'mid'}, optional |
| 64 | + Which CIVET surface to use for geometry of `brainmap`. Default: 'mid' |
| 65 | + version : {'v1', 'v2'}, optional |
| 66 | + Which CIVET version to use for geometry of `brainmap`. Default: 'v1' |
| 67 | + freesurfer : str, optional |
| 68 | + Which version of FreeSurfer space to project data to. Must be one of |
| 69 | + {'fsaverage', 'fsaverage3', 'fsaverage4', 'fsaverage5', 'fsaverage6'}. |
| 70 | + Default: 'fsaverage6' |
| 71 | + method : {'nearest', 'linear'}, optional |
| 72 | + What method of interpolation to use when projecting the data between |
| 73 | + surfaces. Default: 'nearest' |
| 74 | + data_dir : str, optional |
| 75 | + Path to use as data directory. If not specified, will check for |
| 76 | + environmental variable 'NNT_DATA'; if that is not set, will use |
| 77 | + `~/nnt-data` instead. Default: None |
| 78 | +
|
| 79 | + Returns |
| 80 | + ------- |
| 81 | + data : np.ndarray |
| 82 | + Provided `brainmap` mapped to FreeSurfer |
| 83 | + """ |
| 84 | + |
| 85 | + brainmap = np.asarray(brainmap) |
| 86 | + densities = (81924, 327684) |
| 87 | + n_vert = brainmap.shape[0] |
| 88 | + if n_vert not in densities: |
| 89 | + raise ValueError('Unable to interpret `brainmap` space; provided ' |
| 90 | + 'array must have length in {}. Received: {}' |
| 91 | + .format(densities, n_vert)) |
| 92 | + |
| 93 | + n_vert = n_vert // 2 |
| 94 | + icbm = fetch_civet(density='41k' if n_vert == 40962 else '164k', |
| 95 | + version=version, data_dir=data_dir, verbose=0)[surface] |
| 96 | + fsavg = fetch_fsaverage(version=freesurfer, data_dir=data_dir, verbose=0) |
| 97 | + fsavg = fsavg['pial' if surface == 'mid' else surface] |
| 98 | + |
| 99 | + data = [] |
| 100 | + for n, hemi in enumerate(('lh', 'rh')): |
| 101 | + sl = slice(n_vert * n, n_vert * (n + 1)) |
| 102 | + vert_cv, _ = read_civet(getattr(icbm, hemi)) |
| 103 | + vert_fs = nib.affines.apply_affine( |
| 104 | + _MNI305to152, nib.freesurfer.read_geometry(getattr(fsavg, hemi))[0] |
| 105 | + ) |
| 106 | + data.append(griddata(vert_cv, brainmap[sl], vert_fs, method=method)) |
| 107 | + |
| 108 | + return np.hstack(data) |
0 commit comments