Computing source space does not complete

eort · October 18, 2023, 2:57pm

Hi,
I am trying to set up the source space for my MEG data. This worked fine for all but one subject. For that particular subject, the script does not get past the " Computing patch statistics…" step, instead it just keeps on running for hours (I aborted after 30hours).

I use this code:

from mne import setup_source_space
src = setup_source_space(subject='sub-36', subjects_dir=subjects_dir,
                         spacing='oct6', verbose='DEBUG')

and get following output:

Setting up the source space with the following parameters:

SUBJECTS_DIR = /home/user/projects/dd/lorasick/meg_forward/input_anat/freesurfer
Subject      = sub-36
Surface      = white
Octahedron subdivision grade 6

>>> 1. Creating the source space...

Doing the octahedral vertex picking...
Loading /home/user/projects/dd/lorasick/meg_forward/input_anat/freesurfer/sub-36/surf/lh.white...
Mapping lh sub-36 -> oct (6) ...
    Triangle neighbors and vertex normals...
Loading geometry from /home/user/projects/dd/lorasick/meg_forward/input_anat/freesurfer/sub-36/surf/lh.sphere...
Setting up the triangulation for the decimated surface...
loaded lh.white 4098/137840 selected to source space (oct = 6)

Loading /home/user/projects/dd/lorasick/meg_forward/input_anat/freesurfer/sub-36/surf/rh.white...
Mapping rh sub-36 -> oct (6) ...
    Triangle neighbors and vertex normals...
Loading geometry from /home/user/projects/dd/lorasick/meg_forward/input_anat/freesurfer/sub-36/surf/rh.sphere...
Setting up the triangulation for the decimated surface...
loaded rh.white 4098/140304 selected to source space (oct = 6)

Calculating source space distances (limit=inf mm)...

For other subjects, I additionally get:

    Computing patch statistics...
    Patch information added...
    Computing patch statistics...
    Patch information added...
You are now one step closer to computing the gain matrix

Based on that, I suspect the script freezes somewhere at https://github.com/mne-tools/mne-python/blob/4beb8dde7588c3153ee0a240b5e363dc987c95f1/mne/source_space/_source_space.py#L1600C21-L1600C21

I am on Ubuntu with mne version 1.3

Thanks,
Eduard

britta-wstnr · October 23, 2023, 11:36am

Hi Eduard,

mh, this might actually be hard to tackle without debugging it directly - except someone here has come across a similar problem before.

Some thoughts on what you could try to narrow down the problem (maybe you tried this already):

run it with coarser spacing
have you looked at the surfaces? do they look okay? any holes or weird patches?
and, if you want to look at it in more detail and confirm your suspicion: you could put a breakpoint each before and after the call to the suspicious function and see if it reaches either or both of them before it gets stuck.

HTH,
Britta

eort · October 24, 2023, 3:09pm

Hey Britta,
Thanks for the suggestion!

and, if you want to look at it in more detail and confirm your suspicion:

Yeah, I tried that, and it is indeed this call here, that is the problem

github.com

mne-tools/mne-python/blob/89b5a533867012a122fc43dfe57a18ad530c5906/mne/source_space/_source_space.py#L2743-L2745


      
          d = parallel(
              p_fun(adjacency, s["vertno"], r, dist_limit)
              for r in np.array_split(np.arange(len(s["vertno"])), n_jobs)

have you looked at the surfaces? do they look okay? any holes or weird patches?

No holes, nor patches, but it seems like the pial surface (green) penetrates the skull (gray, see fig). Could this be the problem? And if so, is there a way to fix it?

run it with coarser spacing

This worked. I tried ico4 instead of oct6 and there was no problem. This being said, I don’t think it is ideal to have a subject that is differently processed that all other ones. So preferably, I would find a way to fix the reason for this problem, or am I overreacting?

Unrelated question, when I checked the ico4 source space, I realized that we set up the source space with oct6 spacing, but use ico=4 as parameter in the making of the BEM model. Is this allowed? Or are we mixing now spaces in a weird way? To illustrate, would this be a valid set up of BEM and source space?

model = make_bem_model(ico=4, conductivity=[0.3])
src = setup_source_space(spacing='oct6')

Thanks again,
Eduard

drammock · October 24, 2023, 3:28pm

Have you seen the tutorial about fixing BEM surfaces in Blender? Fixing BEM and head surfaces — MNE 1.5.1 documentation

eort · October 25, 2023, 10:16am

Yes, I have. However I am still a little confused what surface to fix. I mean, from the screenshot, it looks as if I should shrink the pial surface a little so that it fits into the inner skull. The log file of creating the watershed bems and scalp surfaces suggest that all basic topological checks have been passed, except the outer skin surface which is not complete (sum of solid angles yielded 1.000149709292878, should be 1.). Finally, going back to the setup_source_space function, the only surfaces that are used there are the lh/rh.white. The problem persists for the left hemisphere, but visually I can’t detect any obvious problem that I could try to blender-fix.

Nevertheless, I tried fixing the pial surface so that it fits nicely into the skull, still the problem with the source space persists.

Any idea?

britta-wstnr · October 25, 2023, 12:01pm

Hi @eort,

the bad (or good?) news is that this needs to be fixed because otherwise subsequent steps will fail as well. For a BEM model, the boundaries cannot be intersecting but need to be contained within each other. From the picture you sent it looks like you might have a similar “unicorn” situation to the video in the tutorial link @drammock sent you. So you would fix this “unicorn horn” by removing it with Blender (= fixing the pial surface).

I am not quite sure if the watershed algo does check the nested surfaces, if so, the scalp check could point to the same issue. In any case, I would inspect the outer skin again carefully.

Have you looked at the surfaces that are used for the source space setup? You could e.g. check out mne.viz.plot_bem — MNE 1.5.1 documentation

Regarding your question about spacing/sampling: yes you can mix those between the BEM model and the source space. In the first case, this determines how many geometrical elements are used to make the BEM model, in the other case it determines how far it is between two grid points of your source space.

So, in summary: check (confirm you checked) all surfaces and then let’s think from there.

Britta

eort · October 31, 2023, 9:22am

Thanks again Britta!

So you would fix this “unicorn horn” by removing it with Blender (= fixing the pial surface).

I did that, but as I said, at the moment the effect is not really observable, because the source space issue does not seem to be related to this particular issue. Interestingly, I noticed that I cannot directly apply the approach in the tutorial because the pial surfaces does not have nodes that I can select and “push back” into the inner skull. Instead I had to “pull out” the inner skull to contain the pretruding pial surface. Do you think this makes sense?

Have you looked at the surfaces that are used for the source space setup?

Yes, going through the slices one by one, I found one peculiar intersection of the left white and pial (left hemisphere, medial side). Could this be it? If so how would I fix it? In blender, neither pial nor white surface have nodes that I could select and shrink?

Here another view in blender, where I colored the white surface pink, and the pial white. So the pink patches would indicate areas where the surfaces intersect, right?

Eduard

Edit: Okay, now I am entirely confused. I just tried the right hemisphere as a comparison, because for that one there is no problem (as in, setup_source_space.py works properly). However, here the blender check looks much worse:

britta-wstnr · November 2, 2023, 5:14pm

Hi Eduard,

yes, just to make this clear again (I think I was not quite precise in my last answer): the fix of the pial surface was anticipating the next step already (for which this would cause trouble).
Regarding the source space setup: I am not sure what is going on. I wonder if the left and right hemisphere (greatly) overlap and that might cause trouble when setting up the source space? But this is just a blind guess from me given your plots.

Couple of thoughts:

maybe @larsoner has an idea or has seen something like this before?
you could try and see if SimNIBS can deal better with this MRI
you could try and get some deeper insights when exactly the computation spirals (maybe during one of our office hours if you are not sure yourself)
you could use the template MRI for this participant

– I think we first wait if #1, asking @larsoner helps - very often it does

Britta

eort · November 7, 2023, 12:11am

Hi Britta,
Thanks for your input.

overlapping L/R surfaces

For the white surfaces I couldn’t fit an overlap, only the medial walls are essentially touching. However for the pial surfaces, there seem to be overlaps:

.

Besides, in setup_source_space the left and right white surfaces are processed individually (in a loop), starting with the left one. If I change the order, the right surface is processed just fine. Would that be possible if an overlap between the hemispheres causes the issue?

you could try and get some deeper insights when exactly the computation spirals (maybe during one of our office hours if you are not sure yourself)

Okay, I tried to dig down the rabbit hole. It seems the script crashes for one particular vertex. during the Dijkstra algorithm at https://github.com/mne-tools/mne-python/blob/89b5a533867012a122fc43dfe57a18ad530c5906/mne/source_space/_source_space.py#L2793

The vertex in question is this one here (see the red little cross hairs). I don’t see anything that is obviously wrong with it (I plot the left white and pial surfaces, the crosshair is on the white)

Unfortunately, I can’t get more specific that this (bc the source code in scipy is compiled?). In the debugger the last lines I get are these:

> <ipython-input-1-254ab182b81c>(22)_do_src_distances()
-> out = func(con, indices=idx)
(Pdb) s
--Call--
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/scipy/sparse/csgraph/_validation.py(9)validate_graph()
-> def validate_graph(csgraph, directed, dtype=DTYPE,
(Pdb) n
--Return--
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/scipy/sparse/csgraph/_validation.py(56)validate_graph()-><137840x13784...se Row format>
-> return csgraph
(Pdb) 
--Call--
> <__array_function__ internals>(177)any()
(Pdb) 
> <__array_function__ internals>(179)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
> <__array_function__ internals>(181)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
--Return--
> <__array_function__ internals>(180)any()->False
(Pdb) 
--Call--
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/scipy/sparse/_base.py(119)get_shape()
-> def get_shape(self):
(Pdb) 
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/scipy/sparse/_base.py(121)get_shape()
-> return self._shape
(Pdb) 
--Return--
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/scipy/sparse/_base.py(121)get_shape()->(137840, 137840)
-> return self._shape
(Pdb) 
--Call--
> <__array_function__ internals>(177)atleast_1d()
(Pdb) 
> <__array_function__ internals>(179)atleast_1d()
(Pdb) 
> <__array_function__ internals>(180)atleast_1d()
(Pdb) 
> <__array_function__ internals>(181)atleast_1d()
(Pdb) 
> <__array_function__ internals>(180)atleast_1d()
(Pdb) 
--Return--
> <__array_function__ internals>(180)atleast_1d()->array([130388], dtype=int32)
(Pdb) 
--Call--
> <__array_function__ internals>(177)any()
(Pdb) 
> <__array_function__ internals>(179)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
> <__array_function__ internals>(181)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
--Return--
> <__array_function__ internals>(180)any()->False
(Pdb) 
--Call--
> <__array_function__ internals>(177)any()
(Pdb) 
> <__array_function__ internals>(179)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
> <__array_function__ internals>(181)any()
(Pdb) 
> <__array_function__ internals>(180)any()
(Pdb) 
--Return--
> <__array_function__ internals>(180)any()->False
(Pdb) 
--Call--
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/numpy/core/numeric.py(289)full()
-> @set_array_function_like_doc
(Pdb) s
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/numpy/core/numeric.py(337)full()
-> if like is not None:
(Pdb) 
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/numpy/core/numeric.py(340)full()
-> if dtype is None:
(Pdb) 
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/numpy/core/numeric.py(343)full()
-> a = empty(shape, dtype, order)
(Pdb) 
> /home/user/tools/miniconda3/envs/mne/lib/python3.10/site-packages/numpy/core/numeric.py(344)full()
-> multiarray.copyto(a, fill_value, casting='unsafe')
(Pdb) 
-> return a
(Pdb) a
shape = (1, 137840)
fill_value = inf
dtype = <class 'numpy.float64'>
order = 'C'
like = None
(Pdb) fill_value

As soon as I continue from here, I get stuck. So, I guess the problem occurs after. When I try this debugging for a vertex that has no issue, I don’t see any difference up to the point before it starts to spiral.

Is this helpful in any way? Or should I just give up that subject?

Thanks,
Eduard

larsoner · November 15, 2023, 4:27pm

I don’t think the overlap with the skull or overlap between left and right hemi pial surfaces will make a difference for your patch stats / distance calculations.

Assuming your individual surface topology is okay (no holes, etc.) then in principle the add_dist=True should only take a couple of minutes to complete. The fact that it gets “stuck” for you suggests to me that there is some bug in scipy.spatial.distance.dijkstra, which is what we call under the hood. Any chance you can share the surf directory for this subject with me? Then I’d be able to replicate locally and open an upstream bug.

In the meantime you could try add_dist='patch', which does far fewer calculations, but would still result in the same improved source normals and thus forward/inv op you’d get with add_dist=True. add_dist=True is really mostly useful for when you need distances between the source points, like in clustering. So usually it’s most useful for fsaverage or whatever you eventually morph all subject data to.

eort · November 15, 2023, 8:19pm

Hi @larsoner,

Thanks for the feedback, indeed, setting add_dist='patch' solved the stuckiness.

Thanks again!
Eduard

edit: Link removed

larsoner · November 16, 2023, 2:50pm

So I downloaded the data and unzipped it in a temp folder and ran:

$ time python -c 'import mne; mne.setup_source_space("temp", subjects_dir=".", verbose="debug")'
...
real	5m57.213s

so it only took about 5 minutes on my system. No idea why it’s broken on yours. SciPy does the heavy lifting here, this is my relevant sys_info:

Platform             Linux-6.5.0-10-generic-x86_64-with-glibc2.38
Python               3.11.6 (main, Oct  8 2023, 05:06:43) [GCC 13.2.0]
Executable           /home/larsoner/python/virtualenvs/base/bin/python3
CPU                  x86_64 (8 cores)
Memory               62.7 GB

Core
├☑ mne               1.6.0.dev206+g239e574ee2.d20231116 (devel, latest release is 1.5.1)
├☑ numpy             1.26.1 (OpenBLAS 0.3.23.dev with 4 threads)
├☑ scipy             1.12.0.dev0+1888.d7f12d3

I doubt it matters that I’m on a dev version of SciPy but it could. If you want to try it you can install one for example with:

pip install --upgrade --index-url "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple" scipy

eort · November 16, 2023, 8:03pm

That’s weird. Your way didn’t work (I got some Freesurfer related problem), but when I installed the latest mne and reran my script, it worked fine. So yeah, it must be to do with the versions (I used mne 1.3, scipy 1.10).

Still it is interesting that this problem only occurred for one out of 60 subjects.
Well thanks for looking into this!
Eduard

Topic		Replies	Views
mne_setup_source_space ✉️ Mailing List Archive (read-only) list-archive	3	146	June 5, 2013
Help with the 'spacing' parameter in mne.setup_source_space 💬 Support & Discussions source-localization , eeg	2	308	October 29, 2021
Understanding source spaces 💬 Support & Discussions source-localization , visualization	1	334	April 14, 2021
plot_label_from_stc.py: Exception: The source space does not appear to be an ico surface ✉️ Mailing List Archive (read-only) list-archive	6	176	January 6, 2014
a junior question on mne_setup_source_space ✉️ Mailing List Archive (read-only) list-archive	32	1037	May 1, 2014

Computing source space does not complete

Related topics