real-time source modeling

ftclient_rt_source_modelling.py

#!/usr/bin/env python
"""
=============================================
Visualize real-time source estimates of
the most recent samples from FieldTrip client
=============================================

Please refer to `ftclient_rt_average.py` for instructions on
how to get the FieldTrip connector working in MNE-Python.

This example demonstrates how to use it for continuous
visualization of power spectra in real-time using buffer_as_epoch method.

"""

# Author: Sebastian Silfverberg <sebastian.silfveberg@aalto.fi>
#
# License: BSD (3-clause)


print(__doc__)

import os.path as op
from tempfile import mkdtemp
import numpy as np
from mayavi import mlab

from mne.realtime import FieldTripClient
import mne
from mne.minimum_norm.inverse import prepare_inverse_operator
#from mne.datasets import sample
from mne.minimum_norm import read_inverse_operator, compute_source_psd_epochs
from surfer import Brain

import datetime as dt

savedir = mkdtemp()

data_path = '/home/mainak/Desktop/parkkonen_lauri'
#subjects_dir = data_path + '/subjects'
subjects_dir = '/home/mainak/Desktop/'
fname_inv = data_path + '/lppassivel01raw_tsss_mc-5-meg-inv.fif'
subject_id = 'parkkonen_lauri'

snr = 3.0
lambda2 = 1.0 / snr ** 2
method = "dSPM"  # use dSPM method (could also be MNE or sLORETA)

# user must provide list of bad channels because
# FieldTrip header object does not provide that

bads = []

# Load data
inverse_operator = read_inverse_operator(fname_inv)
inv = prepare_inverse_operator(inverse_operator, nave=1, lambda2=lambda2,
                               method=method)

with FieldTripClient(host='sinuhe', port=1972,
                     tmax=150, wait_max=10) as rt_client:

    # get measurement info guessed by MNE-Python
    raw_info = rt_client.get_measurement_info()

    # pick MEG channels
    picks = mne.pick_types(raw_info, meg=True, eeg=False, stim=False,
                           eog=False, exclude=bads)

    # plot the initial brain surface
    surface = 'inflated'
    hemi = 'both'

    # Views parameter can be
    # | 'lateral' | 'medial'  | 'rostral' | 'caudal'   |
    # | 'dorsal'  | 'ventral' | 'frontal' | 'parietal' |

    brain = Brain(subject_id=subject_id, hemi=hemi, surf=surface,
                  subjects_dir=subjects_dir, views='dorsal', title='')

    # define frequencies of interest

    fmin, fmax = 8., 11.
    bandwidth = 2.  # bandwidth of the windows in Hz

    n_samples = 512
    n_fft = 512  # number of FFT. Preferably a power of two.
    time_label = '{0:1d}th window of {1:d} latest samples'

    for ii in range(100):

        epoch = rt_client.get_data_as_epoch(n_samples=n_samples, picks=picks)

        tmin = epoch.events[0][0] / raw_info['sfreq']
        tmax = (epoch.events[0][0] + n_samples) / raw_info['sfreq']

        print("Just got buffer %d" % (ii + 1))
        print('%0.2f to %0.2f secs.' % (tmin, tmax))

        # compute source space psd in label

        # Note: By using "return_generator=True" stcs will be a generator
        # object instead of a list. This allows us so to iterate without
        # having to keep everything in memory.

        t1 = dt.datetime.now()

        stcs = compute_source_psd_epochs(epoch, inv,
                                         lambda2=lambda2, method=method,
                                         fmin=fmin, fmax=fmax,
                                         bandwidth=bandwidth,
                                         return_generator=True,
                                         nave=1, prepared=True)

        #print('compute_source_psd_epochs = %f' %
        #      ((dt.datetime.now() - t1).total_seconds()))
        #print((dt.datetime.now() - t1).total_seconds())
        # plotting the source estimates with PySurfer

        for i, stc in enumerate(stcs):

            print('compute_source_psd_epochs = %f' %
                  ((dt.datetime.now() - t1).total_seconds()))
            #print(i)
            # for vertices and array use lh_vertno for hemi='lh'
            # and rh_vertno for hemi='rh'

            vertices = stc.lh_vertno
            array = np.mean(stc.lh_data, axis=1)

            #t1 = dt.datetime.now()

            brain.add_data(array=array, min=1000, max=50000, thresh=None,
                           colormap='hot', vertices=vertices,
                           smoothing_steps=3,
                           time_label=time_label.format(ii + 1, n_fft),
                           hemi='lh', alpha=0.7,
                           remove_existing=True)

            #print('add_data = %f' %
            #      ((dt.datetime.now() - t1).total_seconds()))

            vertices = stc.rh_vertno
            array = np.mean(stc.rh_data, axis=1)
            brain.add_data(array=array, min=1000, max=50000, thresh=None,
                           colormap='hot', vertices=vertices,
                           smoothing_steps=3,
                           time_label=time_label.format(ii + 1, n_fft),
                           hemi='rh', alpha=0.7,
                           remove_existing=True)
            mlab.view(-90, 90)
            mlab.title(text='%d to %d Hz' % (fmin, fmax), size=0.3, height=0.9)
            mlab.text(x=0.1, y=0.1, text='%0.2f to %0.2f secs.' % (tmin, tmax),
                      width=0.25)
        # XXX: force the figure to be shown

        fig = mlab.gcf()
        fig.scene.save_bmp(op.join(savedir, 'test.bmp'))
mlab.close()