Quick little script to get median distance between sensor helmet and head surface points (excluding the face)

mne.sample.helmet-distance.py

import mne
from mne.io.constants import FIFF
import os.path as op
import numpy as np
from scipy.stats import rankdata

#data_path = mne.datasets.sample.data_path(update_path=False)
data_path = '/export/research/analysis/human/jhouck/shared/tools/MNE-sample-data/MNE-sample-data'
subjects_dir = op.join(data_path, 'subjects')
subject='sample'
fname_raw = op.join(data_path, 'MEG/sample/sample_audvis_raw.fif')
fname_trans = op.join(data_path, 'MEG/sample/sample_audvis_raw-trans.fif')

info =  mne.io.read_info(fname_raw)
trans = mne.read_trans(fname_trans)  
helmet = mne.get_meg_helmet_surf(info, trans, verbose=False)
assert helmet['coord_frame'] == FIFF.FIFFV_COORD_MRI
assert info['dig'][7]['coord_frame'] == FIFF.FIFFV_COORD_HEAD
hpi_loc = np.array([d['r'] for d in info['dig']
                if d['kind'] == FIFF.FIFFV_POINT_HPI])
ext_loc = np.array([d['r'] for d in info['dig']
               if d['kind'] == FIFF.FIFFV_POINT_EXTRA])
tmppts = np.concatenate((ext_loc, hpi_loc), axis=0)
# get z & y thresholds from hpi_loc
zrank = rankdata(hpi_loc[:,2])
zthresh = hpi_loc[np.where(zrank==3), 2][0,0] #  second-highest z
ythresh = hpi_loc[np.where(zrank==3), 1][0,0] #  y coord of 2nd-highest z
# identify and remove face points (anything inferior and anterior to the lowest forehead HPI coil)
pts_bad = tmppts[(tmppts[:,2] < zthresh ) & (tmppts[:,1] > ythresh )]
mask = np.isin(tmppts, pts_bad, invert=True).astype(int)
pts_orig = tmppts * mask
pts_new = mne.transforms.apply_trans(trans, pts_orig[(pts_orig[:,0] != 0.)])
nearest, distances = mne.surface._compute_nearest(pts_new, helmet['rr'], method='BallTree', return_dists=True)
dist_helmet = np.median(distances) * 1e3
print('Median distance between helmet surface and head points is %.2f mm' % dist_helmet)
print('Minimum distance between helmet surface and head points is %.2f mm' % (np.min(distances) * 1e3))
mne.viz.plot_alignment(info=info, trans=trans, subject=subject, subjects_dir=subjects_dir, surfaces='seghead', meg='helmet')