To interpolate from a set of scattered points (e.g. gridding them onto regular grids), Delaunay_triangulation_2 is used to build the triangle mesh, natural_neighbor_coordinates_2() and linear_interpolation() are used for interpolation.
A problem I encountered is that when the input points are from some regular grids, the interpolation process could get "stuck" at some output location: the process is occupied by natural_neighbor_coordinates_2() but it never returns. It will run through if random noise is added onto the coordinates of the input points.
Wonder if anyone also had this problem and what is the solution. Adding random noise is OK but affects the accuracy of interpolation.
The scripts for interpolation are as below (I am using Armadillo for matrix)
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT Coord_type;
Delaunay_triangulation T;
arma::fmat points = ... ; //matrix saving point coordinates and function value
float output_x=...,output_y=...; //location for interpolation
std::map<Point, Coord_type, K::Less_xy_2> function_values;
// build mesh
for (long long i=0;i<points.n_cols;i++)
{
K::Point_2 p(points(0,i),points(1,i));
T.insert(p);
function_values.insert(std::make_pair(p,points(2,i)));
}
// interpolate
K::Point_2 p(output_x,output_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)); //res is the interpolation result.
This reminds me of a problem with random_polygon_2() getting stuck if the points are aligned, discussed here: http://cgal-discuss.949826.n4.nabble.com/random-polygon-2-gets-stuck-possible-CGAL-bug-td4659470.html
I'd suggest to try running your code with the set of points from Sebastien's answer. Your problem might also be related to points being aligned (just a guess).