跑完一次肝脏多模态配准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积分之后才是真正的形变。表示形式存储内容单位常见来源直接可视化的可行性位移场每个体素的相对偏移向量毫米ITK、SimpleITK、多数深度学习配准网络可以直接画变形场每个体素的目标物理坐标毫米部分老式工具、某些论文附件减掉恒等映射后可画B 样条系数控制点网格上的系数视实现而定elastix、部分 3D Slicer 模块必须先求值速度场微分同胚的生成元毫米SyN、LDDMM 中间量必须先积分我自己的习惯是拿到任何形变场文件第一件事不是画图而是做一次往返验证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 会自动缩放箭头长度一定要显式设置scale1和scale_unitsxy否则算出来的物理位移和你看到的箭头长度不成比例就没法定量判断了。切片选择有个小技巧不要只画正中间那一层。我会先算整场的位移幅值找出幅值最大的 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 分量分通道 RGBux、uy、uz 各占一通道三维信息完整颜色可解释性差不推荐直接展示箭头场方向 相对大小空间位置明确密集时遮挡严重方向色轮的一个硬伤是只能表达平面方向。对轴状面切片来说z 方向的位移就被丢掉了。我的处理是轴状面看 (ux, uy)同时单独出一张 uz 的灰度图或者干脆在某些关键层上改用箭头场用箭头长度的第三个维度来表达。3.3 箭头场与流线箭头场的实现就是matplotlib的 quiver或者 3D 场景里pyvista的 glyph。它比色轮多一个好处是保留空间位置的直觉你能直接看到“这个位置的向量指向哪里”。实操上有几个参数必须手动调。scale1、scale_unitsxy保证箭头长度等于真实物理长度anglesxy保证箭头方向按数据坐标计算不受屏幕宽高比影响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(vcenter1)蓝色表示压缩、红色表示膨胀然后把 det ≤ 0 的区域单独用高饱和的第三种颜色描边或者填充保证一眼就能定位。一个小细节算梯度的时候 spacing 一定要传对而且是物理单位的 spacing。np.gradient 的第二个参数按轴顺序给轴 0 是 z、轴 1 是 y、轴 2 是 x。这个顺序错一次det(J) 的数值会整体偏掉但图看着还算正常非常阴险。3.5 棋盘格、差值图与闪烁对比前三类方案都是直接看形变场本身棋盘格和差值图属于间接验证——它们看的是配准后图像的吻合程度。棋盘格的做法是把固定图像和变形后的浮动图像按块交替拼接块大小一般取 20 到 40 个体素。如果配准到位块与块之间的解剖结构是连续的如果错位就会在块边界处看到明显的断裂。这个方法的优点是和临床医生的直觉最接近很多放射科同事看棋盘格比看雅可比热图更容易接受。差值图绝对差或带符号差适合看残余误差的分布但要注意它同时混合了配准误差和图像本身的强度差异尤其是多模态配准CT 和 MR 的强度根本没有可比性这时候差值图基本没用只能用棋盘格看结构连续性。我一般的组合是棋盘格看全局吻合雅可比看形变合理性网格看局部形态。这三张图放在一起一个配准结果的好坏基本就无所遁形了。3.6 三维交互查看二维切片再怎么看也只是抽样。真正要把一个形变场看透还是得上交互式三维。我常用的几个工具3D SlicerTransforms 模块可以直接加载 ITK 位移场自带网格显示改色和密度都方便适合和临床同事一起看。ITK-SNAP叠加显示方便网格可视化没那么灵活但胜在轻量。napariPython 生态适合在自己的脚本流程里直接弹出查看配合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(fspacing{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, n21, color#28c8ff, lw0.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], colorcolor, lwlw, alpha0.9) for y in ys: pts np.array([transform.TransformPoint((x, y, z)) for x in xs]) ax.plot(pts[:, 0], pts[:, 1], colorcolor, lwlw, alpha0.9) def show_slice(ax, fixed, z_index, cmapgray, alpha1.0): x0, x1, y0, y1, _ physical_xy_range(fixed, z_index) arr sitk.GetArrayFromImage(fixed)[z_index] ax.imshow(arr, cmapcmap, extent[x0, x1, y0, y1], originlower, alphaalpha) 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, size256): 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_pct99): 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, alpha0.55) im axes[1].imshow(np.ma.masked_where(~m, mag), cmapinferno, extentext, originlower, vmin0, vmaxvmax) axes[1].set_title(displacement magnitude (mm)) fig.colorbar(im, axaxes[1], fraction0.046) show_slice(axes[2], fixed, z_index, alpha0.35) axes[2].imshow(direction_rgb(np.where(m, ux, 0), np.where(m, uy, 0), vmax), extentext, originlower) 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, dpi200, facecolorwhite) plt.close(fig)这段里有两个我踩过坑的地方。一个是np.ma.masked_where——用掩膜数组而不是把掩膜外设成 0因为设成 0 会让背景显示成色图的最低色看起来像是“这个位置没有形变”实际上只是没数据。另一个是vmax用掩膜内的 99 分位而不是全场最大值理由前面讲过了。4.5 雅可比行列式与折叠率统计def jacobian_det(u, spacing): u: (Z, Y, X, 3) 单位 mmspacing: (sx, sy, sz) 单位 mm。返回 (Z, Y, X) 的 det(J)。 sz, sy, sx spacing[2], spacing[1], spacing[0] dz [np.gradient(u[..., c], sz, axis0) for c in range(3)] # d(u_c)/dz dy [np.gradient(u[..., c], sy, axis1) for c in range(3)] # d(u_c)/dy dx [np.gradient(u[..., c], sx, axis2) 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, maskNone): 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_pathNone): os.makedirs(out_dir, exist_okTrue) 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, indent2) 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 不同工具的方向约定差异这是跨工具链最容易出问题的地方。我把自己用过的几种情况整理成一张表但必须强调任何约定差异最终都要用已知点做往返验证来确认不要凭记忆。来源输出形式常见约定使用前必做的检查ANTs*_1Warp.nii.gz物理空间与固定图像同网格用antsApplyTransformsToPoints对一个已知点做正反变换验证SimpleITK 位移场变换DisplacementFieldTransform物理空间毫米定义在固定图像网格确认GetSize()与固定图像一致SimpleITK Warpsitk.Warp的输入位移场定义在输出图像空间与DisplacementFieldTransform互为反向务必区分elastixTransformParameters.*.txtB 样条系数体素相关单位先用transformix -def导出采样后的场再可视化自研网络输出.npz/.h5视实现而定可能是归一化坐标打印幅值分位数确认量纲最稳的做法是统一转换成 SimpleITK 的DisplacementFieldTransform之后所有可视化代码只针对这一种表示写。转换的时候用一个已知的解剖标志点做一次往返验证误差在亚毫米级就认为转换正确。5.6 常见问题速查表现象最可能的原因优先排查动作图像上下颠倒方向矩阵带负元未处理打印GetDirection()加轴翻转网格与解剖错位位移场与图像网格不一致比对GetSize()、GetSpacing()、GetOrigin()位移幅值大得离谱单位是体素不是毫米或背景外推打印幅值分位数加掩膜折叠率异常高梯度计算 spacing 顺序错用刚体场验证 det(J) 是否恒为 1折叠率异常低恒为 0统计范围被掩膜裁得太小打印掩膜体积占比色图几乎全是一个颜色colorbar 量程被背景拉偏改用掩膜内 99 分位截断箭头长度与真实位移不符quiver 自动缩放设scale1、scale_unitsxy三维查看时场显示为空组件数或 dtype 不对转sitkVectorFloat32再读6. 几个我踩坑之后固定下来的习惯写了这么多年可视化脚本有几条经验已经变成了肌肉记忆顺手分享出来。第一条是任何形变场文件到手先算幅值分位数再画图这个五秒钟的操作挡掉的 bug 比任何代码审查都多。第二条是雅可比实现一定要配一个刚体变换的单元测试因为它是整个可视化链路里最容易静默出错的一环——错了不会报异常只会给你一个看起来挺像回事的结论。第三条是任何用于论文的图都先在 25% 缩放下看一眼缩小之后还看得清的结构才值得留在图里。第四条也是我觉得最有价值的一条把调试图和结果图彻底分开做两套脚本不要试图用一套代码同时满足这两种需求最后一定是两边都不好用。至于后续还能怎么扩展我最近在尝试把二维切片热图按关键层自动选层之后拼成一张全景图再配上折叠区域的计数标注做成每个病例一张的“形变场体检报告”这样批量筛查的时候不用切来切去扫一眼就能定位异常病例。