glass_brain.py

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Use glass brain plotting for surface STCs.
"""

import os.path as op
import numpy as np
from scipy import sparse
import mne
import nibabel as nib
from nilearn.plotting import plot_glass_brain

data_path = mne.datasets.sample.data_path()
subjects_dir = data_path + '/subjects'
subject = 'sample'
stc = mne.read_source_estimate(
    data_path + '/MEG/sample/sample_audvis-meg-eeg-lh.stc')
surf_to_mri = 0.
n_vertices = sum(len(v) for v in stc.vertices)
offset = 0
for hi, hemi in enumerate(('lh', 'rh')):
    ribbon = nib.load(op.join(subjects_dir, subject, 'mri',
                              '%s.ribbon.mgz' % hemi))
    xfm = ribbon.header.get_vox2ras_tkr()
    mri_data = ribbon.get_data()
    ijk = np.array(np.where(mri_data)).T
    xyz = mne.transforms.apply_trans(xfm, ijk) / 1000.
    row_ind = np.where(mri_data.ravel())[0]
    data = np.ones(row_ind.size)
    rr = mne.read_surface(subjects_dir + '/sample/surf/%s.white' % hemi)[0]
    rr /= 1000.
    rr = rr[stc.vertices[hi]]
    col_ind = mne.surface._compute_nearest(rr, xyz) + offset
    surf_to_mri = surf_to_mri + sparse.csr_matrix(
        (data, (row_ind, col_ind)), shape=(mri_data.size, n_vertices))
    offset += len(stc.vertices[hi])

stc.crop(0.08, 0.08)
data = surf_to_mri.dot(stc.data[:, 0])
data.shape = mri_data.shape

###############################################################################
# This is taken from mne.source_space._read_talxfm:
fname = op.join(subjects_dir, subject, 'mri', 'transforms',
                'talairach.xfm')
# read the RAS to MNI transform from talairach.xfm
with open(fname, 'r') as fid:

    # read lines until we get the string 'Linear_Transform', which precedes
    # the data transformation matrix
    got_it = False
    comp = 'Linear_Transform'
    for line in fid:
        if line[:len(comp)] == comp:
            # we have the right line, so don't read any more
            got_it = True
            break

    if got_it:
        xfm = list()
        # read the transformation matrix (3x4)
        for ii, line in enumerate(fid):
            digs = [float(s) for s in line.strip('\n;').split()]
            xfm.append(digs)
            if ii == 2:
                break
        xfm.append([0., 0., 0., 1.])
        xfm = np.array(xfm, dtype=float)
    else:
        raise ValueError('failed to find \'Linear_Transform\' string in '
                         'xfm file:\n%s' % fname)
###############################################################################

fmin, fmax = 5., 15.
affine = np.dot(xfm, ribbon.header.get_vox2ras())
img = nib.Nifti1Image(data, affine)
plot_glass_brain(img, threshold=fmin, vmax=fmax)

brain = stc.plot(subject='sample', subjects_dir=subjects_dir, hemi='rh',
                 clim=dict(kind='value',
                           pos_lims=[fmin, (fmin + fmax) / 2., fmax]))