在包围点的网格内找到三角形的快速方法

Question

我运行遇到了我需要完成的任务的性能问题。目前的瓶颈之一是从非结构化网格中获取插值字段值。

慢的部分是，给定一个二维点和一个非结构化的二维网格，找到紧邻该点的网格点。找到它落入的三角形就好了。

我现在正在使用 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

任何人都可以指出可接受的更高性能解决方案的方向，无论是不同的库还是算法？

Answer 1

CGAL 包含一个实现 Triangulation Hierarchy 其中

implements a triangulation augmented with a data structure to efficiently answer point location queries. [...] the data structure remains small and achieves fast point location queries on real data.

它的性能对于 Delaunay 三角剖分是最佳的。

^{图 36.8}

Answer 2

根据我自己的经验，Axis Aligned Bouding Box trees 对于这个问题相当有效（而且我观察到比 "walking in the triangulation" 更高效，即使用 locate()）。

我正在使用自己的实现（用于 3D 网格，但可以很容易地适应 2D）： http://alice.loria.fr/software/geogram/doc/html/classGEO_1_1MeshFacetsAABB.html

AABB 树也在 CGAL 中实现： http://doc.cgal.org/latest/AABB_tree/

注意：如果您的查询点是结构化的（例如按规则网格组织），那么有多种方法可以获得大量性能，例如：

使用 CGAL 中的 Delaunay 三角剖分的 locate() 函数，该函数将提示作为参数，对于提示，使用与前一点关联的三角形之一。这确保了点定位算法不必走得太远。这种方法的收益通常非常显着。如果你不想改变你的 class API，你可以有一个存储提示的可变成员（如果点没有结构化也可以工作，但需要在空间上对它们进行排序，另请参见 Marc Glisse 的评论）。
"drawing" 点上的三角形（每个三角形 "sends" 它的值到它覆盖的点）。这可以通过用于在屏幕上绘制三角形的 "computer graphics" 算法来完成。这种方法的收获更为重要（但需要更多的工作，它将彻底改变你的算法）。

Answer 3

如果网格是不变的，则可以使用哈希 table 在 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 行，几乎每个正方形都与一个三角形相关联。如果网格有不规则的形状或孔洞，许多方块可能不会落在网格上。这对散列 table 没有影响。要设置散列 table，每个正方形最多只能关联一个三角形。如果两个或多个三角形想要与同一个正方形相关联，则正方形太大。要获得更多方块，请设置更小的区域以获得更大的 Nx 和 Ny。现在每个三角形都与某一列 #X 和行 #Y 中的一个正方形相关联，可以生成一百万个关键字。对于列 #X 和行 #Y 中的单元格，设置关键字字符串“±#X±#Y”，其中 ±#X 和 ±#Y 表示整数，如果数字大于 —1，则带有前导“+”符号，并且如果数字小于零，则前导“—”符号。要为所有关键字提供相同的长度，请用前导零填充整数值。典型的唯一关键字可能类似于“-0378+0087”或“+0029-1007”。在这些示例中，每个关键字正好是 10 个字符。按照三角形的存储顺序设置散列 table 中的关键字。

要使用哈希 table，对于给定的 x 和 y 点，设置整数 ix = round(Sx*x) 和整数 iy = round(Sy*y) 并形成关键字“±ix ±iy" 根据 ix 和 iy 的符号。如果 x 和 y 在网格中，关键字“±ix±iy”应该在散列 table 中，它应该 return 包围 x 和 y 的三角形的索引或非常接近的三角形说到点子上。在某些情况下，生成的“±ix±iy”关键字可能不在散列 table 中。在这种情况下，散列 table 函数将 return 一些随机三角形索引。发生这种情况时，最好的解决方法是将 ix and/or iy 扰乱一个单位并生成另一个关键字。当然，如果该点在网格之外，散列 table 可能会发回大量随机三角形索引。因此，请确保该点有效。

本讨论假定三角形形状合理且面积近似均匀。如果不是这种情况，并且几个三角形必须适合同一个正方形或几个正方形必须适合同一个三角形，则解决方法是在侧面生成重定向 table。当几个三角形适合同一个正方形时，将最中心的三角形分配给正方形，并理解随后需要快速搜索所需的三角形。当几个正方形必须适合同一个三角形时，将其中一个正方形指向三角形索引，将所有其他正方形指向一个间接 table，其索引在三角形索引之后开始。间接 table 将用于指向包含这些正方形的三角形索引。

在包围点的网格内找到三角形的快速方法

Fast way to find the triangle inside a mesh that encloses a point

c++

performance

computational-geometry

cgal