[论文解读] A numerical algorithm for $L_2$ semi-discrete optimal transport in 3D

Bruno Lévy|arXiv (Cornell University)|Sep 3, 2014

Markov Chains and Monte Carlo Methods参考文献 16被引用 16

一句话总结

本文提出了一种数值算法，用于计算具有分段线性密度的三维四面体网格与一组狄拉克质量之间的 $L_2$ 半离散最优传输映射。该方法利用重心图（power diagrams）和拟牛顿法最小化一个凸目标函数，通过一种高效的几何算法计算重心图与网格之间的交集，实现在包含多达一百万个狄拉克质量的网格上，计算时间在45分钟内完成，具有鲁棒性和可扩展性。

ABSTRACT

This paper introduces a numerical algorithm to compute the $L_2$ optimal transport map between two measures $μ$ and $ν$, where $μ$ derives from a density $ρ$ defined as a piecewise linear function (supported by a tetrahedral mesh), and where $ν$ is a sum of Dirac masses. I first give an elementary presentation of some known results on optimal transport and then observe a relation with another problem (optimal sampling). This relation gives simple arguments to study the objective functions that characterize both problems. I then propose a practical algorithm to compute the optimal transport map between a piecewise linear density and a sum of Dirac masses in 3D. In this semi-discrete setting, Aurenhammer et.al [\emph{8th Symposium on Computational Geometry conf. proc.}, ACM (1992)] showed that the optimal transport map is determined by the weights of a power diagram. The optimal weights are computed by minimizing a convex objective function with a quasi-Newton method. To evaluate the value and gradient of this objective function, I propose an efficient and robust algorithm, that computes at each iteration the intersection between a power diagram and the tetrahedral mesh that defines the measure $μ$. The numerical algorithm is experimented and evaluated on several datasets, with up to hundred thousands tetrahedra and one million Dirac masses.

研究动机与目标

• 开发一种实用且可扩展的算法，用于计算具有分段线性密度的三维四面体网格与由狄拉克质量组成的离散测度之间的 $L_2$ 最优传输映射。
• 通过计算自然且视觉上令人愉悦的传输映射，支持计算机图形学中的应用，如三维渐变与形状变形。
• 解决在三维空间中高几何保真度下准确计算最优传输映射的挑战，尤其是在拓扑变化发生时。
• 提供一个鲁棒且高效的计算流水线，集成几何判定、重心图构建与多级优化。

提出的方法

• 通过确定最小化凸目标函数的重心图权重来计算最优传输映射，该方法在半离散设置下由 Aurenhammer 等人建立。
• 采用拟牛顿优化方法最小化目标函数，梯度计算基于重心图与四面体网格之间的几何交集。
• 设计一种高效且鲁棒的算法，同时遍历四面体网格与重心图以计算其交集，采用浮点数与任意精度算术，结合算术过滤与符号扰动。
• 集成多级方法，结合 BRIO 预排序与二阶回归，以加速收敛并提升数值稳定性。
• 几何判定与重心图构建在公开的 geogram 库的 PCK 模块中实现，以确保可重现性与可重用性。
• 通过检测 Delaunay 与常规加权三角剖分之间的交集，实现对变形过程中传输映射不连续性的检测，从而支持拓扑变化。

实验结果

研究问题

• RQ1 当源测度为四面体网格上的分段线性密度时，如何高效且鲁棒地求解三维 $L_2$ 半离散最优传输问题？
• RQ2 半离散最优传输中的目标函数结构是什么？其梯度如何通过几何交集高效计算？
• RQ3 在变形过程中，尤其是发生拓扑变化时，最优传输映射中的不连续性能否被准确检测与表示？
• RQ4 当处理包含多达一百万个狄拉克质量的大规模问题时，该算法在计算时间与内存消耗方面的可扩展性如何？

主要发现

• 该算法成功计算了具有分段线性密度的三维四面体网格与最多一百万个狄拉克质量之间的最优传输映射，对于一百万个质量点的计算时间约为45分钟。
• 该方法实现了稳定且视觉上令人愉悦的三维渐变效果，如将犰狳网格变形为球体，即使在显著的拓扑变化下也表现良好。
• 采用结合 BRIO 预排序与二阶回归的多级算法显著提升了收敛速度与数值鲁棒性。
• 采用任意精度算术与过滤技术实现的几何交集算法，确保了在退化构型与近退化情况下的鲁棒性。
• 尽管具有鲁棒性，但由于不连续性采样限制，传输映射中仍存在约一个四面体宽度的微小间隙，尤其在高变形区域较为明显。
• 该算法表现出良好的可扩展性，计算时间从1,000个质量点的1.45秒增长至一百万个质量点的2,649秒（约44分钟），在犰狳到球体的基准测试中表现稳定。