Forward object contains triangulation which uses omitted sources

[makkostya]

I first computed source space for sample subject: src = mne.setup_source_space(subject, spacing='ico5', subjects_dir=subjects_dir, add_dist=False) Then forward solution:

model = mne.make_bem_model(subject=subject, ico=None,
                           subjects_dir=subjects_dir)
bem = mne.make_bem_solution(model)
fwd = mne.make_forward_solution(evoked.info, trans=trans, src=src, bem=bem,
                                meg=True, eeg=True, mindist=3.0, n_jobs=2)

Some source points were out of distance limits, so were deleted from the model:

9 source space points omitted because they are outside the inner skull surface. 114 source space points omitted because of the 3.0-mm distance limit.

But the forward object still contains triangulation which uses these omitted sources

print src[0]['nuse'], src[0]['nuse_tri']
print fwd['src'][0]['nuse'], fwd_fixed['src'][0]['nuse_tri']

10242 [20480] 10119 [20480]

After discussing this issue with @agramfort , a solution could be not to omit these sources, but to interpolate their lead fields using morph_data tool. I can do a PR for this, adding this option to mne.make_forward_solution() function.

[larsoner]

interpolate their lead fields using morph_data tool

Has anyone checked that the results are sensible? It should be pretty straightforward -- compute forward with some dist limit that works for all sources, then compute it with whatever dist limit causes exactly one source to be omitted, and check your interpolated value against the veridical one. I'm not sure it will be all that close.

To decide how close is "close enough" you could look at how similar the MEG lead fields are for using a 1-layer vs 3-layer BEM, which presumably should not be too different (right @agramfort?).

But the deeper question is: why do you need the triangulation information to be complete? If that is all that you really need (and you don't actually need the source point itself), I'm sure there are ways to "re-close" the mesh. For one, you could try using the spherical surface to perform a Delaunay of the points that formerly touched the removed vertex.

[agramfort]

usecase is for example to make sure we have the same vertno for all subjects or if you need a closed mesh for data-driven parcellation computation

in terms of accuracy is will not be correct (think that radial dipoles close to inner skull cannot be guessed by interpolation) but it should good enough for some use cases.

[larsoner]

I'd like to see the result of the suggested test. It can work as a unit test anyway so it's useful for multiple reasons

[makkostya]

think that radial dipoles close to inner skull cannot be guessed by interpolation

I don't think that orientation meters a lot, cause anyway x, y, z directions are estimated, not the fixed orientation lead field. For me personally, I would prefer to not interpolate but to omit these points and then reclose the mesh. But it's harder to implement (at least me, I don't know how to do it), so the proposed solution is just the simplest one.

[larsoner]

I don't think that orientation meters a lot, cause anyway x, y, z directions are estimated, not the fixed orientation lead field.

Yes but think about it differently -- we compute in cardinal XYZ simply for convenience. What we calculate is equivalent (through mult with a rotation matrix) to computing a radial component and two other components orthogonal to it. So @agramfort is saying one of these three degrees of freedom is likely to be pretty bad.

How bad can be established with the little test I suggested above (which we'll need anyway as a unit test). So hopefully you can implement it and see what the quality is like!

I would prefer to not interpolate but to omit these points and then reclose the mesh

This actually is not so difficult a problem because by Freesurfer we have a way to map all our cortical vertices onto a sphere (preserving the existing triangulation). But it does not solve as many problems as interpolation would.

Realistically we might consider having the API on_missing=:

1. 'info': current behavior, inform the user but don't do anything special
2. 'warn': allow it to warn instead
3. 'raise': raise RuntimeError
4. 'interpolate': initial suggestion
5. 'close': just close/complete the mesh after omission

At least the first 4 make for a hopefully simple PR. I can add the "close" option if you'd like to explore how useful (or not) it is.

[agramfort]

just give it a try and look at fwd sensitivity maps to see how it looks like

[makkostya]

ok, i'll start with this