Note

Go to the end to download the full example code

Surface processing#

Create layer surfaces#

Uses FreeSurfer mris-inflate. Based on the number of surfaces you specify. Each will be approximately equidistant. Vertex correspondence is maintained, i.e each layer will have vertices along the vector connecting a white matter surface vertex with the corresponding pial surface vertex. Each surface will then be converted to gifti format and the vertex coordinates will adjusted by the RAS at the center of the volume to put them in native space.

import os
import numpy as np
import nibabel as nib
import k3d
import matplotlib.pyplot as plt
from lameg.viz import show_surface, rgbtoint

# Get name of each mesh that makes up the layers of the multilayer mesh
n_layers = 11
layers = np.linspace(1, 0, n_layers)

surf_path = '../test_data/sub-104/surf'
layer_fnames = []
for l, layer in enumerate(layers):
    if layer == 1:
        layer_fnames.append(os.path.join(surf_path, f'lh.pial.gii'))
    elif layer > 0 and layer < 1:
        layer_name = '{:.3f}'.format(layer)
        layer_fnames.append(os.path.join(surf_path, f'lh.{layer_name}.gii'))
    elif layer == 0:
        layer_fnames.append(os.path.join(surf_path, f'lh.white.gii'))

cam_view=[-143, 14, 31.5,
          -32, 22.5, 38.5,
          0, 0, 1]
col_r = plt.cm.cool(np.linspace(0,1, num=n_layers))

for l_idx,layer_fname in enumerate(layer_fnames):
    mesh = nib.load(layer_fname)
    plot = show_surface(mesh, camera_view=cam_view, color=col_r[l_idx,:3]*255, height=256)

Remove deep vertices#

Freesurfer operates on hemispheres independently, resulting in deep vertices and mesh faces cutting through subcortical structures. These are removed because they are not part of the cortical source space. Any vertex not labelled or labelled as unknown in the Desikan-Killiany atlas are removed. This is applied to each hemisphere of each layer mesh.

Before

surf_fname = os.path.join(surf_path, f'lh.pial.gii')

cam_view=[110, 18, 36,
          -32, 22, 42,
          0, 0, 1]
col_r = plt.cm.cool(np.linspace(0,1, num=n_layers))

mesh = nib.load(surf_fname)
plot = show_surface(mesh, camera_view=cam_view, color=col_r[0,:3]*255)

After

surf_fname = os.path.join(surf_path, f'lh.pial.nodeep.gii')

mesh = nib.load(surf_fname)
plot = show_surface(mesh, camera_view=cam_view, color=col_r[0,:3]*255)

Combine hemispheres#

The left and right hemisphere meshes are combined by concatenation of their vertices and faces (left then right). No new faces are created. This is done for each layer.

surf_fname = os.path.join(surf_path, f'pial.gii')

cam_view=[13, 64, 217,
          21, 21, 33,
          0, 1, 0]
col_r = plt.cm.cool(np.linspace(0,1, num=n_layers))

mesh = nib.load(surf_fname)
plot = show_surface(mesh, camera_view=cam_view, color=col_r[0,:3]*255)

Downsample#

The surfaces are much too dense to use in source reconstruction. They must be downsampled, but we must maintain vertex correspondence between them (see [Bonaiuto et al. 2020, Estimates of cortical column orientation improve MEG source inversion](https://doi.org/10.1016/j.neuroimage.2020.116862)). The pial surface is therefore downsampled by a factor of 10 using the vtkDecimatePro algorithm, which only removes vertices rather than creating new ones. The removed vertices are also removed from each other layer mesh, and the face structure from the pial mesh is copied to them (though the faces are not used in the source reconstruction if link vector orientations are used).

surf_fname = os.path.join(surf_path, f'pial.ds.link_vector.fixed.gii')

cam_view=[13, 64, 217,
          21, 21, 33,
          0, 1, 0]
col_r = plt.cm.cool(np.linspace(0,1, num=n_layers))

mesh = nib.load(surf_fname)
plot = show_surface(mesh, camera_view=cam_view, color=col_r[0,:3]*255)

Combine layers#

The layer meshes are then combined into a single mesh by concatenating their vertices and faces (from pial to white matter). No new faces are created (i.e. there are no edges connecting vertices across layers)

layer_fnames = []
for l, layer in enumerate(layers):
    if layer == 1:
        layer_fnames.append(os.path.join(surf_path, f'pial.ds.link_vector.fixed.gii'))
    elif layer > 0 and layer < 1:
        layer_name = '{:.3f}'.format(layer)
        layer_fnames.append(os.path.join(surf_path, f'{layer_name}.ds.link_vector.fixed.gii'))
    elif layer == 0:
        layer_fnames.append(os.path.join(surf_path, f'white.ds.link_vector.fixed.gii'))

col_r = plt.cm.cool(np.linspace(0,1, num=n_layers+1))

plot = k3d.plot(
    grid_visible=False, menu_visibility=False, camera_auto_fit=False
)

for l_idx,layer_fname in enumerate(layer_fnames):
    surface = nib.load(layer_fname)
    vertices, faces, _ = surface.agg_data()

    mesh = k3d.mesh(vertices, faces, side="double", color=rgbtoint(col_r[l_idx,:3]*255), opacity=l_idx/(n_layers-1))
    plot += mesh

plot.camera=cam_view

plot.display()

Putting it all together#

All of these steps can be run using the function:

> postprocess_freesurfer_surfaces

from lameg.surf import postprocess_freesurfer_surfaces

# Create a 2-layer surface (only pial and white)
postprocess_freesurfer_surfaces(
    'sub-104',
    '../test_data/sub-104/surf',
    'multilayer.2.ds.link_vector.fixed.gii',
    n_surfaces=2,
    ds_factor=0.1,
    orientation='link_vector',
    remove_deep=True
)

# Create an 11-layer surface
postprocess_freesurfer_surfaces(
    'sub-104',
    '../test_data/sub-104/surf',
    'multilayer.11.ds.link_vector.fixed.gii',
    n_surfaces=11,
    ds_factor=0.1,
    orientation='link_vector',
    remove_deep=True
)

# Create a 15-layer surface
postprocess_freesurfer_surfaces(
    'sub-104',
    '../test_data/sub-104/surf',
    'multilayer.15.ds.link_vector.fixed.gii',
    n_surfaces=15,
    ds_factor=0.1,
    orientation='link_vector',
    remove_deep=True
)

Total running time of the script: (0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery