mri.py

"""
export MRI=/mnt/bakraid/larsoner/mri/Infants/Sources/BEM/AVG7-5Months3T_segmented_BEM4.nii.gz
mri_binarize --i $MRI --o inner_skull.mgz --min 1 --max 2
mri_binarize --i $MRI --o outer_skull.mgz --min 1 --max 3
mri_binarize --i $MRI --o outer_skin.mgz --min 1 --max 4
mri_tessellate inner_skull.mgz 1 lh.inner_skull_dense
mri_tessellate outer_skull.mgz 1 lh.outer_skull_dense
mri_tessellate outer_skin.mgz 1 lh.outer_skin_dense
mris_extract_main_component lh.inner_skull_dense lh.inner_skull_dense
mris_extract_main_component lh.outer_skull_dense lh.outer_skull_dense
mris_extract_main_component lh.outer_skin_dense lh.outer_skin_dense
mris_smooth -nw lh.inner_skull_dense lh.inner_skull_dense.smooth
mris_smooth -nw lh.outer_skull_dense lh.outer_skull_dense.smooth
mris_smooth -nw lh.outer_skin_dense lh.outer_skin_dense.smooth
mris_inflate -no-save-sulc lh.inner_skull_dense.smooth lh.inner_skull_dense.inflated
mris_inflate -no-save-sulc lh.outer_skull_dense.smooth lh.outer_skull_dense.inflated
mris_inflate -no-save-sulc lh.outer_skin_dense.smooth lh.outer_skin_dense.inflated
mris_sphere -q lh.inner_skull_dense.inflated lh.inner_skull_dense.qsphere
mris_sphere -q lh.outer_skull_dense.inflated lh.outer_skull_dense.qsphere
mris_sphere -q lh.outer_skin_dense.inflated lh.outer_skin_dense.qsphere
"""  # noqa: E501
import os
import os.path as op
import numpy as np
import mne
ico_surf = mne.surface._get_ico_surface(4)
norms = np.linalg.norm(ico_surf['rr'], axis=1)
np.testing.assert_allclose(norms, 1., atol=1e-4)
os.makedirs('mysubj/bem', exist_ok=True)
for surf in ('inner_skull', 'outer_skull', 'outer_skin'):
    fname_dense = 'lh.' + surf + '_dense.smooth'
    rr, tris, volume_info = mne.read_surface(fname_dense, read_metadata=True)
    fname_sphere = op.splitext(fname_dense)[0] + '.qsphere'
    sphere_rr, _ = mne.read_surface(fname_sphere)
    norms = np.linalg.norm(sphere_rr, axis=1, keepdims=True)
    np.testing.assert_allclose(norms, 100.)  # freesurfer
    sphere_rr /= norms
    idx = mne.surface._compute_nearest(sphere_rr, ico_surf['rr'])
    rr, tris = rr[idx], ico_surf['tris']
    mne.write_surface('mysubj/bem/%s.surf' % surf, rr, tris,
                      volume_info=volume_info, overwrite=True)
bem = mne.make_bem_model('mysubj', subjects_dir='.', verbose=True)