cgal natural neighbor interpolation

Question

I write a module for GRASS GIS that provide natural neighbor interpolation using CGAL. I have measured points saved in vector points. Also I have a vector function_values for these points: function_values.insert(std::make_pair(Point(x,y),z));

I'd like to save into a grid values interpolated within the region containing these points.

So far I have code that works but is quite slow, because I have to go throw the grid step-by-step by cols and rows to define point p K::Point_2 p(coor_x,coor_y) and call the CGAL::natural_neighbor_coordinates_2 for each single step.

/* 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);

Is there any option how to fill the grid in different way?

Thanks in advance Adam

score 1 · Answer 1 · answered Jul 24 '15 at 12:08

Main idea:

In the natural neighbors computation, the most expensive operation is determining which triangle contains the query point (locate() function). This locate() function can be accelerated by giving to it a "hint", that is to say a triangle (Dt::Face_handle) that is not too far away from the triangle that contains the point. Since all your query points are organized in a regular grid, if you use as a hint the triangle that contained the previous point, then it will significantly accelerate the program.

How to do it:

the function natural_neighbor_coordinates_2() has an optional parameter 'start' to specify the hint.

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  
        ...
   }
}

Going even faster:

CGAL is thread-safe, therefore you may use OpenMP to speed-up things (do not forget to declare fh as private so that each thread has its own copy:

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

cgal natural neighbor interpolation

1 Answers1