怀念图像配准先驱，他的算法让医学图像更准确

2023年12月28日，笔者的同事Allen Tannenbaum博士因病医治无效，在纽约长岛逝世，享年70岁。Allen患有癌症，在Covid之后，身体状况一直不稳定，近期急剧恶化。他在自动控制、计算机视觉和医学图像领域都做出了杰出贡献，尤其在医学图像领域，堪称是泰斗级的科学家。

Allen于1976年在哈佛大学获得数学博士学位，师从菲尔兹奖得主广中平佑（ Heisuke Hironaka ），研究领域为代数几何。毕业后，他一直在应用数学领域从事研究工作。Allen早期的工作主要是基于复分析中的Nevanlinna-Pick复插值理论来研究自动控制问题，他的主要突破是用椭圆曲线技巧证明了定义在多项式环上的系统中极点配置问题。后来，Allen应用偏微分方程理论来研究计算机视觉和医学图像问题，比较著名的工作包括他和Guillermo Sapiro发明的基于仿射不变热方程理论的图像增强算法；他和Steven Haker与Sigurd Angenent共同发明的基于最优传输的图像配准算法，等等。

笔者非常熟悉Allen的工作，曾经为Allen写过推荐信，帮助他转到石溪并且申请石溪计算机科学、应用数学和统计的杰出教授职位。笔者常常参加Allen的Seminar，与Allen深入探讨交流，也经常在自己的课上讲解Allen发明的各种算法。Allen的学术生涯非常成功，但是个人生活却并不顺利。Allen的儿子Emmanuel Tannenbaum患有脑瘤，Allen一直努力寻找各种方法来治疗，最终还是无法避免白发人送黑发人的悲剧。因此，笔者一直认为Allen为医学图像领域的研究倾尽了心血绝非为了名利，而是真正为了挽救生命。这里，笔者回顾一下Allen发明的基于最优传输理论的图像注册方法。

图像注册是医学图像领域中非常基本的问题之一，给定两张二维或者三维的CT或者MRI图像，和 ，寻找一个微分同胚，将图像中所有像素之间建立对应关系，准确反映人体器官的形变。Allen将源图像视为一个概率分布，每个像素的灰度值代表了密度函数。同样，目标图像视为目标概率分布，。不失一般性，我们假设密度函数绝对连续。两幅图像的总测度相等，即

图像注册的目的是找到一个微分同胚，使得对应像素的灰度值相近。这一点可以用保测度的性质来描述，即对于任意Borel集合，

保测度的性质记为。同时，根据生理学上的依据，每个原像点与像点之间的距离应当尽量小，即映射的整体几何畸变尽量小，因此我们寻找所有保测度映射中整体畸变最小者，这归结为所谓的蒙日问题（Monge's Problem）：

这里距离可以选为距离。蒙日问题的解被称为从到的最优传输映射。如果我们令，那么根据Brenier定理，最优传输映射存在且唯一，并且存在一个凸函数，其梯度映射等于最优传输映射。保测度条件的微分表示为

由此我们得到经典的蒙日-安培方程，

这里是凸函数的Hessian矩阵。边界条件为：

因此，图像注册问题转化为求解蒙日-安培方程问题。

蒙日-安培方程强烈非线性，求解蒙日-安培方程是非常具有挑战性的问题。在2001年，Allen, Haker和Angenent基于Brenier极分解定理，独树一帜地提出了一种新颖算法，令人耳目一新。极分解定理 给定映射，这里是中的凸集，是上的绝对连续测度，如果也是绝对连续，那么存在凸函数和映射，保测度，使得

并且和几乎处处唯一定义。即可以被唯一分解成两个映射的复合,分别是最优传输映射和保测度映射。

图1. 默比乌斯变换。

如图1所示，是平面圆盘，是圆盘到自身的默比乌斯变换,

这里圆盘上的测度为经典的勒贝格测度。根据极分解定理，可以被分解为最优传输映射和保勒贝格测度映射。在图2中，最优传输映射显示在左帧，每个棋盘格的面积发生变化，棋盘格的面积等于默比乌斯变换后相应棋盘格的面积。保勒贝格测度映射显示在右帧中，每个棋盘格的面积都等于初始面积。

图2. 映射的极分解。

Angenent-Haker-Tannenbaum（AHT）算法的思路如下：给定和，首先计算任意一个保测度映射，满足；其次，通过迭代，找到，，这样最优传输映射

初始保测度映射可以用流体力学方法求出。在我们的问题中为欧氏空间中的紧凸集。我们设计一个时变矢量场，每个粒子在流场中流动，以为速度场，粒子轨迹满足方程：

速度场与边界平行，即

这样流体不会溢出。如此得到时变微分同胚，

流体的密度为，其演化满足连续性方程，

我们设计密度为

同时假设矢量场，代入连续性方程，得到关于的泊松方程：

满足Neumann边界条件:

泊松方程是线性椭圆型PDE，可以用快速傅里叶变换求解，从而得到，流速场，积分得到微分同胚，初始保测度满足。

第二阶段，我们计算Brenier极分解中的保初始测度映射，。这一步计算基于矢量场的Hodge分解。

Hodge分解定理 给定定义在一个欧氏空间区域上的光滑矢量场，可以被唯一地分解为无旋场和无散场

这里是一个函数，为无散场，。Hodge分解也可以通过求解泊松方程求得：

由此得到

假设我们给定一个初始映射，我们构造矢量场

然后计算Hodge分解，得到无散场，构造流场

得到同胚,由连续性方程

由此保持，。我们更新，

这里是步长。重复迭代，直至无旋场的模足够小，算法终止。

图3. 基于最优传输的图像注册。

图3显示了AHT算法得到的大脑MRI图像注册结果，我们可以看到这一算法给出了微分同胚。图4显示了三维大脑MRI图像的注册结果，这证实了AHT算法可以直接向高维推广。

图4. 基于最优传输的3D图像注册。

Allen在2001年提出的求解最优传输问题的算法非常具有创意，迄今为止也是基于极分解理论的唯一算法，非常具有前瞻性。与传统的连续性方法求解蒙日-安培方程，或者基于几何变分法求解最优传输映射相比，AHT方法也是迭代法求解非线性蒙日安培方程，每一迭代步骤求解线性的椭圆型偏微分方程，这里是泊松方程。这种方法用流体力学方法，通过构造流场来计算微分同胚。各种限制条件都可以加入到流场中，因此在医学图像应用中非常灵活。

另一方面，目前AHT算法的理论证明依然缺失。在二维情形，AHT算法收敛于最优传输映射的理论证明已经完成；但是高维情形，这一算法的收敛性依然未知。主要困难在于，算法是收敛于全局最优还是局部最优，目前无法证明。另一方面，AHT算法基于Brenier极分解定理，因此要求的凸性，并且，同时对于解的正则性要求加高。我们知道，如果是凸区域，非凸，那么最优传输映射有可能在某个奇异集合上非连续。这种情形，AHT算法无法求解，而几何变分算法可以求解。

图5. 非连续的最优传输映射。

图5显示了具有奇异集合的最优传输映射，红色区域为凸集，蓝色区域非凸，源概率密度和目标概率密度都是勒贝格测度，右帧Brenier势能函数连续，但是非光滑，中存在奇异集合（黑色Graph），在奇异集合上Brenier势能函数连续不可导（红色曲线），最优传输映射在奇异集合上非连续（这种非连续性质可以解释AI生成模型中的模式坍塌现象）。

每次聚会，Allen都会讲很多笑话，见解独特，非常幽默风趣。但是私下里，Allen喜欢将办公室的门窗遮蔽起来，只留下一盏台灯，在黑暗中长时间的静默思考。每次Allen讲解算法或者理论，总能一针见血，言简意赅地指出问题关键。他的思想深邃独到，从不追波逐流，反而经常引领学术风潮。

Allen具有深厚的数学功底和敏锐深刻的洞察力，他能够将应用领域的基本问题提炼出数学本质，找到相应的现代数学理论，并且为抽象的理论发明出高效实用的计算方法，从而真正解决实际问题。Allen所发明的各种算法将会在工程和医疗领域中继续被发扬光大，造福人类。希望他的灵魂在天堂中得以安息。

本文转载自微信公众号“老顾谈几何”，原标题为《怀念Allen Tannenbaum博士》。

长按下方图片关注「返朴」，查看更多历史文章

微信实行乱序推送，常点“在看”，可防失联