自然邻域插值

Question

我为GRASS GIS编写了一个模块，使用CGAL提供自然邻居插值。我把测量点保存在向量points中。对于这些点我还有一个向量function_values: function_values.insert(std::make_pair(Point(x,y),z));

我想保存到一个网格值内插值的区域包含这些点。

到目前为止，我有工作的代码，但相当慢，因为我必须一步一步地通过颜色和行来定义点p K::Point_2 p(coor_x,coor_y)，并为每一步调用CGAL::natural_neighbor_coordinates_2。

/* perform NN interpolation */
G_message(_("Computing..."));
Delaunay_triangulation T;
typedef CGAL::Data_access< std::map<Point, Coord_type, K::Less_xy_2 > >
Value_access;
T.insert(points.begin(), points.end());
//coordinate computation in grid
double coor_x, coor_y;
coor_x = window.west;
coor_y = window.north;
for (int rows=0 ; rows<window.rows ; rows++) {
    G_percent(rows, window.rows, 2);
    coor_x = window.west;
    for (int cols=0 ; cols<window.cols ; cols++) {
        K::Point_2 p(coor_x,coor_y);
        std::vector< std::pair< Point, Coord_type > > coords;
        Coord_type norm = CGAL::natural_neighbor_coordinates_2(T, p,std::back_inserter(coords)).second;
        Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(), norm, Value_access(function_values));
        G_debug(5, "x: %f y: %f -> res: %f (row=%d; col=%d)",
        coor_x, coor_y, res, rows, cols);
        coor_x += ewres;
        std::cout << res << " ";
    }
    coor_y -= nsres;
    std::cout << std::endl;
}
G_percent(1, 1, 1);

是否有任何选项如何以不同的方式填充网格?

提前感谢亚当

Answer 1

主旨:

在自然邻居计算中，最昂贵的操作是确定哪个三角形包含查询点(locate()函数)。这个locate()函数可以通过给它一个"提示"来加速，也就是说，一个离包含该点的三角形不远的三角形(Dt::Face_handle)。由于所有的查询点都组织在一个规则的网格中，如果您使用包含前一个点的三角形作为提示，那么它将大大加快程序。

怎么做:

函数natural_neighbor_coordinates_2()有一个可选参数'start'来指定提示。

Delaunay_triangulation::Facet_handle fh; // the variable that stores the 'hint'
for(int rows ...) {
   for (int columns ...) {
       K::Point_2 p(...)
       ...
       fh = T.locate(p,fh); // locate p using the previous triangle as a hint
        Coord_type norm = CGAL::natural_neighbor_coordinates_2(T,p,std::back_inserter(coords),fh).second;     
        // Give the hint to the function that computes natural neighbour coordinates  
        ...
   }
}

走得更快:

CGAL是线程安全的，因此您可以使用OpenMP来加速事情(不要忘记将fh声明为私有，以便每个线程都有自己的副本:

#pragma omp parallel for private(fh)
for(int row ...) {
   for(int columns ...) {
   }
}