tSSS on Empty-Room Recording

lkandemir · March 16, 2021, 10:43am

Hi,

In my preprocessing pipeline, I apply tSSS to all of my recordings. I have a question regarding how this is applied on empty-room recordings. In both MaxFilter and MNE implementations of tSSS, the expansion origin of the spherical harmonics is (0, 0, 0) in MEG coordinate system when there is no head-device transformation information. Since this information is omitted in empty-room recordings, it is automatically taken care of. However, in Brainstorm, a default head-device transformation is applied automatically to empty-room recordings, so it is possible to still use head coordinate system (which is the default in Brainstorm). I personally tried both head and meg coordinate systems for my empty-room recordings and I somehow get much better results (i.e. results as expected) when I use head coordinate system. When I try to use meg coordinate system, the variation between the channels becomes so high that it seems tSSS is introducing new noise to the data. I believe the reason for this could be that when automatic head transformation is applied, meg coordinate system is affected too and maybe the new origin is outside the MEG helmet. Is there a way to deactivate this automatic transformation? Also, would you have any further suggestions on this issue?

Sorry for the long post and thank you for your help.

Best regards,
Levent

Francois · March 17, 2021, 10:45am

Indeed, there is a problem here.

Brainstorm adds a default device=>head transformation when it is missing from the FIF file (Z +4cm), in order to homogenize all the default display procedures. Unlike in most other MEG software, the Brainstorm 2D topo plots are projections of the 3D sensors on a 2D plane, therefore we always need valid 3D coordinates (even if fake) for rendering the sensor recordings or power in 2D. Hence this default transformation.
https://github.com/brainstorm-tools/brainstorm3/blob/master/toolbox/io/in_channel_fif.m#L73

When exporting the the contents of the channel file to a MNE-Python structure, all the transformations added by Brainstorm to the sensors locations are reverted to what was present in the initial FIF file. However, when exporting, the default Brainstorm device=>head transformation is exported as an extra dev_head_t field that was initially not present in the FIF file. And its presence seems to mess up with the behavior of the mne.maxwell_filter function (and probably others).

If I understood well the problem, this issue could be solved by keeping everything unchanged in the data import to Brainstorm but preventing from exporting the default Z+4cm transformation as a dev_head_t field in the MNE-Python object.

This is what I tried to do in this commit:
https://github.com/brainstorm-tools/brainstorm3/commit/75a6794ac0d2f648a5b9f9667d634b49fd22a575#diff-12c42dea9a032ee334e1e52881bfd66e3f7329c715c06021ed22c2893729570fR134

Please review the changes to see if it makes sense.
Then update Brainstorm to get this fix and try again filtering your noise recordings.

Good catch - Thanks!

Francois · March 17, 2021, 10:49am

@Alexandre @dengeman FYI

lkandemir · March 19, 2021, 12:38pm

Hi Francois,

Thank you for swift response. I checked the results after the bug-fix and it seems that the results are not better when 'meg' coordinate frame is selected. Moreover, when 'head' coordinate frame is selected, it is worse compared to before the bug-fix. In order to cross-check, I ran the same pipeline in Python. It seems that 'head' coordinate frame was producing similar output to MNE-Python before the bug-fix. However, it is not the case anymore. Also, when 'head' coordinate frame is selected in MNE-Python, MNE throws an error stating that there is no device-head transformation so 'meg' coordinate frame should be selected. However, this is not the case in BST. It does not complain even if I select 'head' coordinate frame. Therefore, I believe somehow the transformation information still remains after the bug-fix. I attached the results both in MNE-Python and Brainstorm to this post for evaluation.

Thank you,
Best regards,
Levent

[tSSS_BST.pdf|attachment](uploatSSS_BST.pdf (1.0 MB) tSSS_MNE.pdf (722.5 KB) d://sBq4Rzeb33g9dQlii8BQF1sif7A.pdf) (1.0 MB)

Francois · March 23, 2021, 8:44am

I was hoping I could find some quick solution, but I didn't, this will require some finer debugging, starting with a comparison of the maxwell_filter input of:

the Python object generated by Brainstorm (out_mne_data.m/out_mne_channel.m),
the Python object read from the MNE-Python script.

I might be able to work on this this week, but it would help me a lot if you could provide me with some material to get started with this comparison. Could you please:

Share the example empty-room recordings you did your tests with (the original .fif file)? Upload it somewhere and post the download link here.
Include any additional file needed for processing it (calibration / cross-talk files):
Copy-paste here your MNE-Python script to load and process this file with tSSS.
Include a screen capture of the options window for the process in Brainstorm (or the corresponding .m generated by Brainstorm)

Thanks!

lkandemir · March 23, 2021, 12:44pm

Hi Francois,

Thank you for looking into this! I uploaded the files to the link below. I also uploaded the notebook file of the tSSS application in MNE-Python including Welch's PSD comparison before and after tSSS.

https://drive.google.com/drive/folders/1Ft3Z2earYHnqWxfX1wIQj0Kk4gnHjibV?usp=sharing

The codes I used in Python are below:

data_path = 'F:/emptyroom_05062018.fif'
raw = mne.io.read_raw_fif(data_path)
raw.info['bads'] = ['MEG0422', 'MEG0643', 'MEG0811', 'MEG1032']
fine_cal_file = 'F:/sss_cal_3042_131004_422.dat'
crosstalk_file = 'F:/ct_sparse.fif'

nraw=mne.preprocessing.maxwell_filter(raw, origin='auto', coord_frame='meg', int_order=8, ext_order=3, 
                            calibration=fine_cal_file, cross_talk=crosstalk_file, st_duration=60, 
                            st_correlation=0.98, regularize='in', ignore_ref=False, st_fixed=True,
                            st_only=False, mag_scale=100.0, verbose=True,
                            skip_by_annotation=('edge', 'bad_acq_skip'))

And the settings in Brainstorm are:

Please let me know if you need anything else.

Best regards,
Levent

Francois · March 25, 2021, 4:23pm

Thanks for the example files.
I fixed one very obvious bug in the channel file export to MNE-Python:
https://github.com/brainstorm-tools/brainstorm3/commit/68d19bcda65a87fe2f036ad192827a448f177b9f

Can you please update Brainstorm and try again?

I tested I was obtaining the same signals from the MNE-Python FIF reader or the Brainstorm import+export as MNE RawArray object with the following script:

% Read Brainstorm file
MatFile = 'C:\Work\Protocols\Test\data\MNE\@rawemptyroom_05062018\data_0raw_emptyroom_05062018.mat';
pyBst = out_mne_data(MatFile, 'Raw');
pyBst.info{'bads'} = py.list({'MEG0422', 'MEG0643', 'MEG0811', 'MEG1032'});

% Read Python file
noise_file = 'C:\Work\RawData\Test\MNE_tSSS\emptyroom_05062018.fif';
pyMne = py.mne.io.read_raw_fif(noise_file);
pyMne.info{'bads'} = py.list({'MEG0422', 'MEG0643', 'MEG0811', 'MEG1032'});

% tSSS options
fine_cal_file = 'C:\Work\RawData\Test\MNE_tSSS\sss_cal_3042_131004_422.dat';
crosstalk_file = 'C:\Work\RawData\Test\MNE_tSSS\ct_sparse.fif';
sss_args = pyargs(...
    'origin', 'auto', ...
    'coord_frame', 'meg', ...
    'int_order', py.int(8), ...
    'ext_order', py.int(3), ...
    'calibration', fine_cal_file, ...
    'cross_talk', crosstalk_file, ...
    'st_duration', 60, ...
    'st_correlation', 0.98, ...
    'regularize', 'in', ...
    'ignore_ref', py.False, ...
    'st_fixed', py.True, ...
    'st_only', py.False, ...
    'mag_scale', 100.0, ...
    'verbose', py.True, ...
    'skip_by_annotation', py.tuple({'edge', 'bad_acq_skip'}));

% Run tSSS
pyBst_sss = py.mne.preprocessing.maxwell_filter(pyBst, sss_args);
pyMne_sss = py.mne.preprocessing.maxwell_filter(pyMne, sss_args);

% Compare visually
pyBst_sss.plot(pyargs('duration',0.5));
pyMne_sss.plot(pyargs('duration',0.5));

lkandemir · March 27, 2021, 10:04pm

Thank you for the update. I am getting the below error now:

bst_error

Could this be related to the Matlab version? I am using 2019b currently.

Francois · March 29, 2021, 12:53pm

The Matlab version should not be a problem. I tested with 2019a and 2018a and this syntax seems to work.
But there is maybe something wrong with the way the Python environment is initialized.

Can you try executing the following commands, in this order, and copy-paste all the output from the command window?

py.NoneType
brainstorm
[isOk, errorMsg, pyVer] = bst_python_init('Initialize', 1)
py.NoneType
pyversion
py.NoneType
PythonExeBst = bst_get('PythonExe')
[pyVer, PythonExeMatlab, isLoaded] = bst_python_ver()

Thanks

lkandemir · March 30, 2021, 7:42am

I got the below output with those commands:

py.NoneType

Unable to resolve the name py.NoneType.

brainstorm

BST> Starting Brainstorm:
BST> =================================
BST> Version: 25-Mar-2021
BST> Compiling main interface files...
BST> Emptying temporary directory...
BST> Deleting old process reports...
BST> Loading configuration file...
BST> Checking internet connectivity... BST> ERROR: webread and urlread failed reading URL: http://neuroimage.usc.edu/bst/getversion.php
failed
BST> Could not check for Brainstorm updates.
BST> Initializing user interface...
BST> Starting OpenGL engine... hardware
BST> Reading process folder...
BST> Loading current protocol...
BST> =================================

[isOk, errorMsg, pyVer]=bst_python_init('Initialize', 1)

BST> Python 3.8 executable: C:\Users\Administrator\AppData\Local\Programs\Python\Python38\python.exe

isOk =

1

errorMsg =

pyVer =

'3.8'

py.NoneType

Unable to resolve the name py.NoneType.

pyversion

   version: '3.8'
executable: 'C:\Users\Administrator\AppData\Local\Programs\Python\Python38\python.exe'
   library: 'C:\Users\Administrator\AppData\Local\Programs\Python\Python38\python38.dll'
      home: 'C:\Users\Administrator\AppData\Local\Programs\Python\Python38'
  isloaded: 1

py.NoneType

Unable to resolve the name py.NoneType.

PythonExeBst = bst_get('PythonExe')

PythonExeBst =

'C:\Users\Administrator\AppData\Local\Programs\Python\Python38\python.exe'

[pyVer, PythonExeMatlab, isLoaded] = bst_python_ver()

pyVer =

'3.8'

PythonExeMatlab =
'C:\Users\Administrator\AppData\Local\Programs\Python\Python38\python.exe'

Thank you for your help.

Levent

Francois · March 30, 2021, 8:21am

BST> Checking internet connectivity... BST> ERROR: webread and urlread failed reading URL: http://neuroimage.usc.edu/bst/getversion.php

Are you behind a firewall that is not configured in Matlab?

Unable to resolve the name py.NoneType.

py.None seems to be available with all Matlab versions, I fixed the code to use it:

Please update Brainstorm again and tell me how it goes.

lkandemir · March 30, 2021, 10:41am

Hi Francois,

Are you behind a firewall that is not configured in Matlab?

Actually the computer I am using is offline. Therefore, it can't check for the newer version.

I updated and it seems to be working now. I also checked with 'head' selected as the coordinate frame and it gives the error just as MNE-Python and suggests me to use 'meg'. On the way back from Python, BST tells me that it applied the default transformation to the file but I think this is a normal behaviour since the file does not contain any transformation information and BST needs that for visualization. I cross checked the MNE and BST results for tSSS and it seems almost identical. There are only really minor (one has to look really careful to spot them) differences at low frequencies (in PSD). I think this could be related to rounding errors during data transformation and probably normal, right?

Thank you once again for all the help.

Best regards,
Levent

Francois · March 31, 2021, 7:37am

This is great news!
Have you checked by computing and displaying the PSD in Brainstorm for both files?

Some minor differences are introduced with the successive changes of coordinates systems in Brainstorm. The positions of all the sensors and digitized head points are first transformed to the Brainstorm SCS coordinates, and then this transformation is inverted and applied again to revert to the original Neuromag device coordinates. These two transformations cause small differences in the sensors coordinates (up to 0.1 mm). I'd consider these differences to be minor enough not to be worrying for MEG head localization, but they are probably sufficient for causing observable differences it the SSS filtering.

lkandemir · March 31, 2021, 9:00am

Have you checked by computing and displaying the PSD in Brainstorm for both files?

Yes, that is how I checked the differences and as I said they are really minor. Great, thank you for the help and information.

Best regards,
Levent