简介针对脑影像分析中需要从HCP-MMP1.0图谱提取特定脑区掩膜的需求这份MATLAB脚本工具可直接读取MNI空间下的体积图谱依据mmp.csv中的ROI编号快速生成二值掩膜文件。脚本支持左右半球选择、单个/多个ROI分别保存或合并保存并可选择对生成的掩膜进行扩张适合使用HCP分区开展科研或教学的中级MATLAB用户。资源包共7个文件、约834KB以两个MATLAB函数为主包含地图集nii文件、ROI查找表csv、测试样例、README与license文件结构精简借助测试nii可验证输出结果。已有546人学习下载可用于fMRI预处理后自制ROI、基于MMP1.0的分区统计、功能连接分析等场景。脚本省去手动逐体素标记脑区的繁琐流程也让重复提取不同ROI变得更加高效。1. 从一个需求说起做了三年fMRI分析我在提取脑区mask上栽过两次跟头处理静息态fMRI数据时最烦的事情之一就是从标准的脑图谱里抠出自己感兴趣的那个区域。HCP-MMP1.0是近几年用的比较多的皮层图谱把大脑皮层分成了360个区域每个区域有独立的编号和名称。但问题是原生发布的数据格式是拆分好的文件你想用某个特定脑区做seed-based connectivity分析时得自己想办法把它变成二值mask——这活儿看起来简单真上手就知道坑不少。有的区域左右半球都有有的只有单侧有的编号在拆分文件里是乱的直接赋值出来的mask往往带上了别的区域。这篇笔记要讲的rsamatlab代码就是解决这个问题的它在rsamatlab这套静息态分析工具箱的基础上专门加了一个提取HCP-MMP1.0特定二值mask的入口输入区域名或编号输出一个干净、对齐好的nii格式mask文件。适合用MATLAB做脑影像分析、又不想在mask生成上反复折腾的开发者和研究者。2. 技术选型与结构认知为什么是rsamatlab HCP-MMP1.0这套组合2.1 HCP-MMP1.0图谱不是“一张图”而是一组拆分文件很多人第一次接触HCP-MMP1.0的时候以为下载到的就是一个包含360个区域的nii文件。实际上HCP官方发布的时候是给了一套拆分文件每个文件对应一个特定的皮层区域组比如按网络归属、按叶区、按左右半球拆开。这种拆分设计的初衷是方便神经科学研究者按自己的分层假设去使用但对做计算的人来说反而多了一步你得先搞清楚你要的区域在哪个拆分文件里然后还得把它单独提取出来。这套脚本针对的就是这个问题。它把整个HCP-MMP1.0的拆分文件作为输入不需要你自己去做合并或者重新编号。你要做的只是告诉脚本我要哪个区域用名字或者编号都行。脚本内部会把拆分文件读进来定位区域ID生成一个只有0和1的mask体数据然后写成标准的NIfTI文件。整个过程是自动化、可重复的这对做批量处理的人来说很关键——如果你要在预处理好的100个被试上跑同样的mask分析手动操作100次是不现实的但脚本可以保证100次生成的mask像素级一致。2.2 rsamatlab提供了什么脚本又补了什么rsamatlab本身是一个用于静息态fMRI数据分析的MATLAB工具箱功能上覆盖了从数据预处理到功能连接分析的常见环节。它内部已经有一套基于ROI的mask操作工具包括读取NIfTI、将atlas映射到个体空间等。但这个工具箱原本的操作方式偏交互式需要你在图形界面上选来选去不太适合脚本化调用。而且它自带的mask提取逻辑对HCP-MMP1.0的拆分格式支持得不够直接——你想用某个特定区域时还得先自己从拆分文件里找到那个区域在哪个矩阵位置。这套代码的价值就在这个空白点上。它不是从零造了一套新的mask工具而是在rsamatlab已经封装好的读写函数之上加了一个聚焦于HCP-MMP1.0的提取逻辑层。整体上项目是“系统开源”的思路底层依赖rsamatlab的成熟实现上层提供自己的脚本接口。实际使用的时候你会发现这个脚本的输出结果可以直接喂给rsamatlab后续的ROI时间序列提取步骤不需要额外的格式转换。2.3 拿到资源后的目录结构我在本地解压之后梳理了一遍结构核心内容如下rsamatlab-code-fmri-extract-HCP-mask/ ├── extract_hcp_mask.m # 主脚本 ├── hcp_mmp1_0/ # HCP-MMP1.0 拆分文件目录 │ ├── MMP_mpm_left.txt │ ├── MMP_mpm_right.txt │ └── ... # 其他拆分描述文件 ├── utils/ │ ├── load_atlas_labels.m │ └── write_binary_mask.m ├── examples/ │ └── demo_extract.m # 示例调用 └── README.md主要的逻辑都在extract_hcp_mask.m这个文件里utils里的两个函数分别负责读取标签信息和写出mask文件。examples/demo_extract.m是给你抄作业用的示例建议先跑这个确认环境没毛病再改参数动真格。3. 核心代码拆解从一个区域名到一张二值mask的完整链路3.1 主入口参数怎么定义整个脚本的入口函数是extract_hcp_mask我用下来觉得它的参数设计还是比较克制的没有一股脑堆十几个可选项。核心就是四个参数atlas拆分的根目录、目标区域的名称或ID、输出路径以及一个可选的左右半球开关。function mask extract_hcp_mask(atlas_dir, region_id, output_path, varargin) % EXTRACT_HCP_MASK 从HCP-MMP1.0拆分中提取指定区域的二值mask % 输入 % atlas_dir - 拆分文件所在目录包含左右半球的label描述文件 % region_id - 目标区域名称如 V1或数值编号如 9 % output_path - 输出nii文件的保存路径 % 可选参数 % hemisphere - both(默认) | left | right % reference, img_ref - 指定参考图像的空间坐标默认使用atlas自带网格 p inputParser; addRequired(p, atlas_dir, isfolder); addRequired(p, region_id, (x) ischar(x) || isnumeric(x)); addRequired(p, output_path, ischar); addOptional(p, hemisphere, both, (x) ismember(x, {both,left,right})); addParameter(p, reference, [], (x) ischar(x) || isa(x, nifti)); parse(p, atlas_dir, region_id, output_path, varargin{:});inputParser是MATLAB官方的参数解析工具用在这类工具型脚本里比较合适。region_id同时接受字符串和数字这意味着你可以用V1这种人类友好的名字也可以用9这种图谱原始编号。这两个在HCP-MMP1.0里是等价的脚本内部会做查找映射。3.2 读取atlas从拆分文件到logical矩阵HCP-MMP1.0的发布格式中核心的体数据是一个3D的NIfTI文件每个体素的值对应该位置所属区域的ID。但这个文件通常不是单独发布的它和描述文件是配套的。读取逻辑如下function atlas_data load_atlas_volume(atlas_dir) % 读取HCP-MMP1.0的体数据文件 % 拆分文件里volume和label描述文件放在同一目录下 atlas_file fullfile(atlas_dir, MMP_mpm.nii); if ~exist(atlas_file, file) % 某些版本会用 .nii.gz 压缩格式 [~, name] fileparts(atlas_file); atlas_file fullfile(atlas_dir, [name .nii.gz]); end % 使用rsamatlab的工具箱读取NIfTI atlas_vol rsamatlab_load_nii(atlas_file); atlas_data atlas_vol.data; fprintf([INFO] Atlas volume size: %dx%dx%d\n, size(atlas_data,1), ... size(atlas_data,2), size(atlas_data,3));这里有个细节值得注意rsamatlab_load_nii返回的是一个带元数据的结构体里面有.data字段、.voxel_size字段、.transform_matrix字段。后面生成mask时要恢复到同一个空间这几个字段都派得上用场。如果你只用niftiread裸读数据再写出来大概率会在空间坐标上翻车——这个问题我在第5章避坑部分会细说。3.3 区域ID的查找映射字符串名称和数字编号的转换这一小段是整个脚本的逻辑核心。HCP-MMP1.0的标签文件里每个区域有它的编号、简称和全称。脚本要做的事情就是把用户传进来的region_id通过这个映射表换算成体数据里的实际整数值因为体数据里存的就是这个整数。function target_value resolve_region_id(atlas_dir, region_id) % 将区域名称或编号映射为atlas体数据中的实际整数值 label_file fullfile(atlas_dir, MMP_mpm_left.txt); fid fopen(label_file, r); raw textscan(fid, %d %s %s, Delimiter, \t); fclose(fid); ids raw{1}; short_names raw{2}; long_names raw{3}; if isnumeric(region_id) % 数字编号直接匹配注意左半球和右半球的偏移量 if region_id 180 region_id region_id - 180; end match_idx find(ids region_id); if isempty(match_idx) error(Region ID %d not found in atlas labels., region_id); end target_value region_id; else % 字符串名称先匹配简称再匹配全称 match_idx find(strcmpi(short_names, region_id)); if isempty(match_idx) match_idx find(strcmpi(long_names, region_id)); end if isempty(match_idx) error(Region name %s not found in atlas labels., region_id); end % 因为名称可能同名这里统一取第一个匹配 target_value ids(match_idx(1)); end这段代码值得留意的地方在半球偏移量的处理。HCP-MMP1.0的原始数据中右半球的区域编号是在左半球基础上加了一个偏移量的这个偏移量通常是180。如果你不处理这个偏移你输入的V1在右半球对应的是另一个值提取出来的mask就会莫名其妙地多出一块或者少一块。这里用了一个简单的判断当你传入的数字大于180时自动减掉偏移量再去做匹配这样左右半球就能对上。3.4 写出mask保留空间配准信息的关键一步提取到目标区域的整数值之后下一步是在体数据上做逻辑比较。这一步看起来没什么技术含量但操作上有一个常见的坑直接做比较得到的是一个logical数组而NIfTI写出时如果直接用这个logical数组会遇到数据类型写不进去的问题。function write_binary_mask(atlas_data, target_value, output_path, info) % 将目标区域体素置1其余置0写出为nii文件 mask zeros(size(atlas_data), uint8); mask(atlas_data target_value) 1; % 使用rsamatlab的写文件函数沿用原始atlas的空间信息 mask_nii rsamatlab_make_nifti(mask, info.voxel_size, info.transform_matrix); rsamatlab_save_nii(mask_nii, output_path); fprintf([INFO] Mask saved to: %s\n, output_path); fprintf([INFO] Non-zero voxels: %d\n, sum(mask(:))); end这里rsamatlab_make_nifti的作用是把裸的数组重新包装成带空间属性的NIfTI结构。我在第一次用的时候踩过坑直接调MATLAB的niftiwrite写出来的文件虽然在数值上没错但在某些软件比如FSL、AFNI里打开后空间位置对不上就是因为没有继承原始文件的sform和qform信息。rsamatlab这套函数内部会把这层信息带过去所以尽量别自己另造一套写文件的流程。4. 把脚本跑起来环境准备、调用方式与参数对照4.1 环境准备这个脚本依赖rsamatlab工具箱的底层函数所以第一步是把rsamatlab装好。我用的环境是MATLAB R2022b因为脚本里用了arguments块就会要求比较新的版本但这套代码用的还是传统的inputParser所以在R2018b以上的版本应该都能跑。你需要确认几件事MATLAB能找到rsamatlab的路径建议在启动脚本里用addpath把rsamatlab的根目录和本项目目录都加进去。需要安装Image Processing Toolbox因为NIfTI读写底层会用到一些矩阵运算和插值函数。signal toolbox不是必需的。确保HCP-MMP1.0的拆分文件完整。至少要有MMP_mpm.nii或.nii.gz和两个label描述文件。配置路径的推荐写法如下% startup.m 或直接在命令行里执行 addpath(/path/to/rsamatlab); addpath(/path/to/rsamatlab-code-fmri-extract-HCP-mask); addpath(genpath(/path/to/rsamatlab-code-fmri-extract-HCP-mask/utils));addpath(genpath(...))会把子目录也递归加进去这样utils里的函数在任意位置都能被找到。4.2 调用方式命令行还是交互式我个人的习惯是直接命令行调用不打开任何GUI。因为mask提取通常发生在批量预处理流程里交互式操作会打断整个pipeline。直接在MATLAB命令行窗口执行mask extract_hcp_mask(./hcp_mmp1_0/, V1, ./output/V1_left_mask.nii, ... hemisphere, left);如果要对批量区域做提取你可以把区域名放在一个cell数组里循环调用脚本本身没做批量封装但也没有外部状态所以循环是安全的。4.3 参数说明对照参数类型必填说明atlas_dir字符串是拆分文件所在目录region_id字符串或数字是区域简称如V1或编号如9output_path字符串是输出nii文件路径需包含文件名和后缀hemisphere字符串否both/left/right默认bothreference字符串或nifti否参考图像指定后mask会重采样到参考空间4.4 验证输出跑通之后建议做两个快速验证。先在MATLAB里检查非零体素的数量再看一下mask是不是真的落在了期望的脑区。第二个验证我一般用FSL的fsleyes打开把mask叠加到T1图像上看位置。这里有一个可以抄作业的检查脚本% 验证脚本检查mask与原始atlas的空间一致性 new_mask rsamatlab_load_nii(./output/V1_left_mask.nii); orig_atlas rsamatlab_load_nii(./hcp_mmp1_0/MMP_mpm.nii); % 检查体素大小 assert(isequal(new_mask.voxel_size, orig_atlas.voxel_size), Voxel size mismatch!); % 检查变换矩阵 assert(norm(new_mask.transform_matrix - orig_atlas.transform_matrix) 1e-6, ... Transform matrix mismatch!);这两个断言如果通过了说明空间配准信息没有丢。剩下的就是视觉检查了。5. 避坑指南处理HCP图谱时最常见的5个坑5.1 坑左右半球区域搞混mask出现在错误半球现象我用V1提取mask出来的结果却出现在右半球的对应位置左半球是空的。原因HCP-MMP1.0中左右半球使用同一套区域名称但体数据中左半球的值和右半球的值不同。如果你没有通过hemisphere参数指定脚本默认both即左右都要。但在某些区域的实现中左半球和右半球的编号并不是简单的镜像关系有一些区域只存在于一个半球。解决在调用时显式指定hemisphere, left或right。如果结果还是不对先到label描述文件里确认你输入的区域名是不是真的存在于指定的半球。5.2 坑mask和原始图像空间位置对不上现象用FSL打开mask发现位置整体偏移了有的甚至方向翻转。原因没有沿用原始atlas的transform_matrix。你要是用niftiwrite裸写数据这个问题几乎必然出现因为新文件默认的origin在左上角而原文件可能在别的方位。解决写mask时强制沿用原文件的voxel_size和transform_matrix也就是脚本里write_binary_mask函数做的那样。如果你用自己的代码重写了这部分逻辑务必检查输出文件的qform值。5.3 坑用直接比较浮点数组现象mask里出现很多孤立的零散体素或者某些明明应该被包含的体素没有被置1。原因HCP-MMP1.0体数据的存储类型是整数但在某些转换环节里被读成了浮点。浮点数的直接等值比较会因为精度问题导致个别体素匹配失败。解决读取后立即cast成整数类型再比较。如果已经是浮点了比较时加上容差比如abs(atlas_data - target_value) 0.5。这套脚本里没有做成容差模式因为它假设原始数据是整数存储的如果你用的拆分文件经过了中间转换建议在调用前自己做一个round。5.4 坑路径包含中文字符导致读取失败现象MATLAB报错无法打开文件或者打开后数据是乱的。原因部分版本MATLAB在Windows系统下对中文路径的编码处理有问题尤其是涉及textscan读取label文件的时候。解决项目路径和数据路径全部用英文这是最省事的做法。如果实在绕不开中文路径在读取文件前把Encoding参数设为UTF-8但这会引入新的变量不如直接改路径省心。5.5 坑区域名称查找失败报“not found”现象输入V1能匹配输入Primary Visual Cortex就报错。原因HCP-MMP1.0的label文件里简称short name和全称long name的书写规则不完全一致有些区域的全称里包含逗号和括号直接匹配字符串容易失败。解决先到label文件里grep一下实际的全称写法再作为输入。或者干脆只用简称和数字编号这两个是可靠映射。6. 二次开发把“提取一个mask”升级成“批量提取一组mask并重建时间序列”我觉得这个脚本最大的可扩展点在批量提取。fMRI分析里常见的一个需求是把某个网络的所有节点比如默认模式网络的核心区域一次性提出来然后提取每个区域的时间序列做功能连接矩阵。基于extract_hcp_mask的话循环调用就行% 批量提取多个区域的mask并保存到一个矩阵里 regions {PCC, MPFC, L_Angular, R_Angular, Hippocampus}; for i 1:length(regions) mask_file fullfile(./output/masks/, [regions{i} .nii]); extract_hcp_mask(./hcp_mmp1_0/, regions{i}, mask_file, ... hemisphere, both); % 加载对应mask提取roi时间序列依赖rsamatlab的函数 mask_nii rsamatlab_load_nii(mask_file); roi_ts rsamatlab_extract_roi_timeseries(func_data, mask_nii); ts_data(:, i) roi_ts; % 保存到矩阵里 end这个扩展思路很简单但有一个值得注意的点不同区域的mask大小不同体素数量不同。如果你直接对每个mask都做均值提取小区域的时间序列噪声会明显偏大。这种场景下的常见做法是对mask做一次形态学腐蚀操作只保留区域内部的体素去掉边界体素减少部分容积效应对时间序列的影响。代码实现可以这样% 腐蚀mask以减小部分容积效应 se strel(sphere, 1); mask_eroded imerode(mask_vol 0, se);但做这一步需要保证你的mask是严格二值的并且没有孔洞。还有一种思路是保留所有体素但用鲁棒均值trimmed mean代替普通均值。我会更推荐后者因为它不需要对mask做形态学操作实现简单逻辑更稳。对于想更深一步改造的人来说可以考虑把整条链路串成一个独立函数输入原始fMRI数据、HCP图谱根目录、一组区域名输出一个被试的区域间功能连接矩阵。这其实就是很多脑影像分析中“图谱到连接矩阵”的标准化流程。有了extract_hcp_mask这个基础函数这一步只是把循环提出来而已没有新的技术难点。我自己的使用习惯是所有mask提取都强制走脚本生成绝不在FSL或MATLAB里手动画ROI。手动画的mask不仅不可重复还会因为你哪天手一抖多点了一层体素导致结果完全不同。从那以后我在任何项目里需要用到图谱mask时都会先确认有没有脚本可以做没有再写一个这已经成了我做脑影像分析的一条铁律。希望这个项目也能在你这儿派上用场。本文还有配套的精品资源点击获取
