Python mne.extract_label_time_course函数代码示例

def cal_labelts(stcs_path, fn_func_list, condition='LLst',
                min_subject='fsaverage', subjects_dir=None):
    '''
    Extract stcs from special ROIs, and store them for funther causality
    analysis.

    Parameter
    ---------
    stcs_path: string
        The path of stc's epochs.
    fn_ana_list: string
        The path of the file including pathes of functional labels.
    condition: string
        The condition for experiments.
    min_subject: the subject for common brain
    '''
    path_list = get_files_from_list(stcs_path)
    minpath = subjects_dir + '/%s' % (min_subject)
    srcpath = minpath + '/bem/fsaverage-ico-5-src.fif'
    src_inv = mne.read_source_spaces(srcpath)
    # loop across all filenames
    for stcs_path in path_list:
        caupath = stcs_path[:stcs_path.rfind('/%s' % condition)]
        fn_stcs_labels = caupath + '/%s_labels_ts.npy' % (condition)
        _, _, files = os.walk(stcs_path).next()
        trials = len(files) / 2
        # Get unfiltered and morphed stcs
        stcs = []
        i = 0
        while i < trials:
            fn_stc = stcs_path + 'trial%s_fsaverage' % (str(i))
            stc = mne.read_source_estimate(fn_stc + '-lh.stc',
                                           subject=min_subject)
            stcs.append(stc)
            i = i + 1
        # Get common labels
        list_file = fn_func_list
        with open(list_file, 'r') as fl:
                  file_list = [line.rstrip('\n') for line in fl]
        fl.close()
        rois = []
        labels = []
        for f in file_list:
            label = mne.read_label(f)
            labels.append(label)
            rois.append(label.name)
        # Extract stcs in common labels
        label_ts = mne.extract_label_time_course(stcs, labels, src_inv,
                                                 mode='pca_flip')
        # make label_ts's shape as (sources, samples, trials)
        label_ts = np.asarray(label_ts).transpose(1, 2, 0)
        np.save(fn_stcs_labels, label_ts)

    picks = mne.fiff.pick_channels_regexp(raw.info["ch_names"], "M..-*")
    raw.filter(l_freq=1, h_freq=50, picks=picks)
    er_raw.filter(l_freq=1, h_freq=50, picks=picks)

    noise_cov = mne.compute_raw_data_covariance(er_raw)
    # note that MNE reads CTF data as magnetometers!
    noise_cov = mne.cov.regularize(noise_cov, raw.info, mag=noise_reg)
    inverse_operator = mne.minimum_norm.make_inverse_operator(raw.info, forward, noise_cov, loose=0.2, depth=0.8)
    data, time = raw[0, :]  #
    events = fg.get_good_events(markers[subj], time, window_length)

    epochs = mne.Epochs(raw, events, None, 0, window_length, preload=True, baseline=None, detrend=0, picks=picks)
    stcs = mne.minimum_norm.apply_inverse_epochs(epochs, inverse_operator, lambda2, "MNE", return_generator=False)

    labels, label_colors = mne.labels_from_parc(subj, parc="aparc")
    label_ts = mne.extract_label_time_course(stcs, labels, forward["src"], mode=label_mode)

    # label_data is nlabels by time, so here we can use whatever connectivity method we fancy
    con, freqs, times, n_epochs, n_tapers = mne.connectivity.spectral_connectivity(
        label_ts,
        method=method,
        mode="multitaper",
        sfreq=raw.info["sfreq"],
        fmin=[1, 4, 8, 13, 30],
        fmax=[4, 8, 13, 30, 50],
        faverage=True,
        n_jobs=3,
        mt_adaptive=False,
    )
    np.save(dir_out + subj + "-" + label_mode + "-" + "-".join(method), con)

snr = 1.
lambda2 = 1. / snr**2

labels = mne.read_labels_from_annot(
    subject=subject, parc="PALS_B12_Brodmann", regexp="Brodmann")

condition = "interupt"

inv = read_inverse_operator(mne_folder + "%s_%s-inv.fif" % (subject, condition
                                                            ))
epochs = mne.read_epochs(epochs_folder + "%s_%s-epo.fif" % (subject, condition
                                                            ))
# epochs.resample(500)

stcs = apply_inverse_epochs(
    epochs["press"], inv, lambda2, method=method, pick_ori=None)
ts = [
    mne.extract_label_time_course(
        stc, labels, inv["src"], mode="mean_flip") for stc in stcs
]

# for h, tc in enumerate(ts):
#     for j, t in enumerate(tc):
#         t *= np.sign(t[np.argmax(np.abs(t))])
#         tc[j, :] = t
#     ts[h] = tc

ts = np.asarray(ts)
stc.save(source_folder + "%s_%s_epo" % (subject, condition))
np.save(source_folder + "ave_ts/%s_%s_ts-epo.npy" % (subject, condition), ts)

#.........这里部分代码省略.........
        else:
            epochs = mne.Epochs(raw, events, event_id, t_min, t_max,
                                picks=picks, baseline=(None, 0), reject=reject)
            stc = apply_inverse_epochs(epochs, inverse_operator, lambda2,
                                       inv_method, pick_ori=None)

            print '***'
            print 'len stc %d' % len(stc)
            print '***'

    elif is_epoched and event_id is None:
        stc = apply_inverse_epochs(epochs, inverse_operator, lambda2,
                                   inv_method, pick_ori=None)

        print '***'
        print 'len stc %d' % len(stc)
        print '***'
    else:
        stc = apply_inverse_raw(raw, inverse_operator, lambda2, inv_method,
                                label=None,
                                start=None, stop=None,
                                buffer_size=1000,
                                pick_ori=None)  # None 'normal'

        print '***'
        print 'stc dim ' + str(stc.shape)
        print '***'

    if save_stc:
        if aseg:
            for i in range(len(stc)):
                try:
                    os.mkdir(op.join(subj_path, 'TS'))
                except OSError:
                    pass
                stc_file = op.join(subj_path, 'TS', basename + '_' +
                                   inv_method + '_stc_' + str(i) + '.npy')

                if not op.isfile(stc_file):
                    np.save(stc_file, stc[i].data)

    labels_cortex = mne.read_labels_from_annot(sbj_id, parc=parc,
                                               subjects_dir=sbj_dir)
    if is_blind:
        for l in labels_cortex:
            if l.name in labels_removed:
                print l.name
                labels_cortex.remove(l)

    print '\n*** %d ***\n' % len(labels_cortex)

    src = inverse_operator['src']

    # allow_empty : bool -> Instead of emitting an error, return all-zero time
    # courses for labels that do not have any vertices in the source estimate
    label_ts = mne.extract_label_time_course(stc, labels_cortex, src,
                                             mode='mean',
                                             allow_empty=True,
                                             return_generator=False)

    # save results in .npy file that will be the input for spectral node
    print '\n*** SAVE ROI TS ***\n'
    print len(label_ts)

    ts_file = op.abspath(basename + '_ROI_ts.npy')
    np.save(ts_file, label_ts)

    if aseg:
        print sbj_id
        labels_aseg = get_volume_labels_from_src(src, sbj_id, sbj_dir)
        labels = labels_cortex + labels_aseg
    else:
        labels = labels_cortex

    print labels[0].pos
    print len(labels)

    labels_file = op.abspath('labels.dat')
    with open(labels_file, "wb") as f:
        pickle.dump(len(labels), f)
        for value in labels:
            pickle.dump(value, f)

    label_names_file = op.abspath('label_names.txt')
    label_coords_file = op.abspath('label_coords.txt')

    label_names = []
    label_coords = []

    for value in labels:
        label_names.append(value.name)
#        label_coords.append(value.pos[0])
        label_coords.append(np.mean(value.pos, axis=0))

    np.savetxt(label_names_file, np.array(label_names, dtype=str),
               fmt="%s")
    np.savetxt(label_coords_file, np.array(label_coords, dtype=float),
               fmt="%f %f %f")

    return ts_file, labels_file, label_names_file, label_coords_file

def test_extract_label_time_course():
    """Test extraction of label time courses from stc
    """
    n_stcs = 3
    n_times = 50

    src = read_inverse_operator(fname_inv)['src']
    vertices = [src[0]['vertno'], src[1]['vertno']]
    n_verts = len(vertices[0]) + len(vertices[1])

    # get some labels
    labels_lh, _ = labels_from_parc('sample', hemi='lh',
                                    subjects_dir=subjects_dir)
    labels_rh, _ = labels_from_parc('sample', hemi='rh',
                                    subjects_dir=subjects_dir)
    labels = list()
    labels.extend(labels_lh[:5])
    labels.extend(labels_rh[:4])

    n_labels = len(labels)

    label_means = np.arange(n_labels)[:, None] * np.ones((n_labels, n_times))

    # compute the mean with sign flip
    label_means_flipped = np.zeros_like(label_means)
    for i, label in enumerate(labels):
        label_means_flipped[i] = i * np.mean(label_sign_flip(label, src))

    # generate some stc's with known data
    stcs = list()
    for i in range(n_stcs):
        data = np.zeros((n_verts, n_times))
        # set the value of the stc within each label
        for j, label in enumerate(labels):
            if label.hemi == 'lh':
                idx = np.intersect1d(vertices[0], label.vertices)
                idx = np.searchsorted(vertices[0], idx)
            elif label.hemi == 'rh':
                idx = np.intersect1d(vertices[1], label.vertices)
                idx = len(vertices[0]) + np.searchsorted(vertices[1], idx)
            data[idx] = label_means[j]

        this_stc = SourceEstimate(data, vertices, 0, 1)
        stcs.append(this_stc)

    # test some invalid inputs
    assert_raises(ValueError, extract_label_time_course, stcs, labels,
                  src, mode='notamode')

    # have an empty label
    empty_label = labels[0].copy()
    empty_label.vertices += 1000000
    assert_raises(ValueError, extract_label_time_course, stcs, empty_label,
                  src, mode='mean')

    # but this works:
    tc = extract_label_time_course(stcs, empty_label, src, mode='mean',
                                   allow_empty=True)
    for arr in tc:
        assert_true(arr.shape == (1, n_times))
        assert_array_equal(arr, np.zeros((1, n_times)))

    # test the different modes
    modes = ['mean', 'mean_flip', 'pca_flip']

    for mode in modes:
        label_tc = extract_label_time_course(stcs, labels, src, mode=mode)
        label_tc_method = [stc.extract_label_time_course(labels, src,
                           mode=mode) for stc in stcs]
        assert_true(len(label_tc) == n_stcs)
        assert_true(len(label_tc_method) == n_stcs)
        for tc1, tc2 in zip(label_tc, label_tc_method):
            assert_true(tc1.shape == (n_labels, n_times))
            assert_true(tc2.shape == (n_labels, n_times))
            assert_true(np.allclose(tc1, tc2, rtol=1e-8, atol=1e-16))
            if mode == 'mean':
                assert_array_almost_equal(tc1, label_means)
            if mode == 'mean_flip':
                assert_array_almost_equal(tc1, label_means_flipped)

                             "%s_trial_start-epo.fif" % subject)
    # epochs.drop_bad_epochs(reject_params)
    # epochs.resample(250, n_jobs=4)

    for condition in conditions:
        stcs = apply_inverse_epochs(epochs[condition],
                                    inverse_operator,
                                    lambda2,
                                    method,
                                    pick_ori=None)

        for label in labels_sel:
            label_ts = []
            for j in range(len(stcs)):
                ts = mne.extract_label_time_course(stcs[j],
                                                   labels=label,
                                                   src=src,
                                                   mode="pca_flip")
                ts = np.squeeze(ts)
                ts *= np.sign(ts[np.argmax(np.abs(ts))])
                label_ts.append(ts)

            label_ts = np.asarray(label_ts)
            tfr = cwt_morlet(label_ts, epochs.info["sfreq"], freqs,
                             use_fft=True, n_cycles=n_cycle)

            np.save(tf_folder + "%s_%s_%s_%s_%s_sf-tfr" % (subject,
                                                           condition[:3],
                                                           condition[4:],
                                                           label.name, method),
                    tfr)
            np.save(tf_folder + "%s_%s_%s_%s_%s_sf-ts" % (subject,

        label = mne.read_label(label_dir + roi)
        roi_pretty = roi.split('.')[0]
    else:
        roi_pretty = roi.split('/')[-1].split('+')[0]
	# right labels are normally in the second index
	if avg_label[0] is not None:
        	label = avg_label[0].morph(subject_to=s)
	else:
        	label = avg_label[1].morph(subject_to=s)

    evoked_fname = '/Volumes/Shaw/MEG_data/analysis/stop/evoked/%s_stop_BP1-35_DS120-ave.fif' % s
    inv_fname = '/Volumes/Shaw/MEG_data/analysis/stop/%s_task-5-meg-inv.fif' % s
    inverse_operator = read_inverse_operator(inv_fname)

    for i, c in enumerate(conds):
        # calculate source estimates for the whole brain
        evoked = mne.read_evokeds(evoked_fname, condition=c)
        stc = apply_inverse(evoked, inverse_operator, lambda2, method,
                            pick_ori=None)
        ts = mne.extract_label_time_course(stc, label, inverse_operator['src'])
        data[i].append(ts)

# export one CSV file for each condition
for i, c in enumerate(conds):
    fname = out_dir + '%s_%s_%s.csv' % (c, roi_pretty, method)
    fid = open(fname, 'w')
    fid.write('time,' + ','.join(['%.3f' % t for t in evoked.times]) + '\n')
    for j, d in enumerate(data[i]):
        fid.write('%s,' % subjs[j] + ','.join(['%e' % t for t in d[0]]) + '\n')
    fid.close()

lambda2 = 1.0 / snr ** 2

# Compute inverse operator
inverse_operator = make_inverse_operator(evoked.info, fwd, noise_cov,
                                         depth=None, fixed=False)

stc = apply_inverse(evoked, inverse_operator, lambda2, inv_method,
                    pick_ori=None)

# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi
labels_parc = mne.read_labels_from_annot(
    subject, parc=parc, subjects_dir=subjects_dir)

###############################################################################
# Average the source estimates within each label of the cortical parcellation
# and each sub structure contained in the src space

src = inverse_operator['src']

label_ts = mne.extract_label_time_course(
    [stc], labels_parc, src, mode='mean', allow_empty=True)

# plot the times series of 2 labels
fig, axes = plt.subplots(1)
axes.plot(1e3 * stc.times, label_ts[0][0, :], 'k', label='bankssts-lh')
axes.plot(1e3 * stc.times, label_ts[0][71, :].T, 'r', label='Brain-stem')
axes.set(xlabel='Time (ms)', ylabel='MNE current (nAm)')
axes.legend()
mne.viz.tight_layout()

# %%
stcsNormal = apply_inverse_epochs(epochs, inverse_operator, lambda2,
                                method, pick_ori="normal",
                                return_generator=True)

# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi
labels = mne.read_labels_from_annot('subject_1', parc='aparc.DKTatlas40',
                                    subjects_dir=subjects_dir)


# Average the source estimates within each label using sign-flips to reduce
# signal cancellations, also here we return a generator
src = inverse_operator['src']
labelTsNormal = mne.extract_label_time_course(stcsNormal, labels, src,
                                            mode='mean_flip',
                                            return_generator=False)


# %%
from nitime import TimeSeries
from nitime.analysis import MTCoherenceAnalyzer
from nitime.viz import drawmatrix_channels

f_up = 13  # upper limit
f_lw = 8  # lower limit

cohMatrixNormal = np.empty([np.shape(labelTsNormal)[1], np.shape(labelTsNormal)[1],
                          np.shape(labelTsNormal)[0]])

labels_name = []

# Load data
fname_inv = mne_folder + "%s-inv.fif" % subject
inv = mne.minimum_norm.read_inverse_operator(fname_inv)
fname_evoked = epochs_folder + "%s_filtered_ica_mc_tsss-ave.fif" % subject
evokeds = mne.read_evokeds(fname_evoked, baseline=(None, 0))
src = mne.read_source_spaces(mne_folder + "%s-oct-6-src.fif" % subject)

for evk in evokeds:
    stc = apply_inverse(evk, inv, lambda2=lambda2,  
                        method=method)
    exec("stc_%s_%s = stc" % (subject, evk.comment))


# src = mne.read_source_spaces(mne_folder + "%s-oct-6-src.fif" % subject)
labels = mne.read_labels_from_annot(subject, parc='PALS_B12_Lobes',
                                    # regexp="Bro",
                                    subjects_dir=subjects_dir)
labels_occ = [labels[9], labels[10], labels[9]+labels[10]]

lbl_ent_left = mne.extract_label_time_course(stc_0006_ent_left,
                                             labels=[labels[9]],
                                             src=src,
                                             mode="pca_flip")

lbl_ctl_left = mne.extract_label_time_course(stc_0006_ctl_left,
                                             labels=[labels[9]],
                                             src=src,
                                                 mode="pca_flip")

    else:
        label = avg_label[1].morph(subject_to=s)

    epochs_fname = home + '/data/meg/stop/parsed/%s_stop_parsed_matched_clean_BP1-100_DS300-epo.fif.gz' % s
    epochs = mne.read_epochs(epochs_fname, proj=True)
    fwd_fname = home + '/data/meg/stop/%s_task-5-fwd.fif' % s
    fwd = mne.read_forward_solution(fwd_fname, surf_ori=True)

    # calculate source power estimates for the whole brain
    # quick hack in tmax ot make it the same length as btmax
    data_csds = compute_epochs_csd(epochs[cond], mode='multitaper',
                                   tmin=tmin, tmax=tmax + btmax,
                                   fmin=band[0], fmax=band[1],
                                   fsum=False)
    noise_csds = compute_epochs_csd(epochs[cond], mode='multitaper',
                                    tmin=btmin, tmax=btmax,
                                    fmin=band[0], fmax=band[1],
                                    fsum=False)
    stc = dics_source_power(epochs.info, fwd, noise_csds, data_csds)
    ts = mne.extract_label_time_course(stc, label, fwd['src'])
    data.append(ts)

# export one CSV file
fname = out_dir + '%s_%s_%02dto%02d_tmin%.2f.csv' % (cond, roi_pretty, band[0],
                                                     band[1], tmin)
fid = open(fname, 'w')
fid.write('subj,power\n')
for j, d in enumerate(data):
    fid.write('%s,' % subjs[j] + ','.join(['%e' % t for t in d[0]]) + '\n')
fid.close()

# labels = mne.read_labels_from_annot('subject_1', parc='aparc.DKTatlas40',
#                                     subjects_dir=subjects_dir)

for cond in epochs.event_id.keys():
    stcs = apply_inverse_epochs(epochs[cond], inverse_operator, lambda2,
                                method, pick_ori="normal")
    exec("stcs_%s = stcs" % cond)

labels_name = [label.name for label in labels_occ]

for label in labels_occ:
    labels_name += [label.name]

# Extract  time series
ts_ctl_left = mne.extract_label_time_course(stcs_ctl_left,
                                            labels_occ,
                                            src=inverse_operator["src"],
                                            mode = "mean_flip")

ts_ent_left = mne.extract_label_time_course(stcs_ent_left,
                                            labels_occ,
                                            src=inverse_operator["src"],
                                            mode = "mean_flip")

stcs_all_left = stcs_ctl_left + stcs_ent_left
ts_all_left = np.asarray(mne.extract_label_time_course(stcs_all_left,
                                            labels_occ,
                                            src=inverse_operator["src"],
                                            mode = "mean_flip"))

number_of_permutations = 2000
index = np.arange(0, len(ts_all_left))

    epochs = mne.read_epochs(epochs_folder +
                             "%s_ds_filtered_ica_mc_tsss-epo.fif" % subject)
    # epochs.resample(250, n_jobs=4)

    for condition in conditions:
        stcs = apply_inverse_epochs(epochs[condition],
                                    inverse_operator,
                                    lambda2,
                                    method,
                                    pick_ori="normal")

        for label in labels_occ:
            label_ts = []
            for j in range(len(stcs)):
                label_ts.append(mne.extract_label_time_course(stcs[j],
                                                              labels=label,
                                                              src=src,
                                                              mode="mean_flip"))

            label_ts = np.squeeze(np.asarray(label_ts))

            tfr = cwt_morlet(label_ts, epochs.info["sfreq"], freqs,
                             use_fft=True, n_cycles=n_cycle)

            np.save(tf_folder + "%s_%s_%s_MNE-tfr" % (subject, condition,
                                                      label.name),
                    tfr)

        del stcs
        del tfr

    del epochs

    for first in range(start, stop, step):
        last = first + step
        if last >= stop:
            last = stop
        raw_segment = raw_data[:, first:last]
        mu += raw_segment.sum(axis=1)
        data += np.dot(raw_segment, raw_segment.T)
        n_samples += raw_segment.shape[1]
    mu /= n_samples
    data -= n_samples * mu[:, None] * mu[None, :]
    data /= (n_samples - 1.0)
    ch_names = [raw.info['ch_names'][k] for k in picks]
    data_cov = mne.Covariance(None)
    data_cov.update(kind=mne.fiff.FIFF.FIFFV_MNE_NOISE_COV, diag=False, 
        dim=len(data), names=ch_names, data=data, 
        projs=cp.deepcopy(raw.info['projs']), bads=raw.info['bads'], 
        nfree=n_samples, eig=None, eigvec=None)

    noise_cov = mne.compute_raw_data_covariance(er_raw)
    # note that MNE reads CTF data as magnetometers!
    noise_cov = mne.cov.regularize(noise_cov, raw.info, mag=noise_reg)
    events = fg.get_good_events(markers[subj], time, window_length)

    epochs = mne.Epochs(raw, events, None, 0, window_length, preload=True, baseline=None, detrend=0, picks=picks)
    stcs = mne.beamformer.lcmv_epochs(epochs, forward, noise_cov.as_diag(), data_cov, reg=data_reg, pick_ori='max-power')

    for net in avg_intersects:
        subj_labels = [label.morph('fsaverage',subj) for label in net]
        label_ts = mne.extract_label_time_course(stcs, subj_labels, forward['src'], mode=label_mode, allow_empty=True)
        con, freqs, times, n_epochs, n_tapers = mne.connectivity.spectral_connectivity(label_ts, method=method, mode='multitaper', sfreq=raw.info['sfreq'], fmin=[1,4,8,13,30], fmax=[4,8,13,30,50], faverage=True, n_jobs=3, mt_adaptive=False)
        np.save(dir_out + subj + '-' + net[0].name + '-' + label_mode +'-' + '-'.join(method), con)

                                         fixed=False)


stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, inv_method,
                            pick_ori=None, return_generator=True)

# Get labels for FreeSurfer 'aparc' cortical parcellation with 34 labels/hemi
labels_parc = mne.read_labels_from_annot(subject, parc=parc,
                                         subjects_dir=subjects_dir)

# Average the source estimates within each label of the cortical parcellation
# and each sub structures contained in the src space
# If mode = 'mean_flip' this option is used only for the cortical label
src = inverse_operator['src']
label_ts = mne.extract_label_time_course(stcs, labels_parc, src,
                                         mode='mean_flip',
                                         allow_empty=True,
                                         return_generator=False)

# We compute the connectivity in the alpha band and plot it using a circular
# graph layout
fmin = 8.
fmax = 13.
sfreq = raw.info['sfreq']  # the sampling frequency
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    label_ts, method='pli', mode='multitaper', sfreq=sfreq, fmin=fmin,
    fmax=fmax, faverage=True, mt_adaptive=True, n_jobs=1)

# We create a list of Label containing also the sub structures
labels_aseg = mne.get_volume_labels_from_src(src, subject, subjects_dir)
labels = labels_parc + labels_aseg