简介：面向医学图像处理初学者与需在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/MAT：Matlab 读取链路的三个决策点

手里有一批 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{end+1} = r; %#ok<AGROW> 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 end

sort按投影坐标升序，比按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); end

findgroups和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和RescaleIntercept；MR 或定量图建议用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)); end

outClass是第一个要决定的参数：int16省内存，适合 CT；single兼容性更好，适合 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 RAS：x、y 取反，z 不变 A = diag([-1, -1, 1, 1]) * A_lps; end

这段A_lps的列顺序容易被写反。记住一点：Matlab 矩阵第一维是行下标，行下标增加的方向是 DICOM 的第二组方向余弦c；第二维是列下标，列下标增加的方向是第一组方向余弦r。ps(1)是行间距，配c；ps(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 保存 MAT：matfile 与 -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 保存 NPZ：Matlab 没有原生 npz，走 Python 桥或手写 zip

NPZ 本质是一个 zip 包，里面放若干.npy文件。Matlab 没有原生savez，常见做法有两个：配置pyenv后直接调用 Python 的 numpy，或者先

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