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

在给定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,Kernel::Less_xy_2 > m_vY;
    std::map< MeshType,Kernel::Less_xy_2 > m_vZ;

    typedef CGAL::Data_access< std::map< MeshType,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(),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(),ValueAccess( m_vY )
        );

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

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

#endif

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