球面最优传输映射的C++实现

曲面的微分几何理论断言给定曲面的第一基本形式（黎曼度量）和第二基本形式，曲面在三维欧氏空间中的嵌入就被确定（相差一个刚体变换）。对于封闭的凸曲面而言，黎曼度量就可以确定曲面的形状。但是如何由黎曼度量来计算凸曲面形状，一直以来没有成熟的算法。为此，我们安排了一次作业，实现从高斯曲率来重构凸曲面的算法。

图1. 由高斯曲率来重构凸曲面。

Minkowski问题

假如我们在球面上定义了高斯曲率函数，则相应的黎曼度量与球面度量之间相差一个共形因子函数，满足经典的Yamabe方程：，这里是球面度量决定的Laplace-Beltrami算子。共形因子可以通过曲面的Ricci流（计算共形几何课程第8次作业）求得，

由此，我们可以计算出来给定高斯曲率相应的黎曼度量。但是这里黎曼度量抽象给出，并没有具体的曲面在三维空间中的嵌入信息。如果我们希望进一步得到曲面的嵌入，则需要求解闵可夫斯基（Minkowski）问题。Minkowski问题提法如下：给定一个凸区域，其边界为凸曲面，其极坐标表示为，

即曲面上任意一点表示为极径向量与极径长度的乘积。高斯曲率函数为定义在球面上的正值函数，在点处的高斯曲率为。简而言之，Minkowski问题就是从给定的来反求。我们定义法向映射（normal mapping），对于凸曲面上的任意一点,

高斯映射（Gauss Mapping）定义为, 高斯映射诱导了球面上的高斯测度（Gauss measure），对于任意Borel集合,

这里H是Hausdorff测度。Minkowski问题：给定单位球面上的Borel测度，求凸区域，使得。Minkowski证明了解的存在性和唯一性。

球面最优传输

图2：广义勒让德变换。

考察凸曲面的一个支撑平面，法向量为，高度为，那么高度满足等式：

,

反之，凸曲面可以视为其所有支撑平面的内包络，

,

我们定义，，传输代价函数为

.

由此，我们定义广义Legendre变换（c-transform），

,

则我们有，并且。那么Minkowski问题可以被转化为最优传输问题：给定球面的Hausdorff测度和高斯测度，求一个球面到自身的传输映射，将Hausdorff测度映射到高斯测度，同时在所有保测度的映射中，使得传输代价最小，

.

这里传输映射由凸曲面的高斯映射给出。由Kantorovich理论，这等价于优化下面的泛函：

由最优传输理论，我们可以得到解的存在性和唯一性。

离散球面最优传输映射

图3. 离散最优传输映射。

在离散情形，我们给定 单位球面上的离散点集，待求径向长度记为，凸曲面是凸包

可以证明，每个顶点对偶的平面为，

凸曲面的Legendre对偶是支撑平面的包络，

凸包的法向映射诱导了球面的一个胞腔分解，这个胞腔分解等价于Legendre对偶向球心的投影，

图4：球面Power Voronoi胞腔分解。这个胞腔分解等价于球面上的Power Voronoi Diagram，其中球面上的Power距离由传输代价所定义。每个胞腔的边界都是球面测地线。

图5：球面最优传输映射能量的Hessian矩阵。

算法的核心是优化下面的凸能量，

能量的梯度为

能量的Hessian矩阵为

Hessian矩阵对角线元素为：

我们可以用牛顿法进行优化，但是在优化过程中应该时刻注意每一步优化都在可容许解空间中进行,

。

算法设计

离散凸曲面用三角网格表示，存储为半边数据结构，并且也用于表示Legendre对偶曲面。凸包算法可以用传统的增量算法（incremental convex hull），或者Lawson Edge Flip算法。实用中，增量算法比较稳定，但是计算费时，边对换算法比较省时，但是数值稳定性较为逊色。

图6. 输入三角网格，普朗克头像。

输入是三角网格，亏格为零的封闭曲面。数据格式是通用的三维数据格式，例如obj，mesh文件。

图7. 球面调和映射的像。

我们先用球面调和映射将三角网格映射到单位球面上，然后每个顶点在调和映射下的像作为球面采样点，每个顶点相邻面的面积和作为此采样点对应的目标测度。

图8. 球面胞腔分解。

球面的Delaunay三角剖分和胞腔分解由同一个三角网格来表示，三角剖分的顶点、边和面对偶于胞腔分解的面、边和顶点。

图9. Minkowski问题的解，由曲率决定的凸曲面。

优化用牛顿法实现，为了保证优化中间过程都在可容许空间中进行，我们采用阻尼法。在优化过程的每一次循环，如果极半径更新后，球面胞腔分解的某个胞腔消失，则我们将步长减半，重新计算极半径，如此循环，直至临时解诱导的胞腔分解满足每个胞腔都非空的条件。

图10. Legendre对偶曲面。

如此我们得到的凸包曲面是Minkowski问题的解，顶点的离散高斯曲率等于给定的目标测度。计算结果既包括Minkowski凸曲面，又包括Legendre对偶曲面，这一对曲面相互对偶，彼此相互决定。

图11. 大脑皮层曲面。图12. 大脑皮层曲面的球面调和映射。图13. 球面上的测地胞腔分解。图14. Minkowski问题的解。图15. Minkowski凸曲面的Legendre对偶曲面。

软件使用和下载

我们将基于几何变分的球面最优传输映射算法设计成课程的一次作业，并且将源码公布，以起到抛砖引玉的作用。算法用通用C++语言实现，可以在Windows，Linux平台上编译运行，图形界面基于OpenGL，线性方程组求解基于Eigen库，底层几何计算基于Shewchuk的自适应算术库。演示程序、源代码、测试数据都可以从【1】下载，或者扫描下方二维码：

解压演示程序后，双击bat文件可以自动演示。程序的源码用CMake可以生成项目文件，再用Visual Studio等编译器进行编译。我们也提供了测试数据，格式为通用的obj文件。

算法的理论基础来自于丘成桐先生和郑绍远教授关于蒙日-安培方程的理论，我们提出的几何变分原理【4】，也参考了汪徐家教授的球面最优传输理论【5】，提出了基于几何变分的球面最优传输映射算法【3】，得到斯杭博士、汪徐家教授、罗锋教授、崔丽教授、苏科华教授的建议和帮助，雷娜教授进行具体指导。开发团队主要成员包括齐鑫、陈伟、安东生、郭洋、赵彤、温成峰、章敏、曾薇等博士、博士后。

未来，我们会进一步公布和开发更加通用计算工具，例如高维和一般密度函数的最优传输映射，同时也欢迎更多的同学加入到开发团队之中。希望大家发现算法漏洞，提出改进方法，更加欢迎深入的学术探讨，宝贵的建议和批评！

---

【1】https://www3.cs.stonybrook.edu/~gu/software/SOT/index.html

【2】Na Lei, Dongsheng An, Yang Guo, Kehua Su, Shixia Liu, Zhongxuan Luo, Shing-Tung Yau and Xianfeng Gu, "Geometric Understanding of Deep Learning", Journal of Engineering, 2020.

【3】Li Cui, Xin Qi, Chengfeng Wen, Na Lei, Xinyuan Li, Min Zhang and Xianfeng Gu, Spherical Optimal Transportation", Computer-Aided Design (CAD), Volume 119, Pages 181-193, 2019.

【4】Xianfeng Gu, Feng Luo, Jian Sun and Shing-Tung Yau, "Variational Principles for Minkowski Type Problems, Discrete Optimal Transport, and Discrete Monge-Ampere Equations",  AJM 20(2), 383-398, 2016.

【5】Xu-Jia Wang, "On the design of a reflector antenna II",  Calc Var Partial Differential Equations, 20(3), pp. 329-341, 2004.

---

请长按下方二维码，选择 “识别图中二维码”，即可关注。

【老顾谈几何】邀请国内国际著名纯粹数学家，应用数学家，理论物理学家和计算机科学家，讲授现代拓扑和几何的理论，算法和应用。