CN108022005A

CN108022005A - 一种高效的全离散最优传输方法

Info

Publication number: CN108022005A
Application number: CN201711195311.8A
Authority: CN
Inventors: 苏科华; 陈彩玲; 焦冲; 顾险峰
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2017-11-24
Filing date: 2017-11-24
Publication date: 2018-05-11

Abstract

本发明公开了一种高效的全离散最优传输方法，首先将源区域和概率测度离散化表示，输入目标离散点，给每个目标离散点赋予一个狄拉克测度；初始化截距向量，将源区域离散点的坐标和截距向量的系数用矩阵存储，将目标离散点的坐标及截距向量用矩阵存储；每次迭代过程中，通过矩阵相乘、根据行最大值将离散点分类的方法构造近似离散加权维诺图；求全离散最优传输的能量函数及其梯度；更新截距向量；循环直到传输映射的凸能量达到极小值，得到最终的全离散最优传输的解。本发明简单高效，易于实现，可通过多线程或CUDA等并行计算技术进一步加速提高计算效率，可用于求解大规模最优传输问题。本发明与源区域的维度无关，适用于求解任意维度的最优传输问题。

Description

一种高效的全离散最优传输方法

技术领域

本发明属于计算机图形图像处理技术领域，涉及一种高效的全离散最优传输新方法。

背景技术

法国数学家Monge最早提出了著名的最优传输问题，即设计一个将某质量分布(X,μ)运输到另一质量分布(Y,ν)的最佳方案，满足质量相等∫_Xμ(x)dx＝∫_Yν(y)dy，同时使得运输代价最小。其中X,Y是两个有界区域，μ,ν表示密度函数。1948年，俄罗斯数学家Kantorovich运用线性规划思想成功求解最优传输问题，并证明了其解的存在性和唯一性，由此获得了1975年的诺贝尔经济奖。如今，Monge-Kantorovich理论已取得大量的研究成果，极大地推动了数学、经济学、医学、人工智能、计算机科学、医学等学科的发展。对于纯粹数学，最优传输理论可用于求解曲面等距嵌入问题等凸几何问题，如Alexandrov问题、Minkowski问题；对于计算机数学，最优传输算法可用于求解非线性的椭圆型偏微分方程，如蒙日-安培方程；对于经济学，最优传输可以应用于生产消耗运输、资源配置与调度等；对于人工智能，Wasserstein距离给出了衡量两概率测度分布距离的方法，在WGAN模型等深度学习领域发挥了显著作用；此外，最优传输理论在计算机视觉、三维模型参数化、曲面注册、医学成像等领域发挥着重大作用，取得很好的效果。

对Monge-Kantorovich的研究主要分为三类：第一类是全离散最优传输，即源质量分布(X,μ)和目标质量分布(Y,ν)都是离散的；第二类是半离散最优传输，即源质量分布(X,μ)是连续的，而目标质量分布(Y,ν)是离散的；第三类是连续最优传输，即源质量分布(X,μ)和目标质量分布(Y,ν)都是连续的。

最优传输问题的求解方法，归纳起来主要有三种：线性规划方法、熵正则化距离近似方法以及基于计算几何的凸优化法。Kantorovich最早采用线性规划思想求解最优传输问题。线性规划的思想是将空间区域X和Y表示成离散点集，将概率测度μ,ν离散成狄拉克测度，进而将最优传输问题转换成线性规划问题来求解。最优传输问题的等效解决方案是Earth Mover’s Distance，简称EMD，也称Wasserstein距离，用于测量两概率测度分布之间的距离。以测地距离为基础内核，并用热核来进行近似，其迭代数值解方案具有线性收敛性，每次迭代只需要求解高斯卷积。与线性规划法相比，该方法极大地提高了计算性能。但由于采用了近似方法简化问题，其解并不是真正意义上的最优传输。法国数学家与1991年Brenier证明，当传输代价为L²距离且传输区域为凸区域时，某个凸函数的梯度即为最优传输映射，由此建立了最优传输与凸几何的联系。近年来，顾险峰对Brenier理论进行离散化，基于变分法通过动态构造加权维诺图(power diagram)以及采用凸优化来求解最优传输问题。

上述求解最优传输问题的方法各具优缺点，线性规划原理简单，但它的空间复杂度高，计算效率低，无法用于数据量大的大规模最优传输问题中。熵正则化距离近似法计算性能有所提高，但其解并不是真正意义上的最优传输。基于计算几何的凸优化方法计算复杂度低，算法性能高，但只适合求解半离散最优传输，不适用全离散最优传输问题的求解。根据顾险峰等学者的研究理论，求解最优传输问题的核心是动态构造加权维诺图，而加权维诺图的构造等价于求解线性函数的上包络。但目前主流的最优传输算法中，源区域都是连续的，加权维诺图的构造较为复杂,每个胞腔的测度、能量函数及其梯度函数的计算繁琐复杂。

发明内容

为了解决上述技术问题，本发明提供了一种高效的全离散最优传输方法。

本发明所采用的技术方案是：一种高效的全离散最优传输方法，其特征在于，包括以下步骤：

步骤1：将源区域和概率测度离散化表示，输入目标离散点，给每个目标离散点赋予一个狄拉克测度；

步骤2：初始化截距向量，将源区域离散点的坐标和截距向量的系数用矩阵存储，将目标离散点的坐标及截距向量用矩阵存储；

步骤3：每次迭代过程中，通过矩阵相乘、根据行最大值将离散点分类的方法构造近似离散加权维诺图；

步骤4：求全离散最优传输的凸能量函数及其梯度；

步骤5：更新截距向量；

步骤6：循环步骤3-步骤5直到传输映射的凸能量函数达到极小值，得到最终的全离散最优传输的结果。

作为优选，步骤1中，源区域和目标区域都离散化表示，而且它们的总测度相等，即：

其中，N表示源区域的离散顶点总数，μ_i表示点x_i上的狄拉克测度；M表示目标离散点个数，ν_j为目标离散点y_j上赋予的狄拉克测度。

作为优选，步骤2中，将源离散点的信息用大小为N×(D+1)的矩阵S存储，其中N表示源区域的离散顶点总数，D表示源离散点的维度，矩阵的前D列存储源离散点的坐标信息，第D+1列存储截距向量的系数，值为1；

将截距向量h初始化为零向量，即h＝0，或者初始化为(0,1)之间的随机值，或者初始化为每个目标离散点y_j到原点距离平方的一半，即其中<,>表示点积；

将目标离散点的坐标信息和截距向量存于一个大小为M×(D+1)的矩阵T中，其中前D列为目标离散点的坐标信息，第D+1列为截距向量；

其中，(a_i,b_i,···)表示源区域离散点x_i的顶点坐标；(c_j,d_j,···)表示目标区域离散点y_j的顶点坐标，h_j为该目标离散点的权重。

作为优选，步骤3中，根据加权维诺图的定义：

推导出近似的离散加权维诺图：

因此，求解近似的离散加权维诺图的具体方法如下：

每次迭代中，矩阵S与T的转置矩阵T'相乘得到大小为N×M的矩阵P，即：

矩阵P表示一系列线性函数，其元素p_ij为p_ij＝<x_i,y_j>+h_j，表示源离散点x_i与目标离散点y_j的点乘与截距向量之和。根据推导的近似离散加权维诺图定义：可知，每次迭代过程中，矩阵P中所有线性函数的上包络即为所要求的加权维诺图。因此可依据矩阵P每一行元素最大值所在列将源离散点x_i进行分类。如矩阵第i行的最大值在第j列，则表示与离散点x_i相关的所有线性函数的上包络位于目标离散点y_j处，即源离散点x_i位于胞腔Pow(y_j)内。类别相同的离散点属于同一个加权维诺图胞腔，将所有离散点标上类别即求得近似离散加权维诺图。

作为优选，步骤4中，全离散最优传输映射的凸能量函数计算公式为：

凸能量函数的梯度计算公式如下：

其中，h表示截距向量，h_i表示目标离散点y_j的权重；M表示目标离散点个数；q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和，m_k表示离散加权维诺图胞腔Pow(y_j)所包含的源区域离散点的个数，μ_i为源离散点x_i上的狄拉克测度；ν_j为目标离散点y_j上赋予的狄拉克测度。

作为优选，步骤5中，截距向量h的更新规则如下：

其中，stepLength为下降步长，步长越大，能量函数收敛越快，但结果越不精确；步长越小，结果越精确，但收敛速度越慢；表示凸能量函数的梯度。

作为优选，步骤6中，每次迭代，重新计算矩阵S与矩阵T'的乘积，矩阵P元素的值被更新，因此源区域的离散点所属的胞腔被重新标记，得到新的离散加权维诺图和凸能量函数；当凸能量函数达到极小值时，算法终止，得到最终的全离散最优传输算法的结果。

作为优选，所述凸能量函数达到极小值的条件是，梯度的二范数小于预设的收敛误差，即下述条件成立：

其中，M表示目标离散点个数，表示离散加权维诺图胞腔Pow(y_j)的能量梯度；q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和，ν_j为目标离散点y_j上赋予的狄拉克测度，m_k表示离散加权维诺图胞腔Pow(y_j)所包含的源区域离散点的个数，μ_k表示源离散点x_k的狄拉克测度。

本发明的有益效果在于：本发明提供了一种高效的全离散最优传输新算法，通过将源区域和目标区域离散化，表示成一簇离散点，将复杂的加权维诺图构造、能量函数及梯度的计算转换成简单的矩阵相乘和求最大值问题。方法如下：将源离散点的坐标及截距向量的系数存储于一个矩阵中，将目标离散点及截距向量存于另一个矩阵中。通过矩阵相乘和求上包络的方法，将源离散点分类，以构造近似离散加权维诺图。再通过迭代法更新截距向量和离散加权维诺图，优化凸能量函数，得到最终的最优传输结果。本发明提供的全离散最优传输新算法，简单高效，易于实现，通过多线程或CUDA等并行计算技术进一步加速提高计算效率，可用于求解大规模最优传输问题。此外，本发明提出的全离散最优传输算法与源区域的维度无关，适用于求解任意维度的最优传输问题，这是现有最优传输算法所不具有的特性。

附图说明

图1为本发明实施例的全离散最优传输问题示意图；

图2为本发明实施例的流程图；

图3为本发明实施例流程图的补充示意图；

图4为本发明实施例中二维最优传输结果示例图；

图5为本发明实施例中三维最优传输结果示例图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明的全离散最优传输问题是：源区域被离散化，表示成一簇离散点，每个离散点x_i上带有一个狄拉克测度μ_i。目标区域也被离散化，表示成一族离散点，每个目标离散点y_j上赋予一个狄拉克测度ν_j。求解将源区域的离散点映射到目标离散点的最优传输映射f，使得映射f满足以下条件：映射前后测度相等，即且传输代价最小，即映射的凸能量取得极小值。

为了解决上述问题，请见图2和图3，本发明提供了一种高效的全离散最优传输方法，包括以下步骤：

在本实施方式中，输入为一簇任意维的源离散点X＝{x₁,x₂,···,x_N}及其上的狄拉克测度μ_i，一簇目标离散点Y＝{y₁,y₂,···,y_M}及其上赋予的狄拉克测度。如图3，图(a)表示连续的源区域，图(b)中绿色点为源区域的离散化表示，红色点为目标离散点。为了方便计算，源狄拉克测度和目标狄拉克测度满足以下条件：

其中，N表示源区域的离散点总数，μ_i表示源点x_i上的狄拉克测度；M表示目标离散点个数，ν_j为目标离散点y_j上赋予的狄拉克测度。

二维图像可离散表示成像素，三维曲面可用三角网格近似表示，进而将其离散化成一簇顶点，三维体网格可离散化成体素表示。

随机生成目标离散点和赋予测度，或者根据实际需要生成目标离散点和相应的狄拉克测度。

步骤2：初始化截距向量，将源区域离散点的坐标和截距向量的系数用矩阵存储，将目标离散点的坐标及截距向量用矩阵存储；；

将源离散点的信息用大小为N×(D+1)的矩阵S存储，其中D表示源离散点的维度，矩阵的前D列存储源离散点的坐标信息，第D+1列存储截距向量系数，值为1。

其中，(a_i,b_i,···)表示离散点x_i的顶点信息，第D+1列表示截距向量的系数全为1。

在二维情况下，

其中，(a_i,b_i)分别表示离散点x_i的X轴坐标和Y轴坐标。

同理，三维情况下，

其中，(a_i,b_i,c_i)分别表示离散点x_i的X轴坐标、Y轴坐标和Z轴坐标。

截距向量h可以初始化为零向量，即h＝0，或者初始化为随机(0,1)之间的随机值，也可以初始化为每个目标离散点y_j到原点距离的一半，即其中<,>表示点积。将目标离散点的信息和截距向量存于一个大小为M×(D+1)的矩阵T，其中前D维为目标离散点的坐标信息，第D+1为权重。

即，在二维情况下，其中(c_j,d_j)分别表示目标离散点y_j的X轴坐标和Y轴坐标。

在三维情况下，其中(c_j,d_j,e_j)分别表示目标离散点y_j的X轴坐标，Y轴坐标和Z轴坐标。

步骤3：通过矩阵相乘、根据行最大值将离散点分类的方法构造近似离散加权维诺图；

每次迭代中，矩阵S与矩阵T'相乘得到大小为N×M的矩阵P。矩阵P的元素p_ij表示源离散点x_i与目标离散点y_j的点乘与截距向量之和，即p_ij＝<x_i,y_j>+h_j。

根据加权维诺图的定义：

可推导出近似的离散加权维诺图为：

因此可以根据矩阵P每一行元素最大值所在的列将源离散点x_i进行分类，例如，若第i行的最大值在第j列，则将源离散点x_i的类别标记为j。每次迭代，x_i的分类都根据最大值进行更新。

将所有的离散点都标记类别之后，即获得了这次迭代的离散加权维诺图。图3(c)显示某次迭代过程的结果，此时离散加权维诺图胞腔Pow(y_j)的测度与目标离散点y_j上赋予的狄拉克测度ν_j并不相等。

步骤4：求全离散最优传输的凸能量函数及其梯度；

全离散最优传输映射的运输代价为目标离散点的加权维诺图中包含的源离散点的总测度与目标离散点给定的测度之差的总和，即：

其中，h表示截距向量，h_i表示目标离散点y_j的权重；M表示目标离散点个数；q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和，m_k表示类别标记为j的离散点的个数，μ_i为源离散点x_i上的狄拉克测度；ν_j为目标离散点y_j上赋予的狄拉克测度。

则能量函数的梯度函数为：

步骤5：更新截距向量；

由加权维诺图的定义可知，截距向量h决定离散加权维诺图，进而决定全离散最优传输的运输代价和最终结果，截距向量h的更新规则如下：

其中，表示凸能量函数的梯度。stepLength为下降步长，步长越大，能量函数收敛越快，但结果越不精确；步长越小，结果越精确，但收敛速度越慢。在本实施例中，为了获得较精确的结果和较快的收敛速度，将迭代步长stepLength的阈值设为0.001。

步骤6：循环步骤3-步骤5直到传输映射的凸能量函数达到极小值，得到最终的全离散最优传输的结果；

每次迭代，重新计算矩阵S与矩阵T'的乘积，矩阵P的元素被更新，因此源区域的离散点所属的胞腔被重新更新，得到新的离散加权维诺图和凸能量函数。当凸能量函数取得极小值时，算法终止，得到最终的全离散最优传输算法的结果。此时加权维诺图胞腔Pow(y_j)的测度与目标离散点y_j上赋予的狄拉克测度ν_j相等。如图3(d)所示，当所有目标离散点上的狄拉克测度相等时，加权离散维诺图每个胞腔的面积也相等。

凸能量函数达到极小值的条件是，梯度的二范数小于预设的收敛误差，即

其中，M表示目标离散点个数，q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和；ν_j为目标离散点y_j上赋予的狄拉克测度，m_k表示离散加权维诺图胞腔Pow(y_j)所包含的源区域离散点的个数，μ_k表示源离散点x_k的狄拉克测度。在本实施例中，收敛阈值设为10^-6。

在本实施例中，通过将源区域和目标区域离散化表示，将加权维诺图的求解转换成矩阵相乘求每行最大值问题，简化了全离散最优传输代价的计算以及能量函数的优化，大大提高了计算效率。此外，本发明提出的全离散最优传输算法适用于任意维度的最优传输问题。图4是二维的最优传输结果示例，其中源区域由1000000个均匀采样的离散点组成，其上的狄拉克测度定义为0.000001。目标离散点是随机生成的100个点，其上赋予的狄拉克测度均为0.01，截距向量初始化为(0,1)之间的随机数。由图4可知，每个胞腔的离散点个数相等，即胞腔的面积相等。图5是三维的最优传输结果，源区域由N＝1000*1000*1000个均匀采样点组成，每个离散点上的狄拉克测度为1/N；随机生成100个目标离散点，其上赋予的狄拉克测度为0.01。其中图(a)为正视图，图(b)为旋转一定角度后的结果图，图(c)为截面图。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。以上优选实施例虽然只是二维和三维的，但是本发明提出的高效的全离散最优传输新算法适用于求解任意维度的最优传输问题。此外，以上优选实施例中源区域的离散点均匀采样，测度均匀分布，目标离散点的测度均匀分布，但本算法同样适用于源狄拉克测度、目标测度不均匀分布等情况。

Claims

1.一种高效的全离散最优传输方法，其特征在于，包括以下步骤：

步骤4：求全离散最优传输的凸能量函数及其梯度；

步骤5：更新截距向量；

2.根据权利要求1所述的高效的全离散最优传输方法，其特征在于：步骤1中，源区域和目标区域都离散化表示，而且它们的总测度相等，即：

<mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>&mu;</mi> <mi>i</mi> </msub> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mi>v</mi> <mi>j</mi> </msub> </mrow>

其中，N表示源区域的离散顶点总数，μ_i表示源离散点x_i上的狄拉克测度；M表示目标离散点个数，ν_j为目标离散点y_j上赋予的狄拉克测度。

3.根据权利要求2所述的高效的全离散最优传输方法，其特征在于：步骤2中，将源离散点的坐标信息及截距向量的系数用大小为N×(D+1)的矩阵S存储，其中N表示源区域的离散顶点总数，D表示源离散点的维度，矩阵的前D列存储源离散点的坐标信息，第D+1列存储截距向量的系数，值为1；

将目标离散点的坐标信息和截距向量存于一个大小为M×(D+1)的矩阵T中，其中前D列为目标离散点的坐标信息，第D+1列为截距向量h；

其中，(a_i,b_i,···)表示源区域离散点x_i的顶点坐标，第D+1列表示截距向量的系数为1；(c_j,d_j,···)表示目标区域离散点y_j的顶点坐标，h_j为该目标离散点的权重。

4.根据权利要求3所述的高效的全离散最优传输方法，其特征在于：步骤3中，根据加权维诺图的定义：

<mrow> <mi>P</mi> <mi>o</mi> <mi>w</mi> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>:</mo> <mo>=</mo> <mo>{</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>|</mo> <mo>|</mo> <mo>|</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>h</mi> <mi>j</mi> </msub> <mo><</mo> <mo>|</mo> <mo>|</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>-</mo> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>,</mo> <mo>&ForAll;</mo> <mi>j</mi> <mo>&NotEqual;</mo> <mi>k</mi> <mo>}</mo> <mo>,</mo> </mrow>

可推导出近似的离散加权维诺图：

<mrow> <mi>P</mi> <mi>o</mi> <mi>w</mi> <mrow> <mo>(</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>:</mo> <mo>=</mo> <mo>{</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>|</mo> <mo><</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>></mo> <mo>+</mo> <msub> <mi>h</mi> <mi>j</mi> </msub> <mo>></mo> <mo><</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>k</mi> </msub> <mo>></mo> <mo>+</mo> <msub> <mi>h</mi> <mi>k</mi> </msub> <mo>,</mo> <mo>&ForAll;</mo> <mi>j</mi> <mo>&NotEqual;</mo> <mi>k</mi> <mo>}</mo> </mrow>

因此，求解近似离散加权维诺图的具体方法如下：

每次迭代过程中，将矩阵S与T的转置矩阵T'相乘得到大小为N×M的矩阵P，即：

<mrow> <mi>P</mi> <mo>=</mo> <mi>S</mi> <mo>&times;</mo> <msup> <mi>T</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <msub> <mi>p</mi> <mn>11</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>12</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>13</mn> </msub> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mn>1</mn> <mi>M</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>p</mi> <mn>21</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>22</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>23</mn> </msub> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mn>2</mn> <mi>M</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <msub> <mi>p</mi> <mn>31</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>32</mn> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mn>33</mn> </msub> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mn>3</mn> <mi>M</mi> </mrow> </msub> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mo>.</mo> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <mo>.</mo> </mtd> </mtr> <mtr> <mtd> <msub> <mi>p</mi> <mrow> <mi>N</mi> <mn>1</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mi>N</mi> <mn>2</mn> </mrow> </msub> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mi>N</mi> <mn>3</mn> </mrow> </msub> </mtd> <mtd> <mn>...</mn> </mtd> <mtd> <msub> <mi>p</mi> <mrow> <mi>N</mi> <mi>M</mi> </mrow> </msub> </mtd> </mtr> </mtable> </mfenced> </mrow>

矩阵P表示一系列线性函数，其元素p_ij为p_ij＝<x_i,y_j>+h_j，表示源离散点x_i与目标离散点y_j的点积与截距向量之和；根据推导的近似离散加权维诺图定义：可知，每次迭代过程中，这一系列线性函数的上包络即为所要求的加权维诺图；因此依据矩阵P每一行元素最大值所在的列将源离散点x_i进行分类；每次迭代，x_i所属胞腔都根据线性函数的最大值进行更新。

5.根据权利要求4所述的高效的全离散最优传输方法，其特征在于：步骤4中，全离散最优传输映射的凸能量函数计算公式为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>E</mi> <mrow> <mo>(</mo> <mi>h</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <mrow> <mo>(</mo> <msub> <mi>q</mi> <mi>j</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>,</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>q</mi> <mi>j</mi> </msub> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>m</mi> <mi>k</mi> </msub> </munderover> <msub> <mi>&mu;</mi> <mi>k</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>

凸能量函数的梯度计算公式如下：

其中，h表示截距向量，决定每个目标离散点的权重；M表示目标离散点个数；q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和，ν_j为目标离散点y_j上赋予的狄拉克测度；m_k表示离散加权维诺图胞腔Pow(y_j)所包含的源区域离散点的个数。

6.根据权利要求3所述的高效的全离散最优传输方法，其特征在于：步骤5中，截距向量h的更新规则如下：

<mrow> <mi>h</mi> <mo>=</mo> <mi>h</mi> <mo>+</mo> <mi>s</mi> <mi>t</mi> <mi>e</mi> <mi>p</mi> <mi>L</mi> <mi>e</mi> <mi>n</mi> <mi>g</mi> <mi>t</mi> <mi>h</mi> <mo>&times;</mo> <mo>&dtri;</mo> <mi>E</mi> <mrow> <mo>(</mo> <mi>h</mi> <mo>)</mo> </mrow> </mrow>

7.根据权利要求4所述的高效的全离散最优传输方法，其特征在于：步骤6中，每次迭代，重新计算矩阵S与矩阵T'的乘积，矩阵P的元素被更新，因此源区域的离散点所属的胞腔被重新更新，得到新的离散加权维诺图和凸能量函数；当凸能量函数取得极小值时，算法终止，得到最终的全离散最优传输算法的结果。

8.根据权利要求7所述的高效的全离散最优传输方法，其特征在于：所述凸能量函数达到极值的条件是，梯度的二范数小于预设的收敛误差,即满足以下条件：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msqrt> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msup> <mrow> <mo>(</mo> <mo>&dtri;</mo> <mi>E</mi> <mo>(</mo> <msub> <mi>h</mi> <mi>i</mi> </msub> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>=</mo> <msqrt> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>q</mi> <mi>j</mi> </msub> <mo>-</mo> <msub> <mi>v</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo><</mo> <mi>t</mi> <mi>h</mi> <mi>r</mi> <mi>e</mi> <mi>s</mi> <mi>h</mi> <mi>o</mi> <mi>l</mi> <mi>d</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>q</mi> <mi>j</mi> </msub> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>m</mi> <mi>k</mi> </msub> </munderover> <msub> <mi>&mu;</mi> <mi>k</mi> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced>

其中，M表示目标离散点个数，q_j表示离散加权维诺图胞腔Pow(y_j)内所包含的源狄拉克测度的总和；ν_j为目标离散点y_j上赋予的狄拉克测度，m_k表示离散加权维诺图胞腔Pow(y_j)所包含的源区域离散点的个数，μ_k表示源离散点x_k的狄拉克测度。