基于python mne库构造自定义fNIRS数据并可视化地形图

在科研中遇到需要可视化fNIRS通道的重要性，参考了网上已有的一篇文章，发现只是导入元数据进行替换来实现的，并不符合自己目标（不是需要可视化原始数据，而是需要可视化通道间重要性，每个通道值为0-1，所有通道加起来合为1），因此花了不少时间参照mne库完成了实现。
实现效果：
1、自定义每个通道的位置，无需元数据位置信息
2、自定义通道数据，可以绘制相对重要性而非元数据电压分布图
3、替换了强度标识，而非具体电压大小
代码如下：

点击查看代码

import numpy as np
import mne
import yaml
import matplotlib.pyplot as plt

from mne.channels import make_dig_montage

def get_1020_positions(fnirs_to_1020_mapping):
    """
    根据10/20国际标准体系获取fNIRS电极的具体三维坐标
    参数:
    ----------
    fnirs_to_1020_mapping : dict
        映射字典，格式: {'S1': 'FC3', 'S2': 'FCz', 'D1': 'FC1', ...}
        其中键是fNIRS电极名，值是10/20系统中的标准电极名
    
    返回:
    ----------
    tuple: (source_positions, detector_positions)
        - source_positions: 光源位置字典 {'S1': array([x,y,z]), ...}
        - detector_positions: 探测器位置字典 {'D1': array([x,y,z]), ...}

    """
    
    # 1. 获取标准10/20系统的所有电极位置
    montage_1020 = mne.channels.make_standard_montage('standard_1020')
    positions_1020 = montage_1020.get_positions()['ch_pos']
    
    # 2. 初始化返回字典
    source_positions = {}
    detector_positions = {}
    missing_electrodes = []
    
    # 3. 遍历映射字典，获取每个电极的坐标
    for fnirs_name, eeg_name in fnirs_to_1020_mapping.items():
        if eeg_name in positions_1020:
            position = positions_1020[eeg_name]
            
            # 根据电极名的首字母判断是光源还是探测器
            if fnirs_name.startswith('S'):
                source_positions[fnirs_name] = position
            elif fnirs_name.startswith('D'):
                detector_positions[fnirs_name] = position
        else:
            missing_electrodes.append((fnirs_name, eeg_name))
    if missing_electrodes:
        print(f"\n⚠ 警告: 以下电极不在标准10/20系统中:")
        for fnirs_name, eeg_name in missing_electrodes:
            print(f"  {fnirs_name} -> {eeg_name}")
    return source_positions, detector_positions

def create_fNIRS_data(data,config_path):

    """根据numpy data创建fNIRS数据 data维度为[channels*2,times]"""
    # 1. 定义通道配置,包含source和detector对应位置，以及通道定义
 
    with open(config_path, 'r', encoding='utf-8') as f:
        config = yaml.safe_load(f)
    channel_config = config['channel_config']
    
    # 2. 获取光源和探测器的3D位置
    fnirs_to_1020_mapping = config['fnirs_to_1020_mapping']
    source_positions, detector_positions = get_1020_positions(fnirs_to_1020_mapping)
    '''
    例子：
        source_positions = {
        'S1': np.array([0.03, 0.08, 0.05]),   # 前额左
        'S2': np.array([-0.03, 0.08, 0.05]),  # 前额右
    }
    detector_positions = {
        'D1': np.array([0.02, 0.05, 0.05]),   # 前额左探测器
        'D2': np.array([0.00, 0.05, 0.05]),   # 前额中探测器
        'D3': np.array([-0.02, 0.05, 0.05]),  # 前额右探测器
    }
    '''

    # 3. 创建通道名
    ch_names = []
    for src, det, _ in channel_config:
        ch_names.append(f"{src}_{det} hbo")  # 760nm波长
        ch_names.append(f"{src}_{det} hbr")  # 850nm波长
    # n_channels = len(ch_names)
    # n_times = 500
    sfreq = 10.0 
    # 5. 创建info对象，这是mne中数据的元数据，包含基本信息，主要是数据每维度的通道名，type，要对应
    ch_types = []
    for ch_name in ch_names:
        if 'hbo' in ch_name:
            ch_types.append('hbo')  
        elif 'hbr' in ch_name:
            ch_types.append('hbr')  
    info = mne.create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)

    # 7. 创建montage，这是电极布局，eeg数据通道直接导入，fnirs数据通过提供3D数据也能创建
    ch_pos = {**source_positions, **detector_positions}
    montage = make_dig_montage(ch_pos=ch_pos, coord_frame='head')
    info.set_montage(montage)
    # 8. 创建Raw对象，数据+信息结构
    raw = mne.io.RawArray(data, info, verbose=False)
    # 创建Evoked对象用于地形图，和raw对象区别是分批次，便于取平均
    evoked = mne.EvokedArray(data, info, tmin=0, comment='fNIRS data')
  
    return evoked,ch_names,raw, channel_config, source_positions, detector_positions

def plot_fNIRS_layout(raw, source_positions, detector_positions):
    """绘制fNIRS电极位置布局"""
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    # 1. 2D传感器布局
    ax1 = axes[0]
    raw.plot_sensors(axes=ax1, show_names=False)
    # 调整电极点大小：找到所有的散点集合并调整大小
    for collection in ax1.collections:
        if hasattr(collection, 'get_sizes'):
            # 获取当前的散点大小
            current_sizes = collection.get_sizes()
            if len(current_sizes) > 0:
                # 设置为指定的小尺寸
                collection.set_sizes([current_sizes/5])
    ax1.set_title("fNIRS sensors")
    # 2. 在坐标轴上绘制光源和探测器位置
    ax2 = axes[1]  
    # 绘制光源
    for src, pos in source_positions.items():
        ax2.scatter(pos[0], pos[1], c='red', s=200, marker='o', label='source' if src == 'S1' else "")
        ax2.text(pos[0], pos[1], src, fontsize=12, ha='center', va='bottom')
    # 绘制探测器
    for det, pos in detector_positions.items():
        ax2.scatter(pos[0], pos[1], c='blue', s=200, marker='s', label='detector' if det == 'D1' else "")
        ax2.text(pos[0], pos[1], det, fontsize=12, ha='center', va='bottom')
    ax2.set_xlabel('X (m)')
    ax2.set_ylabel('Y (m)')
    ax2.set_title('source and detector position')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    ax2.axis('equal') 
    plt.tight_layout()
    return fig

def plot_custom_topmap(evoked, ch_names, time_point=0.0, figsize=(10, 8)):
    # 筛选HBO和HBR通道
    hbo_picks = [i for i, ch in enumerate(ch_names) if 'hbo' in ch]
    hbr_picks = [i for i, ch in enumerate(ch_names) if 'hbr' in ch]
    
    if not hbo_picks or not hbr_picks:
        raise ValueError("未找到HBO或HBR通道")
    
    # 创建图形
    fig, axes = plt.subplots(1, 2, figsize=figsize)
    
    # 绘制HBO地形图
    hbo_data = evoked.data[hbo_picks, 0]  # 取第一个时间点
    hbo_info = mne.pick_info(evoked.info, hbo_picks)
    
    im_hbo, _ = mne.viz.plot_topomap(
        data=hbo_data,
        pos=hbo_info,
        axes=axes[0],
        show=False,
        cmap='RdBu_r',
        vmin=-5,
        vmax=5,
        sensors=True,
        contours=False
    )
    axes[0].set_title('HbO Topography', fontsize=12, fontweight='bold')
    plt.colorbar(im_hbo, ax=axes[0], shrink=0.8)
    
    # 绘制HBR地形图
    hbr_data = evoked.data[hbr_picks, 0]
    hbr_info = mne.pick_info(evoked.info, hbr_picks)
    
    im_hbr, _ = mne.viz.plot_topomap(
        data=hbr_data,
        pos=hbr_info,
        axes=axes[1],
        show=False,
        cmap='RdBu_r',
        vmin=-5,
        vmax=5,
        sensors=True,
        contours=False
    )
    axes[1].set_title('HbR Topography', fontsize=12, fontweight='bold')
    plt.colorbar(im_hbr, ax=axes[1], shrink=0.8)
    
    plt.suptitle(f'fNIRS Topography at t={time_point}s', fontsize=14, fontweight='bold')
    plt.tight_layout() 
    return fig

def plot_fNIRS_topmap(data,config_path):
    evoked, ch_names, raw, channel_config, source_positions, detector_positions = create_fNIRS_data(data,config_path)
    # # 绘制布局
    fig_layout = plot_fNIRS_layout(raw, source_positions, detector_positions)
    fig_layout.savefig("fNIRS_layout.png", dpi=300, bbox_inches='tight')
    
    #绘制地形图
    # 修改后的topomap_args
    topomap_args = dict(
        extrapolate='local',
        sphere=(0, 0, 0, 0.115),
    #   image_interp='bicubic',
        contours=True,        # 禁用等高线
        outlines='head',
        sensors=True,
        colorbar=False,
        res=64,
        show=False,

    )
    
    im = evoked.plot_topomap(ch_type='hbr', 
                                times=0.,          
                            **topomap_args)
    fig = plt.gcf()  # 获取当前图形
    ax = fig.axes[0]
    
    # 从 axes 中找到图像对象
    image_obj = None
    for child in ax.get_children():
        if hasattr(child, 'get_array'):
            image_obj = child
            break

    # 创建颜色条
    if image_obj is not None:
        cbar = plt.colorbar(image_obj, ax=ax, 
                        fraction=0.046,
                        pad=0.04,
                        shrink=0.8)
    # 移除所有数字刻度
    cbar.set_ticks([])  
    
    # 在颜色条顶部添加"强"标签
    cbar.ax.text(0.5, 1.05, 'high',  # 位置：水平居中，垂直稍上方
                transform=cbar.ax.transAxes,  # 使用相对坐标
                ha='center', va='bottom',  # 水平居中，垂直底部对齐
                fontsize=6,  # 字体大小
                fontweight='bold',  # 粗体
                color='black')
    
    # 在颜色条底部添加"弱"标签
    cbar.ax.text(0.5, -0.05, 'low',  # 位置：水平居中，垂直稍下方
                transform=cbar.ax.transAxes,
                ha='center', va='top',  # 水平居中，垂直顶部对齐
                fontsize=6,
                fontweight='bold',
                color='black')

    fig.savefig("fig_topography.png", dpi=300, bbox_inches='tight')

if __name__ == "__main__":
    # 创建数据
    n_hbo = 36  
    n_hbr = 36 
    
    # 生成随机数据
    random_hbo = np.random.rand(n_hbo, 1)  
    random_hbr = np.random.rand(n_hbr, 1)  

    # 归一化，使得数据为0-1
    hbo_normalized = random_hbo / random_hbo.sum() 
    hbr_normalized = random_hbr / random_hbr.sum() 
    
    # 合并数据，按HBO-HBR交替排列
    data = np.zeros((n_hbo+n_hbr, 1))
    for i in range(n_hbo):
        data[2*i, 0] = hbo_normalized[i, 0]      # HBO
        data[2*i + 1, 0] = hbr_normalized[i, 0]  # HBR

    config_path = 'public_dataset_config.yaml'
    plot_fNIRS_topmap(data,config_path)

public_dataset_config.yaml如下：

点击查看代码

channel_config:
  - [S2, D2, C1]  # 通道1
  - [S3, D2, C2]  # 通道2
  - [S3, D1, C3]  # 通道3
  - [S1, D2, C4]  # 通道4
  - [S1, D1, C5]  # 通道5
  - [S1, D3, C6]  # 通道6
  - [S5, D1, C7]  # 通道7
  - [S5, D3, C8]  # 通道8
  - [S4, D3, C9]  # 通道9
  - [S14, D14, C10]  # 通道10
  - [S14, D15, C11]  # 通道11
  - [S14, D16, C12]  # 通道12
  - [S8, D8, C13]  # 通道13
  - [S8, D6, C14]  # 通道14
  - [S8, D4, C15]  # 通道15
  - [S7, D6, C16]  # 通道16
  - [S7, D4, C17]  # 通道17
  - [S7, D5, C18]  # 通道18
  - [S9, D8, C19]  # 通道19
  - [S9, D4, C20]  # 通道20
  - [S9, D7, C21]  # 通道21
  - [S6, D4, C22]  # 通道22
  - [S6, D5, C23]  # 通道23
  - [S6, D7, C24]  # 通道24
  - [S10, D10, C25]  # 通道25
  - [S10, D12, C26]  # 通道26
  - [S10, D9, C27]  # 通道27
  - [S11, D10, C28]  # 通道28
  - [S11, D9, C29]  # 通道29
  - [S11, D11, C30]  # 通道30
  - [S13, D13, C31]  # 通道31
  - [S13, D12, C32]  # 通道32
  - [S13, D9, C33]  # 通道33
  - [S12, D13, C34]  # 通道34
  - [S12, D9, C35]  # 通道35
  - [S12, D11, C36]  # 通道36

fnirs_to_1020_mapping:
  S1: Fpz     
  S2: AF7     
  S3: AF3
  S4: AF8
  S5: AF4
  S6: C1
  S7: FC3
  S8: C5
  S9: CP3     
  S10: C2     
  S11: FC4
  S12: C6
  S13: CP4
  S14: Oz

  D1: AFz     
  D2: Fp1    
  D3: Fp2
  D4: C3  
  D5: FC1   
  D6: FC5
  D7: CP1
  D8: CP5
  D9: C4     
  D10: FC2    
  D11: FC6  
  D12: CP2   
  D13: CP6
  D14: POz
  D15: O1
  D16: O2