医学图像配准形变场可视化：从位移场到雅可比折叠检测

跑完一次肝脏多模态配准，Dice 刷到 0.94，我把形变场画出来一看——肝门附近的网格拧成了麻花，局部雅可比行列式已经是负的。这种“指标好看、形变场上全是坑”的情况，在医学图像配准的日常里出现频率高得离谱。形变场可视化这件事，说白了就是把配准算法吐出来的那个三维位移向量场，翻译成人的眼睛能判断对错的东西：网格有没有被扯断、器官边界有没有被撕裂、位移幅值是不是超出了解剖上合理的范围。它既是算法研发阶段的调试工具，也是临床质控环节的验收手段，还是论文里最能说明问题的那张图。不管你是刚开始做医学图像配准的研究生，还是已经在做多模态融合、放疗剂量累积、纵向随访萎缩率量化的工程同学，把形变场画好这一关都绕不过去。下面我把这十几年里踩过的坑、写过的脚本、以及一套现在还在用的可视化流程，完整地摊开讲一遍。

1. 形变场可视化到底在解决什么问题

1.1 标量指标看不见的局部畸变

配准评价里最常被引用的几个数：Dice、HD95、TRE、互信息。它们有一个共同的毛病——把整个体积压成了一个数。一个体积几百万个体素，只要大结构对得上，Dice 就好看；但真正决定配准能不能用的，往往是那些体积占比很小、却对临床意义极大的局部区域。

举个我自己经手过的场景：腹部 CT-MR 配准，肝脏整体 Dice 0.95，看起来相当漂亮。可是把形变场叠上去之后发现，肝门区靠近血管的地方出现了明显的折叠，负雅可比体素占比 0.7%。0.7% 听起来很小，但折算成体积是好几立方厘米，如果拿这个形变场去做射频消融的术前规划，穿刺路径就可能被映射到错误的位置上。

这就是形变场可视化的第一层价值：它把“平均意义上正确”拆开成“每个位置是否正确”。标量指标告诉你及格了，可视化告诉你哪里不及格。

顺便说一句，Dice 还有个隐蔽的问题——它对大结构有天然偏袒。肺配准里气管、小血管的分支，脑配准里海马、丘脑这类小核团，它们的体积可能只占整体的千分之几，就算全错了 Dice 也不会掉多少。所以只要你的配准结果要进入下游任务，形变场可视化就不是“可选项”，而是“必选项”。

1.2 形变场在数学上到底是什么

很多人画不好形变场，根子不在绘图技巧，而在于没搞清楚手里那个数组到底代表什么。形变场在工程上有好几种完全不同的表示方式，混用一次就够你查一整天 bug。

最常见的是位移场（displacement field）。它对每个体素存一个三维向量 u(x)，实际的空间映射写成 φ(x) = x + u(x)。这种表示最直观，也最常见，缺点是只适合表达“不太离谱”的形变，位移大到一定程度的时候，正向映射可能不再是单射。

另一种是变形场或绝对位置场（warp field / transformation field），它存的是每个体素最终映射到的绝对物理坐标，而不是相对偏移。两者之间差一个恒等映射，转换写起来就一行代码，但方向搞反了整个图会完全错乱。

还有两种在特定工具链里出现得更多：B 样条控制点系数（elastix 的参数文件就是这个），以及速度场（LDDMM、SyN 这类微分同胚算法内部迭代的中间量）。B 样条系数不能直接当位移场画，必须先经过一次基函数求值或者用transformix -def转成采样后的场；速度场也不能直接画，要通过指数映射（scaling-and-squaring）积分之后才是真正的形变。

我自己的习惯是：拿到任何形变场文件，第一件事不是画图，而是做一次往返验证（round-trip）。挑几个已知的解剖标志点，用形变场把它们从固定图像空间映射到浮动图像空间，再映射回来，看误差是不是接近零。如果回不来，说明表示方式或者方向约定理解错了，这时候画出来的图再漂亮也是错的。

1.3 三类读者，三种不同的图

形变场可视化没有“一张万能图”。我在项目里会按读者分三档来做，而且这三档的图长得完全不一样。

第一档是给我自己看的调试图。要的是信息密度最大，好不好看完全无所谓。我会把原始位移场的切片、幅值热图、方向色轮、雅可比行列式、棋盘格拼在一张大图上，一次刷 20 个测试病例，找异常。

第二档是给临床同事或项目评审看的质控图。要的是“一眼能判断有没有问题”。通常是固定图像做灰度背景，上面叠稀疏网格，旁边配一个雅可比行列式的双色热力图，红蓝色标出折叠区，再附一行统计数字：折叠率、位移 95 分位、最大位移。

第三档是给论文用的插图。要的是清晰、克制、可复现。网格线不能太密，颜色不能太艳，字号要能缩小到单栏宽度还看得清，色轮图例必须单独出一张。

把这三档分开，能省下大量返工时间。我见过太多人把调试图直接塞进论文，结果审稿人一句“the deformation grid is too dense to interpret”就得重画。

1.4 影响范围：从算法迭代到临床决策

形变场可视化的影响链条其实很长。在研发阶段，它决定了你调参的效率——没有它，你只能靠 Dice 的零点几个百分点的波动盲猜；有了它，一眼就能看出是正则化太弱还是太强。在验证阶段，它决定了你的结论可不可信，因为平均指标会掩盖失败样本。在应用阶段，它直接决定了下游任务的安全性，比如放疗中把剂量从计划 CT 映射到多次分次的 CBCT 上，折叠区域的剂量映射会产生几倍的系统误差。

还有一类容易被忽略的影响：纵向随访。做脑萎缩率测量的时候，形变场的雅可比行列式直接就是局部体积变化率。如果可视化没做好，一个折叠区就能让某个脑区的萎缩率虚高十几个百分点，而你从最终数字上完全看不出问题。

2. 数据准备与坐标系归一化

2.1 位移场的存储格式与读取要点

常见的形变场文件格式有这么几种：.nii.gz、.mha、.nrrd、.npz、.h5。前三者是带空间信息的医学图像格式，后两者是裸数组，没有空间信息——这点非常重要，因为裸数组丢失了 spacing、origin 和方向矩阵，你没法把位移向量和真实解剖位置对应起来。

判断一个文件是不是位移场，看通道数：矢量图（vector image）的组件数（components per pixel）是 3。ITK 系工具读进来之后，组件顺序是 (x, y, z)，而数组维度顺序是 (z, y, x)。这两个顺序反着，是新手最容易翻车的地方之一。

import SimpleITK as sitk import numpy as np fixed = sitk.ReadImage("fixed.nii.gz", sitk.sitkFloat32) warp = sitk.ReadImage("warp.nii.gz", sitk.sitkVectorFloat32) print("size :", fixed.GetSize()) print("spacing:", fixed.GetSpacing()) print("origin :", fixed.GetOrigin()) print("direction:\n", np.round(np.array(fixed.GetDirection()).reshape(3, 3), 4)) print("components:", warp.GetNumberOfComponentsPerPixel()) u = sitk.GetArrayFromImage(warp) # 形状 (Z, Y, X, 3)，最后一位是 (ux, uy, uz)，单位毫米 fixed_arr = sitk.GetArrayFromImage(fixed) assert u.shape[:3] == fixed_arr.shape, "位移场和固定图像网格不一致"

提示：位移场的 dtype 一定要转成 float32 或 float64 再参与任何运算。有些工具会把它存成 int16 或者 uint8，直接读进来算梯度会得到一堆台阶状的噪声。

还有一个检查项是位移场所在的网格。位移场必须和固定图像同网格，这在 ITK 体系里是硬性要求。如果你拿到一个和浮动图像同网格的场，所有画出来的一切都是错的，而且错得很隐蔽——图看着有形状，就是和背景对不上。

2.2 物理坐标与体素坐标：翻车重灾区

体素索引 (i, j, k) 到物理坐标 p 的转换公式是 p = O + D · S · i，其中 O 是原点，D 是 3×3 的方向矩阵，S 是对角线上放着 spacing 的缩放矩阵。对着公式手写代码看着简单，实际项目里我统计过，形变场可视化画错的原因至少一半出在这上面。

三个高频错误：

把体素单位的位移当成毫米。有些工具（尤其是一些自己写的 PyTorch 实现）输出的位移场是归一化到 [-1, 1] 或者以体素为单位的。如果 spacing 是 (0.7, 0.7, 2.0) 这种各向异性的情况，直接用会让 z 方向的形变看起来比实际小 3 倍。

把 (z, y, x) 当成 (x, y, z)。numpy 数组是 (z, y, x)，ITK 的物理坐标是 (x, y, z)。取切片的时候一次疏忽，画出来的就是冠状面当轴状面用。

忽略了方向矩阵。NIfTI 文件里 radiological 和 neurological 朝向的差别，会让图像左右翻转。更麻烦的是有些扫描协议的方向矩阵带负的对角元，直接 imshow 出来的图像是上下颠倒的，而你拿它当背景叠网格，就会觉得网格和脑沟完全对不上。

最省事的规避办法是：可视化之前先把所有数据重采样到一个规范的 LPS 网格上，方向矩阵设成单位阵，spacing 全部取正数。代价是一次重采样引入的插值误差（用线性插值基本可以忽略），换来的是后面所有代码都只需要处理最简单的情况。我现在几乎所有可视化脚本的第一段都是这个规范化步骤。

def to_canonical(img): """把图像重采样到方向矩阵为单位阵、spacing 为正的规范网格。""" d = np.array(img.GetDirection()).reshape(3, 3) sp = np.array(img.GetSpacing()) # 只处理轴对齐的情况，带旋转的方向矩阵建议直接用原始坐标变换 flip = np.sign(np.diag(d)) new_sp = sp * np.abs(np.diag(d)) new_sp[new_sp == 0] = sp[new_sp == 0] return sitk.DirectionCorrection # 占位，实际用 sitk.Resample + 手工构造参考网格

实际做的时候，我会用sitk.ResampleImageFilter加一个手工构造的参考网格（size 用GetSize()，spacing 用GetSpacing()，direction 用单位阵，origin 用TransformIndexToPhysicalPoint((0,0,0))），把图像和位移场分别重采样过去，这样位移场重采样时用线性插值就够，不要用高阶插值。

注意：位移场重采样千万别用三次以上的插值。高阶插值会在场里造出过冲，肉眼看着平滑，算雅可比的时候会出现假的折叠。

2.3 采样抽稀：网格、箭头、切片怎么定

一个 512×512×300 的体积，位移场有七千多万个向量。全画出来既不可能也没意义。抽稀策略直接决定了图能不能读。

网格图的密度我用“物理间距”来控制，而不是体素间隔。默认每 8 到 10 毫米一条线，脑部配准这种小器官可以收到 5 毫米，腹部可以放到 15 毫米。换算成线条数量，一张 200 毫米宽的切片上放 20 到 25 条线比较舒服。线太密会有摩尔纹，太疏看不出局部形变。

箭头场的采样间隔一般比网格再稀一点，我常用 12 到 16 个体素一个箭头。箭头密度对可视化效果的影响比箭头本身的长度大得多。matplotlib里默认的 quiver 会自动缩放箭头长度，一定要显式设置scale=1和scale_units='xy'，否则算出来的物理位移和你看到的箭头长度不成比例，就没法定量判断了。

切片选择有个小技巧：不要只画正中间那一层。我会先算整场的位移幅值，找出幅值最大的 3 个轴向切片位置，再加上器官的质心层，一共画 4 到 6 层。事实证明，90% 的形变异常都出现在幅值最大的那几层里。

2.4 掩膜与值域裁剪：别让背景毁掉 colorbar

这是我想重点强调的一条。位移场在背景空气区域往往是完全没意义的——那里的图像几乎没有梯度信息，配准算法在那里输出的位移纯粹是正则化项外推的结果，可能会达到几十甚至上百毫米。如果你画幅值热图的时候不做处理，colorbar 的量程会被这些背景值拉到 100 毫米以上，器官内部 2~3 毫米的真实差异就被压成了一片均匀的颜色。

我的处理方式是两步走。第一步，用一个器官掩膜（或者至少是图像强度阈值生成的粗略前景掩膜）把统计范围限制住。第二步，用一个稳健的分位数截断，比如取掩膜内位移幅值的 99 分位作为 colorbar 上限，超出的部分直接饱和处理。这样器官内部的结构能看得非常清楚。

雅可比行列式的可视化同理，但更要注意：背景区域的雅可比数值往往极端，因为它们是在噪声上求梯度。做折叠率统计的时候，一定要说明统计范围是掩膜内还是全体积，这两者报出来的数字可能差十倍。论文里我建议明确写 “folding rate within the organ mask”。

3. 五种常用可视化方案与实现细节

3.1 网格形变图：最直观也最容易骗人

网格形变图的逻辑很朴素：在固定图像空间里画一组规则的正交直线，然后用形变场把这些线的采样点映射到目标空间，连起来。规则性被破坏的程度就是形变的直观读数。

它最大的优点是人眼对“直线变弯”极其敏感，只要有一点点局部剪切或者旋转，立刻就能看出来。缺点是它只反映采样点上的形变。如果你的网格稀到 20 毫米一条线，中间发生的小范围折叠（比如 5 毫米大的区域）可能完全被漏掉。我遇过最坑的一次是：网格图干干净净，但把雅可比热图叠上去，掩膜内有 0.3% 的折叠区，位置正好落在两条网格线中间。

所以我的做法是网格图和雅可比热图必须配着看。网格图负责解释形变的形态，雅可比负责保证不漏检。

另外网格图有个视觉陷阱：有些形变看起来很大，其实只是整体平移造成的网格偏移，并不是局部畸变。判断局部畸变要看网格线的相对形变，也就是间距是否均匀变化、角度是否保持。这一点在写图注的时候最好说明一下，不然审稿人可能会误读。

3.2 位移幅值热图与方向色轮

幅值热图是最容易做也最容易被误读的图。它把三维向量取模变成一个标量，信息损失很大，但胜在直观。我一般把它当作“找异常位置”的索引图，而不是定量图。

真正承载方向信息的是方向色轮编码。做法是把位移向量的平面分量 (ux, uy) 转成极坐标，角度映射到 HSV 的色相 H，幅值映射到明度 V。这样一张图里，颜色告诉你往哪个方向变形，亮度告诉你变形的强度。旁边配一张标准色轮图例，读者就能直接读出方向。

方向色轮的一个硬伤是只能表达平面方向。对轴状面切片来说，z 方向的位移就被丢掉了。我的处理是：轴状面看 (ux, uy)，同时单独出一张 uz 的灰度图；或者干脆在某些关键层上改用箭头场，用箭头长度的第三个维度来表达。

3.3 箭头场与流线

箭头场的实现就是matplotlib的 quiver，或者 3D 场景里pyvista的 glyph。它比色轮多一个好处是保留空间位置的直觉，你能直接看到“这个位置的向量指向哪里”。

实操上有几个参数必须手动调。scale=1、scale_units='xy'保证箭头长度等于真实物理长度；angles='xy'保证箭头方向按数据坐标计算，不受屏幕宽高比影响；width控制在 0.001 到 0.003 之间，太粗会糊成一片；颜色我习惯用纯白或纯黄，叠在灰度解剖图上对比度最好。

流线（streamline）用得少一些，一般只在需要展示大范围位移轨迹的时候用。用matplotlib.streamplot的话要先把位移场转成规则的速度网格，注意它只支持均匀网格，非均匀 spacing 要先插值。

3.4 雅可比行列式：量化折叠的硬指标

如果只能保留一种形变场可视化方案，我选雅可比行列式。因为它是唯一能给出定量、可比较、可自动化筛查结论的方案。

形变映射 φ 的雅可比矩阵 J = ∂φ/∂x，对位移场来说就是 J = I + ∂u/∂x。它的行列式 det(J) 有非常明确的物理意义：局部体积变化率。det(J) = 1 表示局部体积不变，大于 1 表示膨胀，小于 1 表示压缩，小于等于 0 表示发生了折叠（folding），映射不再是一一对应，物理上不可实现。

这几个阈值在项目里可以直接当验收标准用：

• det(J) ≤ 0 的体素占比（折叠率）：腹部配准我一般要求小于 0.1%，脑部要求小于 0.01%。
• det(J) 的 1% 和 99% 分位数：用来判断压缩和膨胀是否过度。压缩到 0.5 以下或者膨胀到 1.5 以上，通常说明正则化参数该调了。
• det(J) 的标准差：全局形变平滑性的一个粗指标。

可视化上，我一般用双色 diverging 色图，以 1 为中心（用TwoSlopeNorm(vcenter=1)），蓝色表示压缩、红色表示膨胀，然后把 det ≤ 0 的区域单独用高饱和的第三种颜色描边或者填充，保证一眼就能定位。

一个小细节：算梯度的时候 spacing 一定要传对，而且是物理单位的 spacing。np.gradient 的第二个参数按轴顺序给，轴 0 是 z、轴 1 是 y、轴 2 是 x。这个顺序错一次，det(J) 的数值会整体偏掉，但图看着还算正常，非常阴险。

3.5 棋盘格、差值图与闪烁对比

前三类方案都是直接看形变场本身，棋盘格和差值图属于间接验证——它们看的是配准后图像的吻合程度。

棋盘格的做法是把固定图像和变形后的浮动图像按块交替拼接，块大小一般取 20 到 40 个体素。如果配准到位，块与块之间的解剖结构是连续的；如果错位，就会在块边界处看到明显的断裂。这个方法的优点是和临床医生的直觉最接近，很多放射科同事看棋盘格比看雅可比热图更容易接受。

差值图（绝对差或带符号差）适合看残余误差的分布，但要注意它同时混合了配准误差和图像本身的强度差异（尤其是多模态配准，CT 和 MR 的强度根本没有可比性，这时候差值图基本没用，只能用棋盘格看结构连续性）。

我一般的组合是：棋盘格看全局吻合，雅可比看形变合理性，网格看局部形态。这三张图放在一起，一个配准结果的好坏基本就无所遁形了。

3.6 三维交互查看

二维切片再怎么看也只是抽样。真正要把一个形变场看透，还是得上交互式三维。我常用的几个工具：

• 3D Slicer：Transforms 模块可以直接加载 ITK 位移场，自带网格显示，改色和密度都方便，适合和临床同事一起看。
• ITK-SNAP：叠加显示方便，网格可视化没那么灵活，但胜在轻量。
• napari：Python 生态，适合在自己的脚本流程里直接弹出查看，配合napari的 vectors layer 可以实时看三维向量场。
• pyvista：做出版级三维渲染最方便，glyph 和 streamlines 都是一行代码的事。

交互查看的正确定位是“抽查”，不是“批量筛查”。几百个病例你还是得靠脚本自动算折叠率和分位数，把可疑病例挑出来，再用交互工具细看。指望人工一个个点开看，做二十个病例你就崩溃了。

4. 代码实操：一套可复用的绘图脚本

4.1 环境与依赖

我现在的标准组合是 SimpleITK 读数据、numpy 算场、matplotlib 出图。不依赖深度学习框架，装起来干净，在服务器上没有图形界面也能跑（matplotlib.use("Agg")）。

pip install SimpleITK numpy scipy matplotlib

如果要做三维渲染再加pyvista和napari。我建议把绘图逻辑和配准逻辑彻底解耦，绘图脚本只接收“固定图像路径 + 位移场路径 + 输出目录”三个参数。这样无论是 ANTs、elastix 还是自己训的网络产出的场，只要格式统一，同一套脚本就能用。

4.2 读取与坐标系校验

import numpy as np import SimpleITK as sitk def load_pair(fixed_path, warp_path): fixed = sitk.ReadImage(fixed_path, sitk.sitkFloat32) warp = sitk.ReadImage(warp_path, sitk.sitkVectorFloat32) assert warp.GetNumberOfComponentsPerPixel() == 3, "不是位移场" assert warp.GetSize() == fixed.GetSize(), "网格不一致，先重采样" u = sitk.GetArrayFromImage(warp).astype(np.float64) # (Z, Y, X, 3) fa = sitk.GetArrayFromImage(fixed).astype(np.float32) return fixed, fa, u def sanity_check(fixed, u): sp = np.array(fixed.GetSpacing()) # (sx, sy, sz) mag = np.linalg.norm(u, axis=-1) # (Z, Y, X) p = np.percentile(mag, [50, 95, 99, 100]) print(f"spacing={sp}, size={fixed.GetSize()}") print(f"位移幅值(mm) 中位={p[0]:.2f} 95%={p[1]:.2f} 99%={p[2]:.2f} 最大={p[3]:.2f}") # 粗略判断单位是不是搞错了：最大位移超过 200mm 基本可以确定有问题 if p[3] > 200: print("警告：位移幅值异常大，检查单位是否为毫米、是否包含背景外推") return mag

这段 sanity check 我建议每次跑都打印。它只要五秒钟，但能挡掉大部分“单位搞错”和“场和图像不匹配”的低级错误。实测下来，最大位移超过 200 毫米的情况，九成以上是单位问题或者背景外推，不是配准真的那么差。

4.3 网格形变图

import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt def physical_xy_range(img, z_index): size = img.GetSize() o = np.array(img.TransformIndexToPhysicalPoint((0, 0, z_index))) e = np.array(img.TransformIndexToPhysicalPoint((size[0] - 1, size[1] - 1, z_index))) x0, x1 = sorted((o[0], e[0])) y0, y1 = sorted((o[1], e[1])) return x0, x1, y0, y1, o[2] def draw_grid(ax, transform, fixed, z_index, n=21, color="#28c8ff", lw=0.7): x0, x1, y0, y1, z = physical_xy_range(fixed, z_index) xs = np.linspace(x0, x1, n) ys = np.linspace(y0, y1, n) for x in xs: pts = np.array([transform.TransformPoint((x, y, z)) for y in ys]) ax.plot(pts[:, 0], pts[:, 1], color=color, lw=lw, alpha=0.9) for y in ys: pts = np.array([transform.TransformPoint((x, y, z)) for x in xs]) ax.plot(pts[:, 0], pts[:, 1], color=color, lw=lw, alpha=0.9) def show_slice(ax, fixed, z_index, cmap="gray", alpha=1.0): x0, x1, y0, y1, _ = physical_xy_range(fixed, z_index) arr = sitk.GetArrayFromImage(fixed)[z_index] ax.imshow(arr, cmap=cmap, extent=[x0, x1, y0, y1], origin="lower", alpha=alpha) ax.set_xlim(x0, x1) ax.set_ylim(y0, y1) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([])

transform这里既可以是sitk.DisplacementFieldTransform(u_img)，也可以是组合变换sitk.CompositeTransform([...])。用TransformPoint而不是自己写矩阵乘法，好处是方向矩阵、origin、spacing 全部由 ITK 处理，你不用操心。

需要注意一点：网格点变形之后可能跑到切片平面之外（因为 z 方向也有位移），所以ax.plot画出来的其实是这些点投影到 z = 常数平面的结果。对轴状面来说这个投影误差通常很小，但如果 z 方向位移很大，最好改用真正的三维网格在三维视图里看。

配色上我试过很多组合，最后固定在深色背景（#101418）+ 高饱和青色网格线 + 灰度解剖图。这个组合在投影仪和打印稿上都还看得清。纯白网格叠在骨窗 CT 上会糊掉，不建议用。

4.4 位移幅值热图与方向色轮

from matplotlib.colors import hsv_to_rgb, Normalize import matplotlib.cm as cm def direction_rgb(ux, uy, vmax): ang = np.arctan2(uy, ux) h = (ang / (2 * np.pi)) % 1.0 v = np.clip(np.hypot(ux, uy) / vmax, 0, 1) s = np.ones_like(h) return hsv_to_rgb(np.stack([h, s, v], axis=-1)) def save_colorwheel(path, size=256): yy, xx = np.mgrid[-1:1:complex(size), -1:1:complex(size)] h = (np.arctan2(yy, xx) / (2 * np.pi)) % 1.0 r = np.clip(np.hypot(xx, yy), 0, 1) s = np.ones_like(h) plt.imsave(path, hsv_to_rgb(np.stack([h, s, r], axis=-1))) def plot_mag_and_dir(fixed, u, z_index, mask, out_path, vmax_pct=99): ux = u[z_index, :, :, 0] uy = u[z_index, :, :, 1] mag = np.hypot(ux, uy) m = mask[z_index] if mask is not None else np.ones_like(mag, bool) vmax = np.percentile(mag[m], vmax_pct) x0, x1, y0, y1, _ = physical_xy_range(fixed, z_index) ext = [x0, x1, y0, y1] fig, axes = plt.subplots(1, 3, figsize=(15, 5)) show_slice(axes[0], fixed, z_index) axes[0].set_title("fixed image") show_slice(axes[1], fixed, z_index, alpha=0.55) im = axes[1].imshow(np.ma.masked_where(~m, mag), cmap="inferno", extent=ext, origin="lower", vmin=0, vmax=vmax) axes[1].set_title("displacement magnitude (mm)") fig.colorbar(im, ax=axes[1], fraction=0.046) show_slice(axes[2], fixed, z_index, alpha=0.35) axes[2].imshow(direction_rgb(np.where(m, ux, 0), np.where(m, uy, 0), vmax), extent=ext, origin="lower") axes[2].set_title("direction (hue) + magnitude (value)") for ax in axes: ax.set_xticks([]); ax.set_yticks([]) fig.tight_layout() fig.savefig(out_path, dpi=200, facecolor="white") plt.close(fig)

这段里有两个我踩过坑的地方。一个是np.ma.masked_where——用掩膜数组而不是把掩膜外设成 0，因为设成 0 会让背景显示成色图的最低色，看起来像是“这个位置没有形变”，实际上只是没数据。另一个是vmax用掩膜内的 99 分位，而不是全场最大值，理由前面讲过了。

4.5 雅可比行列式与折叠率统计

def jacobian_det(u, spacing): """u: (Z, Y, X, 3) 单位 mm；spacing: (sx, sy, sz) 单位 mm。返回 (Z, Y, X) 的 det(J)。""" sz, sy, sx = spacing[2], spacing[1], spacing[0] dz = [np.gradient(u[..., c], sz, axis=0) for c in range(3)] # d(u_c)/dz dy = [np.gradient(u[..., c], sy, axis=1) for c in range(3)] # d(u_c)/dy dx = [np.gradient(u[..., c], sx, axis=2) for c in range(3)] # d(u_c)/dx # 行索引 c 对应输出分量，列索引 j 对应输入坐标 (x, y, z) du_dx = np.stack([dx, dy, dz], axis=-1) # (Z, Y, X, 3, 3) J = du_dx + np.eye(3) return np.linalg.det(J) def folding_report(detJ, mask=None): v = detJ[mask] if mask is not None else detJ.ravel() v = v[np.isfinite(v)] return { "folding_rate": float((v <= 0).mean()), "p01": float(np.percentile(v, 1)), "p99": float(np.percentile(v, 99)), "mean": float(v.mean()), "std": float(v.std()), }

用np.stack([dx, dy, dz], axis=-1)这个写法要特别小心。最后那个 axis=-1 是列索引，也就是输入坐标 x、y、z 的顺序。如果写成axis=-2，得到的矩阵是转置的，行列式虽然还算得出来，但数值完全不对（对称矩阵下甚至看不出来）。我在这个坑里待过半天，最后是靠一个已知的刚体旋转场验证才发现的——单位阵加旋转的雅可比行列式应该恒等于 1，如果算出来不是 1，那一定是矩阵布局错了。这个方法推荐你也在自己的脚本里加一个单元测试。

提示：验证雅可比实现是否正确的黄金标准就是刚体变换。构造一个已知的小角度旋转加平移的位移场，理论上 det(J) 处处为 1。如果你的实现在这个测试上偏离超过 1e-4，先别往下走。

4.6 批量导出与论文排版参数

我现在的流程是：一个病例一个子目录，脚本自动跑一遍，输出grid_z*.png、mag_z*.png、jac_z*.png、checker_z*.png和一份report.json。所有病例跑完，再写一个汇总脚本读所有 report.json，按折叠率排序，把前 10% 的病例挑出来人工复查。

import json, os def run_case(fixed_path, warp_path, out_dir, mask_path=None): os.makedirs(out_dir, exist_ok=True) fixed, fa, u = load_pair(fixed_path, warp_path) mag = sanity_check(fixed, u) mask = sitk.GetArrayFromImage(sitk.ReadImage(mask_path)).astype(bool) if mask_path else None detJ = jacobian_det(u, fixed.GetSpacing()) rep = folding_report(detJ, mask) rep["n_folding_voxels"] = int((detJ <= 0).sum()) with open(os.path.join(out_dir, "report.json"), "w") as f: json.dump(rep, f, indent=2) return rep

论文排版参数我固定成这几条：单栏图宽度 8.5 厘米、双栏 17 厘米；dpi 300 起步（矢量图优先存 PDF）；字体用 Arial 或 Helvetica，字号在最终显示尺寸下不小于 7pt；网格线宽 0.6~0.8pt；色图优先选感知均匀的（inferno、viridis、cividis），绝对不要用 jet——jet 会在中间位置造出一个虚假的亮带，让人误以为那里形变最大。

5. 常见问题排查与避坑清单

5.1 整张图左右反了或者上下颠倒

症状很好认：网格的形状大体合理，但和背景解剖结构对不上，脑室跑到了另一侧，或者网格整体镜像了。

排查顺序是这样的。先打印固定图像的GetDirection()，看看对角线上有没有负值，再看 sform 和 qform 是否一致（用nibabel加载后对比affine和get_qform()）。如果方向矩阵带负元，imshow出来的数组需要按对应轴翻转，翻转代码我前面show_slice里已经给了。还有一种情况是位移场和固定图像的方向矩阵不同（比如一个是 LPS 一个是 RAS），这种直接重采样规范到同一个网格最省事。

5.2 网格看着很正常，但配准其实是失败的

这是最危险的一类问题，因为它会让你误以为配准没问题。典型场景是：网格整体非常平滑，几乎没有畸变，但棋盘格上解剖结构明显错位。

原因通常是配准算法压根没动——可能是优化器早停了，也可能是初始变换太强、形变项被正则化压死了。判断方法很简单：算一下位移幅值的中位数。如果中位数接近 0（比如小于 0.5 毫米），而两个图像之间明显存在需要大形变才能对齐的差异，那基本就是配准失败了。

我的处理是给每个病例加一条硬性检查：位移幅值中位数小于 1 毫米且 Dice 低于阈值的时候，直接标记为失败病例。这个规则救过我很多次。

5.3 雅可比出现大面积负值

大面积负值分两种性质完全不同的情况。

第一种是真折叠，也就是配准算法确实产生了不合理的形变，通常伴随位移场在某些位置出现剧烈跳变。这时候去看网格图，往往能看到网格线交叉。处理办法是加强正则化（B 样条配准里调FinalGridSpacingInPhysicalUnits和SP_a，深度学习配准里加位移场的梯度惩罚项），或者改用微分同胚的算法。

第二种是数值假象。如果负值只出现在掩膜外或者图像边缘，而且数量很稀疏，大概率是背景噪声上求梯度导致的。这时候用掩膜限制统计范围就消掉了。还有一种情况是位移场本身带噪声（比如网络输出的场没有做后处理平滑），梯度会把高频噪声放大好几倍。可以先对位移场做一次小尺度的高斯平滑（sigma 约 1 个体素）再算雅可比，注意平滑只用于可视化和统计，不要把平滑后的场写回去当作配准结果。

5.4 箭头太密、网格太乱、论文里看不清

这个问题几乎每个人都会遇到。我的经验值是：单张切片上，网格线控制在 15 到 25 条之间，箭头控制在 20×20 个以内。如果一定要展示高密度形变，宁可放大一个局部区域（ROI 裁剪），也不要硬塞。

还有一个技巧是分层展示：一张图只讲一件事。展示全局形变就放稀疏网格；展示局部异常就放大 ROI 配雅可比；展示方向就用色轮。想在一张图里把什么都说完，结果就是什么都看不清。

5.5 不同工具的方向约定差异

这是跨工具链最容易出问题的地方。我把自己用过的几种情况整理成一张表，但必须强调：任何约定差异，最终都要用已知点做往返验证来确认，不要凭记忆。

最稳的做法是：统一转换成 SimpleITK 的DisplacementFieldTransform，之后所有可视化代码只针对这一种表示写。转换的时候用一个已知的解剖标志点做一次往返验证，误差在亚毫米级就认为转换正确。

5.6 常见问题速查表

6. 几个我踩坑之后固定下来的习惯

写了这么多年可视化脚本，有几条经验已经变成了肌肉记忆，顺手分享出来。第一条是任何形变场文件到手，先算幅值分位数再画图，这个五秒钟的操作挡掉的 bug 比任何代码审查都多。第二条是雅可比实现一定要配一个刚体变换的单元测试，因为它是整个可视化链路里最容易静默出错的一环——错了不会报异常，只会给你一个看起来挺像回事的结论。第三条是任何用于论文的图，都先在 25% 缩放下看一眼，缩小之后还看得清的结构才值得留在图里。第四条，也是我觉得最有价值的一条：把调试图和结果图彻底分开做两套脚本，不要试图用一套代码同时满足这两种需求，最后一定是两边都不好用。至于后续还能怎么扩展，我最近在尝试把二维切片热图按关键层自动选层之后拼成一张全景图，再配上折叠区域的计数标注，做成每个病例一张的“形变场体检报告”，这样批量筛查的时候不用切来切去，扫一眼就能定位异常病例。