图像配准是什么
图像配准是医学图像领域一个常见的任务，其目标是把待配准图像$m$ (moving image)通过一个变换$T$对应到参考图像$f$(fixed image)上
Sotiras A, Davatzikos C, Paragios N. Deformable medical image registration: A survey. IEEE transactions on medical imaging, 2013, 32(7): 1153-1190.
一个比较直观的示意图如下：
Hoffmann M, Billot B, Greve D N, et al. SynthMorph: learning contrast-invariant registration without acquired images. IEEE transactions on medical imaging, 2021, 41(3): 543-558.
https://colab.research.google.com/drive/1zaDnAJGUokS0knqWttuTgrRJMb6zxukI#scrollTo=V8PIMBl0Idsk
主要用于运动矫正，影响组学分析等领域
在实际应用中判断配准质量一个难点是，$m$和$f$可能属于不同的模态，例如PET/CT/MRI之间的跨模态配准，或MRI之间不同模态（T1/T2/DWI）之间的配准。不同模态采集图像的对比度可能有较大的变化，例如心脏T1成像中不同翻转时间的图像对比度会有较大变化
因此无法通过简单的MSE来判断配准质量。通常需要借助分割标签来辅助判断，理想情况下配准后图像和$f$的同一个分割区域应该是完全重叠的，可以通过分割区域的重叠度来判断配准质量。
DSC(Dice Similarity Coefficient)
说到分割标签，一个常用的Metric是Dice系数：
$$
\mathrm{Dice}(A,B)=\frac{2|A\cap B|}{|A|+|B|}
$$
https://mp.ofweek.com/ai/a956714557177
当两个区域完全重叠时，Dice系数为1，完全无重叠时为0。
Python实现：
import numpy as np
def dice_coefficient(pred: np.ndarray, true: np.ndarray, smooth: float = 1e-5) -> float:
"""
Dice Similarity Coefficient (DSC) between two binary masks.
pred, true : binary arrays of the same shape (any dimensionality)
smooth : avoids division by zero when both masks are empty
"""
pred = pred.astype(bool)
true = true.astype(bool)
intersection = np.sum(pred & true)
denom = np.sum(pred) + np.sum(true)
return (2.0 * intersection + smooth) / (denom + smooth)
Dice系数的问题在于它的本质是面积重叠程度，对于局部边界不敏感，例如对于心肌来说，即使心肌边界有明显偏移，由于其他部分重叠较多，Dice系数可能看起来也还不错。因此需要一个另外的指标来查看边界的偏移程度。
HD(Hausdorf Distance)
豪斯多夫距离用来衡量最坏情况下的边界误差，首先定义单项距离： $$ d(A,B)=\max_{a\in A}\min_{b\in B}\|a-b\| $$ 即A中每个点到B的最短距离，取这些最短距离中的最大值
看起来有点烧脑，我们看示意图，其实就是区域A上离B最远的一点到B的最近距离。
Jungeblut P, Kleist L, Miltzow T. The complexity of the Hausdorff distance. Discrete & Computational Geometry, 2024, 71(1): 177-213.
HD是双向的，取两个方向中比较大的一个： $$ HD(A,B)=\max\left(\max_{a\in A}\min_{b\in B}\|a-b\|,\max_{b\in B}\min_{a\in A}\|b-a\|\right) $$ Python实现：
import numpy as np
from scipy.spatial.distance import cdist
from scipy.ndimage import binary_erosion
def _surface_points(mask: np.ndarray) -> np.ndarray:
"""Coordinates of the boundary pixels/voxels of a binary mask."""
mask = mask.astype(bool)
eroded = binary_erosion(mask)
boundary = mask & ~eroded
return np.argwhere(boundary)
def hausdorff_distance(pred: np.ndarray, true: np.ndarray, voxel_spacing: float = 1.0) -> float:
"""
Symmetric Hausdorff Distance between the boundaries of two binary masks.
Returns NaN if either mask has no foreground (empty boundary).
"""
pred_pts, true_pts = _surface_points(pred), _surface_points(true)
if len(pred_pts) == 0 or len(true_pts) == 0:
return np.nan
d = cdist(pred_pts, true_pts) * voxel_spacing
return max(d.min(axis=1).max(), d.min(axis=0).max())
这里的一个小优化是，把两个区域中所有点都计算一遍的计算量过大，而且也没有意义，因为最坏的情况只会出现在边界上，所以先用_surface_points取出区域边界，再计算两个边界上的点两两之间的最近距离，再取最大值
另外由于HD对边界敏感，容易被极值点带偏，因此在使用中通常使用HD95，也就是取95百分位数，把可能出现的异常点过滤掉，实现如下
def hausdorff_distance_95(pred: np.ndarray, true: np.ndarray, voxel_spacing: float = 1.0) -> float:
"""
95th-percentile Hausdorff Distance (HD95) — more robust to a single
outlier boundary point than the raw HD above.
"""
pred_pts, true_pts = _surface_points(pred), _surface_points(true)
if len(pred_pts) == 0 or len(true_pts) == 0:
return np.nan
d = cdist(pred_pts, true_pts) * voxel_spacing
return max(np.percentile(d.min(axis=1), 95), np.percentile(d.min(axis=0), 95))