简介：面向医学影像处理学习者与科研人员的MATLAB非刚性图像配准代码包，围绕不同时间、设备或成像方式获取的医学图像对齐问题，提供从刚体预处理、二维/三维B样条变形建模，到互信息相似度度量与梯度优化求解的完整实现思路。压缩包共34个文件，约240KB，以22个m脚本、6个C源码为主，辅以png示例图像、fig界面配置和txt说明；m文件承载配准示例、误差计算与结果显示，C源码实现底层变换和互信息直方图计算，便于理解算法与程序之间的衔接。已有1336人学习，适合具备一定MATLAB基础、希望钻研医学图像配准与多模态相似度度量的读者，可据此复现实验流程、观察配准效果，并在此基础上调整变形模型和优化参数。
1. 一张 CT 和一张 MR 叠不上：医学图像配准真正要解决的是什么
影像科把同一患者的 CT 和 MR 摆在一起看，十有八九是错位的：层厚一个 1 mm 一个 5 mm，摆位一个仰头一个平躺，FOV 也不一样。医学图像配准要干的事很朴素——求一个空间变换，把 moving 图搬进 fixed 图的坐标系，让同一块解剖结构落在同一个像素上。MATLAB 里这条链路是现成的：geometric transform 描述变换，metric 给相似度打分，optimizer 负责搜索，imwarp 做重采样，四个零件拼起来就是imregtform和imregister。
它最适合的场景是算法验证和课题复现：手上有几十例数据，想在一两天内跑出对比曲线，不想为了编译某个 C++ 库折腾一下午。反过来，千万级体素、需要 GPU 加速、或者要接深度学习配准网络的活儿，这套工具链就不划算了——它胜在参数透明、每一步都能停下来看中间结果，而不是胜在速度。
2. 变换、度量、优化器：医学图像配准三件套的选型分界线
配准跑不出结果，八成不是代码写错，是这三件套配错了。变换的自由度给多了，优化器会在噪声里乱逛；度量选错了，单模态能收敛的配置换到跨模态直接发散。这一章把三者的边界划清楚，后面写脚本才不会靠试。
2.1 变换模型：刚体、仿射、B 样条各自负责什么
先把自由度算清楚，再谈选型。同一个器官在不同时间点扫，形变主要是体位造成的整体旋转平移，用刚体就够；不同模态之间还存在尺度差异（像素间距不同、重建 FOV 不同），必须上仿射；软组织器官的呼吸形变、脑组织的开颅后形变，只能靠稠密位移场。
imregister家族支持的变换只有'translation'、'rigid'、'similarity'、'affine'四种，B 样条不在里面，形变配准要另走imregdemons。
| 变换类型 | 2D 自由度 | 3D 自由度 | MATLAB 入口 | 典型场景 | 不该用的场景 |
|---|---|---|---|---|---|
| translation | 2 | 3 | imregtform(...,'translation',...) | 同机同时序、只差位移 | 有旋转的摆位差 |
| rigid | 3 | 6 | 'rigid' | 头颅 CT/MR、骨性标志 | 有尺度差的两台设备 |
| similarity | 4 | 7 | 'similarity' | 尺度未知的跨设备 | 各向异性拉伸 |
| affine | 6 | 12 | 'affine' | CT-MR、PET-CT | 软组织局部形变 |
| 稠密位移场 | 2×N | 3×N | imregdemons | 呼吸、形变、图谱映射 | 大位移未预对齐 |
一个常见误用是拿仿射去啃腹部呼吸配准：仿射只能表达全局线性形变，肝脏下缘对上了，膈顶又跑了。正确做法是刚体或仿射先做粗对齐，再用imregdemons收尾，这一点在第 3 章的脚本里会体现。
2.2 相似性度量：为什么跨模态必须换互信息
Image Processing Toolbox 里内置的 metric 只有两个类：registration.metric.MeanSquares和registration.metric.MattesMutualInformation。前者就是逐像素灰度差平方和，本质假设两图灰度线性相关；后者统计两图灰度的联合直方图，衡量的是"知道 moving 的灰度后，fixed 灰度的不确定性下降了多少"。
CT 的灰度是亨斯菲尔德单位，骨头 1000 以上、空气 -1000；MR 的灰度取决于序列和线圈，同一块软组织在 T1 和 T2 上可以一亮一暗。这两个东西之间根本不存在线性映射，用 MeanSquares 就是在优化一个错误的假设，结果往往是"配准完成但完全没对上"——优化器很诚实，它只是把错误的目标函数降到了局部最小。
互信息的关键参数有三个，改错任何一个都会让结果飘：
NumberOfHistogramBins：联合直方图分箱数。8 bit 图像 32～64 够用，12 bit 数据用 32 会粗到分辨不出峰，一般提到 64～100。NumberOfSpatialSamples：每次迭代随机采样的像素数。默认 20000，图像大、噪声重时提到 50000，代价是单次迭代变慢。UseAllPixels：设为true则用全部像素，可复现性好但慢得多；调试阶段建议false，出最终结果时再打开做一次确认。
提示：互信息对重叠区域极度敏感。两图 FOV 差很多时，背景（空气）占了直方图的大部分样本，优化器会把"让背景重叠"当成最优解。裁到共同 ROI 再配，比调任何参数都有效。
2.3 优化器：梯度下降和进化算法不是随便挑的
imregconfig('monomodal')返回RegularStepGradientDescent，imregconfig('multimodal')返回OnePlusOneEvolutionary。这个默认配置背后是有道理的：单模态的相似度曲面相对光滑，梯度可信，梯度下降收敛快；互信息的曲面是一堆锯齿，梯度方向基本是噪声，只能用无导数的随机搜索。
RegularStepGradientDescent的四个参数决定了它会不会"冲过头"：
| 参数 | 默认量级 | 作用 | 调大的后果 |
|---|---|---|---|
MaximumStepLength | 0.0625 | 初始步长 | 直接跳过最优点，来回震荡 |
MinimumStepLength | 1e-5 | 停止阈值 | 迭代次数暴涨，收益趋零 |
RelaxationFactor | 0.5 | 每步衰减系数 | 接近 1 时收敛慢，太小则过早停 |
MaximumIterations | 100 | 迭代上限 | 只影响耗时，不影响精度 |
OnePlusOneEvolutionary看的是另外三个：InitialRadius（初始搜索半径）、GrowthFactor（半径增长因子，默认 1.05）、Epsilon（收敛判据）。跨模态配准跑飞，最常见的原因就是InitialRadius用了默认值而两图初始错位超过了几十毫米——搜索半径还没长到那么大，迭代次数就耗完了。
2.4 预处理与多分辨率金字塔：一半的精度在这里决定
我一般把预处理固定成四步，缺一步后面都要返工：
- 统一数据类型：
im2double把整型转成 [0,1]，避免整型运算溢出。 - 灰度拉伸：
mat2gray做 min-max 归一化，让窄灰度带的图像也有足够对比度参与直方图统计。 - 视野对齐：
imcrop裁到两图共同覆盖的解剖区域。 - 挂物理坐标：用
imref2d/imref3d把 DICOM 的PixelSpacing带进去，否则算出来的平移量是"像素"不是"毫米"，跟临床对不上。
金字塔层数由'PyramidLevels'控制，默认 3 层。它的作用是把大位移拆成"粗层找大体位置、细层抠细节"两件事。层数太少，大旋转会掉进局部极小；层数太多（比如 5 层以上），最粗那层的图像已经被高斯平滑糊成一团，反而找不到方向。2D 数据 3 层、3D 数据 3～4 层是稳妥区间。
3. 从 imregtform 到 imregdemons：四段能直接跑的脚本
参数讲完，落到代码。下面四段可以按顺序拼成一个完整的处理流程，也可以单独摘出来用。所有示例都基于灰度图，RGB 要先转灰度。
3.1 读数据与预处理：让两张图先站在同一起跑线
% 读同一病例的 CT 与 MR 单层，统一成 double 并拉伸灰度 fixedInfo = dicominfo('case01_ct.dcm'); movingInfo = dicominfo('case01_mr.dcm'); fixed = dicomread(fixedInfo); moving = dicomread(movingInfo); % im2double 负责类型归一，mat2gray 负责对比度拉伸，两步不要合并 fixed = mat2gray(im2double(fixed)); moving = mat2gray(im2double(moving)); % 先裁到共同解剖区域，再挂物理坐标，顺序反了 ref 就对不上 fixed = imcrop(fixed, [40 40 200 200]); moving = imcrop(moving, [30 30 210 210]); % PixelSpacing 是 [行间距 列间距]，而 imref2d 要的是 [列 行] fixedRef = imref2d(size(fixed), fixedInfo.PixelSpacing(2), fixedInfo.PixelSpacing(1)); movingRef = imref2d(size(moving), movingInfo.PixelSpacing(2), movingInfo.PixelSpacing(1));这段的关键在最后两行的顺序。DICOM 标签PixelSpacing的存储顺序是行优先，imref2d的构造参数是列优先，直接照抄会导致 X/Y 间距互换，图像各向异性的时候误差能有百分之几。裁剪必须放在建 ref 之前，否则 ref 里记录的尺寸还是裁剪前的，后面imwarp的OutputView会算错。
imcrop的矩形参数是[xmin ymin width height]，以左上角为原点，两个病例的裁剪框不必相同，只要保证覆盖的是同一块解剖结构就行。
3.2 单模态刚体配准：imregconfig 的默认值够不够用
% 单模态：同一次扫描的两个序列，灰度分布接近，用均方差即可 [optimizer, metric] = imregconfig('monomodal'); % 默认 MaximumIterations=100 对临床数据偏少，脑部配准常需要 200~300 optimizer.MaximumIterations = 300; optimizer.MinimumStepLength = 1e-5; optimizer.MaximumStepLength = 0.0625; tformRigid = imregtform(moving, fixed, 'rigid', optimizer, metric, ... 'PyramidLevels', 3, 'DisplayOptimization', true); % spacing 不同的情况下必须走带 ref 的重采样，否则只按像素算 movingRigid = imwarp(moving, movingRef, tformRigid, 'OutputView', fixedRef); figure; subplot(1,2,1); imshowpair(fixed, moving, 'falsecolor'); title('配准前'); subplot(1,2,2); imshowpair(fixed, movingRigid, 'falsecolor'); title('配准后');'DisplayOptimization', true会把每次迭代的相似度值打印到命令行，这是排查"到底收敛没收敛"最直接的手段——曲线如果在最后几十次迭代里还在上下大幅跳动，说明MaximumStepLength给大了，得往下调。
imwarp的两种调用形式差别很大。省略 ref 时，MATLAB 按像素索引处理，变换矩阵里的平移分量单位是像素；带上movingRef和'OutputView', fixedRef之后，平移分量才有物理意义。这一步不做，后面所有"位移了多少毫米"的结论都是错的。
3.3 跨模态 CT-MR：互信息度量配进化优化器
% 跨模态：CT 与 MR 灰度无线性关系，必须换互信息 [optimizer, metric] = imregconfig('multimodal'); metric.NumberOfHistogramBins = 64; % 12 bit 数据 32 太粗 metric.NumberOfSpatialSamples = 50000; % 采样多一点，MI 估计更稳 metric.UseAllPixels = false; % 调试阶段关掉，速度优先 optimizer.InitialRadius = 6.25e-3; % 初始错位大时适当放大 optimizer.GrowthFactor = 1.05; optimizer.Epsilon = 1.5e-6; optimizer.MaximumIterations = 500; % 用刚体结果当初值，别让进化优化器从零开始瞎猜 initAffine = affine2d(tformRigid.T); tformMI = imregtform(moving, fixed, 'affine', optimizer, metric, ... 'InitialTransformation', initAffine, 'PyramidLevels', 3); movingMI = imwarp(moving, movingRef, tformMI, 'OutputView', fixedRef);'InitialTransformation'是这段里最值钱的一行。跨模态配准的互信息曲面局部极多，直接跑仿射经常落在某个"看着还行"的伪解上；先跑刚体拿到大致的旋转平移，再把它转成affine2d作为初值，收敛率和精度都会明显改善。转换用tformRigid.T提齐次矩阵再构造，比手工复制参数可靠。
UseAllPixels = false意味着每次迭代随机采样，同一份数据跑两遍结果会有微小差异。做论文对比实验时，务必在最后统一设成true，或者在调用前rng(0)固定随机种子，否则别人复现不出你的数字。
3.4 形变配准：imregdemons 怎么接在刚体之后
% 形变配准只做最后一公里，输入必须是已经刚体/仿射对齐的图 [D, movingDemons] = imregdemons(movingMI, fixed, [100 50 25], ... 'AccumulatedFieldSmoothing', 1.5, 'DisplayWaitbar', false); % 检查雅可比行列式：detJ <= 0 说明位移场把像素折叠了 [dx_dcol, dx_drow] = gradient(D(:,:,1)); [dy_dcol, dy_drow] = gradient(D(:,:,2)); detJ = (1 + dx_dcol) .* (1 + dy_drow) - dx_drow .* dy_dcol; fprintf('折叠体素占比: %.4f%%\n', 100 * mean(detJ(:) <= 0)); figure; imshowpair(fixed, movingDemons, 'checkerboard');imregdemons的第二个参数[100 50 25]是三个尺度上的迭代次数，从粗到细。AccumulatedFieldSmoothing控制每轮之后对位移场做多大强度的平滑：小于 1.0 时位移场会追着噪声跑，出现大量局部折叠；大于 3.0 时形变被抹平，等于白做。1.0～2.0 是比较稳的区间，解剖形变剧烈的部位取小值。
雅可比行列式这一步别省。折叠体素占比超过千分之一，说明位移场已经不可逆，拿去算剂量分布或者做形变映射会得到物理上讲不通的结果。真出现折叠，先降AccumulatedFieldSmoothing到 1.0，再不行就把最细尺度的迭代次数从 25 降到 10。
4. 配准结果怎么验收：别让"看着对齐"骗了你
配准最容易翻车的地方不是跑不出结果，而是跑出了一个看起来对、实际上错的结果。灰阶图叠在一起，人眼对 2～3 mm 的错位其实很不敏感，必须靠指标和交叉检查。
4.1 定量指标：MI、NCC、Dice 与 TRE 的适用边界
| 指标 | 计算方式 | 适用条件 | 数值解读 |
|---|---|---|---|
| 互信息 MI | 联合直方图熵差 | 跨模态 | 越大越好，无绝对阈值，只看前后对比 |
| 归一化互相关 NCC | normxcorr2或自算 | 单模态、线性灰度 | 0～1，0.9 以上算对齐良好 |
| Dice 系数 | 标签交并比的 2 倍 | 有分割标签 | 0.85 以上可用于统计 |
| TRE | 标志点配对距离 | 有金标准 | 直接是毫米误差，最可信 |
| 雅可比行列式 | 位移场偏导行列式 | 形变配准 | ≤0 的比例应接近 0 |
单模态任务里 NCC 和 MI 会同时上升，可以互相印证；跨模态任务里 NCC 没有意义，只看 MI。Dice 需要人工或自动分割，标注成本高但对临床最直观。TRE 是唯一能给出"误差几毫米"的指标，问题是你得先有可信的标志点——解剖标志、植入物、或者仿真数据里已知的变换。
dice在 MATLAB 里有现成实现（Image Processing Toolbox），输入是两个逻辑矩阵：
% 用配准后图像的分割结果和 fixed 的分割结果算 Dice diceScore = dice(fixedMask, movingMaskWarped); fprintf('Dice = %.4f\n', diceScore);注意movingMaskWarped必须用最近邻插值重采样，标签图用双线性插值会造出 0.37 这种不存在的类别。
4.2 目视检查：棋盘格、差值图与轮廓叠加
imshowpair的几种模式各有用途，别只用'falsecolor'：
'checkerboard'：棋盘格交替显示两图，最适合看边缘是否连续。错位 1 mm 以上，骨皮质和脑室边缘会出现明显的台阶。'diff'：绝对差值图，单模态配准时差值越接近全黑越好，跨模态基本没用。'blend'：半透明叠加，适合看整体轮廓。'montage'：左右并排，适合放在汇报材料里。
figure; subplot(1,3,1); imshowpair(fixed, movingMI, 'checkerboard'); title('棋盘格'); subplot(1,3,2); imshowpair(fixed, movingMI, 'diff'); title('差值图'); subplot(1,3,3); imshowpair(fixed, movingMI, 'blend'); title('半透明叠加');4.3 变换矩阵怎么读：从 T 反推旋转角与平移毫米数
imregtform返回的是affine2d或rigid2d对象，.T是 3×3（2D）或 4×4（3D）的齐次矩阵。会读这个矩阵，才能判断结果是否符合解剖常识。
T = tformMI.T; % 3x3，最后一列恒为 [0;0;1] theta = atan2d(T(2,1), T(1,1)); % 旋转角，单位度 scaleX = hypot(T(1,1), T(2,1)); % X 方向尺度因子 scaleY = hypot(T(1,2), T(2,2)); % Y 方向尺度因子 tx_mm = T(3,1) * movingRef.PixelExtentInWorldX; ty_mm = T(3,2) * movingRef.PixelExtentInWorldY; fprintf('旋转 %.2f 度, 平移 (%.2f, %.2f) mm, 尺度 (%.3f, %.3f)\n', ... theta, tx_mm, ty_mm, scaleX, scaleY);MATLAB 的仿射变换采用行向量右乘约定，也就是[x y 1] * T，所以平移在第三行而不是第三列。这一点跟很多教材上的列向量约定正好转置，照抄公式会得到完全相反的结论。
读出来之后对照常识：头部 CT-MR 的旋转一般不超过 15 度，平移不超过 30 mm，尺度因子应该非常接近 1。如果跑出 90 度旋转、或者尺度因子 0.7，几乎可以断定是优化器跳到了伪解，别急着往下做。
4.4 回放验证：用已知变换测一遍自己的流程
在真实数据上没人知道正确答案，那就自己造一个。取一张图，施加一个已知的刚体变换，再用同一套流程把它配回来，看残余误差有多大。这是唯一能在没有金标准的情况下量化流程精度的办法。
% 人为施加 8 度旋转 + (12,-7) 像素平移 gt = rigid2d('Rotation', 8, 'Translation', [12 -7]); movingWarped = imwarp(moving, gt, 'OutputView', imref2d(size(moving))); [opt, mtr] = imregconfig('monomodal'); tEst = imregtform(movingWarped, moving, 'rigid', opt, mtr); % 复合估计变换与已知变换，看四角残余位移 corners = [1 1; size(moving,2) 1; 1 size(moving,1); size(moving,2) size(moving,1)]; pGT = transformPointsForward(gt, corners); pEst = transformPointsForward(tEst, pGT); residual = sqrt(sum((pEst - corners).^2, 2)); fprintf('四角残余位移(像素): %.2f %.2f %.2f %.2f\n', residual);四角残余位移都在 0.5 像素以内，说明这套参数配置在这个数据尺度上是可信的；某个角跑到 5 像素以上，多半是那个角落的解剖结构对比度太低，优化器没有约束住。这个测试我一般在换数据集之后必跑一遍，比直接上真实数据瞎调参数省时间。
5. 避坑与排查：五条血泪记录，从全黑输出到配准跑飞
下面五条基本都是我在不同数据集上真实撞过的，按"现象 → 原因 → 解决"记，遇到时可以直接对号入座。
5.1 显示与数据类的两个坑
现象一：imwarp输出的图整张全黑，或者一半黑一半正常。
原因通常是三选一：输出视图没指定，导致重采样后的图像落在原图坐标系之外；输入是uint16而变换后的插值结果超出显示范围；或者用了imshow而数据是single类型且最小值不为 0。
解决：显式指定'OutputView'，并且重采样后统一做一次拉伸。
warped = imwarp(moving, movingRef, tformMI, 'OutputView', fixedRef); warped = mat2gray(warped); % 拉回 [0,1]，避免显示全黑现象二：平移量算出来是几百毫米，明显不合理。
原因是没挂物理坐标，imwarp和imregtform全程按像素索引计算，而你按毫米去解读了。另一可能是PixelSpacing的行列顺序取反，各向异性数据（比如 0.7×0.7 vs 1.2×1.2）下会放大成几倍的误差。
解决：建 ref 时严格用PixelSpacing(2)对应 X、PixelSpacing(1)对应 Y，并在读出平移量之后乘以PixelExtentInWorldX/Y做一次单位换算，两头对得上才算完。
5.2 优化与参数类的两个坑
现象三：配准跑完，imshowpair一看，两个器官完全错开，相似度曲线还告诉你收敛了。
原因是相似度度量选错（跨模态用了 MeanSquares）或者初始错位超出了优化器的搜索范围。RegularStepGradientDescent的搜索是局部的，初始偏差超过十几毫米就基本没救；OnePlusOneEvolutionary的InitialRadius如果只有 6.25e-3，搜索半径要长很多轮才够大。
解决：先确认度量类型与模态关系匹配，再给'InitialTransformation'一个粗配准结果，最后才考虑放大InitialRadius。三者的优先级不能反，先调参数往往是白费功夫。
现象四：同一份数据跑两次，结果差了几个像素。
原因是UseAllPixels = false下的随机采样，每次迭代抽的像素不同，MI 的估计值有抖动。数据量小或者噪声重的时候，这个抖动足以让优化器走向不同的局部极小。
解决：调试阶段接受抖动但固定rng(0)；出最终结果时把UseAllPixels设为true，代价是耗时增加几倍，换来的是可复现性。
5.3 形变场类的坑
现象五：imregdemons跑完，某些区域的组织像被"揉皱"了，像素互相穿透。
原因是位移场出现了负雅可比，也就是折叠。常见诱因有三个：输入没有做刚体预对齐，大位移让 demons 的局部假设失效；AccumulatedFieldSmoothing设得太小（比如 0.5），位移场追着噪声跑；图像本身有大面积均匀区域（如膀胱、脑室），梯度接近零，位移场在那里无法定义。
解决：先刚体后形变，这是硬顺序；平滑参数从 1.5 起步往回调；对均匀区域可以先做掩膜，配准后只在掩膜内应用位移场。验收时一定跑第 3.4 节那段雅可比检查，折叠比例超过千分之一就重来。
6. 进阶：多分辨率加 B 样条，把 2 mm 的残差压到 1 mm 以内
单次调用imregtform的天花板很明显，脑部 CT-MR 的 TRE 通常在 2 mm 上下，再想往下压，靠调参数已经没用了，得改流程结构。我常用的做法是三级串联：刚体 → 仿射（带初值）→ 形变，每一级的输出作为下一级的输入，中间用'InitialTransformation'传递变换，而不是把上一级的重采样图像再喂进去。这么做的差别在于：串联变换不会因为多次插值而累积模糊，图像只重采样一次。
% 三级串联：变换逐级复合，图像只在最后重采样一次 [optMono, mtrMono] = imregconfig('monomodal'); t1 = imregtform(moving, fixed, 'rigid', optMono, mtrMono, 'PyramidLevels', 4); [optMulti, mtrMulti] = imregconfig('multimodal'); mtrMulti.NumberOfHistogramBins = 64; t2 = imregtform(moving, fixed, 'affine', optMulti, mtrMulti, ... 'InitialTransformation', affine2d(t1.T), 'PyramidLevels', 3); % 复合前两级，得到形变配准的输入 movingCoarse = imwarp(moving, movingRef, affine2d(t2.T), 'OutputView', fixedRef); movingCoarse = mat2gray(movingCoarse); % 形变场用三级尺度，平滑系数逐级递减 [D, ~] = imregdemons(movingCoarse, fixed, [200 100 50], ... 'AccumulatedFieldSmoothing', 1.2); movingFinal = imwarp(movingCoarse, D);这里有个细节值得单独说：imregdemons的迭代次数[200 100 50]和AccumulatedFieldSmoothing是耦合的。粗尺度迭代多、平滑强，先把大形变抓住；细尺度迭代少、平滑弱，只抠局部。如果反过来（粗尺度迭代少、细尺度迭代多），位移场会在细尺度上反复修正同一块区域，最后折叠一片。
参数配置我一般按这张表起步，再根据雅可比检查结果微调：
| 阶段 | 关键参数 | 起步值 | 调整方向 |
|---|---|---|---|
| 刚体 | PyramidLevels | 4 | 大旋转时加到 4，稳定数据用 3 |
| 仿射 | NumberOfHistogramBins | 64 | 12 bit 数据提到 100 |
| 仿射 | MaximumIterations | 500 | 收敛曲线未平就加到 800 |
| 形变 | 迭代序列 | [200 100 50] | 形变剧烈时[300 150 75] |
| 形变 | AccumulatedFieldSmoothing | 1.2 | 折叠比例高就往上调，最大 2.0 |
验证方法还是 4.4 节的回放测试，但这次要测两级：先测刚体+仿射的残余误差，再测加上形变之后的。如果加了形变 TRE 反而变大，说明形变在"过拟合"图像噪声，把AccumulatedFieldSmoothing往上调一档，或者直接砍掉最细的那一级迭代。
3D 数据还要额外注意内存。一个 512×512×200 的single体积已经占 200 MB，imregdemons内部会存位移场和多个金字塔层，8 GB 内存的机器跑到第三级就可能开始用交换分区，速度断崖式下跌。我的做法是先在降采样一半的数据上跑刚体和仿射，把变换矩阵记下来，再回到原分辨率只跑形变——刚体仿射对分辨率不敏感，形变才需要全分辨率。
最后说个习惯。早期我总想一步到位，调一次参数就指望出结果，结果每次换数据集都要重来一遍，浪费的时间比写代码还多。后来固定成"先跑回放测试定参数，再上真实数据，最后必看雅可比"这三步，虽然每次多花十分钟，但再也没有出现过做完一周才发现配准结果是错的这种情况。前面那套脚本你可以直接抄，但那个验证顺序，建议你也一并抄走。希望帮到你。
本文还有配套的精品资源，点击获取
