rsamatlab提取HCP-MMP1.0脑区mask：从区域名到二值掩码

简介：针对脑影像分析中需要从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 参数说明对照

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时，都会先确认有没有脚本可以做，没有再写一个，这已经成了我做脑影像分析的一条铁律。希望这个项目也能在你这儿派上用场。

本文还有配套的精品资源，点击获取