Last active
July 7, 2016 11:37
-
-
Save cosacog/4e7e79c4538ec7cb14e5 to your computer and use it in GitHub Desktop.
labelをいじるfunction、label内の時系列データのうち活動の多い方から10%のvertexの平均をプロットするfunction
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def half_ydir_label(label_orig): | |
''' | |
要mne等 | |
labelをy軸方向に半分にするスクリプト | |
params: label_orig. 元のlabel | |
return: label_half. y軸方向に半分になったlabel | |
''' | |
ypos_label = label_orig.pos[:,1] | |
idx_half_label = np.where(ypos_label >= np.median(ypos_label)) # example: y positions are anterior | |
pos_label_half = label_orig.pos[idx_half_label,:] | |
vertices_half = label_orig.vertices[idx_half_label] | |
values_half = label_orig.values[idx_half_label] | |
name_half = label_orig.name[:-3] + '_half'+label_orig.name[-3:] | |
label_half = label_orig.copy() | |
label_half.pos = pos_label_half | |
label_half.vertices = vertices_half | |
label_half.values = values_half | |
label_half.name= name_half | |
return label_half | |
def plot_conditions_in_label_10percentile(list_path_stc, list_dict_labels, percentile_vertices, t_range, src, subjects_dir, is_label_smooth): | |
''' | |
list_path_stcのファイルから読み込んだデータを重ね書き(plot)するfunction | |
list_dict_labelsに入ったラベルの領域毎にsubplotでプロット | |
適当(0.1-0.2 sec)の最大振幅の時刻を基準に10 percentileのverticesを取得して平均振幅を描画 | |
<params> list_path_stc: stc fileのlist | |
<params> list_dict_labels: dict(label と title)のlist | |
<params> percentile_vertices: label内verticesから読み込む数. percentile表示(e.g. 0.1で10%) | |
<params> t_range: 最大の時刻, vertexを検出する時間幅. e.g. (0.1, 0.2) | |
<params> src: いわゆるsrc. ただし中身は mne.add_source_space_distances 後のsrcにしとかないと計算が20分くらいかかる | |
<params> subjects_dir: いわゆるsubjects_dir. | |
<params> is_label_smooth: True or False. 通常Trueでよいが、Falseでstcの中のvertexの位置が具体的にわかる | |
return: subplotで各labelの中のstcデータを平均してプロット。それぞれlabel内で大きい方から10 percentileの頂点を選ぶ | |
<out> array_mean_percentile: 平均波形のarray, labelとファイルの数だけ返って来る | |
<out> header_csv: 各列のlabel,ファイル名がわかるようにしたheader | |
<out> list_label: 検出したlabelのlist, labelとファイルの順番は上記header_csvに準じる | |
<out> list_latency: 検出した頂点潜時のlist, 順番は上に同じ | |
''' | |
# settings | |
time_range4peak_detect = t_range | |
colors = ['red','black','blue','green'] # plot colors: ファイルの種類に対応する | |
# main | |
import numpy as np | |
import matplotlib.pyplot as plt | |
import mne | |
header_csv = 'time' | |
idx = 0 | |
list_label = [] | |
list_latency = [] | |
for idx_plot, dict_label_title in enumerate(list_dict_labels): | |
label = dict_label_title['label'] | |
title = dict_label_title['title'] | |
idx_hemi = 0 if label.hemi == 'lh' else 1 | |
plt.subplot(int(np.ceil(len(list_dict_labels)/2)), 2, idx_plot+1) | |
legend_text = [os.path.basename(path)[:-4] for path in list_path_stc] | |
for ii, path_stc in enumerate(list_path_stc): | |
stc = mne.read_source_estimate(path_stc) | |
stc_in_label = stc.in_label(label) | |
n_vertices = stc_in_label.data.shape[0] | |
stc_data = stc_in_label.data | |
idx_toi = np.where((stc_in_label.times >= time_range4peak_detect[0]) & | |
(stc_in_label.times <= time_range4peak_detect[1]))[0] | |
# detect peak latency and vertex | |
# idx_vertex_peak, idx_time_peak = np.where(stc_data[:,idx_toi] == np.max(stc_data[:,idx_toi])) | |
stc_mean_data_toi = np.mean(stc_data[:,idx_toi], axis=0) | |
idx_time_peak = np.where(stc_mean_data_toi == np.max(stc_mean_data_toi))[0] # 全vertex平均の最大時刻のindex | |
idx_vertices_peak = np.argsort(stc_data[:,idx_time_peak[0]],axis=0)[-int(np.ceil(n_vertices * percentile_vertices)):] | |
mean_stc_data_percentile = np.mean(stc_data[idx_vertices_peak,],axis=0) # | |
# detected latency | |
latency_detected = stc_in_label.times[idx_toi[0] + idx_time_peak][0] | |
list_latency.append(latency_detected) | |
# plot | |
plt.plot(stc.times, mean_stc_data_percentile, color = colors[ii]) | |
plt.xlim(stc.times[0],stc.times[-1]) | |
plt.title(title) | |
# data concatenate: 三項演算子 | |
array_mean_percentile = mean_stc_data_percentile if idx == 0 else np.r_['1,2,0', array_mean_percentile, mean_stc_data_percentile] | |
# csv header | |
header_csv += ',' + title + '-' + os.path.basename(path_stc)[:-7] | |
# label | |
stc_in_label.data[:,:] = 0 | |
stc_in_label.data[idx_vertices_peak,0]=1 #ここだけlabelになる | |
label_peak_both = mne.stc_to_label(stc_in_label, src=src, smooth=is_label_smooth, connected=False, subjects_dir=subjects_dir) | |
label_peak_hemi = label_peak_both[idx_hemi] # rh or lh only | |
label_peak_hemi.subject = stc.subject | |
label_peak_hemi.name = os.path.basename(path_stc)[:-7] + ' ' + title | |
list_label.append(label_peak_hemi) | |
idx += 1 | |
plt.legend(legend_text) | |
array_mean_percentile = np.r_['1,2,0', stc.times, array_mean_percentile] | |
return array_mean_percentile, header_csv, list_label, list_latency | |
def find_label_in_annot(label_name, labels): | |
''' | |
mne.read_labels_from_annotで読みだしたlabelsからほしい名前のlabelの番号を取り出す | |
<params> | |
label_name: ほしいlabelの名前の一部'fusiform'とか | |
labels: mne.read_labels_from_annotで読みだしたlabels | |
<return> | |
labels_sub: 取り出したlabel | |
idx_labels: 取り出したlabelの番号.labels[idx_labels]とかするとほしいlabelが出せるはず | |
''' | |
import numpy as np | |
labels_sub = [] | |
idx_labels = np.array([]) | |
ii = 0 | |
for label in labels: | |
idx = label.name.find(label_name) | |
if idx > -1: | |
labels_sub.append(label) | |
idx_labels = np.append(idx_labels, ii) | |
ii += 1 | |
return labels_sub, idx_labels | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
.annotファイルからほしいlabelのindexを探すfunction: find_label_in_annot
要mne, numpy
ほしいエリアにどのような名前のラベルがついてるかはfreeviewあたりで調べましょう.