LNCC：局部归一化互相关损失与 PyTorch 实现

医学图像配准不仅要让两幅图像在视觉上更相似，还要建立正确的解剖对应，并保证形变场在几何上合理。局部归一化互相关（Local Normalized Cross-Correlation，LNCC）只回答其中一个问题：固定图像与配准后的移动图像，在局部灰度变化模式上是否一致？

这个边界很重要。LNCC 可以作为单模态无监督配准的训练损失，也可以作为图像层面的评价指标，但更高的 LNCC 并不必然代表解剖结构已经正确对齐，更不能单独说明形变保持了拓扑。

1. 配准目标中的图像相似性

记固定图像为 $I_f$，移动图像为 $I_m$，变换为 $\phi$。按照本文采用的坐标映射约定，形变后的移动图像写为：

$$I_w(p)=I_m\bigl(\phi(p)\bigr)=(I_m\circ\phi)(p).$$

无监督形变配准通常同时优化图像相似性与形变正则项：

$$\mathcal{L}_{reg} = \mathcal{L}_{sim}(I_f,I_w) + \lambda\mathcal{R}(\phi).$$

$\mathcal{L}_{sim}$ 约束图像内容对齐，$\mathcal{R}(\phi)$ 则抑制剧烈、不连续或不合理的形变。LNCC 属于前者，它并不直接约束位移场的平滑性、可逆性或 Jacobian 行列式。

2. 从 NCC 开始

对于单模态图像，固定图像与移动图像通常具有相近的灰度含义。归一化互相关（Normalized Cross-Correlation，NCC）先对两幅图像去均值，再通过各自的方差归一化：

$$\operatorname{NCC}(I_f,I_w) = \frac{ \sum_{p\in\Omega} \bigl(I_f(p)-\bar I_f\bigr) \bigl(I_w(p)-\bar I_w\bigr) }{ \sqrt{ \sum_{p\in\Omega}\bigl(I_f(p)-\bar I_f\bigr)^2 } \sqrt{ \sum_{p\in\Omega}\bigl(I_w(p)-\bar I_w\bigr)^2 } +\varepsilon }.$$

其中，$\Omega$ 是离散像素或体素网格，$\bar I_f$ 与 $\bar I_w$ 是全局平均强度，$\varepsilon$ 用于避免分母为零。

从统计角度看，NCC 就是标准化协方差：

$$\operatorname{NCC}(I_f,I_w) = \frac{ \operatorname{Cov}(I_f,I_w) }{ \sqrt{ \operatorname{Var}(I_f)\operatorname{Var}(I_w) } +\varepsilon }.$$

如果两幅图像满足近似的正线性强度关系，NCC 会接近 $1$；若局部变化方向相反，则可能接近 $-1$。去均值和方差归一化使 NCC 对整体亮度偏移与正比例缩放相对不敏感，但它仍然使用一组全局统计量。

3. 为什么需要局部 NCC

医学图像中的强度关系往往并非全局一致。线圈敏感度、偏置场、扫描协议以及预处理差异，都可能使同一种组织在不同位置呈现不同的灰度偏移。全局 NCC 会把整幅图像压缩成一组均值、方差和协方差，容易忽略局部结构是否真正对齐。

LNCC 将统计范围限制在每个位置 $p$ 周围的邻域窗口 $\mathcal{N}(p)$。先定义局部均值：

$$\mu_f(p) = \frac{1}{|\mathcal{N}(p)|} \sum_{q\in\mathcal{N}(p)} I_f(q), \qquad \mu_w(p) = \frac{1}{|\mathcal{N}(p)|} \sum_{q\in\mathcal{N}(p)} I_w(q).$$

局部协方差的未归一化形式为：

$$C_{fw}(p) = \sum_{q\in\mathcal{N}(p)} \bigl(I_f(q)-\mu_f(p)\bigr) \bigl(I_w(q)-\mu_w(p)\bigr).$$

对应的局部方差为：

$$V_f(p) = \sum_{q\in\mathcal{N}(p)} \bigl(I_f(q)-\mu_f(p)\bigr)^2, \qquad V_w(p) = \sum_{q\in\mathcal{N}(p)} \bigl(I_w(q)-\mu_w(p)\bigr)^2.$$

于是，位置 $p$ 的局部相关系数可以写为：

$$\rho(p) = \frac{C_{fw}(p)} {\sqrt{V_f(p)V_w(p)+\varepsilon}}.$$

最后对有效位置求平均：

$$\operatorname{LNCC}(I_f,I_w) = \frac{1}{|\Omega'|} \sum_{p\in\Omega'}\rho(p).$$

$\Omega'$ 表示两个局部窗口都具有足够灰度变化的位置集合。均匀背景中的局部方差接近零，此时相关系数缺乏统计意义，实践中通常通过 $\varepsilon$、有效窗口掩膜或前景掩膜处理。

4. 指标与损失的符号约定

LNCC 是“越大越好”的相似性指标，而优化器默认最小化损失。常见转换有两种：

$$\mathcal{L}_{LNCC}=1-\operatorname{LNCC},$$

或：

$$\mathcal{L}_{LNCC}=-\operatorname{LNCC}.$$

两者只相差一个常数，不会改变梯度方向，但日志中的数值范围不同。报告实验结果时，应明确展示的是 LNCC 指标还是 LNCC loss。

一些配准实现使用平方相关形式：

$$\rho^2(p) = \frac{C_{fw}^2(p)}{V_f(p)V_w(p)+\varepsilon}.$$

平方形式数值非负，却会把强负相关也视为高相关。对于预期具有正相关强度关系的单模态任务，这通常不是问题；如果符号本身具有意义，则应保留未平方的相关系数。不同论文中所谓的 “LNCC loss” 可能对应不同版本，复现实验时需要先确认公式和代码。

5. 窗口尺度决定了 LNCC 看见什么

局部窗口不是一个无关紧要的超参数，它决定了相似性度量的观察尺度。

• 窗口过小：局部方差容易接近零，对噪声和插值误差更敏感，优化信号可能不稳定。
• 窗口过大：LNCC 会逐渐接近全局 NCC，对局部强度漂移和细粒度错位的敏感性下降。
• 多尺度配准：可以在不同分辨率上采用不同窗口，或组合多个窗口尺度的 LNCC。
• 各向异性体数据：窗口大小最好结合体素间距理解。相同的 $9\times9\times9$ 体素窗口，在不同扫描间距下对应的物理范围可能完全不同。

因此，窗口应该与目标解剖结构的尺度、图像分辨率和网络当前优化层级共同确定，而不是机械地固定为某个体素数。

6. PyTorch 实现

局部统计量可以通过全一卷积核高效计算。下面的实现同时支持二维与三维张量，输入形状分别为 [B, C, H, W] 和 [B, C, D, H, W]。边界位置使用实际参与统计的像素或体素数量，而不是把补零区域计入窗口。

from collections.abc import Sequence

import torch
import torch.nn as nn
import torch.nn.functional as F


class LNCCLoss(nn.Module):
    """Local normalized cross-correlation loss for 2D or 3D images."""

    def __init__(
        self,
        window_size: int | Sequence[int] = 9,
        squared: bool = False,
        eps: float = 1e-5,
    ) -> None:
        super().__init__()
        self.window_size = window_size
        self.squared = squared
        self.eps = eps

    def _normalize_window(self, spatial_dims: int) -> tuple[int, ...]:
        if isinstance(self.window_size, int):
            window = (self.window_size,) * spatial_dims
        else:
            window = tuple(self.window_size)

        if len(window) != spatial_dims:
            raise ValueError(
                f"window_size must have {spatial_dims} values, got {window}"
            )
        if any(size < 1 or size % 2 == 0 for size in window):
            raise ValueError("window sizes must be positive odd numbers")
        return window

    def forward(
        self,
        fixed: torch.Tensor,
        warped: torch.Tensor,
    ) -> torch.Tensor:
        if fixed.shape != warped.shape:
            raise ValueError(
                f"shape mismatch: fixed={fixed.shape}, warped={warped.shape}"
            )
        if not fixed.is_floating_point() or not warped.is_floating_point():
            raise TypeError("LNCC inputs must be floating-point tensors")

        spatial_dims = fixed.ndim - 2
        if spatial_dims not in (2, 3):
            raise ValueError("expected [B,C,H,W] or [B,C,D,H,W] tensors")

        # Accumulate local statistics in float32 under mixed precision.
        if fixed.dtype in (torch.float16, torch.bfloat16):
            fixed = fixed.float()
            warped = warped.float()

        window = self._normalize_window(spatial_dims)
        padding = tuple(size // 2 for size in window)
        channels = fixed.shape[1]

        kernel = torch.ones(
            (channels, 1, *window),
            dtype=fixed.dtype,
            device=fixed.device,
        )
        conv = F.conv2d if spatial_dims == 2 else F.conv3d

        def local_sum(value: torch.Tensor) -> torch.Tensor:
            return conv(
                value,
                kernel,
                padding=padding,
                groups=channels,
            )

        count = local_sum(torch.ones_like(fixed))
        fixed_sum = local_sum(fixed)
        warped_sum = local_sum(warped)
        fixed_sq_sum = local_sum(fixed.square())
        warped_sq_sum = local_sum(warped.square())
        cross_sum = local_sum(fixed * warped)

        cross = cross_sum - fixed_sum * warped_sum / count
        fixed_var = (fixed_sq_sum - fixed_sum.square() / count).clamp_min(0)
        warped_var = (warped_sq_sum - warped_sum.square() / count).clamp_min(0)

        variance_product = fixed_var * warped_var
        valid = (fixed_var > self.eps) & (warped_var > self.eps)
        if not torch.any(valid):
            raise ValueError("no valid LNCC windows: local variance is too small")

        if self.squared:
            score = cross.square() / (variance_product + self.eps)
            score = score.clamp(0.0, 1.0)
        else:
            score = cross / torch.sqrt(variance_product + self.eps)
            score = score.clamp(-1.0, 1.0)

        return 1.0 - score[valid].mean()

调用时应传入已经完成空间变换的移动图像，而不是原始移动图像或位移场本身：

warped = spatial_transform(moving, displacement)
loss_sim = LNCCLoss(window_size=9)(fixed, warped)

对于三维各向异性数据，可以显式设置不同轴上的窗口大小：

loss_fn = LNCCLoss(window_size=(5, 9, 9))

7. 实现与实验中的常见问题

输入次序与变换方向

LNCC 本身对两个输入对称，但空间变换并不对称。必须先确认 $\phi$ 表示 fixed-to-moving 采样坐标还是 moving-to-fixed 位移，并确保生成的 $I_w$ 与 $I_f$ 位于同一网格。

低方差窗口

背景区域或均匀组织中的局部方差接近零。直接除法会放大数值误差；只增加 $\varepsilon$ 虽然能避免 NaN，却不一定能提供有意义的梯度。前景掩膜或有效窗口筛选通常更稳妥。

边界统计

使用零填充时，如果仍把窗口大小固定为 $k^d$，边缘的局部均值会被补零值拉低。上面的实现额外卷积全一张量，以获得边界窗口实际包含的样本数。

强度预处理

LNCC 对线性亮度变化较稳健，但并不免疫离群值、截断差异和非线性强度变化。训练前仍应统一重采样、强度裁剪和归一化策略。

多模态配准

LNCC 假设局部强度模式具有近似线性相关关系，因此更适合单模态任务。跨模态图像的对应组织可能没有稳定的线性灰度关系，此时通常需要互信息、模态无关描述子或学习式相似性度量。

8. LNCC 不能替代完整的配准评价

LNCC 衡量的是图像模式一致性。弱纹理区域、重复结构或具有相似灰度模式的错误位置，都可能得到较高相关性。因此，完整的医学图像配准评价还应覆盖至少两个层面：

• 解剖结构一致性：使用 Dice、目标配准误差（TRE）或表面距离衡量分割、标志点和边界是否对齐。
• 形变合理性：使用 Jacobian 行列式、折叠比例、位移或梯度统计检查形变是否连续、可逆并保持拓扑。

换句话说，LNCC 可以说明“局部图像看起来是否更一致”，却不能独立证明“解剖关系是否正确”或“形变是否合理”。训练损失与评价指标应围绕这三个问题分别设计，而不是把单一相似性分数当作配准质量的全部。

9. 小结

LNCC 可以看作在滑动窗口内重复计算 NCC，再对有效位置求平均。它通过局部去均值和方差归一化，减弱空间变化的亮度偏移与对比度差异，因而成为单模态医学图像配准中常见的相似性损失。

真正使用时，需要明确窗口的物理尺度、相关系数是否平方、损失的符号约定、低方差区域处理以及形变方向。更重要的是，LNCC 只描述图像层面的局部一致性；它应与解剖结构指标和形变场指标一起使用，才能形成完整的配准评价。