Head model fiducials and landmarks

Cedalion ships with segmentations as well as brain and scalp surfaces for the Colin27 and ICBM-152 heads.

This notebook documents the source of the fiducial landmarks and compares the outputs of the landmark builder which we distribute together with the head models.

[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.dataclasses as cdc
import cedalion.data
import cedalion.geometry.landmarks
import cedalion.geometry.segmentation
import cedalion.io
import cedalion.vis.blocks as vbx
import cedalion.dot
import cedalion.xrutils as xrutils
import matplotlib.pyplot as p
import numpy as np
import pyvista
import xarray as xr

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

Utility function

[3]:
def compare_landmarks(scalp, landmarks, reference):
    fiducial_labels = ["Nz", "Iz", "LPA", "RPA", "Cz"]

    plt = pyvista.Plotter()
    vbx.plot_surface(plt, scalp, opacity=0.7)
    vbx.plot_labeled_points(plt, reference, color="y")
    vbx.plot_labeled_points(plt, landmarks, color="g")
    vbx.plot_labeled_points(plt, landmarks.loc[fiducial_labels], color="r")
    plt.show()

    common_labels = list(set(reference.label.values) & set(landmarks.label.values))

    diffs = xrutils.norm(
        landmarks.sel(label=landmarks.label.isin(common_labels))
        - reference.sel(label=reference.label.isin(common_labels)),
        dim="mni",
    )
    diffs = diffs.pint.to("mm").pint.dequantify()
    median_diff = np.median(diffs).item()

    p.figure()
    p.hist(diffs, np.arange(0, 10.5, 0.5))
    p.axvline(median_diff, c="r", ls=":")
    p.xlabel(r"$|\vec{x}_{ours} - \vec{x}_{ref}|$ / mm")
    p.ylabel(r"# landmarks")
    print(f"The median distance between landmarks is {median_diff:.1f} mm.")

Colin 27

Load Colin27 head model

[4]:
colin_ijk = cedalion.dot.get_standard_headmodel("colin27")
colin_ras = colin_ijk.apply_transform(colin_ijk.t_ijk2ras)

Load fieldtrip 10-5 coordinates:

[5]:
colin_fieldtrip = cedalion.data.get_fieldtrip_colin27_landmarks()
colin_fieldtrip = colin_fieldtrip.points.set_crs("mni")
colin_fieldtrip
Downloading file 'fieldtrip_standard1005.elc' from 'https://raw.githubusercontent.com/fieldtrip/fieldtrip/refs/heads/master/template/electrode/standard_1005.elc' to '/home/runner/.cache/cedalion/dev'.
[5]:
<xarray.DataArray (label: 346, mni: 3)> Size: 8kB
<Quantity([[-8.60761e+01 -1.99897e+01 -4.79860e+01]
 [ 8.57939e+01 -2.00093e+01 -4.80310e+01]
 [ 8.30000e-03  8.68110e+01 -3.99830e+01]
 ...
 [ 8.57939e+01 -4.50093e+01 -6.80310e+01]
 [-8.60761e+01 -2.49897e+01 -6.79860e+01]
 [ 8.57939e+01 -2.50093e+01 -6.80310e+01]], 'millimeter')>
Coordinates:
  * label    (label) <U6 8kB 'LPA' 'RPA' 'Nz' 'Fp1' ... 'M1' 'M2' 'A1' 'A2'
    type     (label) object 3kB PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: mni

Visualize fieldtrip coordinates and the Colin27 scalp surface that ships with Cedalion. Note tha Nz floats above the surface.

[6]:
plt = pyvista.Plotter()
vbx.plot_surface(plt, colin_ras.scalp, opacity=.7)
vbx.plot_labeled_points(plt, colin_fieldtrip, color="y")
plt.show()
../../_images/examples_head_models_48_headmodel_landmarks_verification_11_0.png

Select the five fiducicals, with Nz snapped to the scalp.

[7]:
colin_fiducials = xr.concat(
    (
        colin_fieldtrip.loc[["Iz", "LPA", "RPA", "Cz"]],
        colin_ras.scalp.snap(colin_fieldtrip.loc[["Nz"]]),
    ),
    dim="label",
)
colin_fiducials
[7]:
<xarray.DataArray (label: 5, mni: 3)> Size: 120B
<Quantity([[ 4.5000000e-03 -1.1856500e+02 -2.3078000e+01]
 [-8.6076100e+01 -1.9989700e+01 -4.7986000e+01]
 [ 8.5793900e+01 -2.0009300e+01 -4.8031000e+01]
 [ 4.0090000e-01 -9.1670000e+00  1.0024400e+02]
 [-6.7919500e-01  8.3259347e+01 -3.9408527e+01]], 'millimeter')>
Coordinates:
  * label    (label) <U6 120B 'Iz' 'LPA' 'RPA' 'Cz' 'Nz'
    type     (label) object 40B PointType.LANDMARK ... PointType.LANDMARK
Dimensions without coordinates: mni

Construct remaining landmarks:

[8]:
lmbuilder = cedalion.geometry.landmarks.LandmarksBuilder1010(colin_ras.scalp, colin_fiducials)
colin_all_landmarks = lmbuilder.build()
lmbuilder.plot()
/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:319: UserWarning: WIP: distance calculation around ears
  warnings.warn("WIP: distance calculation around ears")
../../_images/examples_head_models_48_headmodel_landmarks_verification_15_1.png

Compare constructed landmarks to fieldtrip reference

[9]:
compare_landmarks(colin_ras.scalp, colin_all_landmarks, colin_fieldtrip)
../../_images/examples_head_models_48_headmodel_landmarks_verification_17_0.png
The median distance between landmarks is 1.0 mm.
../../_images/examples_head_models_48_headmodel_landmarks_verification_17_2.png

ICBM-152

The fiducials and reference landmarks for the ICBM-152 head were adopted from Cutini et al.

Translate data from the paper to labeled points:

[10]:
def get_cutini_icbm_landmarks():
    # from the supplementary table of Cutini et al., 2010
    # doi: 10.1016/j.neuroimage.2010.09.030

    cutini_icbm_coords = """
    Fp1     -28     83      -5
    Fp2     31      84      -5
    Fz      0       58      66
    F3      -50     55      39
    F4      50      55      41
    F7      -68     41      -12
    F8      69      41      -12
    Cz      0       -9      100
    C3      -66     -11     63
    C4      69      -11     63
    T7      -80     -16     -10
    T8      81      -17     -10
    Pz      0       -82     78
    P3      -52     -80     54
    P4      53      -80     54
    P7      -71     -72     -5
    P8      70      -72     -5
    O1      -26     -111    3
    O2      27      -109    3
    Fpz     0       88      -1
    AFz     0       79      37
    AF3     -40     72      24
    AF4     40      73      24
    AF7     -53     66      -9
    AF8     54      65      -9
    F1      -30     60      56
    F2      29      59      59
    F5      -63     47      13
    F6      64      47      15
    FCz     0       26      87
    FC1     -36     28      76
    FC2     35      27      77
    FC3     -60     25      52
    FC4     62      24      54
    FC5     -74     19      18
    FC6     75      20      21
    FT7     -75     13      -14
    FT8     78      13      -15
    C1      -38     -8      89
    C2      38      -10     89
    C5      -80     -15     26
    C6      81      -14     26
    CPz     0       -49     99
    CP1     -37     -48     91
    CP2     37      -51     90
    CP3     -65     -50     66
    CP4     65      -51     64
    CP5     -79     -47     30
    CP6     78      -48     29
    TP7     -80     -45     -8
    TP8     80      -45     -9
    P1      -28     -83     74
    P2      30      -83     73
    P5      -66     -79     26
    P6      65      -77     26
    POz     0       -103    46
    PO3     -37     -101    37
    PO4     39      -100    37
    PO7     -54     -97     0
    PO8     52      -97     0
    Oz      0       -115    10
    """

    labels = []
    coords = []
    for line in cutini_icbm_coords.splitlines():
        line = line.strip()
        if not line:
            continue

        fields = line.split()
        labels.append(fields[0])
        coords.append(np.asarray([float(i) for i in fields[1:]]))

    return cdc.build_labeled_points(
        np.asarray(coords), crs="mni", units="mm", labels=labels
    )


def get_cutini_fiducials():
    # from Cutini et al., 2010, Table 1
    # doi: 10.1016/j.neuroimage.2010.09.030
    return cdc.build_labeled_points(
        np.asarray(
            [
                [0.0, 84.0, -43.0],
                [0.0, -114.0, -30.0],
                [-75.09, -19.49, -47.98],
                [76, -19.45, -47.7],
                [0.0, -9.0, 100.0],
            ]
        ),
        crs="mni",
        units="mm",
        labels=["Nz", "Iz", "LPA", "RPA", "Cz"],
    )

[11]:
icbm_ijk = cedalion.dot.get_standard_headmodel("icbm152")
icbm_ras = icbm_ijk.apply_transform(icbm_ijk.t_ijk2ras)
[12]:
plt = pyvista.Plotter()
vbx.plot_surface(plt, icbm_ras.scalp, opacity=.7)
vbx.plot_labeled_points(plt, get_cutini_icbm_landmarks(), color="y")
vbx.plot_labeled_points(plt, get_cutini_fiducials(), color="r")
plt.show()
../../_images/examples_head_models_48_headmodel_landmarks_verification_21_0.png
[13]:
icbm_fiducials = icbm_ras.scalp.snap(get_cutini_fiducials())
icbm_fiducials = get_cutini_fiducials()
[14]:
lmbuilder = cedalion.geometry.landmarks.LandmarksBuilder1010(icbm_ras.scalp, icbm_fiducials)
icbm_all_landmarks = lmbuilder.build()
lmbuilder.plot()
/home/runner/work/cedalion/cedalion/src/cedalion/geometry/landmarks.py:319: UserWarning: WIP: distance calculation around ears
  warnings.warn("WIP: distance calculation around ears")
../../_images/examples_head_models_48_headmodel_landmarks_verification_23_1.png
[15]:
compare_landmarks(icbm_ras.scalp, icbm_all_landmarks, get_cutini_icbm_landmarks())
../../_images/examples_head_models_48_headmodel_landmarks_verification_24_0.png
The median distance between landmarks is 2.7 mm.
../../_images/examples_head_models_48_headmodel_landmarks_verification_24_2.png

Export landmarks

[16]:
#cedalion.io.probe_geometry.save_mrk_json("colin.mrk.json", icbm_all_landmarks, "mni")
#cedalion.io.probe_geometry.save_mrk_json("icbm.mrk.json", icbm_all_landmarks, "mni")

References

[17]:
cedalion.bib.dump_to_notebook()

Methods used

[1]Holmes1998cedalion.data.get_colin27_headmodel_filesColin J. Holmes, Rick Hoge, Louis Collins, Roger Woods, Arthur W. Toga, and Alan C. Evans. Enhancement of mr images using registration for signal averaging. Journal of Computer Assisted Tomography, 22(2):324–333, March 1998. doi:10.1097/00004728-199803000-00032.
[2]Fischl2012cedalion.data.get_colin27_headmodel_files, cedalion.data.get_icbm152_headmodel_filesBruce Fischl. FreeSurfer. NeuroImage, 62(2):774–781, 2012. doi:10.1016/j.neuroimage.2012.01.021.
[3]Schaefer2018cedalion.data.get_colin27_headmodel_files, cedalion.data.get_icbm152_headmodel_filesAlexander Schaefer, Ru Kong, Evan M Gordon, Timothy O Laumann, Xi-Nian Zuo, Avram J Holmes, Simon B Eickhoff, and BT Thomas Yeo. Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri. Cerebral cortex, 28(9):3095–3114, 2018. doi:10.1093/cercor/bhx179.
[4]Oostenveld2001cedalion.geometry.landmarks.LandmarksBuilder1010.__init__Robert Oostenveld and Peter Praamstra. The five percent electrode system for high-resolution eeg and erp measurements. Clinical Neurophysiology, 112(4):713–719, 2001. doi:https://doi.org/10.1016/S1388-2457(00)00527-7.
[5]Fonov2011cedalion.data.get_icbm152_headmodel_filesVladimir Fonov, Alan C. Evans, Kelly Botteron, C. Robert Almli, Robert C. McKinstry, and D. Louis Collins. Unbiased average age-appropriate atlases for pediatric studies. NeuroImage, 54(1):313–327, January 2011. doi:10.1016/j.neuroimage.2010.07.033.