label_merger.py

import mne
import os
import numpy as np
import os.path as op
from mne.datasets import sample
from nose.tools import assert_equal


data_path = sample.data_path()

meg_path = data_path + '/MEG/sample'

subject = 'sample'
subjects_dir = data_path + '/subjects'

label_fname  = meg_path + '/labels/Aud-lh.label'
label = mne.read_label(label_fname, subject=subject)



def get_vertices_faces(vertices_ind, subjects_dir, subject, hemi, surf):

    surf = mne.surface.read_surface(op.join(subjects_dir, subject, 'surf', hemi + '.' + surf))
    ind = np.in1d(surf[1][:,0], vertices_ind)*\
        np.in1d(surf[1][:,1], vertices_ind)*\
        np.in1d(surf[1][:,2], vertices_ind)

    vertices = surf[0][vertices_ind]
    faces = surf[1][ind]
    faces = np.array([np.where(vertices_ind == face)[0]
                      for face in faces.ravel()]).reshape(faces.shape)

    return vertices, faces


def triangle_area(vertices, faces):

    r12 = vertices[faces[:,0],:]
    r13 = vertices[faces[:,2],:] - r12
    r12 = vertices[faces[:,1],:] - r12
    return np.sum(np.sqrt(np.sum(np.cross(r12, r13)**2,axis=1))/2.)



def compute_area(label_fname, subject):
    log_file = op.join(subjects_dir, subject, 'label', 'temp.log')
    command = 'mris_anatomical_stats -l ' + label_fname + ' ' + subject + ' ' + label_fname[-8:-6] + ' white >& ' + log_file
    status = os.system(command)

    if not bool(status):
        datafile = file(log_file)
        for line in datafile:
            if 'total surface area' in line:
                datafile.close()
                return int(line.split(' ')[-2])
        os.remove(log_file)


vertices, faces = get_vertices_faces(label.vertices, subjects_dir, label.subject, label.hemi, 'white')

label_area_python = triangle_area(vertices, faces)

label_area_mne = compute_area(label_fname, subject)


def test_triangle_area():
    faces = np.array([[0, 1, 2], [0, 2, 3]])
    vertices = np.array([[0, 0, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0]])
    assert_equal(triangle_area(vertices, faces), 1)