载入fNIRS数据：https://mne.tools/stable/auto_tutorials/io/plot_30_reading_fnirs_data.html#sphx-glr-auto-tutorials-io-plot-30-reading-fnirs-data-py
fNIRS预处理：https://mne.tools/stable/auto_tutorials/preprocessing/plot_70_fnirs_processing.html#sphx-glr-auto-tutorials-preprocessing-plot-70-fnirs-processing-py
1.载入数据，我的数据是snirf格式的，使用
# preload=False不加载数据，preload=True加载数据
raw_od= mne.io.read_raw_snirf("数据路径", preload=True)载入数据。若数据格式是NIRx记录的，可以用mne.io.read_raw_nirx()载入数据，该函数只在NIRScout上测试过可行。
注：mne的0.20.7版本还无法使用read_raw_snirf函数，conda install、update命令也无法升级，直接用pip install -U mne进行安装或升级到最新版本。
读取数据过程中出现TypeError: iteration over a 0-d array错误
原因是，sources和detectors的Labels没有成功读入，可能SNIRF格式的种类太多，我的数据中没有这部分数据。
sources = np.array(dat.get('nirs/probe/sourceLabels'))
detectors = np.array(dat.get('nirs/probe/detectorLabels'))好在所有数据channel都是固定的，手动将Labels设置上
read_raw_snirf()# sources = np.array(dat.get('nirs/probe/sourceLabels'))
# detectors = np.array(dat.get('nirs/probe/detectorLabels'))
# sources = [s.decode('UTF-8') for s in sources]
# detectors = [d.decode('UTF-8') for d in detectors]
# 根据探头数量设置
sources = ['S1','S2','S3','S4','S5','S6','S7','S8','S9','S10','S11','S12','S13','S14','S15','S16']
detectors = ['D1','D2','D3','D4','D5','D6','D7','D8','D9','D10','D11','D12','D13','D14','D15','D16']数据格式如果有异常，可以根据官方提供的数据排查问题，并手动更改数据
2.根据https://mne.tools/stable/auto_tutorials/preprocessing/plot_70_fnirs_processing.html#sphx-glr-auto-tutorials-preprocessing-plot-70-fnirs-processing-py提供代码对数据进行预处理。我已经在MATLAB中预处理过了，直接将数据进行替换
from scipy.io import loadmatraw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)
nirs_HbO = loadmat("HbO路径")
nirs_HbR = loadmat("HbR路径")
raw_haemo._data[list(range(0,96,2)), :] = nirs_HbO['oxyData'].transpose()
raw_haemo._data[list(range(1,96,2)), :] = nirs_HbR['dxyData'].transpose()发现标签数据（raw_haemo._annotations.description）异常，手动读入标签信息，进行替换（参考：https://www.cnblogs.com/hackpig/p/8215786.html）
import numpy as npdef get_ananotation ():filename = 'trigger信息文件路径' labels = []with open(filename, 'r') as file_to_read:while True:lines = file_to_read.readline() # 整行读取数据if not lines:breakpassa, b, c, d, label, e = [(i) for i in lines.split(";")] #我的数据是按“；”分隔的labels.append(label) # 添加新读取的数据passlabels = np.array(labels) # 将数据从list类型转换为array类型。passreturn labelsraw_haemo._annotations.description = get_ananotation()
4.绘图所需数据准备
#event_id根据raw_haemo._annotations.description的格式和事件数进行更改
events, _ = mne.events_from_annotations(raw_haemo, event_id={'1': 1,'2': 2,'3': 3,'4': 4})
event_dict = {'hand': 1, 'wrist': 2, 'shoulder': 3, 'rest': 4}
tmin, tmax = -5, 30 # 所需的时间段在标签打下前5s到标签打下后30sepochs = mne.Epochs(raw_haemo, events, event_id=event_dict,tmin=tmin, tmax=tmax,reject_by_annotation=True,proj=True, baseline=(None, 0), preload=True,detrend=None, verbose=True)# nrows和ncols分别表示行数和列数，列数应算上colorbar的那一列，gridspec_kw最后的0.1为colorbar的比重
fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(9, 5),gridspec_kw=dict(width_ratios=[1, 1, 1, 1, 0.1]))
vmin, vmax, ts = -8, 8, 9.0 # vmin, vmax表示colorbar的范围，ts表示绘制的是9.0s时刻的脑地形图
topomap_args = dict(extrapolate='local')evoked_hand = epochs['hand'].average()
evoked_wrist = epochs['wrist'].average()
evoked_shoulder = epochs['shoulder'].average()
evoked_rest = epochs['rest'].average()
5.绘制脑地形图
evoked_hand.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 0],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_hand.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 0],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_wrist.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 1],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_wrist.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 1],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_shoulder.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 2],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_shoulder.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 2],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
# 最后一列的两张图，colorbar=True，axes后有个‘：’
evoked_rest.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 3:],vmin=vmin, vmax=vmax, colorbar=True,**topomap_args)
evoked_rest.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 3:],vmin=vmin, vmax=vmax, colorbar=True,**topomap_args)
6.设置标题
for column, condition in enumerate(['hand', 'wrist', 'shoulder', 'rest']):for row, chroma in enumerate(['HbO', 'HbR']):axes[row, column].set_title('{}: {}'.format(chroma, condition))
fig.tight_layout()
7.官方数据及全部代码
import os
import numpy as np
import matplotlib.pyplot as plt
from itertools import compress
import mnefnirs_data_folder = mne.datasets.fnirs_motor.data_path()
fnirs_cw_amplitude_dir = os.path.join(fnirs_data_folder, 'Participant-1')
raw_intensity = mne.io.read_raw_nirx(fnirs_cw_amplitude_dir, verbose=True)
raw_intensity.load_data()picks = mne.pick_types(raw_intensity.info, meg=False, fnirs=True)
dists = mne.preprocessing.nirs.source_detector_distances(raw_intensity.info, picks=picks)
raw_intensity.pick(picks[dists > 0.01])raw_od = mne.preprocessing.nirs.optical_density(raw_intensity)sci = mne.preprocessing.nirs.scalp_coupling_index(raw_od)raw_od.info['bads'] = list(compress(raw_od.ch_names, sci < 0.5))raw_haemo = mne.preprocessing.nirs.beer_lambert_law(raw_od)events, _ = mne.events_from_annotations(raw_haemo, event_id={'1.0': 1,'2.0': 2,'3.0': 3})
event_dict = {'Control': 1, 'Tapping/Left': 2, 'Tapping/Right': 3}reject_criteria = dict(hbo=80e-6)
tmin, tmax = -5, 15epochs = mne.Epochs(raw_haemo, events, event_id=event_dict,tmin=tmin, tmax=tmax,reject=reject_criteria, reject_by_annotation=True,proj=True, baseline=(None, 0), preload=True,detrend=None, verbose=True)fig, axes = plt.subplots(nrows=2, ncols=4, figsize=(9, 5),gridspec_kw=dict(width_ratios=[1, 1, 1, 0.1]))
vmin, vmax, ts = -8, 8, 9.0evoked_left = epochs['Tapping/Left'].average()
evoked_right = epochs['Tapping/Right'].average()evoked_left.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 0],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_left.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 0],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_right.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 1],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)
evoked_right.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 1],vmin=vmin, vmax=vmax, colorbar=False,**topomap_args)evoked_diff = mne.combine_evoked([evoked_left, evoked_right], weights=[1, -1])evoked_diff.plot_topomap(ch_type='hbo', times=ts, axes=axes[0, 2:],vmin=vmin, vmax=vmax, colorbar=True,**topomap_args)
evoked_diff.plot_topomap(ch_type='hbr', times=ts, axes=axes[1, 2:],vmin=vmin, vmax=vmax, colorbar=True,**topomap_args)for column, condition in enumerate(['Tapping Left', 'Tapping Right', 'Left-Right']):for row, chroma in enumerate(['HbO', 'HbR']):axes[row, column].set_title('{}: {}'.format(chroma, condition))
fig.tight_layout()
8.官方绘制脑地形图
支付宝扫一扫
微信扫一扫
