Constructing 10-10 coordinates on segmented MRI scans

[1]:

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.

[2]:

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

[3]:

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’)

[4]:

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')>)

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

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

[5]:

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

Construct landmarks

[6]:

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

[7]:

lmbuilder.plot()
../../_images/examples_head_models_42_1010_system_13_0.png 
[8]:

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.69114113e+00 -2.57564888e+01  1.00475578e+02]
 [-8.25655518e+01 -1.83361855e+01 -2.66997890e+01]
 [ 8.50526962e+01 -1.57870302e+01 -1.28135719e+01]
 [-1.54080927e+00  8.74485703e+01  1.28976452e+00]
 [-3.02458239e+00  7.78819656e+01  4.12112923e+01]
 [-4.10502481e+00  5.12564507e+01  7.33338242e+01]
 [-4.62859392e+00  1.48468704e+01  9.25502930e+01]
 [-4.17725658e+00 -6.60458374e+01  9.34200821e+01]
 [-2.85223293e+00 -9.46670227e+01  6.35407906e+01]
 [-1.28595889e+00 -1.10956093e+02  2.54971142e+01]
 [ 3.18488479e-01 -1.16504837e+02 -1.52032127e+01]
 [-3.06368866e+01  8.28258667e+01 -5.99113083e+00]
 [-5.46991692e+01  6.75067520e+01 -1.29668846e+01]
 [-6.86669846e+01  4.15411911e+01 -1.84337406e+01]
 [-7.93556747e+01  1.36666317e+01 -2.32817936e+01]
 [-8.35859604e+01 -4.53174896e+01 -2.91818581e+01]
...
 [ 3.40733261e+01  7.58159790e+01  2.89205475e+01]
 [ 4.56429977e+01  7.29144363e+01  1.29440489e+01]
 [-7.77949448e+01 -5.24241524e+01  1.06570177e+01]
 [-6.71673355e+01 -5.92329941e+01  4.93465042e+01]
 [-4.14317245e+01 -6.42502289e+01  7.96890640e+01]
 [ 3.60008316e+01 -6.41902618e+01  8.72871475e+01]
 [ 6.66431046e+01 -5.88866768e+01  6.11364479e+01]
 [ 8.08320007e+01 -5.19773712e+01  2.44334316e+01]
 [-6.63746567e+01 -8.01964035e+01  2.91434383e+00]
 [-5.58452682e+01 -8.74079132e+01  3.20668449e+01]
 [-3.33801575e+01 -9.26997681e+01  5.41993027e+01]
 [ 2.88454914e+01 -9.34314880e+01  6.02286034e+01]
 [ 5.44201660e+01 -8.85375137e+01  4.20879364e+01]
 [ 6.82210312e+01 -8.12486954e+01  1.38504372e+01]
 [-4.87957840e+01 -1.02225250e+02 -6.95061302e+00]
 [-3.82177429e+01 -1.06838860e+02  9.88006401e+00]
 [-2.25627651e+01 -1.10200142e+02  2.23148594e+01]
 [ 1.79007702e+01 -1.10253891e+02  2.33871746e+01]
 [ 3.62451439e+01 -1.07865463e+02  1.51914787e+01]
 [ 5.00070915e+01 -1.03793289e+02  8.37729096e-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