Constructing 10-10 coordinates on segmented MRI scans

[1]:

# This cells setups the environment when executed in Google Colab.
try:
    import google.colab
    !curl -s https://raw.githubusercontent.com/ibs-lab/cedalion/dev/scripts/colab_setup.py -o colab_setup.py
    # Select branch with --branch "branch name" (default is "dev")
    %run colab_setup.py
except ImportError:
    pass

[2]:

import cedalion
import cedalion.io
import cedalion.geometry.segmentation
import cedalion.geometry.landmarks
from cedalion.imagereco.forward_model import TwoSurfaceHeadModel
import cedalion.datasets
import os.path
import pyvista

#pyvista.set_jupyter_backend("html")
pyvista.set_jupyter_backend("static")

Load segmentation masks

This example constructs the 10-10 system on the Colin27 average brain.

[3]:

SEG_DATADIR, mask_files, landmarks_file = cedalion.datasets.get_colin27_segmentation()
masks, t_ijk2ras = cedalion.io.read_segmentation_masks(SEG_DATADIR, mask_files)

Wrap the segmented head with derived surfaces in a TwoSurfaceHeadModel

[4]:

head = TwoSurfaceHeadModel.from_surfaces(
    segmentation_dir=SEG_DATADIR,
    mask_files = mask_files,
    brain_surface_file= os.path.join(SEG_DATADIR, "mask_brain.obj"),
    scalp_surface_file= os.path.join(SEG_DATADIR, "mask_scalp.obj"),
    landmarks_ras_file=landmarks_file,
    brain_face_count=None,
    scalp_face_count=None
)

Transform the scalp surface from voxel space (‘ijk’) to RAS space (‘aligned’)

[5]:

scalp_surface = head.scalp
display(scalp_surface)
scalp_surface = scalp_surface.apply_transform(t_ijk2ras)
display(scalp_surface)

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(10050, 3), faces.shape=(20096, 3))>, crs='ijk', units=<Unit('dimensionless')>, vertex_coords={})

TrimeshSurface(mesh=<trimesh.Trimesh(vertices.shape=(10050, 3), faces.shape=(20096, 3))>, crs='aligned', units=<Unit('millimeter')>, vertex_coords={})

Transform initial landmarks from voxel space (‘ijk’) to RAS space (‘aligned’)

[6]:

landmarks_ras = head.landmarks.points.apply_transform(t_ijk2ras)

Construct landmarks

[7]:

lmbuilder = cedalion.geometry.landmarks.LandmarksBuilder1010(scalp_surface, landmarks_ras)
all_landmarks = lmbuilder.build()

/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:242: UserWarning: WIP: distance calculation around ears
  warnings.warn("WIP: distance calculation around ears")

Visualize

[8]:

lmbuilder.plot()

../../_images/examples_head_models_42_1010_system_14_0.png

[9]:

display(all_landmarks)

<xarray.DataArray (label: 73, aligned: 3)> Size: 2kB
<Quantity([[-4.99239750e-02  7.97894669e+01 -3.61411209e+01]
 [ 1.75606728e+00 -1.00949997e+02 -5.48579865e+01]
 [-7.19499969e+01 -1.61194057e+01 -5.46569786e+01]
 [ 7.59499969e+01 -1.30380936e+01 -5.40778427e+01]
 [-4.50463189e+00 -2.57601506e+01  1.00486147e+02]
 [-8.25398134e+01 -1.83963801e+01 -2.58655098e+01]
 [ 8.50524916e+01 -1.57862304e+01 -1.28503999e+01]
 [-1.49455525e+00  8.74494123e+01  1.29367116e+00]
 [-2.92668956e+00  7.78851914e+01  4.12196329e+01]
 [-3.96320511e+00  5.12599604e+01  7.33488565e+01]
 [-4.45736100e+00  1.48428770e+01  9.25585733e+01]
 [-3.99464688e+00 -6.60469443e+01  9.34368697e+01]
 [-2.70350958e+00 -9.46608241e+01  6.35509013e+01]
 [-1.18320783e+00 -1.10959783e+02  2.55067853e+01]
 [ 3.70594105e-01 -1.16504924e+02 -1.51988858e+01]
 [-3.06872504e+01  8.27403384e+01 -5.71760589e+00]
 [-5.47414793e+01  6.72495662e+01 -1.24577241e+01]
 [-6.88353407e+01  4.12741251e+01 -1.78084721e+01]
 [-7.92158628e+01  1.37193703e+01 -2.24480436e+01]
 [-8.34870294e+01 -4.53348542e+01 -2.83145550e+01]
...
 [ 3.40164052e+01  7.58534693e+01  2.88817762e+01]
 [ 4.56434960e+01  7.29187153e+01  1.29169121e+01]
 [-7.76712839e+01 -5.24810720e+01  1.16422583e+01]
 [-6.69674751e+01 -5.92159293e+01  4.97946055e+01]
 [-4.10877778e+01 -6.42529363e+01  8.00248409e+01]
 [ 3.64243104e+01 -6.41920426e+01  8.70440208e+01]
 [ 6.66232577e+01 -5.89685026e+01  6.11064040e+01]
 [ 8.08032929e+01 -5.20668567e+01  2.43955413e+01]
 [-6.62975121e+01 -8.02922082e+01  3.58543120e+00]
 [-5.55636584e+01 -8.74852436e+01  3.26403019e+01]
 [-3.28938756e+01 -9.27391912e+01  5.45402688e+01]
 [ 2.88678402e+01 -9.34487937e+01  6.01805054e+01]
 [ 5.47290596e+01 -8.84850186e+01  4.16505668e+01]
 [ 6.81948269e+01 -8.13164070e+01  1.37985355e+01]
 [-4.88957402e+01 -1.02145995e+02 -6.45460392e+00]
 [-3.81615526e+01 -1.06806089e+02  1.02908121e+01]
 [-2.23346802e+01 -1.10185812e+02  2.25205415e+01]
 [ 1.81919090e+01 -1.10274054e+02  2.32645868e+01]
 [ 3.64615246e+01 -1.07861961e+02  1.48493294e+01]
 [ 5.01652160e+01 -1.03775222e+02  4.08961569e-01]], 'millimeter')>
Coordinates:
  * label    (label) <U3 876B 'Nz' 'Iz' 'LPA' 'RPA' ... 'PO1' 'PO2' 'PO4' 'PO6'
    type     (label) object 584B PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: aligned