agramfort
commented
3 years ago Did you track down what files are used to see where it can go wrong? it's going to be hard to debug without access to the data...
Open yh-luo opened 3 years ago
agramfort
commented
3 years ago Did you track down what files are used to see where it can go wrong? it's going to be hard to debug without access to the data...
yh-luo
commented
3 years ago @agramfort I think all the files are fine. They look correct in FreeView. That's why I suspect the problem is from the rendering process. I updated the information for clarity (hopefully).
If I remember correctly, bem/subject-head.fif (head in coregistration gui) is converted from bem/subject_outer_skin_surface, and bem/subject-head-dense.fif (high-resolution head) is converted from subject/mri/lh.seghead
Both subject/mri/lh.seghead and subject_outer_skin_surface looks correct in FreeView (the first image, subject_outer_skin_surface not showing).
From the first line in the error messages, I think the problem is gui related. subject-346 was the specified subject but subject-102's file (another subject) was mentioned, which is weird.
Exception occurred in traits notification handler for object: <mne.gui._file_traits.SurfaceSource object at 0x7f3db86bcb30>, trait: file, old value: subjects//subject-102/bem/subject-102-head-dense.fif, new value: subjects//subject-346/bem/subject-346-head-dense.fifI would look into the fif files you have to see if the information on subject they contain is correct
yh-luo
commented
3 years ago I'm not that familiar with the fif files sadly... subjects/subject-346/bem can be downloaded via the link.
agramfort
commented
3 years ago you can use:
In [25]: import mne
In [26]: pos, tri = mne.read_surface('outer_skin.surf')
In [27]: from mayavi import mlab
In [28]: mlab.triangular_mesh(*pos.T, tri) Out[28]: <mayavi.modules.surface.Surface at 0x7fe64be30720>
to plot the .surf files. For the fif bem files you can use:
In [37]: surf = mne.bem.read_bem_surfaces('subject-346-head.fif')[0]
In [38]: mlab.triangular_mesh(*surf['rr'].T, surf['tris'])
I let you dig in
larsoner
commented
3 years ago I think all the files are fine. They look correct in FreeView
Based on this output:
RuntimeError: Surface outer skin has topological defects: 13 / 383510 vertices have fewer than three neighboring triangles [44385, 44387, 44389, 45198, 46031, 46033, 46861, 46862, 47708, 48546, 49428, 50338, 53081]I don't think the surfaces are currently usable. They might look fine in FreeView, but mne coreg (and other parts of MNE code) rely on the surface being topologically correct, and it's not. A topologically incorrect surface (with holes, missing neighbors, etc.) can look just fine by eye.
The coreg error is that when it tries to switch surfaces between subjects it fails, leaving things in an inconsistent state. That's a bug, but if we fix your surfaces it'll get you going again and avoid it. And even if we make the failure more graceful then you'll still need to fix the surface for that subject in order to be able to process them, so let's work on that first.
As a first step, you can try using the code here:
It has some dependencies like pymeshfix and trimesh that you'll have to install, but basically I would 1) backup the original files, 2) read them with read_surface or read_bem_surfaces, 3) fix them using this function, and then 4) write them back out with write_surface or some code hacked together from our private bem functions (note: using private functions is not a stable long-term solution as the API can change at any time).
In the meantime I could work on a PR to make these mesh fixing functions public (albeit requiring these packages) and then add a write_head_bem or so to allow people to write out custom head.fif files.
yh-luo
commented
3 years ago @larsoner thanks for the explanation. Yes, there are holes in the head eyeing from FreeView.
I understand the importance of BEM surfaces and follow the amazing Blender tutorial to fix a lot of surfaces, but only for brain, inner_skull and outer_skull which I think will be used in MEG source reconstruction (inner_skull to be specific).
I wonder if head.fif and head-dense.fif need to be fixed too. That would be a painful job because nearly all my subjects have topological defects in the head-dense.fif due to wearing goggles and headphones in the scan. The data were collected before I'm in the loop so 🤦
Another thing intrigued me is that mne coreg only complains this subject, but not all the subjects.
If I understand correctly, these head files are only used for coregistration? My current workaround is to use FreeView to view the T1 images and decide the possible MNI coordinates for fiducial points and use that in the gui. Then, the usual alignment process is followed. Do you think that's good enough for coregistration?
The code and the paper look great (infant template brains are life-saving!), I will try that!
larsoner
commented
3 years ago I wonder if head.fif and head-dense.fif need to be fixed too.
Some of the computations even in mne coreg assume surfaces are topologically correct, such as the inside/outside coloring of points. Really it's best to fix the issues, and if we can make a function / use the functions from the infant MRI templates, then you can loop over all of your subjects and fix them all in one go without needing any manual intervention.
If I understand correctly, these head files are only used for coregistration? My current workaround is to use FreeView to view the T1 images and decide the possible MNI coordinates for fiducial points and use that in the gui. Then, the usual alignment process is followed. Do you think that's good enough for coregistration?
I would ensure that the actual correct subject surface is being used, otherwise it seems less reliable.
yh-luo
commented
3 years ago It is easier than I previously thought thanks to the well-organized code.
I was able to fix(?) the high-resolution head using fix_all_defects. mne coreg accepts the head-dense.fif though the holes are still obvious. Not quite sure why...
Since it's more like a real-world problem than a bug, I'm closing the issue.
Before Unable to load.
After (head-dense.fif in mne coreg)
Codes
import os.path as op
import shutil
import mne
from mne.bem import _read_bem_surfaces_fif, _surfaces_to_bem, read_surface
from mne.io.constants import FIFF
from config import subjects_dir
from template_building import fix_all_defects
subject = 'subject-346'
bem_dir = op.join(subjects_dir, subject, 'bem')
# Load the surface
surf = _read_bem_surfaces_fif(op.join(bem_dir, f'{subject}-head-dense.fif'),
None)[0]
coords = surf['rr']
faces = surf['tris']
fix_coords, fixed_faces = fix_all_defects(coords, faces)
# get info
# _, _, volume_info = read_surface(op.join(bem_dir, 'watershed',
# f'{subject}_outer_skin_surface'),
# read_metadata=True)
# surf to bem for .fif
surf = _surfaces_to_bem([dict(rr=fix_coords, tris=fixed_faces)],
[FIFF.FIFFV_BEM_SURF_ID_HEAD], [1],
rescale=False)
fname_dense_head = op.join(bem_dir, f'{subject}-head-dense.fif')
# copy the original file
if op.isfile(fname_dense_head):
shutil.copy(fname_dense_head,
fname_dense_head.replace('dense', 'dense_orig'))
# Overwrite the surface
mne.write_bem_surfaces(fname_dense_head, surf, overwrite=True)Essential template_building.py (functions that are not used or relying on unpublished codes are removed)
import numpy as np
import pymeshfix
import trimesh
from mne.bem import _get_solids
from mne.surface import _triangle_neighbors
def fix_all_defects(vertices, faces_):
if has_degenerated_faces(vertices, faces_):
vertices, faces_ = remove_degenerated_faces(vertices, faces_)
assert (not has_degenerated_faces(vertices, faces_))
if has_topological_defects(vertices, faces_):
print("The decimated mesh has topological defects. Fixing it.")
vertices, faces_ = fix_topological_defects(vertices, faces_)
if has_degenerated_faces(vertices, faces_):
vertices, faces_ = remove_degenerated_faces(vertices, faces_)
assert (not has_topological_defects(vertices, faces_))
if not surface_is_complete(vertices, faces_) or not trimesh.Trimesh(
vertices, faces_).is_watertight:
print("The decimated mesh has holes. Fixing it.")
vertices, faces_ = repair_holes(vertices, faces_)
check_mesh(vertices, faces_)
return vertices, faces_
def has_degenerated_faces(vertices, faces_):
return not np.all(
trimesh.Trimesh(vertices, faces_).remove_degenerate_faces())
def remove_degenerated_faces(vertices, faces_):
mesh_ = trimesh.Trimesh(vertices, faces_)
mesh_.remove_degenerate_faces()
return mesh_.vertices, mesh_.faces
def has_topological_defects(vertices, faces_):
zero, one, two = get_topological_defects(vertices, faces_)
return len(zero) or len(one) or len(two)
# Code extracted and slighly modified from mne.surface.complete_surface_info
# for compactness of the example
def fix_topological_defects(vertices, faces_):
zero, one, two = get_topological_defects(vertices, faces_)
ind_faces_to_remove = []
if len(zero) > 0:
print(' Vertices do not have any neighboring '
'triangles: [%s]' % ', '.join(str(z) for z in zero))
print(' Correcting by removing these vertices.')
if len(one) > 0:
print(' Vertices have only one neighboring '
'triangles: [%s]' % ', '.join(str(tri) for tri in one))
print(
' Correcting by removing these vertices and their neighboring triangles.'
)
ind_faces_to_remove.extend(np.where(faces_ == one)[0].tolist())
if len(two) > 0:
print(' Vertices have only two neighboring '
'triangles, removing neighbors: [%s]' %
', '.join(str(tri) for tri in two))
print(' Correcting by merging the two neighboring '
'triangles and removing the faulty vertices.')
ind_faces, faces_to_add = correction_two_neighboring_tri(
vertices, faces_, two)
ind_faces_to_remove.extend(ind_faces)
faces_ = np.concatenate((np.delete(faces_,
np.array(ind_faces_to_remove,
dtype=int),
axis=0), faces_to_add))
vertices_to_remove = np.concatenate((zero, one, two)).astype(int)
if len(vertices_to_remove):
vertices, faces_ = reindex_vertices(vertices, faces_,
vertices_to_remove)
else:
print("No issue found with the mesh.")
return vertices, faces_
def surface_is_complete(vertices, faces_):
"""Check the sum of solid angles as seen from inside."""
cm = vertices.mean(axis=0)
tot_angle = _get_solids(vertices[faces_], cm[np.newaxis, :])[0]
prop = tot_angle / (2 * np.pi)
return np.abs(prop - 1.0) < 1e-5
def repair_holes(vertices, faces_):
# trimesh has a hole fixing function, but it just deals with
# 3 or 4 vertices holes.
meshfix = pymeshfix.MeshFix(vertices, faces_)
meshfix.repair()
vertices = meshfix.v # numpy np.float array
faces_ = meshfix.f # numpy np.int32 array
# The return mesh has a solid angle of -1 instead of 1.
# Correcting this.
mesh_ = trimesh.Trimesh(vertices=vertices, faces=faces_)
trimesh.repair.fix_normals(mesh_, multibody=False)
return mesh_.vertices, mesh_.faces
def check_mesh(vertices, faces_):
assert (surface_is_complete(vertices, faces_))
assert (not has_topological_defects(vertices, faces_))
assert (not has_degenerated_faces(vertices, faces_))
assert trimesh.Trimesh(vertices, faces_).is_watertight
def get_topological_defects(vertices, faces_):
# Find neighboring triangles, accumulate vertex normals, normalize
neighbor_tri = _triangle_neighbors(faces_, len(vertices))
# Check for topological defects
zero, one, two = list(), list(), list()
for ni, n in enumerate(neighbor_tri):
if len(n) < 3:
if len(n) == 0:
zero.append(ni)
elif len(n) == 1:
one.append(ni)
else:
two.append(ni)
return zero, one, two
def correction_two_neighboring_tri(vertices, faces_, faulty_vert_ind):
ind_faces_to_remove = []
new_faces = []
for ind in faulty_vert_ind:
ind_faces = np.where(faces_ == ind)[0]
ind_faces_to_remove.extend(ind_faces)
face1, face2 = faces_[ind_faces]
new_face = np.unique(np.concatenate((face1, face2)))
new_face = np.delete(new_face, np.where(new_face == ind))
assert (
len(new_face) == 3
) # If == 4, it means that face1 and face2 do not share a common edge
new_det = np.linalg.det(vertices[new_face])
assert new_det # If zero, the three points are colinear
# Align the normals
det1 = np.linalg.det(vertices[face1])
if np.sign(det1) == np.sign(new_det):
new_face = new_face[[1, 0, 2]]
new_faces.append(new_face)
return np.array(ind_faces_to_remove, dtype=int), new_faces
def reindex_vertices(vertices, faces_, ind_vertices_to_remove):
decrement = np.cumsum(
np.zeros(vertices.shape[0], dtype=int) +
np.in1d(np.arange(vertices.shape[0]), ind_vertices_to_remove))
vertices = np.delete(vertices, ind_vertices_to_remove, axis=0)
faces_ = faces_ - decrement[faces_]
return vertices, faces_@yh-luo I'll reopen until we get a variant of this code in MNE-Python. I've also needed to use it so that's at least three people, and these topology questions come up a lot so I think it's worth pursuing proper inclusion as a function.
@christian-oreilly since you did all the heavy lifting on this, would you be up for opening a EDIT: rough PR to add a mne.surface.fix_topology(rr, tris) function? Then I'm happy to hack away at the code to add tests, etc. to get it working, but the first commit will already get you code credit for the changes when we squash+merge, which I think is helpful since you did the heavy lifting!
I think in the end we probably also want a write_head_bem function since the code to do this currently requires poking around in private attributes. Then the fix could have maybe just been for you @yh-luo :
surf = mne.read_bem_surfaces(fname_head)[0]
surf['rr'], surf['tris'] = mne.surface.fix_topology(surf['rr'], surf['tris'])
mne.write_head_bem(fname_head, surf)@larsoner Another small thing is that mne.read_bem_surfaces complains about the topological defects and raise RuntimeError, that's why I used the private _read_bem_surfaces_fif to get around this.
To make your proposal work, an argument like on_missing in mne.Epochs (maybe on_defects) to let the user proceed to fix_topology would be nice? If so:
surf = mne.read_bem_surfaces(fname_head, on_defects='warn')[0]
surf['rr'], surf['tris'] = mne.surface.fix_topology(surf['rr'], surf['tris'])
mne.write_head_bem(fname_head, surf)RuntimeWarning: Surface outer skin has topological defects: 13 / 383510 vertices have fewer than three neighboring triangles [44385, 44387, 44389, 45198, 46031, 46033, 46861, 46862, 47708, 48546, 49428, 50338, 53081]Both write_head_bem and this are easy to implement, I can make a PR for it if no one is already up to.
larsoner
commented
3 years ago That would be great @yh-luo, please do!
christian-oreilly
commented
3 years ago Yes, that would be @yh-luo ; I could have a look at this issue at some point but not in a near future because I need to do a few other things in priority at the moment. So if you can tackle it, that is perfect!
tillhabersetzer
commented
2 years ago Thank you very much for providing the code. I had the same problem and it worked for me too! I would really appreciate to see this function mne.surface.fix_topology(rr, tris) implemented in mne python
Describe the bug
Someone else's high-resolution head was used for the selected subject.
lh.segheadin FreeView looks like:But the high-resolution head rendered in
mne coreggui looks like: (Notice the differences of the top-left corner of the goggle, the ears, and the nose bridge.) I can confirm that it is from other subject'shead-dense.fif(lh.seghead) via FreeView.The high-resolution head also did not align to the head. No rotations were made, the head looks like:
FreeView showed that
lh.segheadandouter_skin.surfwere aligned. From the above images, it looks like someone else's dense head was used. The specified subject is subject-346, the error messages indicated that subject-102 dense head was used instead?subject-102 is the first subject in the subject dropdown dialog. It seems to me that
mne coreggui renders the first subject and then update with the selected subject. Maybe the bug result from that process?Steps to reproduce
I can't provide a MWE because it only happened to this subject.
Expected results
subject-346's
lh.segheadActual results
subject-102(?)'s
lh.segheadAdditional information
The above images in this post will be removed in the future.
The high-resolution head was created with
mne make_scalp_surfaces.Logs in mne make_scalp_surfaces
``` 1. Creating a dense scalp tessellation with mkheadsurf... Running subprocess: mkheadsurf -subjid subject-346 -srcvol T1.mgz INFO: log file is /mnt/sdb1/MEG_Prosody/subjects/subject-346/scripts/mkheadsurf.log -------------------------------- Tue 02 Feb 2021 10:22:20 AM CST /mnt/sdb1/MEG_Prosody mri_seghead --invol /mnt/sdb1/MEG_Prosody/subjects/subject-346/mri/T1.mgz --outvol /mnt/sdb1/MEG_Prosody/subjects/subject-346/mri/seghead.mgz --fill 1 --thresh1 20 --thresh2 20 --nhitsmin 2 -------------------------------- input volume: /mnt/sdb1/MEG_Prosody/subjects/subject-346/mri/T1.mgz output volume: /mnt/sdb1/MEG_Prosody/subjects/subject-346/mri/seghead.mgz threshold1: 20 threshold2: 20 nhitsmin: 2 fill value: 1 Loading input volume Filling Columns Filling Rows Filling Slices Merging and Inverting Growing Counting N Head Voxels = 2763508 N Back Voxels = 14013708 Avg. Back Intensity = 1.099117 Max. Back Intensity = 235.000000 Writing output Done -------------------------------- Tue 02 Feb 2021 10:22:22 AM CST /mnt/sdb1/MEG_Prosody mri_tessellate /mnt/sdb1/MEG_Prosody/subjects/subject-346/mri/seghead.mgz 1 /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead -------------------------------- 7.1.1 7.1.1 slice 40: 4930 vertices, 5067 faces slice 50: 13321 vertices, 13533 faces slice 60: 23339 vertices, 23642 faces slice 70: 32759 vertices, 33068 faces slice 80: 41073 vertices, 41395 faces slice 90: 49432 vertices, 49789 faces slice 100: 58965 vertices, 59458 faces slice 110: 68826 vertices, 69328 faces slice 120: 77977 vertices, 78498 faces slice 130: 86895 vertices, 87425 faces slice 140: 95446 vertices, 96008 faces slice 150: 104054 vertices, 104633 faces slice 160: 112518 vertices, 113104 faces slice 170: 122539 vertices, 123275 faces slice 180: 135046 vertices, 136220 faces slice 190: 149553 vertices, 150907 faces slice 200: 164959 vertices, 166490 faces slice 210: 181548 vertices, 183291 faces slice 220: 189480 vertices, 191117 faces slice 230: 190186 vertices, 191755 faces slice 240: 190186 vertices, 191755 faces slice 250: 190186 vertices, 191755 faces using the conformed surface RAS to save vertex points... writing /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead using vox2ras matrix: -1.00000 0.00000 0.00000 128.00000; 0.00000 0.00000 1.00000 -128.00000; 0.00000 -1.00000 0.00000 128.00000; 0.00000 0.00000 0.00000 1.00000; -------------------------------- Tue 02 Feb 2021 10:22:24 AM CST /mnt/sdb1/MEG_Prosody mris_smooth -n 10 -b area.seghead -c curv.seghead /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead -------------------------------- smoothing for 10 iterations smoothing surface tessellation for 10 iterations... smoothing complete - recomputing first and second fundamental forms... writing smoothed curvature to /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.curv.seghead writing smoothed area to /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.area.seghead -------------------------------- Tue 02 Feb 2021 10:22:27 AM CST /mnt/sdb1/MEG_Prosody mris_inflate -n 10 -sulc sulc.seghead /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead.inflated -------------------------------- niterations = 10 sulc name = sulc.seghead Reading /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead avg radius = 87.9 mm, total surface area = 117737 mm^2 step 000: RMS=0.148 (target=0.015) step 005: RMS=0.114 (target=0.015) step 010: RMS=0.103 (target=0.015) step 015: RMS=0.097 (target=0.015) step 020: RMS=0.092 (target=0.015) step 025: RMS=0.088 (target=0.015) step 030: RMS=0.085 (target=0.015) step 035: RMS=0.082 (target=0.015) step 040: RMS=0.081 (target=0.015) step 045: RMS=0.080 (target=0.015) step 050: RMS=0.079 (target=0.015) step 055: RMS=0.080 (target=0.015) step 060: RMS=0.080 (target=0.015) writing inflated surface to /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.seghead.inflated writing sulcal depths to /mnt/sdb1/MEG_Prosody/subjects/subject-346/surf/lh.sulc.seghead inflation complete. inflation took 0.3 minutes mris_inflate utimesec 16.246170 mris_inflate stimesec 0.531678 mris_inflate ru_maxrss 299472 mris_inflate ru_ixrss 0 mris_inflate ru_idrss 0 mris_inflate ru_isrss 0 mris_inflate ru_minflt 743476 mris_inflate ru_majflt 0 mris_inflate ru_nswap 0 mris_inflate ru_inblock 0 mris_inflate ru_oublock 10448 mris_inflate ru_msgsnd 0 mris_inflate ru_msgrcv 0 mris_inflate ru_nsignals 0 mris_inflate ru_nvcsw 0 mris_inflate ru_nivcsw 1765 Started at: Tue 02 Feb 2021 10:22:20 AM CST Ended at: Tue 02 Feb 2021 10:22:43 AM CST mkheadsurf done 2. Creating /mnt/sdb1/MEG_Prosody/subjects/subject-346/bem/subject-346-head-dense.fif ... outer skin CM is 10.24 5.95 -47.79 mm Surfaces passed the basic topology checks. 3. Creating medium tessellation... 3.1 Decimating the dense tessellation... 3.2 Creating /mnt/sdb1/MEG_Prosody/subjects/subject-346/bem/subject-346-head-medium.fif outer skin CM is 18.27 27.64 -63.59 mm Surfaces passed the basic topology checks. 4. Creating sparse tessellation... 4.1 Decimating the dense tessellation... 4.2 Creating /mnt/sdb1/MEG_Prosody/subjects/subject-346/bem/subject-346-head-sparse.fif outer skin CM is 17.12 24.57 -62.60 mm Surfaces passed the basic topology checks. /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin has topological defects: 13 / 383510 vertices have fewer than three neighboring triangles [44385, 44387, 44389, 45198, 46031, 46033, 46861, 46862, 47708, 48546, 49428, 50338, 53081] Consider using --force as an additional input parameter. sys.exit(main()) /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin is not complete (sum of solid angles yielded 0.999793, should be 1.) sys.exit(main()) /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin has topological defects: 53 / 30000 vertices have fewer than three neighboring triangles [1659, 1663, 1690, 1732, 1734, 1767, 1812, 1813, 1815, 1849, 1850, 1887, 2012, 2127, 2212, 2272, 2507, 3476, 5047, 5435, 6169, 6296, 6297, 6532, 6540, 6587, 6976, 7295, 7855, 8106, 8199, 8528, 8533, 8538, 8727, 8795, 8953, 9039, 10498, 10546, 11537, 11560, 11755, 12293, 12587, 13594, 13779, 14204, 14207, 14300, 14313, 14463, 14532] Consider using --force as an additional input parameter. sys.exit(main()) /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin is not complete (sum of solid angles yielded 0.999567, should be 1.) sys.exit(main()) /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin has topological defects: 17 / 2500 vertices have fewer than three neighboring triangles [176, 180, 184, 189, 271, 272, 677, 775, 991, 1060, 1061, 1107, 1128, 1160, 1161, 1162, 1190] Consider using --force as an additional input parameter. sys.exit(main()) /home/ntubabylab/anaconda3/envs/mne/bin/mne:8: RuntimeWarning: Surface outer skin is not complete (sum of solid angles yielded 1.00046, should be 1.) sys.exit(main()) ```