Integration with Brainstorm

[christian-oreilly]

Initial work on this topic has been discussed here: https://neuroimage.usc.edu/forums/t/atlas-for-1-year-old-infants/15475/71

There are currently some issues with the co-registration of the montage and the surfaces in Brainstorm. From previous email exchange with @ftadel and @larsoner

Hi Eric & Christian,

As I wrote in the forum thread (https://neuroimage.usc.edu/forums/t/atlas-for-1-year-old-infants/15475/70):
The coregistration in Brainstorm is independent from the original coordinate system of either the surfaces 
or the electrodes. Everything imported to Brainstorm is converted to CTF/SCS coordinate system using 
the NAS/LPA/RPA (https://neuroimage.usc.edu/brainstorm/CoordinateSystems), so as long as these points 
are correctly defined in the MRI and in the montage file (ie. in the same coordinate system as the rest of 
the 3D points in the same file), there is no possible coordinate system issue once imported in Brainstorm.
=> You can use the coordinate system you want in the montage.fif files.

The current issues (documented in the forum thread):

1. The bem/outer_skin surface is too large wrt to the other bem surfaces and the MRI
2. The coordinates in the .fif files are not labelled (impossible to know which electrode they 
correspond to)  => It would be much easier and efficient to transfer these files between programs 
if they were simple text files, one 3D point per line, "Name X Y Z" (and you can add an extra 5th column 
Type if you can't derive it from the name of the point). Like the Polhemus files... all the complexity of the 
FIF files does not bring anything good here, especially if the coordinate system is not even the one 
documented in the FIF structure.
    If you want something truly portable, use the BIDS electrodes.tsv specification, it would be good 
for you to have an export function to this format in any case:
    https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/03-electroencephalography.html#electrodes-description-_electrodestsv
3.    The NAS/LPA/RPA positions for the import of the MRI+surface must be defined in 
Brainstorm MRI coordinates (this is only annoying part that might require extra manual work)

Could we keep this discussion going somewhere it's archived, organized and accessible 
by others - instead of via  emails? (either the Brainstorm forum thread referenced above or a
new brainstorm or MNE github issue)

Cheers
Francois

[christian-oreilly]

I have noticed some compatibility issues between the code exporting the montage and the code importing it for Brainstorm. I am looking into it.

Regarding BIDS compatibility, I am all for in theory. In practice, it is part of the BIDS ecosystem and I am not convinced it is not interdependent with the other expected files for BIDS projects such as the _coordsystem.json and _T1w.json for the fiducials. I think BIDS is relatively new an I am not sure I want to get into writing import/export functions for software packages like Brainstorm, MNE, and FieldTrip if these functionalities are not already in place. So, for the time being, I'll fix the incompatibilities I said I have noticed and check if it solves the issue at hand.

[christian-oreilly]

@larsoner The .fif export function mne.io.meas_info.write_dig(montage_fname, montage.dig) takes directly a list of mne.io._digitization.DigPoint objects which are essentially dictionaries with keys ['r', 'ident', 'kind', 'coord_frame'] and with corresponding values e.g. like [array([0.02494444, 0.03210928, 0.05049407]), 1 (FIFFV_POINT_LPA), 1 (FIFFV_POINT_CARDINAL), 5 (FIFFV_COORD_MRI)]. I see no clear way to save an actual name for the individual channels/fiducials. Am I right to understand that this information cannot be saved in the FIF format (at least with the actual MNE-Python functions)?

[larsoner]

Yes I would only use that function for MRI fiducials SUBJECT/bem/SUBJECT-fiducials.fif for MNE. I like @ftadel's suggestion of just using some suitable text format for the EEG locations (probably Polhemus-like?) for the EEG montages, knowing that they are in FreeSurfer surface RAS coordinates (if relevant to people).

[christian-oreilly]

OK. I'll go with the popular vote and prepare a Polhemus-like format then.

[christian-oreilly]

@ftadel @larsoner Does one of you know where I can find a description of that format and/or example files? I know that it is a simple format but there are still things that need to be specified, e.g., is there a title row? What is the unit? What is the coordinate system, e.g. LPI? What is the columns order? What is the spacing separator (space vs tabs)? And so on...

[ftadel]

The BIDS electrodes.tsv and coordsystem.json are a robust way of documenting both the electrodes positions and the referential. These files are (or will soon be) readable by all the programs dealing with EEG. They are simple text files, easy to edit with any text editor, human readable and editable. Using these formats doesn't mean you have to use a full BIDS structure, or that you need to document your MRI in BIDS either.

https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/03-electroencephalography.html#electrodes-description-_electrodestsv https://bids-specification.readthedocs.io/en/stable/04-modality-specific-files/03-electroencephalography.html#coordinate-system-json-_coordsystemjson

You can find some examples here: https://github.com/bids-standard/bids-examples

[ftadel]

(skip the coordsystem.json if it doesn't make much sense to you, and just add to the filename of the tsv a tag indicating in which coordinate systems they are supposed to be)

[christian-oreilly]

Here is what I have right now for e.g. the 10-20 montage (10-20_electrodes.tsv)

name    x   y   z
FZ  -0.003055893452212003   -0.032064090002355834   0.06451582182251714
CZ  -0.00296771730663196    -0.06977075333711903    0.04008209482506764
PZ  -0.0021336846531962396  -0.07582308854450227    -0.006966221879654919
OZ  -0.017715934930783665   0.04168709100177952 0.030925180428014887
FP1 0.014960283043851674    0.041850759295829036    0.03303760429864649
FP2 -0.030231570180668153   0.014258093537300384    0.047310049616561836
F3  0.027036921474096004    0.015378651880977595    0.05032452361776425
F4  -0.040706601439877865   0.022890903430381613    0.015435878672161086
F7  0.04102759027159819 0.024366703496549222    0.02106096110564212
F8  -0.04276038505091253    -0.02259439469485351    0.04209152095672312
C3  0.040038080118693   -0.024764628668414924   0.046762038050382405
C4  -0.04753591853871812    -0.004965693658377632   -0.0031449540422044207
T7  0.048257535492658744    -0.0068005557623915655  0.0011775356013593006
T8  -0.03434671963946244    -0.060219449082566344   0.024608599920027285
P3  0.03082875147042082 -0.06353352752464983    0.02870305745408268
P4  -0.04232946702432042    -0.040302583452643655   -0.0064429931472090125
P7  0.04340974338795373 -0.043749768222233316   -0.0032022100664471412
P8  -0.02027696559388767    -0.07056118236440255    -0.007539028015980339
O1  0.017546746740611636    -0.07282044712298018    -0.006041763414346149

and the accompanying fiducial file (fiducials.tsv)

name    x   y   z
LPA 0.023432455652875273    0.02951177467802661 0.04414493514701666
Nz  -0.04705270613411553    -0.025418123431664483   -0.0006206020544484739
RPA -0.0014942770523743206  0.04826398113209235 0.011539915674043048
AC  -0.0005 -0.001  0.0125
PC  -0.0415 0.003   -0.0235
Iz  0.04177321890348855 0.0019660082100442173   -0.021722263403512024
LMA -0.041832591573434816   -0.021128485096939095   -0.02466985490888543
RMA 0.04216618983303278 -0.021879278120898776   -0.021371300356866553
LEOG    -0.0014970128735639107  0.05121324257406698 0.004655698015611161
REOG    -0.034005117804820816   0.03742933393647815 0.002670292182298529
LEMG    0.030949635606162908    0.03944971691571252 0.002546296565711347
REMG    -0.021617247115082894   0.0495608026706042  -0.007285594184332787
mid_EOG_nose    0.018263009394794506    0.048814659668375976    -0.004625029842968209
Cz  -0.001549808042629331   -0.014931175460064486   0.06593069103913642
IH  2.4995  -20.001 35.5125

FYI, the apparent difference in precision for AC, PC and IH with respect to the other fiducials is not a bug. This is due to the other fiducial being slightly moved from their NMD position so that they fall precisely on the BEM surface, whereas the "internal" fiducials (AC, PC, IH) are obviously not corrected in that fashion.

@ftadel @larsoner Are we good to go with these file formats? These will supplement, not replace, the montage.fif and fiducial.fif files.

[larsoner]

Are we good to go with these file formats?

Looks reasonable to me

These will supplement, not replace, the montage.fif and fiducial.fif

Okay but please do replace these files with the coord_frame set properly. It should be doable in MNE-Python master following https://github.com/mne-tools/mne-python/pull/8532

[ftadel]

Are we good to go with these file formats?

This looks good.

To read them in Brainstorm easily, we will need to merge them. If you import positions from a file that contains a list of labelled 3D points including NAS, LPA and RPA: these points are used to convert automatically all the other coordinates to the Brainstorm CTF/SCS coordinate system, and removed from the list.

Putting the other points with the electrodes would make your electrodes.tsv file not BIDS-compliant anymore (minus the file name). Which is probably not a problem, since the other file fiducials.tsv is not BIDS-compliant anyway... Since merging them is very easy, I'm not going to push you to do it if it feels more natural to you like this :)

[christian-oreilly]

@ftadel The only reason I did not merge these is to keep the electrodes.tsv file closer to BIDS-compliance (minus the file name as you mentioned). I looked at going all the way BIDS-compliant by adding the coordsystem.json file, but reading at https://bids-specification.readthedocs.io/en/stable/99-appendices/08-coordinate-systems.html it seems like this case would fall into the "individual" type of "nonstandard coordinate system identifiers" which requires "specifying an additional, participant-specific file to be fully defined"... so it feels like a rabbit hole and I am not sure there is real need for this for now. I'll be happy to revisit this if it becomes needed at some point.

[ftadel]

Indeed, a coordsystem.json might not be needed now. Any text file with labelled 3D points easily editable (hence mergeable) by humans is good. Your solution is perfectly fine, everybody wanting to use your atlases will be able to easily.

[christian-oreilly]

@larsoner Yep, the coordinate frame has been fixed. There is currently a bug MNE so it works only after locally correcting this bug in MNE code. I reported this bug in mne-tools/mne-python#8532 .

[christian-oreilly]

@ftadel In what coordinate system ('voxel', 'mri', 'scs', 'mni', 'world') are the sensors when we look at them in BS? scs?

[christian-oreilly]

@ftadel Actually, I think I will expend my question... here is my updated attempt at importing the channel nets:

    % Create folder for channel file
    for iMontage = 1:5
        switch iMontage
            case 1
                MontageName = '10-5';
                RefChannelFile = 'channel_ASA_10-05_343.mat';
            case 2
                MontageName = '10-10';
                RefChannelFile = 'channel_10-10_65.mat';
            case 3
                MontageName = '10-20';
                RefChannelFile = 'channel_10-20_19.mat';
            case 4
                MontageName = 'HGSN128';
                RefChannelFile = 'channel_GSN_HydroCel_128_E1.mat';
            case 5
                MontageName = 'HGSN128';
                RefChannelFile = 'channel_GSN_HydroCel_128_E001.mat';
        end
        % Load reference channel file
        ChannelMat = in_bst_channel(bst_fullfile(RefChannelPath, RefChannelFile));
        FolderName = file_standardize(str_remove_parenth(ChannelMat.Comment));

        % Create new folder
        iStudy = db_add_condition(subjectName, FolderName);
        sStudy = bst_get('Study', iStudy);

        montage = [MontageName '_electrodes.tsv'];

        % Read channel file
        channels = tdfread([FsDir '/montages/' montage], '\t');  
        channels_mri = [channels.x channels.y channels.z] ./ 1000;  

        for iChannel=1:size(channels_mri, 1)
            if iMontage >= 4
                assert(str2num(ChannelMat.Channel(iChannel).Name(2:end)) == ...
                       str2num(channels.name(iChannel, 2:end)));   
                ChannelMat.Channel(iChannel).Loc = transpose(channels_mri(iChannel, :));
            else                
                ChannelMat.Channel(iChannel) = ChannelMat.Channel(1);
                ChannelMat.Channel(iChannel).Loc = transpose(channels_mri(iChannel, :));
                ChannelMat.Channel(iChannel).Name = strtrim(channels.name(iChannel, :));
            end
        end
        ChannelMat.SCS.LPA = sFid.LPA;
        ChannelMat.SCS.NAS = sFid.NAS;
        ChannelMat.SCS.RPA = sFid.RPA;           

        %ChannelMat.HeadPoints = FifMat.HeadPoints;
        % Re-Register based on NAS/LPA/RPA
        % ChannelMat = channel_detect_type(ChannelMat, 1, 0);        

        % Save channel file
        ChannelFile = bst_fullfile(bst_fileparts(file_fullpath(sStudy.FileName)), RefChannelFile);
        bst_save(ChannelFile, ChannelMat, 'v7');
        db_reload_studies(iStudy);
    end    

Here is how the surfaces are imported:

    % Import head surface
    BemFiles = {...
        bst_fullfile(FsDir, 'bem', 'inner_skull.surf'), ...
        bst_fullfile(FsDir, 'bem', 'outer_skull.surf'), ...
        bst_fullfile(FsDir, 'bem', 'outer_skin.surf')};
    [iNewSurfaces, BstBemFiles] = import_surfaces(iSubject, BemFiles, 'FS', 0);    

And here is the current result:

The vertices and the channel position, as they are saved in the file that are here imported, are already co-registered (except for a 1000 multiplicative factor if I remember well). So, apparently, some transform is applied under the hood either to the surfaces or to the sensor position when they are imported. So, I guess the question is what is that transform that is applied on what (the surfaces or the sensors?) and how should I bring them back in co-registration?

[ftadel]

In what coordinate system ('voxel', 'mri', 'scs', 'mni', 'world') are the sensors when we look at them in BS? scs?

SCS, indeed: https://neuroimage.usc.edu/brainstorm/CoordinateSystems#Subject_Coordinate_System_.28SCS_.2F_CTF.29

[ftadel]

The surfaces are always going to be converted to SCS, based on the NAS/LPA/RPA points defined in the MRI imported before the surfaces.

For the electrodes positions, I think this is optional to convert to SCS, based on the NAS/LPA/RPA points defined in the positions file (.tsv?) imported. If you disable the call to channel_detect_type, it might keep all the 3D positions as they are in the file (not sure, maybe it does it anyway).

But since you MUST have everything in SCS in Brainstorm, you MUST convert everything correctly to SCS. And this is done independently for the electrodes and for the surfaces, you need to check both procedures. When you import the MRI (before importing the surfaces), the NAS/LPA/RPA must be set to points that are matching EXACTLY the positions of the NAS/LPA/RPA defined in the .tsv positions file.

[christian-oreilly]

the NAS/LPA/RPA must be set to points that are matching EXACTLY the positions of the NAS/LPA/RPA defined in the .tsv positions file.

Except for a factor 1000 according to the code in channel_detect_type.m

Actually, as I dive a bit deeper in BS code I find it a bit confusing because there is a lot of these *1000 and /1000 around. It seems clear to me that the MRI space is in mm. Can you confirm if the SCS space is in meter or in mm? Just to make sure that everything is coherent in my code...

[christian-oreilly]

Actually, I take that back, even for the MRI space it is not that clear because there is code like...

ChannelMat.Channel(i).Loc = cs_convert(ChannelMat, 'mri', 'scs', ChannelMat.Channel(i).Loc' ./ 1000)' .* 1000;

This transform doing MRI --> SCS its input (ChannelMat.Channel(i).Loc' ./ 1000) should be in proper units for the MRI space... which in this case is meter given the /1000... So the MRI space is in MM but the saved fiducials in structures like sMri are expected to be in meters?

[ftadel]

You're right, some things are not very obvious. Some values were kept in millimeters for backward compatibility... This is the case of the MRI coordinate system, as documented here: https://neuroimage.usc.edu/brainstorm/CoordinateSystems#MRI_coordinates

In the volumes, the fields MRI.SCS.NAS / .LPA / .RPA are in MRI coordinates (millimeters) In the channel files, fields ChannelMat.SCS.NAS / .LPA / .RPA are in SCS coordinates (meters)

Using the same cs_convert for both requires these *1000 and /1000 hacks, as documented here: https://github.com/brainstorm-tools/brainstorm3/blob/master/toolbox/sensors/channel_detect_type.m#L192

[christian-oreilly]

Maybe I am a bit thick here but if ChannelMat.SCS.* are in meters (e.g., == 1 [m]), then ChannelMat.SCS.NAS ./ 1000 (i.e., 0.001 [km]) are in km and cs_convert expect to work with mm ("cs_convert() is intended to work on sMri structures, in which NAS/LPA/RPA/T fields are in millimeters"). Am I crazy or this operation is inverted? Or it is my internal semantic parser that fails me badly? :)

[ftadel]

In cs_convert, the T field (translation) is divided by 1000, but not the coordinates of the points, as it is designed for the MRI coordinates in millimeters. If you want the transformation to be correct for input coordinates in meters, you can either multiple the .T field by 1000, or divide the input coordinates by 1000. And the second is easier to do as a one-liner. Does it make sense?

[ftadel]

A Brainstorm user reminded us that this still requires some work :) https://neuroimage.usc.edu/forums/t/one-click-in-bst/27268/4

[christian-oreilly]

Solved by #9 thanks to the work of @ftadel!