Bad results after removing digital points under Nasion

Support & Discussions
, , ,
1

I am trying to reproduce this paper: https://doi.org/10.1016/j.neuroimage.2021.117822. In it, they used the old MNE version which may be 0.18 to do the MEG coregistration. They removed the digital points below Nasion because the subjects in the dataset have been removed the face for privacy protection. They calculated the distance between MRI anatomical points and MEG digital points to ensure it is small enough. Below is the original code:
def coreg_head2mri(subjects_dir, subject, native_fid, raw_path, raw_NosePtsOut, trans_dst, flag_fid=False):
import scipy.io
import numpy as np
from mne.gui._coreg_gui import CoregModel

model = CoregModel()
model.mri.subjects_dir = subjects_dir
model.mri.subject = subject

# Remove Polhemus points around the nose (y>0, z<0)
model.hsp.file = raw_path
head_pts = model.hsp.points
raw  = read_raw_fif(raw_path, preload=True)  # !
pos  = np.where((head_pts[:,1] <= 0) | (head_pts[:,2] >= 0))                                   
dig  = raw.info['dig']
dig2 = dig[0:8]
dig3 = [dig[p+7] for p in pos[0]]
dig_yeah = dig2+dig3
raw.info['dig'] = dig_yeah
raw.save(raw_NosePtsOut, overwrite=True) # !

model.hsp.file = raw_NosePtsOut

# Load CamCAN fiducials from matlab file
if flag_fid:
    fid = scipy.io.loadmat(native_fid, struct_as_record=False, squeeze_me=True) # !
    fid = fid['fid']

    model.mri.lpa = np.reshape(fid.native.mm.lpa*0.001,(1,3))
    model.mri.nasion = np.reshape(fid.native.mm.nas*0.001,(1,3))
    model.mri.rpa = np.reshape(fid.native.mm.rpa*0.001,(1,3))
    assert (model.mri.fid_ok)
    
    lpa_distance = model.lpa_distance
    nasion_distance = model.nasion_distance
    rpa_distance = model.rpa_distance
          
    model.nasion_weight = 1.
    model.fit_fiducials(0)
    old_x = lpa_distance ** 2 + rpa_distance ** 2 + nasion_distance ** 2
    new_x = (model.lpa_distance ** 2 + model.rpa_distance ** 2 +
             model.nasion_distance ** 2)
    assert new_x < old_x
    
avg_point_distance = np.mean(model.point_distance)

while True:
    model.fit_icp(0)
    new_dist = np.mean(model.point_distance)
    assert new_dist < avg_point_distance
    if model.status_text.endswith('converged)'):
        break

model.save_trans(trans_dst) # !
trans = mne.read_trans(trans_dst)
np.testing.assert_allclose(trans['trans'], model.head_mri_t)

But now, I am using MNE version 1.7 to reproduce it. So I try to use the automated co-registration function to do the coregistration. Below is my code:
def coreg_head2mri(subjects_dir, subject, raw_path, raw_NosePtsOut, trans_dst):
import numpy as np
from mne.coreg import Coregistration

raw  = read_raw_fif(raw_path, preload=True)

dig = raw.info['dig']
dig_points_list = []
for i,dig_obj in enumerate(raw.info['dig']):
    # print(i,dig_obj['r'])
    dig_points_list.append(dig_obj['r'])
dig_points_array = np.array(dig_points_list)

pos  = np.where((dig_points_array[:,1] <= 0) | (dig_points_array[:,2] >= 0))
pos_np = pos[0]


dig2 = dig_points_array[0:8]

# 将原始dig当中pos中7以后的数据保存在dig3当中
dig3_list = []
for p in pos_np[8:]:
    # print(p)
    dig3_list.append(dig[p])
    # print('\n',dig3_list)


dig3_points_list = []
for i,dig3_obj in enumerate(dig3_list):
    print(i,dig3_obj['r'])
    dig3_points_list.append(dig3_obj['r'])
dig3 = np.array(dig3_points_list)


dig_yeah =np.concatenate((dig2, dig3), axis=0, dtype='float32')


new_montage = make_dig_montage(ch_pos=None, nasion=dig_yeah[1], lpa=dig_yeah[0], rpa=dig_yeah[2], hpi=dig_yeah[3:])

raw.set_montage(montage=new_montage)

raw.save(raw_NosePtsOut, overwrite=True)


raw_NosePtsOut_info = read_info(raw_NosePtsOut_dir)
fiducials = raw_NosePtsOut_info['dig']


model = Coregistration(raw_NosePtsOut, subject, subjects_dir, fiducials=fiducials)

model.fit_fiducials(verbose=True)

model.fit_icp(n_iterations=6, nasion_weight=2.0, verbose=True)
model.fit_icp(n_iterations=30,lpa_weight=2.0, nasion_weight=10.0,rpa_weight=2.0, verbose=True)

dists = model.compute_dig_mri_distances() * 1e3  # in mm
print(
f" {subject} Distance between HSP and MRI (mean/min/max):\n{np.mean(dists):.2f} mm "
f"/ {np.min(dists):.2f} mm / {np.max(dists):.2f} mm"
)

mne.write_trans(trans_dst, model.trans)

But I got very bad results from my code:

image
image1920×999 45.9 KB

Now, I’m trying to use Brainstorm to do the co-registration and export the result after finishing the co-registration. Is it possible? Or should I continue to work on the old MNE version?

2

Hi,

I haven’t run your code, but I noticed you’re using ICP for automated co-registration. In my experience, including the nose (part below the nasion) in both the MRI and MEG improves the ICP alignment. The nose from the MEG “locks” to the nose from the MRI and improves the ICP. Without the nose, the headshape generated from the MEG digital points is pretty round and the ICP might rotate it in relation to the MRI.

Best,
Konstantinos

1 Like