快速查找网格内围绕一个点的三角形的方法

5
我遇到了一个性能问题,需要完成一项任务。目前其中一个瓶颈是从非结构化网格中获取插值字段值。慢的部分是给定一个二维点和一个非结构化的二维网格,找到立即环绕该点的网格点。最好只需找到它所在的三角形。
目前我正在使用CGAL,但速度太慢了。使用当前实现,整个任务将在高端CPU上并行运行数天才能完成。
我认为慢的部分是CGAL::natural_neighbor_coordinates_2函数。
#ifndef FIELD_INTERPOLATOR_H
#define FIELD_INTERPOLATOR_H

#include "Vec.h"

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>

#include <map>
#include <vector>

typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
typedef CGAL::Delaunay_triangulation_2< Kernel > Delaunay_triangulation;

typedef Kernel::FT FieldType;
typedef Kernel::Point_2 MeshType;

struct FieldInterpolator23 {

    Delaunay_triangulation m_triangulation;

    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vX;
    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vY;
    std::map< MeshType, FieldType, Kernel::Less_xy_2 > m_vZ;

    typedef CGAL::Data_access< std::map< MeshType, FieldType, Kernel::Less_xy_2 > > ValueAccess;

    FieldInterpolator23() {}

    FieldInterpolator23( const std::vector< TN::Vec2 > & mesh, const std::vector< TN::Vec3 > & field )
    {
        const int N = mesh.size();
        for ( int i = 0; i < N; ++i ) {

            MeshType p( mesh[i].x(), mesh[i].y() );

            m_triangulation.insert( p );
            m_vX.insert( std::make_pair( p, field[i].x() ) );
            m_vY.insert( std::make_pair( p, field[i].y() ) ); 
            m_vZ.insert( std::make_pair( p, field[i].z() ) );                        
        }       
    }

    void set( const std::vector< TN::Vec2 > & mesh, const std::vector< TN::Vec3 > & field ) {

        m_triangulation.clear();
        m_vX.clear();
        m_vY.clear();
        m_vZ.clear();

        const int N = mesh.size();
        for ( int i = 0; i < N; ++i ) {

            MeshType p( mesh[i].x(), mesh[i].y() );

            m_triangulation.insert( p );
            m_vX.insert( std::make_pair( p, field[i].x() ) );
            m_vY.insert( std::make_pair( p, field[i].y() ) );
            m_vZ.insert( std::make_pair( p, field[i].z() ) );
        }
    }

    TN::Vec3 operator() ( TN::Vec2 p ) {

        MeshType pos( p.x(), p.y() );

        std::vector< std::pair< MeshType, FieldType > > coords;

        FieldType norm =
            CGAL::natural_neighbor_coordinates_2( m_triangulation, pos, std::back_inserter( coords ) ).second;

        FieldType resX =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vX )
        );

        FieldType resY =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vY )
        );

        FieldType resZ =
            CGAL::linear_interpolation(
                coords.begin(), 
                coords.end(),
                norm,
                ValueAccess( m_vZ )
        );

        return TN::Vec3( resX, resY, resZ );
    }
};

#endif

有没有人能够指导我一个更高效的解决方案,无论是不同的库还是算法?


1
你的三角剖分中有多少个顶点?你需要在三角剖分中定位多少个查询点? - Sneftel
1
大约有一百万个网格点,我需要进行大约20亿次查询。 - MVTC
5
如果您提前知道查询,将它们排序会有所帮助(比如说使用空间填充曲线,CGAL有相关函数),然后按照这个顺序进行点定位查询,并将之前的位置作为每个查询的起始点。这样大多数查询只需检查起始点是否已经是正确的三角形,其他查询最多只需要移动1或2步即可到达正确的三角形。 - Marc Glisse
如果您的查询点位于规则网格上,一种快速的方法是通过三角形填充:创建一个图像,其中每个像素对应一个查询点,并使用三角形ID作为填充颜色来填充每个三角形。这将使查询时间恒定。(对于如此大量的查询,这可能在您的情况下不切实际 - 除非您可以按子窗口进行操作。) - user1196549
3个回答

3
CGAL包含了一个三角剖分层次结构的实现,它实现了一个带有数据结构的三角剖分,以便高效地回答点定位查询。该数据结构保持小巧,并在真实数据上实现快速的点定位查询。它在Delaunay三角剖分中的性能是最优的。
Triangulation
图36.8

我改用了分层三角剖分,但仅导致了轻微的性能提升。作为一个测试,我尝试了天真地迭代整个网格节点集以找到最近的节点并返回其值。使用此方法,整个任务完成速度比使用CGAL三角剖分快约10倍。 - MVTC
@MVTC 听起来很奇怪,但我不太明白你在比较什么。你有一个基准测试可以分享吗?它只做点定位并显示朴素迭代比三角测量方法更快的结果吗? - Marc Glisse
@MVTC:很奇怪,暴力线性搜索比对数搜索表现更好。我猜您是把预处理放在一边,只看查询时间了吧? - Joseph O'Rourke
1
预处理被搁置了。此外,在我的应用程序中,使用locate函数查找包含该点的面比使用自然邻居函数快得多,这很奇怪,因为natural_neighbors_coord_2函数似乎只是在内部使用locate。时间差异如此之大,以至于似乎三角剖分出了问题?我需要更多时间来调查发生了什么。我将把这个作为接受的答案,并在弄清楚原因后更新。 - MVTC

2
如果网格是不变的,可以使用哈希表以O(1)时间访问任何条目。假设数据通过水平x和垂直y值访问,将一个边界框放在网格周围,并将其垂直和水平地切成大致相等于平均网格元素面积的正方形。将Nx设置为垂直切片的数量,将Ny设置为水平切片的数量。如果网格数据从X_min到X_max和Y_min到Y_max延伸,则将Sx = Nx / (X_max—Xmin)和Sy = Ny / (Y_max—Y_min)设置为。然后从左到右递增地对垂直列进行编号,从round(Sx*X_min)到round(Sx*X_max)开始。同样,从底部向上逐渐编号水平行,从round(Sy*Y_min)到round(Sy*Y_max)开始。
如果网格大致呈正方形,则大约有1000列和1000行,几乎每个正方形都与一个三角形相关联。如果网格具有不规则的形状或空洞,则许多正方形可能不会落在网格上。这对哈希表没有影响。要设置哈希表,每个正方形最多只能与一个三角形相关联。如果两个或多个三角形想要与同一个正方形关联,则正方形太大。要获得更多正方形,请设置较小的面积以获得较大的Nx和Ny。现在,每个三角形都与某一列#X和行#Y中的一个正方形相关联,可以生成100万个关键字。对于列#X和行#Y中的单元格,设置关键字字符串“±#X±#Y”,其中±#X和±#Y是具有前导“+”符号的整数,如果数字大于-1,则为前导“-”符号。如果数字小于零。为了使所有关键字具有相同的长度,请使用前导零填充整数值。典型的唯一关键字可能看起来像“-0378 + 0087”或“+0029-1007”。在这些示例中,每个关键字都恰好为10个字符。按存储三角形的顺序将关键字设置在哈希表中。
要使用哈希表,对于给定的x和y点,请设置整数ix = round(Sx * x)和整数iy = round(Sy * y),并根据ix和iy的符号形成关键字“±ix±iy”。如果x和y在网格中,则关键字“±ix±iy”应该在哈希表中,并且应该返回包围x和y或非常接近该点的三角形的索引。在某些情况下,生成的“±ix±iy”关键字可能不在哈希表中。在这种情况下,哈希表函数将返回一些随机三角形索引。当发生这种情况时,最好的解决方法是扰动ix和/或iy一个单位并生成另一个关键字。当然,如果该点在网格之外,则哈希表可能会发送回大量随机三角形索引。因此,请确保该点有效以开始。

本讨论假设三角形形状合理且面积大致相同。如果不是这种情况,需要将多个三角形适配到同一个正方形或多个正方形适配到同一个三角形中,则需要在侧边生成重定向表。当多个三角形适配到同一个正方形中时,将最中央的三角形分配给该正方形,并理解需要快速搜索所需的三角形。当多个正方形必须适配到同一个三角形中时,将一个正方形指向三角形索引,将其余所有正方形指向一个间接表,其索引从三角形索引之后开始。然后可以使用间接表来指向包含这些正方形的三角形索引。


2
根据我的经验,轴对齐边界框树在解决这个问题时是相当高效的(我观察到比“在三角剖分中行走”即使用locate()方法更快)。我正在使用自己的实现(针对3D网格,但可以很容易地适应2D):http://alice.loria.fr/software/geogram/doc/html/classGEO_1.1MeshFacetsAABB.html。CGAL也实现了AABB树:http://doc.cgal.org/latest/AABB_tree/
注意:如果您的查询点是有结构的(例如组织在规则网格上),那么有几种方式可以获得很高的性能,例如:
使用CGAL中Delaunay三角剖分的locate()函数,该函数接受一个提示参数,对于提示,使用前一个点相邻的三角形之一。这确保了点定位算法不必走得太远。这种方法通常非常有效。如果您不想更改类API,可以有一个可变成员来存储提示(如果点不是结构化的,则也可以工作,但需要在空间上进行排序,请参见Marc Glisse的评论)。
将三角形“绘制”到点上(每个三角形向其覆盖的点“发送”其值)。这可以通过用于在屏幕上绘制三角形的“计算机图形学”算法来完成。这种方法的收益更加重要(但需要更多的工作,它将完全改变您的算法)。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接