Hi all,
Is there any formula to find out the surface area corresponding to a dipole within a source space? When I am looking at the MNE documentation regarding ico-5 source space, its written as:
Source space: ‘ico5’
Sources per hemisphere: 10242
Source spacing: 3.1 mm
Surface area per source : 9.8 mm2 (is there any general formula to calculate this?)
Thanks,
@larsoner do you have any idea about it? thanks!
There is this info:
>>> import mne
>>> src = mne.read_source_spaces(mne.datasets.sample.data_path() / 'subjects' / 'fsaverage' / 'bem' / 'fsaverage-ico-5-src.fif', patch_stats=True)
>>> src[0].keys()
dict_keys(['id', 'type', 'np', 'ntri', 'coord_frame', 'rr', 'nn', 'tris', 'nuse', 'inuse', 'vertno', 'nuse_tri', 'use_tris', 'nearest', 'nearest_dist', 'pinfo', 'patch_inds', 'dist', 'dist_limit', 'subject_his_id', 'tri_area', 'tri_cent', 'tri_nn', 'use_tri_cent', 'use_tri_nn', 'use_tri_area'])
>>> src[0]["use_tri_area"].shape
(20480,)
>>> src[0]["use_tris"].shape
(20480, 3)
So this could get you the surface area per triangle in the decimated space. But you want it per vertex not per triangle.
Really the “surface area per source” can be thought of as “total surface area / total number of sources”. So looking at it that way, you could calculate it as:
sum(s["tri_area"].sum() for s in src) / sum(s['nuse'] for s in src) * 1e6
6.367790682029399
The * 1e6 is a conversion from “sources per square m” to sources per square mm". For an oct-6 for example you’d get:
>>> src = mne.read_source_spaces(mne.datasets.sample.data_path() / 'subjects' / 'sample' / 'bem' / 'sample-oct-6-src.fif', patch_stats=True)
>>> sum(s["use_tri_area"].sum() for s in src) / sum(s['nuse'] for s in src) * 1e6
20.69581503362706
I’m not sure why the value for the ico-5 doesn’t match, assuming that manual entry is for fsaverage (because it will vary by subject). Even that for sample wouldn’t be 9.8:
>>> sum(s["use_tri_area"].sum() for s in src) / 20484 * 1e6
8.280750830677961
… oh in those last two I errantly used use_tri_area – if we use tri_area we do get a value that matches (if you floor instead of round):
>>> sum(s["tri_area"].sum() for s in src) / 20484 * 1e6
9.882712190552494
@larsoner Thanks a lot for the explanation.
It works with sample data but not the 'fsaverge'
In [38]: src= read_source_spaces('/home/dip_meg/mne_data/MNE-sample-data/subjects/fsaverage/bem/fsaverage-ico-5-src.fif')
Reading a source space...
[done]
Reading a source space...
[done]
2 source spaces read
In [39]: src
Out[39]: <SourceSpaces: [<surface (lh), n_vertices=163842, n_used=10242>, <surface (rh), n_vertices=163842, n_used=10242>] MRI (surface RAS) coords, subject 'fsaverage', ~21.9 MB>
In [40]: sum(s["tri_area"].sum() for s in src) / 20484 * 1e6
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[40], line 1
----> 1 sum(s["tri_area"].sum() for s in src) / 20484 * 1e6
Cell In[40], line 1, in <genexpr>(.0)
----> 1 sum(s["tri_area"].sum() for s in src) / 20484 * 1e6
KeyError: 'tri_area'
In [41]: sum(s["use_tri_area"].sum() for s in src) / 20484 * 1e6
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[41], line 1
----> 1 sum(s["use_tri_area"].sum() for s in src) / 20484 * 1e6
Cell In[41], line 1, in <genexpr>(.0)
----> 1 sum(s["use_tri_area"].sum() for s in src) / 20484 * 1e6
KeyError: 'use_tri_area'
It seems like use_tri_area or tri_area attributes are missing in the src. Am I missing something?
thanks!
You have to pass patch_stats=True to the read_source_spaces call
I can print the values now but for ico-5 for fsaverage, it doesn’t give the same value of 9.8 as documented in MNE. Would that mean the value is rather based on the sample subject src and not for fsaverage?
In [3]: from mne import read_source_spaces
In [4]: src= read_source_spaces('/home/dip_meg/mne_data/MNE-fsaverage-data/fsave
...: rage/bem/fsaverage-ico-5-src.fif', patch_stats=True)
Reading a source space...
[done]
Completing triangulation info...
[done]
Completing selection triangulation info...
[done]
Reading a source space...
[done]
Completing triangulation info...
[done]
Completing selection triangulation info...
[done]
2 source spaces read
In [5]: src
Out[5]: <SourceSpaces: [<surface (lh), n_vertices=163842, n_used=10242>, <surface (rh), n_vertices=163842, n_used=10242>] MRI (surface RAS) coords, subject 'fsaverage', ~59.1 MB>
In [6]: src.info
Out[6]:
{'fname': '/home/dip_meg/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-ico-5-src.fif',
'working_dir': '/autofs/space/megmix_002/users/mluessi/data/code/mne-scripts',
'command_line': 'mne_make_source_space --surf lh.white:rh.white --ico 5 --src /cluster/fusion/mluessi/MNE-sample-data-build/MNE-sample-data/subjects/fsaverage/bem/fsaverage-ico-5-src.fif'}
In [7]: sum(s["tri_area"].sum() for s in src) / 20484 * 1e6
Out[7]: 6.367790682029399
In [8]: sum(s["use_tri_area"].sum() for s in src) / 20484 * 1e6
Out[8]: 6.246565416014217
I guess it must be. It’ll depend on the brain surface area since the number of vertices for ico-5 is always the same…
Yes, you are right. I looked it up with my subjects data.
thanks for helping.
@larsoner sorry, I have a follow up question.
If I look at the mne.setup_source_space — MNE 1.5.0 documentation , MNE creates (at-least with default setup) source space at the interface between gray and white matter, i.e. at the white surface (decimated surface, more suitable for morphological feature)? and not at the middle surface, i.e. a surface that runs at the mid-distance between white and pial??
I am looking at the mne code to compute Cortical surface area for the triangulation.
If I understand the code properly, then this means (I put the comments at the side):
def _complete_source_space_info(this, verbose=None):
"""Add more info on surface."""
# Main triangulation
logger.info(" Completing triangulation info...")
this["tri_area"] = np.zeros(this["ntri"])
r1 = this["rr"][this["tris"][:, 0], :]
r2 = this["rr"][this["tris"][:, 1], :]
r3 = this["rr"][this["tris"][:, 2], :]
this["tri_cent"] = (r1 + r2 + r3) / 3.0
this["tri_nn"] = fast_cross_3d((r2 - r1), (r3 - r1))
this["tri_area"] = _normalize_vectors(this["tri_nn"]) / 2.0
logger.info("[done]")
# Selected triangles
logger.info(" Completing selection triangulation info...")
if this["nuse_tri"] > 0:
r1 = this["rr"][this["use_tris"][:, 0], :] >>>>>>>>> a triangular face say, r1r2r3; vertex coordinates
r2 = this["rr"][this["use_tris"][:, 1], :] >>>>>>>>> vertex coordinates
r3 = this["rr"][this["use_tris"][:, 2], :] >>>>>>>>> vertex coordinates
this["use_tri_cent"] = (r1 + r2 + r3) / 3.0 >>>>>>>>> Conversion from facewise to vertexwise each vertex one-third of the sum of the areas of all faces that meet at that vertex
this["use_tri_nn"] = fast_cross_3d((r2 - r1), (r3 - r1)) >>>>>>>>>>> computes the triangle vector normal
this["use_tri_area"] = np.linalg.norm(this["use_tri_nn"], axis=1) / 2.0 >>>>>>>>>>> computes The triangle areas from vector normal in m2
logger.info("[done]")
Am I correct? I am following this paper : https://s3.us-east-2.amazonaws.com/brainder/publications/2012/winkler2012_facewise_surface_area.pdf to understand the logic.
Highly appreciate your help.
Dip
That sounds right, the white surface should be the white/gray interface
I don’t have the bandwidth to “check your work” so-to-speak on this – it’s probably better to convince yourself using more reading and/or checking against other libraries (like vtk) anyway!
Thanks! I will do that.