简介：这是一套面向医学影像分析学习者与科研人员的MATLAB图像配准代码包，重点解决不同时间、不同设备或不同成像方式获取图像的对齐、比较与融合问题，尤其适合研究非刚性配准及生物组织形变适应。压缩包共34个文件，约240KB，以22个.m脚本和6个.c源码为核心，辅以4张png示例图、1个fig界面文件及说明文本。其中.m文件承担主流程与示例演示，.c源码负责计算密集的变换与相似度度量，示例图与界面文件用于直观核对配准效果。内容覆盖二维与三维B样条变换、刚体变换预处理、互信息相似度度量、基于梯度的优化求解以及GUI可视化，并提供多组配准示例，便于读者理解从图像预处理、模型选择到参数优化的完整流程。已有1336人学习，适合具备MATLAB基础、希望深入掌握经典配准算法实现细节并用于课程实验或课题复现的读者参考。

1. 两张片子叠不上的时候，医学图像配准这套 MATLAB 工程能干什么

早上查房前，某影像科要把同一患者的术前 CT 和术中 MR 叠在一张图上，滑块拖了半天：脑室对上了，颅骨又偏两个像素。医学图像配准做的就是这件事——在 MATLAB 里求一个空间变换，把浮动图映射到参考图的坐标系，让两幅图在解剖结构上严丝合缝。这份工程把互信息测度、仿射与非刚体变换、多分辨率优化和一套能点着看的 GUI 串成完整链路，适合两类人：手里有 CT/MR/US 数据要做术前术后对比的研究生，以及想把配准当成预处理步骤塞进自己分割、重建流水线的工程师。它不追求跑分，胜在每一步都能打开函数看清楚。

2. 工程结构拆解：从 main 入口到 imregconfig 的六个目录各管什么

拿到一份 MATLAB 工程，我第一件事不是双击 main.m，而是先看目录树。配准类代码有个通病：算法本身不到一百行，剩下的几百行全在读图、转格式、坐标换算、画图上，结构不清楚的话，改一个参数能被调用关系绕晕。这份工程的目录划分走的是典型的三层路子——入口层、算法层、工具层，配置单独抽出来，数据和代码分开。

2.1 目录树与文件分工

常见做法是把代码按职责切成下面这个样子，实际下载下来即使文件命名有出入，职责划分基本跑不出这个框架：

medical_registration/ ├── main.m % 单组配准入口，改两行路径就能跑 ├── main_batch.m % 批量遍历 data/ 下的病例目录 ├── gui_registration.m % GUI 回调逻辑，与同名 .fig 成对出现 ├── gui_registration.fig ├── config/ │ └── cfg_default.m % 金字塔层数、bins、优化器参数集中在这里 ├── core/ │ ├── mi_measure.m % 互信息相似性测度 │ ├── ncc_measure.m % 归一化互相关，单模态场景用 │ ├── objective.m % 参数向量 -> tform -> 测度值 │ ├── pyramid_reg.m % 由粗到精的多分辨率主循环 │ └── demons_reg.m % 非刚体位移场估计 ├── utils/ │ ├── load_pair.m % 读图、转灰度、归一化到 [0,1] │ ├── dice_score.m % 分割掩膜重叠度评估 │ └── checkerboard.m % 棋盘格叠加，肉眼验收用 └── data/ ├── case01/ └── case02/

入口层三个文件对应三种使用姿势：命令行单跑、批量跑、界面交互跑。算法层全部是无状态函数（core/下不读文件、不画图、不弹窗），输入图像矩阵和参数，输出变换对象或位移场。这么做的好处是单元测试好写，我在自己的改造里给mi_measure单独喂一对人工平移过的图像，就能验证测度是不是在零位移处取极值。工具层是纯搬砖的，load_pair负责把 uint8、uint16、三通道、不同尺寸这些脏活挡在算法之外，checkerboard只管可视化。

config/cfg_default.m是个返回结构体的函数，不是脚本。这点很关键——脚本形式的配置会污染 base 工作区，批量跑的时候外层循环的变量被内层脚本覆盖，是血泪经验级别的坑。函数形式每次调用返回一份干净副本，cfg = cfg_default();之后想改哪个字段改哪个字段。

2.2 一条配准链路的三段式骨架

抛开 GUI 和文件 IO，整个配准过程就是三件事：给一个空间变换，把浮动图搬过去，用某个指标衡量搬完之后像不像，然后让优化器去搜让指标最好的变换。

第一段是空间变换。工程里默认用affine2d，一个 3×3 矩阵同时表达平移、旋转、缩放和剪切，能覆盖大部分刚体加近似缩放的场景。需要更强的形变能力时切到imregdemons输出的稠密位移场。

第二段是图像重采样。imwarp拿到变换后，内部按逆映射把参考图网格上的每个像素点映射回浮动图坐标，再插值取值。这里有个方向问题必须先想清楚：affine2d定义的是浮动图坐标到参考图坐标的前向映射，imwarp自己负责求逆。写反了的话，配准照样能跑完，只是结果永远差一截，属于最难查的那类 bug。

第三段是优化。目标函数是「参数向量 → 变换矩阵 → 重采样 → 测度值」，优化器在里面反复试参数。因为测度对参数不一定可导（互信息涉及直方图统计，是分段常值函数），工程里走的是无导数优化路线。

2.3 为什么用 MATLAB，以及它撑不住的场景

选 MATLAB 做医学图像配准，核心理由是 Image Processing Toolbox 把链路铺好了：imregconfig直接给出单模态和多模态两套优化器与测度预设，imregtform一行就能跑出仿射变换，imregdemons把非刚体光流做成了黑盒，imref2d管理世界坐标系省掉大量手写坐标换算。自己从零写互信息加 Powell，最少两天；调通 MATLAB 的imregister，半小时。

它的边界也很明确。一是内存，三维体积数据（比如 512×512×300 的 CT）转 double 就是 600MB 起，imwarp重采样再翻一倍，普通笔记本跑两个病例就开始交换内存。二是速度，无导数优化每次迭代都要做一次完整重采样，一次配准几百次迭代很正常，四维或者大体积数据基本得靠降采样或者转 C 混合编译。三是部署，想让临床科室直接双击就用，MATLAB Runtime 打出来的包动辄几百兆，还挑版本。所以我的判断是：做算法验证、写论文、算指标，MATLAB 是效率最高的选择；做成日处理上百例的产线服务，最后还是要把核心循环用 MEX 或者别的语言重写一遍。这份工程的价值在前者，它把每一步都拆开摆着，方便你对照着移植。

3. 核心算法落到代码：互信息测度、仿射矩阵与优化器怎么接

这一章是我拆这类工程时花时间最多的地方。配准结果不准，八成不是优化器不行，而是测度写错了、坐标搞反了、直方图被补零区域污染了。下面按调用顺序把核心几个函数拆开，代码是按工程里最常见的写法整理的，你自己对着源码看，逻辑应该能一一对上。

3.1 相似性测度：互信息为什么能扛住跨模态

单模态配准用归一化互相关就够了——同一台机器两次扫描，灰度值基本是线性相关的，NCC 在正确配准位置有尖峰。但 CT 和 MR 之间不存在这种线性关系：骨在 CT 上是高亮，在 MR 上可能一片黑；软组织反过来。互信息不假设灰度之间的具体函数形式，只看两幅图的联合直方图是否足够「集中」——对齐之后，一个灰度值大概率只对应另一个特定的灰度值，联合分布窄；没对齐时联合分布摊开，熵大。互信息就是边缘熵之和减去联合熵。

function mi = mi_measure(fixed, moving, nbins) % 输入：同尺寸 double 图像，已归一化到 [0,1] % 输出：互信息值，越大表示对齐程度越高 fixed = double(fixed(:)); moving = double(moving(:)); % 只在两幅图都非零的重叠区统计，避免 imwarp 补零的边缘拉低测度 valid = fixed > 0 & moving > 0; fixed = fixed(valid); moving = moving(valid); if numel(fixed) < 100 mi = -inf; % 重叠区太小，直接判定为坏解 return; end edges = linspace(0, 1, nbins + 1); joint = histcounts2(moving, fixed, edges, edges) / numel(fixed); pM = sum(joint, 2); % 浮动图边缘分布 pF = sum(joint, 1); % 参考图边缘分布 pMF = pM * pF; % 外积，与 joint 同尺寸 nz = joint > 0 & pMF > 0; % 去掉零概率项，否则 log2(0) 出 NaN mi = sum(joint(nz) .* log2(joint(nz) ./ pMF(nz))); end

nbins一般取 32 到 64，太小了测度曲线很平滑但分辨不出细微偏移，太大了直方图变成稀疏矩阵，测度曲线抖得优化器原地打转。我一般从 32 起步，发现配准精度不够再往上加。

代码里valid那一行是整个函数最容易被忽略也最影响结果的地方。imwarp重采样后，落在参考图视野外的像素会被填成 0，如果把这些 0 也算进联合直方图，那么每变换一次就凭空多出一大块「灰度 0 对应灰度 0」的统计量，测度在平移量变大时会虚高，优化器反而往边界跑。工程里如果没做这个过滤，表现就是浮动图越移越远，最后贴到角上。

3.2 空间变换：3×3 仿射矩阵里的平移、旋转与缩放

affine2d的矩阵是列向量右乘形式，第三行是平移量：

% 参数向量 p = [tx, ty, theta, scale] t = p(3); s = p(4); A = [ s*cos(t) s*sin(t) 0; -s*sin(t) s*cos(t) 0; p(1) p(2) 1 ]; tform = affine2d(A);

前两行是旋转与缩放的组合，第三行是平移。注意affine2d要求最后一行是[0 0 1]，如果手写成[tx ty 1]放在最后一行、把旋转写在下面的行，MATLAB 会直接报矩阵非法。这一点和不少教材里用 OpenCV 的写法正好反过来——OpenCV 的warpAffine前两行是变换、最后一列是平移，两边对照着改代码时百分之百会栽一次，我第一次移植就翻了车。

变换打通之后，重采样这一步写成：

moved = imwarp(moving, tform, ... 'OutputView', imref2d(size(fixed)), ... % 输出网格对齐参考图 'FillValues', 0, ... % 视野外补 0，配合 valid 掩膜 'Interp', 'linear'); % 线性插值，够用且快

OutputView一定要显式给，不写的话imwarp会按变换后图像的外接矩形输出，尺寸和参考图对不上，后面的测度算出来全是错位。Interp在粗配准阶段可以降成'nearest'提速，精配准阶段再换'linear'或者'cubic'。

3.3 优化器与多分辨率金字塔：参数从粗到精怎么传

把上面两段拼成目标函数，剩下的交给优化器：

function loss = objective(p, fixed, moving, cfg) t = p(3); s = p(4); A = [ s*cos(t) s*sin(t) 0; -s*sin(t) s*cos(t) 0; p(1) p(2) 1]; moved = imwarp(moving, affine2d(A), ... 'OutputView', imref2d(size(fixed)), 'FillValues', 0); loss = -mi_measure(fixed, moved, cfg.nbins); % 取负，把最大化转成最小化 end

只跑单尺度会翻车：初始错位几十像素时，互信息曲面在搜索空间里有一堆局部极值，优化器很容易卡在某个「看上去还行」的位置。所以主循环做成金字塔，先降采样把大位移吃掉，再逐层细化：

p = cfg.initParams; % [0 0 0 1] opts = optimset('MaxIter', cfg.maxIter, 'TolX', 1e-4, 'TolFun', 1e-6, ... 'Display', 'off'); for level = cfg.levels:-1:1 % 层数从粗到精，比如 3 -> 2 -> 1 sc = 2^(level - 1); fLow = imresize(fixed, 1/sc, 'Antialiasing', true); mLow = imresize(moving, 1/sc, 'Antialiasing', true); p = fminsearch(@(q) objective(q, fLow, mLow, cfg), p, opts); p(1:2) = p(1:2) * sc; % 平移量还原到下一层尺度 end

imresize必须开Antialiasing，否则降采样产生混叠，互信息曲面全是毛刺。p(1:2)乘回尺度这一行最容易被漏掉：在 4 倍降采样层求出平移 3 像素，实际原图是 12 像素，不还原的话下一层从 3 起步，等于白算了上面一层。我见过有实现干脆每层都从零初值重新搜，那金字塔就退化成了三倍耗时，一点精度收益都没有。

fminsearch用的是 Nelder-Mead 单纯形，不需要梯度，对这个不可导的测度很合适。有 Optimization Toolbox 的话换成patternsearch，可以给每个参数设上下界，把旋转限制在 ±30°、缩放限制在 0.8 到 1.2，防止优化器搜出个把图像翻面的解——那种解互信息值经常也不低，因为你把图像镜像过去之后，灰度分布的对齐关系还在。

3.4 非刚体延伸：Demons 位移场与 B 样条自由形变

仿射只能处理整体平移旋转，软组织形变（呼吸、组织挤压）得靠非刚体。工程里如果带demons_reg.m，大概率是包了一层imregdemons：

% 多分辨率级数 [200 100 50] 表示从粗到细的迭代次数 % AccumulatedFieldSmoothing 越大位移场越平滑，形变越保守 [D, ~] = imregdemons(moving, fixed, [200 100 50], ... 'AccumulatedFieldSmoothing', 1.5, ... 'DisplayWaitbar', false); registered = imwarp(moving, D); % D 是 H×W×2 的位移场

D(:,:,1)是每个像素的行方向位移，D(:,:,2)是列方向位移。AccumulatedFieldSmoothing是最重要的参数，取值 1.0 左右形变很激进，容易把噪声也拟合进去；取 2.0 以上位移场被过度平滑，形变基本退化成刚体。我一般先用 1.0 出一版看趋势，看形变集中在不该动的地方就往上调到 1.5 附近。

更灵活的方案是 B 样条自由形变（FFD），把控制网格铺在图像上，每个控制点的位移作为待优化参数。这是配准论文里最常见的非刚体模型，网格越密能表达的形变越细，参数也越多，优化代价随控制点数立方增长。工程源码里如果没有现成实现，可以用bsplinegrid之类的第三方函数配合ppval自己搭，思路和仿射一致：控制点位移 → 插值出稠密位移场 →imwarp重采样 → 算测度。

4. 跑通第一组结果：从 addpath 到棋盘格叠加图落盘

拆完结构就该动手了。这一章按我第一次运行这类工程的顺序走一遍，中间该注意的地方都标出来，照着做应该能在半小时内看到第一张叠加图。

4.1 环境与路径：三个工具箱和一个容易忽略的编译前提

依赖三个东西：MATLAB 主体、Image Processing Toolbox、Optimization Toolbox（只有用到patternsearch或fminunc时才必须）。检查一下直接敲：

ver % 看工具箱装没装 which imregdemons % 空结果说明 IPT 缺失或版本太老

路径设置不要用鼠标点「添加到路径」然后关窗口就忘，写成脚本放工程根目录，每次开 MATLAB 先跑一遍：

root = fileparts(mfilename('fullpath')); addpath(fullfile(root, 'core')); addpath(fullfile(root, 'utils')); addpath(fullfile(root, 'config'));

fileparts(mfilename('fullpath'))拿到的是当前脚本所在目录，比写绝对路径可靠——换台机器、换个盘符都不用改。用genpath一把梭也行，但会把data/里的子目录都加进去，如果里面有同名的旧版函数，MATLAB 按路径顺序调用的老版本，改代码不生效就是它闹的。

如果工程里带 MEX 文件（有些实现为了提速把互信息循环编译成 C），要先跑mex -setup选编译器。没装编译器的机器上直接Undefined function或者Invalid MEX-file，报错信息不会告诉你缺编译器。

4.2 单模态与多模态两条调用入口

有两条路。快的一条直接用imregconfig：

fixed = im2double(imread('data/case01/fixed.png')); moving = im2double(imread('data/case01/moving.png')); % 关键一步：多模态用 multimodal，单模态用 monomodal [optimizer, metric] = imregconfig('multimodal'); optimizer.InitialRadius = 0.006; % 初始步长，太大容易一步跨过头 optimizer.MaximumIterations = 300; tform = imregtform(moving, fixed, 'affine', optimizer, metric); moved = imwarp(moving, tform, 'OutputView', imref2d(size(fixed)));

imregconfig('multimodal')内部用的是 Mattes 互信息加 1+1 进化策略优化器，参数已经调过一轮，直接跑通常能出结果。InitialRadius是它最容易调坏的地方——默认值 0.01 对大幅错位够用，对精细配准嫌大，会看到输出图在正确位置附近来回抖；改成 0.004 到 0.006 之间比较稳。

慢的一条是自己写的pyramid_reg，可控性高，每一步都能打日志。两条路的差别在于：前者省事、内部是个黑匣子；后者参数全暴露、能加约束、能换测度。我第一次跑建议先用imregtform拿个基准结果，再用自写版本对照，两者结果偏差超过两三个像素，说明自写版本的某一步出了问题。

4.3 结果怎么验：互信息值、Dice 与棋盘格叠加

配准跑完只看一张叠图是不够的，我一般从三个角度验。第一是数值，把配准前后的互信息值打出来，正常情况应该有明显提升；提升不到 5%，要么本来就对齐了，要么算法没起效。第二是定性叠加，棋盘格最直观：

function out = checkerboard(A, B, tile) % A、B 同尺寸，按 tile×tile 的方块交替取像素 [h, w] = size(A); [X, Y] = meshgrid(1:w, 1:h); mask = mod(floor((X-1)/tile) + floor((Y-1)/tile), 2) == 0; out = A; out(mask) = B(mask); end

方块取 32 像素左右，看方块边界上的血管、骨缘有没有错位断线，一眼就能定生死。第三是定量指标，如果手头有分割掩膜（肿瘤、器官），直接算 Dice：

dice = 2 * nnz(fixedMask & movedMask) / (nnz(fixedMask) + nnz(movedMask));

Dice 单一数值报出来，比互信息值更能说服审稿人或临床同事。没有掩膜的话，可以人工在两张图上各点几个对应解剖标志点，算两点间距离（TRE），取均值和中位数。

5. 避坑与排查：配准跑飞、结果全黑、一夜跑不完的五种现场

下面这几条都是我自己或者身边同学真实撞过的，不分先后，按出现频率排。

5.1 输出图全黑或四周一圈黑边

现象：imwarp出来一片黑，或者主体挤在角落，外面一圈黑框。

原因有两个。一是输入图像还是 uint8，变换后插值溢出回到 0 或 255；二是没指定OutputView，输出尺寸按外接矩形算，结果图和参考图对不上。

解决：读图之后立刻im2double，并且统一归一化到[0,1]；imwarp显式传'OutputView', imref2d(size(fixed))。如果黑边确实来自视野外，那是正常现象，测度计算里用掩膜排掉就行，别去改FillValues里面填什么——填 255 只是把问题换个颜色。

5.2 优化器迭代满次数，矩阵越跑越歪

现象：exitflag是 0（到达最大迭代次数），输出变换的旋转角跑到三四十度，图像已经翻得不成样子。

原因：初始错位太大，单尺度下互信息曲面的局部极值把它吸走了；或者参数没有约束，缩放和剪切在搜索时失控。

解决：先把金字塔层数开起来，粗层降采样到 1/8，大位移在这一层就能吃掉；再给参数加上下界，用patternsearch的LowerBound/UpperBound，旋转卡在 ±30°、缩放卡在 0.8～1.2。如果仍然发散，检查目标函数是不是返回了NaN——重叠区被valid掩膜清空时就会这样，优化器收到NaN之后的行为完全不可预测，必须提前用-inf或者一个大正数兜住。

5.3 跨模态用 NCC，怎么调都对不准

现象：CT 配 MR，结果整体偏几个像素，肉眼看骨骼对上了软组织又错位，换参数只是错位方向变。

原因：ncc_measure假设两图灰度线性相关，跨模态这个假设直接不成立。优化器在错误的曲面上找到了「最优」，当然怎么调都差一截。

解决：跨模态一律换互信息。工程里判断模态的办法很简单——同设备同序列用 NCC，不同设备或不同成像原理用 MI。同时也别把 MI 无脑用在单模态上，同模态数据 NCC 的峰值更尖锐、收敛更快，混用是给自己添麻烦。

5.4 互信息算得慢，一晚上跑不完一组

现象：单次配准跑十几分钟，批量跑一晚上还没出结果。

原因：mi_measure里对每个像素做循环、每次迭代都全分辨率重采样、histcounts2每个候选参数都重新分箱。三件事叠起来就是灾难。

解决：先把像素循环改成histcounts2向量化；再用金字塔，粗层用 1/4 分辨率算测度，成本直接降到 1/16；最后用parfor把病例级并行跑起来，MATLAB 默认的并行池开满核，四五个病例同时跑没压力。要是还嫌慢，把mi_measure编译成 MEX，循环里的开销能再砍掉一大半。

5.5 坐标 (x,y) 与 (row,col) 反了

现象：配准结果整体转了 90 度，或者平移方向完全相反，但互信息值居然也不低。

原因：MATLAB 图像矩阵是(row, col)，也就是(y, x)；而affine2d的平移量tx, ty是按(x, y)世界坐标定义的。手写矩阵时把两者混用，平移就斜着去了。

解决：所有跟矩阵索引打交道的地方用(row, col)，所有跟tform、imref2d打交道的地方用(x, y)，中间只在一处转换。更省事的办法是干脆不手写矩阵，用affine2d加imtranslate、imrotate组合，让工具箱替你管坐标顺序。

6. 二次开发：把配准函数塞进自己的批处理流水线

改造的第一步是把入口脚本里的逻辑抽成一个纯函数，签名固定下来，后面所有改动都围绕它：

function [tform, moved, metrics] = reg_one_pair(fixedPath, movingPath, cfg) % fixedPath / movingPath: 图像路径 % cfg: 配置结构体，来自 cfg_default() 再按需覆盖 % tform: affine2d 或位移场；moved: 配准后图像；metrics: 互信息、耗时等 [fixed, moving] = load_pair(fixedPath, movingPath, cfg); switch lower(cfg.mode) case 'rigid', tform = pyramid_reg(fixed, moving, cfg); case 'demons', tform = demons_reg(fixed, moving, cfg); otherwise, error('未知配准模式: %s', cfg.mode); end moved = imwarp(moving, tform, 'OutputView', imref2d(size(fixed))); metrics = struct('mi', mi_measure(fixed, moved, cfg.nbins), ... 'dice', dice_score(fixed, moved, cfg)); end

批量层只做三件事：遍历目录、调reg_one_pair、把metrics写进表格。我一般用parfor加上DataQueue做进度回传，跑的时候能看到「第 7/40 例，平均 MI 提升 23%」，不用干等。

换测度或换优化器时，改动边界要划清楚：

有一点值得单独拎出来说：改测度函数时，返回值的方向必须和原来的约定一致。原来mi_measure是越大越好，objective里取负；你新写一个 SSD，它是越小越好，如果不统一，优化器会朝着最差的方向狂奔，跑出来的图看起来还挺「正常」，只是完全没对齐，属于典型的黑匣子事故。

从那以后我每接一个配准需求，都强制自己先跑通一组肉眼可验的样例数据、把互信息和 Dice 打印出来存进表里，再动代码——没有这个基准，改了什么、改好没改好，全靠感觉，最后一定是白干。希望帮到你。

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