Matlab读取DICOM并转NIfTI/NPZ/MAT的完整链路与避坑指南

简介面向医学图像处理初学者与需在Matlab与Python间迁移数据的科研人员这套代码示例聚焦DCM医学影像的读取与格式转换覆盖从原始DCM到NII、NPZ、MAT三种常用存储格式的流程。包内共9个文件以8个.m函数脚本和1个.txt说明文档为主m文件承担DICOM读取、NIfTI头信息解析与保存等核心功能txt提供使用说明压缩包约20KB便于直接嵌入现有工程调试。已有3181人学习下载说明该流程在医学影像预处理中通用性较高。读者可获得可直接运行的转换脚本与函数调用范例理解dicomread、dicominfo等接口用法并借助load_nii、save_nii等函数完成NIfTI结构读写掌握元数据保留与格式适配思路减少搭建管线时的试错成本。1. 从 DICOM 到 NIfTI/NPZ/MATMatlab 读取链路的三个决策点手里有一批 DICOM 序列要在 Matlab 里读成体数据再分别保存为 NII、NPZ 或 MAT难点通常不在dicomread那一行而在序列分组、层序、方向、仿射和数据类型。这份资源解决的是完整链路Matlab 读取 DCM 文件按SeriesInstanceUID重建三维体数据保存为 NIfTI 供影像工具查看或保存为 NPZ/MAT 供训练脚本和数值分析使用。适合做医学影像预处理、算法验证、跨语言数据交换的从业者也适合刚接手一批 DICOM 目录、需要先跑通最小闭环的新手。真正会翻车的地方往往是把“能读出来”当成“读对了”。2. 目录扫描与 DICOM 元数据解析先把序列、层厚和方向搞对DICOM 目录很少是理想状态一个文件夹里可能混着定位像、增强前后序列、不同重建核的薄层和厚层甚至还有导出时留下的重复文件。直接对所有文件循环dicomread再按文件名排序最后得到的体数据大概率层序错乱。我的习惯是先把“像素读取”放一边只读头信息建索引把SeriesInstanceUID、InstanceNumber、ImagePositionPatient、ImageOrientationPatient、PixelSpacing这些字段抓出来再决定怎么堆体数据。2.1 用 dicominfo 建索引从文件名猜序列是下策dicominfo只解析头信息速度比dicomread快很多适合先扫一遍目录。下面这个函数把常见字段拍平成表避免后面反复读头。注意代码里把三维向量拆成IPP1、IPP2这类标量列是为了让struct2table不把向量字段展开成不可控的列。function idx build_dicom_index(rootDir) % 扫描目录只读 DICOM 头信息不读像素 files dir(fullfile(rootDir, **, *.dcm)); if isempty(files) % 有些导出工具不带 .dcm 后缀退化为扫全部文件再试读 files dir(fullfile(rootDir, **, *)); end rows {}; for k 1:numel(files) if files(k).isdir continue; end f fullfile(files(k).folder, files(k).name); try info dicominfo(f); catch continue; % 不是 DICOM、权限不足或文件损坏跳过 end r struct(); r.File string(f); r.SeriesUID string(getfield_or(info, SeriesInstanceUID, )); r.InstanceNumber getfield_or(info, InstanceNumber, NaN); r.Rows getfield_or(info, Rows, NaN); r.Columns getfield_or(info, Columns, NaN); % 位置和方向拆成标量方便后面做表操作 r.IPP1 getvec(info, ImagePositionPatient, 1, NaN); r.IPP2 getvec(info, ImagePositionPatient, 2, NaN); r.IPP3 getvec(info, ImagePositionPatient, 3, NaN); r.IOP1 getvec(info, ImageOrientationPatient, 1, NaN); r.IOP2 getvec(info, ImageOrientationPatient, 2, NaN); r.IOP3 getvec(info, ImageOrientationPatient, 3, NaN); r.IOP4 getvec(info, ImageOrientationPatient, 4, NaN); r.IOP5 getvec(info, ImageOrientationPatient, 5, NaN); r.IOP6 getvec(info, ImageOrientationPatient, 6, NaN); ps getfield_or(info, PixelSpacing, [NaN; NaN]); r.RowSpacing ps(1); % DICOM PixelSpacing 第一个是行间距 r.ColSpacing ps(2); r.SliceThickness getfield_or(info, SliceThickness, NaN); r.SpacingBetweenSlices getfield_or(info, SpacingBetweenSlices, NaN); rows{end1} r; %#okAGROW end if isempty(rows) idx table(); return; end idx struct2table([rows{:}]); end function v getfield_or(s, name, default) if isfield(s, name) ~isempty(s.(name)) v s.(name); else v default; end end function v getvec(s, name, k, default) if isfield(s, name) numel(s.(name)) k v double(s.(name)(k)); else v default; end end这段代码的关键参数是rootDir和递归通配符。dir(fullfile(rootDir, **, *.dcm))在较新的 Matlab 里能递归子目录老版本如果不认**换成genpath逐目录扫或者先在外面用文件清单喂进来。第二步的“扫全部文件”是兜底代价是可能对大量非 DICOM 文件做try/catch目录特别大时建议先按扩展名和文件头过滤。getfield_or和getvec只是防止某些设备缺字段导致整个索引中断缺失值我一般保留NaN后面排序时再显式处理。索引表建好后不要急着取“文件数最多的序列”就开工。常见做法是先看每个SeriesUID的切片数、行列数、ImageType和层间距确认它是原始横断位、增强序列还是重建图像。文件名里的T1、T2、ADC只能当辅助线索真正可靠的还是 UID 和空间信息。2.2 判断体素间距、层间距与缺失切片DICOM 的PixelSpacing是[rowSpacing, colSpacing]也就是行方向间距和列方向间距不是简单的[dx, dy]。层间距则有两个常见来源SliceThickness和SpacingBetweenSlices。前者是层厚后者才是相邻层中心之间的实际间距。增强薄层重建时两者可能不相等如果代码里默认dz SliceThickness体素尺寸就会偏。更稳妥的做法是用ImagePositionPatient在层面法向上的投影差来算实际层间距。法向量由ImageOrientationPatient的两组方向余弦叉乘得到function [sub, dz] sort_series_by_position(sub) % 按层面法向投影排序并估计实际层间距 iop1 [sub.IOP1, sub.IOP2, sub.IOP3]; % 行方向余弦 iop2 [sub.IOP4, sub.IOP5, sub.IOP6]; % 列方向余弦 n cross(iop1, iop2, 2); % 每层的法向量 pos [sub.IPP1, sub.IPP2, sub.IPP3]; d sum(pos .* n, 2); % 投影到法向 [~, ord] sort(d); sub sub(ord, :); d d(ord); dz median(diff(d)); if isempty(dz) || isnan(dz) dz median(sub.SpacingBetweenSlices, omitnan); end if isnan(dz) dz median(sub.SliceThickness, omitnan); end % 间距突变通常意味着缺层、重复层或混入了其他序列 if numel(d) 2 dev abs(diff(d) - dz); if max(dev) 0.2 * abs(dz) warning(层间距不一致最大偏差 %.3f检查是否缺层或混层, max(dev)); end end endsort按投影坐标升序比按InstanceNumber更贴近空间真实顺序。InstanceNumber在部分导出工具里会被重排或者不同序列共用编号单靠它排序很容易出现“中间某一层插到最前面”的情况。median(diff(d))对少量异常层不敏感但如果偏差超过 20%就值得停下来看可能是缺层也可能是同一序列里混了不同层厚。warning不是摆设批量处理时把它记进日志比后面看到体数据被拉伸再回头查要省时间。2.3 按 SeriesInstanceUID 分组并排序一个目录里多序列是常态。分组时以SeriesUID为主键不要用文件夹名硬编码。分组后先看每个序列的切片数是否和Rows、Columns一致再决定用哪一组。function [sub, dz] pick_series(idx, targetUID) % 按 SeriesUID 取子表targetUID 为空时取切片数最多的序列 if nargin 2 || strlength(string(targetUID)) 0 [g, uid] findgroups(idx.SeriesUID); counts splitapply(numel, idx.SeriesUID, g); [~, imax] max(counts); targetUID uid(imax); end sub idx(idx.SeriesUID string(targetUID), :); if isempty(sub) error(没有找到指定 SeriesUID 的切片); end % 同一序列内行列数必须一致 if numel(unique(sub.Rows)) 1 || numel(unique(sub.Columns)) 1 error(同一序列内行列数不一致先检查是否混入定位像或重建图); end [sub, dz] sort_series_by_position(sub); endfindgroups和splitapply适合快速统计targetUID为空时取切片数最多的序列只是启动策略不是诊断结论。真实项目里我一般会把候选序列的 UID、切片数、层间距、行列数打印成一张小表人工确认一次再把确认过的 UID 写进配置。这样下次换数据时至少知道该改哪一行而不是重新猜。unique(sub.Rows)检查能挡住一类常见错误同一 UID 下混入了不同矩阵尺寸的图像堆体数据时会在vol(:,:,k)处报维度不匹配。3. 体数据重建与 NIfTI 写出niftiwrite 的 affine 与数据类型体数据重建的核心不是“把切片塞进一个三维数组”而是让数组下标和病人空间坐标一一对应。Matlab 的niftiwrite能把三维数组写成 NIfTI但头信息里的Affine如果不自己构造查看器只能按默认方向猜。NIfTI 用 RAS 坐标DICOM 的ImagePositionPatient和方向余弦是 LPS 坐标中间差两个符号。这个转换不处理左右翻转就会像玄学一样出现。3.1 把像素矩阵堆成三维体数据先确定输出数据类型。CT 的 HU 值范围通常能塞进int16但前提是先应用RescaleSlope和RescaleInterceptMR 或定量图建议用single避免小数被截断。下面函数按排序后的子表逐层读取预分配三维数组读一层写一层。function [vol, meta] load_series(sub, outClass) % sub 必须已经按空间位置排序outClass 常用 int16 或 single n height(sub); rows sub.Rows(1); cols sub.Columns(1); vol zeros(rows, cols, n, outClass); firstInfo dicominfo(char(sub.File(1))); for k 1:n info dicominfo(char(sub.File(k))); raw dicomread(info); % 单帧 DICOM 正常是二维出现三维/四维时先检查是不是多帧 if ndims(raw) 2 raw squeeze(raw); if ndims(raw) ~ 2 error(第 %d 个文件是多帧 DICOM需要按帧拆分后再堆体数据, k); end end raw double(raw); slope getfield_or(info, RescaleSlope, 1); inter getfield_or(info, RescaleIntercept, 0); val raw * slope inter; if contains(outClass, int) lo double(intmin(outClass)); hi double(intmax(outClass)); if any(val(:) lo) || any(val(:) hi) warning(第 %d 层像素值超出 %s 范围可能被截断, k, outClass); end end vol(:, :, k) cast(val, outClass); end % 用第一层和实际层间距构造仿射 dz median(diff([sub.IPP1, sub.IPP2, sub.IPP3] * ... cross([sub.IOP1(1) sub.IOP2(1) sub.IOP3(1)], ... [sub.IOP4(1) sub.IOP5(1) sub.IOP6(1)]))); if isnan(dz) dz median(sub.SpacingBetweenSlices, omitnan); end meta struct(); meta.Affine dicom_affine(firstInfo, abs(dz)); meta.VoxelSize [sub.RowSpacing(1), sub.ColSpacing(1), abs(dz)]; meta.SeriesUID char(sub.SeriesUID(1)); meta.Source char(sub.File(1)); endoutClass是第一个要决定的参数int16省内存适合 CTsingle兼容性更好适合 MR 和后续归一化。dicomread(info)返回的是原始像素矩阵RescaleSlope和RescaleIntercept要自己乘加CT 的-1024空气值才会正确出现。cast之前先做范围检查是因为整数类型转换会截断等 NIfTI 写出去再发现窗宽窗位不对排查成本更高。多帧 DICOM 这里直接报错是故意的多帧文件的帧间距、帧起始位置可能和单帧序列不同硬塞进vol(:,:,k)只会得到一个看似成功、实际错位的体数据。3.2 用 niftiwrite 保存 NII头信息如何对齐NIfTI 的Affine是 4×4 矩阵前 3 列是体素三个方向在 RAS 空间中的向量第 4 列是原点。DICOM 的ImageOrientationPatient给的是行、列方向余弦ImagePositionPatient给的是第一层第一像素的位置。构造时要把行列顺序和 Matlab 数组的row, col, slice对齐。function A dicom_affine(info, dz) % 输入单层 DICOM 头信息输出 NIfTI 用的 4x4 RAS affine iop double(info.ImageOrientationPatient(:)); r iop(1:3); % 行方向余弦对应列下标增加方向 c iop(4:6); % 列方向余弦对应行下标增加方向 n cross(r, c); % 层面法向DICOM LPS ps double(info.PixelSpacing(:)); % [rowSpacing; colSpacing] ipp double(info.ImagePositionPatient(:)); % vol(row, col, slice) 的下标顺序 % 第 1 维 row - c * rowSpacing % 第 2 维 col - r * colSpacing % 第 3 维 slice - n * dz A_lps [c * ps(1), r * ps(2), n * dz, ipp]; % DICOM LPS 转 NIfTI RASx、y 取反z 不变 A diag([-1, -1, 1, 1]) * A_lps; end这段A_lps的列顺序容易被写反。记住一点Matlab 矩阵第一维是行下标行下标增加的方向是 DICOM 的第二组方向余弦c第二维是列下标列下标增加的方向是第一组方向余弦r。ps(1)是行间距配cps(2)是列间距配r。n * dz的符号由叉乘决定dz取正排序时已经按投影升序所以第三列方向与层序一致。LPS 转 RAS 只改前两行符号原点也跟着变不要把ipp单独取反后拼回去。写出时先写一次拿到niftiinfo模板再替换Affine重写。这样比手工拼 NIfTI 结构稳也能让niftiwrite处理压缩和头字段。niiPath fullfile(outDir, series.nii.gz); niftiwrite(vol, char(niiPath), Compressed, true); nii niftiinfo(char(niiPath)); nii.Affine meta.Affine; nii.Description DICOM series converted in Matlab; niftiwrite(vol, char(niiPath), nii, Compressed, true);Compressed为true时输出.nii.gz路径后缀最好也跟着写.nii.gz不要用.nii骗自己。nii.Affine必须是 4×4 双精度矩阵如果从别处拿来的是 3×4先补[0 0 0 1]。Description只是备注不参与几何定位。部分 Matlab 版本对niftiwrite的info参数和Compressed组合支持不一致如果第二次写出报参数错误就退回先写未压缩.nii再用外部工具压缩别在版本兼容上耗太久。3.3 验证 NII 读回与方向一致性写完必须读回验证不要只看文件大小。至少检查三件事数组尺寸是否一致、Affine是否非奇异、读回后的数据类型是否符合预期。V niftiread(char(niiPath)); infoBack niftiinfo(char(niiPath)); assert(isequal(size(V), size(vol)), NIfTI 尺寸和源体数据不一致); assert(abs(det(infoBack.Affine(1:3, 1:3))) 0, Affine 奇异方向矩阵有问题); assert(strcmp(class(V), class(vol)), 读回数据类型发生变化); fprintf(NIfTI 尺寸: %s\n, mat2str(size(V))); fprintf(体素尺寸: %.3f %.3f %.3f\n, infoBack.PixelDimensions(1:3));det检查只能说明矩阵可逆不能证明方向正确方向还要靠查看器或已知标记点确认。常见做法是挑一个明显不对称的层面比如带肝脏或带病灶的一侧和 DICOM 查看器对比左右。niftiread读回的数据类型可能被 NIfTI 数据类型字段影响如果原来是int16读回变成double先看nii.Datatype是否被模板覆盖。PixelDimensions和Affine的体素尺寸应当一致不一致时以Affine为准因为查看器定位靠的是它。4. NPZ 与 MAT 的取舍压缩、字段结构和跨语言读取NII 适合给影像查看器和分割/配准工具MAT 适合 Matlab 内部继续处理NPZ 适合交给 Python 训练脚本。三种格式没有绝对优劣关键是字段命名、数据类型和压缩方式要提前约定。批量转换时最怕“Matlab 里读得出来Python 里读出来是 object 数组”那通常是保存方式或字段结构埋的雷。4.1 保存 MATmatfile 与 -v7.3 的边界MAT 是 Matlab 原生格式保存体数据和元信息最直接。小文件用默认save就够超过 2 GB 或需要 HDF5 兼容时用-v7.3。但-v7.3不是万能钥匙老版本 Matlab 打不开Python 侧scipy.io.loadmat也读不了需要改用 HDF5 读取方式。matPath fullfile(outDir, series.mat); % 小数据、需要给旧版 Matlab 或 scipy 读取时用 -v7 % 大数据、单变量超过 2GB 时用 -v7.3 save(char(matPath), vol, meta, -v7.3); % 需要追加变量或分块写入时用 matfile m matfile(char(matPath), Writable, true); m.vol vol; m.meta meta;-v7.3的代价是写入慢、文件结构变成 HDF5部分工具链不认。我的习惯是如果下游是 Python 且数据量不大优先存-v7如果只是 Matlab 内部用且体数据很大才上-v7.3。matfile适合“先建文件再分块写”的场景但不要把它当成无限内存的替代品写入过程中仍然要注意变量类型和维度顺序。meta里放Affine、VoxelSize、SeriesUID、Source这些字段就够了不要把完整 DICOM 头原样塞进去字段太多会让下游读取变慢。4.2 保存 NPZMatlab 没有原生 npz走 Python 桥或手写 zipNPZ 本质是一个 zip 包里面放若干.npy文件。Matlab 没有原生savez常见做法有两个配置pyenv后直接调用 Python 的 numpy或者先本文还有配套的精品资源点击获取