Codecov Report

Merging #7639 into master will increase coverage by 0.04%.
The diff coverage is 95.38%.

@@            Coverage Diff             @@
##           master    #7639      +/-   ##
==========================================
+ Coverage   90.31%   90.35%   +0.04%     
==========================================
  Files         455      455              
  Lines       83999    84427     +428     
  Branches    13302    13370      +68     
==========================================
+ Hits        75860    76285     +425     
- Misses       5270     5272       +2     
- Partials     2869     2870       +1     

Okay @agramfort this is ready for review/merge as well. Might be good for you to take a first pass, then once we're happy with the API ask others who would be interested to look

The new tutorial image showing extract_label_time_course on a volume STC:

I was trying to read the volume stc which were computed previously using mne.read_source_estimate in mne.v0.21.dev0.
I met with an error
TypeError: vertices must be an instance of list, got <class 'numpy.ndarray'> instead.
without reading stcs i cannot check the extract_label_time_course ,correct @larsoner ?

@poojaprabhu9 I pushed a fix that should fix that bug hopefully

I tried the same by installing mne.v0.21.dev0.
I error still exist.

that helped, thanks.
I had computed stc after morphing.
Thus, now when i am extracting the time course, i have to setup the volume source space for fsaverage to get ''src" (same src was used while doing morphing)

src_fs = mne.setup_volume_source_space( 'fsaverage', pos=5.,mri=aseg_fname, 
             subjects_dir=subjects_dir,verbose=True)

then i tried to extract the time course from the stc

tc=stc_morph.extract_label_time_course(labels='/path/fsaverage/mri/aparc+aseg.mgz',src=src_fs,mode='mean',allow_empty=True)

i am encountering the unexpected error:
ValueError: 9900/20377 left hemisphere stc vertices missing from the source space, likely mismatch.

How is that possible when i am using the same src to extract the time course as well as to compute the morphed stc?

I tried to replicate the code here
two observations are done,

1. when i execute the code in mne.v.0.20, the stc obtained after morphing has 20377 vertices whereas when executed in development version(pip install -U https://github.com/larsoner/mne-python/archive/atlas.zip ) it has 24303 vertices. I can extract the time course from the morphed stc computed from development version. In my case i have the morphed stc computed using mne.v.0.20, when i read this in development version and try to extract the time course i get the ValueError: 9900/20377 left hemisphere stc vertices missing from the source space, likely mismatch.
2. when i try to execute (line 49 of gist)
   src = mne.setup_volume_source_space('sample', mri=aseg_fname,pos=5., sphere=sphere, mindist=3.,subjects_dir=subjects_dir)
   in development version i get
   MemoryError: Unable to allocate 402. TiB for an array with shape (48003, 48003, 48003) and data type int32
   whereas it works fine in mne.v.0.40.
   so you can't run the whole code in development version. So i computed the forward solution in mne.v.0.20 and then i read the forward solution in development version.

You have to be careful about sphere_units, they changed between versions. You probably need to set this to mm, the default is now m.

In general if you stick to using dev versions you need to follow changes or read the release notes very closely because deprecation cycles can be effectively days or weeks in duration instead of 6 months.

You have to be careful about sphere_units, they changed between versions. You probably need to set this to mm, the default is now m.

In general if you stick to using dev versions you need to follow changes or read the release notes very closely because deprecation cycles can be effectively days or weeks in duration instead of 6 months.

that solved the second point.
first problem is solved by adding sphere=(0,0,0,90) to

src_fs = mne.setup_volume_source_space( 'fsaverage', pos=5.,mri=aseg_fname, 
             subjects_dir=subjects_dir,verbose=True)

Then compute the extract_time_course.
I could extract the time course. Thanks!

Yes right now, because everything is done in Freesurfer surface RAS coordinates (vox_mri_t in MNE, get_vox2ras_tkr in nibabel) and this doesn't exist AFAIK for NIfTI. Right now you have to make sure that whatever atlas image you use has been transformed properly to the subject's MRI space so that things like mri_resolution=True work properly, and that a given voxel in your atlas corresponds to the same voxel in your T1.mgz. The builtin Freesurfer aseg.mgz meets these requirements.

Have you already transformed the HCPMMP1 parcellation to your subject's MRI space such that a given voxel in your HCPMMP1 volume corresponds to the same voxel in their aseg.mgz / T1.mgz? If so, how did you do this transformation (we should document it)? If on the other hand you just have HCPMMP1 in NIfTI for fsaverage, we'll need to figure out how to transform it to the subject's space in the first place. It should be doable, but we need to figure it out. Once we figure that out, we can document it and then figure out what needs to be done at the MNE end to facilitate using other atlases.

Check this link: https://cjneurolab.org/2016/11/22/hcp-mmp1-0-volumetric-nifti-masks-in-native-structural-space/
As per the link, the code to create volume parcellation using HCPMMP1 is here (this code is from the same link). the final volume will be in .nii format.
We just need it for fsaverage brain.
We need to apply this parcellation on the already morphed volume estimates. For this purpose we have to extract the time courses for each labels of HCPMMP1 parcellation. Currently, to extract the time courses we have the atlas in nifti. So the previous question.

So fsaverage is pretty nice:

>>> import nibabel as nib
>>> img = nib.load('/mnt/bakraid/data/structurals/fsaverage/mri/aseg.mgz')
>>> import numpy as np
>>> np.array_equal(img.affine, img.header.get_vox2ras_tkr())
True

You don't have to worry about these differences. So in theory if you just save your nii to mgz, the MNE code should work, because (I think) nibabel will write the vox2ras_tkr as being equivalent to vox2ras / affine.

... and personally I would probably update that script you used to output .mgz files with the vox2ras_tkr set properly. Then it would even work for non-fsaverage subjects.

You mean to say I can use atlas in .nii format to extract time courses from stc?
(Without updating the MNE)

Or you mean I need to convert .nii to .mgz? Then try to extract the time course

Or you mean I need to convert .nii to .mgz? Then try to extract the time course

You almost certainly need to convert it to mgz. Even just doing mri_convert /path/whatever.nii /path/whatever.mgz might be enough for things to work.

We can probably ship some HCPMMP1 + aseg volumes for fsaverage at least. I think that doing the conversion EDIT: for other subjects is out of scope for us. But it sounds like having the volumes for fsaverage is already enough, assuming you're morphing all subject data there anyway (?).

... looking into this a bit, the solution they posted is a start but it's not exactly straightforward to get it working. I trimmed the script a bit just to use with fsaverage and output an mgz file. However, it does not output the actual region LUT that is to be used with the volume, so this needs to be created based on the math it actually does. And doing these things in Bash is not exactly straightforward.

I think it would actually be easier to do most of this in Python with nibabel. All that we really need FreeSurfer for (I think) is to generate a aseg volume from the annotation:

mri_aparc2aseg --s fsaverage --o HCPMMP1.mgz  --annot HCPMMP1

From there in Python we can:

1. Load HCPMMP1.mgz
2. Remove hippocampal voxels (they do this in the script)
3. Shift the values by some amount (they do this in the script)
4. Load the standard aseg.mgz for the subject, and add it to the HCPMMP one (including only whatever regions should be included)
5. Keep track of which vertex values correspond to which regions, and then write out a custom LUT.

Eventually this could be done for other subjects, too, for example by using read_labels_from_annot followed by morph_labels then write_labels_to_annot for each subject (they do this sort of operation in their bash script).

For now, perhaps try using aparc.a2009s+aseg.mgz. Each FreeSurfer-reconstructed subject should already have this. It's a much less challenging path, and you should already have the files you need. In principle later HCPMMP1 can be substituted as the atlas+LUT combination, but in the meantime you can set up your scripts to work without this being a bottleneck.

@poojaprabhu9 @virvw can you try this FreeSurfer command (after running mne.datasets.fetch_hcpmmp_parcellation):

$ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1
$ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1_combined

Then run this gist, which will create SUBJECTS_DIR/fsaverage/mri/HCPMMP1ColorLUT.txt and HCPMMP1_combinedColorLUT.txt. You can then view these with FreeSurfer (from SUBJECTS_DIR/fsaverage/mri directory:

$ freeview -v HCPMMP1+aseg.mgz:colormap=lut:lut=HCPMMP1ColorLUT.txt

freeview -v HCPMMP1_combined+aseg.mgz:colormap=lut:lut=HCPMMP1_combinedColorLUT.txt

Then in MNE you can do something like the following to extract label time courses:

mri_dir = op.join(mne.get_config('SUBJECTS_DIR'), 'fsaverage', 'mri')
labels = mne.read_labels_from_annot('fsaverage', 'HCPMMP1', 'both')
lut = mne.read_freesurfer_lut(op.join(mri_dir, 'HCPMMP1ColorLUT.txt'))
# limit to just the HCPMMP labels with this, or omit it to keep subcortical ones:
lut = {label.name: lut[label.name] for label in labels}
atlas = op.join(mri_dir, 'HCPMMP1+aseg.mgz')
label_tc = mne.extract_label_time_course(stc, (atlas, lut), fwd['src'])

I was trying to use the above code but unfortunately i could not read the stc file (i updated the mne using, pip install https://api.github.com/repos/mne-tools/mne-python/zipball/master) because of this error.
TypeError: vertices must be an instance of list, got <class 'numpy.ndarray'> instead
I remember this problem was solved (when i was using, pip install -U https://github.com/larsoner/mne-python/archive/atlas.zip)
Did i miss something?

pip install https://api.github.com/repos/mne-tools/mne-python/zipball/master

This might not have done anything, try doing it again with --upgrade in there. If it still fails, upload the STC somewhere and I can see if I can load it

I tried upgrading by
pip --upgrade install https://api.github.com/repos/mne-tools/mne-python/zipball/master
Still the error exist.

You can find the stc file here
Please try stc.in_label and stc.extract_time_courses.

Can you upload the source space that you used? I can't try those functions without it

just run this line

src = mne.setup_volume_source_space(
                        'fsaverage', pos=5., subjects_dir=subjects_dir, sphere=(0,0,0,90),
                        verbose=True,sphere_units='mm')

Okay -- it's better to use the bem/brain to define the set of sources to keep rather than sphere

Okay -- it's better to use the bem/brain to define the set of sources to keep rather than sphere

will try that for the further computations.

@larsoner @poojaprabhu9
This worked:
$ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1
$ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1_combined

I am stuck with the gist for several reasons independent of the code itself
@poojaprabhu9 if you manage to run that bit on the fsaverage and forward the outcome, I'll finish up...

@larsoner @poojaprabhu9

The gist was ran removing the initial option "sort" because not available in our version.

Note that the created file "FreeSurferColorLUT.txt" was found by the script "create_subj_volume_parcellation.sh" when placed in the same folder as the script and not in the "fsaverage/label/"

Now I am encountering issue when running "create_subj_volume_parcellation.sh" since we are missing "wm.mgz" for fsaverage

FreeSurfer/vv221713@is152869:/neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI$ bash create_subj_volume_parcellation.sh -L subject_list.txt -a HCPMMP1 -d HCPMMP1_parcellation -s YES -t YES

     >>>>         Current FreeSurfer subjects folder is /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI

reading colortable from annotation file...
colortable with 181 entries read (originally ./left.fsaverage164.label.gii)
reading colortable from annotation file...
colortable with 181 entries read (originally ./left.fsaverage164.label.gii)
reading colortable from annotation file...
colortable with 181 entries read (originally ./right.fsaverage164.label.gii)
reading colortable from annotation file...
colortable with 181 entries read (originally ./right.fsaverage164.label.gii)

     >>>>         PREPROCESSING fsaverage         <<<<
     >>>>         START TIME: Sun May 10 16:49:56 CEST 2020         <<<<

Annotation files lh.fsaverage_HCPMMP1.annot and rh.fsaverage_HCPMMP1.annot already exist in fsaverage/label. Won't perform transformations
reading colortable from annotation file...
colortable with 182 entries read (originally HCPMMP1_parcellation/temp_1_0_1246/colortab_HCPMMP1_L)
reading colortable from annotation file...
colortable with 182 entries read (originally HCPMMP1_parcellation/temp_1_0_1246/colortab_HCPMMP1_R)
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
Creating subcortical aseg masks for subject fsaverage
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
fslmaths: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
computing statistics for each annotation in /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/label/lh.fsaverage_HCPMMP1.annot.
reading volume /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz...
mghRead(/neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz, -1): could not open file
mris_anatomical_stats: could not read input volume /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz
computing statistics for each annotation in /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/label/rh.fsaverage_HCPMMP1.annot.
reading volume /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz...
mghRead(/neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz, -1): could not open file
mris_anatomical_stats: could not read input volume /neurospin/meg/meg_tmp/2015_MiniManip_rsMEG/POLTI_Ignacio_RSretroTIME/MRI/fsaverage/mri/wm.mgz
(standard_in) 2: syntax error
sed: -e expression n°1, caractère 3: ,' inattendue sed: -e expression n°1, caractère 3: ,' inattendue
sed: -e expression n°1, caractère 3: ,' inattendue sed: -e expression n°1, caractère 3: ,' inattendue

     >>>>         fsaverage STARTED AT Sun May 10 16:49:56 CEST 2020, ENDED AT: Sun May 10 16:51:04 CEST 2020

The gist was ran removing the initial option "sort" because not available in our version.

It will not work properly without that option. The regions will be mislabeled. It was added to master recently, so if you update it should be available.

Note that the created file "FreeSurferColorLUT.txt" was found by the script

You should not run this at all anymore, just the two steps mentioned in the comment above:

1. $ mri_aparc2aseg --s fsaverage --volmask --annot HCPMMP1
2. python hcpmmp_lut.py from here

This should create all the files for fsaverage that you need to be able to extract label time courses (i.e., run the Python code from this comment).

There is something missing in the line 17 of hcpmmp_lut.py. is it so?

That's the start of a multi-line string, looks okay to me. Did you try and and get some error?

considering thet the sort option is set to FALSE in your script,
does that matter at all?

In master the default is True (and anything before master implicitly had sort=True), so yes unfortunately you need it.

Also, if you plan to use the HCPMMP1_combined version, you probably need to (for now at least) regenerate it using latest master, which you can do by deleting the lh/rh.HCPMMP1_combined.annot files from fsaverage/label the doing mne.datasets.fetch_hcpmmp_parcellation(). And if you already did the FreeSurfer mri_aparc2aseg command, redo it as well, and then re-run my gist. This also has to do with the sort option defauting to True, which causes problems when generating the low-res/combined parc...

That's the start of a multi-line string, looks okay to me. Did you try and and get some error?

Everything is fine. It was a mistake from my side.

@larsoner @virvw
It is working!!
Thanks a lot @larsoner