Python | 非线性配准基本原理

《Python | ANTs 多模态医学影像配准》

《Python | PyElastix 非刚性配准》

但这些只是学会了使用现成的工具，对于非线性配准的基本原理并没有深入介绍。

这篇文章会直观地介绍：当两张图像存在非线性形变时，算法怎样从图像内容中估计出一个连续的非线性形变场？

一、非线性配准到底想求什么？

配准里通常有两张图：

fixed image:  Fmoving image: M

目标是找一个空间变换，让 moving 经过变换后尽量像 fixed。

在线性配准里，这个变换可能只是平移、旋转、缩放：

x' = A x + b

但非线性配准允许不同位置有不同位移：

φ(x) = x + u(x)

这里的 u(x) 就是形变场。在二维图像里，每个像素有一个二维位移：

u(x, y) = (dy, dx)

所以非线性配准的核心可以写成：

找到一个形变场 u(x)，让 M(x + u(x)) 尽量接近 F(x)

也就是：

warped moving = moving(x + u(x))

但这里有一个重要问题：如果每个像素都能动，怎样在这么高的自由度下，找到合适的形变场？

答案是：非线性配准不是让每个像素独立乱动，而是在优化一个“整体平滑的空间变换”。常见目标函数可以粗略写成：

总损失 = 图像不相似程度 + 形变正则化

也就是：

E(u) = D(F, M(x + u(x))) + λ R(u)

其中：

• D 衡量 warped moving 和 fixed 是否相似；

• R 惩罚形变场太毛躁、太不连续；

• λ 控制“贴得更紧”和“形变更平滑”之间的权衡。

这就是非线性配准最核心的思想：既要对齐图像，也要让形变合理。

二、非线性形变示例

下面这个例子里，fixed 是 Cameraman 原图，moving 在fixed的基础上有扭曲。

• 第一张是 fixed；

• 第二张是 moving；

• 第三张是二者的初始差异，也叫做residual；

• 第四张是 fixed 和 moving 的 overlay。

人、相机、三脚架的边缘都有错位。我们的目标就是找到一个形变场，把 moving 拉回 fixed。

三、通过梯度得到形变场

假设当前已经有一个形变场 u，我们可以把 moving warp 到 fixed 空间：

warped = moving(x + u(x))

然后计算 residual：

residual = fixed - warped

一轮更新里会看到这些中间结果：

C图中的箭头。表示 warped moving 中“往哪里采样会让灰度变亮最快”。

D图代表每个像素根据自己的 residual 和 gradient 算出来的局部移动建议。

Residual 很重要，但它只告诉我们“差多少”：

这个位置 warped moving 比 fixed 太亮，或者太暗

它本身并不告诉我们采样坐标应该往左、右、上、下哪个方向挪。方向来自 warped moving 的图像梯度。

图像梯度不是猜出来的，而是邻近像素灰度差。对当前 warped moving 图像 W，中心差分可以写成：

∂W/∂x ≈ (右边像素 - 左边像素) / 2∂W/∂y ≈ (下边像素 - 上边像素) / 2

Python 里可以直接写：

gy, gx = np.gradient(warped)

其中 gx 是左右方向灰度变化，gy 是上下方向灰度变化。

梯度回答的是： 如果我把采样坐标稍微挪一点，warped moving 的灰度会怎么变？

所以 gradient 不是位移本身，而是“位移会怎样影响灰度”。

现在我们有两个量：

residual：当前还差多少，是需要变亮还是变暗gradient：往哪个方向采样，灰度会变亮最快

于是局部更新可以理解为：

δu ∝ residual × gradient

• 如果 residual > 0，说明 warped moving 太暗，需要变亮，就沿 gradient 方向采样。 

• 如果 residual < 0，说明 warped moving 太亮，需要变暗，就沿 gradient 反方向采样。 

• 如果 residual 接近 0，这个位置基本不用动。

二维图像里，梯度由左右方向和上下方向两个分量组成，所以它本身就是一个二维向量：

梯度向量指向灰度增加最快的方向，而 residual 决定最终更新是顺着梯度还是逆着梯度。

四、关键代码：warp、更新场、平滑和累积

先看 warp。配准里的 warp 本质上是采样：对 fixed 空间里的每个位置 (y, x)，去 moving 里的 (y + dy, x + dx) 取灰度值。

import numpy as npfrom scipy.ndimage import map_coordinatesdef warp_image(image, disp):    rows, cols = np.indices(image.shape)    coords = [        rows + disp[..., 0],  # y + dy        cols + disp[..., 1],  # x + dx    ]    return map_coordinates(image, coords, order=1, mode="nearest")

接下来写一个简化的局部更，展示residual 和 gradient 如何组合成 update field：

def local_update(fixed, warped, step=0.25):    residual = fixed - warped    gy, gx = np.gradient(warped)    denom = gx * gx + gy * gy + 0.45 * residual * residual + 1e-5    update = np.zeros((*fixed.shape, 2), dtype=float)    update[..., 0] = step * residual * gy / denom    update[..., 1] = step * residual * gx / denom    return update, residualfrom scipy.ndimage import gaussian_filterdef smooth_field(field, sigma):    smoothed = np.empty_like(field)    for k in range(2):        smoothed[..., k] = gaussian_filter(field[..., k], sigma=sigma)    return smoothed

直接算出来的局部更新场通常会很毛躁，因为每个像素只看自己附近的 residual 和 gradient。比如一条边缘附近，左边像素可能建议“往左一点”，右边像素可能建议“往右一点”；纹理和噪声附近还会给出很多很碎的方向。真实配准不能让相邻像素各走各的路，所以要平滑更新场：

这也解释了为什么图里的 raw update 和 smoothed update 看起来差别很大。Raw update 保留了很多高频、零散、边缘驱动的局部建议；smooth 以后，这些局部建议被周围邻域平均，孤立的小方向会被削弱，大片区域共同的移动趋势会被保留下来。

Smooth本质是在表达配准里的正则化假设：真实形变应该连续、平滑，相邻位置应该大体一起移动。

最后，小更新要累积成总形变场。严格来说，形变场应该通过变换组合来更新，而不是简单相加。一个简化的 displacement field 组合可以写成：

def warp_vector(field, disp):    rows, cols = np.indices(field.shape[:2])    coords = [rows + disp[..., 0], cols + disp[..., 1]]    warped = np.empty_like(field)    for k in range(2):        warped[..., k] = map_coordinates(            field[..., k], coords, order=1, mode="nearest"        )    return warpeddef compose_displacements(current, update):    return update + warp_vector(current, update)

把这些合在一起，主循环就是：

u = np.zeros((*fixed.shape, 2), dtype=float)for iteration in range(num_iterations):    warped = warp_image(moving, u)    update, residual = local_update(fixed, warped)    update = smooth_field(update, sigma=update_sigma)    u = compose_displacements(u, update)    u = smooth_field(u, sigma=total_sigma)

这段循环就是“形变场逐渐长出来”的核心。 每一轮只做一点点局部更新，然后平滑、累积。很多轮之后，dense deformation field 就形成了。

五、迭代更新形变场

可以把形变场想象成一张坐标纸被拉弯。规则网格变成弯曲网格，说明空间坐标被局部拉伸、压缩或弯曲；网格仍然连续，说明形变不是每个像素各走各的。

下面这张图展示了形变场逐步发挥作用的过程：

最终结果如下：

residual 相关指标随 iteration 的变化：

这里包括：

MAE：平均绝对误差RMSE：均方根误差Corr：fixed 和 warped moving 的相关系数

这个例子里：

MAE:  0.0507 -> 0.0105Corr: 0.9246 -> 0.9954

说明 warped moving 确实越来越接近 fixed。

warped moving 与变形网格随 iteration 的变化：

六、实际非线性配准

本文用 residual × gradient 生成局部更新，是为了把直观展示怎样“将图像差异怎样变成移动方向”。不代表所有非线性配准都必须使用像素差。

实际方法通常还会选择不同的组成部分：

刚性或仿射初始配准多分辨率金字塔不同的图像相似性度量不同的形变参数化方式不同的平滑或正则化约束不同的优化与变换组合方式

不同软件和算法会用不同方式回答这些问题。例如：

• ANTs 是一套常用的医学图像配准工具，SyN 是其中著名的对称可微分同胚配准方法；SyN 通常使用局部互相关 CC，或者 MI 等metric

• Demons 是另一类经典的稠密配准思路。

• Elastix 非线性配准会使用 B-spline transform。此时算法优化的不是每个像素的二维位移，而是一组较稀疏的 B-spline 控制点参数

最后总结一下：

1. 形变场 u(x) 告诉我们在 moving 里去哪里采样。2. residual 告诉我们当前 warped moving 和 fixed 差多少。3. gradient 告诉我们采样坐标往哪移动会改变灰度。4. residual × gradient 给出局部更新方向。5. 平滑和正则化让相邻位置一起动。6. 多轮小更新累积成 dense deformation field。

所以形变场不是逐像素 random search 出来的，也不是算法凭空猜出来的。它是通过“图像差异 + 图像梯度 + 平滑正则化 + 多轮累积”一步步得到的。

这就是非线性配准最核心的直觉。

本篇教程的所有Python代码，可以在GitHub上下载：

https://link.zhihu.com/?target=https%3A//github.com/ethanzhao9/Medical-Image-Processing

希望对大家有帮助~

创作不易，点个关注再走吧

如果有任何问题，欢迎在文章下方留言

作者 | Treasure琛

排版 | 帽子妈妈

往期回顾

Python | 实现 3D Slice Viewer

Python | 读取显微镜图像格式

Python | pyStackReg 图像配准

Python | 小鼠标准脑图谱(CCFv3)下载和使用

Python | PyElastix 非刚性配准