顾险峰教授：解读医学影像配准的基本算法

本文来自微信公众号“老顾谈几何”，作者顾险峰教授，雷锋网 (公众号：雷锋网) 经顾险峰授权转载。

顾险峰教授，美国纽约州立大学石溪分校计算机系和应用数学系的终身教授，也是清华大学丘成桐数学科学中心访问教授。曾获美国国家自然科学基金CAREER奖，中国国家自然科学基金海外杰出青年奖（与胡事民教授合作），“华人菲尔茨奖”：晨兴应用数学金奖。

丘成桐先生和顾险峰博士团队，将微分几何，代数拓扑，黎曼面理论，偏微分方程与计算机科学相结合，创立跨领域学科“计算共形几何”，并广泛应用于计算机图形学，计算机视觉，几何建模，无线传感器网络，医学图像等领域。目前已经发表二百篇余篇国际论文，学术专著包括“Computational Conformal Geometry”（计算共形几何）， “Ricci Flow for Surface Registration and Shape Analysis”等。

以下是顾险峰教授著作的原文：

纽约长岛有着得天独厚的自然条件，这里气候温润，海岸线漫长。美丽的石溪大学坐落在长岛中部，坐拥一片私人海滩。这片海滩不为人知，空旷寂寥，深沉静谧。学校师生经常来这里沐浴海风，远眺大洋，特别在黄昏时分，在落日余晖中参省冥思，体悟自然真理。

1970年代初期，年轻的化学系助理教授Paul Lauterbur博士经常带着他的女儿来这里散步思考。有一天，三岁的小女孩在岸边捡到了一只微小的黑色蛤蜊，Paul为这只蛤蜊拍摄了人类历史上第一张活体断层图像。很可惜，短视的石溪大学拒绝为Paul的发明申请专利，专家们一致认为“这项发明可能带来的转让费不会弥补专利申请费。”最终，石溪失去了Lauterbur博士，也失去了经济上腾飞的一次机遇。历史是公正的，因为发明了核磁共振断层扫描技术，Paul Lauterbur博士于2003年获得诺贝尔医学奖。

数十年来，核磁共振技术和CT断层扫描技术彻底地革命了医学，医学影像技术使得医生可以直接看到病人体内，从而精准地进行诊断，制定治疗方案，检验治疗效果。很多病变都会诱发器官组织的变形，或者由器官变形所诱发，例如大脑皮层的萎缩退化诱导老年失智，各种肿瘤会在器官表面形成凸起，骨质流失会引起骨骼的变形。因此，医生可以通过精确比对器官的几何形状，来判断脏器是否反常；通过分析肿瘤的几何特征，来判断肿瘤的良性恶性。这些可以归结为医学影像的配准（registration）和分析（analysis）问题。

医学影像是一个非常庞大的学科，有种类繁多的计算方法。图像配准的方法也是异常丰富，理论严谨，富有实效。从基础的方法论角度而言，比较普适的主要有基于微分几何的方法，基于流体力学的方法和基于概率测度理论的方法。虽然，这些途径中问题提法、数学 工具 、计算机算法非常不同，但是最终都可以归结为几何偏微分方程，进而转换为变分优化问题。这三种方法各有优缺点，相辅相成，都有很大的实用价值。每种方法都在迅猛发展，分支众多。下面，我们简略讨论这三种方法的核心思想和最简单的算法。每一种方法深入下去，都是一片浩瀚的海洋。

微分几何方法

首先，我们从医学图像中提取我们感兴趣的器官，表示成二维曲面或者三维实体。例如，我们希望研究大脑皮层曲面。首先得到整个颅部的核磁共振断层扫描图像；然后对每个断层的图像进行分割（segmentation），将不同的组织进行区分，将骨骼、脑灰质、脑白质分开；其次，提取我们感兴趣的组织轮廓（contour），主要是脑灰质；然后将每个断层图像中的轮廓线摞在一起，组合成封闭曲面。当然，实际算法远比理想描述复杂得多，每一步都会引入大量的几何、拓扑误差，需要精细而严密的方法以提高稳定性和精度。

图1. 从医学图像中重建的大脑皮层曲面（王雅琳）

例如，我们比较同一个人不同时期扫描得到的大脑皮层曲面，来定量分析大脑萎缩情况。从根本层面上而言，我们希望将第一个大脑皮层上的每一个细胞映到第二个大脑皮层上相同的细胞位置。当然，在目前技术现实中，这是无法真正做到的，我们可以通过设计各种能量，加上一些限制条件来逼近理想情形。大脑皮层上面有解剖学特征，例如沟回皱褶，我们要求微分同胚将主要的沟回映射到相应的沟回上面。同时，我们要求曲面映射的几何畸变最小。如果我们假设大脑皮层曲面是具有弹性的，那么我们要求映射带来的弹性形变势能最小。弹性势能一般用调和能量来表述。根据不同的医学应用，我们可以设计不同的能量。假设我们已经得到同一器官的不同曲面，或者两个三维实体，注册问题的微分几何提法如下：给定三维欧氏空间中的两个嵌入流形， ，给定初始映射  ，求两个流形间的微分同胚 ，这个微分同胚和初始微分同胚同伦，，满足一定的限制条件，例如将特征点映到相应的特征点，；同时优化某种能量，例如调和能量。形式上，我们可以把注册问题表示成：

从理论层面来讲，我们需要明确知晓这一问题解的存在性、唯一性、稳定性和正则性。从算法角度而言，我们需要求出能量的一阶变分和二阶变分。理论上的困难有很多，通常情况下，我们考虑的泛函空间-流形间的微分同胚空间 本身并没有紧性，即便我们找到了能量极小化序列，其极限不见得是微分同胚，因此解的存在性证明比较困难；如果能量没有凸性，我们很难保证解的唯一性，能量的凸性很多时候依赖于流形本身的几何特性，例如高斯曲率；各种直观的几何限制条件，例如特征点的对应，很难转换成解析形式，往往用偏微分方程来表述，或者用变分的局部线性逼近来显示表达。算法角度的困难也很多，计算机所能处理的计算都是离散的，如何将微积分表述的概念来离散逼近，如何加速收敛，如何用单纯复形的组合结构（离散多面体）来描述流形结构，如何寻找适合数值偏微分方程的三角剖分或者样条表示，这些问题也都具有本质困难。

基于微分几何的注册方法也有多种，比较常用的是基于共形几何（conformal geometry）的内蕴方法，这种方法只依赖曲面的黎曼度量，不需要曲面的嵌入。如图2所示，其核心方法如下：首先，我们用单值化定理（uniformization）将曲面映射到标准空间中，；其次，我们在标准空间中寻找最优微分同胚， ；最后，曲面间的微分同胚由复合映射给出， 。

图2. 基于共形几何的曲面配准方法

第一步归结为寻找一个函数 ，变换曲面本身的黎曼度量，使得新的黎曼度量诱导常值高斯曲率。根据曲面的拓扑，高斯曲率有三种选择+1,0，-1，对应的标准空间为单位球面，欧氏平面和双曲平面，如图3所示。曲面单值化可以用曲面黎奇流方法得到。这样，我们将三维空间中曲面间的映射归结为标准空间的自同胚，减小了搜索范围，降低了计算难度。

图3. 曲面单值化

第二步，我们首先将标准空间的自同胚用拟共形映射（quasi-conformal map）来表述。所有的K-拟共形映射构成的空间具有紧性，可以保证解的存在性。每个自映射可以用Beltrami微分来表示，因此我们在Beltrami微分空间中对相应能量进行变分，构造能量极小序列，从而求得极限。这里问题的关键在于如何将限制条件用Beltrami微分来表示，和如何求得欧拉-拉格朗日方程。如果我们要求解是同胚，则Beltrami 微分的模小于1，如果我们要求解保持共形结构，则Beltrami微分和源曲面上的全纯二次微分作用为0，等等。很多时候，解的唯一性和正则性依赖于目标曲面上的高斯曲率。例如，曲面间度为1的调和映射，如果目标曲面上的高斯曲率处处为负，则必为微分同胚。再如，曲面间映射的任意同伦类中存在唯一的Teichmuller映射，使得曲面共形结构畸变最小。调和映照和Teichmuller映射在医学图像配准问题中被经常使用。在求解这些映射过程中，曲面的拓扑、黎曼度量和共形结构经常被灵活变化，从而达到最优的目的。

图4. 腰椎骨质流失测量，共形几何方法（雷诺明）

流场方法

大形变微分同胚度量映射（LDDMM）方法考虑空间中的流场，每个粒子在空中流动，每一刹那的流场用粒子的速度向量场来描述。固定时刻，粒子从起始位置到达终点位置，这给出了空间到自身的微分同胚。

流场的速度场可以被视作是有限支集的光滑映射，所有的流速场记为 。给定时变流速场 ， ，给出了空间中每一点每一刻运动的速度向量，我们得到依赖于时间的微分同胚满足微分方程：

如果我们固定一点， ，给出了此点的运动轨迹。我们定义矢量场的模为

欧氏空间的微分同胚变换群记为：

微分同胚群中的右不变距离为：

可以证明，是的一个微分同胚群，在度量下是完备的，任意两个微分同胚 之间存在一条极短测地线。

给定两幅医学图像，记为和，我们寻找一个微分同胚，使得图像和图像尽量接近，同时和恒同映射尽量接近。由此，我们得到能量：

,

等价地，我们可以用流速场来表示同样的能量：

,

图像注册问题归结为极小化能量，求得流速场。

我们可以用梯度下降法来优化能量。微分同胚对路径非线性依赖，但是的变分可以表示成：

,

这里 。由此，我们可以得到能量关于流速场的梯度，从而可以进行优化。

图5.孕期婴儿大脑皮层36周到43周的生长情况（Laurent Risser， LDDMM 方法）

概率测度方法

从根本层面上而言，CT图像和MRI图像反映的是人体内部组织的密度。由此，我们可以将医学图像的灰度值看成是概率测度。两个图像之间的匹配应该满足如下条件：密度相近的点彼此对应，每个原像点和其像点的距离尽量接近。从这个角度出发，人们提出了基于最优传输理论（optimal transportation theory）的图像匹配方法。

例如我们匹配CT体图像，设是三维欧氏空间中的凸区域，（例如单位立方体），两幅图像代表定义其上的概率测度， ，满足总测度相同，

,

我们寻找自同胚 ，将推前到，这意味着对于一切波莱尔集,

,

记成，对应的微分方程为雅克比方程：

。

雅克比方程的解有无穷多个。每个映射的传输代价定义为：

我们寻找最优传输映射，极小化传输代价， 。当p为2时，最优传输映射是某个凸函数的梯度映射， ，。代入雅克比方程，

我们的蒙日-安培方程：

.

图6. 大脑白质图像匹配（Allen Tannenbaum，最优传输方法）

小结

这三类方法基于不同的物理和几何的理解，所用的理论工具相距甚远，风格迥异，各有千秋。

基于微分几何的方法有着强烈的几何直观，只用黎曼度量，不需要流形在欧氏空间的嵌入，洗练简洁，严密普适，计算效率较高。但是需要比较抽象的黎曼几何概念，需要先从图像中提取器官曲面。基于流场的方法直接灵活，适用于各种数据类型，但是需要流形的嵌入，形变过程中曲面不允许自交。同时为了二维曲面匹配，需要计算三维背景空间的自同胚，计算量较大。基于最优传输的方法兼顾物理直觉，可向高维直接推广，但是所解方程高度非线性。另一方面，这几种方法具有内在的密切联系：最优传输的蒙日-安培方程和微分几何中的Alexandrov定理等价；最优传输存在流体力学的解释。总之，所有方法都是在微分同胚群中进行优化，问题具有本质的难度。

雷锋网版权文章，未经授权禁止转载。详情见 转载须知 。