CN111080509B - 三维裁剪Voronoi图的多线程并行计算方法、系统 - Google Patents

三维裁剪Voronoi图的多线程并行计算方法、系统 Download PDF

Info

Publication number
CN111080509B
CN111080509B CN201911255557.9A CN201911255557A CN111080509B CN 111080509 B CN111080509 B CN 111080509B CN 201911255557 A CN201911255557 A CN 201911255557A CN 111080509 B CN111080509 B CN 111080509B
Authority
CN
China
Prior art keywords
voronoi
tetrahedron
sites
neighbor
clipping
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911255557.9A
Other languages
English (en)
Other versions
CN111080509A (zh
Inventor
严冬明
刘晓寒
郭建伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Automation of Chinese Academy of Science
Original Assignee
Institute of Automation of Chinese Academy of Science
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Automation of Chinese Academy of Science filed Critical Institute of Automation of Chinese Academy of Science
Priority to CN201911255557.9A priority Critical patent/CN111080509B/zh
Publication of CN111080509A publication Critical patent/CN111080509A/zh
Application granted granted Critical
Publication of CN111080509B publication Critical patent/CN111080509B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T1/00—General purpose image data processing
    • G06T1/20—Processor architectures; Processor configuration, e.g. pipelining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F9/00—Arrangements for program control, e.g. control units
    • G06F9/06—Arrangements for program control, e.g. control units using stored programs, i.e. using an internal store of processing equipment to receive or retain programs
    • G06F9/30—Arrangements for executing machine instructions, e.g. instruction decode
    • G06F9/38—Concurrent instruction execution, e.g. pipeline or look ahead
    • G06F9/3885—Concurrent instruction execution, e.g. pipeline or look ahead using a plurality of independent parallel functional units

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明属于数字图像处理领域，具体涉及一种三维裁剪Voronoi图的多线程并行计算方法、系统，旨在为了解决现有三维裁剪Voronoi图计算技术在效率上的不足的问题。本发明方法包括获取每个采样站点的k近邻站点；基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；利用相交的部分更新对应的Voronoi单元Ωj。本发明本发明的极大提高了三维裁剪Voronoi图计算处理效率。

Description

三维裁剪Voronoi图的多线程并行计算方法、系统
技术领域
本发明属于数字图像处理领域，具体涉及一种三维裁剪Voronoi图的多线程并行计算方法、系统。
背景技术
Voronoi图(Voronoi Diagram，也被称为Dirichlet Tessellation，即迪利克雷镶嵌)是俄国数学家格奥尔吉·沃罗诺伊(Georgii Voronoi)于1908年提出的一种空间分割算法。该算法依据一组有限点集对空间进行分解，这些点被称为站点(site)。在分解中，每个站点对应一个单元(cell)，该单元由在空间中相比其他站点距离该站点更近的所有点组成。Voronoi图在自然界和日常生活中十分常见，如今已成为计算几何中的一项基础工具，被国内外学者广泛地应用于各领域的实际问题中，在网格优化、三维重建、地理信息系统、气象预测、城市规划等领域均取得了出色的应用成果。由于在普通的Voronoi图中，位于其对偶的Delaunay三角剖分凸包上的站点对应无限大的单元，因此在实际使用时，往往需要将Voronoi图限制在一个有限域内，即使用域的边界对普通Voronoi图进行裁剪。
另一方面，随着地理信息系统的迅速发展以及建设数字城市的需要，对大规模三维信息进行处理的需求变得日益迫切。以目前的计算技术来看，要想在消费级设备上处理如此大规模的数据，离不开基于图形处理器(Graphics Processing Unit,GPU)的并行计算。近年来，在人工智能与游戏市场双重需求的推动下，GPU的计算能力已经有了飞速的发展，具备了实现大规模并行计算的能力。其中，NVIDIA公司提出的统一计算架构(ComputeUnified Device Architecture,CUDA)就是该技术的代表。CUDA凭借其高性价比、低通信开销、卓越的并行计算能力，使得针对海量化且计算复杂度高的三维模型数据的快速处理成为可能。CUDA采用C语言作为编程语言，提供了大量高性能计算指令的开发能力，使开发者能够在GPU的强大计算能力的基础上建立起一种效率更高的密集数据计算解决方案。利用CUDA编程，我们可以一次并行地开启大量线程进行运算，从而大幅度提升运算效率且不会显著降低结果质量。
目前，尽管国内外的研究学者提出了多种高效计算普通Voronoi图的方法，并且存在一些GPU上的实现，但裁剪Voronoi图沉重的计算成本仍在阻碍其在许多应用中的广泛使用。Rong等人(RONG G,LIU Y,WANG W,et al.Gpu-assisted computation of centroidalvoronoi tessellation[J].IEEE Transactions on Visualization and ComputerGraphics,2011,17:345-356.)通过几何图像参数化计算曲面Voronoi图，并将跳跃泛洪算法(RONG G,TAN T S.Jump flooding in gpu with applications to voronoi diagramand distance transform[C]//SI3D.2006.)用于二维Voronoi图的计算。Han等人(HAN J,YAN D M,WANG L,et al.Computing restricted voronoi diagram on graphicshardware[C]//Pacific Graphics 2017.Eurographics Association,2017:23-26.)通过在GPU上应用多边形裁剪算法来计算曲面网格上的受限Voronoi图。最近，Ray等人(RAY N,SOKOLOV D,LEFEBVRE S,et al.Meshless voronoi on the gpu[J].ACM Trans.onGraphics,2018,37(6):265:1-265:12.)在GPU上将Voronoi单元以对偶的形式存储为简单的三角形网格，从而实现了大规模Voronoi图的并行计算。但是，该方法只能计算无边界的普通Voronoi图，这也是目前裁剪Voronoi图的瓶颈所在。
发明内容
为了解决现有技术中的上述问题，即为了解决现有三维裁剪Voronoi图计算技术在效率上的不足的问题，本发明的一方面，提出了一种三维裁剪Voronoi图的多线程并行计算方法，该方法包括以下步骤：
步骤S100，获取每个采样站点的k近邻站点；
步骤S200，基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
步骤S300，将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
步骤S400，利用相交的部分更新对应的Voronoi单元Ωj。
在一些优选实时方式中，所述计算方法中，
步骤S100还包括获取三维网格中每个四面体的k近邻站点；
步骤S300中“将每个Voronoi单元Ωj分别与三维网格中的四面体配对”其方法为：将k近邻站点包含采样站点xj的四面体ti与Voronoi单元Ωj进行配对，获得多个“四面体-单元”对。
在一些优选实时方式中，所述并行计算，其方法为：将独立计算的各部分内容分配至多线程同时计算。
在一些优选实时方式中，步骤S100中“获取三维网格中每个四面体的k近邻站点”，其方法为：
对于四面体ti，计算其质心ci并利用k近邻查询得到距其质心最近的k个站点。
在一些优选实时方式中，步骤S100中每个四面体的k近邻站点的数量k值，其计算方法为
Figure BDA0002310149880000031
其中，λ、α、β、γ为预设参数，
Figure BDA0002310149880000042
为三维网格中四面体的数量，|X|为站点的数量。
在一些优选实时方式中，λ＝20,α＝10,β＝8,γ＝7。
在一些优选实时方式中，步骤S300中“获取对应的相交部分”，其方法为：基于每个Voronoi单元Ωj、及其配对的四面体，
利用四面体的四个面对Voronoi单元Ωj进行裁剪，获取对应的相交部分；和/或
利用Voronoi单元Ωj的面，对四面体进行裁剪，获取对应的相交部分。
在一些优选实时方式中，步骤S400中“利用相交的部分更新对应的Voronoi单元Ωj”，其方法为：
将相交的部分的凸多面体P分解为由若干四面体组成的集合，然后计算每个四面体的质心及体积，以体积作为权重计算质心坐标的加权和，同时对体积进行累加；
利用atomic_Add原子操作，将质心坐标加权和与体积分别其累加到裁剪Voronoi单元
Figure BDA0002310149880000041
的加权和与总体积中。
在一些优选实时方式中，k近邻站点的计算方法为：
基于预设的站点数阈值，若k近邻站点数大于该阈值，采用k近邻搜索算法获取k近邻站点，否则采用蛮力法获取k近邻站点。
本发明的第二方面，提出了一种三维裁剪Voronoi图的多线程并行计算系统，该系统包括：第一模块、第二模块、第三模块、第四模块；
所述第一模块，配置为获取每个采样站点的k近邻站点；
所述第二模块，配置为基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
所述第三模块，配置为将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
所述第四模块，配置为利用相交的部分更新对应的Voronoi单元Ωj。
本发明的第三方面，提出了一种存储装置，其中存储有多条程序，其特征在于，所述程序适于由处理器加载并执行以实现上述的三维裁剪Voronoi图的多线程并行计算方法。
本发明的第四方面，提出了一种处理装置，包括处理器、存储装置；处理器，适于执行各条程序；存储装置，适于存储多条程序；其特征在于，所述程序适于由处理器加载并执行以实现上述的三维裁剪Voronoi图的多线程并行计算方法。
本发明的有益效果：
本发明采用多线程并行计算技术，使得裁剪Voronoi图的每个相交区域都可以在一个单独的线程上进行计算，本发明的极大提高了三维裁剪Voronoi图计算处理效率，在网格质量优化、三维重建等领域有很好的应用价值。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：
图1是本发明一种实施例的三维裁剪Voronoi图的多线程并行计算方法流程示意图；
图2是本发明一种实施例中“单元-四面体”策略的二维示意图；
图3是本发明一种实施例中“四面体-单元”策略的二维示意图；
图4是本发明一种实施例中构造普通Voronoi单元的示意图；
图5是本发明一种实施例中使用四面体对Voronoi单元进行裁剪的示意图；
图6是本发明一种实施例与Geogram的方法在Joint模型(四面体数为31K)上，不断增加站点数得到的计算时间曲线图；
图7是本发明一种实施例与Geogram的方法在Joint模型上，固定站点数为15K，不断增加网格模型四面体数得到的计算时间曲线图；
图8是本发明一种实施例中进行Lloyd迭代算法计算质心Voronoi图的结果示例图；
图9是适于用来实现本申请实施例的电子设备的计算机系统的结构示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。
本发明的一种三维裁剪Voronoi图的多线程并行计算方法，如图1所示，该方法包括以下步骤：
步骤S100，获取每个采样站点的k近邻站点；
步骤S200，基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
步骤S300，将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
步骤S400，利用相交的部分更新对应的Voronoi单元Ωj。
为了更清晰地对本发明三维裁剪Voronoi图的多线程并行计算方法进行说明，下面结合附图对本方发明方法一种实施例中各步骤进行展开详述。
本发明一种实施例的三维裁剪Voronoi图的多线程并行计算方法，包括步骤S100-步骤S400。本实施例中，通过多个GPU进行多线程并行计算。
步骤S100，k近邻站点的获取。
获取每个采样站点的k近邻站点：利用k近邻查询得到距其最近的ksite个站点；获取三维网格中每个四面体的k近邻站点：对于四面体ti，计算其质心ci并利用k近邻查询得到距其质心最近的k个站点。同时，基于预设的站点数阈值，若k近邻站点数大于该阈值，采用k近邻搜索算法获取k近邻站点，否则采用蛮力法获取k近邻站点。本实施例中站点数阈值为14，当然还可以为其他阈值。
本发明实施例是一种用于计算限制于三维域空间
Figure BDA0002310149880000071
中的裁剪Voronoi图的高效算法。为便于计算，输入的三维域空间被离散化为四面体网格。该算法以一组站点X和一个四面体网格
Figure BDA0002310149880000072
作为输入，其中
Figure BDA0002310149880000073
是网格顶点的集合，
Figure BDA0002310149880000074
是四面体元素的集合。在四面体网格中，每个四面体ti由四个有序顶点表示，索引值分别为0、1、2和3。任意三个顶点的组合都可以组成一个三角形面，该面的索引值与其对侧的顶点索引值相同。
k近邻搜索是一种非常实用的空间搜索技术，在众多领域都有着广泛的应用，如三维几何处理、模式识别等。本发明提出的方法也是基于k近邻搜索技术。目前，这一技术已有多种可用的CPU和GPU实现，且计算效率能够满足实际应用的需求。由于本发明主要是在GPU上进行计算，为了避免在CPU和GPU之间进行频繁的数据传输，本发明选择了两种基于GPU的k近邻搜索算法，分别适用于不同场景下的查询操作。考虑到需要处理的数据仅为三维数据，本发明采用了一种基于网格划分的k近邻搜索算法用于计算每个站点的k近邻。然而，该方法在站点数量较少时(约为14K)无法正常进行计算，并且要求查询点只能为输入点集内的点。因此，本发明同时采用了一种简单的蛮力法来计算每个四面体的k近邻站点，以及当站点数量少于14K时每个站点的k近邻。得益于GPU高度并行的架构，这种简单的蛮力法在本发明的应用场景下也能取得与网格划分法相当的性能。下面对这两种方法分别进行简要介绍。
蛮力法的算法思想非常简单，直接计算每两点之间的距离并选出距离最小的k个即可。给定d维空间中的一组参考的
Figure BDA0002310149880000081
以及在同一空间中一组查询点
Figure BDA0002310149880000082
距离查询点最近的k个参考点可以通过如下步骤计算：
计算当前查询点与每个参考点直接的距离di＝d(q,ri),1≤i≤m；
对得到的距离
Figure BDA0002310149880000083
按照升序排序；
输出距离值最小的k个参考点的索引值。
上述每个步骤都是可以高度并行化的，适合在GPU上进行实现。
基于网格划分的方法以一组三维空间点集
Figure BDA0002310149880000091
为输入，对每个点集中的点x∈X，输出距离其最近的k个点的索引值。首先，通过四舍五入将输入点嵌入到由若干单元组成的三维立方体网格中，平均每个单元约包含5个点。然后，按照单元格索引的升序对所有点进行排序，以便快速检索给定单元格内的所有点。
网格结构划分完成后，对于每个查询点，都能够以同心圆的形状快速找到其附近的网格单元，从而访问包含在这些单元内的邻居。在遍历访问其邻居时，需要维护一个最大二叉堆来存储k个候选点。如果查询点与其邻居之间的距离小于堆中的最大距离，则将最大元素替换为邻居并调整堆的结构。在访问以查询点所在单元为中心的周围若干层单元的过程中，如果堆中有k个元素，并且最大距离小于到下一层的最小距离，则结束查询并返回结果。
对于本发明提出的算法而言，至关重要的一点是如何确定一个四面体需要裁剪的近邻Voronoi单元数量，即参数k的值。一般来说，采样站点可以看作是大致均匀的分布，因此可以根据站点数量与四面体数量的比值来确定k的值。为了提高算法的实用性，本发明给出了一个公式来计算k的参考值：
Figure BDA0002310149880000092
其中，其中，λ、α、β、γ为预设参数，
Figure BDA0002310149880000093
为三维网格中四面体的数量，|X|为站点的数量。本实施例中，λ＝20,α＝10,β＝8,γ＝7。当然，用户也可以手动地对参数k进行调整从而在效率与准确度之间进行权衡。
步骤S200，基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj。
根据Voronoi图的定义，每个站点对应的Voronoi单元都是三维空间
Figure BDA0002310149880000101
的一个子集，包含了相比于其他站点到该站点距离更近的所有点。显然，站点xi对应的Voronoi单元Ωi可以由一系列半空间求交得到，即Ωi＝∩j≠iΠ+(i,j)，其中Π+(i,j)表示由两站点(xi,xj)垂直平分面定义的且包含xi一侧的半空间。然而，由于距离较远的站点实际上不会对Voronoi单元有任何贡献，因此无需用其余所有站点对应的半空间对Voronoi单元进行裁剪。基于这一性质，对于每次计算得到的Voronoi单元，可以考虑一个以站点xi为球心的单元包围球，记其半径为R。如果当前站点与下一站点(按照距离升序排序)的距离d(xi,xj)>2R，那么下一站点对应的裁剪平面将无法接触到当前的包围球，也就不会对Voronoi单元产生影响。此时的半径R被称为“安全半径”。如图4所示，令x1,x2,…,x10表示站点xi的10个邻近站点，按照距离升序排序。可以看出，对Voronoi单元有贡献的平分线会裁剪掉一个半空间，并成为Voronoi图的一条边，而与站点xi的距离大于2R的站点x7,x8,…,x10则不会对Voronoi单元做出贡献。因此，我们可以由近到远地使用站点xi与其他站点定义的半空间对整个空间依次进行裁剪，并在达到安全半径后停止裁剪，从而得到站点xi对应的普通Voronoi单元Ωi。同时，由于每个Voronoi单元的构造只需要其若干邻居，且构造过程是独立的，不会影响其他单元，故每个单元都可以被并行计算。
步骤S300，将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分。
本实施例中，“将每个Voronoi单元Ωj分别与三维网格中的四面体配对”其方法为：将k近邻站点包含采样站点xj的四面体ti与Voronoi单元Ωj进行配对，获得多个“四面体-单元”对。
本实施例中，“获取对应的相交部分”，其方法为：基于每个Voronoi单元Ωj、及其配对的四面体，利用四面体的四个面对Voronoi单元Ωj进行裁剪，获取对应的相交部分；和/或利用Voronoi单元Ωj的面，对四面体进行裁剪，获取对应的相交部分。
计算相交区域两种策略可以进一步描述为：第一种是先计算普通Voronoi单元，然后用由四面体的四个面定义的半空间分别去裁剪相应单元，从而获得Voronoi单元位于四面体内部的部分，即“单元-四面体”策略(如图2所示，图中xl为采样站点，x1-x5为其对应的k近邻站点)。另一种策略与现有方法类似，以四面体作为初始结构，用由一个站点及其最近的k个站点的等分线所界定的半空间来裁剪该四面体，即“四面体-单元策略”(如图3所示，图中xl为采样站点，x1-x5为其对应的k近邻站点)。接下来详细介绍并比较以上两种策略。
对于“单元-四面体”策略，首先需要构造普通Voronoi图。得到普通Voronoi单元后，关键问题就在于如何将四面体和Voronoi单元进行配对，并计算出它们的相交区域。对于每个四面体ti，首先计算其质心ci并利用k近邻查询得到距其质心最近的k个站点。然后，对于每个邻居站点xj，四面体ti和站点xj对应的Voronoi单元Ωj可以构成一个“四面体-单元”配对(ti,Ωj)。对每个配对，都可以创建一个线程来计算Voronoi单元在四面体内部的部分。线程首先将预先计算的Voronoi单元从全局内存加载到其共享内存，然后使用由四面体的四个面定义的四个半空间对单元进行裁剪。令{v1,v2,v3}为四面体某个面的三个顶点，法线方向n＝(nx,ny,nz)指向四面体内侧，其序号按照右手定则排序。则半空间的方程可以根据一个顶点(如v1)和其法向n得到：
n＝(v2-v1)×(v3-v1)
Figure BDA0002310149880000111
图5展示了整个裁剪过程，普通Voronoi单元(凸多面体)依次被由四面体的四个面所定义的半空间裁剪。如果四面体和Voronoi单元之间的相交区域不为空，则将使用步骤S400的方法并借助atomicAdd原子操作更新裁剪Voronoi单元
Figure BDA0002310149880000121
的总体积和质心坐标。值得注意的是，以上所有计算操作都是在GPU上进行实现的，避免了在CPU和GPU之间较为耗时的数据传输。
在“单元-四面体”策略中，普通Voronoi单元需要先存储在GPU的全局内存中，然后再将其加载到共享内存中。众所周知，共享内存的存取速度比全局内存要快很多。因此，本发明同时提出了“四面体-单元”策略来直接计算裁剪Voronoi图，避免了此处对全局内存的存取。图3展示了“四面体-单元”策略的裁剪流程。
“四面体-单元”策略使用了与“单元-四面体”策略相同的方法来对四面体和其近邻Voronoi单元进行配对，不同之处在于此处使用四面体作为初始区域，并且直接使用由当前站点与其近邻站点的垂直平分面界定的半空间对四面体初始区域进行裁剪。
步骤S400，利用相交的部分更新对应的Voronoi单元Ωj。
在该步骤中，基于GPU上的原子操作，制定了一种维护Voronoi单元的质心和体积的算法。该算法避免了在计算完裁剪Voronoi单元的每个局部区域后，更新对应Voronoi单元时发生的内存写入冲突，保证了计算的正确性，且不会对计算效率产生较大影响，充分利用了现代GPU的并行计算能力。具体算法包括：
若步骤S300计算得到的结果不为空，则裁剪后所剩区域应为某个普通Voronoi单元的一部分，且是一个凸多面体，记为P，同时记其对应的Voronoi单元的索引为j。要将区域P更新到索引为j的裁剪Voronoi单元中，需要计算该区域的体积和质心。首先将区域P分解为由若干四面体组成的集合，然后计算每个四面体的质心及体积，以体积作为权重计算质心坐标的加权和，同时对体积进行累加。最后利用atomicAdd原子操作，将质心坐标加权和与体积分别其累加到裁剪Voronoi单元
Figure BDA0002310149880000131
的加权和与总体积中。(此处以
Figure BDA0002310149880000132
为裁剪Voronoi单元，表示限制于有限域
Figure BDA0002310149880000133
中的Voronoi单元，是本发明的计算目标，以区别于普通Voronoi单元Ωj。Ωj没有有限域的限制，可能存在无限大的单元。)
图6展示了在Joint模型(四面体数为31K)上，不断增加站点数量，对本发明算法与Geogram 1.7.1(project ALICE-GRAPHYS Inria.Geogram:a programming library ofgeometric algorithms[EB/OL].2019.http://alice.loria.fr/software/geogram/doc/html/index.html.)的比较结果。随着采样站点的数量从1K增加到31K，本发明提出的算法在计算效率上的优势越来越明显，并且两种策略均可以达到相同的性能水平。在图7中，站点数量被固定为15K，同时输入网格的四面体数量从10K开始逐渐增加到100K。在这一过程中，可以发现本发明算法的性能几乎没有变化。这是因为当站点数量与四面体数量的比值减小时，参数k的取值也随之减小，从而节约了计算成本，且在大部分情况下不会影响计算结果的准确性。
质心Voronoi图(CVT)是一种特殊的Voronoi图，它的每个站点都与其对应的Voronoi单元的质心重合。CVT在计算机图形学中有许多有用的应用，例如网格划分、采样、点彩画等等。在CVT的计算方法中，目前最流行的方法之一是Lloyd算法，它是一种迭代算法。Lloyd算法首先计算一组站点的裁剪Voronoi图，然后将每个站点移动至其Voronoi单元的质心位置，再重复这一过程。在Lloyd算法中，最困难且最耗时的步骤就是裁剪Voronoi图的计算。因此，本发明提出的基于GPU的算法可以显着提升Lloyd算法的计算性能，同时体积误差也可以被保持在可接受的范围内。如图8所示，我们从输入四面体网格内的3K个随机站点开始执行Lloyd算法。可以发现，经过若干次迭代后，这些站点的分布逐渐变得均匀。我们使用本发明提出的算法分别在Torus(四面体数为8K)、Bunny(四面体数为10K)、Joint(四面体数为31K)和Gargoyle(四面体数为199K)模型上进行了Lloyd迭代，对应图中从上到下的顺序。图中从左到右依次为：输入的四面体网格，3K个初始采样站点的裁剪Voronoi图，以及5、20和120次迭代后的裁剪Voronoi图。在四个模型上完成120次迭代分别需要的时间为10.51秒、9.50秒、10.43秒和32.65秒，与CPU方法相比大约快了一个数量级。
本发明的方法的特色和创新在于，提出了一种简单而有效的三维裁剪Voronoi图的GPU并行计算方法。该算法设计了两种并行策略，使得裁剪Voronoi图的每个相交区域都可以在一个单独的线程上进行计算。与CPU上的方法相比，本发明提出的算法大大提高了三维裁剪Voronoi图的计算效率，为后续工作及应用带来了更多便利。
上述实验结果和计算方法，可以用于网格质量优化和三维重建等多个领域，具有较高的实际应用价值。
本发明第二实施例的一种三维裁剪Voronoi图的多线程并行计算系统，包括：第一模块、第二模块、第三模块、第四模块；
所述第一模块，配置为获取每个采样站点的k近邻站点；
所述第二模块，配置为基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
所述第三模块，配置为将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
所述第四模块，配置为利用相交的部分更新对应的Voronoi单元Ωj。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
需要说明的是，上述实施例提供的三维裁剪Voronoi图的多线程并行计算系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。
本发明第三实施例的一种存储装置，其中存储有多条程序，所述程序适于由处理器加载并执行以实现上述的三维裁剪Voronoi图的多线程并行计算方法。
本发明第四实施例的一种处理装置，包括处理器、存储装置；处理器，适于执行各条程序；存储装置，适于存储多条程序；所述程序适于由处理器加载并执行以实现上述的三维裁剪Voronoi图的多线程并行计算方法。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
下面参考图9，其示出了适于用来实现本申请方法、系统、装置实施例的服务器的计算机系统的结构示意图。图9示出的服务器仅仅是一个示例，不应对本申请实施例的功能和使用范围带来任何限制。
如图9所示，计算机系统包括中央处理单元(CPU，Central Processing Unit)901，其可以根据存储在只读存储器(ROM，Read Only Memory)902中的程序或者从存储部分908加载到随机访问存储器(RAM，Random Access Memory)903中的程序而执行各种适当的动作和处理。在RAM 903中，还存储有系统操作所需的各种程序和数据。CPU 901、ROM 902以及RAM 903通过总线904彼此相连。输入/输出(I/O，Input/Output)接口905也连接至总线904。
以下部件连接至I/O接口905：包括键盘、鼠标等的输入部分906；包括诸如阴极射线管(CRT，Cathode Ray Tube)、液晶显示器(LCD，Liquid Crystal Display)等以及扬声器等的输出部分907；包括硬盘等的存储部分908；以及包括诸如LAN(局域网，Local AreaNetwork)卡、调制解调器等的网络接口卡的通信部分909。通信部分909经由诸如因特网的网络执行通信处理。驱动器910也根据需要连接至I/O接口905。可拆卸介质911，诸如磁盘、光盘、磁光盘、半导体存储器等等，根据需要安装在驱动器910上，以便于从其上读出的计算机程序根据需要被安装入存储部分908。
特别地，根据本公开的实施例，上文参考流程图描述的过程可以被实现为计算机软件程序。例如，本公开的实施例包括一种计算机程序产品，其包括承载在计算机可读介质上的计算机程序，该计算机程序包含用于执行流程图所示的方法的程序代码。在这样的实施例中，该计算机程序可以通过通信部分909从网络上被下载和安装，和/或从可拆卸介质911被安装。在该计算机程序被中央处理单元(CPU)901执行时，执行本申请的方法中限定的上述功能。需要说明的是，本申请上述的计算机可读介质可以是计算机可读信号介质或者计算机可读存储介质或者是上述两者的任意组合。计算机可读存储介质例如可以是——但不限于——电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。计算机可读存储介质的更具体的例子可以包括但不限于：具有一个或多个导线的电连接、便携式计算机磁盘、硬盘、随机访问存储器(RAM)、只读存储器(ROM)、可擦式可编程只读存储器(EPROM或闪存)、光纤、便携式紧凑磁盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。在本申请中，计算机可读存储介质可以是任何包含或存储程序的有形介质，该程序可以被指令执行系统、装置或者器件使用或者与其结合使用。而在本申请中，计算机可读的信号介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了计算机可读的程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。计算机可读的信号介质还可以是计算机可读存储介质以外的任何计算机可读介质，该计算机可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。计算机可读介质上包含的程序代码可以用任何适当的介质传输，包括但不限于：无线、电线、光缆、RF等等，或者上述的任意合适的组合。
可以以一种或多种程序设计语言或其组合来编写用于执行本申请的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
附图中的流程图和框图，图示了按照本申请各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims (10)

1.一种三维裁剪Voronoi图的多线程并行计算方法，其特征在于，该方法包括以下步骤：
步骤S100，获取每个采样站点的k近邻站点；
还包括获取三维网格中每个四面体的k近邻站点；每个采样站点的近邻站点和每个四面体的k近邻站点的数量k值，其计算方法为：
Figure FDA0003486013020000011
其中，λ、α、β、γ为预设参数，
Figure FDA0003486013020000012
为三维网格中四面体的数量，|X|为站点的数量；
步骤S200，基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
由近到远地使用站点xi与其他站点定义的半空间对整个空间依次进行裁剪，并在达到安全半径后停止裁剪，从而得到站点xi对应的普通Voronoi单元Ωi；
步骤S300，将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
具体包括：
将k近邻站点包含采样站点xj的四面体ti与Voronoi单元Ωj进行配对，获得多个“四面体-单元”对；
单元-四面体策略，先计算普通Voronoi单元，然后用由四面体的四个面定义的半空间分别去裁剪相应单元，从而获得Voronoi单元位于四面体内部的部分；
四面体-单元策略，以四面体作为初始结构，用由一个站点及其最近的k个站点的等分线所界定的半空间来裁剪该四面体；
步骤S400，利用相交的部分更新对应的Voronoi单元Ωj。
2.根据权利要求1所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，所述并行计算，其方法为：将独立计算的各部分内容分配至多线程同时计算。
3.根据权利要求1所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，步骤S100中“获取三维网格中每个四面体的k近邻站点”，其方法为：
对于四面体ti，计算其质心ci并利用k近邻查询得到距其质心最近的个站点。
4.根据权利要求1所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，λ＝20，α＝10，β＝8，γ＝7。
5.根据权利要求1-4任一项所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，步骤S300中“获取对应的相交部分”，其方法为：基于每个Voronoi单元Ωj、及其配对的四面体，
利用四面体的四个面对Voronoi单元Ωj进行裁剪，获取对应的相交部分；和/或
利用Voronoi单元Ωj的面，对四面体进行裁剪，获取对应的相交部分。
6.根据权利要求1-4任一项所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，步骤S400中“利用相交的部分更新对应的Voronoi单元Ωj”，其方法为：
将相交的部分的凸多面体P分解为由若干四面体组成的集合，然后计算每个四面体的质心及体积，以体积作为权重计算质心坐标的加权和，同时对体积进行累加；
利用atomic_Add原子操作，将质心坐标加权和与体积分别其累加到裁剪Voronoi单元
Figure FDA0003486013020000033
的加权和与总体积中。
7.根据权利要求1-4任一项所述的三维裁剪Voronoi图的多线程并行计算方法，其特征在于，k近邻站点的计算方法为：
基于预设的站点数阈值，若k近邻站点数大于该阈值，采用k近邻搜索算法获取k近邻站点，否则采用蛮力法获取k近邻站点。
8.一种三维裁剪Voronoi图的多线程并行计算系统，其特征在于，该系统包括：第一模块、第二模块、第三模块、第四模块；
所述第一模块，配置为获取每个采样站点的k近邻站点；
还包括获取三维网格中每个四面体的k近邻站点；
每个采样站点的k近邻站点和每个四面体的k近邻站点的数量k值，其计算方法为：
Figure FDA0003486013020000031
其中，λ、α、β、γ为预设参数，
Figure FDA0003486013020000032
为三维网格中四面体的数量，|X|为站点的数量；
所述第二模块，配置为基于各采样站点及其对应的k近邻站点，通过并行计算分别获取各采样站点xj对应的Voronoi单元Ωj；
由近到远地使用站点xi与其他站点定义的半空间对整个空间依次进行裁剪，并在达到安全半径后停止裁剪，从而得到站点xi对应的普通Voronoi单元Ωi；
所述第三模块，配置为将每个Voronoi单元Ωj分别与三维网格中的四面体配对，通过并行计算分别获取对应的相交部分；
具体包括：
将k近邻站点包含采样站点xj的四面体ti与Voronoi单元Ωj进行配对，获得多个“四面体-单元”对；
单元-四面体策略，先计算普通Voronoi单元，然后用由四面体的四个面定义的半空间分别去裁剪相应单元，从而获得Voronoi单元位于四面体内部的部分；
四面体-单元策略，以四面体作为初始结构，用由一个站点及其最近的k个站点的等分线所界定的半空间来裁剪该四面体；
所述第四模块，配置为利用相交的部分更新对应的Voronoi单元Ωj。
9.一种存储装置，其中存储有多条程序，其特征在于，所述程序适于由处理器加载并执行以实现权利要求1-7任一项所述的三维裁剪Voronoi图的多线程并行计算方法。
10.一种处理装置，包括处理器、存储装置；处理器，适于执行各条程序；存储装置，适于存储多条程序；其特征在于，所述程序适于由处理器加载并执行以实现权利要求1-7任一项所述的三维裁剪Voronoi图的多线程并行计算方法。
CN201911255557.9A 2019-12-10 2019-12-10 三维裁剪Voronoi图的多线程并行计算方法、系统 Active CN111080509B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911255557.9A CN111080509B (zh) 2019-12-10 2019-12-10 三维裁剪Voronoi图的多线程并行计算方法、系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911255557.9A CN111080509B (zh) 2019-12-10 2019-12-10 三维裁剪Voronoi图的多线程并行计算方法、系统

Publications (2)

Publication Number Publication Date
CN111080509A CN111080509A (zh) 2020-04-28
CN111080509B true CN111080509B (zh) 2022-03-08

Family

ID=70313484

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911255557.9A Active CN111080509B (zh) 2019-12-10 2019-12-10 三维裁剪Voronoi图的多线程并行计算方法、系统

Country Status (1)

Country Link
CN (1) CN111080509B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113627019B (zh) * 2021-08-12 2022-03-18 宁波市气象安全技术中心 一种自动气象站覆盖优化方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6384826B1 (en) * 1998-08-14 2002-05-07 Xerox Corporation Method, apparatus and computer medium for surface reconstruction by Voronoi filtering
CN102074052A (zh) * 2011-01-20 2011-05-25 山东理工大学 基于样点拓扑近邻的散乱点云曲面拓扑重建方法
CN103226804A (zh) * 2013-04-12 2013-07-31 山东大学 一种基于流线重心Voronoi图的流场可视化方法
CN104915991A (zh) * 2015-06-04 2015-09-16 南京师范大学 一种针对相交线状地理要素的空间剖分方法
CN108053483A (zh) * 2017-11-03 2018-05-18 北京航空航天大学 一种基于gpu加速的维诺图三维网格重构方法
CN109063332A (zh) * 2018-08-02 2018-12-21 中国科学院自动化研究所 重心Voronoi图的规整性提升方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6384826B1 (en) * 1998-08-14 2002-05-07 Xerox Corporation Method, apparatus and computer medium for surface reconstruction by Voronoi filtering
CN102074052A (zh) * 2011-01-20 2011-05-25 山东理工大学 基于样点拓扑近邻的散乱点云曲面拓扑重建方法
CN103226804A (zh) * 2013-04-12 2013-07-31 山东大学 一种基于流线重心Voronoi图的流场可视化方法
CN104915991A (zh) * 2015-06-04 2015-09-16 南京师范大学 一种针对相交线状地理要素的空间剖分方法
CN108053483A (zh) * 2017-11-03 2018-05-18 北京航空航天大学 一种基于gpu加速的维诺图三维网格重构方法
CN109063332A (zh) * 2018-08-02 2018-12-21 中国科学院自动化研究所 重心Voronoi图的规整性提升方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Computing 3D Clipped Voronoi Diagrams on GPU;Xiaohan Liu等;《In SIGGRAPH Asia 2019 Posters(SA’19 Posters)》;20191117;摘要、正文第1-2节、图2 *
Voronoi图k阶邻近并行矩阵迭代算法;余婧等;《计算机工程与应用》;20141231;第50卷(第6期);第102-105、131页 *
基于Voronoi邻近的移动地图自适应裁剪模型;陈军等;《测绘学报》;20090430;第38卷(第2期);第152-155、161页 *

Also Published As

Publication number Publication date
CN111080509A (zh) 2020-04-28

Similar Documents

Publication Publication Date Title
Huang et al. Explorations of the implementation of a parallel IDW interpolation algorithm in a Linux cluster-based parallel GIS
Bose et al. A survey of geodesic paths on 3D surfaces
US7786991B2 (en) Applications of interval arithmetic for reduction of number of computations in ray tracing problems
Lu et al. Massive point cloud space management method based on octree-like encoding
CN108520142B (zh) 一种城市群边界识别方法、装置、设备及存储介质
US20220011136A1 (en) Road data processing method, apparatus, device, and storage medium
CN108053483A (zh) 一种基于gpu加速的维诺图三维网格重构方法
Biniaz et al. A faster circle-sweep Delaunay triangulation algorithm
WO2024138878A1 (zh) 一种基于空间划分和模型体素化的快速检索方法
Coelho et al. Intersecting and trimming parametric meshes on finite element shells
CN110647596A (zh) 地图数据处理方法和装置
CN108364331A (zh) 一种等值线生成方法、系统和存储介质
CN116030218A (zh) 四面体网格划分方法、装置、系统及存储介质
CN111695281B (zh) 一种四面体网格划分有限元粒子模拟的粒子快速定位方法
Stojanovic et al. High performance processing and analysis of geospatial data using CUDA on GPU
Wang et al. Spatial query based virtual reality GIS analysis platform
CN111080509B (zh) 三维裁剪Voronoi图的多线程并行计算方法、系统
Stojanovic et al. High–performance computing in GIS: Techniques and applications
CN114417075B (zh) 一种建立寻路网格数据索引的方法、装置、介质和设备
Xu et al. Feature-preserving simplification framework for 3D point cloud
Fellegara et al. Terrain trees: a framework for representing, analyzing and visualizing triangulated terrains
Stojanović et al. Performance improvement of viewshed analysis using GPU
CN113342999A (zh) 一种基于多层跳序树结构的变分辨率点云简化方法
Mateo et al. Hierarchical, Dense and Dynamic 3D Reconstruction Based on VDB Data Structure for Robotic Manipulation Tasks
Steinlechner et al. Adaptive pointcloud segmentation for assisted interactions

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant