test_coreg.py

# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

import os
from functools import reduce
from glob import glob
from shutil import copyfile

import numpy as np
import pytest
from numpy.testing import (
    assert_allclose,
    assert_array_almost_equal,
    assert_array_equal,
    assert_array_less,
)

import mne
from mne._fiff.constants import FIFF
from mne.channels import DigMontage
from mne.coreg import (
    Coregistration,
    _is_mri_subject,
    coregister_fiducials,
    create_default_subject,
    fit_matched_points,
    get_mni_fiducials,
    scale_labels,
    scale_mri,
    scale_source_space,
)
from mne.datasets import testing
from mne.io import read_fiducials, read_info
from mne.source_space import write_source_spaces
from mne.transforms import (
    Transform,
    _angle_between_quats,
    apply_trans,
    invert_transform,
    read_trans,
    rot_to_quat,
    rotation,
    scaling,
    translation,
)
from mne.utils import catch_logging

data_path = testing.data_path(download=False)
subjects_dir = data_path / "subjects"
fid_fname = subjects_dir / "sample" / "bem" / "sample-fiducials.fif"
raw_fname = data_path / "MEG" / "sample" / "sample_audvis_trunc_raw.fif"
trans_fname = data_path / "MEG" / "sample" / "sample_audvis_trunc-trans.fif"


@pytest.fixture
def few_surfaces(monkeypatch):
    """Set the _MNE_FEW_SURFACES env var."""
    monkeypatch.setenv("_MNE_FEW_SURFACES", "true")
    yield


def test_coregister_fiducials():
    """Test coreg.coregister_fiducials()."""
    # prepare head and MRI fiducials
    trans = Transform(
        "head", "mri", rotation(0.4, 0.1, 0).dot(translation(0.1, -0.1, 0.1))
    )
    coords_orig = np.array(
        [
            [-0.08061612, -0.02908875, -0.04131077],
            [0.00146763, 0.08506715, -0.03483611],
            [0.08436285, -0.02850276, -0.04127743],
        ]
    )
    coords_trans = apply_trans(trans, coords_orig)

    def make_dig(coords, cf):
        return (
            {"coord_frame": cf, "ident": 1, "kind": 1, "r": coords[0]},
            {"coord_frame": cf, "ident": 2, "kind": 1, "r": coords[1]},
            {"coord_frame": cf, "ident": 3, "kind": 1, "r": coords[2]},
        )

    mri_fiducials = make_dig(coords_trans, FIFF.FIFFV_COORD_MRI)
    info = {"dig": make_dig(coords_orig, FIFF.FIFFV_COORD_HEAD)}

    # test coregister_fiducials()
    trans_est = coregister_fiducials(info, mri_fiducials)
    assert trans_est.from_str == trans.from_str
    assert trans_est.to_str == trans.to_str
    assert_array_almost_equal(trans_est["trans"], trans["trans"])


@pytest.mark.slowtest  # can take forever on OSX Travis
@testing.requires_testing_data
@pytest.mark.parametrize("scale", (0.9, [1, 0.2, 0.8]))
def test_scale_mri(tmp_path, few_surfaces, scale):
    """Test creating fsaverage and scaling it."""
    pytest.importorskip("nibabel")
    # create fsaverage using the testing "fsaverage" instead of the FreeSurfer
    # one
    fake_home = data_path
    create_default_subject(subjects_dir=tmp_path, fs_home=fake_home, verbose=True)
    assert _is_mri_subject("fsaverage", tmp_path), "Creating fsaverage failed"

    fid_path = tmp_path / "fsaverage" / "bem" / "fsaverage-fiducials.fif"
    os.remove(fid_path)
    create_default_subject(update=True, subjects_dir=tmp_path, fs_home=fake_home)
    assert fid_path.exists(), "Updating fsaverage"

    # copy MRI file from sample data (shouldn't matter that it's incorrect,
    # so here choose a small one)
    path_from = fake_home / "subjects" / "sample" / "mri" / "T1.mgz"
    path_to = tmp_path / "fsaverage" / "mri" / "orig.mgz"
    copyfile(path_from, path_to)

    # remove redundant label files
    label_temp = tmp_path / "fsaverage" / "label" / "*.label"
    label_paths = glob(str(label_temp))
    for label_path in label_paths[1:]:
        os.remove(label_path)

    # create source space
    bem_path = tmp_path / "fsaverage" / "bem"
    bem_fname = "fsaverage-%s-src.fif"
    src = mne.setup_source_space(
        "fsaverage", "ico0", subjects_dir=tmp_path, add_dist=False
    )
    mri = tmp_path / "fsaverage" / "mri" / "orig.mgz"
    vsrc = mne.setup_volume_source_space(
        "fsaverage",
        pos=50,
        mri=mri,
        subjects_dir=tmp_path,
        add_interpolator=False,
    )
    write_source_spaces(bem_path / (bem_fname % "vol-50"), vsrc)

    # scale fsaverage
    write_source_spaces(bem_path / (bem_fname % "ico-0"), src, overwrite=True)
    scale_mri(
        "fsaverage",
        "flachkopf",
        scale,
        True,
        subjects_dir=tmp_path,
        verbose="debug",
    )
    assert _is_mri_subject("flachkopf", tmp_path), "Scaling failed"
    spath = tmp_path / "flachkopf" / "bem"
    spath_fname = "flachkopf-%s-src.fif"

    assert (spath / (spath_fname % "ico-0")).exists()
    assert (tmp_path / "flachkopf" / "surf" / "lh.sphere.reg").is_file()
    vsrc_s = mne.read_source_spaces(spath / (spath_fname % "vol-50"))
    for vox in ([0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 2, 3]):
        idx = np.ravel_multi_index(vox, vsrc[0]["shape"], order="F")
        err_msg = f"idx={idx} @ {vox}, scale={scale}"
        assert_allclose(
            apply_trans(vsrc[0]["src_mri_t"], vox), vsrc[0]["rr"][idx], err_msg=err_msg
        )
        assert_allclose(
            apply_trans(vsrc_s[0]["src_mri_t"], vox),
            vsrc_s[0]["rr"][idx],
            err_msg=err_msg,
        )
    scale_labels("flachkopf", subjects_dir=tmp_path)

    # add distances to source space after hacking the properties to make
    # it run *much* faster
    src_dist = src.copy()
    for s in src_dist:
        s.update(rr=s["rr"][s["vertno"]], nn=s["nn"][s["vertno"]], tris=s["use_tris"])
        s.update(
            np=len(s["rr"]),
            ntri=len(s["tris"]),
            vertno=np.arange(len(s["rr"])),
            inuse=np.ones(len(s["rr"]), int),
        )
    mne.add_source_space_distances(src_dist)
    write_source_spaces(bem_path / (bem_fname % "ico-0"), src_dist, overwrite=True)

    # scale with distances
    os.remove(spath / (spath_fname % "ico-0"))
    scale_source_space("flachkopf", "ico-0", subjects_dir=tmp_path)
    ssrc = mne.read_source_spaces(spath / (spath_fname % "ico-0"))
    assert ssrc[0]["dist"] is not None
    assert ssrc[0]["nearest"] is not None

    # check patch info computation (only if SciPy is new enough to be fast)
    for s in src_dist:
        for key in ("dist", "dist_limit"):
            s[key] = None
    write_source_spaces(bem_path / (bem_fname % "ico-0"), src_dist, overwrite=True)

    # scale with distances
    os.remove(spath / (spath_fname % "ico-0"))
    scale_source_space("flachkopf", "ico-0", subjects_dir=tmp_path)
    ssrc = mne.read_source_spaces(spath / (spath_fname % "ico-0"))
    assert ssrc[0]["dist"] is None
    assert ssrc[0]["nearest"] is not None


@pytest.mark.slowtest  # can take forever on OSX Travis
@testing.requires_testing_data
def test_scale_mri_xfm(tmp_path, few_surfaces, subjects_dir_tmp_few):
    """Test scale_mri transforms and MRI scaling."""
    pytest.importorskip("nibabel")
    # scale fsaverage
    sample_dir = subjects_dir_tmp_few / "sample"
    subject_to = "flachkopf"
    spacing = "oct2"
    for subject_from in ("fsaverage", "sample"):
        if subject_from == "fsaverage":
            scale = 1.0  # single dim
        else:
            scale = [0.9, 2, 0.8]  # separate
        src_from_fname = (
            subjects_dir_tmp_few
            / subject_from
            / "bem"
            / (f"{subject_from}-{spacing}-src.fif")
        )
        src_from = mne.setup_source_space(
            subject_from,
            spacing,
            subjects_dir=subjects_dir_tmp_few,
            add_dist=False,
        )
        write_source_spaces(src_from_fname, src_from)
        vertices_from = np.concatenate([s["vertno"] for s in src_from])
        assert len(vertices_from) == 36
        hemis = [0] * len(src_from[0]["vertno"]) + [1] * len(src_from[0]["vertno"])
        mni_from = mne.vertex_to_mni(
            vertices_from, hemis, subject_from, subjects_dir=subjects_dir_tmp_few
        )
        if subject_from == "fsaverage":  # identity transform
            source_rr = np.concatenate([s["rr"][s["vertno"]] for s in src_from]) * 1e3
            assert_allclose(mni_from, source_rr)
        if subject_from == "fsaverage":
            overwrite = skip_fiducials = False
        else:
            with pytest.raises(OSError, match="No fiducials file"):
                scale_mri(
                    subject_from,
                    subject_to,
                    scale,
                    subjects_dir=subjects_dir_tmp_few,
                )
            skip_fiducials = True
            with pytest.raises(OSError, match="already exists"):
                scale_mri(
                    subject_from,
                    subject_to,
                    scale,
                    subjects_dir=subjects_dir_tmp_few,
                    skip_fiducials=skip_fiducials,
                )
            overwrite = True
        if subject_from == "sample":  # support for not needing all surf files
            os.remove(sample_dir / "surf" / "lh.curv")
        scale_mri(
            subject_from,
            subject_to,
            scale,
            subjects_dir=subjects_dir_tmp_few,
            verbose="debug",
            overwrite=overwrite,
            skip_fiducials=skip_fiducials,
        )
        if subject_from == "fsaverage":
            assert _is_mri_subject(subject_to, subjects_dir_tmp_few)
        src_to_fname = (
            subjects_dir_tmp_few
            / subject_to
            / "bem"
            / (f"{subject_to}-{spacing}-src.fif")
        )
        assert src_to_fname.exists(), "Source space was not scaled"
        # Check MRI scaling
        fname_mri = subjects_dir_tmp_few / subject_to / "mri" / "T1.mgz"
        assert fname_mri.exists(), "MRI was not scaled"
        # Check MNI transform
        src = mne.read_source_spaces(src_to_fname)
        vertices = np.concatenate([s["vertno"] for s in src])
        assert_array_equal(vertices, vertices_from)
        mni = mne.vertex_to_mni(
            vertices, hemis, subject_to, subjects_dir=subjects_dir_tmp_few
        )
        assert_allclose(mni, mni_from, atol=1e-3)  # 0.001 mm
        # Check head_to_mni (the `trans` here does not really matter)
        trans = rotation(0.001, 0.002, 0.003) @ translation(0.01, 0.02, 0.03)
        trans = Transform("head", "mri", trans)
        pos_head_from = np.random.RandomState(0).randn(4, 3)
        pos_mni_from = mne.head_to_mni(
            pos_head_from, subject_from, trans, subjects_dir_tmp_few
        )
        pos_mri_from = apply_trans(trans, pos_head_from)
        pos_mri = pos_mri_from * scale
        pos_head = apply_trans(invert_transform(trans), pos_mri)
        pos_mni = mne.head_to_mni(pos_head, subject_to, trans, subjects_dir_tmp_few)
        assert_allclose(pos_mni, pos_mni_from, atol=1e-3)
        # another way
        pos_mri_from_2 = mne.head_to_mri(
            pos_head_from, subject_from, trans, subjects_dir_tmp_few
        )
        pos_mri_from_ras = mne.head_to_mri(
            pos_head_from,
            subject_from,
            trans,
            subjects_dir_tmp_few,
            kind="ras",
        )
        mri_eq_ras = np.allclose(pos_mri_from_2, pos_mri_from_ras, atol=1e-1)
        if subject_from == "fsaverage":
            assert mri_eq_ras  # fsaverage is special this way
        else:
            assert not mri_eq_ras  # sample is not
        assert_allclose(pos_mri_from_2, 1e3 * pos_mri_from, atol=1e-3)
        with pytest.raises(OSError, match=r"parameters\.cfg"):
            mne.head_to_mri(
                pos_head_from,
                subject_from,
                trans,
                subjects_dir_tmp_few,
                unscale=True,
                kind="mri",
            )
        # yet another way
        pos_mri_from_3 = mne.head_to_mri(
            pos_head,
            subject_to,
            trans,
            subjects_dir_tmp_few,
            kind="mri",
            unscale=True,
        )
        assert_allclose(pos_mri_from_3, 1e3 * pos_mri_from, atol=1e-3)


def test_fit_matched_points():
    """Test fit_matched_points: fitting two matching sets of points."""
    tgt_pts = np.random.RandomState(42).uniform(size=(6, 3))

    # rotation only
    trans = rotation(2, 6, 3)
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, translate=False, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(tgt_pts, est_pts, 2, "fit_matched_points with rotation")

    # rotation & translation
    trans = np.dot(translation(2, -6, 3), rotation(2, 6, 3))
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(
        tgt_pts, est_pts, 2, "fit_matched_points with rotation and translation."
    )

    # rotation & translation & scaling
    trans = reduce(
        np.dot, (translation(2, -6, 3), rotation(1.5, 0.3, 1.4), scaling(0.5, 0.5, 0.5))
    )
    src_pts = apply_trans(trans, tgt_pts)
    trans_est = fit_matched_points(src_pts, tgt_pts, scale=1, out="trans")
    est_pts = apply_trans(trans_est, src_pts)
    assert_array_almost_equal(
        tgt_pts,
        est_pts,
        2,
        "fit_matched_points with rotation, translation and scaling.",
    )

    # test exceeding tolerance
    tgt_pts[0, :] += 20
    pytest.raises(RuntimeError, fit_matched_points, tgt_pts, src_pts, tol=10)


@testing.requires_testing_data
def test_get_mni_fiducials():
    """Test get_mni_fiducials."""
    pytest.importorskip("nibabel")
    fids, coord_frame = read_fiducials(fid_fname)
    assert coord_frame == FIFF.FIFFV_COORD_MRI
    assert [f["ident"] for f in fids] == list(range(1, 4))
    fids = np.array([f["r"] for f in fids])
    fids_est = get_mni_fiducials("sample", subjects_dir)
    fids_est = np.array([f["r"] for f in fids_est])
    dists = np.linalg.norm(fids - fids_est, axis=-1) * 1000.0  # -> mm
    assert (dists < 8).all(), dists


@pytest.mark.slowtest
@testing.requires_testing_data
@pytest.mark.parametrize(
    "scale_mode,ref_scale,grow_hair,fiducials,fid_match",
    [
        (None, [1.0, 1.0, 1.0], 0.0, None, "nearest"),
        (None, [1.0, 1.0, 1.0], 0.0, "estimated", "nearest"),
        (None, [1.0, 1.0, 1.0], 2.0, "auto", "nearest"),
        ("uniform", [1.0, 1.0, 1.0], 0.0, None, "nearest"),
        ("3-axis", [1.0, 1.0, 1.0], 0.0, "auto", "nearest"),
        ("uniform", [0.8, 0.8, 0.8], 0.0, "auto", "nearest"),
        ("3-axis", [0.8, 1.2, 1.2], 0.0, "auto", "matched"),
    ],
)
def test_coregistration(scale_mode, ref_scale, grow_hair, fiducials, fid_match):
    """Test automated coregistration."""
    pytest.importorskip("nibabel")
    subject = "sample"
    if fiducials is None:
        fiducials, coord_frame = read_fiducials(fid_fname)
        assert coord_frame == FIFF.FIFFV_COORD_MRI
    info = read_info(raw_fname)
    for d in info["dig"]:
        d["r"] = d["r"] * ref_scale
    trans = read_trans(trans_fname)
    coreg = Coregistration(
        info, subject=subject, subjects_dir=subjects_dir, fiducials=fiducials
    )
    assert np.allclose(coreg._last_parameters, coreg._parameters)
    assert len(coreg.fiducials.dig) == 3
    for dig_point in coreg.fiducials.dig:
        assert dig_point["coord_frame"] == FIFF.FIFFV_COORD_MRI
        assert dig_point["kind"] == FIFF.FIFFV_POINT_CARDINAL

    coreg.set_fid_match(fid_match)
    default_params = list(coreg._default_parameters)
    coreg.set_rotation(default_params[:3])
    coreg.set_translation(default_params[3:6])
    coreg.set_scale(default_params[6:9])
    coreg.set_grow_hair(grow_hair)
    coreg.set_scale_mode(scale_mode)
    # Identity transform
    errs_id = coreg.compute_dig_mri_distances()
    is_scaled = ref_scale != [1.0, 1.0, 1.0]
    id_max = 0.03 if is_scaled and scale_mode == "3-axis" else 0.02
    assert 0.005 < np.median(errs_id) < id_max
    # Fiducial transform + scale
    coreg.fit_fiducials(verbose=True)
    assert coreg._extra_points_filter is None
    coreg.omit_head_shape_points(distance=0.02)
    assert coreg._extra_points_filter is not None
    errs_fid = coreg.compute_dig_mri_distances()
    assert_array_less(0, errs_fid)
    if is_scaled or scale_mode is not None:
        fid_max = 0.05
        fid_med = 0.02
    else:
        fid_max = 0.03
        fid_med = 0.01
    assert_array_less(errs_fid, fid_max)
    assert 0.001 < np.median(errs_fid) < fid_med
    assert not np.allclose(coreg._parameters, default_params)
    coreg.omit_head_shape_points(distance=-1)
    coreg.omit_head_shape_points(distance=5.0 / 1000)
    assert coreg._extra_points_filter is not None
    # ICP transform + scale
    coreg.fit_icp(verbose=True)
    assert isinstance(coreg.trans, Transform)
    errs_icp = coreg.compute_dig_mri_distances()
    assert_array_less(0, errs_icp)
    if is_scaled or scale_mode == "3-axis":
        icp_max = 0.015
    else:
        icp_max = 0.01
    assert_array_less(errs_icp, icp_max)
    assert 0.001 < np.median(errs_icp) < 0.004
    assert (
        np.rad2deg(
            _angle_between_quats(
                rot_to_quat(coreg.trans["trans"][:3, :3]),
                rot_to_quat(trans["trans"][:3, :3]),
            )
        )
        < 13
    )
    if scale_mode is None:
        atol = 1e-7
    else:
        atol = 0.35
    assert_allclose(coreg._scale, ref_scale, atol=atol)
    coreg.reset()
    assert_allclose(coreg._parameters, default_params)


@pytest.mark.slowtest
@testing.requires_testing_data
def test_coreg_class_gui_match():
    """Test that using Coregistration matches mne coreg."""
    pytest.importorskip("nibabel")
    fiducials, _ = read_fiducials(fid_fname)
    info = read_info(raw_fname)
    coreg = Coregistration(
        info, subject="sample", subjects_dir=subjects_dir, fiducials=fiducials
    )
    assert_allclose(coreg.trans["trans"], np.eye(4), atol=1e-6)
    # mne coreg -s sample -d subjects -f MEG/sample/sample_audvis_trunc_raw.fif
    # then "Fit Fid.", Save... to get trans, read_trans:
    want_trans = [
        [9.99428809e-01, 2.94733196e-02, 1.65350307e-02, -8.76054692e-04],
        [-1.92420650e-02, 8.98512006e-01, -4.38526988e-01, 9.39774036e-04],
        [-2.77817696e-02, 4.37958330e-01, 8.98565888e-01, -8.29207990e-03],
        [0, 0, 0, 1],
    ]
    coreg.set_fid_match("matched")
    coreg.fit_fiducials(verbose=True)
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)
    # Set ICP iterations to one, click "Fit ICP"
    want_trans = [
        [9.99512792e-01, 2.80128177e-02, 1.37659665e-02, 6.08855276e-04],
        [-1.91694051e-02, 8.98992002e-01, -4.37545270e-01, 9.66848747e-04],
        [-2.46323701e-02, 4.37068194e-01, 8.99091005e-01, -1.44129358e-02],
        [0, 0, 0, 1],
    ]
    coreg.fit_icp(1, verbose=True)
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)
    # Set ICP iterations to 20, click "Fit ICP"
    with catch_logging() as log:
        coreg.fit_icp(20, verbose=True)
    log = log.getvalue()
    want_trans = [
        [9.97582495e-01, 2.12266613e-02, 6.61706254e-02, -5.07694029e-04],
        [1.81089472e-02, 8.39900672e-01, -5.42437911e-01, 7.81218382e-03],
        [-6.70908988e-02, 5.42324841e-01, 8.37485850e-01, -2.50057746e-02],
        [0, 0, 0, 1],
    ]
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)
    assert "ICP 19" in log
    assert "ICP 20" not in log  # converged on 19
    # Change to uniform scale mode, "Fit Fiducials" in scale UI
    coreg.set_scale_mode("uniform")
    coreg.fit_fiducials()
    want_scale = [0.975] * 3
    want_trans = [
        [9.99428809e-01, 2.94733196e-02, 1.65350307e-02, -9.25998494e-04],
        [-1.92420650e-02, 8.98512006e-01, -4.38526988e-01, -1.03350170e-03],
        [-2.77817696e-02, 4.37958330e-01, 8.98565888e-01, -9.03170835e-03],
        [0, 0, 0, 1],
    ]
    assert_allclose(coreg.scale, want_scale, atol=5e-4)
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)
    # Click "Fit ICP" in scale UI
    with catch_logging() as log:
        coreg.fit_icp(20, verbose=True)
    log = log.getvalue()
    assert "ICP 18" in log
    assert "ICP 19" not in log
    want_scale = [1.036] * 3
    want_trans = [
        [9.98992383e-01, 1.72388796e-02, 4.14364934e-02, 6.19427126e-04],
        [6.80460501e-03, 8.54430079e-01, -5.19521892e-01, 5.58008114e-03],
        [-4.43605632e-02, 5.19280374e-01, 8.53451848e-01, -2.03358755e-02],
        [0, 0, 0, 1],
    ]
    assert_allclose(coreg.scale, want_scale, atol=5e-4)
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)
    # Change scale mode to 3-axis, click "Fit ICP" in scale UI
    coreg.set_scale_mode("3-axis")
    with catch_logging() as log:
        coreg.fit_icp(20, verbose=True)
    log = log.getvalue()
    assert "ICP  7" in log
    assert "ICP  8" not in log
    want_scale = [1.025, 1.010, 1.121]
    want_trans = [
        [9.98387098e-01, 2.04762165e-02, 5.29526398e-02, 4.97257097e-05],
        [1.13287698e-02, 8.42087150e-01, -5.39222538e-01, 7.09863892e-03],
        [-5.56319728e-02, 5.38952649e-01, 8.40496957e-01, -1.46372067e-02],
        [0, 0, 0, 1],
    ]
    assert_allclose(coreg.scale, want_scale, atol=5e-4)
    assert_allclose(coreg.trans["trans"], want_trans, atol=1e-6)


@testing.requires_testing_data
@pytest.mark.parametrize(
    "drop_point_kind",
    (
        FIFF.FIFFV_POINT_CARDINAL,
        FIFF.FIFFV_POINT_HPI,
        FIFF.FIFFV_POINT_EXTRA,
        FIFF.FIFFV_POINT_EEG,
    ),
)
def test_coreg_class_init(drop_point_kind):
    """Test that Coregistration can be instantiated with various digs."""
    pytest.importorskip("nibabel")
    fiducials, _ = read_fiducials(fid_fname)
    info = read_info(raw_fname)

    dig_list = []
    eeg_chans = []
    for pt in info["dig"]:
        if pt["kind"] != drop_point_kind:
            dig_list.append(pt)
            if pt["kind"] == FIFF.FIFFV_POINT_EEG:
                eeg_chans.append(f"EEG {pt['ident']:03d}")

    this_info = info.copy()
    this_info.set_montage(
        DigMontage(dig=dig_list, ch_names=eeg_chans), on_missing="ignore"
    )
    Coregistration(
        this_info, subject="sample", subjects_dir=subjects_dir, fiducials=fiducials
    )