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

Question

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

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

现在我正在使用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 包含三角剖分层次结构的实现， 其中

实现用数据结构增强的三角测量，以有效地回答点位置查询。[...]数据结构保持较小，并实现对真实数据的快速点位置查询。

其性能对于 Delaunay 三角剖分而言是最佳的。

         
          ^图36.8