球面最优传输映射的C++实现
曲面的微分几何理论断言给定曲面的第一基本形式（黎曼度量）和第二基本形式，曲面在三维欧氏空间中的嵌入就被确定（相差一个刚体变换）。对于封闭的凸曲面而言，黎曼度量就可以确定曲面的形状。但是如何由黎曼度量来计算凸曲面形状，一直以来没有成熟的算法。为此，我们安排了一次作业，实现从高斯曲率来重构凸曲面的算法。
图1. 由高斯曲率来重构凸曲面。
由此，我们可以计算出来给定高斯曲率相应的黎曼度量。但是这里黎曼度量抽象给出，并没有具体的曲面在三维空间中的嵌入信息。如果我们希望进一步得到曲面的嵌入，则需要求解闵可夫斯基（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对偶向球心的投影，
图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对偶曲面。
【1】https://www3.cs.stonybrook.edu/~gu/software/SOT/index.html
请长按下方二维码，选择 “识别图中二维码”，即可关注。
【老顾谈几何】邀请国内国际著名纯粹数学家，应用数学家，理论物理学家和计算机科学家，讲授现代拓扑和几何的理论，算法和应用。