@@ -0,0 +1,2255 @@
+from __future__ import division%0A%0Aimport os%0Aimport os.path as op%0A%0Aimport numpy as np%0Afrom scipy.spatial import KDTree%0A%0Aimport nibabel as nib%0Aimport nibabel.freesurfer as nifs%0Afrom nibabel.affines import apply_affine%0A%0A%0Adef vol_to_surf_xfm(vol_fname, reg_fname):%0A    %22%22%22Obtain a transformation from vol voxels -%3E Freesurfer surf coords.%0A    Parameters%0A    ----------%0A    vol_fname : string%0A        Filename pointing at image defining the vol space.%0A    reg_fname : string%0A        Filename pointing at registration file (from bbregister) that maps%0A        %60%60vol_fname%60%60 to the Freesurfer anatomy.%0A%0A    Returns%0A    -------%0A    xfm : 4 x 4 numpy array%0A        Transformation matrix that can be applied to surf coords.%0A%0A    %22%22%22%0A    # Load the Freesurfer %22tkreg%22 style transform file%0A    # Confusingly, this file actually encodes the anat-to-func transform%0A    anat2func_xfm = np.genfromtxt(reg_fname, skip_header=4, skip_footer=1)%0A    func2anat_xfm = np.linalg.inv(anat2func_xfm)%0A%0A    # Get a tkreg-compatibile mapping from IJK to RAS%0A    vol_img = nib.load(vol_fname)%0A    mgh_img = nib.MGHImage(np.zeros(vol_img.shape%5B:3%5D),%0A                           vol_img.affine,%0A                           vol_img.header)%0A    vox2ras_tkr = mgh_img.header.get_vox2ras_tkr()%0A%0A    # Combine the two transformations%0A    xfm = np.dot(func2anat_xfm, vox2ras_tkr)%0A%0A    return xfm%0A%0A%0Adef surf_to_voxel_coords(subj, hemi, xfm, surf=%22graymid%22,%0A                         subjects_dir=None):%0A    %22%22%22Obtain voxel coordinates of surface vertices in the EPI volume.%0A%0A    Parameters%0A    ----------%0A    subj : string%0A        Freesurfer subject ID.%0A    hemi : lh %7C rh%0A        Hemisphere of surface to map.%0A    xfm : 4 x 4 array%0A        Linear transformation matrix between spaces.%0A    surf : string%0A        Freesurfer surface name defining coords.%0A%0A    Returns%0A    -------%0A    i, j, k : 1d int arrays%0A        Arrays of voxel indices.%0A%0A    %22%22%22%0A    # Load the surface geometry%0A    if subjects_dir is None:%0A        subjects_dir = os.environ%5B%22SUBJECTS_DIR%22%5D%0A    surf_fname = os.path.join(subjects_dir, subj, %22surf%22,%0A                              %22%7B%7D.%7B%7D%22.format(hemi, surf))%0A    coords, _ = nib.freesurfer.read_geometry(surf_fname)%0A%0A    return apply_affine(xfm, coords).round().astype(np.int).T%0A%0A%0A%0A

Add two basic freesurfer utilities